particletracer.cpp

Go to the documentation of this file.

/******************************************************************************
 *
 * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
 * Copyright 2003-2006 Nils Thuerey
 *
 * Particle Viewer/Tracer
 *
 *****************************************************************************/

#include <stdio.h>
//#include "../libs/my_gl.h"
//#include "../libs/my_glu.h"

/* own lib's */
#include "particletracer.h"
#include "ntl_matrices.h"
#include "ntl_ray.h"
#include "ntl_matrices.h"

#include <zlib.h>


// particle object id counter
int ParticleObjectIdCnt = 1;

/******************************************************************************
 * Standard constructor
 *****************************************************************************/
ParticleTracer::ParticleTracer() :
    ntlGeometryObject(),
    mParts(),
    //mTrailLength(1), mTrailInterval(1),mTrailIntervalCounter(0),
    mPartSize(0.01),
    mStart(-1.0), mEnd(1.0),
    mSimStart(-1.0), mSimEnd(1.0),
    mPartScale(0.1) , mPartHeadDist( 0.1 ), mPartTailDist( -0.1 ), mPartSegments( 4 ),
    mValueScale(0),
    mValueCutoffTop(0.0), mValueCutoffBottom(0.0),
    mDumpParts(0), //mDumpText(0), 
    mDumpTextFile(""), 
    mDumpTextInterval(0.), mDumpTextLastTime(0.), mDumpTextCount(0),
    mShowOnly(0), 
    mNumInitialParts(0), mpTrafo(NULL),
    mInitStart(-1.), mInitEnd(-1.),
    mPrevs(), mTrailTimeLast(0.), mTrailInterval(-1.), mTrailLength(0)
{
    debMsgStd("ParticleTracer::ParticleTracer",DM_MSG,"inited",10);
};

ParticleTracer::~ParticleTracer() {
    debMsgStd("ParticleTracer::~ParticleTracer",DM_MSG,"destroyed",10);
}

/*****************************************************************************/
/*****************************************************************************/
void ParticleTracer::parseAttrList(AttributeList *att) 
{
    AttributeList *tempAtt = mpAttrs; 
    mpAttrs = att;

    mNumInitialParts = mpAttrs->readInt("particles",mNumInitialParts, "ParticleTracer","mNumInitialParts", false);
    //errMsg(" NUMP"," "<<mNumInitialParts);
    mPartScale    = mpAttrs->readFloat("part_scale",mPartScale, "ParticleTracer","mPartScale", false);
    mPartHeadDist = mpAttrs->readFloat("part_headdist",mPartHeadDist, "ParticleTracer","mPartHeadDist", false);
    mPartTailDist = mpAttrs->readFloat("part_taildist",mPartTailDist, "ParticleTracer","mPartTailDist", false);
    mPartSegments = mpAttrs->readInt  ("part_segments",mPartSegments, "ParticleTracer","mPartSegments", false);
    mValueScale   = mpAttrs->readInt  ("part_valscale",mValueScale, "ParticleTracer","mValueScale", false);
    mValueCutoffTop = mpAttrs->readFloat("part_valcutofftop",mValueCutoffTop, "ParticleTracer","mValueCutoffTop", false);
    mValueCutoffBottom = mpAttrs->readFloat("part_valcutoffbottom",mValueCutoffBottom, "ParticleTracer","mValueCutoffBottom", false);

    mDumpParts   = mpAttrs->readInt  ("part_dump",mDumpParts, "ParticleTracer","mDumpParts", false);
    // mDumpText deprecatd, use mDumpTextInterval>0. instead
    mShowOnly    = mpAttrs->readInt  ("part_showonly",mShowOnly, "ParticleTracer","mShowOnly", false);
    mDumpTextFile= mpAttrs->readString("part_textdumpfile",mDumpTextFile, "ParticleTracer","mDumpTextFile", false);
    mDumpTextInterval= mpAttrs->readFloat("part_textdumpinterval",mDumpTextInterval, "ParticleTracer","mDumpTextInterval", false);

    string matPart;
    matPart = mpAttrs->readString("material_part", "default", "ParticleTracer","material", false);
    setMaterialName( matPart );

    mInitStart = mpAttrs->readFloat("part_initstart",mInitStart, "ParticleTracer","mInitStart", false);
    mInitEnd   = mpAttrs->readFloat("part_initend",  mInitEnd, "ParticleTracer","mInitEnd", false);

    // unused...
    //int mTrailLength  = 0; // UNUSED
    //int mTrailInterval= 0; // UNUSED
    mTrailLength  = mpAttrs->readInt("traillength",mTrailLength, "ParticleTracer","mTrailLength", false);
    mTrailInterval= mpAttrs->readFloat("trailinterval",mTrailInterval, "ParticleTracer","mTrailInterval", false);

    // restore old list
    mpAttrs = tempAtt;
}

