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