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