/******************************************************************************
 * draw the particle array
 *****************************************************************************/
void ParticleTracer::draw()
{
}

/******************************************************************************
 * init trafo matrix
 *****************************************************************************/
void ParticleTracer::initTrafoMatrix() {
    ntlVec3Gfx scale = ntlVec3Gfx(
            (mEnd[0]-mStart[0])/(mSimEnd[0]-mSimStart[0]),
            (mEnd[1]-mStart[1])/(mSimEnd[1]-mSimStart[1]),
            (mEnd[2]-mStart[2])/(mSimEnd[2]-mSimStart[2])
            );
    ntlVec3Gfx trans = mStart;
    if(!mpTrafo) mpTrafo = new ntlMat4Gfx(0.0);
    mpTrafo->initId();
    for(int i=0; i<3; i++) { mpTrafo->value[i][i] = scale[i]; }
    for(int i=0; i<3; i++) { mpTrafo->value[i][3] = trans[i]; }
}

/******************************************************************************
 * adapt time step by rescaling velocities
 *****************************************************************************/
void ParticleTracer::adaptPartTimestep(float factor) {
    for(size_t i=0; i<mParts.size(); i++) {
        mParts[i].setVel( mParts[i].getVel() * factor );
    }
} 


/******************************************************************************
 * add a particle at this position
 *****************************************************************************/
void ParticleTracer::addParticle(float x, float y, float z) {
    ntlVec3Gfx p(x,y,z);
    ParticleObject part( p );
    mParts.push_back( part );
}
void ParticleTracer::addFullParticle(ParticleObject &np) {
    mParts.push_back( np );
}


void ParticleTracer::cleanup() {
    // cleanup
    int last = (int)mParts.size()-1;
    if(mDumpTextInterval>0.) { errMsg("ParticleTracer::cleanup","Skipping cleanup due to text dump..."); return; }

    for(int i=0; i<=last; i++) {
        if( mParts[i].getActive()==false ) {
            ParticleObject *p = &mParts[i];
            ParticleObject *p2 = &mParts[last];
            *p = *p2; last--; mParts.pop_back();
        }
    }
}
        
extern bool glob_mpactive;
extern int glob_mpindex,glob_mpnum;

/******************************************************************************
 *! dump particles if desired 
 *****************************************************************************/
