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