isosurface.h

Go to the documentation of this file.

/******************************************************************************
 *
 * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
 * Copyright 2003-2006 Nils Thuerey
 *
 * Marching Cubes "displayer"
 *
 *****************************************************************************/

#ifndef ISOSURFACE_H

#include "ntl_geometryobject.h"
#include "ntl_bsptree.h"

#define ISO_STRICT_DEBUG 0
#define ISOSTRICT_EXIT *((int *)0)=0;

/* access some 3d array */
#define ISOLEVEL_INDEX(ii,ij,ik) ((mSizex*mSizey*(ik))+(mSizex*(ij))+((ii)))

class ParticleTracer;

/* struct for a small cube in the scalar field */
typedef struct {
  ntlVec3Gfx   pos[8];
  double  value[8];
  int i,j,k;
} IsoLevelCube;


typedef struct {
  ntlVec3Gfx v; // vertex
  ntlVec3Gfx n; // vertex normal
} IsoLevelVertex;

// the fluid surface, templated by scalar field access object 
class IsoSurface : 
    public ntlGeometryObject //, public S
{

    public:

        IsoSurface(double iso);
        virtual ~IsoSurface();

        virtual void initializeIsosurface(int setx, int sety, int setz, ntlVec3Gfx extent);

        void resetAll(gfxReal val);

        void triangulate( void );

        void setParticles(ParticleTracer *pnt,float psize){ mpIsoParts = pnt; mPartSize=psize; };
        void setSubdivs(int s) { 
            if(mInitDone) errFatal("IsoSurface::setSubdivs","Changing subdivs after init!", SIMWORLD_INITERROR);
            if(s<1) s=1; if(s>10) s=10;
            mSubdivs = s; }
        int  getSubdivs() { return mSubdivs;}
        void setUseFulledgeArrays(bool set) { 
            if(mInitDone) errFatal("IsoSurface::setUseFulledgeArrays","Changing usefulledge after init!", SIMWORLD_INITERROR);
            mUseFullEdgeArrays = set;}

    protected:

        /* variables ... */

        int mSizex, mSizey, mSizez;

        float *mpData;

        double mIsoValue;

        vector<IsoLevelVertex> mPoints;

        bool mUseFullEdgeArrays;
        int *mpEdgeVerticesX;
        int *mpEdgeVerticesY;
        int *mpEdgeVerticesZ;
        int mEdgeArSize;


        vector<unsigned int> mIndices;

        ntlVec3Gfx mStart, mEnd;

        ntlVec3Gfx mDomainExtent;

        bool mInitDone;

        float mSmoothSurface;
        float mSmoothNormals;
        
        vector<int> mAcrossEdge;
        vector< vector<int> > mAdjacentFaces;

        int mCutoff;
        int *mCutArray;
        ParticleTracer *mpIsoParts;
        float mPartSize;
        int mSubdivs;
        
        vector<int> flags;
        int mFlagCnt;
        vector<ntlVec3Gfx> cornerareas;
        vector<float> pointareas;
        vector< vector<int> > neighbors;

    public:
        // miscelleanous access functions 

        void setStart(ntlVec3Gfx set) { mStart = set; };
        ntlVec3Gfx getStart() { return mStart; };
        void setEnd(ntlVec3Gfx set) { mEnd = set; };
        ntlVec3Gfx getEnd() { return mEnd; };
        inline void setIsolevel(double set) { mIsoValue = set; };
        inline void setSmoothSurface(float set) { mSmoothSurface = set; };
        inline void setSmoothNormals(float set) { mSmoothNormals = set; };
        inline float getSmoothSurface() { return mSmoothSurface; }
        inline float getSmoothNormals() { return mSmoothNormals; }

        // geometry object functions 
        virtual void getTriangles(double t, vector<ntlTriangle> *triangles, 
                vector<ntlVec3Gfx> *vertices, 
                vector<ntlVec3Gfx> *normals, int objectId );

        virtual inline ntlVec3Gfx *getBBStart()     { return &mStart; }
        virtual inline ntlVec3Gfx *getBBEnd()       { return &mEnd; }

        inline float* getData(){ return mpData; }
        inline float* getData(int ii, int jj, int kk){ 
#if ISO_STRICT_DEBUG==1
            if(ii<0){ errMsg("IsoStrict"," invX- |"<<ii<<","<<jj<<","<<kk); ISOSTRICT_EXIT; }
            if(jj<0){ errMsg("IsoStrict"," invY- |"<<ii<<","<<jj<<","<<kk); ISOSTRICT_EXIT; }
            if(kk<0){ errMsg("IsoStrict"," invZ- |"<<ii<<","<<jj<<","<<kk); ISOSTRICT_EXIT; }
            if(ii>mSizex-1){ errMsg("IsoStrict"," invX+ |"<<ii<<","<<jj<<","<<kk); ISOSTRICT_EXIT; }
            if(jj>mSizey-1){ errMsg("IsoStrict"," invY+ |"<<ii<<","<<jj<<","<<kk); ISOSTRICT_EXIT; }
            if(kk>mSizez-1){ errMsg("IsoStrict"," invZ+ |"<<ii<<","<<jj<<","<<kk); ISOSTRICT_EXIT; }
            return mpData + ISOLEVEL_INDEX(ii, jj, kk); 
#else //ISO_STRICT_DEBUG==1
            return mpData + ISOLEVEL_INDEX(ii, jj, kk); 
#endif
        }
        inline float* lbmGetData(int ii, int jj, int kk){ 
#if ISO_STRICT_DEBUG==1
            ii++; jj++; kk++;
            if(ii<0){ errMsg("IsoStrict"," invX- |"<<ii<<","<<jj<<","<<kk); ISOSTRICT_EXIT; }
            if(jj<0){ errMsg("IsoStrict"," invY- |"<<ii<<","<<jj<<","<<kk); ISOSTRICT_EXIT; }
            if(kk<0){ errMsg("IsoStrict"," invZ- |"<<ii<<","<<jj<<","<<kk); ISOSTRICT_EXIT; }
            if(ii>mSizex-1){ errMsg("IsoStrict"," invX+ |"<<ii<<","<<jj<<","<<kk); ISOSTRICT_EXIT; }
            if(jj>mSizey-1){ errMsg("IsoStrict"," invY+ |"<<ii<<","<<jj<<","<<kk); ISOSTRICT_EXIT; }
            if(kk>mSizez-1){ errMsg("IsoStrict"," invZ+ |"<<ii<<","<<jj<<","<<kk); ISOSTRICT_EXIT; }
            return mpData + ISOLEVEL_INDEX(ii, jj, kk); 
#else //ISO_STRICT_DEBUG==1
            return mpData + ISOLEVEL_INDEX(ii+1,jj+1,kk+1); 
#endif
        }
        inline void setCutoff(int set) { mCutoff = set; };
        inline void setCutArray(int *set) { mCutArray = set; };

        unsigned int getIsoVertexCount() {
            return mPoints.size();
        }
        unsigned int getIsoIndexCount() {
            return mIndices.size();
        }
        char* getIsoVertexArray() {
            return (char *) &(mPoints[0]);
        }
        unsigned int *getIsoIndexArray() {
            return &(mIndices[0]);
        }
        
        // surface smoothing functions
        void setSmoothRad(float radi1, float radi2, ntlVec3Gfx mscc);
        void smoothSurface(float val, bool smoothNorm);
        void smoothNormals(float val);
        void computeNormals();

    protected:

        inline ntlVec3Gfx getNormal(int i, int j,int k);
        bool diffuseVertexField(ntlVec3Gfx *field, int pointerScale, int v, float invsigma2, ntlVec3Gfx &flt);
        vector<int> mDboundary;
        float mSCrad1, mSCrad2;
        ntlVec3Gfx mSCcenter;
};


#define ISOSURFACE_H
#endif