void ParticleTracer::notifyOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename, double simtime) {
    debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"obj:"<<this->getName()<<" frame:"<<frameNrStr<<" dumpp"<<mDumpParts<<" t"<<simtime, 10); // DEBUG

    if(
            (dumptype==DUMP_FULLGEOMETRY)&&
            (mDumpParts>0)) {
        // dump to binary file
        std::ostringstream boutfilename("");
        boutfilename << outfilename <<"_particles_" << frameNrStr;
        if(glob_mpactive) {
            if(glob_mpindex>0) { boutfilename << "mp"<<glob_mpindex; }
        }
        boutfilename << ".gz";
        debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"B-Dumping: "<< this->getName() <<", particles:"<<mParts.size()<<" "<< " to "<<boutfilename.str()<<" #"<<frameNr , 7);
        //debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"B-Dumping: partgeodeb sim:"<<mSimStart<<","<<mSimEnd<<" geosize:"<<mStart<<","<<mEnd,2 );

        // output to zipped file
        gzFile gzf;
        gzf = gzopen(boutfilename.str().c_str(), "wb1");
        if(gzf) {
            int numParts;
            if(sizeof(numParts)!=4) { errMsg("ParticleTracer::notifyOfDump","Invalid int size"); return; }
            // only dump active particles
            numParts = 0;
            for(size_t i=0; i<mParts.size(); i++) {
                if(!mParts[i].getActive()) continue;
                numParts++;
            }
            gzwrite(gzf, &numParts, sizeof(numParts));
            for(size_t i=0; i<mParts.size(); i++) {
                if(!mParts[i].getActive()) { continue; }
                ParticleObject *p = &mParts[i];
                //int type = p->getType();  // export whole type info
                int type = p->getFlags(); // debug export whole type & status info
                ntlVec3Gfx pos = p->getPos();
                float size = p->getSize();

                if(type&PART_FLOAT) { // WARNING same handling for dump!
                    // add one gridcell offset
                    //pos[2] += 1.0; 
                } 
                // display as drop for now externally
                //else if(type&PART_TRACER) { type |= PART_DROP; }

                pos = (*mpTrafo) * pos;

                ntlVec3Gfx v = p->getVel();
                v[0] *= mpTrafo->value[0][0];
                v[1] *= mpTrafo->value[1][1];
                v[2] *= mpTrafo->value[2][2];
                // FIXME check: pos = (*mpTrafo) * pos;
                gzwrite(gzf, &type, sizeof(type)); 
                gzwrite(gzf, &size, sizeof(size)); 
                for(int j=0; j<3; j++) { gzwrite(gzf, &pos[j], sizeof(float)); }
                for(int j=0; j<3; j++) { gzwrite(gzf, &v[j], sizeof(float)); }
            }
            gzclose( gzf );
        }
    } // dump?
}

void ParticleTracer::checkDumpTextPositions(double simtime) {
    // dfor partial & full dump
    if(mDumpTextInterval>0.) {
        debMsgStd("ParticleTracer::checkDumpTextPositions",DM_MSG,"t="<<simtime<<" last:"<<mDumpTextLastTime<<" inter:"<<mDumpTextInterval,7);
    }

    if((mDumpTextInterval>0.) && (simtime>mDumpTextLastTime+mDumpTextInterval)) {
        // dump to binary file
        std::ostringstream boutfilename("");
        if(mDumpTextFile.length()>1) {   
            boutfilename << mDumpTextFile <<  ".cpart2"; 
        } else {                           
            boutfilename << boutfilename <<"_particles" <<  ".cpart2"; 
        }
        debMsgStd("ParticleTracer::checkDumpTextPositions",DM_MSG,"T-Dumping: "<< this->getName() <<", particles:"<<mParts.size()<<" "<< " to "<<boutfilename.str()<<" " , 7);

        int numParts = 0;
        // only dump bubble particles
        for(size_t i=0; i<mParts.size(); i++) {
            //if(!mParts[i].getActive()) continue;
            //if(!(mParts[i].getType()&PART_BUBBLE)) continue;
            numParts++;
        }

        // output to text file
        //gzFile gzf;
        FILE *stf;
        if(mDumpTextCount==0) {
            //gzf = gzopen(boutfilename.str().c_str(), "w0");
            stf = fopen(boutfilename.str().c_str(), "w");

            fprintf( stf, "\n\n# cparts generated by elbeem \n# no. of parts \nN %d \n\n",numParts);
            // fixed time scale for now
            fprintf( stf, "T %f \n\n", 1.0);
        } else {
            //gzf = gzopen(boutfilename.str().c_str(), "a+0");
            stf = fopen(boutfilename.str().c_str(), "a+");
        }

        // add current set
        if(stf) {
            fprintf( stf, "\n\n# new set at frame %d,t%f,p%d --------------------------------- \n\n", mDumpTextCount, simtime, numParts );
            fprintf( stf, "S %f \n\n", simtime );
            
            for(size_t i=0; i<mParts.size(); i++) {
                ParticleObject *p = &mParts[i];
                ntlVec3Gfx pos = p->getPos();
                float size = p->getSize();
                float infl = 1.;
                //if(!mParts[i].getActive()) { size=0.; } // switch "off"
                if(!mParts[i].getActive()) { infl=0.; } // switch "off"
                if(!mParts[i].getInFluid()) { infl=0.; } // switch "off"
                if(mParts[i].getLifeTime()<0.) { infl=0.; } // not yet active...

                pos = (*mpTrafo) * pos;
                ntlVec3Gfx v = p->getVel();
                v[0] *= mpTrafo->value[0][0];
                v[1] *= mpTrafo->value[1][1];
                v[2] *= mpTrafo->value[2][2];
                
                fprintf( stf, "P %f %f %f \n", pos[0],pos[1],pos[2] );
                if(size!=1.0) fprintf( stf, "s %f \n", size );
                if(infl!=1.0) fprintf( stf, "i %f \n", infl );
                fprintf( stf, "\n" );
            }

            fprintf( stf, "# %d end  ", mDumpTextCount );
            //gzclose( gzf );
            fclose( stf );

            mDumpTextCount++;
        }

        mDumpTextLastTime += mDumpTextInterval;
    }

}


