particle.h

Go to the documentation of this file.

/******************************************************************************
 *
 * MantaFlow fluid solver framework
 * Copyright 2011 Tobias Pfaff, Nils Thuerey 
 *
 * This program is free software, distributed under the terms of the
 * GNU General Public License (GPL) 
 * http://www.gnu.org/licenses
 *
 * Base class for particle systems
 *
 ******************************************************************************/

#ifndef _PARTICLE_H
#define _PARTICLE_H

#include <vector>
#include "grid.h"
#include "vectorbase.h"
#include "integrator.h"
#include "randomstream.h"
namespace Manta {

// fwd decl
template<class T> class Grid;
class ParticleDataBase;
template<class T> class ParticleDataImpl;

PYTHON class ParticleBase : public PbClass {
public:
    enum SystemType { BASE=0, PARTICLE, VORTEX, FILAMENT, FLIP, TURBULENCE, INDEX };
    
    enum ParticleStatus {
        PNONE         = 0,
        PNEW          = (1<<1),  // particles newly created in this step
        PDELETE       = (1<<10), // mark as deleted, will be deleted in next compress() step
        PINVALID      = (1<<30), // unused
    };

    PYTHON() ParticleBase(FluidSolver* parent);
    virtual ~ParticleBase();

    virtual void cloneParticleData(ParticleBase* nm);

    virtual SystemType getType() const { return BASE; }
    virtual std::string infoString() const; 
    virtual ParticleBase* clone() { assertMsg( false , "Dont use, override..."); return NULL; } 

    virtual IndexInt getSizeSlow() const { assertMsg( false , "Dont use, override..."); return 0; } 

    inline void addBuffered(const Vec3& pos);


    PYTHON() PbClass* create(PbType type, PbTypeVec T=PbTypeVec(), const std::string& name = "");
    void registerPdata(ParticleDataBase* pdata);
    void registerPdataReal(ParticleDataImpl<Real>* pdata);
    void registerPdataVec3(ParticleDataImpl<Vec3>* pdata);
    void registerPdataInt (ParticleDataImpl<int >* pdata);
    void deregister(ParticleDataBase* pdata);
    void addAllPdata();
    // note - deletion of pdata is handled in compress function

    IndexInt getNumPdata() const { return mPartData.size(); }
    ParticleDataBase* getPdata(int i) { return mPartData[i]; }

protected:  
    std::vector<Vec3> mNewBuffer;

    bool mAllowCompress;

    std::vector<ParticleDataBase*> mPartData;
    std::vector< ParticleDataImpl<Real> *> mPdataReal;
    std::vector< ParticleDataImpl<Vec3> *> mPdataVec3;
    std::vector< ParticleDataImpl<int> *>  mPdataInt;
    bool mFreePdata;
};



PYTHON template<class S> class ParticleSystem : public ParticleBase {
public:    
    PYTHON() ParticleSystem(FluidSolver* parent) : ParticleBase(parent), mDeletes(0), mDeleteChunk(0) {}
    virtual ~ParticleSystem() {};
    
    virtual SystemType getType() const { return S::getType(); };
    
    inline S& operator[](IndexInt idx)             { DEBUG_ONLY(checkPartIndex(idx)); return mData[idx]; }
    inline const S& operator[](IndexInt idx) const { DEBUG_ONLY(checkPartIndex(idx)); return mData[idx]; }
    inline IndexInt size() const { return mData.size(); }
    virtual IndexInt getSizeSlow() const { return size(); }
    PYTHON() int pySize() const { return (int)mData.size(); }

    inline int  getStatus(IndexInt idx) const { DEBUG_ONLY(checkPartIndex(idx)); return mData[idx].flag; }
    inline bool isActive(IndexInt idx) const  { DEBUG_ONLY(checkPartIndex(idx)); return (mData[idx].flag & PDELETE) == 0; }
    
    PYTHON() void setPos(const IndexInt idx, const Vec3& pos) { DEBUG_ONLY(checkPartIndex(idx)); mData[idx].pos = pos; }
    PYTHON() Vec3 getPos(IndexInt idx) const                  { DEBUG_ONLY(checkPartIndex(idx)); return mData[idx].pos; }
    PYTHON() void getPosPdata(ParticleDataImpl<Vec3>& target) const;
    PYTHON() void setPosPdata(const ParticleDataImpl<Vec3>& source);
    void transformPositions( Vec3i dimOld, Vec3i dimNew );

    void doCompress() { if ( mDeletes > mDeleteChunk) compress(); }
    void insertBufferedParticles();
    void resizeAll(IndexInt newsize);
    
    inline void kill(IndexInt idx);
    IndexInt add(const S& data);
    PYTHON() void clear();
            
    PYTHON() void advectInGrid( const FlagGrid& flags, const MACGrid& vel, const int integrationMode,
        const bool deleteInObstacle=true, const bool stopInObstacle=true,
        const ParticleDataImpl<int> *ptype=NULL, const int exclude=0);
    
    PYTHON() void projectOutside(Grid<Vec3>& gradient);
    PYTHON() void projectOutOfBnd(const FlagGrid &flags, const Real bnd, const std::string& plane="xXyYzZ", const ParticleDataImpl<int> *ptype=NULL, const int exclude=0);
    
    virtual ParticleBase* clone();
    virtual std::string infoString() const;

    inline void checkPartIndex(IndexInt idx) const;
    
protected:  
    IndexInt mDeletes, mDeleteChunk;    
    std::vector<S> mData;    

    virtual void compress(); 
};

//******************************************************************************

struct BasicParticleData {
public:
    BasicParticleData() : pos(0.), flag(0) {}
    BasicParticleData(const Vec3& p) : pos(p), flag(0) {}
    static ParticleBase::SystemType getType() { return ParticleBase::PARTICLE; }

