2.5: warning fixes
[blender.git] / source / blender / blenkernel / intern / smoke.c
1 /**
2  * smoke.c
3  *
4  * $Id$
5  *
6  * ***** BEGIN GPL LICENSE BLOCK *****
7  *
8  * This program is free software; you can redistribute it and/or
9  * modify it under the terms of the GNU General Public License
10  * as published by the Free Software Foundation; either version 2
11  * of the License, or (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with this program; if not, write to the Free Software Foundation,
20  * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
21  *
22  * The Original Code is Copyright (C) Blender Foundation.
23  * All rights reserved.
24  *
25  * The Original Code is: all of this file.
26  *
27  * Contributor(s): Daniel Genrich
28  *
29  * ***** END GPL LICENSE BLOCK *****
30  */
31
32 /* Part of the code copied from elbeem fluid library, copyright by Nils Thuerey */
33
34 #include <GL/glew.h>
35
36 #include "MEM_guardedalloc.h"
37
38 #include <float.h>
39 #include <math.h>
40 #include "stdio.h"
41
42 #include "BLI_linklist.h"
43 #include "BLI_rand.h"
44 #include "BLI_jitter.h"
45 #include "BLI_blenlib.h"
46 #include "BLI_arithb.h"
47 #include "BLI_edgehash.h"
48 #include "BLI_kdtree.h"
49 #include "BLI_kdopbvh.h"
50
51 #include "BKE_bvhutils.h"
52 #include "BKE_cdderivedmesh.h"
53 #include "BKE_customdata.h"
54 #include "BKE_DerivedMesh.h"
55 #include "BKE_modifier.h"
56 #include "BKE_particle.h"
57 #include "BKE_utildefines.h"
58
59 #include "DNA_customdata_types.h"
60 #include "DNA_group_types.h"
61 #include "DNA_mesh_types.h"
62 #include "DNA_meshdata_types.h"
63 #include "DNA_modifier_types.h"
64 #include "DNA_object_types.h"
65 #include "DNA_particle_types.h"
66 #include "DNA_scene_types.h"
67 #include "DNA_smoke_types.h"
68
69 #include "smoke_API.h"
70
71 #include "BKE_smoke.h"
72
73 #ifdef _WIN32
74 #include <time.h>
75 #include <stdio.h>
76 #include <conio.h>
77 #include <windows.h>
78
79 static LARGE_INTEGER liFrequency;
80 static LARGE_INTEGER liStartTime;
81 static LARGE_INTEGER liCurrentTime;
82
83 static void tstart ( void )
84 {
85         QueryPerformanceFrequency ( &liFrequency );
86         QueryPerformanceCounter ( &liStartTime );
87 }
88 static void tend ( void )
89 {
90         QueryPerformanceCounter ( &liCurrentTime );
91 }
92 static double tval()
93 {
94         return ((double)( (liCurrentTime.QuadPart - liStartTime.QuadPart)* (double)1000.0/(double)liFrequency.QuadPart ));
95 }
96 #else
97 #include <sys/time.h>
98 static struct timeval _tstart, _tend;
99 static struct timezone tz;
100 static void tstart ( void )
101 {
102         gettimeofday ( &_tstart, &tz );
103 }
104 static void tend ( void )
105 {
106         gettimeofday ( &_tend,&tz );
107 }
108 static double tval()
109 {
110         double t1, t2;
111         t1 = ( double ) _tstart.tv_sec*1000 + ( double ) _tstart.tv_usec/ ( 1000 );
112         t2 = ( double ) _tend.tv_sec*1000 + ( double ) _tend.tv_usec/ ( 1000 );
113         return t2-t1;
114 }
115 #endif
116
117 struct Object;
118 struct Scene;
119 struct DerivedMesh;
120 struct SmokeModifierData;
121
122 // forward declerations
123 static void get_cell(struct SmokeModifierData *smd, float *pos, int *cell, int correct);
124 static void get_bigcell(struct SmokeModifierData *smd, float *pos, int *cell, int correct);
125 void calcTriangleDivs(Object *ob, MVert *verts, int numverts, MFace *tris, int numfaces, int numtris, int **tridivs, float cell_len);
126
127 #define TRI_UVOFFSET (1./4.)
128
129 int smokeModifier_init (SmokeModifierData *smd, Object *ob, Scene *scene, DerivedMesh *dm)
130 {
131         if((smd->type & MOD_SMOKE_TYPE_DOMAIN) && smd->domain && !smd->domain->fluid)
132         {
133                 size_t i;
134                 float min[3] = {FLT_MAX, FLT_MAX, FLT_MAX}, max[3] = {-FLT_MAX, -FLT_MAX, -FLT_MAX};
135                 float size[3];
136                 MVert *verts = dm->getVertArray(dm);
137                 float scale = 0.0;
138                 int res;                
139
140                 res = smd->domain->maxres;
141
142                 // get BB of domain
143                 for(i = 0; i < dm->getNumVerts(dm); i++)
144                 {
145                         float tmp[3];
146
147                         VECCOPY(tmp, verts[i].co);
148                         Mat4MulVecfl(ob->obmat, tmp);
149
150                         // min BB
151                         min[0] = MIN2(min[0], tmp[0]);
152                         min[1] = MIN2(min[1], tmp[1]);
153                         min[2] = MIN2(min[2], tmp[2]);
154
155                         // max BB
156                         max[0] = MAX2(max[0], tmp[0]);
157                         max[1] = MAX2(max[1], tmp[1]);
158                         max[2] = MAX2(max[2], tmp[2]);
159                 }
160
161                 VECCOPY(smd->domain->p0, min);
162                 VECCOPY(smd->domain->p1, max);
163
164                 // calc other res with max_res provided
165                 VECSUB(size, max, min);
166                 if(size[0] > size[1])
167                 {
168                         if(size[0] > size[1])
169                         {
170                                 scale = res / size[0];
171                                 smd->domain->dx = size[0] / res;
172                                 smd->domain->res[0] = res;
173                                 smd->domain->res[1] = (int)(size[1] * scale + 0.5);
174                                 smd->domain->res[2] = (int)(size[2] * scale + 0.5);
175                         }
176                         else
177                         {
178                                 scale = res / size[1];
179                                 smd->domain->dx = size[1] / res;
180                                 smd->domain->res[1] = res;
181                                 smd->domain->res[0] = (int)(size[0] * scale + 0.5);
182                                 smd->domain->res[2] = (int)(size[2] * scale + 0.5);
183                         }
184                 }
185                 else
186                 {
187                         if(size[1] > size[2])
188                         {
189                                 scale = res / size[1];
190                                 smd->domain->dx = size[1] / res;
191                                 smd->domain->res[1] = res;
192                                 smd->domain->res[0] = (int)(size[0] * scale + 0.5);
193                                 smd->domain->res[2] = (int)(size[2] * scale + 0.5);
194                         }
195                         else
196                         {
197                                 scale = res / size[2];
198                                 smd->domain->dx = size[2] / res;
199                                 smd->domain->res[2] = res;
200                                 smd->domain->res[0] = (int)(size[0] * scale + 0.5);
201                                 smd->domain->res[1] = (int)(size[1] * scale + 0.5);
202                         }
203                 }
204
205                 // TODO: put in failsafe if res<=0 - dg
206
207                 // printf("res[0]: %d, res[1]: %d, res[2]: %d\n", smd->domain->res[0], smd->domain->res[1], smd->domain->res[2]);
208
209                 // dt max is 0.1
210                 smd->domain->fluid = smoke_init(smd->domain->res, smd->domain->p0, 0.1);
211                 smd->domain->wt = smoke_turbulence_init(smd->domain->res,  (smd->domain->flags & MOD_SMOKE_HIGHRES) ? (smd->domain->amplify + 1) : 0, smd->domain->noise);
212                 smd->time = scene->r.cfra;
213                 smd->domain->firstframe = smd->time;
214                 
215                 smoke_initBlenderRNA(smd->domain->fluid, &(smd->domain->alpha), &(smd->domain->beta));
216
217                 if(smd->domain->wt)
218                         smoke_initWaveletBlenderRNA(smd->domain->wt, &(smd->domain->strength));
219
220                 return 1;
221         }
222         else if((smd->type & MOD_SMOKE_TYPE_FLOW) && smd->flow)
223         {
224                 // handle flow object here
225                 // XXX TODO
226
227                 smd->time = scene->r.cfra;
228
229                 // update particle lifetime to be one frame
230                 // smd->flow->psys->part->lifetime = scene->r.efra + 1;
231 /*
232                 if(!smd->flow->bvh)
233                 {
234                         // smd->flow->bvh = MEM_callocN(sizeof(BVHTreeFromMesh), "smoke_bvhfromfaces");
235                         // bvhtree_from_mesh_faces(smd->flow->bvh, dm, 0.0, 2, 6);
236
237                         // copy obmat
238                         // Mat4CpyMat4(smd->flow->mat, ob->obmat);
239                         // Mat4CpyMat4(smd->flow->mat_old, ob->obmat);
240                 }
241 */
242
243                 return 1;
244         }
245         else if((smd->type & MOD_SMOKE_TYPE_COLL))
246         {
247                 smd->time = scene->r.cfra;
248
249                 // todo: delete this when loading colls work -dg
250                 if(!smd->coll)
251                         smokeModifier_createType(smd);
252
253                 if(!smd->coll->points)
254                 {
255                         // init collision points
256                         SmokeCollSettings *scs = smd->coll;
257                         MVert *mvert = dm->getVertArray(dm);
258                         MFace *mface = dm->getFaceArray(dm);
259                         size_t i = 0, divs = 0;
260                         int *tridivs = NULL;
261                         float cell_len = 1.0 / 50.0; // for res = 50
262                         size_t newdivs = 0;
263                         //size_t max_points = 0;
264                         size_t quads = 0, facecounter = 0;
265
266                         // copy obmat
267                         Mat4CpyMat4(scs->mat, ob->obmat);
268                         Mat4CpyMat4(scs->mat_old, ob->obmat);
269
270                         // count quads
271                         for(i = 0; i < dm->getNumFaces(dm); i++)
272                         {
273                                 if(mface[i].v4)
274                                         quads++;
275                         }
276
277                         calcTriangleDivs(ob, mvert, dm->getNumVerts(dm), mface,  dm->getNumFaces(dm), dm->getNumFaces(dm) + quads, &tridivs, cell_len);
278
279                         // count triangle divisions
280                         for(i = 0; i < dm->getNumFaces(dm) + quads; i++)
281                         {
282                                 divs += (tridivs[3 * i] + 1) * (tridivs[3 * i + 1] + 1) * (tridivs[3 * i + 2] + 1);
283                         }
284
285                         // printf("divs: %d\n", divs);
286
287                         scs->points = MEM_callocN(sizeof(float) * (dm->getNumVerts(dm) + divs) * 3, "SmokeCollPoints");
288
289                         for(i = 0; i < dm->getNumVerts(dm); i++)
290                         {
291                                 float tmpvec[3];
292                                 VECCOPY(tmpvec, mvert[i].co);
293                                 Mat4MulVecfl (ob->obmat, tmpvec);
294                                 VECCOPY(&scs->points[i * 3], tmpvec);
295                         }
296                         
297                         for(i = 0, facecounter = 0; i < dm->getNumFaces(dm); i++)
298                         {
299                                 int again = 0;
300                                 do
301                                 {
302                                         size_t j, k;
303                                         int divs1 = tridivs[3 * facecounter + 0];
304                                         int divs2 = tridivs[3 * facecounter + 1];
305                                         //int divs3 = tridivs[3 * facecounter + 2];
306                                         float side1[3], side2[3], trinormorg[3], trinorm[3];
307                                         
308                                         if(again == 1 && mface[i].v4)
309                                         {
310                                                 VECSUB(side1,  mvert[ mface[i].v3 ].co, mvert[ mface[i].v1 ].co);
311                                                 VECSUB(side2,  mvert[ mface[i].v4 ].co, mvert[ mface[i].v1 ].co);
312                                         }
313                                         else
314                                         {
315                                                 VECSUB(side1,  mvert[ mface[i].v2 ].co, mvert[ mface[i].v1 ].co);
316                                                 VECSUB(side2,  mvert[ mface[i].v3 ].co, mvert[ mface[i].v1 ].co);
317                                         }
318
319                                         Crossf(trinormorg, side1, side2);
320                                         Normalize(trinormorg);
321                                         VECCOPY(trinorm, trinormorg);
322                                         VecMulf(trinorm, 0.25 * cell_len);
323
324                                         for(j = 0; j <= divs1; j++)
325                                         {
326                                                 for(k = 0; k <= divs2; k++)
327                                                 {
328                                                         float p1[3], p2[3], p3[3], p[3]={0,0,0}; 
329                                                         const float uf = (float)(j + TRI_UVOFFSET) / (float)(divs1 + 0.0);
330                                                         const float vf = (float)(k + TRI_UVOFFSET) / (float)(divs2 + 0.0);
331                                                         float tmpvec[3];
332                                                         
333                                                         if(uf+vf > 1.0) 
334                                                         {
335                                                                 // printf("bigger - divs1: %d, divs2: %d\n", divs1, divs2);
336                                                                 continue;
337                                                         }
338
339                                                         VECCOPY(p1, mvert[ mface[i].v1 ].co);
340                                                         if(again == 1 && mface[i].v4)
341                                                         {
342                                                                 VECCOPY(p2, mvert[ mface[i].v3 ].co);
343                                                                 VECCOPY(p3, mvert[ mface[i].v4 ].co);
344                                                         }
345                                                         else
346                                                         {
347                                                                 VECCOPY(p2, mvert[ mface[i].v2 ].co);
348                                                                 VECCOPY(p3, mvert[ mface[i].v3 ].co);
349                                                         }
350
351                                                         VecMulf(p1, (1.0-uf-vf));
352                                                         VecMulf(p2, uf);
353                                                         VecMulf(p3, vf);
354                                                         
355                                                         VECADD(p, p1, p2);
356                                                         VECADD(p, p, p3);
357
358                                                         if(newdivs > divs)
359                                                                 printf("mem problem\n");
360
361                                                         // mMovPoints.push_back(p + trinorm);
362                                                         VECCOPY(tmpvec, p);
363                                                         VECADD(tmpvec, tmpvec, trinorm);
364                                                         Mat4MulVecfl (ob->obmat, tmpvec);
365                                                         VECCOPY(&scs->points[3 * (dm->getNumVerts(dm) + newdivs)], tmpvec);
366                                                         newdivs++;
367
368                                                         if(newdivs > divs)
369                                                                 printf("mem problem\n");
370
371                                                         // mMovPoints.push_back(p - trinorm);
372                                                         VECCOPY(tmpvec, p);
373                                                         VECSUB(tmpvec, tmpvec, trinorm);
374                                                         Mat4MulVecfl (ob->obmat, tmpvec);
375                                                         VECCOPY(&scs->points[3 * (dm->getNumVerts(dm) + newdivs)], tmpvec);
376                                                         newdivs++;
377                                                 }
378                                         }
379
380                                         if(again == 0 && mface[i].v4)
381                                                 again++;
382                                         else
383                                                 again = 0;
384
385                                         facecounter++;
386
387                                 } while(again!=0);
388                         }
389
390                         scs->numpoints = dm->getNumVerts(dm) + newdivs;
391
392                         MEM_freeN(tridivs);
393                 }
394
395                 if(!smd->coll->bvhtree)
396                 {
397                         smd->coll->bvhtree = NULL; // bvhtree_build_from_smoke ( ob->obmat, dm->getFaceArray(dm), dm->getNumFaces(dm), dm->getVertArray(dm), dm->getNumVerts(dm), 0.0 );
398                 }
399
400         }
401
402         return 0;
403 }
404
405 /*! init triangle divisions */
406 void calcTriangleDivs(Object *ob, MVert *verts, int numverts, MFace *faces, int numfaces, int numtris, int **tridivs, float cell_len) 
407 {
408         // mTriangleDivs1.resize( faces.size() );
409         // mTriangleDivs2.resize( faces.size() );
410         // mTriangleDivs3.resize( faces.size() );
411
412         size_t i = 0, facecounter = 0;
413         float maxscale[3] = {1,1,1}; // = channelFindMaxVf(mcScale);
414         float maxpart = ABS(maxscale[0]);
415         float scaleFac = 0;
416         float fsTri = 0;
417         if(ABS(maxscale[1])>maxpart) maxpart = ABS(maxscale[1]);
418         if(ABS(maxscale[2])>maxpart) maxpart = ABS(maxscale[2]);
419         scaleFac = 1.0 / maxpart;
420         // featureSize = mLevel[mMaxRefine].nodeSize
421         fsTri = cell_len * 0.5 * scaleFac;
422
423         if(*tridivs)
424                 MEM_freeN(*tridivs);
425
426         *tridivs = MEM_callocN(sizeof(int) * numtris * 3, "Smoke_Tridivs");
427
428         for(i = 0, facecounter = 0; i < numfaces; i++) 
429         {
430                 float p0[3], p1[3], p2[3];
431                 float side1[3];
432                 float side2[3];
433                 float side3[3];
434                 int divs1=0, divs2=0, divs3=0;
435
436                 VECCOPY(p0, verts[faces[i].v1].co);
437                 Mat4MulVecfl (ob->obmat, p0);
438                 VECCOPY(p1, verts[faces[i].v2].co);
439                 Mat4MulVecfl (ob->obmat, p1);
440                 VECCOPY(p2, verts[faces[i].v3].co);
441                 Mat4MulVecfl (ob->obmat, p2);
442
443                 VECSUB(side1, p1, p0);
444                 VECSUB(side2, p2, p0);
445                 VECSUB(side3, p1, p2);
446
447                 if(INPR(side1, side1) > fsTri*fsTri) 
448                 { 
449                         float tmp = Normalize(side1);
450                         divs1 = (int)ceil(tmp/fsTri); 
451                 }
452                 if(INPR(side2, side2) > fsTri*fsTri) 
453                 { 
454                         float tmp = Normalize(side2);
455                         divs2 = (int)ceil(tmp/fsTri); 
456                         
457                         /*
458                         // debug
459                         if(i==0)
460                                 printf("b tmp: %f, fsTri: %f, divs2: %d\n", tmp, fsTri, divs2);
461                         */
462                 }
463
464                 (*tridivs)[3 * facecounter + 0] = divs1;
465                 (*tridivs)[3 * facecounter + 1] = divs2;
466                 (*tridivs)[3 * facecounter + 2] = divs3;
467
468                 // TODO quad case
469                 if(faces[i].v4)
470                 {
471                         divs1=0, divs2=0, divs3=0;
472
473                         facecounter++;
474                         
475                         VECCOPY(p0, verts[faces[i].v3].co);
476                         Mat4MulVecfl (ob->obmat, p0);
477                         VECCOPY(p1, verts[faces[i].v4].co);
478                         Mat4MulVecfl (ob->obmat, p1);
479                         VECCOPY(p2, verts[faces[i].v1].co);
480                         Mat4MulVecfl (ob->obmat, p2);
481
482                         VECSUB(side1, p1, p0);
483                         VECSUB(side2, p2, p0);
484                         VECSUB(side3, p1, p2);
485
486                         if(INPR(side1, side1) > fsTri*fsTri) 
487                         { 
488                                 float tmp = Normalize(side1);
489                                 divs1 = (int)ceil(tmp/fsTri); 
490                         }
491                         if(INPR(side2, side2) > fsTri*fsTri) 
492                         { 
493                                 float tmp = Normalize(side2);
494                                 divs2 = (int)ceil(tmp/fsTri); 
495                         }
496
497                         (*tridivs)[3 * facecounter + 0] = divs1;
498                         (*tridivs)[3 * facecounter + 1] = divs2;
499                         (*tridivs)[3 * facecounter + 2] = divs3;
500                 }
501                 facecounter++;
502         }
503 }
504
505 void smokeModifier_freeDomain(SmokeModifierData *smd)
506 {
507         if(smd->domain)
508         {
509                 // free visualisation buffers
510                 if(smd->domain->bind)
511                 {
512                         glDeleteTextures(smd->domain->max_textures, (GLuint *)smd->domain->bind);
513                         MEM_freeN(smd->domain->bind);
514                 }
515                 smd->domain->max_textures = 0; // unnecessary but let's be sure
516
517                 if(smd->domain->tray)
518                         MEM_freeN(smd->domain->tray);
519                 if(smd->domain->tvox)
520                         MEM_freeN(smd->domain->tvox);
521                 if(smd->domain->traybig)
522                         MEM_freeN(smd->domain->traybig);
523                 if(smd->domain->tvoxbig)
524                         MEM_freeN(smd->domain->tvoxbig);
525
526                 if(smd->domain->fluid)
527                         smoke_free(smd->domain->fluid);
528
529                 if(smd->domain->wt)
530                         smoke_turbulence_free(smd->domain->wt);
531
532                 MEM_freeN(smd->domain);
533                 smd->domain = NULL;
534         }
535 }
536
537 void smokeModifier_freeFlow(SmokeModifierData *smd)
538 {
539         if(smd->flow)
540         {
541 /*
542                 if(smd->flow->bvh)
543                 {
544                         free_bvhtree_from_mesh(smd->flow->bvh);
545                         MEM_freeN(smd->flow->bvh);
546                 }
547                 smd->flow->bvh = NULL;
548 */
549                 MEM_freeN(smd->flow);
550                 smd->flow = NULL;
551         }
552 }
553
554 void smokeModifier_freeCollision(SmokeModifierData *smd)
555 {
556         if(smd->coll)
557         {
558                 if(smd->coll->points)
559                 {
560                         MEM_freeN(smd->coll->points);
561                         smd->coll->points = NULL;
562                 }
563
564                 if(smd->coll->bvhtree)
565                 {
566                         BLI_bvhtree_free(smd->coll->bvhtree);
567                         smd->coll->bvhtree = NULL;
568                 }
569
570                 if(smd->coll->dm)
571                         smd->coll->dm->release(smd->coll->dm);
572                 smd->coll->dm = NULL;
573
574                 MEM_freeN(smd->coll);
575                 smd->coll = NULL;
576         }
577 }
578
579 void smokeModifier_reset(struct SmokeModifierData *smd)
580 {
581         if(smd)
582         {
583                 if(smd->domain)
584                 {
585                         // free visualisation buffers
586                         if(smd->domain->bind)
587                         {
588                                 glDeleteTextures(smd->domain->max_textures, (GLuint *)smd->domain->bind);
589                                 MEM_freeN(smd->domain->bind);
590                                 smd->domain->bind = NULL;
591                         }
592                         smd->domain->max_textures = 0;
593                         if(smd->domain->viewsettings < MOD_SMOKE_VIEW_USEBIG)
594                                 smd->domain->viewsettings = 0;
595                         else
596                                 smd->domain->viewsettings = MOD_SMOKE_VIEW_USEBIG;
597
598                         if(smd->domain->tray)
599                                 MEM_freeN(smd->domain->tray);
600                         if(smd->domain->tvox)
601                                 MEM_freeN(smd->domain->tvox);
602                         if(smd->domain->traybig)
603                                 MEM_freeN(smd->domain->traybig);
604                         if(smd->domain->tvoxbig)
605                                 MEM_freeN(smd->domain->tvoxbig);
606
607                         smd->domain->tvox = NULL;
608                         smd->domain->tray = NULL;
609                         smd->domain->tvoxbig = NULL;
610                         smd->domain->traybig = NULL;
611
612                         if(smd->domain->fluid)
613                         {
614                                 smoke_free(smd->domain->fluid);
615                                 smd->domain->fluid = NULL;
616                         }
617                         
618                         if(smd->domain->wt)
619                         {
620                                 smoke_turbulence_free(smd->domain->wt);
621                                 smd->domain->wt = NULL;
622                         }
623                 }
624                 else if(smd->flow)
625                 {
626                         /*
627                         if(smd->flow->bvh)
628                         {
629                                 free_bvhtree_from_mesh(smd->flow->bvh);
630                                 MEM_freeN(smd->flow->bvh);
631                         }
632                         smd->flow->bvh = NULL;
633                         */
634                 }
635                 else if(smd->coll)
636                 {
637                         if(smd->coll->points)
638                         {
639                                 MEM_freeN(smd->coll->points);
640                                 smd->coll->points = NULL;
641                         }
642
643                         if(smd->coll->bvhtree)
644                         {
645                                 BLI_bvhtree_free(smd->coll->bvhtree);
646                                 smd->coll->bvhtree = NULL;
647                         }
648
649                         if(smd->coll->dm)
650                                 smd->coll->dm->release(smd->coll->dm);
651                         smd->coll->dm = NULL;
652
653                 }
654         }
655 }
656
657 void smokeModifier_free (SmokeModifierData *smd)
658 {
659         if(smd)
660         {
661                 smokeModifier_freeDomain(smd);
662                 smokeModifier_freeFlow(smd);
663                 smokeModifier_freeCollision(smd);
664         }
665 }
666
667 void smokeModifier_createType(struct SmokeModifierData *smd)
668 {
669         if(smd)
670         {
671                 if(smd->type & MOD_SMOKE_TYPE_DOMAIN)
672                 {
673                         if(smd->domain)
674                                 smokeModifier_freeDomain(smd);
675
676                         smd->domain = MEM_callocN(sizeof(SmokeDomainSettings), "SmokeDomain");
677
678                         smd->domain->smd = smd;
679
680                         /* set some standard values */
681                         smd->domain->fluid = NULL;
682                         smd->domain->wt = NULL;
683                         smd->domain->eff_group = NULL;
684                         smd->domain->fluid_group = NULL;
685                         smd->domain->coll_group = NULL;
686                         smd->domain->maxres = 32;
687                         smd->domain->amplify = 1;
688                         smd->domain->omega = 1.0;
689                         smd->domain->alpha = -0.001;
690                         smd->domain->beta = 0.1;
691                         smd->domain->flags = MOD_SMOKE_DISSOLVE_LOG;
692                         smd->domain->strength = 2.0;
693                         smd->domain->noise = MOD_SMOKE_NOISEWAVE;
694                         smd->domain->visibility = 1;
695                         smd->domain->diss_speed = 5;
696
697                         // init 3dview buffer
698                         smd->domain->tvox = NULL;
699                         smd->domain->tray = NULL;
700                         smd->domain->tvoxbig = NULL;
701                         smd->domain->traybig = NULL;
702                         smd->domain->viewsettings = 0;
703                         smd->domain->bind = NULL;
704                         smd->domain->max_textures = 0;
705                 }
706                 else if(smd->type & MOD_SMOKE_TYPE_FLOW)
707                 {
708                         if(smd->flow)
709                                 smokeModifier_freeFlow(smd);
710
711                         smd->flow = MEM_callocN(sizeof(SmokeFlowSettings), "SmokeFlow");
712
713                         smd->flow->smd = smd;
714
715                         /* set some standard values */
716                         smd->flow->density = 1.0;
717                         smd->flow->temp = 1.0;
718
719                         smd->flow->psys = NULL;
720
721                 }
722                 else if(smd->type & MOD_SMOKE_TYPE_COLL)
723                 {
724                         if(smd->coll)
725                                 smokeModifier_freeCollision(smd);
726
727                         smd->coll = MEM_callocN(sizeof(SmokeCollSettings), "SmokeColl");
728
729                         smd->coll->smd = smd;
730                         smd->coll->points = NULL;
731                         smd->coll->numpoints = 0;
732                         smd->coll->bvhtree = NULL;
733                         smd->coll->dm = NULL;
734                 }
735         }
736 }
737
738 // forward declaration
739 void smoke_calc_transparency(struct SmokeModifierData *smd, float *light, int big);
740
741 void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedMesh *dm, int useRenderParams, int isFinalCalc)
742 {       
743         if(scene->r.cfra >= smd->time)
744                 smokeModifier_init(smd, ob, scene, dm);
745
746         if((smd->type & MOD_SMOKE_TYPE_FLOW))
747         {
748                 if(scene->r.cfra > smd->time)
749                 {
750                         // XXX TODO
751                         smd->time = scene->r.cfra;
752
753                         // rigid movement support
754                         /*
755                         Mat4CpyMat4(smd->flow->mat_old, smd->flow->mat);
756                         Mat4CpyMat4(smd->flow->mat, ob->obmat);
757                         */
758                 }
759                 else if(scene->r.cfra < smd->time)
760                 {
761                         smd->time = scene->r.cfra;
762                         smokeModifier_reset(smd);
763                 }
764         }
765         else if(smd->type & MOD_SMOKE_TYPE_COLL)
766         {
767                 if(scene->r.cfra > smd->time)
768                 {
769                         // XXX TODO
770                         smd->time = scene->r.cfra;
771                         
772                         if(smd->coll->dm)
773                                 smd->coll->dm->release(smd->coll->dm);
774
775                         smd->coll->dm = CDDM_copy(dm);
776
777                         // rigid movement support
778                         Mat4CpyMat4(smd->coll->mat_old, smd->coll->mat);
779                         Mat4CpyMat4(smd->coll->mat, ob->obmat);
780                 }
781                 else if(scene->r.cfra < smd->time)
782                 {
783                         smd->time = scene->r.cfra;
784                         smokeModifier_reset(smd);
785                 }
786         }
787         else if(smd->type & MOD_SMOKE_TYPE_DOMAIN)
788         {
789                 SmokeDomainSettings *sds = smd->domain;
790                 
791                 if(scene->r.cfra > smd->time)
792                 {
793                         GroupObject *go = NULL;
794                         Base *base = NULL;
795                         
796                         tstart();
797                         
798                         if(sds->flags & MOD_SMOKE_DISSOLVE)
799                         {
800                                 smoke_dissolve(sds->fluid, sds->diss_speed, sds->flags & MOD_SMOKE_DISSOLVE_LOG);
801                                 
802                                 if(sds->wt)
803                                         smoke_dissolve_wavelet(sds->wt, sds->diss_speed, sds->flags & MOD_SMOKE_DISSOLVE_LOG);
804                         }
805
806                         /* reset view for new frame */
807                         if(sds->viewsettings < MOD_SMOKE_VIEW_USEBIG)
808                                 sds->viewsettings = 0;
809                         else
810                                 sds->viewsettings = MOD_SMOKE_VIEW_USEBIG;
811
812                         // do flows and fluids
813                         if(1)
814                         {
815                                 Object *otherobj = NULL;
816                                 ModifierData *md = NULL;
817
818                                 if(sds->fluid_group) // we use groups since we have 2 domains
819                                         go = sds->fluid_group->gobject.first;
820                                 else
821                                         base = scene->base.first;
822
823                                 while(base || go)
824                                 {
825                                         otherobj = NULL;
826
827                                         if(sds->fluid_group) 
828                                         {
829                                                 if(go->ob)
830                                                         otherobj = go->ob;
831                                         }
832                                         else
833                                                 otherobj = base->object;
834
835                                         if(!otherobj)
836                                         {
837                                                 if(sds->fluid_group)
838                                                         go = go->next;
839                                                 else
840                                                         base= base->next;
841
842                                                 continue;
843                                         }
844
845                                         md = modifiers_findByType(otherobj, eModifierType_Smoke);
846                                         
847                                         // check for active smoke modifier
848                                         if(md && md->mode & (eModifierMode_Realtime | eModifierMode_Render))
849                                         {
850                                                 SmokeModifierData *smd2 = (SmokeModifierData *)md;
851                                                 
852                                                 // check for initialized smoke object
853                                                 if((smd2->type & MOD_SMOKE_TYPE_FLOW) && smd2->flow)
854                                                 {
855                                                         // we got nice flow object
856                                                         SmokeFlowSettings *sfs = smd2->flow;
857                                                         
858                                                         if(sfs->psys && sfs->psys->part && sfs->psys->part->type==PART_EMITTER) // is particle system selected
859                                                         {
860                                                                 ParticleSystem *psys = sfs->psys;
861                                                                 ParticleSettings *part=psys->part;
862                                                                 ParticleData *pa = NULL;
863                                                                 int p = 0;
864                                                                 float *density = smoke_get_density(sds->fluid);
865                                                                 float *bigdensity = smoke_turbulence_get_density(sds->wt);
866                                                                 float *heat = smoke_get_heat(sds->fluid);
867                                                                 float *velocity_x = smoke_get_velocity_x(sds->fluid);
868                                                                 float *velocity_y = smoke_get_velocity_y(sds->fluid);
869                                                                 float *velocity_z = smoke_get_velocity_z(sds->fluid);
870                                                                 unsigned char *obstacle = smoke_get_obstacle(sds->fluid);
871                                                                 int bigres[3];  
872
873                                                                 printf("found flow psys\n");
874                                                                 
875                                                                 // mostly copied from particle code
876                                                                 for(p=0, pa=psys->particles; p<psys->totpart; p++, pa++)
877                                                                 {
878                                                                         int cell[3];
879                                                                         size_t i = 0;
880                                                                         size_t index = 0;
881                                                                         int badcell = 0;
882                                                                         
883                                                                         if(pa->alive == PARS_KILLED) continue;
884                                                                         else if(pa->alive == PARS_UNBORN && (part->flag & PART_UNBORN)==0) continue;
885                                                                         else if(pa->alive == PARS_DEAD && (part->flag & PART_DIED)==0) continue;
886                                                                         else if(pa->flag & (PARS_UNEXIST+PARS_NO_DISP)) continue;
887                                                                         
888                                                                         // VECCOPY(pos, pa->state.co);
889                                                                         // Mat4MulVecfl (ob->imat, pos);
890                                                                         
891                                                                         // 1. get corresponding cell
892                                                                         get_cell(smd, pa->state.co, cell, 0);
893                                                                 
894                                                                         // check if cell is valid (in the domain boundary)
895                                                                         for(i = 0; i < 3; i++)
896                                                                         {
897                                                                                 if((cell[i] > sds->res[i] - 1) || (cell[i] < 0))
898                                                                                 {
899                                                                                         badcell = 1;
900                                                                                         break;
901                                                                                 }
902                                                                         }
903                                                                                 
904                                                                         if(badcell)
905                                                                                 continue;
906                                                                         
907                                                                         // 2. set cell values (heat, density and velocity)
908                                                                         index = smoke_get_index(cell[0], sds->res[0], cell[1], sds->res[1], cell[2]);
909                                                                         
910                                                                         if(!(sfs->type & MOD_SMOKE_FLOW_TYPE_OUTFLOW) && !(obstacle[index] & 2)) // this is inflow
911                                                                         {
912                                                                                 // heat[index] += sfs->temp * 0.1;
913                                                                                 // density[index] += sfs->density * 0.1;
914
915                                                                                 heat[index] = sfs->temp;
916                                                                                 density[index] = sfs->density;
917
918                                                                                 /*
919                                                                                 velocity_x[index] = pa->state.vel[0];
920                                                                                 velocity_y[index] = pa->state.vel[1];
921                                                                                 velocity_z[index] = pa->state.vel[2];
922                                                                                 */
923                                                                                 obstacle[index] |= 2;
924
925                                                                                 // we need different handling for the high-res feature
926                                                                                 if(bigdensity)
927                                                                                 {
928                                                                                         // init all surrounding cells according to amplification, too
929                                                                                         int i, j, k;
930
931                                                                                         smoke_turbulence_get_res(smd->domain->wt, bigres);
932
933                                                                                         for(i = 0; i < smd->domain->amplify + 1; i++)
934                                                                                                 for(j = 0; j < smd->domain->amplify + 1; j++)
935                                                                                                         for(k = 0; k < smd->domain->amplify + 1; k++)
936                                                                                                         {
937                                                                                                                 index = smoke_get_index((smd->domain->amplify + 1)* cell[0] + i, bigres[0], (smd->domain->amplify + 1)* cell[1] + j, bigres[1], (smd->domain->amplify + 1)* cell[2] + k);
938                                                                                                                 bigdensity[index] = sfs->density;
939                                                                                                         }
940                                                                                 }
941                                                                         }
942                                                                         else if(sfs->type & MOD_SMOKE_FLOW_TYPE_OUTFLOW) // outflow
943                                                                         {
944                                                                                 heat[index] = 0.f;
945                                                                                 density[index] = 0.f;
946                                                                                 velocity_x[index] = 0.f;
947                                                                                 velocity_y[index] = 0.f;
948                                                                                 velocity_z[index] = 0.f;
949
950                                                                                 // we need different handling for the high-res feature
951                                                                                 if(bigdensity)
952                                                                                 {
953                                                                                         // init all surrounding cells according to amplification, too
954                                                                                         int i, j, k;
955
956                                                                                         smoke_turbulence_get_res(smd->domain->wt, bigres);
957
958                                                                                         for(i = 0; i < smd->domain->amplify + 1; i++)
959                                                                                                 for(j = 0; j < smd->domain->amplify + 1; j++)
960                                                                                                         for(k = 0; k < smd->domain->amplify + 1; k++)
961                                                                                                         {
962                                                                                                                 index = smoke_get_index((smd->domain->amplify + 1)* cell[0] + i, bigres[0], (smd->domain->amplify + 1)* cell[1] + j, bigres[1], (smd->domain->amplify + 1)* cell[2] + k);
963                                                                                                                 bigdensity[index] = 0.f;
964                                                                                                         }
965                                                                                 }
966                                                                         }
967                                                                 }
968                                                         }
969                                                         else
970                                                         {
971                                                                 /*
972                                                                 for()
973                                                                 {
974                                                                         // no psys
975                                                                         BVHTreeNearest nearest;
976
977                                                                         nearest.index = -1;
978                                                                         nearest.dist = FLT_MAX;
979
980                                                                         BLI_bvhtree_find_nearest(sfs->bvh->tree, pco, &nearest, sfs->bvh->nearest_callback, sfs->bvh);
981                                                                 }*/
982                                                         }
983                                                 }       
984                                         }
985
986                                         if(sds->fluid_group)
987                                                 go = go->next;
988                                         else
989                                                 base= base->next;
990                                 }
991                         }
992
993                         // do effectors
994                         /*
995                         if(sds->eff_group)
996                         {
997                                 for(go = sds->eff_group->gobject.first; go; go = go->next) 
998                                 {
999                                         if(go->ob)
1000                                         {
1001                                                 if(ob->pd)
1002                                                 {
1003                                                         
1004                                                 }
1005                                         }
1006                                 }
1007                         }
1008                         */
1009
1010                         // do collisions        
1011                         if(1)
1012                         {
1013                                 Object *otherobj = NULL;
1014                                 ModifierData *md = NULL;
1015
1016                                 if(sds->coll_group) // we use groups since we have 2 domains
1017                                         go = sds->coll_group->gobject.first;
1018                                 else
1019                                         base = scene->base.first;
1020
1021                                 while(base || go)
1022                                 {
1023                                         otherobj = NULL;
1024
1025                                         if(sds->coll_group) 
1026                                         {
1027                                                 if(go->ob)
1028                                                         otherobj = go->ob;
1029                                         }
1030                                         else
1031                                                 otherobj = base->object;
1032
1033                                         if(!otherobj)
1034                                         {
1035                                                 if(sds->coll_group)
1036                                                         go = go->next;
1037                                                 else
1038                                                         base= base->next;
1039
1040                                                 continue;
1041                                         }
1042                         
1043                                         md = modifiers_findByType(otherobj, eModifierType_Smoke);
1044                                         
1045                                         // check for active smoke modifier
1046                                         if(md && md->mode & (eModifierMode_Realtime | eModifierMode_Render))
1047                                         {
1048                                                 SmokeModifierData *smd2 = (SmokeModifierData *)md;
1049
1050                                                 if((smd2->type & MOD_SMOKE_TYPE_COLL) && smd2->coll)
1051                                                 {
1052                                                         // we got nice collision object
1053                                                         SmokeCollSettings *scs = smd2->coll;
1054                                                         size_t i, j;
1055                                                         unsigned char *obstacles = smoke_get_obstacle(smd->domain->fluid);
1056
1057                                                         for(i = 0; i < scs->numpoints; i++)
1058                                                         {
1059                                                                 int badcell = 0;
1060                                                                 size_t index = 0;
1061                                                                 int cell[3];
1062
1063                                                                 // 1. get corresponding cell
1064                                                                 get_cell(smd, &scs->points[3 * i], cell, 0);
1065                                                         
1066                                                                 // check if cell is valid (in the domain boundary)
1067                                                                 for(j = 0; j < 3; j++)
1068                                                                         if((cell[j] > sds->res[j] - 1) || (cell[j] < 0))
1069                                                                         {
1070                                                                                 badcell = 1;
1071                                                                                 break;
1072                                                                         }
1073                                                                                 
1074                                                                 if(badcell)
1075                                                                         continue;
1076
1077                                                                 // 2. set cell values (heat, density and velocity)
1078                                                                 index = smoke_get_index(cell[0], sds->res[0], cell[1], sds->res[1], cell[2]);
1079                                                                 
1080                                                                 // printf("cell[0]: %d, cell[1]: %d, cell[2]: %d\n", cell[0], cell[1], cell[2]);
1081                                                                 // printf("res[0]: %d, res[1]: %d, res[2]: %d, index: %d\n\n", sds->res[0], sds->res[1], sds->res[2], index);
1082                                                                         
1083                                                                 obstacles[index] = 1;
1084
1085                                                                 // for moving gobstacles
1086                                                                 /*
1087                                                                 const LbmFloat maxVelVal = 0.1666;
1088                                                                 const LbmFloat maxusqr = maxVelVal*maxVelVal*3. *1.5;
1089
1090                                                                 LbmVec objvel = vec2L((mMOIVertices[n]-mMOIVerticesOld[n]) /dvec); { 
1091                                                                 const LbmFloat usqr = (objvel[0]*objvel[0]+objvel[1]*objvel[1]+objvel[2]*objvel[2])*1.5; 
1092                                                                 USQRMAXCHECK(usqr, objvel[0],objvel[1],objvel[2], mMaxVlen, mMxvx,mMxvy,mMxvz); 
1093                                                                 if(usqr>maxusqr) { 
1094                                                                         // cutoff at maxVelVal 
1095                                                                         for(int jj=0; jj<3; jj++) { 
1096                                                                                 if(objvel[jj]>0.) objvel[jj] =  maxVelVal;  
1097                                                                                 if(objvel[jj]<0.) objvel[jj] = -maxVelVal; 
1098                                                                         } 
1099                                                                 } } 
1100
1101                                                                 const LbmFloat dp=dot(objvel, vec2L((*pNormals)[n]) ); 
1102                                                                 const LbmVec oldov=objvel; // debug
1103                                                                 objvel = vec2L((*pNormals)[n]) *dp;
1104                                                                 */
1105                                                         }
1106                                                 }
1107                                         }
1108
1109                                         if(sds->coll_group)
1110                                                 go = go->next;
1111                                         else
1112                                                 base= base->next;
1113                                 }
1114                         }
1115                         
1116                         // set new time
1117                         smd->time = scene->r.cfra;
1118
1119                         // simulate the actual smoke (c++ code in intern/smoke)
1120                         smoke_step(sds->fluid, smd->time);
1121                         if(sds->wt)
1122                                 smoke_turbulence_step(sds->wt, sds->fluid);
1123
1124                         tend();
1125                         printf ( "Frame: %d, Time: %f\n", (int)smd->time, ( float ) tval() );
1126                 }
1127                 else if(scene->r.cfra < smd->time)
1128                 {
1129                         // we got back in time, reset smoke in this case (TODO: use cache later)
1130                         smd->time = scene->r.cfra;
1131                         smokeModifier_reset(smd);
1132                 }
1133         }
1134 }
1135
1136 // update necessary information for 3dview
1137 void smoke_prepare_View(SmokeModifierData *smd, float *light)
1138 {
1139         float *density = NULL;
1140         int x, y, z;
1141
1142         if(!smd->domain->tray)
1143         {
1144                 // TRay is for self shadowing
1145                 smd->domain->tray = MEM_callocN(sizeof(float)*smd->domain->res[0]*smd->domain->res[1]*smd->domain->res[2], "Smoke_tRay");
1146         }
1147         if(!smd->domain->tvox)
1148         {
1149                 // TVox is for tranaparency
1150                 smd->domain->tvox = MEM_callocN(sizeof(float)*smd->domain->res[0]*smd->domain->res[1]*smd->domain->res[2], "Smoke_tVox");
1151         }
1152
1153         // update 3dview
1154         density = smoke_get_density(smd->domain->fluid);
1155         for(x = 0; x < smd->domain->res[0]; x++)
1156                         for(y = 0; y < smd->domain->res[1]; y++)
1157                                 for(z = 0; z < smd->domain->res[2]; z++)
1158                                 {
1159                                         size_t index;
1160
1161                                         index = smoke_get_index(x, smd->domain->res[0], y, smd->domain->res[1], z);
1162                                         // Transparency computation
1163                                         // formula taken from "Visual Simulation of Smoke" / Fedkiw et al. pg. 4
1164                                         // T_vox = exp(-C_ext * h)
1165                                         // C_ext/sigma_t = density * C_ext
1166                                         smoke_set_tvox(smd, index, exp(-density[index] * 7.0 * smd->domain->dx));
1167         }
1168         smoke_calc_transparency(smd, light, 0);
1169 }
1170
1171 // update necessary information for 3dview ("high res" option)
1172 void smoke_prepare_bigView(SmokeModifierData *smd, float *light)
1173 {
1174         float *density = NULL;
1175         size_t i = 0;
1176         int bigres[3];
1177
1178         smoke_turbulence_get_res(smd->domain->wt, bigres);
1179
1180         if(!smd->domain->traybig)
1181         {
1182                 // TRay is for self shadowing
1183                 smd->domain->traybig = MEM_callocN(sizeof(float)*bigres[0]*bigres[1]*bigres[2], "Smoke_tRayBig");
1184         }
1185         if(!smd->domain->tvoxbig)
1186         {
1187                 // TVox is for tranaparency
1188                 smd->domain->tvoxbig = MEM_callocN(sizeof(float)*bigres[0]*bigres[1]*bigres[2], "Smoke_tVoxBig");
1189         }
1190
1191         density = smoke_turbulence_get_density(smd->domain->wt);
1192         for (i = 0; i < bigres[0] * bigres[1] * bigres[2]; i++)
1193         {
1194                 // Transparency computation
1195                 // formula taken from "Visual Simulation of Smoke" / Fedkiw et al. pg. 4
1196                 // T_vox = exp(-C_ext * h)
1197                 // C_ext/sigma_t = density * C_ext
1198                 smoke_set_bigtvox(smd, i, exp(-density[i] * 7.0 * smd->domain->dx / (smd->domain->amplify + 1)) );
1199         }
1200         smoke_calc_transparency(smd, light, 1);
1201 }
1202
1203
1204 float smoke_get_tvox(SmokeModifierData *smd, size_t index)
1205 {
1206         return smd->domain->tvox[index];
1207 }
1208
1209 void smoke_set_tvox(SmokeModifierData *smd, size_t index, float tvox)
1210 {
1211         smd->domain->tvox[index] = tvox;
1212 }
1213
1214 float smoke_get_tray(SmokeModifierData *smd, size_t index)
1215 {
1216         return smd->domain->tray[index];
1217 }
1218
1219 void smoke_set_tray(SmokeModifierData *smd, size_t index, float transparency)
1220 {
1221         smd->domain->tray[index] = transparency;
1222 }
1223
1224 float smoke_get_bigtvox(SmokeModifierData *smd, size_t index)
1225 {
1226         return smd->domain->tvoxbig[index];
1227 }
1228
1229 void smoke_set_bigtvox(SmokeModifierData *smd, size_t index, float tvox)
1230 {
1231         smd->domain->tvoxbig[index] = tvox;
1232 }
1233
1234 float smoke_get_bigtray(SmokeModifierData *smd, size_t index)
1235 {
1236         return smd->domain->traybig[index];
1237 }
1238
1239 void smoke_set_bigtray(SmokeModifierData *smd, size_t index, float transparency)
1240 {
1241         smd->domain->traybig[index] = transparency;
1242 }
1243
1244 long long smoke_get_mem_req(int xres, int yres, int zres, int amplify)
1245 {
1246           int totalCells = xres * yres * zres;
1247           int amplifiedCells = totalCells * amplify * amplify * amplify;
1248
1249           // print out memory requirements
1250           long long int coarseSize = sizeof(float) * totalCells * 22 +
1251                            sizeof(unsigned char) * totalCells;
1252
1253           long long int fineSize = sizeof(float) * amplifiedCells * 7 + // big grids
1254                          sizeof(float) * totalCells * 8 +     // small grids
1255                          sizeof(float) * 128 * 128 * 128;     // noise tile
1256
1257           long long int totalMB = (coarseSize + fineSize) / (1024 * 1024);
1258
1259           return totalMB;
1260 }
1261
1262
1263 static void calc_voxel_transp(SmokeModifierData *smd, int *pixel, float *tRay)
1264 {
1265         // printf("Pixel(%d, %d, %d)\n", pixel[0], pixel[1], pixel[2]);
1266         const size_t index = smoke_get_index(pixel[0], smd->domain->res[0], pixel[1], smd->domain->res[1], pixel[2]);
1267
1268         // T_ray *= T_vox
1269         *tRay *= smoke_get_tvox(smd, index);
1270 }
1271
1272 static void calc_voxel_transp_big(SmokeModifierData *smd, int *pixel, float *tRay)
1273 {
1274         int bigres[3];
1275         size_t index;
1276
1277         smoke_turbulence_get_res(smd->domain->wt, bigres);
1278         index = smoke_get_index(pixel[0], bigres[0], pixel[1], bigres[1], pixel[2]);
1279
1280         /*
1281         if(index > bigres[0]*bigres[1]*bigres[2])
1282                 printf("pixel[0]: %d, [1]: %d, [2]: %d\n", pixel[0], pixel[1], pixel[2]);
1283         */
1284
1285         // T_ray *= T_vox
1286         *tRay *= smoke_get_bigtvox(smd, index);
1287 }
1288
1289 static void bresenham_linie_3D(SmokeModifierData *smd, int x1, int y1, int z1, int x2, int y2, int z2, float *tRay, int big)
1290 {
1291     int dx, dy, dz, i, l, m, n, x_inc, y_inc, z_inc, err_1, err_2, dx2, dy2, dz2;
1292     int pixel[3];
1293
1294     pixel[0] = x1;
1295     pixel[1] = y1;
1296     pixel[2] = z1;
1297
1298     dx = x2 - x1;
1299     dy = y2 - y1;
1300     dz = z2 - z1;
1301
1302     x_inc = (dx < 0) ? -1 : 1;
1303     l = abs(dx);
1304     y_inc = (dy < 0) ? -1 : 1;
1305     m = abs(dy);
1306     z_inc = (dz < 0) ? -1 : 1;
1307     n = abs(dz);
1308     dx2 = l << 1;
1309     dy2 = m << 1;
1310     dz2 = n << 1;
1311
1312     if ((l >= m) && (l >= n)) {
1313         err_1 = dy2 - l;
1314         err_2 = dz2 - l;
1315         for (i = 0; i < l; i++) {
1316                 if(!big)
1317                                 calc_voxel_transp(smd, pixel, tRay);
1318                         else
1319                                 calc_voxel_transp_big(smd, pixel, tRay);
1320                 if(*tRay < 0.0f)
1321                         return;
1322             if (err_1 > 0) {
1323                 pixel[1] += y_inc;
1324                 err_1 -= dx2;
1325             }
1326             if (err_2 > 0) {
1327                 pixel[2] += z_inc;
1328                 err_2 -= dx2;
1329             }
1330             err_1 += dy2;
1331             err_2 += dz2;
1332             pixel[0] += x_inc;
1333         }
1334     } else if ((m >= l) && (m >= n)) {
1335         err_1 = dx2 - m;
1336         err_2 = dz2 - m;
1337         for (i = 0; i < m; i++) {
1338                 if(!big)
1339                                 calc_voxel_transp(smd, pixel, tRay);
1340                         else
1341                                 calc_voxel_transp_big(smd, pixel, tRay);
1342                 if(*tRay < 0.0f)
1343                         return;
1344             if (err_1 > 0) {
1345                 pixel[0] += x_inc;
1346                 err_1 -= dy2;
1347             }
1348             if (err_2 > 0) {
1349                 pixel[2] += z_inc;
1350                 err_2 -= dy2;
1351             }
1352             err_1 += dx2;
1353             err_2 += dz2;
1354             pixel[1] += y_inc;
1355         }
1356     } else {
1357         err_1 = dy2 - n;
1358         err_2 = dx2 - n;
1359         for (i = 0; i < n; i++) {
1360                 if(!big)
1361                                 calc_voxel_transp(smd, pixel, tRay);
1362                         else
1363                                 calc_voxel_transp_big(smd, pixel, tRay);
1364                 if(*tRay < 0.0f)
1365                         return;
1366             if (err_1 > 0) {
1367                 pixel[1] += y_inc;
1368                 err_1 -= dz2;
1369             }
1370             if (err_2 > 0) {
1371                 pixel[0] += x_inc;
1372                 err_2 -= dz2;
1373             }
1374             err_1 += dy2;
1375             err_2 += dx2;
1376             pixel[2] += z_inc;
1377         }
1378     }
1379     if(!big)
1380         calc_voxel_transp(smd, pixel, tRay);
1381     else
1382         calc_voxel_transp_big(smd, pixel, tRay);
1383 }
1384
1385 static void get_cell(struct SmokeModifierData *smd, float *pos, int *cell, int correct)
1386 {
1387         float tmp[3];
1388
1389         VECSUB(tmp, pos, smd->domain->p0);
1390         VecMulf(tmp, 1.0 / smd->domain->dx);
1391
1392         if(correct)
1393         {
1394                 cell[0] = MIN2(smd->domain->res[0] - 1, MAX2(0, (int)floor(tmp[0])));
1395                 cell[1] = MIN2(smd->domain->res[1] - 1, MAX2(0, (int)floor(tmp[1])));
1396                 cell[2] = MIN2(smd->domain->res[2] - 1, MAX2(0, (int)floor(tmp[2])));
1397         }
1398         else
1399         {
1400                 cell[0] = (int)floor(tmp[0]);
1401                 cell[1] = (int)floor(tmp[1]);
1402                 cell[2] = (int)floor(tmp[2]);
1403         }
1404 }
1405 static void get_bigcell(struct SmokeModifierData *smd, float *pos, int *cell, int correct)
1406 {
1407         float tmp[3];
1408         int res[3];
1409         smoke_turbulence_get_res(smd->domain->wt, res);
1410
1411         VECSUB(tmp, pos, smd->domain->p0);
1412
1413         VecMulf(tmp, (smd->domain->amplify + 1)/ smd->domain->dx );
1414
1415         if(correct)
1416         {
1417                 cell[0] = MIN2(res[0] - 1, MAX2(0, (int)floor(tmp[0])));
1418                 cell[1] = MIN2(res[1] - 1, MAX2(0, (int)floor(tmp[1])));
1419                 cell[2] = MIN2(res[2] - 1, MAX2(0, (int)floor(tmp[2])));
1420         }
1421         else
1422         {
1423                 cell[0] = (int)floor(tmp[0]);
1424                 cell[1] = (int)floor(tmp[1]);
1425                 cell[2] = (int)floor(tmp[2]);
1426         }
1427 }
1428
1429
1430 void smoke_calc_transparency(struct SmokeModifierData *smd, float *light, int big)
1431 {
1432         int x, y, z;
1433         float bv[6];
1434         int res[3];
1435         float bigfactor = 1.0;
1436
1437         // x
1438         bv[0] = smd->domain->p0[0];
1439         bv[1] = smd->domain->p1[0];
1440         // y
1441         bv[2] = smd->domain->p0[1];
1442         bv[3] = smd->domain->p1[1];
1443         // z
1444         bv[4] = smd->domain->p0[2];
1445         bv[5] = smd->domain->p1[2];
1446 /*
1447         printf("bv[0]: %f, [1]: %f, [2]: %f, [3]: %f, [4]: %f, [5]: %f\n", bv[0], bv[1], bv[2], bv[3], bv[4], bv[5]);
1448
1449         printf("p0[0]: %f, p0[1]: %f, p0[2]: %f\n", smd->domain->p0[0], smd->domain->p0[1], smd->domain->p0[2]);
1450         printf("p1[0]: %f, p1[1]: %f, p1[2]: %f\n", smd->domain->p1[0], smd->domain->p1[1], smd->domain->p1[2]);
1451         printf("dx: %f, amp: %d\n", smd->domain->dx, smd->domain->amplify);
1452 */
1453         if(!big)
1454         {
1455                 res[0] = smd->domain->res[0];
1456                 res[1] = smd->domain->res[1];
1457                 res[2] = smd->domain->res[2];
1458         }
1459         else
1460         {
1461                 smoke_turbulence_get_res(smd->domain->wt, res);
1462                 bigfactor = 1.0 / (smd->domain->amplify + 1);
1463         }
1464
1465 #pragma omp parallel for schedule(static) private(y, z) shared(big, smd, light, res, bigfactor)
1466         for(x = 0; x < res[0]; x++)
1467                 for(y = 0; y < res[1]; y++)
1468                         for(z = 0; z < res[2]; z++)
1469                         {
1470                                 float voxelCenter[3];
1471                                 size_t index;
1472                                 float pos[3];
1473                                 int cell[3];
1474                                 float tRay = 1.0;
1475
1476                                 index = smoke_get_index(x, res[0], y, res[1], z);
1477
1478                                 // voxelCenter = m_voxelarray[i].GetCenter();
1479                                 voxelCenter[0] = smd->domain->p0[0] + smd->domain->dx * bigfactor * x + smd->domain->dx * bigfactor * 0.5;
1480                                 voxelCenter[1] = smd->domain->p0[1] + smd->domain->dx * bigfactor * y + smd->domain->dx * bigfactor * 0.5;
1481                                 voxelCenter[2] = smd->domain->p0[2] + smd->domain->dx * bigfactor * z + smd->domain->dx * bigfactor * 0.5;
1482
1483                                 // printf("vc[0]: %f, vc[1]: %f, vc[2]: %f\n", voxelCenter[0], voxelCenter[1], voxelCenter[2]);
1484                                 // printf("light[0]: %f, light[1]: %f, light[2]: %f\n", light[0], light[1], light[2]);
1485
1486                                 // get starting position (in voxel coords)
1487                                 if(BLI_bvhtree_bb_raycast(bv, light, voxelCenter, pos) > FLT_EPSILON)
1488                                 {
1489                                         // we're ouside
1490                                         // printf("out: pos[0]: %f, pos[1]: %f, pos[2]: %f\n", pos[0], pos[1], pos[2]);
1491                                         if(!big)
1492                                                 get_cell(smd, pos, cell, 1);
1493                                         else
1494                                                 get_bigcell(smd, pos, cell, 1);
1495                                 }
1496                                 else
1497                                 {
1498                                         // printf("in: pos[0]: %f, pos[1]: %f, pos[2]: %f\n", light[0], light[1], light[2]);
1499                                         // we're inside
1500                                         if(!big)
1501                                                 get_cell(smd, light, cell, 1);
1502                                         else
1503                                                 get_bigcell(smd, light, cell, 1);
1504                                 }
1505
1506                                 // printf("cell - [0]: %d, [1]: %d, [2]: %d\n", cell[0], cell[1], cell[2]);
1507                                 bresenham_linie_3D(smd, cell[0], cell[1], cell[2], x, y, z, &tRay, big);
1508
1509                                 if(!big)
1510                                         smoke_set_tray(smd, index, tRay);
1511                                 else
1512                                         smoke_set_bigtray(smd, index, tRay);
1513                         }
1514 }
1515