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