doxygen: add newline after \file
[blender.git] / intern / elbeem / intern / particletracer.cpp
1 /** \file
2  * \ingroup elbeem
3  */
4 /******************************************************************************
5  *
6  * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
7  * Copyright 2003-2006 Nils Thuerey
8  *
9  * Particle Viewer/Tracer
10  *
11  *****************************************************************************/
12
13 #include <stdio.h>
14 //#include "../libs/my_gl.h"
15 //#include "../libs/my_glu.h"
16
17 /* own lib's */
18 #include "particletracer.h"
19 #include "ntl_matrices.h"
20 #include "ntl_ray.h"
21 #include "ntl_matrices.h"
22 #include "globals.h"
23
24 #include <zlib.h>
25
26
27 // particle object id counter
28 int ParticleObjectIdCnt = 1;
29
30 /******************************************************************************
31  * Standard constructor
32  *****************************************************************************/
33 ParticleTracer::ParticleTracer() :
34         ntlGeometryObject(),
35         mParts(),
36         //mTrailLength(1), mTrailInterval(1),mTrailIntervalCounter(0),
37         mPartSize(0.01),
38         mStart(-1.0), mEnd(1.0),
39         mSimStart(-1.0), mSimEnd(1.0),
40         mPartScale(0.1) , mPartHeadDist( 0.1 ), mPartTailDist( -0.1 ), mPartSegments( 4 ),
41         mValueScale(0),
42         mValueCutoffTop(0.0), mValueCutoffBottom(0.0),
43         mDumpParts(0), //mDumpText(0), 
44         mDumpTextFile(""), 
45         mDumpTextInterval(0.), mDumpTextLastTime(0.), mDumpTextCount(0),
46         mShowOnly(0), 
47         mNumInitialParts(0), mpTrafo(NULL),
48         mInitStart(-1.), mInitEnd(-1.),
49         mPrevs(), mTrailTimeLast(0.), mTrailInterval(-1.), mTrailLength(0)
50 {
51         debMsgStd("ParticleTracer::ParticleTracer",DM_MSG,"inited",10);
52 };
53
54 ParticleTracer::~ParticleTracer() {
55         debMsgStd("ParticleTracer::~ParticleTracer",DM_MSG,"destroyed",10);
56         if(mpTrafo) delete mpTrafo;
57 }
58
59 /*****************************************************************************/
60 //! parse settings from attributes (dont use own list!)
61 /*****************************************************************************/
62 void ParticleTracer::parseAttrList(AttributeList *att) 
63 {
64         AttributeList *tempAtt = mpAttrs; 
65         mpAttrs = att;
66
67         mNumInitialParts = mpAttrs->readInt("particles",mNumInitialParts, "ParticleTracer","mNumInitialParts", false);
68         //errMsg(" NUMP"," "<<mNumInitialParts);
69         mPartScale    = mpAttrs->readFloat("part_scale",mPartScale, "ParticleTracer","mPartScale", false);
70         mPartHeadDist = mpAttrs->readFloat("part_headdist",mPartHeadDist, "ParticleTracer","mPartHeadDist", false);
71         mPartTailDist = mpAttrs->readFloat("part_taildist",mPartTailDist, "ParticleTracer","mPartTailDist", false);
72         mPartSegments = mpAttrs->readInt  ("part_segments",mPartSegments, "ParticleTracer","mPartSegments", false);
73         mValueScale   = mpAttrs->readInt  ("part_valscale",mValueScale, "ParticleTracer","mValueScale", false);
74         mValueCutoffTop = mpAttrs->readFloat("part_valcutofftop",mValueCutoffTop, "ParticleTracer","mValueCutoffTop", false);
75         mValueCutoffBottom = mpAttrs->readFloat("part_valcutoffbottom",mValueCutoffBottom, "ParticleTracer","mValueCutoffBottom", false);
76
77         mDumpParts   = mpAttrs->readInt  ("part_dump",mDumpParts, "ParticleTracer","mDumpParts", false);
78         // mDumpText deprecatd, use mDumpTextInterval>0. instead
79         mShowOnly    = mpAttrs->readInt  ("part_showonly",mShowOnly, "ParticleTracer","mShowOnly", false);
80         mDumpTextFile= mpAttrs->readString("part_textdumpfile",mDumpTextFile, "ParticleTracer","mDumpTextFile", false);
81         mDumpTextInterval= mpAttrs->readFloat("part_textdumpinterval",mDumpTextInterval, "ParticleTracer","mDumpTextInterval", false);
82
83         string matPart;
84         matPart = mpAttrs->readString("material_part", "default", "ParticleTracer","material", false);
85         setMaterialName( matPart );
86
87         mInitStart = mpAttrs->readFloat("part_initstart",mInitStart, "ParticleTracer","mInitStart", false);
88         mInitEnd   = mpAttrs->readFloat("part_initend",  mInitEnd, "ParticleTracer","mInitEnd", false);
89
90         // unused...
91         //int mTrailLength  = 0; // UNUSED
92         //int mTrailInterval= 0; // UNUSED
93         mTrailLength  = mpAttrs->readInt("traillength",mTrailLength, "ParticleTracer","mTrailLength", false);
94         mTrailInterval= mpAttrs->readFloat("trailinterval",mTrailInterval, "ParticleTracer","mTrailInterval", false);
95
96         // restore old list
97         mpAttrs = tempAtt;
98 }
99
100 /******************************************************************************
101  * draw the particle array
102  *****************************************************************************/
103 void ParticleTracer::draw()
104 {
105 }
106
107 /******************************************************************************
108  * init trafo matrix
109  *****************************************************************************/
110 void ParticleTracer::initTrafoMatrix() {
111         ntlVec3Gfx scale = ntlVec3Gfx(
112                         (mEnd[0]-mStart[0])/(mSimEnd[0]-mSimStart[0]),
113                         (mEnd[1]-mStart[1])/(mSimEnd[1]-mSimStart[1]),
114                         (mEnd[2]-mStart[2])/(mSimEnd[2]-mSimStart[2])
115                         );
116         ntlVec3Gfx trans = mStart;
117         if(!mpTrafo) mpTrafo = new ntlMat4Gfx(0.0);
118         mpTrafo->initId();
119         for(int i=0; i<3; i++) { mpTrafo->value[i][i] = scale[i]; }
120         for(int i=0; i<3; i++) { mpTrafo->value[i][3] = trans[i]; }
121 }
122
123 /******************************************************************************
124  * adapt time step by rescaling velocities
125  *****************************************************************************/
126 void ParticleTracer::adaptPartTimestep(float factor) {
127         for(size_t i=0; i<mParts.size(); i++) {
128                 mParts[i].setVel( mParts[i].getVel() * factor );
129         }
130
131
132
133 /******************************************************************************
134  * add a particle at this position
135  *****************************************************************************/
136 void ParticleTracer::addParticle(float x, float y, float z) {
137         ntlVec3Gfx p(x,y,z);
138         ParticleObject part( p );
139         mParts.push_back( part );
140 }
141 void ParticleTracer::addFullParticle(ParticleObject &np) {
142         mParts.push_back( np );
143 }
144
145
146 void ParticleTracer::cleanup() {
147         // cleanup
148         int last = (int)mParts.size()-1;
149         if(mDumpTextInterval>0.) { errMsg("ParticleTracer::cleanup","Skipping cleanup due to text dump..."); return; }
150
151         for(int i=0; i<=last; i++) {
152                 if( mParts[i].getActive()==false ) {
153                         ParticleObject *p = &mParts[i];
154                         ParticleObject *p2 = &mParts[last];
155                         *p = *p2; last--; mParts.pop_back();
156                 }
157         }
158 }
159                 
160 /******************************************************************************
161  *! dump particles if desired 
162  *****************************************************************************/
163 void ParticleTracer::notifyOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename, double simtime) {
164         debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"obj:"<<this->getName()<<" frame:"<<frameNrStr<<" dumpp"<<mDumpParts<<" t"<<simtime, 10); // DEBUG
165
166         if(
167                         (dumptype==DUMP_FULLGEOMETRY)&&
168                         (mDumpParts>0)) {
169                 // dump to binary file
170                 std::ostringstream boutfilename("");
171                 boutfilename << outfilename <<"_particles_" << frameNrStr;
172                 if(glob_mpactive) {
173                         if(glob_mpindex>0) { boutfilename << "mp"<<glob_mpindex; }
174                 }
175                 boutfilename << ".gz";
176                 debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"B-Dumping: "<< this->getName() <<", particles:"<<mParts.size()<<" "<< " to "<<boutfilename.str()<<" #"<<frameNr , 7);
177                 //debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"B-Dumping: partgeodeb sim:"<<mSimStart<<","<<mSimEnd<<" geosize:"<<mStart<<","<<mEnd,2 );
178
179                 // output to zipped file
180                 gzFile gzf;
181                 gzf = gzopen(boutfilename.str().c_str(), "wb1");
182                 if(gzf) {
183                         int numParts;
184                         if(sizeof(numParts)!=4) { errMsg("ParticleTracer::notifyOfDump","Invalid int size"); return; }
185                         // only dump active particles
186                         numParts = 0;
187                         for(size_t i=0; i<mParts.size(); i++) {
188                                 if(!mParts[i].getActive()) continue;
189                                 numParts++;
190                         }
191                         gzwrite(gzf, &numParts, sizeof(numParts));
192                         for(size_t i=0; i<mParts.size(); i++) {
193                                 if(!mParts[i].getActive()) { continue; }
194                                 ParticleObject *p = &mParts[i];
195                                 //int type = p->getType();  // export whole type info
196                                 int type = p->getFlags(); // debug export whole type & status info
197                                 ntlVec3Gfx pos = p->getPos();
198                                 float size = p->getSize();
199
200                                 if(type&PART_FLOAT) { // WARNING same handling for dump!
201                                         // add one gridcell offset
202                                         //pos[2] += 1.0; 
203                                 } 
204                                 // display as drop for now externally
205                                 //else if(type&PART_TRACER) { type |= PART_DROP; }
206
207                                 pos = (*mpTrafo) * pos;
208
209                                 ntlVec3Gfx v = p->getVel();
210                                 v[0] *= mpTrafo->value[0][0];
211                                 v[1] *= mpTrafo->value[1][1];
212                                 v[2] *= mpTrafo->value[2][2];
213                                 // FIXME check: pos = (*mpTrafo) * pos;
214                                 gzwrite(gzf, &type, sizeof(type)); 
215                                 gzwrite(gzf, &size, sizeof(size)); 
216                                 for(int j=0; j<3; j++) { gzwrite(gzf, &pos[j], sizeof(float)); }
217                                 for(int j=0; j<3; j++) { gzwrite(gzf, &v[j], sizeof(float)); }
218                         }
219                         gzclose( gzf );
220                 }
221         } // dump?
222 }
223
224 void ParticleTracer::checkDumpTextPositions(double simtime) {
225         // dfor partial & full dump
226         if(mDumpTextInterval>0.) {
227                 debMsgStd("ParticleTracer::checkDumpTextPositions",DM_MSG,"t="<<simtime<<" last:"<<mDumpTextLastTime<<" inter:"<<mDumpTextInterval,7);
228         }
229
230         if((mDumpTextInterval>0.) && (simtime>mDumpTextLastTime+mDumpTextInterval)) {
231                 // dump to binary file
232                 std::ostringstream boutfilename("");
233                 if(mDumpTextFile.length()>1) {   
234                         boutfilename << mDumpTextFile <<  ".cpart2"; 
235                 } else {                           
236                         boutfilename << "_particles" <<  ".cpart2"; 
237                 }
238                 debMsgStd("ParticleTracer::checkDumpTextPositions",DM_MSG,"T-Dumping: "<< this->getName() <<", particles:"<<mParts.size()<<" "<< " to "<<boutfilename.str()<<" " , 7);
239
240                 int numParts = 0;
241                 // only dump bubble particles
242                 for(size_t i=0; i<mParts.size(); i++) {
243                         //if(!mParts[i].getActive()) continue;
244                         //if(!(mParts[i].getType()&PART_BUBBLE)) continue;
245                         numParts++;
246                 }
247
248                 // output to text file
249                 //gzFile gzf;
250                 FILE *stf;
251                 if(mDumpTextCount==0) {
252                         //gzf = gzopen(boutfilename.str().c_str(), "w0");
253                         stf = fopen(boutfilename.str().c_str(), "w");
254
255                         fprintf( stf, "\n\n# cparts generated by elbeem \n# no. of parts \nN %d \n\n",numParts);
256                         // fixed time scale for now
257                         fprintf( stf, "T %f \n\n", 1.0);
258                 } else {
259                         //gzf = gzopen(boutfilename.str().c_str(), "a+0");
260                         stf = fopen(boutfilename.str().c_str(), "a+");
261                 }
262
263                 // add current set
264                 if(stf) {
265                         fprintf( stf, "\n\n# new set at frame %d,t%f,p%d --------------------------------- \n\n", mDumpTextCount, simtime, numParts );
266                         fprintf( stf, "S %f \n\n", simtime );
267                         
268                         for(size_t i=0; i<mParts.size(); i++) {
269                                 ParticleObject *p = &mParts[i];
270                                 ntlVec3Gfx pos = p->getPos();
271                                 float size = p->getSize();
272                                 float infl = 1.;
273                                 //if(!mParts[i].getActive()) { size=0.; } // switch "off"
274                                 if(!mParts[i].getActive()) { infl=0.; } // switch "off"
275                                 if(!mParts[i].getInFluid()) { infl=0.; } // switch "off"
276                                 if(mParts[i].getLifeTime()<0.) { infl=0.; } // not yet active...
277
278                                 pos = (*mpTrafo) * pos;
279                                 ntlVec3Gfx v = p->getVel();
280                                 v[0] *= mpTrafo->value[0][0];
281                                 v[1] *= mpTrafo->value[1][1];
282                                 v[2] *= mpTrafo->value[2][2];
283                                 
284                                 fprintf( stf, "P %f %f %f \n", pos[0],pos[1],pos[2] );
285                                 if(size!=1.0) fprintf( stf, "s %f \n", size );
286                                 if(infl!=1.0) fprintf( stf, "i %f \n", infl );
287                                 fprintf( stf, "\n" );
288                         }
289
290                         fprintf( stf, "# %d end  ", mDumpTextCount );
291                         //gzclose( gzf );
292                         fclose( stf );
293
294                         mDumpTextCount++;
295                 }
296
297                 mDumpTextLastTime += mDumpTextInterval;
298         }
299
300 }
301
302
303 void ParticleTracer::checkTrails(double time) {
304         if(mTrailLength<1) return;
305         if(time-mTrailTimeLast > mTrailInterval) {
306
307                 if( (int)mPrevs.size() < mTrailLength) mPrevs.resize( mTrailLength );
308                 for(int i=mPrevs.size()-1; i>0; i--) {
309                         mPrevs[i] = mPrevs[i-1];
310                         //errMsg("TRAIL"," from "<<i<<" to "<<(i-1) );
311                 }
312                 mPrevs[0] = mParts;
313
314                 mTrailTimeLast += mTrailInterval;
315         }
316 }
317
318
319 /******************************************************************************
320  * Get triangles for rendering
321  *****************************************************************************/
322 void ParticleTracer::getTriangles(double time, vector<ntlTriangle> *triangles, 
323                                                                                                          vector<ntlVec3Gfx> *vertices, 
324                                                                                                          vector<ntlVec3Gfx> *normals, int objectId )
325 {
326 #ifdef ELBEEM_PLUGIN
327         // suppress warnings...
328         vertices = NULL; triangles = NULL;
329         normals = NULL; objectId = 0;
330         time = 0.;
331 #else // ELBEEM_PLUGIN
332         int pcnt = 0;
333         // currently not used in blender
334         objectId = 0; // remove, deprecated
335         if(mDumpParts>1) { 
336                 return; // only dump, no tri-gen
337         }
338
339         const bool debugParts = false;
340         int tris = 0;
341         int segments = mPartSegments;
342         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]));
343         ntlVec3Gfx trans = mStart;
344         time = 0.; // doesnt matter
345
346         for(size_t t=0; t<mPrevs.size()+1; t++) {
347                 vector<ParticleObject> *dparts;
348                 if(t==0) {
349                         dparts = &mParts;
350                 } else {
351                         dparts = &mPrevs[t-1];
352                 }
353                 //errMsg("TRAILT","prevs"<<t<<"/"<<mPrevs.size()<<" parts:"<<dparts->size() );
354
355         gfxReal partscale = mPartScale;
356         if(t>1) { 
357                 partscale *= (gfxReal)(mPrevs.size()+1-t) / (gfxReal)(mPrevs.size()+1); 
358         }
359         gfxReal partNormSize = 0.01 * partscale;
360         //for(size_t i=0; i<mParts.size(); i++) {
361         for(size_t i=0; i<dparts->size(); i++) {
362                 ParticleObject *p = &( (*dparts)[i] ); //  mParts[i];
363
364                 if(mShowOnly!=10) {
365                         // 10=show only deleted
366                         if( p->getActive()==false ) continue;
367                 } else {
368                         if( p->getActive()==true ) continue;
369                 }
370                 int type = p->getType();
371                 if(mShowOnly>0) {
372                         switch(mShowOnly) {
373                         case 1: if(!(type&PART_BUBBLE)) continue; break;
374                         case 2: if(!(type&PART_DROP))   continue; break;
375                         case 3: if(!(type&PART_INTER))  continue; break;
376                         case 4: if(!(type&PART_FLOAT))  continue; break;
377                         case 5: if(!(type&PART_TRACER))  continue; break;
378                         }
379                 } else {
380                         // by default dont display inter
381                         if(type&PART_INTER) continue;
382                 }
383
384                 pcnt++;
385                 ntlVec3Gfx pnew = p->getPos();
386                 if(type&PART_FLOAT) { // WARNING same handling for dump!
387                         if(p->getStatus()&PART_IN) { pnew[2] += 0.8; } // offset for display
388                         // add one gridcell offset
389                         //pnew[2] += 1.0; 
390                 }
391 #if LBMDIM==2
392                 pnew[2] += 0.001; // DEBUG
393                 pnew[2] += 0.009; // DEBUG
394 #endif 
395
396                 ntlVec3Gfx pdir = p->getVel();
397                 gfxReal plen = normalize( pdir );
398                 if( plen < 1e-05) pdir = ntlVec3Gfx(-1.0 ,0.0 ,0.0);
399                 ntlVec3Gfx pos = (*mpTrafo) * pnew;
400                 gfxReal partsize = 0.0;
401                 if(debugParts) errMsg("DebugParts"," i"<<i<<" new"<<pnew<<" vel"<<pdir<<"   pos="<<pos );
402                 //if(i==0 &&(debugParts)) errMsg("DebugParts"," i"<<i<<" new"<<pnew[0]<<" pos="<<pos[0]<<" scale="<<scale[0]<<"  t="<<trans[0] );
403                 
404                 // value length scaling?
405                 if(mValueScale==1) {
406                         partsize = partscale * plen;
407                 } else if(mValueScale==2) {
408                         // cut off scaling
409                         if(plen > mValueCutoffTop) continue;
410                         if(plen < mValueCutoffBottom) continue;
411                         partsize = partscale * plen;
412                 } else {
413                         partsize = partscale; // no length scaling
414                 }
415                 //if(type&(PART_DROP|PART_BUBBLE)) 
416                 partsize *= p->getSize()/5.0;
417
418                 ntlVec3Gfx pstart( mPartHeadDist *partsize, 0.0, 0.0 );
419                 ntlVec3Gfx pend  ( mPartTailDist *partsize, 0.0, 0.0 );
420                 gfxReal phi = 0.0;
421                 gfxReal phiD = 2.0*M_PI / (gfxReal)segments;
422
423                 ntlMat4Gfx cvmat; 
424                 cvmat.initId();
425                 pdir *= -1.0;
426                 ntlVec3Gfx cv1 = pdir;
427                 ntlVec3Gfx cv2 = ntlVec3Gfx(pdir[1], -pdir[0], 0.0);
428                 ntlVec3Gfx cv3 = cross( cv1, cv2);
429                 //? 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]; }
430                 pstart = (cvmat * pstart);
431                 pend = (cvmat * pend);
432
433                 for(int s=0; s<segments; s++) {
434                         ntlVec3Gfx p1( 0.0 );
435                         ntlVec3Gfx p2( 0.0 );
436
437                         gfxReal radscale = partNormSize;
438                         radscale = (partsize+partNormSize)*0.5;
439                         p1[1] += cos(phi) * radscale;
440                         p1[2] += sin(phi) * radscale;
441                         p2[1] += cos(phi + phiD) * radscale;
442                         p2[2] += sin(phi + phiD) * radscale;
443                         ntlVec3Gfx n1 = ntlVec3Gfx( 0.0, cos(phi), sin(phi) );
444                         ntlVec3Gfx n2 = ntlVec3Gfx( 0.0, cos(phi + phiD), sin(phi + phiD) );
445                         ntlVec3Gfx ns = n1*0.5 + n2*0.5;
446
447                         p1 = (cvmat * p1);
448                         p2 = (cvmat * p2);
449
450                         sceneAddTriangle( pos+pstart, pos+p1, pos+p2, 
451                                         ns,n1,n2, ntlVec3Gfx(0.0), 1, triangles,vertices,normals ); 
452                         sceneAddTriangle( pos+pend  , pos+p2, pos+p1, 
453                                         ns,n2,n1, ntlVec3Gfx(0.0), 1, triangles,vertices,normals ); 
454
455                         phi += phiD;
456                         tris += 2;
457                 }
458         }
459
460         } // t
461
462         debMsgStd("ParticleTracer::getTriangles",DM_MSG,"Dumped "<<pcnt<<"/"<<mParts.size()<<" parts, tris:"<<tris<<", showonly:"<<mShowOnly,10);
463         return; // DEBUG
464
465 #endif // ELBEEM_PLUGIN
466 }
467
468
469
470