828674675e236eb9d45b31be5c626afb50ed6bd8
[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(1.0) , mPartHeadDist( 0.5 ), mPartTailDist( -4.5 ), mPartSegments( 4 ),
37         mValueScale(0),
38         mValueCutoffTop(0.0), mValueCutoffBottom(0.0),
39         mDumpParts(0), mDumpText(0), mDumpTextFile(""), mShowOnly(0), 
40         mNumInitialParts(0), mpTrafo(NULL)
41 {
42         // check env var
43 //#ifdef ELBEEM_PLUGIN
44         //mDumpParts=1; // default on
45 //#endif // ELBEEM_PLUGIN
46         //if(getenv("ELBEEM_DUMPPARTICLE")) { // DEBUG!
47                 //int set = atoi( getenv("ELBEEM_DUMPPARTICLE") );
48                 //if((set>=0)&&(set!=mDumpParts)) {
49                         //mDumpParts=set;
50                         //debMsgStd("ParticleTrace",DM_NOTIFY,"Using envvar ELBEEM_DUMPPARTICLE to set mDumpParts to "<<set<<","<<mDumpParts,8);
51                 //}
52         //}
53         debMsgStd("ParticleTracer::ParticleTracer",DM_MSG,"inited",10);
54 };
55
56 ParticleTracer::~ParticleTracer() {
57         debMsgStd("ParticleTracer::~ParticleTracer",DM_MSG,"destroyed",10);
58 }
59
60 /*****************************************************************************/
61 //! parse settings from attributes (dont use own list!)
62 /*****************************************************************************/
63 void ParticleTracer::parseAttrList(AttributeList *att) 
64 {
65         AttributeList *tempAtt = mpAttrs; 
66         mpAttrs = att;
67
68         mNumInitialParts = mpAttrs->readInt("particles",mNumInitialParts, "ParticleTracer","mNumInitialParts", false);
69         errMsg(" NUMP"," "<<mNumInitialParts);
70         mPartScale    = mpAttrs->readFloat("part_scale",mPartScale, "ParticleTracer","mPartScale", false);
71         mPartHeadDist = mpAttrs->readFloat("part_headdist",mPartHeadDist, "ParticleTracer","mPartHeadDist", false);
72         mPartTailDist = mpAttrs->readFloat("part_taildist",mPartTailDist, "ParticleTracer","mPartTailDist", false);
73         mPartSegments = mpAttrs->readInt  ("part_segments",mPartSegments, "ParticleTracer","mPartSegments", false);
74         mValueScale   = mpAttrs->readInt  ("part_valscale",mValueScale, "ParticleTracer","mValueScale", false);
75         mValueCutoffTop = mpAttrs->readFloat("part_valcutofftop",mValueCutoffTop, "ParticleTracer","mValueCutoffTop", false);
76         mValueCutoffBottom = mpAttrs->readFloat("part_valcutoffbottom",mValueCutoffBottom, "ParticleTracer","mValueCutoffBottom", false);
77
78         mDumpParts   = mpAttrs->readInt  ("part_dump",mDumpParts, "ParticleTracer","mDumpParts", false);
79         mDumpText    = mpAttrs->readInt  ("part_textdump",mDumpText, "ParticleTracer","mDumpText", false);
80         mShowOnly    = mpAttrs->readInt  ("part_showonly",mShowOnly, "ParticleTracer","mShowOnly", false);
81         mDumpTextFile= mpAttrs->readString("part_textdumpfile",mDumpTextFile, "ParticleTracer","mDumpTextFile", false);
82
83         string matPart;
84         matPart = mpAttrs->readString("material_part", "default", "ParticleTracer","material", false);
85         setMaterialName( matPart );
86
87         // unused...
88         int mTrailLength  = 0; // UNUSED
89         int mTrailInterval= 0; // UNUSED
90         mTrailLength  = mpAttrs->readInt("traillength",mTrailLength, "ParticleTracer","mTrailLength", false);
91         mTrailInterval= mpAttrs->readInt("trailinterval",mTrailInterval, "ParticleTracer","mTrailInterval", false);
92
93         // restore old list
94         mpAttrs = tempAtt;
95 }
96
97 /******************************************************************************
98  * draw the particle array
99  *****************************************************************************/
100 void ParticleTracer::draw()
101 {
102 }
103
104 /******************************************************************************
105  * init trafo matrix
106  *****************************************************************************/
107 void ParticleTracer::initTrafoMatrix() {
108         ntlVec3Gfx scale = ntlVec3Gfx(
109                         (mEnd[0]-mStart[0])/(mSimEnd[0]-mSimStart[0]),
110                         (mEnd[1]-mStart[1])/(mSimEnd[1]-mSimStart[1]),
111                         (mEnd[2]-mStart[2])/(mSimEnd[2]-mSimStart[2])
112                         );
113         ntlVec3Gfx trans = mStart;
114         if(!mpTrafo) mpTrafo = new ntlMat4Gfx(0.0);
115         mpTrafo->initId();
116         for(int i=0; i<3; i++) { mpTrafo->value[i][i] = scale[i]; }
117         for(int i=0; i<3; i++) { mpTrafo->value[i][3] = trans[i]; }
118 }
119
120 /******************************************************************************
121  * adapt time step by rescaling velocities
122  *****************************************************************************/
123 void ParticleTracer::adaptPartTimestep(float factor) {
124         for(size_t i=0; i<mParts.size(); i++) {
125                 mParts[i].setVel( mParts[i].getVel() * factor );
126         }
127
128
129
130 /******************************************************************************
131  * add a particle at this position
132  *****************************************************************************/
133 void ParticleTracer::addParticle(float x, float y, float z)
134 {
135         ntlVec3Gfx p(x,y,z);
136         ParticleObject part( p );
137         //mParts.push_back( part );
138         // TODO handle other arrays?
139         //part.setActive( false );
140         mParts.push_back( part );
141         //for(size_t l=0; l<mParts.size(); l++) {
142                 // add deactivated particles to other arrays
143                 // deactivate further particles
144                 //if(l>1) {
145                         //mParts[l][ mParts.size()-1 ].setActive( false );
146                 //}
147         //}
148 }
149
150
151 void ParticleTracer::cleanup() {
152         // cleanup
153         int last = (int)mParts.size()-1;
154         if(mDumpText>0) { errMsg("ParticleTracer::cleanup","Skipping cleanup due to text dump..."); return; }
155
156         for(int i=0; i<=last; i++) {
157                 if( mParts[i].getActive()==false ) {
158                         ParticleObject *p = &mParts[i];
159                         ParticleObject *p2 = &mParts[last];
160                         *p = *p2; last--; mParts.pop_back();
161                 }
162         }
163 }
164                 
165 /******************************************************************************
166  *! dump particles if desired 
167  *****************************************************************************/
168 void ParticleTracer::notifyOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename, double simtime) {
169         debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"obj:"<<this->getName()<<" frame:"<<frameNrStr, 10); // DEBUG
170
171         if(
172                         (dumptype==DUMP_FULLGEOMETRY)&&
173                         (mDumpParts>0)) {
174                 // dump to binary file
175                 std::ostringstream boutfilename("");
176                 boutfilename << outfilename <<"_particles_" << frameNrStr<< ".gz";
177                 debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"B-Dumping: "<< this->getName() <<", particles:"<<mParts.size()<<" "<< " to "<<boutfilename.str()<<" #"<<frameNr , 7);
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                                 //if(mParts[i].getLifeTime()<30) { continue; } //? CHECK
190                                 numParts++;
191                         }
192                         gzwrite(gzf, &numParts, sizeof(numParts));
193                         for(size_t i=0; i<mParts.size(); i++) {
194                                 if(!mParts[i].getActive()) { continue; }
195                                 //if(mParts[i].getLifeTime()<30) { continue; } //? CHECK
196                                 ParticleObject *p = &mParts[i];
197                                 //int type = p->getType();  // export whole type info
198                                 int type = p->getFlags(); // debug export whole type & status info
199                                 ntlVec3Gfx pos = p->getPos();
200                                 float size = p->getSize();
201
202                                 if(type&PART_FLOAT) { // WARNING same handling for dump!
203                                         // add one gridcell offset
204                                         //pos[2] += 1.0; 
205                                 } 
206                                 // display as drop for now externally
207                                 //else if(type&PART_TRACER) { type |= PART_DROP; }
208
209                                 pos = (*mpTrafo) * pos;
210
211                                 ntlVec3Gfx v = p->getVel();
212                                 v[0] *= mpTrafo->value[0][0];
213                                 v[1] *= mpTrafo->value[1][1];
214                                 v[2] *= mpTrafo->value[2][2];
215                                 // FIXME check: pos = (*mpTrafo) * pos;
216                                 gzwrite(gzf, &type, sizeof(type)); 
217                                 gzwrite(gzf, &size, sizeof(size)); 
218                                 for(int j=0; j<3; j++) { gzwrite(gzf, &pos[j], sizeof(float)); }
219                                 for(int j=0; j<3; j++) { gzwrite(gzf, &v[j], sizeof(float)); }
220                         }
221                         gzclose( gzf );
222                 }
223         } // dump?
224
225         // dfor partial & full dump
226         if(mDumpText>0) {
227                 // dump to binary file
228                 std::ostringstream boutfilename("");
229                 if(mDumpTextFile.length()>1) {   
230                         boutfilename << mDumpTextFile <<  ".cpart2"; 
231                 } else {                           
232                         boutfilename << outfilename <<"_particles" <<  ".cpart2"; 
233                 }
234                 debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"T-Dumping: "<< this->getName() <<", particles:"<<mParts.size()<<" "<< " to "<<boutfilename.str()<<" #"<<frameNr , 7);
235
236                 int numParts = 0;
237                 // only dump bubble particles
238                 for(size_t i=0; i<mParts.size(); i++) {
239                         //if(!mParts[i].getActive()) continue;
240                         //if(!(mParts[i].getType()&PART_BUBBLE)) continue;
241                         numParts++;
242                 }
243
244                 // output to text file
245                 gzFile gzf;
246                 if(frameNr==0) {
247                         gzf = gzopen(boutfilename.str().c_str(), "w0");
248
249                         gzprintf( gzf, "\n\n# cparts generated by elbeem \n# no. of parts \nN %d \n\n",numParts);
250                         // fixed time scale for now
251                         gzprintf( gzf, "T %f \n\n", 1.0);
252                 } else {
253                         gzf = gzopen(boutfilename.str().c_str(), "a+0");
254                 }
255
256                 // add current set
257                 if(gzf) {
258                         gzprintf( gzf, "\n\n# new set at frame %d,t%f,p%d --------------------------------- \n\n", frameNr, simtime, numParts );
259                         gzprintf( gzf, "S %f \n\n", simtime );
260                         
261                         for(size_t i=0; i<mParts.size(); i++) {
262                                 ParticleObject *p = &mParts[i];
263                                 ntlVec3Gfx pos = p->getPos();
264                                 float size = p->getSize();
265                                 if(!mParts[i].getActive()) { size=0.; } // switch "off"
266
267                                 pos = (*mpTrafo) * pos;
268                                 ntlVec3Gfx v = p->getVel();
269                                 v[0] *= mpTrafo->value[0][0];
270                                 v[1] *= mpTrafo->value[1][1];
271                                 v[2] *= mpTrafo->value[2][2];
272                                 
273                                 gzprintf( gzf, "P %f %f %f \n", pos[0],pos[1],pos[2] );
274                                 gzprintf( gzf, "s %f \n", size );
275                                 gzprintf( gzf, "\n", size );
276                         }
277
278                         gzprintf( gzf, "# %d end  ", frameNr );
279                         gzclose( gzf );
280                 }
281         }
282 }
283
284
285
286 /******************************************************************************
287  * Get triangles for rendering
288  *****************************************************************************/
289 void ParticleTracer::getTriangles( vector<ntlTriangle> *triangles, 
290                                                                                                          vector<ntlVec3Gfx> *vertices, 
291                                                                                                          vector<ntlVec3Gfx> *normals, int objectId )
292 {
293 #ifdef ELBEEM_PLUGIN
294         // suppress warnings...
295         vertices = NULL; triangles = NULL;
296         normals = NULL; objectId = 0;
297 #else // ELBEEM_PLUGIN
298         // currently not used in blender
299         objectId = 0; // remove, deprecated
300         if(mDumpParts>1) { 
301                 return; // only dump, no tri-gen
302         }
303
304         const bool debugParts = false;
305         int tris = 0;
306         gfxReal partNormSize = 0.01 * mPartScale;
307         int segments = mPartSegments;
308
309         //int lnewst = mTrailLength-1;
310         // TODO get rid of mPart[X] array
311         //? int loldst = mTrailLength-2;
312         // trails gehen nicht so richtig mit der
313         // richtung der partikel...
314
315         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]));
316         ntlVec3Gfx trans = mStart;
317
318         for(size_t i=0; i<mParts.size(); i++) {
319                 //mParts[0][i].setActive(true);
320                 ParticleObject *p = &mParts[i];
321
322                 if(mShowOnly!=10) {
323                         // 10=show only deleted
324                         if( p->getActive()==false ) continue;
325                 } else {
326                         if( p->getActive()==true ) continue;
327                 }
328                 int type = p->getType();
329                 if(mShowOnly>0) {
330                         switch(mShowOnly) {
331                         case 1: if(!(type&PART_BUBBLE)) continue; break;
332                         case 2: if(!(type&PART_DROP))   continue; break;
333                         case 3: if(!(type&PART_INTER))  continue; break;
334                         case 4: if(!(type&PART_FLOAT))  continue; break;
335                         case 5: if(!(type&PART_TRACER))  continue; break;
336                         }
337                 } else {
338                         // by default dont display inter
339                         if(type&PART_INTER) continue;
340                 }
341
342                 ntlVec3Gfx pnew = p->getPos();
343                 if(type&PART_FLOAT) { // WARNING same handling for dump!
344                         // add one gridcell offset
345                         //pnew[2] += 1.0; 
346                 }
347
348                 ntlVec3Gfx pdir = p->getVel();
349                 gfxReal plen = normalize( pdir );
350                 if( plen < 1e-05) pdir = ntlVec3Gfx(-1.0 ,0.0 ,0.0);
351                 ntlVec3Gfx pos = (*mpTrafo) * pnew;
352                 //ntlVec3Gfx pos = pnew; // T
353                 //pos[0] = pos[0]*scale[0]+trans[0]; // T
354                 //pos[1] = pos[1]*scale[1]+trans[1]; // T
355                 //pos[2] = pos[2]*scale[2]+trans[2]; // T
356                 gfxReal partsize = 0.0;
357                 if(debugParts) errMsg("DebugParts"," i"<<i<<" new"<<pnew<<" vel"<<pdir<<"   pos="<<pos );
358                 //if(i==0 &&(debugParts)) errMsg("DebugParts"," i"<<i<<" new"<<pnew[0]<<" pos="<<pos[0]<<" scale="<<scale[0]<<"  t="<<trans[0] );
359                 
360                 // value length scaling?
361                 if(mValueScale==1) {
362                         partsize = mPartScale * plen;
363                 } else if(mValueScale==2) {
364                         // cut off scaling
365                         if(plen > mValueCutoffTop) continue;
366                         if(plen < mValueCutoffBottom) continue;
367                         partsize = mPartScale * plen;
368                 } else {
369                         partsize = mPartScale; // no length scaling
370                 }
371                 //if(type&(PART_DROP|PART_BUBBLE)) 
372                 partsize *= p->getSize()/5.0;
373
374                 ntlVec3Gfx pstart( mPartHeadDist *partsize, 0.0, 0.0 );
375                 ntlVec3Gfx pend  ( mPartTailDist *partsize, 0.0, 0.0 );
376                 gfxReal phi = 0.0;
377                 gfxReal phiD = 2.0*M_PI / (gfxReal)segments;
378
379                 ntlMat4Gfx cvmat; 
380                 cvmat.initId();
381                 pdir *= -1.0;
382                 ntlVec3Gfx cv1 = pdir;
383                 ntlVec3Gfx cv2 = ntlVec3Gfx(pdir[1], -pdir[0], 0.0);
384                 ntlVec3Gfx cv3 = cross( cv1, cv2);
385                 for(int l=0; l<3; l++) {
386                         //? cvmat.value[l][0] = cv1[l];
387                         //? cvmat.value[l][1] = cv2[l];
388                         //? cvmat.value[l][2] = cv3[l];
389                 }
390                 pstart = (cvmat * pstart);
391                 pend = (cvmat * pend);
392
393                 for(int s=0; s<segments; s++) {
394                         ntlVec3Gfx p1( 0.0 );
395                         ntlVec3Gfx p2( 0.0 );
396
397                         gfxReal radscale = partNormSize;
398                         radscale = (partsize+partNormSize)*0.5;
399                         p1[1] += cos(phi) * radscale;
400                         p1[2] += sin(phi) * radscale;
401                         p2[1] += cos(phi + phiD) * radscale;
402                         p2[2] += sin(phi + phiD) * radscale;
403                         ntlVec3Gfx n1 = ntlVec3Gfx( 0.0, cos(phi), sin(phi) );
404                         ntlVec3Gfx n2 = ntlVec3Gfx( 0.0, cos(phi + phiD), sin(phi + phiD) );
405                         ntlVec3Gfx ns = n1*0.5 + n2*0.5;
406
407                         p1 = (cvmat * p1);
408                         p2 = (cvmat * p2);
409
410                         sceneAddTriangle( pos+pstart, pos+p1, pos+p2, 
411                                         ns,n1,n2, ntlVec3Gfx(0.0), 1, triangles,vertices,normals ); 
412                         sceneAddTriangle( pos+pend  , pos+p2, pos+p1, 
413                                         ns,n2,n1, ntlVec3Gfx(0.0), 1, triangles,vertices,normals ); 
414
415                         phi += phiD;
416                         tris += 2;
417                 }
418         }
419
420         //} // trail
421         return; // DEBUG
422
423 #endif // ELBEEM_PLUGIN
424 }
425
426
427
428