void ParticleTracer::checkTrails(double time) {
    if(mTrailLength<1) return;
    if(time-mTrailTimeLast > mTrailInterval) {

        if( (int)mPrevs.size() < mTrailLength) mPrevs.resize( mTrailLength );
        for(int i=mPrevs.size()-1; i>0; i--) {
            mPrevs[i] = mPrevs[i-1];
            //errMsg("TRAIL"," from "<<i<<" to "<<(i-1) );
        }
        mPrevs[0] = mParts;

        mTrailTimeLast += mTrailInterval;
    }
}


/******************************************************************************
 * Get triangles for rendering
 *****************************************************************************/
void ParticleTracer::getTriangles(double time, vector<ntlTriangle> *triangles, 
                                                     vector<ntlVec3Gfx> *vertices, 
                                                     vector<ntlVec3Gfx> *normals, int objectId )
{
#ifdef ELBEEM_PLUGIN
    // suppress warnings...
    vertices = NULL; triangles = NULL;
    normals = NULL; objectId = 0;
    time = 0.;
#else // ELBEEM_PLUGIN
    int pcnt = 0;
    // currently not used in blender
    objectId = 0; // remove, deprecated
    if(mDumpParts>1) { 
        return; // only dump, no tri-gen
    }

    const bool debugParts = false;
    int tris = 0;
    int segments = mPartSegments;
    ntlVec3Gfx scale = ntlVec3Gfx( (mEnd[0]-mStart[0])/(mSimEnd[0]-mSimStart[0]), (mEnd[1]-mStart[1])/(mSimEnd[1]-mSimStart[1]), (mEnd[2]-mStart[2])/(mSimEnd[2]-mSimStart[2]));
    ntlVec3Gfx trans = mStart;
    time = 0.; // doesnt matter

    for(size_t t=0; t<mPrevs.size()+1; t++) {
        vector<ParticleObject> *dparts;
        if(t==0) {
            dparts = &mParts;
        } else {
            dparts = &mPrevs[t-1];
        }
        //errMsg("TRAILT","prevs"<<t<<"/"<<mPrevs.size()<<" parts:"<<dparts->size() );

    gfxReal partscale = mPartScale;
    if(t>1) { 
        partscale *= (gfxReal)(mPrevs.size()+1-t) / (gfxReal)(mPrevs.size()+1); 
    }
    gfxReal partNormSize = 0.01 * partscale;
    //for(size_t i=0; i<mParts.size(); i++) {
    for(size_t i=0; i<dparts->size(); i++) {
        ParticleObject *p = &( (*dparts)[i] ); //  mParts[i];

        if(mShowOnly!=10) {
            // 10=show only deleted
            if( p->getActive()==false ) continue;
        } else {
            if( p->getActive()==true ) continue;
        }
        int type = p->getType();
        if(mShowOnly>0) {
            switch(mShowOnly) {
            case 1: if(!(type&PART_BUBBLE)) continue; break;
            case 2: if(!(type&PART_DROP))   continue; break;
            case 3: if(!(type&PART_INTER))  continue; break;
            case 4: if(!(type&PART_FLOAT))  continue; break;
            case 5: if(!(type&PART_TRACER))  continue; break;
            }
        } else {
            // by default dont display inter
            if(type&PART_INTER) continue;
        }

        pcnt++;
        ntlVec3Gfx pnew = p->getPos();
        if(type&PART_FLOAT) { // WARNING same handling for dump!
            if(p->getStatus()&PART_IN) { pnew[2] += 0.8; } // offset for display
            // add one gridcell offset
            //pnew[2] += 1.0; 
        }
#if LBMDIM==2
        pnew[2] += 0.001; // DEBUG
        pnew[2] += 0.009; // DEBUG
#endif 

        ntlVec3Gfx pdir = p->getVel();
        gfxReal plen = normalize( pdir );
        if( plen < 1e-05) pdir = ntlVec3Gfx(-1.0 ,0.0 ,0.0);
        ntlVec3Gfx pos = (*mpTrafo) * pnew;
        gfxReal partsize = 0.0;
        if(debugParts) errMsg("DebugParts"," i"<<i<<" new"<<pnew<<" vel"<<pdir<<"   pos="<<pos );
        //if(i==0 &&(debugParts)) errMsg("DebugParts"," i"<<i<<" new"<<pnew[0]<<" pos="<<pos[0]<<" scale="<<scale[0]<<"  t="<<trans[0] );
        
        // value length scaling?
        if(mValueScale==1) {
            partsize = partscale * plen;
        } else if(mValueScale==2) {
            // cut off scaling
            if(plen > mValueCutoffTop) continue;
            if(plen < mValueCutoffBottom) continue;
            partsize = partscale * plen;
        } else {
            partsize = partscale; // no length scaling
        }
        //if(type&(PART_DROP|PART_BUBBLE)) 
        partsize *= p->getSize()/5.0;

        ntlVec3Gfx pstart( mPartHeadDist *partsize, 0.0, 0.0 );
        ntlVec3Gfx pend  ( mPartTailDist *partsize, 0.0, 0.0 );
        gfxReal phi = 0.0;
        gfxReal phiD = 2.0*M_PI / (gfxReal)segments;

        ntlMat4Gfx cvmat; 
        cvmat.initId();
        pdir *= -1.0;
        ntlVec3Gfx cv1 = pdir;
        ntlVec3Gfx cv2 = ntlVec3Gfx(pdir[1], -pdir[0], 0.0);
        ntlVec3Gfx cv3 = cross( cv1, cv2);
        //? for(int l=0; l<3; l++) { cvmat.value[l][0] = cv1[l]; cvmat.value[l][1] = cv2[l]; cvmat.value[l][2] = cv3[l]; }
        pstart = (cvmat * pstart);
        pend = (cvmat * pend);

        for(int s=0; s<segments; s++) {
            ntlVec3Gfx p1( 0.0 );
            ntlVec3Gfx p2( 0.0 );

            gfxReal radscale = partNormSize;
            radscale = (partsize+partNormSize)*0.5;
            p1[1] += cos(phi) * radscale;
            p1[2] += sin(phi) * radscale;
            p2[1] += cos(phi + phiD) * radscale;
            p2[2] += sin(phi + phiD) * radscale;
            ntlVec3Gfx n1 = ntlVec3Gfx( 0.0, cos(phi), sin(phi) );
            ntlVec3Gfx n2 = ntlVec3Gfx( 0.0, cos(phi + phiD), sin(phi + phiD) );
            ntlVec3Gfx ns = n1*0.5 + n2*0.5;

            p1 = (cvmat * p1);
            p2 = (cvmat * p2);

            sceneAddTriangle( pos+pstart, pos+p1, pos+p2, 
                    ns,n1,n2, ntlVec3Gfx(0.0), 1, triangles,vertices,normals ); 
            sceneAddTriangle( pos+pend  , pos+p2, pos+p1, 
                    ns,n2,n1, ntlVec3Gfx(0.0), 1, triangles,vertices,normals ); 

            phi += phiD;
            tris += 2;
        }
    }

    } // t

    debMsgStd("ParticleTracer::getTriangles",DM_MSG,"Dumped "<<pcnt<<"/"<<mParts.size()<<" parts, tris:"<<tris<<", showonly:"<<mShowOnly,10);
    return; // DEBUG

#endif // ELBEEM_PLUGIN
}