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