68d055f9eda202299efc7d2e9cc9235e650c7f0e
[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, 2.5 / FPS);
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                 return 1;
233         }
234         else if((smd->type & MOD_SMOKE_TYPE_COLL))
235         {
236                 smd->time = scene->r.cfra;
237
238                 // todo: delete this when loading colls work -dg
239                 if(!smd->coll)
240                         smokeModifier_createType(smd);
241
242                 if(!smd->coll->points)
243                 {
244                         // init collision points
245                         SmokeCollSettings *scs = smd->coll;
246                         MVert *mvert = dm->getVertArray(dm);
247                         MFace *mface = dm->getFaceArray(dm);
248                         size_t i = 0, divs = 0;
249                         int *tridivs = NULL;
250                         float cell_len = 1.0 / 50.0; // for res = 50
251                         size_t newdivs = 0;
252                         size_t max_points = 0;
253                         size_t quads = 0, facecounter = 0;
254
255                         // copy obmat
256                         Mat4CpyMat4(scs->mat, ob->obmat);
257                         Mat4CpyMat4(scs->mat_old, ob->obmat);
258
259                         // count quads
260                         for(i = 0; i < dm->getNumFaces(dm); i++)
261                         {
262                                 if(mface[i].v4)
263                                         quads++;
264                         }
265
266                         calcTriangleDivs(ob, mvert, dm->getNumVerts(dm), mface,  dm->getNumFaces(dm), dm->getNumFaces(dm) + quads, &tridivs, cell_len);
267
268                         // count triangle divisions
269                         for(i = 0; i < dm->getNumFaces(dm) + quads; i++)
270                         {
271                                 divs += (tridivs[3 * i] + 1) * (tridivs[3 * i + 1] + 1) * (tridivs[3 * i + 2] + 1);
272                         }
273
274                         // printf("divs: %d\n", divs);
275
276                         scs->points = MEM_callocN(sizeof(float) * (dm->getNumVerts(dm) + divs) * 3, "SmokeCollPoints");
277
278                         for(i = 0; i < dm->getNumVerts(dm); i++)
279                         {
280                                 float tmpvec[3];
281                                 VECCOPY(tmpvec, mvert[i].co);
282                                 Mat4MulVecfl (ob->obmat, tmpvec);
283                                 VECCOPY(&scs->points[i * 3], tmpvec);
284                         }
285                         
286                         for(i = 0, facecounter = 0; i < dm->getNumFaces(dm); i++)
287                         {
288                                 int again = 0;
289                                 do
290                                 {
291                                         size_t j, k;
292                                         int divs1 = tridivs[3 * facecounter + 0];
293                                         int divs2 = tridivs[3 * facecounter + 1];
294                                         int divs3 = tridivs[3 * facecounter + 2];
295                                         float side1[3], side2[3], trinormorg[3], trinorm[3];
296                                         
297                                         if(again == 1 && mface[i].v4)
298                                         {
299                                                 VECSUB(side1,  mvert[ mface[i].v3 ].co, mvert[ mface[i].v1 ].co);
300                                                 VECSUB(side2,  mvert[ mface[i].v4 ].co, mvert[ mface[i].v1 ].co);
301                                         }
302                                         else
303                                         {
304                                                 VECSUB(side1,  mvert[ mface[i].v2 ].co, mvert[ mface[i].v1 ].co);
305                                                 VECSUB(side2,  mvert[ mface[i].v3 ].co, mvert[ mface[i].v1 ].co);
306                                         }
307
308                                         Crossf(trinormorg, side1, side2);
309                                         Normalize(trinormorg);
310                                         VECCOPY(trinorm, trinormorg);
311                                         VecMulf(trinorm, 0.25 * cell_len);
312
313                                         for(j = 0; j <= divs1; j++)
314                                         {
315                                                 for(k = 0; k <= divs2; k++)
316                                                 {
317                                                         float p1[3], p2[3], p3[3], p[3]={0,0,0}; 
318                                                         const float uf = (float)(j + TRI_UVOFFSET) / (float)(divs1 + 0.0);
319                                                         const float vf = (float)(k + TRI_UVOFFSET) / (float)(divs2 + 0.0);
320                                                         float tmpvec[3];
321                                                         
322                                                         if(uf+vf > 1.0) 
323                                                         {
324                                                                 // printf("bigger - divs1: %d, divs2: %d\n", divs1, divs2);
325                                                                 continue;
326                                                         }
327
328                                                         VECCOPY(p1, mvert[ mface[i].v1 ].co);
329                                                         if(again == 1 && mface[i].v4)
330                                                         {
331                                                                 VECCOPY(p2, mvert[ mface[i].v3 ].co);
332                                                                 VECCOPY(p3, mvert[ mface[i].v4 ].co);
333                                                         }
334                                                         else
335                                                         {
336                                                                 VECCOPY(p2, mvert[ mface[i].v2 ].co);
337                                                                 VECCOPY(p3, mvert[ mface[i].v3 ].co);
338                                                         }
339
340                                                         VecMulf(p1, (1.0-uf-vf));
341                                                         VecMulf(p2, uf);
342                                                         VecMulf(p3, vf);
343                                                         
344                                                         VECADD(p, p1, p2);
345                                                         VECADD(p, p, p3);
346
347                                                         if(newdivs > divs)
348                                                                 printf("mem problem\n");
349
350                                                         // mMovPoints.push_back(p + trinorm);
351                                                         VECCOPY(tmpvec, p);
352                                                         VECADD(tmpvec, tmpvec, trinorm);
353                                                         Mat4MulVecfl (ob->obmat, tmpvec);
354                                                         VECCOPY(&scs->points[3 * (dm->getNumVerts(dm) + newdivs)], tmpvec);
355                                                         newdivs++;
356
357                                                         if(newdivs > divs)
358                                                                 printf("mem problem\n");
359
360                                                         // mMovPoints.push_back(p - trinorm);
361                                                         VECCOPY(tmpvec, p);
362                                                         VECSUB(tmpvec, tmpvec, trinorm);
363                                                         Mat4MulVecfl (ob->obmat, tmpvec);
364                                                         VECCOPY(&scs->points[3 * (dm->getNumVerts(dm) + newdivs)], tmpvec);
365                                                         newdivs++;
366                                                 }
367                                         }
368
369                                         if(again == 0 && mface[i].v4)
370                                                 again++;
371                                         else
372                                                 again = 0;
373
374                                         facecounter++;
375
376                                 } while(again!=0);
377                         }
378
379                         scs->numpoints = dm->getNumVerts(dm) + newdivs;
380
381                         MEM_freeN(tridivs);
382                 }
383
384                 if(!smd->coll->bvhtree)
385                 {
386                         smd->coll->bvhtree = NULL; // bvhtree_build_from_smoke ( ob->obmat, dm->getFaceArray(dm), dm->getNumFaces(dm), dm->getVertArray(dm), dm->getNumVerts(dm), 0.0 );
387                 }
388
389         }
390
391         return 0;
392 }
393
394 /*! init triangle divisions */
395 void calcTriangleDivs(Object *ob, MVert *verts, int numverts, MFace *faces, int numfaces, int numtris, int **tridivs, float cell_len) 
396 {
397         // mTriangleDivs1.resize( faces.size() );
398         // mTriangleDivs2.resize( faces.size() );
399         // mTriangleDivs3.resize( faces.size() );
400
401         size_t i = 0, facecounter = 0;
402         float maxscale[3] = {1,1,1}; // = channelFindMaxVf(mcScale);
403         float maxpart = ABS(maxscale[0]);
404         float scaleFac = 0;
405         float fsTri = 0;
406         if(ABS(maxscale[1])>maxpart) maxpart = ABS(maxscale[1]);
407         if(ABS(maxscale[2])>maxpart) maxpart = ABS(maxscale[2]);
408         scaleFac = 1.0 / maxpart;
409         // featureSize = mLevel[mMaxRefine].nodeSize
410         fsTri = cell_len * 0.5 * scaleFac;
411
412         if(*tridivs)
413                 MEM_freeN(*tridivs);
414
415         *tridivs = MEM_callocN(sizeof(int) * numtris * 3, "Smoke_Tridivs");
416
417         for(i = 0, facecounter = 0; i < numfaces; i++) 
418         {
419                 float p0[3], p1[3], p2[3];
420                 float side1[3];
421                 float side2[3];
422                 float side3[3];
423                 int divs1=0, divs2=0, divs3=0;
424
425                 VECCOPY(p0, verts[faces[i].v1].co);
426                 Mat4MulVecfl (ob->obmat, p0);
427                 VECCOPY(p1, verts[faces[i].v2].co);
428                 Mat4MulVecfl (ob->obmat, p1);
429                 VECCOPY(p2, verts[faces[i].v3].co);
430                 Mat4MulVecfl (ob->obmat, p2);
431
432                 VECSUB(side1, p1, p0);
433                 VECSUB(side2, p2, p0);
434                 VECSUB(side3, p1, p2);
435
436                 if(INPR(side1, side1) > fsTri*fsTri) 
437                 { 
438                         float tmp = Normalize(side1);
439                         divs1 = (int)ceil(tmp/fsTri); 
440                 }
441                 if(INPR(side2, side2) > fsTri*fsTri) 
442                 { 
443                         float tmp = Normalize(side2);
444                         divs2 = (int)ceil(tmp/fsTri); 
445                         
446                         /*
447                         // debug
448                         if(i==0)
449                                 printf("b tmp: %f, fsTri: %f, divs2: %d\n", tmp, fsTri, divs2);
450                         */
451                 }
452
453                 (*tridivs)[3 * facecounter + 0] = divs1;
454                 (*tridivs)[3 * facecounter + 1] = divs2;
455                 (*tridivs)[3 * facecounter + 2] = divs3;
456
457                 // TODO quad case
458                 if(faces[i].v4)
459                 {
460                         divs1=0, divs2=0, divs3=0;
461
462                         facecounter++;
463                         
464                         VECCOPY(p0, verts[faces[i].v3].co);
465                         Mat4MulVecfl (ob->obmat, p0);
466                         VECCOPY(p1, verts[faces[i].v4].co);
467                         Mat4MulVecfl (ob->obmat, p1);
468                         VECCOPY(p2, verts[faces[i].v1].co);
469                         Mat4MulVecfl (ob->obmat, p2);
470
471                         VECSUB(side1, p1, p0);
472                         VECSUB(side2, p2, p0);
473                         VECSUB(side3, p1, p2);
474
475                         if(INPR(side1, side1) > fsTri*fsTri) 
476                         { 
477                                 float tmp = Normalize(side1);
478                                 divs1 = (int)ceil(tmp/fsTri); 
479                         }
480                         if(INPR(side2, side2) > fsTri*fsTri) 
481                         { 
482                                 float tmp = Normalize(side2);
483                                 divs2 = (int)ceil(tmp/fsTri); 
484                         }
485
486                         (*tridivs)[3 * facecounter + 0] = divs1;
487                         (*tridivs)[3 * facecounter + 1] = divs2;
488                         (*tridivs)[3 * facecounter + 2] = divs3;
489                 }
490                 facecounter++;
491         }
492 }
493
494 void smokeModifier_freeDomain(SmokeModifierData *smd)
495 {
496         if(smd->domain)
497         {
498                 // free visualisation buffers
499                 if(smd->domain->bind)
500                 {
501                         glDeleteTextures(smd->domain->max_textures, (GLuint *)smd->domain->bind);
502                         MEM_freeN(smd->domain->bind);
503                 }
504                 smd->domain->max_textures = 0; // unnecessary but let's be sure
505
506                 if(smd->domain->tray)
507                         MEM_freeN(smd->domain->tray);
508                 if(smd->domain->tvox)
509                         MEM_freeN(smd->domain->tvox);
510                 if(smd->domain->traybig)
511                         MEM_freeN(smd->domain->traybig);
512                 if(smd->domain->tvoxbig)
513                         MEM_freeN(smd->domain->tvoxbig);
514
515                 if(smd->domain->fluid)
516                         smoke_free(smd->domain->fluid);
517
518                 if(smd->domain->wt)
519                         smoke_turbulence_free(smd->domain->wt);
520
521                 MEM_freeN(smd->domain);
522                 smd->domain = NULL;
523         }
524 }
525
526 void smokeModifier_freeFlow(SmokeModifierData *smd)
527 {
528         if(smd->flow)
529         {
530                 MEM_freeN(smd->flow);
531                 smd->flow = NULL;
532         }
533 }
534
535 void smokeModifier_freeCollision(SmokeModifierData *smd)
536 {
537         if(smd->coll)
538         {
539                 if(smd->coll->points)
540                 {
541                         MEM_freeN(smd->coll->points);
542                         smd->coll->points = NULL;
543                 }
544
545                 if(smd->coll->bvhtree)
546                 {
547                         BLI_bvhtree_free(smd->coll->bvhtree);
548                         smd->coll->bvhtree = NULL;
549                 }
550
551                 if(smd->coll->dm)
552                         smd->coll->dm->release(smd->coll->dm);
553                 smd->coll->dm = NULL;
554
555                 MEM_freeN(smd->coll);
556                 smd->coll = NULL;
557         }
558 }
559
560 void smokeModifier_reset(struct SmokeModifierData *smd)
561 {
562         if(smd)
563         {
564                 if(smd->domain)
565                 {
566                         // free visualisation buffers
567                         if(smd->domain->bind)
568                         {
569                                 glDeleteTextures(smd->domain->max_textures, (GLuint *)smd->domain->bind);
570                                 MEM_freeN(smd->domain->bind);
571                                 smd->domain->bind = NULL;
572                         }
573                         smd->domain->max_textures = 0;
574                         if(smd->domain->viewsettings < MOD_SMOKE_VIEW_USEBIG)
575                                 smd->domain->viewsettings = 0;
576                         else
577                                 smd->domain->viewsettings = MOD_SMOKE_VIEW_USEBIG;
578
579                         if(smd->domain->tray)
580                                 MEM_freeN(smd->domain->tray);
581                         if(smd->domain->tvox)
582                                 MEM_freeN(smd->domain->tvox);
583                         if(smd->domain->traybig)
584                                 MEM_freeN(smd->domain->traybig);
585                         if(smd->domain->tvoxbig)
586                                 MEM_freeN(smd->domain->tvoxbig);
587
588                         smd->domain->tvox = NULL;
589                         smd->domain->tray = NULL;
590                         smd->domain->tvoxbig = NULL;
591                         smd->domain->traybig = NULL;
592
593                         if(smd->domain->fluid)
594                         {
595                                 smoke_free(smd->domain->fluid);
596                                 smd->domain->fluid = NULL;
597                         }
598                         
599                         if(smd->domain->wt)
600                         {
601                                 smoke_turbulence_free(smd->domain->wt);
602                                 smd->domain->wt = NULL;
603                         }
604                 }
605                 else if(smd->flow)
606                 {
607                                                 
608                 }
609                 else if(smd->coll)
610                 {
611                         if(smd->coll->points)
612                         {
613                                 MEM_freeN(smd->coll->points);
614                                 smd->coll->points = NULL;
615                         }
616
617                         if(smd->coll->bvhtree)
618                         {
619                                 BLI_bvhtree_free(smd->coll->bvhtree);
620                                 smd->coll->bvhtree = NULL;
621                         }
622
623                         if(smd->coll->dm)
624                                 smd->coll->dm->release(smd->coll->dm);
625                         smd->coll->dm = NULL;
626
627                 }
628         }
629 }
630
631 void smokeModifier_free (SmokeModifierData *smd)
632 {
633         if(smd)
634         {
635                 smokeModifier_freeDomain(smd);
636                 smokeModifier_freeFlow(smd);
637                 smokeModifier_freeCollision(smd);
638         }
639 }
640
641 void smokeModifier_createType(struct SmokeModifierData *smd)
642 {
643         if(smd)
644         {
645                 if(smd->type & MOD_SMOKE_TYPE_DOMAIN)
646                 {
647                         if(smd->domain)
648                                 smokeModifier_freeDomain(smd);
649
650                         smd->domain = MEM_callocN(sizeof(SmokeDomainSettings), "SmokeDomain");
651
652                         smd->domain->smd = smd;
653
654                         /* set some standard values */
655                         smd->domain->fluid = NULL;
656                         smd->domain->wt = NULL;
657                         smd->domain->eff_group = NULL;
658                         smd->domain->fluid_group = NULL;
659                         smd->domain->coll_group = NULL;
660                         smd->domain->maxres = 32;
661                         smd->domain->amplify = 1;
662                         smd->domain->omega = 1.0;
663                         smd->domain->alpha = -0.001;
664                         smd->domain->beta = 0.1;
665                         smd->domain->flags = 0;
666                         smd->domain->strength = 2.0;
667                         smd->domain->noise = MOD_SMOKE_NOISEWAVE;
668                         smd->domain->visibility = 1;
669
670                         // init 3dview buffer
671                         smd->domain->tvox = NULL;
672                         smd->domain->tray = NULL;
673                         smd->domain->tvoxbig = NULL;
674                         smd->domain->traybig = NULL;
675                         smd->domain->viewsettings = 0;
676                         smd->domain->bind = NULL;
677                         smd->domain->max_textures = 0;
678                 }
679                 else if(smd->type & MOD_SMOKE_TYPE_FLOW)
680                 {
681                         if(smd->flow)
682                                 smokeModifier_freeFlow(smd);
683
684                         smd->flow = MEM_callocN(sizeof(SmokeFlowSettings), "SmokeFlow");
685
686                         smd->flow->smd = smd;
687
688                         /* set some standard values */
689                         smd->flow->density = 1.0;
690                         smd->flow->temp = 1.0;
691
692                         smd->flow->psys = NULL;
693
694                 }
695                 else if(smd->type & MOD_SMOKE_TYPE_COLL)
696                 {
697                         if(smd->coll)
698                                 smokeModifier_freeCollision(smd);
699
700                         smd->coll = MEM_callocN(sizeof(SmokeCollSettings), "SmokeColl");
701
702                         smd->coll->smd = smd;
703                         smd->coll->points = NULL;
704                         smd->coll->numpoints = 0;
705                         smd->coll->bvhtree = NULL;
706                         smd->coll->dm = NULL;
707                 }
708         }
709 }
710
711 // forward declaration
712 void smoke_calc_transparency(struct SmokeModifierData *smd, float *light, int big);
713
714 void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedMesh *dm, int useRenderParams, int isFinalCalc)
715 {       
716         if(scene->r.cfra >= smd->time)
717                 smokeModifier_init(smd, ob, scene, dm);
718
719         if((smd->type & MOD_SMOKE_TYPE_FLOW))
720         {
721                 if(scene->r.cfra > smd->time)
722                 {
723                         // XXX TODO
724                         smd->time = scene->r.cfra;
725                 }
726                 else if(scene->r.cfra < smd->time)
727                 {
728                         smd->time = scene->r.cfra;
729                         smokeModifier_reset(smd);
730                 }
731         }
732         else if(smd->type & MOD_SMOKE_TYPE_COLL)
733         {
734                 if(scene->r.cfra > smd->time)
735                 {
736                         // XXX TODO
737                         smd->time = scene->r.cfra;
738                         
739                         if(smd->coll->dm)
740                                 smd->coll->dm->release(smd->coll->dm);
741
742                         smd->coll->dm = CDDM_copy(dm);
743
744                         // rigid movement support
745                         Mat4CpyMat4(smd->coll->mat_old, smd->coll->mat);
746                         Mat4CpyMat4(smd->coll->mat, ob->obmat);
747                 }
748                 else if(scene->r.cfra < smd->time)
749                 {
750                         smd->time = scene->r.cfra;
751                         smokeModifier_reset(smd);
752                 }
753         }
754         else if(smd->type & MOD_SMOKE_TYPE_DOMAIN)
755         {
756                 SmokeDomainSettings *sds = smd->domain;
757                 
758                 if(scene->r.cfra > smd->time)
759                 {
760                         GroupObject *go = NULL;
761                         Base *base = NULL;
762                         int cnt_domain = 0;
763                         
764                         tstart();
765
766                         /* reset view for new frame */
767                         if(sds->viewsettings < MOD_SMOKE_VIEW_USEBIG)
768                                 sds->viewsettings = 0;
769                         else
770                                 sds->viewsettings = MOD_SMOKE_VIEW_USEBIG;
771
772                         /* check for 2nd domain, if not there -> no groups are necessary */
773                         for(base = scene->base.first; base; base= base->next) 
774                         {
775                                 Object *ob1= base->object;
776                                 
777                                 if(ob1 && ob1 != ob)
778                                 {
779                                         ModifierData *tmd = modifiers_findByType(ob1, eModifierType_Smoke);
780
781                                         if(tmd && tmd->mode & (eModifierMode_Realtime | eModifierMode_Render))
782                                         {
783                                                 SmokeModifierData *tsmd = (SmokeModifierData *)tmd;
784
785                                                 if((tsmd->type & MOD_SMOKE_TYPE_DOMAIN))
786                                                 {
787                                                         cnt_domain++;
788                                                 }
789                                         }
790                                 }
791                         }
792
793                         // do flows and fluids
794                         if(sds->fluid_group || !cnt_domain)
795                         {
796                                 Object *otherobj = NULL;
797                                 ModifierData *md = NULL;
798
799                                 if(cnt_domain && !sds->fluid_group) // we use groups since we have 2 domains
800                                         go = sds->fluid_group->gobject.first;
801                                 else
802                                         base = scene->base.first;
803
804                                 while(base || go)
805                                 {
806                                         otherobj = NULL;
807
808                                         if(cnt_domain && !sds->fluid_group) 
809                                         {
810                                                 if(go->ob)
811                                                         otherobj = go->ob;
812                                         }
813                                         else
814                                                 otherobj = base->object;
815
816                                         if(!otherobj)
817                                         {
818                                                 if(cnt_domain && !sds->fluid_group)
819                                                         go = go->next;
820                                                 else
821                                                         base= base->next;
822
823                                                 continue;
824                                         }
825
826                                         md = modifiers_findByType(otherobj, eModifierType_Smoke);
827                                         
828                                         // check for active smoke modifier
829                                         if(md && md->mode & (eModifierMode_Realtime | eModifierMode_Render))
830                                         {
831                                                 SmokeModifierData *smd2 = (SmokeModifierData *)md;
832                                                 
833                                                 // check for initialized smoke object
834                                                 if((smd2->type & MOD_SMOKE_TYPE_FLOW) && smd2->flow)
835                                                 {
836                                                         // we got nice flow object
837                                                         SmokeFlowSettings *sfs = smd2->flow;
838                                                         
839                                                         if(sfs->psys && sfs->psys->part && sfs->psys->part->type==PART_EMITTER) // is particle system selected
840                                                         {
841                                                                 ParticleSystem *psys = sfs->psys;
842                                                                 ParticleSettings *part=psys->part;
843                                                                 ParticleData *pa = NULL;
844                                                                 int p = 0;
845                                                                 float *density = smoke_get_density(sds->fluid);
846                                                                 float *bigdensity = smoke_turbulence_get_density(sds->wt);
847                                                                 float *heat = smoke_get_heat(sds->fluid);
848                                                                 float *velocity_x = smoke_get_velocity_x(sds->fluid);
849                                                                 float *velocity_y = smoke_get_velocity_y(sds->fluid);
850                                                                 float *velocity_z = smoke_get_velocity_z(sds->fluid);
851                                                                 int bigres[3];                                                          
852                                                                 
853                                                                 // mostly copied from particle code
854                                                                 for(p=0, pa=psys->particles; p<psys->totpart; p++, pa++)
855                                                                 {
856                                                                         int cell[3];
857                                                                         size_t i = 0;
858                                                                         size_t index = 0;
859                                                                         int badcell = 0;
860                                                                         
861                                                                         if(pa->alive == PARS_KILLED) continue;
862                                                                         else if(pa->alive == PARS_UNBORN && (part->flag & PART_UNBORN)==0) continue;
863                                                                         else if(pa->alive == PARS_DEAD && (part->flag & PART_DIED)==0) continue;
864                                                                         else if(pa->flag & (PARS_UNEXIST+PARS_NO_DISP)) continue;
865                                                                         
866                                                                         // VECCOPY(pos, pa->state.co);
867                                                                         // Mat4MulVecfl (ob->imat, pos);
868                                                                         
869                                                                         // 1. get corresponding cell
870                                                                         get_cell(smd, pa->state.co, cell, 0);
871                                                                 
872                                                                         // check if cell is valid (in the domain boundary)
873                                                                         for(i = 0; i < 3; i++)
874                                                                         {
875                                                                                 if((cell[i] > sds->res[i] - 1) || (cell[i] < 0))
876                                                                                 {
877                                                                                         badcell = 1;
878                                                                                         break;
879                                                                                 }
880                                                                         }
881                                                                                 
882                                                                         if(badcell)
883                                                                                 continue;
884                                                                         
885                                                                         // 2. set cell values (heat, density and velocity)
886                                                                         index = smoke_get_index(cell[0], sds->res[0], cell[1], sds->res[1], cell[2]);
887                                                                         
888                                                                         if(!(sfs->type & MOD_SMOKE_FLOW_TYPE_OUTFLOW)) // this is inflow
889                                                                         {
890                                                                                 heat[index] = sfs->temp;
891                                                                                 density[index] = sfs->density;
892                                                                                 /*
893                                                                                 velocity_x[index] = pa->state.vel[0];
894                                                                                 velocity_y[index] = pa->state.vel[1];
895                                                                                 velocity_z[index] = pa->state.vel[2];
896                                                                                 */
897
898                                                                                 // we need different handling for the high-res feature
899                                                                                 if(bigdensity)
900                                                                                 {
901                                                                                         // init all surrounding cells according to amplification, too
902                                                                                         int i, j, k;
903
904                                                                                         smoke_turbulence_get_res(smd->domain->wt, bigres);
905
906                                                                                         for(i = 0; i < smd->domain->amplify + 1; i++)
907                                                                                                 for(j = 0; j < smd->domain->amplify + 1; j++)
908                                                                                                         for(k = 0; k < smd->domain->amplify + 1; k++)
909                                                                                                         {
910                                                                                                                 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);
911                                                                                                                 bigdensity[index] = sfs->density;
912                                                                                                         }
913                                                                                 }
914                                                                         }
915                                                                         else // outflow
916                                                                         {
917                                                                                 heat[index] = 0.f;
918                                                                                 density[index] = 0.f;
919                                                                                 velocity_x[index] = 0.f;
920                                                                                 velocity_y[index] = 0.f;
921                                                                                 velocity_z[index] = 0.f;
922
923                                                                                 // we need different handling for the high-res feature
924                                                                                 if(bigdensity)
925                                                                                 {
926                                                                                         // init all surrounding cells according to amplification, too
927                                                                                         int i, j, k;
928
929                                                                                         smoke_turbulence_get_res(smd->domain->wt, bigres);
930
931                                                                                         for(i = 0; i < smd->domain->amplify + 1; i++)
932                                                                                                 for(j = 0; j < smd->domain->amplify + 1; j++)
933                                                                                                         for(k = 0; k < smd->domain->amplify + 1; k++)
934                                                                                                         {
935                                                                                                                 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);
936                                                                                                                 bigdensity[index] = 0.f;
937                                                                                                         }
938                                                                                 }
939                                                                         }
940                                                                 }
941                                                         }       
942                                                 }       
943                                         }
944
945                                         if(cnt_domain && !sds->fluid_group)
946                                                 go = go->next;
947                                         else
948                                                 base= base->next;
949                                 }
950                         }
951
952                         // do effectors
953                         /*
954                         if(sds->eff_group)
955                         {
956                                 for(go = sds->eff_group->gobject.first; go; go = go->next) 
957                                 {
958                                         if(go->ob)
959                                         {
960                                                 if(ob->pd)
961                                                 {
962                                                         
963                                                 }
964                                         }
965                                 }
966                         }
967                         */
968
969                         // do collisions        
970                         if(sds->coll_group || !cnt_domain)
971                         {
972                                 Object *otherobj = NULL;
973                                 ModifierData *md = NULL;
974
975                                 if(cnt_domain && !sds->coll_group) // we use groups since we have 2 domains
976                                         go = sds->coll_group->gobject.first;
977                                 else
978                                         base = scene->base.first;
979
980                                 while(base || go)
981                                 {
982                                         otherobj = NULL;
983
984                                         if(cnt_domain && !sds->coll_group) 
985                                         {
986                                                 if(go->ob)
987                                                         otherobj = go->ob;
988                                         }
989                                         else
990                                                 otherobj = base->object;
991
992                                         if(!otherobj)
993                                         {
994                                                 if(cnt_domain && !sds->coll_group)
995                                                         go = go->next;
996                                                 else
997                                                         base= base->next;
998
999                                                 continue;
1000                                         }
1001                         
1002                                         md = modifiers_findByType(otherobj, eModifierType_Smoke);
1003                                         
1004                                         // check for active smoke modifier
1005                                         if(md && md->mode & (eModifierMode_Realtime | eModifierMode_Render))
1006                                         {
1007                                                 SmokeModifierData *smd2 = (SmokeModifierData *)md;
1008
1009                                                 if((smd2->type & MOD_SMOKE_TYPE_COLL) && smd2->coll)
1010                                                 {
1011                                                         // we got nice collision object
1012                                                         SmokeCollSettings *scs = smd2->coll;
1013                                                         size_t i, j;
1014                                                         unsigned char *obstacles = smoke_get_obstacle(smd->domain->fluid);
1015
1016                                                         for(i = 0; i < scs->numpoints; i++)
1017                                                         {
1018                                                                 int badcell = 0;
1019                                                                 size_t index = 0;
1020                                                                 int cell[3];
1021
1022                                                                 // 1. get corresponding cell
1023                                                                 get_cell(smd, &scs->points[3 * i], cell, 0);
1024                                                         
1025                                                                 // check if cell is valid (in the domain boundary)
1026                                                                 for(j = 0; j < 3; j++)
1027                                                                         if((cell[j] > sds->res[j] - 1) || (cell[j] < 0))
1028                                                                         {
1029                                                                                 badcell = 1;
1030                                                                                 break;
1031                                                                         }
1032                                                                                 
1033                                                                 if(badcell)
1034                                                                         continue;
1035
1036                                                                 // 2. set cell values (heat, density and velocity)
1037                                                                 index = smoke_get_index(cell[0], sds->res[0], cell[1], sds->res[1], cell[2]);
1038                                                                 
1039                                                                 // printf("cell[0]: %d, cell[1]: %d, cell[2]: %d\n", cell[0], cell[1], cell[2]);
1040                                                                 // printf("res[0]: %d, res[1]: %d, res[2]: %d, index: %d\n\n", sds->res[0], sds->res[1], sds->res[2], index);
1041                                                                         
1042                                                                 obstacles[index] = 1;
1043
1044                                                                 // for moving gobstacles
1045                                                                 /*
1046                                                                 const LbmFloat maxVelVal = 0.1666;
1047                                                                 const LbmFloat maxusqr = maxVelVal*maxVelVal*3. *1.5;
1048
1049                                                                 LbmVec objvel = vec2L((mMOIVertices[n]-mMOIVerticesOld[n]) /dvec); { 
1050                                                                 const LbmFloat usqr = (objvel[0]*objvel[0]+objvel[1]*objvel[1]+objvel[2]*objvel[2])*1.5; 
1051                                                                 USQRMAXCHECK(usqr, objvel[0],objvel[1],objvel[2], mMaxVlen, mMxvx,mMxvy,mMxvz); 
1052                                                                 if(usqr>maxusqr) { 
1053                                                                         // cutoff at maxVelVal 
1054                                                                         for(int jj=0; jj<3; jj++) { 
1055                                                                                 if(objvel[jj]>0.) objvel[jj] =  maxVelVal;  
1056                                                                                 if(objvel[jj]<0.) objvel[jj] = -maxVelVal; 
1057                                                                         } 
1058                                                                 } } 
1059
1060                                                                 const LbmFloat dp=dot(objvel, vec2L((*pNormals)[n]) ); 
1061                                                                 const LbmVec oldov=objvel; // debug
1062                                                                 objvel = vec2L((*pNormals)[n]) *dp;
1063                                                                 */
1064                                                         }
1065                                                 }
1066                                         }
1067
1068                                         if(cnt_domain && !sds->coll_group)
1069                                                 go = go->next;
1070                                         else
1071                                                 base= base->next;
1072                                 }
1073                         }
1074                         
1075                         // set new time
1076                         smd->time = scene->r.cfra;
1077
1078                         // simulate the actual smoke (c++ code in intern/smoke)
1079                         smoke_step(sds->fluid);
1080                         if(sds->wt)
1081                                 smoke_turbulence_step(sds->wt, sds->fluid);
1082
1083                         tend();
1084                         printf ( "Frame: %d, Time: %f\n", (int)smd->time, ( float ) tval() );
1085                 }
1086                 else if(scene->r.cfra < smd->time)
1087                 {
1088                         // we got back in time, reset smoke in this case (TODO: use cache later)
1089                         smd->time = scene->r.cfra;
1090                         smokeModifier_reset(smd);
1091                 }
1092         }
1093 }
1094
1095 // update necessary information for 3dview
1096 void smoke_prepare_View(SmokeModifierData *smd, float *light)
1097 {
1098         float *density = NULL;
1099         int x, y, z;
1100
1101         if(!smd->domain->tray)
1102         {
1103                 // TRay is for self shadowing
1104                 smd->domain->tray = MEM_callocN(sizeof(float)*smd->domain->res[0]*smd->domain->res[1]*smd->domain->res[2], "Smoke_tRay");
1105         }
1106         if(!smd->domain->tvox)
1107         {
1108                 // TVox is for tranaparency
1109                 smd->domain->tvox = MEM_callocN(sizeof(float)*smd->domain->res[0]*smd->domain->res[1]*smd->domain->res[2], "Smoke_tVox");
1110         }
1111
1112         // update 3dview
1113         density = smoke_get_density(smd->domain->fluid);
1114         for(x = 0; x < smd->domain->res[0]; x++)
1115                         for(y = 0; y < smd->domain->res[1]; y++)
1116                                 for(z = 0; z < smd->domain->res[2]; z++)
1117                                 {
1118                                         size_t index;
1119
1120                                         index = smoke_get_index(x, smd->domain->res[0], y, smd->domain->res[1], z);
1121                                         // Transparency computation
1122                                         // formula taken from "Visual Simulation of Smoke" / Fedkiw et al. pg. 4
1123                                         // T_vox = exp(-C_ext * h)
1124                                         // C_ext/sigma_t = density * C_ext
1125                                         smoke_set_tvox(smd, index, exp(-density[index] * 7.0 * smd->domain->dx));
1126         }
1127         smoke_calc_transparency(smd, light, 0);
1128 }
1129
1130 // update necessary information for 3dview ("high res" option)
1131 void smoke_prepare_bigView(SmokeModifierData *smd, float *light)
1132 {
1133         float *density = NULL;
1134         size_t i = 0;
1135         int bigres[3];
1136
1137         smoke_turbulence_get_res(smd->domain->wt, bigres);
1138
1139         if(!smd->domain->traybig)
1140         {
1141                 // TRay is for self shadowing
1142                 smd->domain->traybig = MEM_callocN(sizeof(float)*bigres[0]*bigres[1]*bigres[2], "Smoke_tRayBig");
1143         }
1144         if(!smd->domain->tvoxbig)
1145         {
1146                 // TVox is for tranaparency
1147                 smd->domain->tvoxbig = MEM_callocN(sizeof(float)*bigres[0]*bigres[1]*bigres[2], "Smoke_tVoxBig");
1148         }
1149
1150         density = smoke_turbulence_get_density(smd->domain->wt);
1151         for (i = 0; i < bigres[0] * bigres[1] * bigres[2]; i++)
1152         {
1153                 // Transparency computation
1154                 // formula taken from "Visual Simulation of Smoke" / Fedkiw et al. pg. 4
1155                 // T_vox = exp(-C_ext * h)
1156                 // C_ext/sigma_t = density * C_ext
1157                 smoke_set_bigtvox(smd, i, exp(-density[i] * 7.0 * smd->domain->dx / (smd->domain->amplify + 1)) );
1158         }
1159         smoke_calc_transparency(smd, light, 1);
1160 }
1161
1162
1163 float smoke_get_tvox(SmokeModifierData *smd, size_t index)
1164 {
1165         return smd->domain->tvox[index];
1166 }
1167
1168 void smoke_set_tvox(SmokeModifierData *smd, size_t index, float tvox)
1169 {
1170         smd->domain->tvox[index] = tvox;
1171 }
1172
1173 float smoke_get_tray(SmokeModifierData *smd, size_t index)
1174 {
1175         return smd->domain->tray[index];
1176 }
1177
1178 void smoke_set_tray(SmokeModifierData *smd, size_t index, float transparency)
1179 {
1180         smd->domain->tray[index] = transparency;
1181 }
1182
1183 float smoke_get_bigtvox(SmokeModifierData *smd, size_t index)
1184 {
1185         return smd->domain->tvoxbig[index];
1186 }
1187
1188 void smoke_set_bigtvox(SmokeModifierData *smd, size_t index, float tvox)
1189 {
1190         smd->domain->tvoxbig[index] = tvox;
1191 }
1192
1193 float smoke_get_bigtray(SmokeModifierData *smd, size_t index)
1194 {
1195         return smd->domain->traybig[index];
1196 }
1197
1198 void smoke_set_bigtray(SmokeModifierData *smd, size_t index, float transparency)
1199 {
1200         smd->domain->traybig[index] = transparency;
1201 }
1202
1203 long long smoke_get_mem_req(int xres, int yres, int zres, int amplify)
1204 {
1205           int totalCells = xres * yres * zres;
1206           int amplifiedCells = totalCells * amplify * amplify * amplify;
1207
1208           // print out memory requirements
1209           long long int coarseSize = sizeof(float) * totalCells * 22 +
1210                            sizeof(unsigned char) * totalCells;
1211
1212           long long int fineSize = sizeof(float) * amplifiedCells * 7 + // big grids
1213                          sizeof(float) * totalCells * 8 +     // small grids
1214                          sizeof(float) * 128 * 128 * 128;     // noise tile
1215
1216           long long int totalMB = (coarseSize + fineSize) / (1024 * 1024);
1217
1218           return totalMB;
1219 }
1220
1221
1222 static void calc_voxel_transp(SmokeModifierData *smd, int *pixel, float *tRay)
1223 {
1224         // printf("Pixel(%d, %d, %d)\n", pixel[0], pixel[1], pixel[2]);
1225         const size_t index = smoke_get_index(pixel[0], smd->domain->res[0], pixel[1], smd->domain->res[1], pixel[2]);
1226
1227         // T_ray *= T_vox
1228         *tRay *= smoke_get_tvox(smd, index);
1229 }
1230
1231 static void calc_voxel_transp_big(SmokeModifierData *smd, int *pixel, float *tRay)
1232 {
1233         int bigres[3];
1234         size_t index;
1235
1236         smoke_turbulence_get_res(smd->domain->wt, bigres);
1237         index = smoke_get_index(pixel[0], bigres[0], pixel[1], bigres[1], pixel[2]);
1238
1239         /*
1240         if(index > bigres[0]*bigres[1]*bigres[2])
1241                 printf("pixel[0]: %d, [1]: %d, [2]: %d\n", pixel[0], pixel[1], pixel[2]);
1242         */
1243
1244         // T_ray *= T_vox
1245         *tRay *= smoke_get_bigtvox(smd, index);
1246 }
1247
1248 static void bresenham_linie_3D(SmokeModifierData *smd, int x1, int y1, int z1, int x2, int y2, int z2, float *tRay, int big)
1249 {
1250     int dx, dy, dz, i, l, m, n, x_inc, y_inc, z_inc, err_1, err_2, dx2, dy2, dz2;
1251     int pixel[3];
1252
1253     pixel[0] = x1;
1254     pixel[1] = y1;
1255     pixel[2] = z1;
1256
1257     dx = x2 - x1;
1258     dy = y2 - y1;
1259     dz = z2 - z1;
1260
1261     x_inc = (dx < 0) ? -1 : 1;
1262     l = abs(dx);
1263     y_inc = (dy < 0) ? -1 : 1;
1264     m = abs(dy);
1265     z_inc = (dz < 0) ? -1 : 1;
1266     n = abs(dz);
1267     dx2 = l << 1;
1268     dy2 = m << 1;
1269     dz2 = n << 1;
1270
1271     if ((l >= m) && (l >= n)) {
1272         err_1 = dy2 - l;
1273         err_2 = dz2 - l;
1274         for (i = 0; i < l; i++) {
1275                 if(!big)
1276                                 calc_voxel_transp(smd, pixel, tRay);
1277                         else
1278                                 calc_voxel_transp_big(smd, pixel, tRay);
1279                 if(*tRay < 0.0f)
1280                         return;
1281             if (err_1 > 0) {
1282                 pixel[1] += y_inc;
1283                 err_1 -= dx2;
1284             }
1285             if (err_2 > 0) {
1286                 pixel[2] += z_inc;
1287                 err_2 -= dx2;
1288             }
1289             err_1 += dy2;
1290             err_2 += dz2;
1291             pixel[0] += x_inc;
1292         }
1293     } else if ((m >= l) && (m >= n)) {
1294         err_1 = dx2 - m;
1295         err_2 = dz2 - m;
1296         for (i = 0; i < m; i++) {
1297                 if(!big)
1298                                 calc_voxel_transp(smd, pixel, tRay);
1299                         else
1300                                 calc_voxel_transp_big(smd, pixel, tRay);
1301                 if(*tRay < 0.0f)
1302                         return;
1303             if (err_1 > 0) {
1304                 pixel[0] += x_inc;
1305                 err_1 -= dy2;
1306             }
1307             if (err_2 > 0) {
1308                 pixel[2] += z_inc;
1309                 err_2 -= dy2;
1310             }
1311             err_1 += dx2;
1312             err_2 += dz2;
1313             pixel[1] += y_inc;
1314         }
1315     } else {
1316         err_1 = dy2 - n;
1317         err_2 = dx2 - n;
1318         for (i = 0; i < n; i++) {
1319                 if(!big)
1320                                 calc_voxel_transp(smd, pixel, tRay);
1321                         else
1322                                 calc_voxel_transp_big(smd, pixel, tRay);
1323                 if(*tRay < 0.0f)
1324                         return;
1325             if (err_1 > 0) {
1326                 pixel[1] += y_inc;
1327                 err_1 -= dz2;
1328             }
1329             if (err_2 > 0) {
1330                 pixel[0] += x_inc;
1331                 err_2 -= dz2;
1332             }
1333             err_1 += dy2;
1334             err_2 += dx2;
1335             pixel[2] += z_inc;
1336         }
1337     }
1338     if(!big)
1339         calc_voxel_transp(smd, pixel, tRay);
1340     else
1341         calc_voxel_transp_big(smd, pixel, tRay);
1342 }
1343
1344 static void get_cell(struct SmokeModifierData *smd, float *pos, int *cell, int correct)
1345 {
1346         float tmp[3];
1347
1348         VECSUB(tmp, pos, smd->domain->p0);
1349         VecMulf(tmp, 1.0 / smd->domain->dx);
1350
1351         if(correct)
1352         {
1353                 cell[0] = MIN2(smd->domain->res[0] - 1, MAX2(0, (int)floor(tmp[0])));
1354                 cell[1] = MIN2(smd->domain->res[1] - 1, MAX2(0, (int)floor(tmp[1])));
1355                 cell[2] = MIN2(smd->domain->res[2] - 1, MAX2(0, (int)floor(tmp[2])));
1356         }
1357         else
1358         {
1359                 cell[0] = (int)floor(tmp[0]);
1360                 cell[1] = (int)floor(tmp[1]);
1361                 cell[2] = (int)floor(tmp[2]);
1362         }
1363 }
1364 static void get_bigcell(struct SmokeModifierData *smd, float *pos, int *cell, int correct)
1365 {
1366         float tmp[3];
1367         int res[3];
1368         smoke_turbulence_get_res(smd->domain->wt, res);
1369
1370         VECSUB(tmp, pos, smd->domain->p0);
1371
1372         VecMulf(tmp, (smd->domain->amplify + 1)/ smd->domain->dx );
1373
1374         if(correct)
1375         {
1376                 cell[0] = MIN2(res[0] - 1, MAX2(0, (int)floor(tmp[0])));
1377                 cell[1] = MIN2(res[1] - 1, MAX2(0, (int)floor(tmp[1])));
1378                 cell[2] = MIN2(res[2] - 1, MAX2(0, (int)floor(tmp[2])));
1379         }
1380         else
1381         {
1382                 cell[0] = (int)floor(tmp[0]);
1383                 cell[1] = (int)floor(tmp[1]);
1384                 cell[2] = (int)floor(tmp[2]);
1385         }
1386 }
1387
1388
1389 void smoke_calc_transparency(struct SmokeModifierData *smd, float *light, int big)
1390 {
1391         int x, y, z;
1392         float bv[6];
1393         int res[3];
1394         float bigfactor = 1.0;
1395
1396         // x
1397         bv[0] = smd->domain->p0[0];
1398         bv[1] = smd->domain->p1[0];
1399         // y
1400         bv[2] = smd->domain->p0[1];
1401         bv[3] = smd->domain->p1[1];
1402         // z
1403         bv[4] = smd->domain->p0[2];
1404         bv[5] = smd->domain->p1[2];
1405 /*
1406         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]);
1407
1408         printf("p0[0]: %f, p0[1]: %f, p0[2]: %f\n", smd->domain->p0[0], smd->domain->p0[1], smd->domain->p0[2]);
1409         printf("p1[0]: %f, p1[1]: %f, p1[2]: %f\n", smd->domain->p1[0], smd->domain->p1[1], smd->domain->p1[2]);
1410         printf("dx: %f, amp: %d\n", smd->domain->dx, smd->domain->amplify);
1411 */
1412         if(!big)
1413         {
1414                 res[0] = smd->domain->res[0];
1415                 res[1] = smd->domain->res[1];
1416                 res[2] = smd->domain->res[2];
1417         }
1418         else
1419         {
1420                 smoke_turbulence_get_res(smd->domain->wt, res);
1421                 bigfactor = 1.0 / (smd->domain->amplify + 1);
1422         }
1423
1424 #pragma omp parallel for schedule(static) private(y, z) shared(big, smd, light, res, bigfactor)
1425         for(x = 0; x < res[0]; x++)
1426                 for(y = 0; y < res[1]; y++)
1427                         for(z = 0; z < res[2]; z++)
1428                         {
1429                                 float voxelCenter[3];
1430                                 size_t index;
1431                                 float pos[3];
1432                                 int cell[3];
1433                                 float tRay = 1.0;
1434
1435                                 index = smoke_get_index(x, res[0], y, res[1], z);
1436
1437                                 // voxelCenter = m_voxelarray[i].GetCenter();
1438                                 voxelCenter[0] = smd->domain->p0[0] + smd->domain->dx * bigfactor * x + smd->domain->dx * bigfactor * 0.5;
1439                                 voxelCenter[1] = smd->domain->p0[1] + smd->domain->dx * bigfactor * y + smd->domain->dx * bigfactor * 0.5;
1440                                 voxelCenter[2] = smd->domain->p0[2] + smd->domain->dx * bigfactor * z + smd->domain->dx * bigfactor * 0.5;
1441
1442                                 // printf("vc[0]: %f, vc[1]: %f, vc[2]: %f\n", voxelCenter[0], voxelCenter[1], voxelCenter[2]);
1443                                 // printf("light[0]: %f, light[1]: %f, light[2]: %f\n", light[0], light[1], light[2]);
1444
1445                                 // get starting position (in voxel coords)
1446                                 if(BLI_bvhtree_bb_raycast(bv, light, voxelCenter, pos) > FLT_EPSILON)
1447                                 {
1448                                         // we're ouside
1449                                         // printf("out: pos[0]: %f, pos[1]: %f, pos[2]: %f\n", pos[0], pos[1], pos[2]);
1450                                         if(!big)
1451                                                 get_cell(smd, pos, cell, 1);
1452                                         else
1453                                                 get_bigcell(smd, pos, cell, 1);
1454                                 }
1455                                 else
1456                                 {
1457                                         // printf("in: pos[0]: %f, pos[1]: %f, pos[2]: %f\n", light[0], light[1], light[2]);
1458                                         // we're inside
1459                                         if(!big)
1460                                                 get_cell(smd, light, cell, 1);
1461                                         else
1462                                                 get_bigcell(smd, light, cell, 1);
1463                                 }
1464
1465                                 // printf("cell - [0]: %d, [1]: %d, [2]: %d\n", cell[0], cell[1], cell[2]);
1466                                 bresenham_linie_3D(smd, cell[0], cell[1], cell[2], x, y, z, &tRay, big);
1467
1468                                 if(!big)
1469                                         smoke_set_tray(smd, index, tRay);
1470                                 else
1471                                         smoke_set_bigtray(smd, index, tRay);
1472                         }
1473 }
1474