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