    Vec3 pos;
    int  flag;
};

PYTHON class BasicParticleSystem : public ParticleSystem<BasicParticleData> {
public:
    PYTHON() BasicParticleSystem(FluidSolver* parent);
    
    PYTHON() void save(const std::string name) const;
    PYTHON() void load(const std::string name);

    void writeParticlesText(const std::string name) const;
    void writeParticlesRawPositionsGz(const std::string name) const;
    void writeParticlesRawVelocityGz(const std::string name) const;

    PYTHON() void readParticles(BasicParticleSystem* from);

    PYTHON() void addParticle(Vec3 pos) { add(BasicParticleData(pos)); }

    std::vector<BasicParticleData>& getData() { return mData; }

    PYTHON() void printParts(IndexInt start=-1, IndexInt stop=-1, bool printIndex=false); 
};


//******************************************************************************

//  used for grid based neighborhood searches on generic particle systems (stores
//  only active particles, and reduces copied data)
//  note - pos & flag are disabled here, do not use!
struct ParticleIndexData {
public:
    ParticleIndexData() : sourceIndex(0) {}
    static ParticleBase::SystemType getType() { return ParticleBase::INDEX; }

    IndexInt  sourceIndex; // index of this particle in the original particle system
    static Vec3 pos;  // do not use... 
    static int  flag; // not needed usally 
    //Vec3 pos; // enable for debugging
};

PYTHON class ParticleIndexSystem : public ParticleSystem<ParticleIndexData> {
public:
    PYTHON() ParticleIndexSystem(FluidSolver* parent) : ParticleSystem<ParticleIndexData>(parent) {};
    
    void resize(IndexInt size) { mData.resize(size); }
};



//******************************************************************************

PYTHON template<class DATA, class CON> class ConnectedParticleSystem : public ParticleSystem<DATA> {
public:
    PYTHON() ConnectedParticleSystem(FluidSolver* parent) : ParticleSystem<DATA>(parent) {}
    
    inline bool isSegActive(int i) { return (mSegments[i].flag & ParticleBase::PDELETE) == 0; }    
    inline int segSize() const { return mSegments.size(); }    
    inline CON& seg(int i) { return mSegments[i]; }
    inline const CON& seg(int i) const { return mSegments[i]; }
        
    virtual ParticleBase* clone();
    
protected:
    std::vector<CON> mSegments;
    virtual void compress();    
};

//******************************************************************************

PYTHON class ParticleDataBase : public PbClass {
public:
    PYTHON() ParticleDataBase(FluidSolver* parent);
    virtual ~ParticleDataBase(); 

    enum PdataType { TypeNone = 0, TypeReal = 1, TypeInt = 2, TypeVec3 = 4 };

    virtual IndexInt  getSizeSlow() const { assertMsg( false , "Dont use, override..."); return 0; } 
    virtual void addEntry()   { assertMsg( false , "Dont use, override..."); return;   }
    virtual ParticleDataBase* clone() { assertMsg( false , "Dont use, override..."); return NULL; }
    virtual PdataType getType() const { assertMsg( false , "Dont use, override..."); return TypeNone; } 
    virtual void resize(IndexInt size)     { assertMsg( false , "Dont use, override..."); return;  }
    virtual void copyValueSlow(IndexInt from, IndexInt to) { assertMsg( false , "Dont use, override..."); return;  }

    void setParticleSys(ParticleBase* set) { mpParticleSys = set; }

    inline void checkPartIndex(IndexInt idx) const;

protected:
    ParticleBase* mpParticleSys;
};


PYTHON template<class T> class ParticleDataImpl : public ParticleDataBase {
public:
    PYTHON() ParticleDataImpl(FluidSolver* parent);
    ParticleDataImpl(FluidSolver* parent, ParticleDataImpl<T>* other);
    virtual ~ParticleDataImpl();

    inline       T& get(IndexInt idx)              { DEBUG_ONLY(checkPartIndex(idx)); return mData[idx]; }
    inline const T& get(IndexInt idx) const        { DEBUG_ONLY(checkPartIndex(idx)); return mData[idx]; }
    inline       T& operator[](IndexInt idx)       { DEBUG_ONLY(checkPartIndex(idx)); return mData[idx]; }
    inline const T& operator[](IndexInt idx) const { DEBUG_ONLY(checkPartIndex(idx)); return mData[idx]; }

    PYTHON() void clear();

    PYTHON() void setSource(Grid<T>* grid, bool isMAC=false );

    virtual IndexInt  getSizeSlow() const;
    virtual void addEntry();
    virtual ParticleDataBase* clone();
    virtual PdataType getType() const;
    virtual void resize(IndexInt s);
    virtual void copyValueSlow(IndexInt from, IndexInt to);

    IndexInt  size() const { return mData.size(); }

    inline void copyValue(IndexInt from, IndexInt to) { get(to) = get(from); } 
    void initNewValue(IndexInt idx, Vec3 pos);

    PYTHON() void setConst(T s);
    PYTHON() void setConstRange(T s, const int begin, const int end);
    PYTHON() ParticleDataImpl<T>& copyFrom(const ParticleDataImpl<T>& a);
    PYTHON() void add(const ParticleDataImpl<T>& a);
    PYTHON() void sub(const ParticleDataImpl<T>& a);
    PYTHON() void addConst(T s);
    PYTHON() void addScaled(const ParticleDataImpl<T>& a, const T& factor); 
    PYTHON() void mult( const ParticleDataImpl<T>& a);
    PYTHON() void multConst(T s);
    PYTHON() void safeDiv(const ParticleDataImpl<T>& a);
    PYTHON() void clamp(Real min, Real max);
    PYTHON() void clampMin(Real vmin);
    PYTHON() void clampMax(Real vmax);

    PYTHON() Real getMaxAbs();
    PYTHON() Real getMax();
    PYTHON() Real getMin();

    PYTHON() T    sum(const ParticleDataImpl<int> *t=NULL, const int itype=0) const;
    PYTHON() Real sumSquare() const;
    PYTHON() Real sumMagnitude() const;

    PYTHON() void setConstIntFlag(T s, const ParticleDataImpl<int>& t, const int flag);

    PYTHON() void printPdata(IndexInt start=-1, IndexInt stop=-1, bool printIndex=false); 
    
    PYTHON() void save(const std::string name);
    PYTHON() void load(const std::string name);
protected:
    std::vector<T> mData; 

    Grid<T>* mpGridSource;
    bool mGridSourceMAC;
};






//******************************************************************************
// Implementation
//******************************************************************************

const int DELETE_PART = 20; // chunk size for compression

void ParticleBase::addBuffered(const Vec3& pos) {
    mNewBuffer.push_back(pos);
}
   
template<class S>
void ParticleSystem<S>::clear() {
    mDeleteChunk = mDeletes = 0;
    this->resizeAll(0); // instead of mData.clear
}

template<class S>
IndexInt ParticleSystem<S>::add(const S& data) {
    mData.push_back(data); 
    mDeleteChunk = mData.size() / DELETE_PART;
    this->addAllPdata();
    return mData.size()-1;
}

template<class S>
inline void ParticleSystem<S>::kill(IndexInt idx)     { 
    assertMsg(idx>=0 && idx<size(), "Index out of bounds");
    mData[idx].flag |= PDELETE; 
    if ( (++mDeletes > mDeleteChunk) && (mAllowCompress) ) compress(); 
}

template<class S>
void ParticleSystem<S>::getPosPdata(ParticleDataImpl<Vec3>& target) const {
    for(IndexInt i=0; i<(IndexInt)this->size(); ++i) {
        target[i] = this->getPos(i);
    }
}
template<class S>
void ParticleSystem<S>::setPosPdata(const ParticleDataImpl<Vec3>& source) {
    for(IndexInt i=0; i<(IndexInt)this->size(); ++i) {
        this->setPos(i, source[i]);
    }
}

template<class S>
void ParticleSystem<S>::transformPositions( Vec3i dimOld, Vec3i dimNew )
{
    Vec3 factor = calcGridSizeFactor( dimNew, dimOld );
    for(IndexInt i=0; i<(IndexInt)this->size(); ++i) {
        this->setPos(i, this->getPos(i) * factor );
    }
}

// check for deletion/invalid position, otherwise return velocity
template<class S> std::vector<Vec3> GridAdvectKernel (std::vector<S>& p, const MACGrid& vel, const FlagGrid& flags, Real dt, bool deleteInObstacle, bool stopInObstacle , const ParticleDataImpl<int> *ptype, const int exclude) {}
;

// final check after advection to make sure particles haven't escaped
// (similar to particle advection kernel)
void KnDeleteInObstacle(std::vector<S>& p, const FlagGrid& flags) {}


// try to get closer to actual obstacle boundary
static inline Vec3 bisectBacktracePos(const FlagGrid& flags, const Vec3& oldp, const Vec3& newp)
{
    Real s = 0.;
    for(int i=1; i<5; ++i) {
        Real ds = 1./(Real)(1<<i);
        if (!flags.isObstacle( oldp*(1.-(s+ds)) + newp*(s+ds) )) {
            s += ds;
        }
    }
    return( oldp*(1.-(s)) + newp*(s) );
}

// at least make sure all particles are inside domain
void KnClampPositions(std::vector<S>& p, const FlagGrid& flags, ParticleDataImpl<Vec3> *posOld = NULL, bool stopInObstacle=true, const ParticleDataImpl<int> *ptype=NULL, const int exclude=0) {}


// advection plugin
template<class S>
void ParticleSystem<S>::advectInGrid(const FlagGrid& flags, const MACGrid& vel, const int integrationMode,
                     const bool deleteInObstacle, const bool stopInObstacle,
                     const ParticleDataImpl<int> *ptype, const int exclude) {
    // position clamp requires old positions, backup
    ParticleDataImpl<Vec3> *posOld = NULL;
    if(!deleteInObstacle) {
        posOld = new ParticleDataImpl<Vec3>(this->getParent());
        posOld->resize(mData.size());
        for(IndexInt i=0; i<(IndexInt)mData.size();++i) (*posOld)[i] = mData[i].pos;
    }

    // update positions
    GridAdvectKernel<S> kernel(mData, vel, flags, getParent()->getDt(), deleteInObstacle, stopInObstacle, ptype, exclude );
    integratePointSet(kernel, integrationMode);

    if(!deleteInObstacle) {
        KnClampPositions<S>   (mData, flags, posOld , stopInObstacle, ptype, exclude );
        delete posOld;
    } else {
        KnDeleteInObstacle<S> (mData, flags);
    }
}

void KnProjectParticles(ParticleSystem<S>& part, Grid<Vec3>& gradient) {}


template<class S>
void ParticleSystem<S>::projectOutside(Grid<Vec3>& gradient) {
    KnProjectParticles<S>(*this, gradient);
}

void KnProjectOutOfBnd(ParticleSystem<S> &part, const FlagGrid &flags, const Real bnd, const bool *axis, const ParticleDataImpl<int> *ptype, const int exclude) {}


template<class S>
void ParticleSystem<S>::projectOutOfBnd(const FlagGrid &flags, const Real bnd, const std::string& plane, const ParticleDataImpl<int> *ptype, const int exclude) {
    bool axis[6] = { false };
    for(std::string::const_iterator it=plane.begin(); it!=plane.end(); ++it) {
        if(*it=='x') axis[0] = true;
        if(*it=='X') axis[1] = true;
        if(*it=='y') axis[2] = true;
        if(*it=='Y') axis[3] = true;
        if(*it=='z') axis[4] = true;
        if(*it=='Z') axis[5] = true;
    }
    KnProjectOutOfBnd<S>(*this, flags, bnd, axis, ptype, exclude);
}

template<class S>
void ParticleSystem<S>::resizeAll(IndexInt size) {
    // resize all buffers to target size in 1 go
    mData.resize(size);
    for(IndexInt i=0; i<(IndexInt)mPartData.size(); ++i)
        mPartData[i]->resize(size);
}

template<class S>
void ParticleSystem<S>::compress() {
    IndexInt nextRead = mData.size();
    for (IndexInt i=0; i<(IndexInt)mData.size(); i++) {
        while ((mData[i].flag & PDELETE) != 0) {
            nextRead--;
            mData[i] = mData[nextRead];
            // ugly, but prevent virtual function calls here:
            for(IndexInt pd=0; pd<(IndexInt)mPdataReal.size(); ++pd) mPdataReal[pd]->copyValue(nextRead, i);
            for(IndexInt pd=0; pd<(IndexInt)mPdataVec3.size(); ++pd) mPdataVec3[pd]->copyValue(nextRead, i);
            for(IndexInt pd=0; pd<(IndexInt)mPdataInt .size(); ++pd) mPdataInt [pd]->copyValue(nextRead, i);
            mData[nextRead].flag = PINVALID;
        }
    }
    if(nextRead<(IndexInt)mData.size()) debMsg("Deleted "<<((IndexInt)mData.size() - nextRead)<<" particles", 1); // debug info

    resizeAll(nextRead);
    mDeletes = 0;
    mDeleteChunk = mData.size() / DELETE_PART;
}

template<class S>
void ParticleSystem<S>::insertBufferedParticles() {
    if(mNewBuffer.size()==0) return;
    IndexInt newCnt = mData.size();
    resizeAll(newCnt + mNewBuffer.size());

    // clear new flag everywhere
    for(IndexInt i=0; i<(IndexInt)mData.size(); ++i) mData[i].flag &= ~PNEW;

    for(IndexInt i=0; i<(IndexInt)mNewBuffer.size(); ++i) {
        // note, other fields are not initialized here...
        mData[newCnt].pos  = mNewBuffer[i];
        mData[newCnt].flag = PNEW;
        // now init pdata fields from associated grids...
        for(IndexInt pd=0; pd<(IndexInt)mPdataReal.size(); ++pd) 
            mPdataReal[pd]->initNewValue(newCnt, mNewBuffer[i] );
        for(IndexInt pd=0; pd<(IndexInt)mPdataVec3.size(); ++pd) 
            mPdataVec3[pd]->initNewValue(newCnt, mNewBuffer[i] );
        for(IndexInt pd=0; pd<(IndexInt)mPdataInt.size(); ++pd) 
            mPdataInt[pd]->initNewValue(newCnt, mNewBuffer[i] );
        newCnt++;
    }
    if(mNewBuffer.size()>0) debMsg("Added & initialized "<<(IndexInt)mNewBuffer.size()<<" particles", 2); // debug info
    mNewBuffer.clear();
}


template<class DATA, class CON>
void ConnectedParticleSystem<DATA,CON>::compress() {
    const IndexInt sz = ParticleSystem<DATA>::size();
    IndexInt *renumber_back = new IndexInt[sz];
    IndexInt *renumber = new IndexInt[sz];
    for (IndexInt i=0; i<sz; i++)
        renumber[i] = renumber_back[i] = -1;
        
    // reorder elements
    std::vector<DATA>& data = ParticleSystem<DATA>::mData;
    IndexInt nextRead = sz;
    for (IndexInt i=0; i<nextRead; i++) {
        if ((data[i].flag & ParticleBase::PDELETE) != 0) {
            nextRead--;
            data[i] = data[nextRead];
            data[nextRead].flag = 0;           
            renumber_back[i] = nextRead;
        } else 
            renumber_back[i] = i;
    }
    
    // acceleration structure
    for (IndexInt i=0; i<nextRead; i++)
        renumber[renumber_back[i]] = i;
    
    // rename indices in filaments
    for (IndexInt i=0; i<(IndexInt)mSegments.size(); i++)
        mSegments[i].renumber(renumber);
        
    ParticleSystem<DATA>::mData.resize(nextRead);
    ParticleSystem<DATA>::mDeletes = 0;
    ParticleSystem<DATA>::mDeleteChunk = ParticleSystem<DATA>::size() / DELETE_PART;
    
    delete[] renumber;
    delete[] renumber_back;
}

template<class S>
ParticleBase* ParticleSystem<S>::clone() {
    ParticleSystem<S>* nm = new ParticleSystem<S>(getParent());
    if(this->mAllowCompress) compress();
    
    nm->mData = mData;
    nm->setName(getName());
    this->cloneParticleData(nm);
    return nm;
}

template<class DATA,class CON>
ParticleBase* ConnectedParticleSystem<DATA,CON>::clone() {
    ConnectedParticleSystem<DATA,CON>* nm = new ConnectedParticleSystem<DATA,CON>(this->getParent());
    if(this->mAllowCompress) compress();
    
    nm->mData = this->mData;
    nm->mSegments = mSegments;
    nm->setName(this->getName());
    this->cloneParticleData(nm);
    return nm;
}

template<class S>  
std::string ParticleSystem<S>::infoString() const { 
    std::stringstream s;
    s << "ParticleSys '" << getName() << "'\n-> ";
    if(this->getNumPdata()>0) s<< "pdata: "<< this->getNumPdata();
    s << "parts: " << size();
    //for(IndexInt i=0; i<(IndexInt)mPartData.size(); ++i) { sstr << i<<":" << mPartData[i]->size() <<" "; } 
    return s.str();
}
    
template<class S>  
inline void ParticleSystem<S>::checkPartIndex(IndexInt idx) const {
    IndexInt mySize = this->size();
    if (idx<0 || idx > mySize ) {
        errMsg( "ParticleBase " << " size " << mySize << " : index " << idx << " out of bound " );
    }
}
    
inline void ParticleDataBase::checkPartIndex(IndexInt idx) const {
    IndexInt mySize = this->getSizeSlow();
    if (idx<0 || idx > mySize ) {
        errMsg( "ParticleData " << " size " << mySize << " : index " << idx << " out of bound " );
    }
    if ( mpParticleSys && mpParticleSys->getSizeSlow()!=mySize ) {
        errMsg( "ParticleData " << " size " << mySize << " does not match parent! (" << mpParticleSys->getSizeSlow() << ") " );
    }
}

// set contents to zero, as for a grid
template<class T>
void ParticleDataImpl<T>::clear() {
    for(IndexInt i=0; i<(IndexInt)mData.size(); ++i) mData[i] = 0.;
}

int countParticles(const ParticleDataImpl<int> &t, const int flag);

} // namespace

#endif