grid4d.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
 *
 * Grid representation
 *
 ******************************************************************************/

#ifndef _GRID4D_H
#define _GRID4D_H

#include "manta.h"
#include "vectorbase.h"
#include "vector4d.h"
#include "kernel.h"


namespace Manta {

    
PYTHON class Grid4dBase : public PbClass {
public:
    enum Grid4dType { TypeNone = 0, TypeReal = 1, TypeInt = 2, TypeVec3 = 4, TypeVec4 = 8 };
        
    PYTHON() Grid4dBase(FluidSolver* parent);
    
    inline int getSizeX() const { return mSize.x; }
    inline int getSizeY() const { return mSize.y; }
    inline int getSizeZ() const { return mSize.z; }
    inline int getSizeT() const { return mSize.t; }
    inline Vec4i getSize() const { return mSize; }
    
    inline IndexInt getStrideX() const { return 1; }
    inline IndexInt getStrideY() const { return mSize.x; }
    inline IndexInt getStrideZ() const { return mStrideZ; }
    inline IndexInt getStrideT() const { return mStrideT; }
    
    inline Real getDx() { return mDx; }
    
    inline void checkIndex(int i, int j, int k, int t) const;
    inline void checkIndex(IndexInt idx) const;
    inline bool isInBounds(const Vec4i& p, int bnd) const;
    inline bool isInBounds(const Vec4i& p) const;
    inline bool isInBounds(const Vec4& p, int bnd = 0) const { return isInBounds(toVec4i(p), bnd); }
    inline bool isInBounds(IndexInt idx) const;
    
    inline Grid4dType getType() const { return mType; }
    inline bool is3D() const { return true; }
    inline bool is4D() const { return true; }

    inline bool isInBounds(int i,int j, int k, int t, int bnd) const { return isInBounds( Vec4i(i,j,k,t), bnd ); }
    
    inline IndexInt index(int i, int j, int k, int t) const { DEBUG_ONLY(checkIndex(i,j,k,t)); return (IndexInt)i + (IndexInt)mSize.x * j + (IndexInt)mStrideZ * k + (IndexInt)mStrideT * t; }
    inline IndexInt index(const Vec4i& pos) const    { DEBUG_ONLY(checkIndex(pos.x,pos.y,pos.z,pos.t)); return (IndexInt)pos.x + (IndexInt)mSize.x * pos.y + (IndexInt)mStrideZ * pos.z + (IndexInt)mStrideT * pos.t; }
protected:
    
    Grid4dType mType;
    Vec4i      mSize;
    Real       mDx;
    // precomputed Z,T shift: to ensure 2D compatibility, always use this instead of sx*sy !
    IndexInt   mStrideZ; 
    IndexInt   mStrideT; 
};

PYTHON template<class T> class Grid4d : public Grid4dBase {
public:
    PYTHON() Grid4d(FluidSolver* parent, bool show = true);
    Grid4d(const Grid4d<T>& a);
    virtual ~Grid4d();
    
    typedef T BASETYPE;
    typedef Grid4dBase BASETYPE_GRID;
    
    PYTHON() void save(std::string name);
    PYTHON() void load(std::string name);
    
    PYTHON() void clear();
    
    inline T get(int i,int j, int k, int t) const         { return mData[index(i,j,k,t)]; }
    inline T& get(int i,int j, int k, int t)              { return mData[index(i,j,k,t)]; }
    inline T get(IndexInt idx) const                           { DEBUG_ONLY(checkIndex(idx)); return mData[idx]; }
    inline T get(const Vec4i& pos) const                  { return mData[index(pos)]; }
    inline T& operator()(int i, int j, int k, int t)      { return mData[index(i, j, k,t)]; }
    inline T operator()(int i, int j, int k, int t) const { return mData[index(i, j, k,t)]; }
    inline T& operator()(IndexInt idx)                  { DEBUG_ONLY(checkIndex(idx)); return mData[idx]; }
    inline T operator()(IndexInt idx) const             { DEBUG_ONLY(checkIndex(idx)); return mData[idx]; }
    inline T& operator()(const Vec4i& pos)         { return mData[index(pos)]; }
    inline T operator()(const Vec4i& pos) const    { return mData[index(pos)]; }
    inline T& operator[](IndexInt idx)                  { DEBUG_ONLY(checkIndex(idx)); return mData[idx]; }
    inline const T operator[](IndexInt idx) const       { DEBUG_ONLY(checkIndex(idx)); return mData[idx]; }
    
    // interpolated access
    inline T    getInterpolated(const Vec4& pos) const { return interpol4d<T>(mData, mSize, mStrideZ, mStrideT, pos); }
    
    // assignment / copy

    //Grid4d<T>& operator=(const Grid4d<T>& a);
    PYTHON() Grid4d<T>& copyFrom(const Grid4d<T>& a, bool copyType=true ); // old: { *this = a; }

    // helper functions to work with grids in scene files 

    PYTHON() void add(const Grid4d<T>& a);
    PYTHON() void sub(const Grid4d<T>& a);
    PYTHON() void setConst(T s);
    PYTHON() void addConst(T s);
    PYTHON() void addScaled(const Grid4d<T>& a, const T& factor); 
    PYTHON() void mult( const Grid4d<T>& a);
    PYTHON() void multConst(T s);
    PYTHON() void clamp(Real min, Real max);
    
    // common compound operators
    PYTHON() Real getMaxAbs();
    PYTHON() Real getMax();
    PYTHON() Real getMin();    
    PYTHON() void setBound(T value, int boundaryWidth=1);
    PYTHON() void setBoundNeumann(int boundaryWidth=1);

    PYTHON() void printGrid(int zSlice=-1, int tSlice=-1,  bool printIndex=false, int bnd=0); 

    // c++ only operators
    template<class S> Grid4d<T>& operator+=(const Grid4d<S>& a);
    template<class S> Grid4d<T>& operator+=(const S& a);
    template<class S> Grid4d<T>& operator-=(const Grid4d<S>& a);
    template<class S> Grid4d<T>& operator-=(const S& a);
    template<class S> Grid4d<T>& operator*=(const Grid4d<S>& a);
    template<class S> Grid4d<T>& operator*=(const S& a);
    template<class S> Grid4d<T>& operator/=(const Grid4d<S>& a);
    template<class S> Grid4d<T>& operator/=(const S& a);
    Grid4d<T>& safeDivide(const Grid4d<T>& a);    
    
    void swap(Grid4d<T>& other);

protected:
    T* mData;
};

// Python doesn't know about templates: explicit aliases needed






inline Vec4 calcGridSizeFactor4d(Vec4i s1, Vec4i s2) {
    return Vec4( Real(s1[0])/s2[0], Real(s1[1])/s2[1], Real(s1[2])/s2[2] , Real(s1[3])/s2[3] );
}
inline Vec4 calcGridSizeFactor4d(Vec4 s1, Vec4 s2) {
    return Vec4( s1[0]/s2[0], s1[1]/s2[1], s1[2]/s2[2] , s1[3]/s2[3] );
}

// prototypes for grid plugins
void getComponent4d(const Grid4d<Vec4>& src, Grid4d<Real>& dst, int c);
void setComponent4d(const Grid4d<Real>& src, Grid4d<Vec4>& dst, int c);


//******************************************************************************
// Implementation of inline functions

inline void Grid4dBase::checkIndex(int i, int j, int k, int t) const {
    if ( i<0 || j<0  || i>=mSize.x || j>=mSize.y || k<0|| k>= mSize.z ||
         t<0|| t>= mSize.t ) {
        std::ostringstream s;
        s << "Grid4d " << mName << " dim " << mSize << " : index " << i << "," << j << "," << k << ","<<t<<" out of bound ";
        errMsg(s.str());
    }
}

inline void Grid4dBase::checkIndex(IndexInt idx) const {
    if (idx<0 || idx >= mSize.x * mSize.y * mSize.z * mSize.t) {
        std::ostringstream s;
        s << "Grid4d " << mName << " dim " << mSize << " : index " << idx << " out of bound ";
        errMsg(s.str());
    }
}

bool Grid4dBase::isInBounds(const Vec4i& p) const { 
    return (p.x >= 0 && p.y >= 0 && p.z >= 0 && p.t >= 0 && 
            p.x < mSize.x && p.y < mSize.y && p.z < mSize.z && p.t < mSize.t); 
}

bool Grid4dBase::isInBounds(const Vec4i& p, int bnd) const { 
    bool ret = (p.x >= bnd && p.y >= bnd && p.x < mSize.x-bnd && p.y < mSize.y-bnd);
    ret &= (p.z >= bnd && p.z < mSize.z-bnd); 
    ret &= (p.t >= bnd && p.t < mSize.t-bnd); 
    return ret;
}
bool Grid4dBase::isInBounds(IndexInt idx) const {
    if (idx<0 || idx >= mSize.x * mSize.y * mSize.z * mSize.t) {
        return false;
    }
    return true;
}

// note - ugly, mostly copied from normal GRID!

void Grid4dAdd (Grid4d<T>& me, const Grid4d<S>& other) {}

void Grid4dSub (Grid4d<T>& me, const Grid4d<S>& other) {}

void Grid4dMult (Grid4d<T>& me, const Grid4d<S>& other) {}

void Grid4dDiv (Grid4d<T>& me, const Grid4d<S>& other) {}

void Grid4dAddScalar (Grid4d<T>& me, const S& other) {}

void Grid4dMultScalar(Grid4d<T>& me, const S& other) {}

void Grid4dScaledAdd (Grid4d<T>& me, const Grid4d<T>& other, const S& factor) {}


void Grid4dSafeDiv (Grid4d<T>& me, const Grid4d<T>& other) {}

void Grid4dSetConst(Grid4d<T>& me, T value) {}


template<class T> template<class S> Grid4d<T>& Grid4d<T>::operator+= (const Grid4d<S>& a) {
    Grid4dAdd<T,S> (*this, a);
    return *this;
}
template<class T> template<class S> Grid4d<T>& Grid4d<T>::operator+= (const S& a) {
    Grid4dAddScalar<T,S> (*this, a);
    return *this;
}
template<class T> template<class S> Grid4d<T>& Grid4d<T>::operator-= (const Grid4d<S>& a) {
    Grid4dSub<T,S> (*this, a);
    return *this;
}
template<class T> template<class S> Grid4d<T>& Grid4d<T>::operator-= (const S& a) {
    Grid4dAddScalar<T,S> (*this, -a);
    return *this;
}
template<class T> template<class S> Grid4d<T>& Grid4d<T>::operator*= (const Grid4d<S>& a) {
    Grid4dMult<T,S> (*this, a);
    return *this;
}
template<class T> template<class S> Grid4d<T>& Grid4d<T>::operator*= (const S& a) {
    Grid4dMultScalar<T,S> (*this, a);
    return *this;
}
template<class T> template<class S> Grid4d<T>& Grid4d<T>::operator/= (const Grid4d<S>& a) {
    Grid4dDiv<T,S> (*this, a);
    return *this;
}
template<class T> template<class S> Grid4d<T>& Grid4d<T>::operator/= (const S& a) {
    S rez((S)1.0 / a);
    Grid4dMultScalar<T,S> (*this, rez);
    return *this;
}


//******************************************************************************
// Other helper functions

inline Vec4 getGradient4d(const Grid4d<Real>& data, int i, int j, int k, int t) {
    Vec4 v;
    if (i > data.getSizeX()-2) i= data.getSizeX()-2;
    if (j > data.getSizeY()-2) j= data.getSizeY()-2;
    if (k > data.getSizeZ()-2) k= data.getSizeZ()-2;
    if (t > data.getSizeT()-2) t= data.getSizeT()-2;
    if (i < 1) i = 1;
    if (j < 1) j = 1;
    if (k < 1) k = 1;
    if (t < 1) t = 1;
    v = Vec4( data(i+1,j  ,k  ,t  ) - data(i-1,j  ,k  ,t  ) ,
              data(i  ,j+1,k  ,t  ) - data(i  ,j-1,k  ,t  ) , 
              data(i  ,j  ,k+1,t  ) - data(i  ,j  ,k-1,t  ) , 
              data(i  ,j  ,k  ,t+1) - data(i  ,j  ,k  ,t-1) );
    return v;
}


void KnInterpolateGrid4dTempl(Grid4d<S>& target, Grid4d<S>& source, const Vec4& sourceFactor , Vec4 offset) {}
 

} //namespace
#endif