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