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