61d9ea99c2f25f14b17bc2aa30b2a59d052dc902
[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_threads.h"
54 #include "BLI_utildefines.h"
55 #include "BLI_voxel.h"
56
57 #include "DNA_anim_types.h"
58 #include "DNA_armature_types.h"
59 #include "DNA_constraint_types.h"
60 #include "DNA_customdata_types.h"
61 #include "DNA_group_types.h"
62 #include "DNA_lamp_types.h"
63 #include "DNA_mesh_types.h"
64 #include "DNA_meshdata_types.h"
65 #include "DNA_modifier_types.h"
66 #include "DNA_object_types.h"
67 #include "DNA_particle_types.h"
68 #include "DNA_scene_types.h"
69 #include "DNA_smoke_types.h"
70
71 #include "BKE_animsys.h"
72 #include "BKE_armature.h"
73 #include "BKE_bvhutils.h"
74 #include "BKE_cdderivedmesh.h"
75 #include "BKE_collision.h"
76 #include "BKE_constraint.h"
77 #include "BKE_customdata.h"
78 #include "BKE_deform.h"
79 #include "BKE_DerivedMesh.h"
80 #include "BKE_global.h"
81 #include "BKE_effect.h"
82 #include "BKE_main.h"
83 #include "BKE_modifier.h"
84 #include "BKE_object.h"
85 #include "BKE_particle.h"
86 #include "BKE_pointcache.h"
87 #include "BKE_scene.h"
88 #include "BKE_smoke.h"
89 #include "BKE_texture.h"
90
91 #include "RE_shader_ext.h"
92
93 /* UNUSED so far, may be enabled later */
94 /* #define USE_SMOKE_COLLISION_DM */
95
96 #include "smoke_API.h"
97
98 #ifdef WITH_SMOKE
99
100 #ifdef _WIN32
101 #include <time.h>
102 #include <stdio.h>
103 #include <conio.h>
104 #include <windows.h>
105
106 static LARGE_INTEGER liFrequency;
107 static LARGE_INTEGER liStartTime;
108 static LARGE_INTEGER liCurrentTime;
109
110 static void tstart(void)
111 {
112         QueryPerformanceFrequency(&liFrequency);
113         QueryPerformanceCounter(&liStartTime);
114 }
115 static void tend(void)
116 {
117         QueryPerformanceCounter(&liCurrentTime);
118 }
119 static double UNUSED_FUNCTION(tval) (void)
120 {
121         return ((double)( (liCurrentTime.QuadPart - liStartTime.QuadPart) * (double)1000.0 / (double)liFrequency.QuadPart));
122 }
123 #else
124 #include <sys/time.h>
125 static struct timeval _tstart, _tend;
126 static struct timezone tz;
127 static void tstart(void)
128 {
129         gettimeofday(&_tstart, &tz);
130 }
131 static void tend(void)
132 {
133         gettimeofday(&_tend, &tz);
134 }
135
136 static double UNUSED_FUNCTION(tval) (void)
137 {
138         double t1, t2;
139         t1 = ( double ) _tstart.tv_sec * 1000 + ( double ) _tstart.tv_usec / (1000);
140         t2 = ( double ) _tend.tv_sec * 1000 + ( double ) _tend.tv_usec / (1000);
141         return t2 - t1;
142 }
143 #endif
144
145 struct Object;
146 struct Scene;
147 struct DerivedMesh;
148 struct SmokeModifierData;
149
150 // timestep default value for nice appearance 0.1f
151 #define DT_DEFAULT 0.1f
152
153 #define ADD_IF_LOWER_POS(a, b) (MIN2((a) + (b), MAX2((a), (b))))
154 #define ADD_IF_LOWER_NEG(a, b) (MAX2((a) + (b), MIN2((a), (b))))
155 #define ADD_IF_LOWER(a, b) (((b) > 0) ? ADD_IF_LOWER_POS((a), (b)) : ADD_IF_LOWER_NEG((a), (b)))
156
157 #else /* WITH_SMOKE */
158
159 /* Stubs to use when smoke is disabled */
160 struct WTURBULENCE *smoke_turbulence_init(int *UNUSED(res), int UNUSED(amplify), int UNUSED(noisetype), const char *UNUSED(noisefile_path), int UNUSED(use_fire), int UNUSED(use_colors)) { return NULL; }
161 //struct FLUID_3D *smoke_init(int *UNUSED(res), float *UNUSED(dx), float *UNUSED(dtdef), int UNUSED(use_heat), int UNUSED(use_fire), int UNUSED(use_colors)) { return NULL; }
162 void smoke_free(struct FLUID_3D *UNUSED(fluid)) {}
163 float *smoke_get_density(struct FLUID_3D *UNUSED(fluid)) { return NULL; }
164 void smoke_turbulence_free(struct WTURBULENCE *UNUSED(wt)) {}
165 void smoke_initWaveletBlenderRNA(struct WTURBULENCE *UNUSED(wt), float *UNUSED(strength)) {}
166 void smoke_initBlenderRNA(struct FLUID_3D *UNUSED(fluid), float *UNUSED(alpha), float *UNUSED(beta), float *UNUSED(dt_factor), float *UNUSED(vorticity),
167                           int *UNUSED(border_colli), float *UNUSED(burning_rate), float *UNUSED(flame_smoke), float *UNUSED(flame_smoke_color),
168                           float *UNUSED(flame_vorticity), float *UNUSED(flame_ignition_temp), float *UNUSED(flame_max_temp)) {}
169 struct DerivedMesh *smokeModifier_do(SmokeModifierData *UNUSED(smd), Scene *UNUSED(scene), Object *UNUSED(ob), DerivedMesh *UNUSED(dm), bool UNUSED(for_render)) { return NULL; }
170 float smoke_get_velocity_at(struct Object *UNUSED(ob), float UNUSED(position[3]), float UNUSED(velocity[3])) { return 0.0f; }
171 void flame_get_spectrum(unsigned char *UNUSED(spec), int UNUSED(width), float UNUSED(t1), float UNUSED(t2)) {}
172
173 #endif /* WITH_SMOKE */
174
175 #ifdef WITH_SMOKE
176
177 void smoke_reallocate_fluid(SmokeDomainSettings *sds, float dx, int res[3], int free_old)
178 {
179         int use_heat = (sds->active_fields & SM_ACTIVE_HEAT);
180         int use_fire = (sds->active_fields & SM_ACTIVE_FIRE);
181         int use_colors = (sds->active_fields & SM_ACTIVE_COLORS);
182
183         if (free_old && sds->fluid)
184                 smoke_free(sds->fluid);
185         if (!min_iii(res[0], res[1], res[2])) {
186                 sds->fluid = NULL;
187                 return;
188         }
189         sds->fluid = smoke_init(res, dx, DT_DEFAULT, use_heat, use_fire, use_colors);
190         smoke_initBlenderRNA(sds->fluid, &(sds->alpha), &(sds->beta), &(sds->time_scale), &(sds->vorticity), &(sds->border_collisions),
191                              &(sds->burning_rate), &(sds->flame_smoke), sds->flame_smoke_color, &(sds->flame_vorticity), &(sds->flame_ignition), &(sds->flame_max_temp));
192
193         /* reallocate shadow buffer */
194         if (sds->shadow)
195                 MEM_freeN(sds->shadow);
196         sds->shadow = MEM_callocN(sizeof(float) * res[0] * res[1] * res[2], "SmokeDomainShadow");
197 }
198
199 void smoke_reallocate_highres_fluid(SmokeDomainSettings *sds, float dx, int res[3], int free_old)
200 {
201         int use_fire = (sds->active_fields & (SM_ACTIVE_HEAT | SM_ACTIVE_FIRE));
202         int use_colors = (sds->active_fields & SM_ACTIVE_COLORS);
203
204         if (free_old && sds->wt)
205                 smoke_turbulence_free(sds->wt);
206         if (!min_iii(res[0], res[1], res[2])) {
207                 sds->wt = NULL;
208                 return;
209         }
210         sds->wt = smoke_turbulence_init(res, sds->amplify + 1, sds->noise, BLI_temporary_dir(), use_fire, use_colors);
211         sds->res_wt[0] = res[0] * (sds->amplify + 1);
212         sds->res_wt[1] = res[1] * (sds->amplify + 1);
213         sds->res_wt[2] = res[2] * (sds->amplify + 1);
214         sds->dx_wt = dx / (sds->amplify + 1);
215         smoke_initWaveletBlenderRNA(sds->wt, &(sds->strength));
216 }
217
218 /* convert global position to domain cell space */
219 static void smoke_pos_to_cell(SmokeDomainSettings *sds, float pos[3])
220 {
221         mul_m4_v3(sds->imat, pos);
222         sub_v3_v3(pos, sds->p0);
223         pos[0] *= 1.0f / sds->cell_size[0];
224         pos[1] *= 1.0f / sds->cell_size[1];
225         pos[2] *= 1.0f / sds->cell_size[2];
226 }
227
228 /* set domain transformations and base resolution from object derivedmesh */
229 static void smoke_set_domain_from_derivedmesh(SmokeDomainSettings *sds, Object *ob, DerivedMesh *dm, int init_resolution)
230 {
231         size_t i;
232         float min[3] = {FLT_MAX, FLT_MAX, FLT_MAX}, max[3] = {-FLT_MAX, -FLT_MAX, -FLT_MAX};
233         float size[3];
234         MVert *verts = dm->getVertArray(dm);
235         float scale = 0.0;
236         int res;
237
238         res = sds->maxres;
239
240         // get BB of domain
241         for (i = 0; i < dm->getNumVerts(dm); i++)
242         {
243                 // min BB
244                 min[0] = MIN2(min[0], verts[i].co[0]);
245                 min[1] = MIN2(min[1], verts[i].co[1]);
246                 min[2] = MIN2(min[2], verts[i].co[2]);
247
248                 // max BB
249                 max[0] = MAX2(max[0], verts[i].co[0]);
250                 max[1] = MAX2(max[1], verts[i].co[1]);
251                 max[2] = MAX2(max[2], verts[i].co[2]);
252         }
253
254         /* set domain bounds */
255         copy_v3_v3(sds->p0, min);
256         copy_v3_v3(sds->p1, max);
257         sds->dx = 1.0f / res;
258
259         /* calculate domain dimensions */
260         sub_v3_v3v3(size, max, min);
261         if (init_resolution) {
262                 zero_v3_int(sds->base_res);
263                 copy_v3_v3(sds->cell_size, size);
264         }
265         /* apply object scale */
266         for (i = 0; i < 3; i++) {
267                 size[i] = fabs(size[i] * ob->size[i]);
268         }
269         copy_v3_v3(sds->global_size, size);
270         copy_v3_v3(sds->dp0, min);
271
272         invert_m4_m4(sds->imat, ob->obmat);
273
274         // prevent crash when initializing a plane as domain
275         if (!init_resolution || (size[0] < FLT_EPSILON) || (size[1] < FLT_EPSILON) || (size[2] < FLT_EPSILON))
276                 return;
277
278         /* define grid resolutions from longest domain side */
279         if (size[0] >= MAX2(size[1], size[2])) {
280                 scale = res / size[0];
281                 sds->scale = size[0] / fabsf(ob->size[0]);
282                 sds->base_res[0] = res;
283                 sds->base_res[1] = (int)(size[1] * scale + 0.5f);
284                 sds->base_res[2] = (int)(size[2] * scale + 0.5f);
285         }
286         else if (size[1] >= MAX2(size[0], size[2])) {
287                 scale = res / size[1];
288                 sds->scale = size[1] / fabsf(ob->size[1]);
289                 sds->base_res[0] = (int)(size[0] * scale + 0.5f);
290                 sds->base_res[1] = res;
291                 sds->base_res[2] = (int)(size[2] * scale + 0.5f);
292         }
293         else {
294                 scale = res / size[2];
295                 sds->scale = size[2] / fabsf(ob->size[2]);
296                 sds->base_res[0] = (int)(size[0] * scale + 0.5f);
297                 sds->base_res[1] = (int)(size[1] * scale + 0.5f);
298                 sds->base_res[2] = res;
299         }
300
301         /* set cell size */
302         sds->cell_size[0] /= (float)sds->base_res[0];
303         sds->cell_size[1] /= (float)sds->base_res[1];
304         sds->cell_size[2] /= (float)sds->base_res[2];
305 }
306
307 static int smokeModifier_init(SmokeModifierData *smd, Object *ob, Scene *scene, DerivedMesh *dm)
308 {
309         if ((smd->type & MOD_SMOKE_TYPE_DOMAIN) && smd->domain && !smd->domain->fluid)
310         {
311                 SmokeDomainSettings *sds = smd->domain;
312                 int res[3];
313                 /* set domain dimensions from derivedmesh */
314                 smoke_set_domain_from_derivedmesh(sds, ob, dm, TRUE);
315                 /* reset domain values */
316                 zero_v3_int(sds->shift);
317                 zero_v3(sds->shift_f);
318                 add_v3_fl(sds->shift_f, 0.5f);
319                 zero_v3(sds->prev_loc);
320                 mul_m4_v3(ob->obmat, sds->prev_loc);
321
322                 /* set resolutions */
323                 if (smd->domain->flags & MOD_SMOKE_ADAPTIVE_DOMAIN) {
324                         res[0] = res[1] = res[2] = 1; /* use minimum res for adaptive init */
325                 }
326                 else {
327                         VECCOPY(res, sds->base_res);
328                 }
329                 VECCOPY(sds->res, res);
330                 sds->total_cells = sds->res[0] * sds->res[1] * sds->res[2];
331                 sds->res_min[0] = sds->res_min[1] = sds->res_min[2] = 0;
332                 VECCOPY(sds->res_max, res);
333
334                 /* allocate fluid */
335                 smoke_reallocate_fluid(sds, sds->dx, sds->res, 0);
336
337                 smd->time = scene->r.cfra;
338
339                 /* allocate highres fluid */
340                 if (sds->flags & MOD_SMOKE_HIGHRES) {
341                         smoke_reallocate_highres_fluid(sds, sds->dx, sds->res, 0);
342                 }
343                 /* allocate shadow buffer */
344                 if (!sds->shadow)
345                         sds->shadow = MEM_callocN(sizeof(float) * sds->res[0] * sds->res[1] * sds->res[2], "SmokeDomainShadow");
346
347                 return 1;
348         }
349         else if ((smd->type & MOD_SMOKE_TYPE_FLOW) && smd->flow)
350         {
351                 smd->time = scene->r.cfra;
352
353                 return 1;
354         }
355         else if ((smd->type & MOD_SMOKE_TYPE_COLL))
356         {
357                 if (!smd->coll)
358                 {
359                         smokeModifier_createType(smd);
360                 }
361
362                 smd->time = scene->r.cfra;
363
364                 return 1;
365         }
366
367         return 2;
368 }
369
370 #endif /* WITH_SMOKE */
371
372 static void smokeModifier_freeDomain(SmokeModifierData *smd)
373 {
374         if (smd->domain)
375         {
376                 if (smd->domain->shadow)
377                         MEM_freeN(smd->domain->shadow);
378                 smd->domain->shadow = NULL;
379
380                 if (smd->domain->fluid)
381                         smoke_free(smd->domain->fluid);
382
383                 if (smd->domain->fluid_mutex)
384                         BLI_rw_mutex_free(smd->domain->fluid_mutex);
385
386                 if (smd->domain->wt)
387                         smoke_turbulence_free(smd->domain->wt);
388
389                 if (smd->domain->effector_weights)
390                         MEM_freeN(smd->domain->effector_weights);
391                 smd->domain->effector_weights = NULL;
392
393                 BKE_ptcache_free_list(&(smd->domain->ptcaches[0]));
394                 smd->domain->point_cache[0] = NULL;
395
396                 MEM_freeN(smd->domain);
397                 smd->domain = NULL;
398         }
399 }
400
401 static void smokeModifier_freeFlow(SmokeModifierData *smd)
402 {
403         if (smd->flow)
404         {
405                 if (smd->flow->dm) smd->flow->dm->release(smd->flow->dm);
406                 if (smd->flow->verts_old) MEM_freeN(smd->flow->verts_old);
407                 MEM_freeN(smd->flow);
408                 smd->flow = NULL;
409         }
410 }
411
412 static void smokeModifier_freeCollision(SmokeModifierData *smd)
413 {
414         if (smd->coll)
415         {
416                 SmokeCollSettings *scs = smd->coll;
417
418                 if (scs->numverts)
419                 {
420                         if (scs->verts_old)
421                         {
422                                 MEM_freeN(scs->verts_old);
423                                 scs->verts_old = NULL;
424                         }
425                 }
426
427                 if (smd->coll->dm)
428                         smd->coll->dm->release(smd->coll->dm);
429                 smd->coll->dm = NULL;
430
431                 MEM_freeN(smd->coll);
432                 smd->coll = NULL;
433         }
434 }
435
436 void smokeModifier_reset_turbulence(struct SmokeModifierData *smd)
437 {
438         if (smd && smd->domain && smd->domain->wt)
439         {
440                 smoke_turbulence_free(smd->domain->wt);
441                 smd->domain->wt = NULL;
442         }
443 }
444
445 static void smokeModifier_reset_ex(struct SmokeModifierData *smd, bool need_lock)
446 {
447         if (smd)
448         {
449                 if (smd->domain)
450                 {
451                         if (smd->domain->shadow)
452                                 MEM_freeN(smd->domain->shadow);
453                         smd->domain->shadow = NULL;
454
455                         if (smd->domain->fluid)
456                         {
457                                 if (need_lock)
458                                         BLI_rw_mutex_lock(smd->domain->fluid_mutex, THREAD_LOCK_WRITE);
459
460                                 smoke_free(smd->domain->fluid);
461                                 smd->domain->fluid = NULL;
462
463                                 if (need_lock)
464                                         BLI_rw_mutex_unlock(smd->domain->fluid_mutex);
465                         }
466
467                         smokeModifier_reset_turbulence(smd);
468
469                         smd->time = -1;
470                         smd->domain->total_cells = 0;
471                         smd->domain->active_fields = 0;
472                 }
473                 else if (smd->flow)
474                 {
475                         if (smd->flow->verts_old) MEM_freeN(smd->flow->verts_old);
476                         smd->flow->verts_old = NULL;
477                         smd->flow->numverts = 0;
478                 }
479                 else if (smd->coll)
480                 {
481                         SmokeCollSettings *scs = smd->coll;
482
483                         if (scs->numverts && scs->verts_old)
484                         {
485                                 MEM_freeN(scs->verts_old);
486                                 scs->verts_old = NULL;
487                         }
488                 }
489         }
490 }
491
492 void smokeModifier_reset(struct SmokeModifierData *smd)
493 {
494         smokeModifier_reset_ex(smd, true);
495 }
496
497 void smokeModifier_free(SmokeModifierData *smd)
498 {
499         if (smd)
500         {
501                 smokeModifier_freeDomain(smd);
502                 smokeModifier_freeFlow(smd);
503                 smokeModifier_freeCollision(smd);
504         }
505 }
506
507 void smokeModifier_createType(struct SmokeModifierData *smd)
508 {
509         if (smd)
510         {
511                 if (smd->type & MOD_SMOKE_TYPE_DOMAIN)
512                 {
513                         if (smd->domain)
514                                 smokeModifier_freeDomain(smd);
515
516                         smd->domain = MEM_callocN(sizeof(SmokeDomainSettings), "SmokeDomain");
517
518                         smd->domain->smd = smd;
519
520                         smd->domain->point_cache[0] = BKE_ptcache_add(&(smd->domain->ptcaches[0]));
521                         smd->domain->point_cache[0]->flag |= PTCACHE_DISK_CACHE;
522                         smd->domain->point_cache[0]->step = 1;
523
524                         /* Deprecated */
525                         smd->domain->point_cache[1] = NULL;
526                         BLI_listbase_clear(&smd->domain->ptcaches[1]);
527                         /* set some standard values */
528                         smd->domain->fluid = NULL;
529                         smd->domain->fluid_mutex = BLI_rw_mutex_alloc();
530                         smd->domain->wt = NULL;
531                         smd->domain->eff_group = NULL;
532                         smd->domain->fluid_group = NULL;
533                         smd->domain->coll_group = NULL;
534                         smd->domain->maxres = 32;
535                         smd->domain->amplify = 1;
536                         smd->domain->alpha = -0.001;
537                         smd->domain->beta = 0.1;
538                         smd->domain->time_scale = 1.0;
539                         smd->domain->vorticity = 2.0;
540                         smd->domain->border_collisions = SM_BORDER_OPEN; // open domain
541                         smd->domain->flags = MOD_SMOKE_DISSOLVE_LOG;
542                         smd->domain->highres_sampling = SM_HRES_FULLSAMPLE;
543                         smd->domain->strength = 2.0;
544                         smd->domain->noise = MOD_SMOKE_NOISEWAVE;
545                         smd->domain->diss_speed = 5;
546                         smd->domain->active_fields = 0;
547
548                         smd->domain->adapt_margin = 4;
549                         smd->domain->adapt_res = 0;
550                         smd->domain->adapt_threshold = 0.02f;
551
552                         smd->domain->burning_rate = 0.75f;
553                         smd->domain->flame_smoke = 1.0f;
554                         smd->domain->flame_vorticity = 0.5f;
555                         smd->domain->flame_ignition = 1.25f;
556                         smd->domain->flame_max_temp = 1.75f;
557                         /* color */
558                         smd->domain->flame_smoke_color[0] = 0.7f;
559                         smd->domain->flame_smoke_color[1] = 0.7f;
560                         smd->domain->flame_smoke_color[2] = 0.7f;
561
562                         smd->domain->viewsettings = MOD_SMOKE_VIEW_SHOWBIG;
563                         smd->domain->effector_weights = BKE_add_effector_weights(NULL);
564                 }
565                 else if (smd->type & MOD_SMOKE_TYPE_FLOW)
566                 {
567                         if (smd->flow)
568                                 smokeModifier_freeFlow(smd);
569
570                         smd->flow = MEM_callocN(sizeof(SmokeFlowSettings), "SmokeFlow");
571
572                         smd->flow->smd = smd;
573
574                         /* set some standard values */
575                         smd->flow->density = 1.0f;
576                         smd->flow->fuel_amount = 1.0f;
577                         smd->flow->temp = 1.0f;
578                         smd->flow->flags = MOD_SMOKE_FLOW_ABSOLUTE | MOD_SMOKE_FLOW_USE_PART_SIZE;
579                         smd->flow->vel_multi = 1.0f;
580                         smd->flow->volume_density = 0.0f;
581                         smd->flow->surface_distance = 1.5f;
582                         smd->flow->source = MOD_SMOKE_FLOW_SOURCE_MESH;
583                         smd->flow->texture_size = 1.0f;
584                         smd->flow->particle_size = 1.0f;
585                         smd->flow->subframes = 0;
586
587                         smd->flow->color[0] = 0.7f;
588                         smd->flow->color[1] = 0.7f;
589                         smd->flow->color[2] = 0.7f;
590
591                         smd->flow->dm = NULL;
592                         smd->flow->psys = NULL;
593
594                 }
595                 else if (smd->type & MOD_SMOKE_TYPE_COLL)
596                 {
597                         if (smd->coll)
598                                 smokeModifier_freeCollision(smd);
599
600                         smd->coll = MEM_callocN(sizeof(SmokeCollSettings), "SmokeColl");
601
602                         smd->coll->smd = smd;
603                         smd->coll->verts_old = NULL;
604                         smd->coll->numverts = 0;
605                         smd->coll->type = 0; // static obstacle
606                         smd->coll->dm = NULL;
607
608 #ifdef USE_SMOKE_COLLISION_DM
609                         smd->coll->dm = NULL;
610 #endif
611                 }
612         }
613 }
614
615 void smokeModifier_copy(struct SmokeModifierData *smd, struct SmokeModifierData *tsmd)
616 {
617         tsmd->type = smd->type;
618         tsmd->time = smd->time;
619
620         smokeModifier_createType(tsmd);
621
622         if (tsmd->domain) {
623                 tsmd->domain->fluid_group = smd->domain->fluid_group;
624                 tsmd->domain->coll_group = smd->domain->coll_group;
625
626                 tsmd->domain->adapt_margin = smd->domain->adapt_margin;
627                 tsmd->domain->adapt_res = smd->domain->adapt_res;
628                 tsmd->domain->adapt_threshold = smd->domain->adapt_threshold;
629
630                 tsmd->domain->alpha = smd->domain->alpha;
631                 tsmd->domain->beta = smd->domain->beta;
632                 tsmd->domain->amplify = smd->domain->amplify;
633                 tsmd->domain->maxres = smd->domain->maxres;
634                 tsmd->domain->flags = smd->domain->flags;
635                 tsmd->domain->highres_sampling = smd->domain->highres_sampling;
636                 tsmd->domain->viewsettings = smd->domain->viewsettings;
637                 tsmd->domain->noise = smd->domain->noise;
638                 tsmd->domain->diss_speed = smd->domain->diss_speed;
639                 tsmd->domain->strength = smd->domain->strength;
640
641                 tsmd->domain->border_collisions = smd->domain->border_collisions;
642                 tsmd->domain->vorticity = smd->domain->vorticity;
643                 tsmd->domain->time_scale = smd->domain->time_scale;
644
645                 tsmd->domain->burning_rate = smd->domain->burning_rate;
646                 tsmd->domain->flame_smoke = smd->domain->flame_smoke;
647                 tsmd->domain->flame_vorticity = smd->domain->flame_vorticity;
648                 tsmd->domain->flame_ignition = smd->domain->flame_ignition;
649                 tsmd->domain->flame_max_temp = smd->domain->flame_max_temp;
650                 copy_v3_v3(tsmd->domain->flame_smoke_color, smd->domain->flame_smoke_color);
651
652                 MEM_freeN(tsmd->domain->effector_weights);
653                 tsmd->domain->effector_weights = MEM_dupallocN(smd->domain->effector_weights);
654         }
655         else if (tsmd->flow) {
656                 tsmd->flow->psys = smd->flow->psys;
657                 tsmd->flow->noise_texture = smd->flow->noise_texture;
658
659                 tsmd->flow->vel_multi = smd->flow->vel_multi;
660                 tsmd->flow->vel_normal = smd->flow->vel_normal;
661                 tsmd->flow->vel_random = smd->flow->vel_random;
662
663                 tsmd->flow->density = smd->flow->density;
664                 copy_v3_v3(tsmd->flow->color, smd->flow->color);
665                 tsmd->flow->fuel_amount = smd->flow->fuel_amount;
666                 tsmd->flow->temp = smd->flow->temp;
667                 tsmd->flow->volume_density = smd->flow->volume_density;
668                 tsmd->flow->surface_distance = smd->flow->surface_distance;
669                 tsmd->flow->particle_size = smd->flow->particle_size;
670                 tsmd->flow->subframes = smd->flow->subframes;
671
672                 tsmd->flow->texture_size = smd->flow->texture_size;
673                 tsmd->flow->texture_offset = smd->flow->texture_offset;
674                 BLI_strncpy(tsmd->flow->uvlayer_name, tsmd->flow->uvlayer_name, sizeof(tsmd->flow->uvlayer_name));
675                 tsmd->flow->vgroup_density = smd->flow->vgroup_density;
676
677                 tsmd->flow->type = smd->flow->type;
678                 tsmd->flow->source = smd->flow->source;
679                 tsmd->flow->texture_type = smd->flow->texture_type;
680                 tsmd->flow->flags = smd->flow->flags;
681         }
682         else if (tsmd->coll) {
683                 /* leave it as initialized, collision settings is mostly caches */
684         }
685 }
686
687 #ifdef WITH_SMOKE
688
689 // forward decleration
690 static void smoke_calc_transparency(SmokeDomainSettings *sds, Scene *scene);
691 static float calc_voxel_transp(float *result, float *input, int res[3], int *pixel, float *tRay, float correct);
692
693 static int get_lamp(Scene *scene, float *light)
694 {
695         Base *base_tmp = NULL;
696         int found_lamp = 0;
697
698         // try to find a lamp, preferably local
699         for (base_tmp = scene->base.first; base_tmp; base_tmp = base_tmp->next) {
700                 if (base_tmp->object->type == OB_LAMP) {
701                         Lamp *la = base_tmp->object->data;
702
703                         if (la->type == LA_LOCAL) {
704                                 copy_v3_v3(light, base_tmp->object->obmat[3]);
705                                 return 1;
706                         }
707                         else if (!found_lamp) {
708                                 copy_v3_v3(light, base_tmp->object->obmat[3]);
709                                 found_lamp = 1;
710                         }
711                 }
712         }
713
714         return found_lamp;
715 }
716
717 /**********************************************************
718  *      Obstacles
719  **********************************************************/
720
721 static void obstacles_from_derivedmesh(Object *coll_ob, SmokeDomainSettings *sds, SmokeCollSettings *scs, unsigned char *obstacle_map, float *velocityX, float *velocityY, float *velocityZ, float dt)
722 {
723         if (!scs->dm) return;
724         {
725                 DerivedMesh *dm = NULL;
726                 MVert *mvert = NULL;
727                 MFace *mface = NULL;
728                 BVHTreeFromMesh treeData = {NULL};
729                 int numverts, i, z;
730
731                 float surface_distance = 0.6;
732
733                 float *vert_vel = NULL;
734                 int has_velocity = 0;
735
736                 tstart();
737
738                 dm = CDDM_copy(scs->dm);
739                 CDDM_calc_normals(dm);
740                 mvert = dm->getVertArray(dm);
741                 mface = dm->getTessFaceArray(dm);
742                 numverts = dm->getNumVerts(dm);
743
744                 // DG TODO
745                 // if (scs->type > SM_COLL_STATIC)
746                 // if line above is used, the code is in trouble if the object moves but is declared as "does not move"
747
748                 {
749                         vert_vel = MEM_callocN(sizeof(float) * numverts * 3, "smoke_obs_velocity");
750
751                         if (scs->numverts != numverts || !scs->verts_old) {
752                                 if (scs->verts_old) MEM_freeN(scs->verts_old);
753
754                                 scs->verts_old = MEM_callocN(sizeof(float) * numverts * 3, "smoke_obs_verts_old");
755                                 scs->numverts = numverts;
756                         }
757                         else {
758                                 has_velocity = 1;
759                         }
760                 }
761
762                 /*      Transform collider vertices to
763                  *   domain grid space for fast lookups */
764                 for (i = 0; i < numverts; i++) {
765                         float n[3];
766                         float co[3];
767
768                         /* vert pos */
769                         mul_m4_v3(coll_ob->obmat, mvert[i].co);
770                         smoke_pos_to_cell(sds, mvert[i].co);
771
772                         /* vert normal */
773                         normal_short_to_float_v3(n, mvert[i].no);
774                         mul_mat3_m4_v3(coll_ob->obmat, n);
775                         mul_mat3_m4_v3(sds->imat, n);
776                         normalize_v3(n);
777                         normal_float_to_short_v3(mvert[i].no, n);
778
779                         /* vert velocity */
780                         VECADD(co, mvert[i].co, sds->shift);
781                         if (has_velocity)
782                         {
783                                 sub_v3_v3v3(&vert_vel[i * 3], co, &scs->verts_old[i * 3]);
784                                 mul_v3_fl(&vert_vel[i * 3], sds->dx / dt);
785                         }
786                         copy_v3_v3(&scs->verts_old[i * 3], co);
787                 }
788
789                 if (bvhtree_from_mesh_faces(&treeData, dm, 0.0f, 4, 6)) {
790 #pragma omp parallel for schedule(static)
791                         for (z = sds->res_min[2]; z < sds->res_max[2]; z++) {
792                                 int x, y;
793                                 for (x = sds->res_min[0]; x < sds->res_max[0]; x++)
794                                         for (y = sds->res_min[1]; y < sds->res_max[1]; y++) {
795                                                 int index = smoke_get_index(x - sds->res_min[0], sds->res[0], y - sds->res_min[1], sds->res[1], z - sds->res_min[2]);
796
797                                                 float ray_start[3] = {(float)x + 0.5f, (float)y + 0.5f, (float)z + 0.5f};
798                                                 BVHTreeNearest nearest = {0};
799                                                 nearest.index = -1;
800                                                 nearest.dist_sq = surface_distance * surface_distance; /* find_nearest uses squared distance */
801
802                                                 /* find the nearest point on the mesh */
803                                                 if (BLI_bvhtree_find_nearest(treeData.tree, ray_start, &nearest, treeData.nearest_callback, &treeData) != -1) {
804                                                         float weights[4];
805                                                         int v1, v2, v3, f_index = nearest.index;
806
807                                                         /* calculate barycentric weights for nearest point */
808                                                         v1 = mface[f_index].v1;
809                                                         v2 = (nearest.flags & BVH_ONQUAD) ? mface[f_index].v3 : mface[f_index].v2;
810                                                         v3 = (nearest.flags & BVH_ONQUAD) ? mface[f_index].v4 : mface[f_index].v3;
811                                                         interp_weights_face_v3(weights, mvert[v1].co, mvert[v2].co, mvert[v3].co, NULL, nearest.co);
812
813                                                         // DG TODO
814                                                         if (has_velocity)
815                                                         {
816                                                                 /* apply object velocity */
817                                                                 {
818                                                                         float hit_vel[3];
819                                                                         interp_v3_v3v3v3(hit_vel, &vert_vel[v1 * 3], &vert_vel[v2 * 3], &vert_vel[v3 * 3], weights);
820                                                                         velocityX[index] += hit_vel[0];
821                                                                         velocityY[index] += hit_vel[1];
822                                                                         velocityZ[index] += hit_vel[2];
823                                                                 }
824                                                         }
825
826                                                         /* tag obstacle cells */
827                                                         obstacle_map[index] = 1;
828
829                                                         if (has_velocity)
830                                                                 obstacle_map[index] |= 8;
831                                                 }
832                                         }
833                         }
834                 }
835                 /* free bvh tree */
836                 free_bvhtree_from_mesh(&treeData);
837                 dm->release(dm);
838
839                 if (vert_vel) MEM_freeN(vert_vel);
840         }
841 }
842
843 /* Animated obstacles: dx_step = ((x_new - x_old) / totalsteps) * substep */
844 static void update_obstacles(Scene *scene, Object *ob, SmokeDomainSettings *sds, float dt,
845                              int UNUSED(substep), int UNUSED(totalsteps))
846 {
847         Object **collobjs = NULL;
848         unsigned int numcollobj = 0;
849
850         unsigned int collIndex;
851         unsigned char *obstacles = smoke_get_obstacle(sds->fluid);
852         float *velx = NULL;
853         float *vely = NULL;
854         float *velz = NULL;
855         float *velxOrig = smoke_get_velocity_x(sds->fluid);
856         float *velyOrig = smoke_get_velocity_y(sds->fluid);
857         float *velzOrig = smoke_get_velocity_z(sds->fluid);
858         float *density = smoke_get_density(sds->fluid);
859         float *fuel = smoke_get_fuel(sds->fluid);
860         float *flame = smoke_get_flame(sds->fluid);
861         float *r = smoke_get_color_r(sds->fluid);
862         float *g = smoke_get_color_g(sds->fluid);
863         float *b = smoke_get_color_b(sds->fluid);
864         unsigned int z;
865
866         smoke_get_ob_velocity(sds->fluid, &velx, &vely, &velz);
867
868         // TODO: delete old obstacle flags
869         for (z = 0; z < sds->res[0] * sds->res[1] * sds->res[2]; z++)
870         {
871                 if (obstacles[z] & 8) // Do not delete static obstacles
872                 {
873                         obstacles[z] = 0;
874                 }
875
876                 velx[z] = 0;
877                 vely[z] = 0;
878                 velz[z] = 0;
879         }
880
881
882         collobjs = get_collisionobjects(scene, ob, sds->coll_group, &numcollobj, eModifierType_Smoke);
883
884         // update obstacle tags in cells
885         for (collIndex = 0; collIndex < numcollobj; collIndex++)
886         {
887                 Object *collob = collobjs[collIndex];
888                 SmokeModifierData *smd2 = (SmokeModifierData *)modifiers_findByType(collob, eModifierType_Smoke);
889
890                 // DG TODO: check if modifier is active?
891
892                 if ((smd2->type & MOD_SMOKE_TYPE_COLL) && smd2->coll)
893                 {
894                         SmokeCollSettings *scs = smd2->coll;
895                         obstacles_from_derivedmesh(collob, sds, scs, obstacles, velx, vely, velz, dt);
896                 }
897         }
898
899         if (collobjs)
900                 MEM_freeN(collobjs);
901
902         /* obstacle cells should not contain any velocity from the smoke simulation */
903         for (z = 0; z < sds->res[0] * sds->res[1] * sds->res[2]; z++)
904         {
905                 if (obstacles[z])
906                 {
907                         velxOrig[z] = 0;
908                         velyOrig[z] = 0;
909                         velzOrig[z] = 0;
910                         density[z] = 0;
911                         if (fuel) {
912                                 fuel[z] = 0;
913                                 flame[z] = 0;
914                         }
915                         if (r) {
916                                 r[z] = 0;
917                                 g[z] = 0;
918                                 b[z] = 0;
919                         }
920                 }
921         }
922 }
923
924
925 /**********************************************************
926  *      Object subframe update method from dynamicpaint.c
927  **********************************************************/
928
929 /* set "ignore cache" flag for all caches on this object */
930 static void object_cacheIgnoreClear(Object *ob, int state)
931 {
932         ListBase pidlist;
933         PTCacheID *pid;
934         BKE_ptcache_ids_from_object(&pidlist, ob, NULL, 0);
935
936         for (pid = pidlist.first; pid; pid = pid->next) {
937                 if (pid->cache) {
938                         if (state)
939                                 pid->cache->flag |= PTCACHE_IGNORE_CLEAR;
940                         else
941                                 pid->cache->flag &= ~PTCACHE_IGNORE_CLEAR;
942                 }
943         }
944
945         BLI_freelistN(&pidlist);
946 }
947
948 static int subframe_updateObject(Scene *scene, Object *ob, int update_mesh, int parent_recursion, float frame, bool for_render)
949 {
950         SmokeModifierData *smd = (SmokeModifierData *)modifiers_findByType(ob, eModifierType_Smoke);
951         bConstraint *con;
952
953         /* if other is dynamic paint canvas, don't update */
954         if (smd && (smd->type & MOD_SMOKE_TYPE_DOMAIN))
955                 return 1;
956
957         /* if object has parents, update them too */
958         if (parent_recursion) {
959                 int recursion = parent_recursion - 1;
960                 int is_domain = 0;
961                 if (ob->parent) is_domain += subframe_updateObject(scene, ob->parent, 0, recursion, frame, for_render);
962                 if (ob->track) is_domain += subframe_updateObject(scene, ob->track, 0, recursion, frame, for_render);
963
964                 /* skip subframe if object is parented
965                  *  to vertex of a dynamic paint canvas */
966                 if (is_domain && (ob->partype == PARVERT1 || ob->partype == PARVERT3))
967                         return 0;
968
969                 /* also update constraint targets */
970                 for (con = ob->constraints.first; con; con = con->next) {
971                         bConstraintTypeInfo *cti = BKE_constraint_get_typeinfo(con);
972                         ListBase targets = {NULL, NULL};
973
974                         if (cti && cti->get_constraint_targets) {
975                                 bConstraintTarget *ct;
976                                 cti->get_constraint_targets(con, &targets);
977                                 for (ct = targets.first; ct; ct = ct->next) {
978                                         if (ct->tar)
979                                                 subframe_updateObject(scene, ct->tar, 0, recursion, frame, for_render);
980                                 }
981                                 /* free temp targets */
982                                 if (cti->flush_constraint_targets)
983                                         cti->flush_constraint_targets(con, &targets, 0);
984                         }
985                 }
986         }
987
988         /* was originally OB_RECALC_ALL - TODO - which flags are really needed??? */
989         ob->recalc |= OB_RECALC_OB | OB_RECALC_DATA | OB_RECALC_TIME;
990         BKE_animsys_evaluate_animdata(scene, &ob->id, ob->adt, frame, ADT_RECALC_ANIM);
991         if (update_mesh) {
992                 /* ignore cache clear during subframe updates
993                  *  to not mess up cache validity */
994                 object_cacheIgnoreClear(ob, 1);
995                 BKE_object_handle_update(G.main->eval_ctx, scene, ob);
996                 object_cacheIgnoreClear(ob, 0);
997         }
998         else
999                 BKE_object_where_is_calc_time(scene, ob, frame);
1000
1001         /* for curve following objects, parented curve has to be updated too */
1002         if (ob->type == OB_CURVE) {
1003                 Curve *cu = ob->data;
1004                 BKE_animsys_evaluate_animdata(scene, &cu->id, cu->adt, frame, ADT_RECALC_ANIM);
1005         }
1006         /* and armatures... */
1007         if (ob->type == OB_ARMATURE) {
1008                 bArmature *arm = ob->data;
1009                 BKE_animsys_evaluate_animdata(scene, &arm->id, arm->adt, frame, ADT_RECALC_ANIM);
1010                 BKE_pose_where_is(scene, ob);
1011         }
1012
1013         return 0;
1014 }
1015
1016 /**********************************************************
1017  *      Flow emission code
1018  **********************************************************/
1019
1020 typedef struct EmissionMap {
1021         float *influence;
1022         float *influence_high;
1023         float *velocity;
1024         int min[3], max[3], res[3];
1025         int hmin[3], hmax[3], hres[3];
1026         int total_cells, valid;
1027 } EmissionMap;
1028
1029 static void em_boundInsert(EmissionMap *em, float point[3])
1030 {
1031         int i = 0;
1032         if (!em->valid) {
1033                 for (; i < 3; i++) {
1034                         em->min[i] = (int)floor(point[i]);
1035                         em->max[i] = (int)ceil(point[i]);
1036                 }
1037                 em->valid = 1;
1038         }
1039         else {
1040                 for (; i < 3; i++) {
1041                         if (point[i] < em->min[i]) em->min[i] = (int)floor(point[i]);
1042                         if (point[i] > em->max[i]) em->max[i] = (int)ceil(point[i]);
1043                 }
1044         }
1045 }
1046
1047 static void clampBoundsInDomain(SmokeDomainSettings *sds, int min[3], int max[3], float *min_vel, float *max_vel, int margin, float dt)
1048 {
1049         int i;
1050         for (i = 0; i < 3; i++) {
1051                 int adapt = (sds->flags & MOD_SMOKE_ADAPTIVE_DOMAIN) ? sds->adapt_res : 0;
1052                 /* add margin */
1053                 min[i] -= margin;
1054                 max[i] += margin;
1055
1056                 /* adapt to velocity */
1057                 if (min_vel && min_vel[i] < 0.0f) {
1058                         min[i] += (int)floor(min_vel[i] * dt);
1059                 }
1060                 if (max_vel && max_vel[i] > 0.0f) {
1061                         max[i] += (int)ceil(max_vel[i] * dt);
1062                 }
1063
1064                 /* clamp within domain max size */
1065                 CLAMP(min[i], -adapt, sds->base_res[i] + adapt);
1066                 CLAMP(max[i], -adapt, sds->base_res[i] + adapt);
1067         }
1068 }
1069
1070 static void em_allocateData(EmissionMap *em, int use_velocity, int hires_mul)
1071 {
1072         int i, res[3];
1073
1074         for (i = 0; i < 3; i++) {
1075                 res[i] = em->max[i] - em->min[i];
1076                 if (res[i] <= 0)
1077                         return;
1078         }
1079         em->total_cells = res[0] * res[1] * res[2];
1080         copy_v3_v3_int(em->res, res);
1081
1082
1083         em->influence = MEM_callocN(sizeof(float) * em->total_cells, "smoke_flow_influence");
1084         if (use_velocity)
1085                 em->velocity = MEM_callocN(sizeof(float) * em->total_cells * 3, "smoke_flow_velocity");
1086
1087         /* allocate high resolution map if required */
1088         if (hires_mul > 1) {
1089                 int total_cells_high = em->total_cells * (hires_mul * hires_mul * hires_mul);
1090
1091                 for (i = 0; i < 3; i++) {
1092                         em->hmin[i] = em->min[i] * hires_mul;
1093                         em->hmax[i] = em->max[i] * hires_mul;
1094                         em->hres[i] = em->res[i] * hires_mul;
1095                 }
1096
1097                 em->influence_high = MEM_callocN(sizeof(float) * total_cells_high, "smoke_flow_influence_high");
1098         }
1099         em->valid = 1;
1100 }
1101
1102 static void em_freeData(EmissionMap *em)
1103 {
1104         if (em->influence)
1105                 MEM_freeN(em->influence);
1106         if (em->influence_high)
1107                 MEM_freeN(em->influence_high);
1108         if (em->velocity)
1109                 MEM_freeN(em->velocity);
1110 }
1111
1112 static void em_combineMaps(EmissionMap *output, EmissionMap *em2, int hires_multiplier, int additive, float sample_size)
1113 {
1114         int i, x, y, z;
1115
1116         /* copyfill input 1 struct and clear output for new allocation */
1117         EmissionMap em1;
1118         memcpy(&em1, output, sizeof(EmissionMap));
1119         memset(output, 0, sizeof(EmissionMap));
1120
1121         for (i = 0; i < 3; i++) {
1122                 if (em1.valid) {
1123                         output->min[i] = MIN2(em1.min[i], em2->min[i]);
1124                         output->max[i] = MAX2(em1.max[i], em2->max[i]);
1125                 }
1126                 else {
1127                         output->min[i] = em2->min[i];
1128                         output->max[i] = em2->max[i];
1129                 }
1130         }
1131         /* allocate output map */
1132         em_allocateData(output, (em1.velocity || em2->velocity), hires_multiplier);
1133
1134         /* base resolution inputs */
1135         for (x = output->min[0]; x < output->max[0]; x++)
1136                 for (y = output->min[1]; y < output->max[1]; y++)
1137                         for (z = output->min[2]; z < output->max[2]; z++) {
1138                                 int index_out = smoke_get_index(x - output->min[0], output->res[0], y - output->min[1], output->res[1], z - output->min[2]);
1139
1140                                 /* initialize with first input if in range */
1141                                 if (x >= em1.min[0] && x < em1.max[0] &&
1142                                         y >= em1.min[1] && y < em1.max[1] &&
1143                                         z >= em1.min[2] && z < em1.max[2]) {
1144                                         int index_in = smoke_get_index(x - em1.min[0], em1.res[0], y - em1.min[1], em1.res[1], z - em1.min[2]);
1145
1146                                         /* values */
1147                                         output->influence[index_out] = em1.influence[index_in];
1148                                         if (output->velocity && em1.velocity) {
1149                                                 copy_v3_v3(&output->velocity[index_out * 3], &em1.velocity[index_in * 3]);
1150                                         }
1151                                 }
1152
1153                                 /* apply second input if in range */
1154                                 if (x >= em2->min[0] && x < em2->max[0] &&
1155                                         y >= em2->min[1] && y < em2->max[1] &&
1156                                         z >= em2->min[2] && z < em2->max[2]) {
1157                                         int index_in = smoke_get_index(x - em2->min[0], em2->res[0], y - em2->min[1], em2->res[1], z - em2->min[2]);
1158
1159                                         /* values */
1160                                         if (additive) {
1161                                                 output->influence[index_out] += em2->influence[index_in] * sample_size;
1162                                         }
1163                                         else {
1164                                                 output->influence[index_out] = MAX2(em2->influence[index_in], output->influence[index_out]);
1165                                         }
1166                                         if (output->velocity && em2->velocity) {
1167                                                 /* last sample replaces the velocity */
1168                                                 output->velocity[index_out * 3]         = ADD_IF_LOWER(output->velocity[index_out * 3], em2->velocity[index_in * 3]);
1169                                                 output->velocity[index_out * 3 + 1] = ADD_IF_LOWER(output->velocity[index_out * 3 + 1], em2->velocity[index_in * 3 + 1]);
1170                                                 output->velocity[index_out * 3 + 2] = ADD_IF_LOWER(output->velocity[index_out * 3 + 2], em2->velocity[index_in * 3 + 2]);
1171                                         }
1172                                 }
1173         } // low res loop
1174
1175
1176
1177         /* initialize high resolution input if available */
1178         if (output->influence_high) {
1179                 for (x = output->hmin[0]; x < output->hmax[0]; x++)
1180                         for (y = output->hmin[1]; y < output->hmax[1]; y++)
1181                                 for (z = output->hmin[2]; z < output->hmax[2]; z++) {
1182                                         int index_out = smoke_get_index(x - output->hmin[0], output->hres[0], y - output->hmin[1], output->hres[1], z - output->hmin[2]);
1183
1184                                         /* initialize with first input if in range */
1185                                         if (x >= em1.hmin[0] && x < em1.hmax[0] &&
1186                                                 y >= em1.hmin[1] && y < em1.hmax[1] &&
1187                                                 z >= em1.hmin[2] && z < em1.hmax[2]) {
1188                                                 int index_in = smoke_get_index(x - em1.hmin[0], em1.hres[0], y - em1.hmin[1], em1.hres[1], z - em1.hmin[2]);
1189                                                 /* values */
1190                                                 output->influence_high[index_out] = em1.influence_high[index_in];
1191                                         }
1192
1193                                         /* apply second input if in range */
1194                                         if (x >= em2->hmin[0] && x < em2->hmax[0] &&
1195                                                 y >= em2->hmin[1] && y < em2->hmax[1] &&
1196                                                 z >= em2->hmin[2] && z < em2->hmax[2]) {
1197                                                 int index_in = smoke_get_index(x - em2->hmin[0], em2->hres[0], y - em2->hmin[1], em2->hres[1], z - em2->hmin[2]);
1198
1199                                                 /* values */
1200                                                 if (additive) {
1201                                                         output->influence_high[index_out] += em2->influence_high[index_in] * sample_size;
1202                                                 }
1203                                                 else {
1204                                                         output->influence_high[index_out] = MAX2(em2->influence_high[index_in], output->influence_high[index_out]);
1205                                                 }
1206                                         }
1207                 } // high res loop
1208         }
1209
1210         /* free original data */
1211         em_freeData(&em1);
1212 }
1213
1214
1215 static void emit_from_particles(Object *flow_ob, SmokeDomainSettings *sds, SmokeFlowSettings *sfs, EmissionMap *em, Scene *scene, float dt)
1216 {
1217         if (sfs && sfs->psys && sfs->psys->part && ELEM(sfs->psys->part->type, PART_EMITTER, PART_FLUID)) // is particle system selected
1218         {
1219                 ParticleSimulationData sim;
1220                 ParticleSystem *psys = sfs->psys;
1221                 float *particle_pos;
1222                 float *particle_vel;
1223                 int totpart = psys->totpart, totchild;
1224                 int p = 0;
1225                 int valid_particles = 0;
1226                 int bounds_margin = 1;
1227
1228                 /* radius based flow */
1229                 float solid = sfs->particle_size * 0.5f;
1230                 float smooth = 0.5f; /* add 0.5 cells of linear falloff to reduce aliasing */
1231                 int hires_multiplier = 1;
1232                 int i, z;
1233                 KDTree *tree;
1234
1235                 sim.scene = scene;
1236                 sim.ob = flow_ob;
1237                 sim.psys = psys;
1238
1239                 /* initialize particle cache */
1240                 if (psys->part->type == PART_HAIR) {
1241                         // TODO: PART_HAIR not supported whatsoever
1242                         totchild = 0;
1243                 }
1244                 else {
1245                         totchild = psys->totchild * psys->part->disp / 100;
1246                 }
1247
1248                 particle_pos = MEM_callocN(sizeof(float) * (totpart + totchild) * 3, "smoke_flow_particles");
1249                 particle_vel = MEM_callocN(sizeof(float) * (totpart + totchild) * 3, "smoke_flow_particles");
1250
1251                 /* setup particle radius emission if enabled */
1252                 if (sfs->flags & MOD_SMOKE_FLOW_USE_PART_SIZE) {
1253                         tree = BLI_kdtree_new(psys->totpart + psys->totchild);
1254
1255                         /* check need for high resolution map */
1256                         if ((sds->flags & MOD_SMOKE_HIGHRES) && (sds->highres_sampling == SM_HRES_FULLSAMPLE)) {
1257                                 hires_multiplier = sds->amplify + 1;
1258                         }
1259
1260                         bounds_margin = (int)ceil(solid + smooth);
1261                 }
1262
1263                 /* calculate local position for each particle */
1264                 for (p = 0; p < totpart + totchild; p++)
1265                 {
1266                         ParticleKey state;
1267                         float *pos;
1268                         if (p < totpart) {
1269                                 if (psys->particles[p].flag & (PARS_NO_DISP | PARS_UNEXIST))
1270                                         continue;
1271                         }
1272                         else {
1273                                 /* handle child particle */
1274                                 ChildParticle *cpa = &psys->child[p - totpart];
1275                                 if (psys->particles[cpa->parent].flag & (PARS_NO_DISP | PARS_UNEXIST))
1276                                         continue;
1277                         }
1278
1279                         state.time = BKE_scene_frame_get(scene); /* use scene time */
1280                         if (psys_get_particle_state(&sim, p, &state, 0) == 0)
1281                                 continue;
1282
1283                         /* location */
1284                         pos = &particle_pos[valid_particles * 3];
1285                         copy_v3_v3(pos, state.co);
1286                         smoke_pos_to_cell(sds, pos);
1287
1288                         /* velocity */
1289                         copy_v3_v3(&particle_vel[valid_particles * 3], state.vel);
1290                         mul_mat3_m4_v3(sds->imat, &particle_vel[valid_particles * 3]);
1291
1292                         if (sfs->flags & MOD_SMOKE_FLOW_USE_PART_SIZE) {
1293                                 BLI_kdtree_insert(tree, valid_particles, pos, NULL);
1294                         }
1295
1296                         /* calculate emission map bounds */
1297                         em_boundInsert(em, pos);
1298                         valid_particles++;
1299                 }
1300
1301                 /* set emission map */
1302                 clampBoundsInDomain(sds, em->min, em->max, NULL, NULL, bounds_margin, dt);
1303                 em_allocateData(em, sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY, hires_multiplier);
1304
1305                 if (!(sfs->flags & MOD_SMOKE_FLOW_USE_PART_SIZE)) {
1306                         for (p = 0; p < valid_particles; p++)
1307                         {
1308                                 int cell[3];
1309                                 size_t i = 0;
1310                                 size_t index = 0;
1311                                 int badcell = 0;
1312
1313                                 /* 1. get corresponding cell */
1314                                 cell[0] = floor(particle_pos[p * 3]) - em->min[0];
1315                                 cell[1] = floor(particle_pos[p * 3 + 1]) - em->min[1];
1316                                 cell[2] = floor(particle_pos[p * 3 + 2]) - em->min[2];
1317                                 /* check if cell is valid (in the domain boundary) */
1318                                 for (i = 0; i < 3; i++) {
1319                                         if ((cell[i] > em->res[i] - 1) || (cell[i] < 0)) {
1320                                                 badcell = 1;
1321                                                 break;
1322                                         }
1323                                 }
1324                                 if (badcell)
1325                                         continue;
1326                                 /* get cell index */
1327                                 index = smoke_get_index(cell[0], em->res[0], cell[1], em->res[1], cell[2]);
1328                                 /* Add influence to emission map */
1329                                 em->influence[index] = 1.0f;
1330                                 /* Uses particle velocity as initial velocity for smoke */
1331                                 if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY && (psys->part->phystype != PART_PHYS_NO))
1332                                 {
1333                                         VECADDFAC(&em->velocity[index * 3], &em->velocity[index * 3], &particle_vel[p * 3], sfs->vel_multi);
1334                                 }
1335                         }   // particles loop
1336                 }
1337                 else if (valid_particles > 0) { // MOD_SMOKE_FLOW_USE_PART_SIZE
1338
1339                         int min[3], max[3], res[3];
1340                         float hr = 1.0f / ((float)hires_multiplier);
1341                         /* slightly adjust high res antialias smoothness based on number of divisions
1342                          * to allow smaller details but yet not differing too much from the low res size */
1343                         const float hr_smooth = smooth * powf(hr, 1.0f / 3.0f);
1344
1345                         /* setup loop bounds */
1346                         for (i = 0; i < 3; i++) {
1347                                 min[i] = em->min[i] * hires_multiplier;
1348                                 max[i] = em->max[i] * hires_multiplier;
1349                                 res[i] = em->res[i] * hires_multiplier;
1350                         }
1351
1352                         BLI_kdtree_balance(tree);
1353
1354                         /* begin thread safe malloc */
1355                         BLI_begin_threaded_malloc();
1356
1357 #pragma omp parallel for schedule(static)
1358                         for (z = min[2]; z < max[2]; z++) {
1359                                 int x, y;
1360                                 for (x = min[0]; x < max[0]; x++)
1361                                         for (y = min[1]; y < max[1]; y++) {
1362                                                 /* take low res samples where possible */
1363                                                 if (hires_multiplier <= 1 || !(x % hires_multiplier || y % hires_multiplier || z % hires_multiplier)) {
1364                                                         /* get low res space coordinates */
1365                                                         int lx = x / hires_multiplier;
1366                                                         int ly = y / hires_multiplier;
1367                                                         int lz = z / hires_multiplier;
1368
1369                                                         int index = smoke_get_index(lx - em->min[0], em->res[0], ly - em->min[1], em->res[1], lz - em->min[2]);
1370                                                         float ray_start[3] = {((float)lx) + 0.5f, ((float)ly) + 0.5f, ((float)lz) + 0.5f};
1371
1372                                                         /* find particle distance from the kdtree */
1373                                                         KDTreeNearest nearest;
1374                                                         float range = solid + smooth;
1375                                                         BLI_kdtree_find_nearest(tree, ray_start, NULL, &nearest);
1376
1377                                                         if (nearest.dist < range) {
1378                                                                 em->influence[index] = (nearest.dist < solid) ? 1.0f : (1.0f - (nearest.dist-solid) / smooth);
1379                                                                 /* Uses particle velocity as initial velocity for smoke */
1380                                                                 if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY && (psys->part->phystype != PART_PHYS_NO))
1381                                                                 {
1382                                                                         VECADDFAC(&em->velocity[index * 3], &em->velocity[index * 3], &particle_vel[nearest.index * 3], sfs->vel_multi);
1383                                                                 }
1384                                                         }
1385                                                 }
1386
1387                                                 /* take high res samples if required */
1388                                                 if (hires_multiplier > 1) {
1389                                                         /* get low res space coordinates */
1390                                                         float lx = ((float)x) * hr;
1391                                                         float ly = ((float)y) * hr;
1392                                                         float lz = ((float)z) * hr;
1393
1394                                                         int index = smoke_get_index(x - min[0], res[0], y - min[1], res[1], z - min[2]);
1395                                                         float ray_start[3] = {lx + 0.5f*hr, ly + 0.5f*hr, lz + 0.5f*hr};
1396
1397                                                         /* find particle distance from the kdtree */
1398                                                         KDTreeNearest nearest;
1399                                                         float range = solid + hr_smooth;
1400                                                         BLI_kdtree_find_nearest(tree, ray_start, NULL, &nearest);
1401
1402                                                         if (nearest.dist < range) {
1403                                                                 em->influence_high[index] = (nearest.dist < solid) ? 1.0f : (1.0f - (nearest.dist-solid) / smooth);
1404                                                         }
1405                                                 }
1406
1407                                         }
1408                         }
1409                         BLI_end_threaded_malloc();
1410                 }
1411
1412                 if (sfs->flags & MOD_SMOKE_FLOW_USE_PART_SIZE) {
1413                         BLI_kdtree_free(tree);
1414                 }
1415
1416                 /* free data */
1417                 if (particle_pos)
1418                         MEM_freeN(particle_pos);
1419                 if (particle_vel)
1420                         MEM_freeN(particle_vel);
1421         }
1422 }
1423
1424 static void sample_derivedmesh(SmokeFlowSettings *sfs, MVert *mvert, MTFace *tface, MFace *mface, float *influence_map, float *velocity_map, int index, int base_res[3], float flow_center[3], BVHTreeFromMesh *treeData, float ray_start[3],
1425                                                                 float *vert_vel, int has_velocity, int defgrp_index, MDeformVert *dvert, float x, float y, float z)
1426 {
1427         float ray_dir[3] = {1.0f, 0.0f, 0.0f};
1428         BVHTreeRayHit hit = {0};
1429         BVHTreeNearest nearest = {0};
1430
1431         float volume_factor = 0.0f;
1432         float sample_str = 0.0f;
1433
1434         hit.index = -1;
1435         hit.dist = 9999;
1436         nearest.index = -1;
1437         nearest.dist_sq = sfs->surface_distance * sfs->surface_distance; /* find_nearest uses squared distance */
1438
1439         /* Check volume collision */
1440         if (sfs->volume_density) {
1441                 if (BLI_bvhtree_ray_cast(treeData->tree, ray_start, ray_dir, 0.0f, &hit, treeData->raycast_callback, treeData) != -1) {
1442                         float dot = ray_dir[0] * hit.no[0] + ray_dir[1] * hit.no[1] + ray_dir[2] * hit.no[2];
1443                         /*  If ray and hit face normal are facing same direction
1444                          *      hit point is inside a closed mesh. */
1445                         if (dot >= 0) {
1446                                 /* Also cast a ray in opposite direction to make sure
1447                                  * point is at least surrounded by two faces */
1448                                 negate_v3(ray_dir);
1449                                 hit.index = -1;
1450                                 hit.dist = 9999;
1451
1452                                 BLI_bvhtree_ray_cast(treeData->tree, ray_start, ray_dir, 0.0f, &hit, treeData->raycast_callback, treeData);
1453                                 if (hit.index != -1) {
1454                                         volume_factor = sfs->volume_density;
1455                                 }
1456                         }
1457                 }
1458         }
1459
1460         /* find the nearest point on the mesh */
1461         if (BLI_bvhtree_find_nearest(treeData->tree, ray_start, &nearest, treeData->nearest_callback, treeData) != -1) {
1462                 float weights[4];
1463                 int v1, v2, v3, f_index = nearest.index;
1464                 float n1[3], n2[3], n3[3], hit_normal[3];
1465
1466                 /* emit from surface based on distance */
1467                 if (sfs->surface_distance) {
1468                         sample_str = sqrtf(nearest.dist_sq) / sfs->surface_distance;
1469                         CLAMP(sample_str, 0.0f, 1.0f);
1470                         sample_str = pow(1.0f - sample_str, 0.5f);
1471                 }
1472                 else
1473                         sample_str = 0.0f;
1474
1475                 /* calculate barycentric weights for nearest point */
1476                 v1 = mface[f_index].v1;
1477                 v2 = (nearest.flags & BVH_ONQUAD) ? mface[f_index].v3 : mface[f_index].v2;
1478                 v3 = (nearest.flags & BVH_ONQUAD) ? mface[f_index].v4 : mface[f_index].v3;
1479                 interp_weights_face_v3(weights, mvert[v1].co, mvert[v2].co, mvert[v3].co, NULL, nearest.co);
1480
1481                 if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY && velocity_map) {
1482                         /* apply normal directional velocity */
1483                         if (sfs->vel_normal) {
1484                                 /* interpolate vertex normal vectors to get nearest point normal */
1485                                 normal_short_to_float_v3(n1, mvert[v1].no);
1486                                 normal_short_to_float_v3(n2, mvert[v2].no);
1487                                 normal_short_to_float_v3(n3, mvert[v3].no);
1488                                 interp_v3_v3v3v3(hit_normal, n1, n2, n3, weights);
1489                                 normalize_v3(hit_normal);
1490                                 /* apply normal directional and random velocity
1491                                  * - TODO: random disabled for now since it doesnt really work well as pressure calc smoothens it out... */
1492                                 velocity_map[index * 3]   += hit_normal[0] * sfs->vel_normal * 0.25f;
1493                                 velocity_map[index * 3 + 1] += hit_normal[1] * sfs->vel_normal * 0.25f;
1494                                 velocity_map[index * 3 + 2] += hit_normal[2] * sfs->vel_normal * 0.25f;
1495                                 /* TODO: for fire emitted from mesh surface we can use
1496                                  *  Vf = Vs + (Ps/Pf - 1)*S to model gaseous expansion from solid to fuel */
1497                         }
1498                         /* apply object velocity */
1499                         if (has_velocity && sfs->vel_multi) {
1500                                 float hit_vel[3];
1501                                 interp_v3_v3v3v3(hit_vel, &vert_vel[v1 * 3], &vert_vel[v2 * 3], &vert_vel[v3 * 3], weights);
1502                                 velocity_map[index * 3]   += hit_vel[0] * sfs->vel_multi;
1503                                 velocity_map[index * 3 + 1] += hit_vel[1] * sfs->vel_multi;
1504                                 velocity_map[index * 3 + 2] += hit_vel[2] * sfs->vel_multi;
1505                         }
1506                 }
1507
1508                 /* apply vertex group influence if used */
1509                 if (defgrp_index != -1 && dvert) {
1510                         float weight_mask = defvert_find_weight(&dvert[v1], defgrp_index) * weights[0] +
1511                                             defvert_find_weight(&dvert[v2], defgrp_index) * weights[1] +
1512                                             defvert_find_weight(&dvert[v3], defgrp_index) * weights[2];
1513                         sample_str *= weight_mask;
1514                 }
1515
1516                 /* apply emission texture */
1517                 if ((sfs->flags & MOD_SMOKE_FLOW_TEXTUREEMIT) && sfs->noise_texture) {
1518                         float tex_co[3] = {0};
1519                         TexResult texres;
1520
1521                         if (sfs->texture_type == MOD_SMOKE_FLOW_TEXTURE_MAP_AUTO) {
1522                                 tex_co[0] = ((x - flow_center[0]) / base_res[0]) / sfs->texture_size;
1523                                 tex_co[1] = ((y - flow_center[1]) / base_res[1]) / sfs->texture_size;
1524                                 tex_co[2] = ((z - flow_center[2]) / base_res[2] - sfs->texture_offset) / sfs->texture_size;
1525                         }
1526                         else if (tface) {
1527                                 interp_v2_v2v2v2(tex_co, tface[f_index].uv[0], tface[f_index].uv[(nearest.flags & BVH_ONQUAD) ? 2 : 1],
1528                                                  tface[f_index].uv[(nearest.flags & BVH_ONQUAD) ? 3 : 2], weights);
1529                                 /* map between -1.0f and 1.0f */
1530                                 tex_co[0] = tex_co[0] * 2.0f - 1.0f;
1531                                 tex_co[1] = tex_co[1] * 2.0f - 1.0f;
1532                                 tex_co[2] = sfs->texture_offset;
1533                         }
1534                         texres.nor = NULL;
1535                         BKE_texture_get_value(NULL, sfs->noise_texture, tex_co, &texres, false);
1536                         sample_str *= texres.tin;
1537                 }
1538         }
1539
1540         /* multiply initial velocity by emitter influence */
1541         if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY && velocity_map) {
1542                 mul_v3_fl(&velocity_map[index * 3], sample_str);
1543         }
1544
1545         /* apply final influence based on volume factor */
1546         influence_map[index] = MAX2(volume_factor, sample_str);
1547 }
1548
1549 static void emit_from_derivedmesh(Object *flow_ob, SmokeDomainSettings *sds, SmokeFlowSettings *sfs, EmissionMap *em, float dt)
1550 {
1551         if (!sfs->dm) return;
1552         {
1553                 DerivedMesh *dm = sfs->dm;
1554                 int defgrp_index = sfs->vgroup_density - 1;
1555                 MDeformVert *dvert = NULL;
1556                 MVert *mvert = NULL;
1557                 MVert *mvert_orig = NULL;
1558                 MFace *mface = NULL;
1559                 MTFace *tface = NULL;
1560                 BVHTreeFromMesh treeData = {NULL};
1561                 int numOfVerts, i, z;
1562                 float flow_center[3] = {0};
1563
1564                 float *vert_vel = NULL;
1565                 int has_velocity = 0;
1566                 int min[3], max[3], res[3];
1567                 int hires_multiplier = 1;
1568
1569                 CDDM_calc_normals(dm);
1570                 mvert = dm->getVertArray(dm);
1571                 mvert_orig = dm->dupVertArray(dm);  /* copy original mvert and restore when done */
1572                 mface = dm->getTessFaceArray(dm);
1573                 numOfVerts = dm->getNumVerts(dm);
1574                 dvert = dm->getVertDataArray(dm, CD_MDEFORMVERT);
1575                 tface = CustomData_get_layer_named(&dm->faceData, CD_MTFACE, sfs->uvlayer_name);
1576
1577                 if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY) {
1578                         vert_vel = MEM_callocN(sizeof(float) * numOfVerts * 3, "smoke_flow_velocity");
1579
1580                         if (sfs->numverts != numOfVerts || !sfs->verts_old) {
1581                                 if (sfs->verts_old) MEM_freeN(sfs->verts_old);
1582                                 sfs->verts_old = MEM_callocN(sizeof(float) * numOfVerts * 3, "smoke_flow_verts_old");
1583                                 sfs->numverts = numOfVerts;
1584                         }
1585                         else {
1586                                 has_velocity = 1;
1587                         }
1588                 }
1589
1590                 /*      Transform dm vertices to
1591                  *   domain grid space for fast lookups */
1592                 for (i = 0; i < numOfVerts; i++) {
1593                         float n[3];
1594                         /* vert pos */
1595                         mul_m4_v3(flow_ob->obmat, mvert[i].co);
1596                         smoke_pos_to_cell(sds, mvert[i].co);
1597                         /* vert normal */
1598                         normal_short_to_float_v3(n, mvert[i].no);
1599                         mul_mat3_m4_v3(flow_ob->obmat, n);
1600                         mul_mat3_m4_v3(sds->imat, n);
1601                         normalize_v3(n);
1602                         normal_float_to_short_v3(mvert[i].no, n);
1603                         /* vert velocity */
1604                         if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY) {
1605                                 float co[3];
1606                                 VECADD(co, mvert[i].co, sds->shift);
1607                                 if (has_velocity) {
1608                                         sub_v3_v3v3(&vert_vel[i * 3], co, &sfs->verts_old[i * 3]);
1609                                         mul_v3_fl(&vert_vel[i * 3], sds->dx / dt);
1610                                 }
1611                                 copy_v3_v3(&sfs->verts_old[i * 3], co);
1612                         }
1613
1614                         /* calculate emission map bounds */
1615                         em_boundInsert(em, mvert[i].co);
1616                 }
1617                 mul_m4_v3(flow_ob->obmat, flow_center);
1618                 smoke_pos_to_cell(sds, flow_center);
1619
1620                 /* check need for high resolution map */
1621                 if ((sds->flags & MOD_SMOKE_HIGHRES) && (sds->highres_sampling == SM_HRES_FULLSAMPLE)) {
1622                         hires_multiplier = sds->amplify + 1;
1623                 }
1624
1625                 /* set emission map */
1626                 clampBoundsInDomain(sds, em->min, em->max, NULL, NULL, (int)ceil(sfs->surface_distance), dt);
1627                 em_allocateData(em, sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY, hires_multiplier);
1628
1629                 /* setup loop bounds */
1630                 for (i = 0; i < 3; i++) {
1631                         min[i] = em->min[i] * hires_multiplier;
1632                         max[i] = em->max[i] * hires_multiplier;
1633                         res[i] = em->res[i] * hires_multiplier;
1634                 }
1635
1636                 if (bvhtree_from_mesh_faces(&treeData, dm, 0.0f, 4, 6)) {
1637 #pragma omp parallel for schedule(static)
1638                         for (z = min[2]; z < max[2]; z++) {
1639                                 int x, y;
1640                                 for (x = min[0]; x < max[0]; x++)
1641                                         for (y = min[1]; y < max[1]; y++) {
1642                                                 /* take low res samples where possible */
1643                                                 if (hires_multiplier <= 1 || !(x % hires_multiplier || y % hires_multiplier || z % hires_multiplier)) {
1644                                                         /* get low res space coordinates */
1645                                                         int lx = x / hires_multiplier;
1646                                                         int ly = y / hires_multiplier;
1647                                                         int lz = z / hires_multiplier;
1648
1649                                                         int index = smoke_get_index(lx - em->min[0], em->res[0], ly - em->min[1], em->res[1], lz - em->min[2]);
1650                                                         float ray_start[3] = {((float)lx) + 0.5f, ((float)ly) + 0.5f, ((float)lz) + 0.5f};
1651
1652                                                         sample_derivedmesh(sfs, mvert, tface, mface, em->influence, em->velocity, index, sds->base_res, flow_center, &treeData, ray_start,
1653                                                                                                 vert_vel, has_velocity, defgrp_index, dvert, (float)lx, (float)ly, (float)lz);
1654                                                 }
1655
1656                                                 /* take high res samples if required */
1657                                                 if (hires_multiplier > 1) {
1658                                                         /* get low res space coordinates */
1659                                                         float hr = 1.0f / ((float)hires_multiplier);
1660                                                         float lx = ((float)x) * hr;
1661                                                         float ly = ((float)y) * hr;
1662                                                         float lz = ((float)z) * hr;
1663
1664                                                         int index = smoke_get_index(x - min[0], res[0], y - min[1], res[1], z - min[2]);
1665                                                         float ray_start[3] = {lx + 0.5f*hr, ly + 0.5f*hr, lz + 0.5f*hr};
1666
1667                                                         sample_derivedmesh(sfs, mvert, tface, mface, em->influence_high, NULL, index, sds->base_res, flow_center, &treeData, ray_start,
1668                                                                                                 vert_vel, has_velocity, defgrp_index, dvert, lx, ly, lz); /* x,y,z needs to be always lowres */
1669                                                 }
1670
1671                                         }
1672                         }
1673                 }
1674                 /* free bvh tree */
1675                 free_bvhtree_from_mesh(&treeData);
1676                 /* restore original mverts */
1677                 CustomData_set_layer(&dm->vertData, CD_MVERT, mvert_orig);
1678                 if (mvert)
1679                         MEM_freeN(mvert);
1680
1681                 if (vert_vel) MEM_freeN(vert_vel);
1682         }
1683 }
1684
1685 /**********************************************************
1686  *      Smoke step
1687  **********************************************************/
1688
1689 static void adjustDomainResolution(SmokeDomainSettings *sds, int new_shift[3], EmissionMap *emaps, unsigned int numflowobj, float dt)
1690 {
1691         int min[3] = {32767, 32767, 32767}, max[3] = {-32767, -32767, -32767}, res[3];
1692         int total_cells = 1, res_changed = 0, shift_changed = 0;
1693         float min_vel[3], max_vel[3];
1694         int x, y, z, i;
1695         float *density = smoke_get_density(sds->fluid);
1696         float *fuel = smoke_get_fuel(sds->fluid);
1697         float *bigdensity = smoke_turbulence_get_density(sds->wt);
1698         float *bigfuel = smoke_turbulence_get_fuel(sds->wt);
1699         float *vx = smoke_get_velocity_x(sds->fluid);
1700         float *vy = smoke_get_velocity_y(sds->fluid);
1701         float *vz = smoke_get_velocity_z(sds->fluid);
1702         int block_size = sds->amplify + 1;
1703         int wt_res[3];
1704
1705         if (sds->flags & MOD_SMOKE_HIGHRES && sds->wt) {
1706                 smoke_turbulence_get_res(sds->wt, wt_res);
1707         }
1708
1709         INIT_MINMAX(min_vel, max_vel);
1710
1711         /* Calculate bounds for current domain content */
1712         for (x = sds->res_min[0]; x <  sds->res_max[0]; x++)
1713                 for (y =  sds->res_min[1]; y <  sds->res_max[1]; y++)
1714                         for (z =  sds->res_min[2]; z <  sds->res_max[2]; z++)
1715                         {
1716                                 int xn = x - new_shift[0];
1717                                 int yn = y - new_shift[1];
1718                                 int zn = z - new_shift[2];
1719                                 int index;
1720                                 float max_den;
1721                                 
1722                                 /* skip if cell already belongs to new area */
1723                                 if (xn >= min[0] && xn <= max[0] && yn >= min[1] && yn <= max[1] && zn >= min[2] && zn <= max[2])
1724                                         continue;
1725
1726                                 index = smoke_get_index(x - sds->res_min[0], sds->res[0], y - sds->res_min[1], sds->res[1], z - sds->res_min[2]);
1727                                 max_den = (fuel) ? MAX2(density[index], fuel[index]) : density[index];
1728
1729                                 /* check high resolution bounds if max density isnt already high enough */
1730                                 if (max_den < sds->adapt_threshold && sds->flags & MOD_SMOKE_HIGHRES && sds->wt) {
1731                                         int i, j, k;
1732                                         /* high res grid index */
1733                                         int xx = (x - sds->res_min[0]) * block_size;
1734                                         int yy = (y - sds->res_min[1]) * block_size;
1735                                         int zz = (z - sds->res_min[2]) * block_size;
1736
1737                                         for (i = 0; i < block_size; i++)
1738                                                 for (j = 0; j < block_size; j++)
1739                                                         for (k = 0; k < block_size; k++)
1740                                                         {
1741                                                                 int big_index = smoke_get_index(xx + i, wt_res[0], yy + j, wt_res[1], zz + k);
1742                                                                 float den = (bigfuel) ? MAX2(bigdensity[big_index], bigfuel[big_index]) : bigdensity[big_index];
1743                                                                 if (den > max_den) {
1744                                                                         max_den = den;
1745                                                                 }
1746                                                         }
1747                                 }
1748
1749                                 /* content bounds (use shifted coordinates) */
1750                                 if (max_den >= sds->adapt_threshold) {
1751                                         if (min[0] > xn) min[0] = xn;
1752                                         if (min[1] > yn) min[1] = yn;
1753                                         if (min[2] > zn) min[2] = zn;
1754                                         if (max[0] < xn) max[0] = xn;
1755                                         if (max[1] < yn) max[1] = yn;
1756                                         if (max[2] < zn) max[2] = zn;
1757                                 }
1758
1759                                 /* velocity bounds */
1760                                 if (min_vel[0] > vx[index]) min_vel[0] = vx[index];
1761                                 if (min_vel[1] > vy[index]) min_vel[1] = vy[index];
1762                                 if (min_vel[2] > vz[index]) min_vel[2] = vz[index];
1763                                 if (max_vel[0] < vx[index]) max_vel[0] = vx[index];
1764                                 if (max_vel[1] < vy[index]) max_vel[1] = vy[index];
1765                                 if (max_vel[2] < vz[index]) max_vel[2] = vz[index];
1766                         }
1767
1768         /* also apply emission maps */
1769         for (i = 0; i < numflowobj; i++)
1770         {
1771                 EmissionMap *em = &emaps[i];
1772
1773                 for (x = em->min[0]; x < em->max[0]; x++)
1774                         for (y = em->min[1]; y < em->max[1]; y++)
1775                                 for (z = em->min[2]; z < em->max[2]; z++)
1776                                 {
1777                                         int index = smoke_get_index(x - em->min[0], em->res[0], y - em->min[1], em->res[1], z - em->min[2]);
1778                                         float max_den = em->influence[index];
1779
1780                                         /* density bounds */
1781                                         if (max_den >= sds->adapt_threshold) {
1782                                                 if (min[0] > x) min[0] = x;
1783                                                 if (min[1] > y) min[1] = y;
1784                                                 if (min[2] > z) min[2] = z;
1785                                                 if (max[0] < x) max[0] = x;
1786                                                 if (max[1] < y) max[1] = y;
1787                                                 if (max[2] < z) max[2] = z;
1788                                         }
1789                                 }
1790         }
1791
1792         /* calculate new bounds based on these values */
1793         mul_v3_fl(min_vel, 1.0f / sds->dx);
1794         mul_v3_fl(max_vel, 1.0f / sds->dx);
1795         clampBoundsInDomain(sds, min, max, min_vel, max_vel, sds->adapt_margin + 1, dt);
1796
1797         for (i = 0; i < 3; i++) {
1798                 /* calculate new resolution */
1799                 res[i] = max[i] - min[i];
1800                 total_cells *= res[i];
1801
1802                 if (new_shift[i])
1803                         shift_changed = 1;
1804
1805                 /* if no content set minimum dimensions */
1806                 if (res[i] <= 0) {
1807                         int j;
1808                         for (j = 0; j < 3; j++) {
1809                                 min[j] = 0;
1810                                 max[j] = 1;
1811                                 res[j] = 1;
1812                         }
1813                         res_changed = 1;
1814                         total_cells = 1;
1815                         break;
1816                 }
1817                 if (min[i] != sds->res_min[i] || max[i] != sds->res_max[i])
1818                         res_changed = 1;
1819         }
1820
1821         if (res_changed || shift_changed) {
1822                 struct FLUID_3D *fluid_old = sds->fluid;
1823                 struct WTURBULENCE *turb_old = sds->wt;
1824                 /* allocate new fluid data */
1825                 smoke_reallocate_fluid(sds, sds->dx, res, 0);
1826                 if (sds->flags & MOD_SMOKE_HIGHRES) {
1827                         smoke_reallocate_highres_fluid(sds, sds->dx, res, 0);
1828                 }
1829
1830                 /* copy values from old fluid to new */
1831                 if (sds->total_cells > 1 && total_cells > 1) {
1832                         /* low res smoke */
1833                         float *o_dens, *o_react, *o_flame, *o_fuel, *o_heat, *o_heatold, *o_vx, *o_vy, *o_vz, *o_r, *o_g, *o_b;
1834                         float *n_dens, *n_react, *n_flame, *n_fuel, *n_heat, *n_heatold, *n_vx, *n_vy, *n_vz, *n_r, *n_g, *n_b;
1835                         float dummy;
1836                         unsigned char *dummy_p;
1837                         /* high res smoke */
1838                         int wt_res_old[3];
1839                         float *o_wt_dens, *o_wt_react, *o_wt_flame, *o_wt_fuel, *o_wt_tcu, *o_wt_tcv, *o_wt_tcw, *o_wt_r, *o_wt_g, *o_wt_b;
1840                         float *n_wt_dens, *n_wt_react, *n_wt_flame, *n_wt_fuel, *n_wt_tcu, *n_wt_tcv, *n_wt_tcw, *n_wt_r, *n_wt_g, *n_wt_b;
1841
1842                         smoke_export(fluid_old, &dummy, &dummy, &o_dens, &o_react, &o_flame, &o_fuel, &o_heat, &o_heatold, &o_vx, &o_vy, &o_vz, &o_r, &o_g, &o_b, &dummy_p);
1843                         smoke_export(sds->fluid, &dummy, &dummy, &n_dens, &n_react, &n_flame, &n_fuel, &n_heat, &n_heatold, &n_vx, &n_vy, &n_vz, &n_r, &n_g, &n_b, &dummy_p);
1844
1845                         if (sds->flags & MOD_SMOKE_HIGHRES) {
1846                                 smoke_turbulence_export(turb_old, &o_wt_dens, &o_wt_react, &o_wt_flame, &o_wt_fuel, &o_wt_r, &o_wt_g, &o_wt_b, &o_wt_tcu, &o_wt_tcv, &o_wt_tcw);
1847                                 smoke_turbulence_get_res(turb_old, wt_res_old);
1848                                 smoke_turbulence_export(sds->wt, &n_wt_dens, &n_wt_react, &n_wt_flame, &n_wt_fuel, &n_wt_r, &n_wt_g, &n_wt_b, &n_wt_tcu, &n_wt_tcv, &n_wt_tcw);
1849                         }
1850
1851
1852                         for (x = sds->res_min[0]; x < sds->res_max[0]; x++)
1853                                 for (y = sds->res_min[1]; y < sds->res_max[1]; y++)
1854                                         for (z = sds->res_min[2]; z < sds->res_max[2]; z++)
1855                                         {
1856                                                 /* old grid index */
1857                                                 int xo = x - sds->res_min[0];
1858                                                 int yo = y - sds->res_min[1];
1859                                                 int zo = z - sds->res_min[2];
1860                                                 int index_old = smoke_get_index(xo, sds->res[0], yo, sds->res[1], zo);
1861                                                 /* new grid index */
1862                                                 int xn = x - min[0] - new_shift[0];
1863                                                 int yn = y - min[1] - new_shift[1];
1864                                                 int zn = z - min[2] - new_shift[2];
1865                                                 int index_new = smoke_get_index(xn, res[0], yn, res[1], zn);
1866
1867                                                 /* skip if outside new domain */
1868                                                 if (xn < 0 || xn >= res[0] ||
1869                                                     yn < 0 || yn >= res[1] ||
1870                                                     zn < 0 || zn >= res[2])
1871                                                         continue;
1872
1873                                                 /* copy data */
1874                                                 n_dens[index_new] = o_dens[index_old];
1875                                                 /* heat */
1876                                                 if (n_heat && o_heat) {
1877                                                         n_heat[index_new] = o_heat[index_old];
1878                                                         n_heatold[index_new] = o_heatold[index_old];
1879                                                 }
1880                                                 /* fuel */
1881                                                 if (n_fuel && o_fuel) {
1882                                                         n_flame[index_new] = o_flame[index_old];
1883                                                         n_fuel[index_new] = o_fuel[index_old];
1884                                                         n_react[index_new] = o_react[index_old];
1885                                                 }
1886                                                 /* color */
1887                                                 if (o_r && n_r) {
1888                                                         n_r[index_new] = o_r[index_old];
1889                                                         n_g[index_new] = o_g[index_old];
1890                                                         n_b[index_new] = o_b[index_old];
1891                                                 }
1892                                                 n_vx[index_new] = o_vx[index_old];
1893                                                 n_vy[index_new] = o_vy[index_old];
1894                                                 n_vz[index_new] = o_vz[index_old];
1895
1896                                                 if (sds->flags & MOD_SMOKE_HIGHRES && turb_old) {
1897                                                         int block_size = sds->amplify + 1;
1898                                                         int i, j, k;
1899                                                         /* old grid index */
1900                                                         int xx_o = xo * block_size;
1901                                                         int yy_o = yo * block_size;
1902                                                         int zz_o = zo * block_size;
1903                                                         /* new grid index */
1904                                                         int xx_n = xn * block_size;
1905                                                         int yy_n = yn * block_size;
1906                                                         int zz_n = zn * block_size;
1907
1908                                                         n_wt_tcu[index_new] = o_wt_tcu[index_old];
1909                                                         n_wt_tcv[index_new] = o_wt_tcv[index_old];
1910                                                         n_wt_tcw[index_new] = o_wt_tcw[index_old];
1911
1912                                                         for (i = 0; i < block_size; i++)
1913                                                                 for (j = 0; j < block_size; j++)
1914                                                                         for (k = 0; k < block_size; k++)
1915                                                                         {
1916                                                                                 int big_index_old = smoke_get_index(xx_o + i, wt_res_old[0], yy_o + j, wt_res_old[1], zz_o + k);
1917                                                                                 int big_index_new = smoke_get_index(xx_n + i, sds->res_wt[0], yy_n + j, sds->res_wt[1], zz_n + k);
1918                                                                                 /* copy data */
1919                                                                                 n_wt_dens[big_index_new] = o_wt_dens[big_index_old];
1920                                                                                 if (n_wt_flame && o_wt_flame) {
1921                                                                                         n_wt_flame[big_index_new] = o_wt_flame[big_index_old];
1922                                                                                         n_wt_fuel[big_index_new] = o_wt_fuel[big_index_old];
1923                                                                                         n_wt_react[big_index_new] = o_wt_react[big_index_old];
1924                                                                                 }
1925                                                                                 if (n_wt_r && o_wt_r) {
1926                                                                                         n_wt_r[big_index_new] = o_wt_r[big_index_old];
1927                                                                                         n_wt_g[big_index_new] = o_wt_g[big_index_old];
1928                                                                                         n_wt_b[big_index_new] = o_wt_b[big_index_old];
1929                                                                                 }
1930                                                                         }
1931                                                 }
1932                                         }
1933                 }
1934                 smoke_free(fluid_old);
1935                 if (turb_old)
1936                         smoke_turbulence_free(turb_old);
1937
1938                 /* set new domain dimensions */
1939                 VECCOPY(sds->res_min, min);
1940                 VECCOPY(sds->res_max, max);
1941                 VECCOPY(sds->res, res);
1942                 sds->total_cells = total_cells;
1943         }
1944 }
1945
1946 BLI_INLINE void apply_outflow_fields(int index, float *density, float *heat, float *fuel, float *react, float *color_r, float *color_g, float *color_b)
1947 {
1948         density[index] = 0.f;
1949         if (heat) {
1950                 heat[index] = 0.f;
1951         }
1952         if (fuel) {
1953                 fuel[index] = 0.f;
1954                 react[index] = 0.f;
1955         }
1956         if (color_r) {
1957                 color_r[index] = 0.f;
1958                 color_g[index] = 0.f;
1959                 color_b[index] = 0.f;
1960         }
1961 }
1962
1963 BLI_INLINE void apply_inflow_fields(SmokeFlowSettings *sfs, float emission_value, int index, float *density, float *heat, float *fuel, float *react, float *color_r, float *color_g, float *color_b)
1964 {
1965         int absolute_flow = (sfs->flags & MOD_SMOKE_FLOW_ABSOLUTE);
1966         float dens_old = density[index];
1967         // float fuel_old = (fuel) ? fuel[index] : 0.0f;  /* UNUSED */
1968         float dens_flow = (sfs->type == MOD_SMOKE_FLOW_TYPE_FIRE) ? 0.0f : emission_value * sfs->density;
1969         float fuel_flow = emission_value * sfs->fuel_amount;
1970         /* add heat */
1971         if (heat && emission_value > 0.0f) {
1972                 heat[index] = ADD_IF_LOWER(heat[index], sfs->temp);
1973         }
1974         /* absolute */
1975         if (absolute_flow) {
1976                 if (sfs->type != MOD_SMOKE_FLOW_TYPE_FIRE) {
1977                         if (dens_flow > density[index])
1978                                 density[index] = dens_flow;
1979                 }
1980                 if (sfs->type != MOD_SMOKE_FLOW_TYPE_SMOKE && fuel && fuel_flow) {
1981                         if (fuel_flow > fuel[index])
1982                                 fuel[index] = fuel_flow;
1983                 }
1984         }
1985         /* additive */
1986         else {
1987                 if (sfs->type != MOD_SMOKE_FLOW_TYPE_FIRE) {
1988                         density[index] += dens_flow;
1989                         CLAMP(density[index], 0.0f, 1.0f);
1990                 }
1991                 if (sfs->type != MOD_SMOKE_FLOW_TYPE_SMOKE && fuel && sfs->fuel_amount) {
1992                         fuel[index] += fuel_flow;
1993                         CLAMP(fuel[index], 0.0f, 10.0f);
1994                 }
1995         }
1996
1997         /* set color */
1998         if (color_r && dens_flow) {
1999                 float total_dens = density[index] / (dens_old + dens_flow);
2000                 color_r[index] = (color_r[index] + sfs->color[0] * dens_flow) * total_dens;
2001                 color_g[index] = (color_g[index] + sfs->color[1] * dens_flow) * total_dens;
2002                 color_b[index] = (color_b[index] + sfs->color[2] * dens_flow) * total_dens;
2003         }
2004
2005         /* set fire reaction coordinate */
2006         if (fuel && fuel[index] > FLT_EPSILON) {
2007                 /* instead of using 1.0 for all new fuel add slight falloff
2008                  * to reduce flow blockiness */
2009                 float value = 1.0f - powf(1.0f - emission_value, 2.0f);
2010
2011                 if (value > react[index]) {
2012                         float f = fuel_flow / fuel[index];
2013                         react[index] = value * f + (1.0f - f) * react[index];
2014                         CLAMP(react[index], 0.0f, value);
2015                 }
2016         }
2017 }
2018
2019 static void update_flowsfluids(Scene *scene, Object *ob, SmokeDomainSettings *sds, float dt, bool for_render)
2020 {
2021         Object **flowobjs = NULL;
2022         EmissionMap *emaps = NULL;
2023         unsigned int numflowobj = 0;
2024         unsigned int flowIndex;
2025         int new_shift[3] = {0};
2026         int active_fields = sds->active_fields;
2027
2028         /* calculate domain shift for current frame if using adaptive domain */
2029         if (sds->flags & MOD_SMOKE_ADAPTIVE_DOMAIN) {
2030                 int total_shift[3];
2031                 float frame_shift_f[3];
2032                 float ob_loc[3] = {0};
2033
2034                 mul_m4_v3(ob->obmat, ob_loc);
2035
2036                 VECSUB(frame_shift_f, ob_loc, sds->prev_loc);
2037                 copy_v3_v3(sds->prev_loc, ob_loc);
2038                 /* convert global space shift to local "cell" space */
2039                 mul_mat3_m4_v3(sds->imat, frame_shift_f);
2040                 frame_shift_f[0] = frame_shift_f[0] / sds->cell_size[0];
2041                 frame_shift_f[1] = frame_shift_f[1] / sds->cell_size[1];
2042                 frame_shift_f[2] = frame_shift_f[2] / sds->cell_size[2];
2043                 /* add to total shift */
2044                 VECADD(sds->shift_f, sds->shift_f, frame_shift_f);
2045                 /* convert to integer */
2046                 total_shift[0] = floor(sds->shift_f[0]);
2047                 total_shift[1] = floor(sds->shift_f[1]);
2048                 total_shift[2] = floor(sds->shift_f[2]);
2049                 VECSUB(new_shift, total_shift, sds->shift);
2050                 copy_v3_v3_int(sds->shift, total_shift);
2051
2052                 /* calculate new domain boundary points so that smoke doesnt slide on sub-cell movement */
2053                 sds->p0[0] = sds->dp0[0] - sds->cell_size[0] * (sds->shift_f[0] - total_shift[0] - 0.5f);
2054                 sds->p0[1] = sds->dp0[1] - sds->cell_size[1] * (sds->shift_f[1] - total_shift[1] - 0.5f);
2055                 sds->p0[2] = sds->dp0[2] - sds->cell_size[2] * (sds->shift_f[2] - total_shift[2] - 0.5f);
2056                 sds->p1[0] = sds->p0[0] + sds->cell_size[0] * sds->base_res[0];
2057                 sds->p1[1] = sds->p0[1] + sds->cell_size[1] * sds->base_res[1];
2058                 sds->p1[2] = sds->p0[2] + sds->cell_size[2] * sds->base_res[2];
2059         }
2060
2061         flowobjs = get_collisionobjects(scene, ob, sds->fluid_group, &numflowobj, eModifierType_Smoke);
2062
2063         /* init emission maps for each flow */
2064         emaps = MEM_callocN(sizeof(struct EmissionMap) * numflowobj, "smoke_flow_maps");
2065
2066         /* Prepare flow emission maps */
2067         for (flowIndex = 0; flowIndex < numflowobj; flowIndex++)
2068         {
2069                 Object *collob = flowobjs[flowIndex];
2070                 SmokeModifierData *smd2 = (SmokeModifierData *)modifiers_findByType(collob, eModifierType_Smoke);
2071
2072                 // check for initialized smoke object
2073                 if ((smd2->type & MOD_SMOKE_TYPE_FLOW) && smd2->flow)
2074                 {
2075                         // we got nice flow object
2076                         SmokeFlowSettings *sfs = smd2->flow;
2077                         int subframes = sfs->subframes;
2078                         EmissionMap *em = &emaps[flowIndex];
2079
2080                         /* just sample flow directly to emission map if no subframes */
2081                         if (!subframes) {
2082                                 if (sfs->source == MOD_SMOKE_FLOW_SOURCE_PARTICLES) {
2083                                         emit_from_particles(collob, sds, sfs, em, scene, dt);
2084                                 }
2085                                 else {
2086                                         emit_from_derivedmesh(collob, sds, sfs, em, dt);
2087                                 }
2088                         }
2089                         /* sample subframes */
2090                         else {
2091                                 int scene_frame = scene->r.cfra;
2092                                 // float scene_subframe = scene->r.subframe;  // UNUSED
2093                                 int subframe;
2094                                 for (subframe = 0; subframe <= subframes; subframe++) {
2095                                         EmissionMap em_temp = {NULL};
2096                                         float sample_size = 1.0f / (float)(subframes+1);
2097                                         float prev_frame_pos = sample_size * (float)(subframe+1);
2098                                         float sdt = dt * sample_size;
2099                                         int hires_multiplier = 1;
2100
2101                                         if ((sds->flags & MOD_SMOKE_HIGHRES) && (sds->highres_sampling == SM_HRES_FULLSAMPLE)) {
2102                                                 hires_multiplier = sds->amplify + 1;
2103                                         }
2104
2105                                         /* set scene frame to match previous frame + subframe
2106                                          * or use current frame for last sample */
2107                                         if (subframe < subframes) {
2108                                                 scene->r.cfra = scene_frame - 1;
2109                                                 scene->r.subframe = prev_frame_pos;
2110                                         }
2111                                         else {
2112                                                 scene->r.cfra = scene_frame;
2113                                                 scene->r.subframe = 0.0f;
2114                                         }
2115
2116                                         if (sfs->source == MOD_SMOKE_FLOW_SOURCE_PARTICLES) {
2117                                                 /* emit_from_particles() updates timestep internally */
2118                                                 emit_from_particles(collob, sds, sfs, &em_temp, scene, sdt);
2119                                         }
2120                                         else { /* MOD_SMOKE_FLOW_SOURCE_MESH */
2121                                                 /* update flow object frame */
2122                                                 subframe_updateObject(scene, collob, 1, 5, BKE_scene_frame_get(scene), for_render);
2123
2124                                                 /* apply flow */
2125                                                 emit_from_derivedmesh(collob, sds, sfs, &em_temp, sdt);
2126                                         }
2127
2128                                         /* combine emission maps */
2129                                         em_combineMaps(em, &em_temp, hires_multiplier, !(sfs->flags & MOD_SMOKE_FLOW_ABSOLUTE), sample_size);
2130                                         em_freeData(&em_temp);
2131                                 }
2132                         }
2133
2134                         /* update required data fields */
2135                         if (em->total_cells && sfs->type != MOD_SMOKE_FLOW_TYPE_OUTFLOW) {
2136                                 /* activate heat field if flow produces any heat */
2137                                 if (sfs->temp) {
2138                                         active_fields |= SM_ACTIVE_HEAT;
2139                                 }
2140                                 /* activate fuel field if flow adds any fuel */
2141                                 if (sfs->type != MOD_SMOKE_FLOW_TYPE_SMOKE && sfs->fuel_amount) {
2142                                         active_fields |= SM_ACTIVE_FIRE;
2143                                 }
2144                                 /* activate color field if flows add smoke with varying colors */
2145                                 if (sfs->type != MOD_SMOKE_FLOW_TYPE_FIRE && sfs->density) {
2146                                         if (!(active_fields & SM_ACTIVE_COLOR_SET)) {
2147                                                 copy_v3_v3(sds->active_color, sfs->color);
2148                                                 active_fields |= SM_ACTIVE_COLOR_SET;
2149                                         }
2150                                         else if (!equals_v3v3(sds->active_color, sfs->color)) {
2151                                                 active_fields |= SM_ACTIVE_COLORS;
2152                                         }
2153                                 }
2154                         }
2155                 }
2156         }
2157
2158         /* monitor active fields based on domain settings */
2159         /* if domain has fire, activate new fields if required */
2160         if (active_fields & SM_ACTIVE_FIRE) {
2161                 /* heat is always needed for fire */
2162                 active_fields |= SM_ACTIVE_HEAT;
2163                 /* also activate colors if domain smoke color differs from active color */
2164                 if (!(active_fields & SM_ACTIVE_COLOR_SET)) {
2165                         copy_v3_v3(sds->active_color, sds->flame_smoke_color);
2166                         active_fields |= SM_ACTIVE_COLOR_SET;
2167                 }
2168                 else if (!equals_v3v3(sds->active_color, sds->flame_smoke_color)) {
2169                         active_fields |= SM_ACTIVE_COLORS;
2170                 }
2171         }
2172
2173         /* Adjust domain size if needed */
2174         if (sds->flags & MOD_SMOKE_ADAPTIVE_DOMAIN) {
2175                 adjustDomainResolution(sds, new_shift, emaps, numflowobj, dt);
2176         }
2177
2178         /* Initialize new data fields if any */
2179         if (active_fields & SM_ACTIVE_HEAT) {
2180                 smoke_ensure_heat(sds->fluid);
2181         }
2182         if (active_fields & SM_ACTIVE_FIRE) {
2183                 smoke_ensure_fire(sds->fluid, sds->wt);
2184         }
2185         if (active_fields & SM_ACTIVE_COLORS) {
2186                 /* initialize all smoke with "active_color" */
2187                 smoke_ensure_colors(sds->fluid, sds->wt, sds->active_color[0], sds->active_color[1], sds->active_color[2]);
2188         }
2189         sds->active_fields = active_fields;
2190
2191         /* Apply emission data */
2192         if (sds->fluid) {
2193                 for (flowIndex = 0; flowIndex < numflowobj; flowIndex++)
2194                 {
2195                         Object *collob = flowobjs[flowIndex];
2196                         SmokeModifierData *smd2 = (SmokeModifierData *)modifiers_findByType(collob, eModifierType_Smoke);
2197
2198                         // check for initialized smoke object
2199                         if ((smd2->type & MOD_SMOKE_TYPE_FLOW) && smd2->flow)
2200                         {
2201                                 // we got nice flow object
2202                                 SmokeFlowSettings *sfs = smd2->flow;
2203                                 EmissionMap *em = &emaps[flowIndex];
2204
2205                                 float *density = smoke_get_density(sds->fluid);
2206                                 float *color_r = smoke_get_color_r(sds->fluid);
2207                                 float *color_g = smoke_get_color_g(sds->fluid);
2208                                 float *color_b = smoke_get_color_b(sds->fluid);
2209                                 float *fuel = smoke_get_fuel(sds->fluid);
2210                                 float *react = smoke_get_react(sds->fluid);
2211                                 float *bigdensity = smoke_turbulence_get_density(sds->wt);
2212                                 float *bigfuel = smoke_turbulence_get_fuel(sds->wt);
2213                                 float *bigreact = smoke_turbulence_get_react(sds->wt);
2214                                 float *bigcolor_r = smoke_turbulence_get_color_r(sds->wt);
2215                                 float *bigcolor_g = smoke_turbulence_get_color_g(sds->wt);
2216                                 float *bigcolor_b = smoke_turbulence_get_color_b(sds->wt);
2217                                 float *heat = smoke_get_heat(sds->fluid);
2218                                 float *velocity_x = smoke_get_velocity_x(sds->fluid);
2219                                 float *velocity_y = smoke_get_velocity_y(sds->fluid);
2220                                 float *velocity_z = smoke_get_velocity_z(sds->fluid);
2221                                 //unsigned char *obstacle = smoke_get_obstacle(sds->fluid);
2222                                 // DG TODO UNUSED unsigned char *obstacleAnim = smoke_get_obstacle_anim(sds->fluid);
2223                                 int bigres[3];
2224                                 float *velocity_map = em->velocity;
2225                                 float *emission_map = em->influence;
2226                                 float *emission_map_high = em->influence_high;
2227
2228                                 int ii, jj, kk, gx, gy, gz, ex, ey, ez, dx, dy, dz, block_size;
2229                                 size_t e_index, d_index, index_big;
2230
2231                                 // loop through every emission map cell
2232                                 for (gx = em->min[0]; gx < em->max[0]; gx++)
2233                                         for (gy = em->min[1]; gy < em->max[1]; gy++)
2234                                                 for (gz = em->min[2]; gz < em->max[2]; gz++)
2235                                                 {
2236                                                         /* get emission map index */
2237                                                         ex = gx - em->min[0];
2238                                                         ey = gy - em->min[1];
2239                                                         ez = gz - em->min[2];
2240                                                         e_index = smoke_get_index(ex, em->res[0], ey, em->res[1], ez);
2241
2242                                                         /* get domain index */
2243                                                         dx = gx - sds->res_min[0];
2244                                                         dy = gy - sds->res_min[1];
2245                                                         dz = gz - sds->res_min[2];
2246                                                         d_index = smoke_get_index(dx, sds->res[0], dy, sds->res[1], dz);
2247                                                         /* make sure emission cell is inside the new domain boundary */
2248                                                         if (dx < 0 || dy < 0 || dz < 0 || dx >= sds->res[0] || dy >= sds->res[1] || dz >= sds->res[2]) continue;
2249
2250                                                         if (sfs->type == MOD_SMOKE_FLOW_TYPE_OUTFLOW) { // outflow
2251                                                                 apply_outflow_fields(d_index, density, heat, fuel, react, color_r, color_g, color_b);
2252                                                         }
2253                                                         else { // inflow
2254                                                                 apply_inflow_fields(sfs, emission_map[e_index], d_index, density, heat, fuel, react, color_r, color_g, color_b);
2255
2256                                                                 /* initial velocity */
2257                                                                 if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY) {
2258                                                                         velocity_x[d_index] = ADD_IF_LOWER(velocity_x[d_index], velocity_map[e_index * 3]);
2259                                                                         velocity_y[d_index] = ADD_IF_LOWER(velocity_y[d_index], velocity_map[e_index * 3 + 1]);
2260                                                                         velocity_z[d_index] = ADD_IF_LOWER(velocity_z[d_index], velocity_map[e_index * 3 + 2]);
2261                                                                 }
2262                                                         }
2263
2264                                                         /* loop through high res blocks if high res enabled */
2265                                                         if (bigdensity) {
2266                                                                 // neighbor cell emission densities (for high resolution smoke smooth interpolation)
2267                                                                 float c000, c001, c010, c011,  c100, c101, c110, c111;
2268
2269                                                                 smoke_turbulence_get_res(sds->wt, bigres);
2270                                                                 block_size = sds->amplify + 1;  // high res block size
2271
2272                                                                 c000 = (ex > 0 && ey > 0 && ez > 0) ? emission_map[smoke_get_index(ex - 1, em->res[0], ey - 1, em->res[1], ez - 1)] : 0;
2273                                                                 c001 = (ex > 0 && ey > 0) ? emission_map[smoke_get_index(ex - 1, em->res[0], ey - 1, em->res[1], ez)] : 0;
2274                                                                 c010 = (ex > 0 && ez > 0) ? emission_map[smoke_get_index(ex - 1, em->res[0], ey, em->res[1], ez - 1)] : 0;
2275                                                                 c011 = (ex > 0) ? emission_map[smoke_get_index(ex - 1, em->res[0], ey, em->res[1], ez)] : 0;
2276
2277                                                                 c100 = (ey > 0 && ez > 0) ? emission_map[smoke_get_index(ex, em->res[0], ey - 1, em->res[1], ez - 1)] : 0;
2278                                                                 c101 = (ey > 0) ? emission_map[smoke_get_index(ex, em->res[0], ey - 1, em->res[1], ez)] : 0;
2279                                                                 c110 = (ez > 0) ? emission_map[smoke_get_index(ex, em->res[0], ey, em->res[1], ez - 1)] : 0;
2280                                                                 c111 = emission_map[smoke_get_index(ex, em->res[0], ey, em->res[1], ez)]; // this cell
2281
2282                                                                 for (ii = 0; ii < block_size; ii++)
2283                                                                         for (jj = 0; jj < block_size; jj++)
2284                                                                                 for (kk = 0; kk < block_size; kk++)
2285                                                                                 {
2286
2287                                                                                         float fx, fy, fz, interpolated_value;
2288                                                                                         int shift_x = 0, shift_y = 0, shift_z = 0;
2289
2290
2291                                                                                         /* Use full sample emission map if enabled and available */
2292                                                                                         if ((sds->highres_sampling == SM_HRES_FULLSAMPLE) && emission_map_high) {
2293                                                                                                 interpolated_value = emission_map_high[smoke_get_index(ex * block_size + ii, em->res[0] * block_size, ey * block_size + jj, em->res[1] * block_size, ez * block_size + kk)]; // this cell
2294                                                                                         }
2295                                                                                         else if (sds->highres_sampling == SM_HRES_NEAREST) {
2296                                                                                                 /* without interpolation use same low resolution
2297                                                                                                  * block value for all hi-res blocks */
2298                                                                                                 interpolated_value = c111;
2299                                                                                         }
2300                                                                                         /* Fall back to interpolated */
2301                                                                                         else
2302                                                                                         {
2303                                                                                                 /* get relative block position
2304                                                                                                  * for interpolation smoothing */
2305                                                                                                 fx = (float)ii / block_size + 0.5f / block_size;
2306                                                                                                 fy = (float)jj / block_size + 0.5f / block_size;
2307                                                                                                 fz = (float)kk / block_size + 0.5f / block_size;
2308
2309                                                                                                 /* calculate trilinear interpolation */
2310                                                                                                 interpolated_value = c000 * (1 - fx) * (1 - fy) * (1 - fz) +
2311                                                                                                                      c100 * fx * (1 - fy) * (1 - fz) +
2312                                                                                                                      c010 * (1 - fx) * fy * (1 - fz) +
2313                                                                                                                      c001 * (1 - fx) * (1 - fy) * fz +
2314                                                                                                                      c101 * fx * (1 - fy) * fz +
2315                                                                                                                      c011 * (1 - fx) * fy * fz +
2316                                                                                                                      c110 * fx * fy * (1 - fz) +
2317                                                                                                                      c111 * fx * fy * fz;
2318
2319
2320                                                                                                 /* add some contrast / sharpness
2321                                                                                                  * depending on hi-res block size */
2322                                                                                                 interpolated_value = (interpolated_value - 0.4f) * (block_size / 2) + 0.4f;
2323                                                                                                 CLAMP(interpolated_value, 0.0f, 1.0f);
2324
2325                                                                                                 /* shift smoke block index
2326                                                                                                  * (because pixel center is actually
2327                                                                                                  * in halfway of the low res block) */
2328                                                                                                 shift_x = (dx < 1) ? 0 : block_size / 2;
2329                                                                                                 shift_y = (dy < 1) ? 0 : block_size / 2;
2330                                                                                                 shift_z = (dz < 1) ? 0 : block_size / 2;
2331                                                                                         }
2332
2333                                                                                         /* get shifted index for current high resolution block */
2334                                                                                         index_big = smoke_get_index(block_size * dx + ii - shift_x, bigres[0], block_size * dy + jj - shift_y, bigres[1], block_size * dz + kk - shift_z);
2335
2336                                                                                         if (sfs->type == MOD_SMOKE_FLOW_TYPE_OUTFLOW) { // outflow
2337                                                                                                 if (interpolated_value) {
2338                                                                                                         apply_outflow_fields(index_big, bigdensity, NULL, bigfuel, bigreact, bigcolor_r, bigcolor_g, bigcolor_b);
2339                                                                                                 }
2340                                                                                         }
2341                                                                                         else { // inflow
2342                                                                                                 apply_inflow_fields(sfs, interpolated_value, index_big, bigdensity, NULL, bigfuel, bigreact, bigcolor_r, bigcolor_g, bigcolor_b);
2343                                                                                         }
2344                                                                                 } // hires loop
2345                                                         }  // bigdensity
2346                                                 } // low res loop
2347
2348                                 // free emission maps
2349                                 em_freeData(em);
2350
2351                         } // end emission
2352                 }
2353         }
2354
2355         if (flowobjs)
2356                 MEM_freeN(flowobjs);
2357         if (emaps)
2358                 MEM_freeN(emaps);
2359 }
2360
2361 static void update_effectors(Scene *scene, Object *ob, SmokeDomainSettings *sds, float UNUSED(dt))
2362 {
2363         ListBase *effectors;
2364         /* make sure smoke flow influence is 0.0f */
2365         sds->effector_weights->weight[PFIELD_SMOKEFLOW] = 0.0f;
2366         effectors = pdInitEffectors(scene, ob, NULL, sds->effector_weights, true);
2367
2368         if (effectors)
2369         {
2370                 float *density = smoke_get_density(sds->fluid);
2371                 float *fuel = smoke_get_fuel(sds->fluid);
2372                 float *force_x = smoke_get_force_x(sds->fluid);
2373                 float *force_y = smoke_get_force_y(sds->fluid);
2374                 float *force_z = smoke_get_force_z(sds->fluid);
2375                 float *velocity_x = smoke_get_velocity_x(sds->fluid);
2376                 float *velocity_y = smoke_get_velocity_y(sds->fluid);
2377                 float *velocity_z = smoke_get_velocity_z(sds->fluid);
2378                 unsigned char *obstacle = smoke_get_obstacle(sds->fluid);
2379                 int x;
2380
2381                 // precalculate wind forces
2382 #pragma omp parallel for schedule(static)
2383                 for (x = 0; x < sds->res[0]; x++)
2384                 {
2385                         int y, z;
2386                         for (y = 0; y < sds->res[1]; y++)
2387                                 for (z = 0; z < sds->res[2]; z++)
2388                                 {
2389                                         EffectedPoint epoint;
2390                                         float mag;
2391                                         float voxelCenter[3] = {0, 0, 0}, vel[3] = {0, 0, 0}, retvel[3] = {0, 0, 0};
2392                                         unsigned int index = smoke_get_index(x, sds->res[0], y, sds->res[1], z);
2393
2394                                         if (((fuel ? MAX2(density[index], fuel[index]) : density[index]) < FLT_EPSILON) || obstacle[index])
2395                                                 continue;
2396
2397                                         vel[0] = velocity_x[index];
2398                                         vel[1] = velocity_y[index];
2399                                         vel[2] = velocity_z[index];
2400
2401                                         /* convert vel to global space */
2402                                         mag = len_v3(vel);
2403                                         mul_mat3_m4_v3(sds->obmat, vel);
2404                                         normalize_v3(vel);
2405                                         mul_v3_fl(vel, mag);
2406
2407                                         voxelCenter[0] = sds->p0[0] + sds->cell_size[0] * ((float)(x + sds->res_min[0]) + 0.5f);
2408                                         voxelCenter[1] = sds->p0[1] + sds->cell_size[1] * ((float)(y + sds->res_min[1]) + 0.5f);
2409                                         voxelCenter[2] = sds->p0[2] + sds->cell_size[2] * ((float)(z + sds->res_min[2]) + 0.5f);
2410                                         mul_m4_v3(sds->obmat, voxelCenter);
2411
2412                                         pd_point_from_loc(scene, voxelCenter, vel, index, &epoint);
2413                                         pdDoEffectors(effectors, NULL, sds->effector_weights, &epoint, retvel, NULL);
2414
2415                                         /* convert retvel to local space */
2416                                         mag = len_v3(retvel);
2417                                         mul_mat3_m4_v3(sds->imat, retvel);
2418                                         normalize_v3(retvel);
2419                                         mul_v3_fl(retvel, mag);
2420
2421                                         // TODO dg - do in force!
2422                                         force_x[index] = min_ff(max_ff(-1.0f, retvel[0] * 0.2f), 1.0f);
2423                                         force_y[index] = min_ff(max_ff(-1.0f, retvel[1] * 0.2f), 1.0f);
2424                                         force_z[index] = min_ff(max_ff(-1.0f, retvel[2] * 0.2f), 1.0f);
2425                                 }
2426                 }
2427         }
2428
2429         pdEndEffectors(&effectors);
2430 }
2431
2432 static void step(Scene *scene, Object *ob, SmokeModifierData *smd, DerivedMesh *domain_dm, float fps, bool for_render)
2433 {
2434         SmokeDomainSettings *sds = smd->domain;
2435         /* stability values copied from wturbulence.cpp */
2436         const int maxSubSteps = 25;
2437         float maxVel;
2438         // maxVel should be 1.5 (1.5 cell max movement) * dx (cell size)
2439
2440         float dt;
2441         float maxVelMag = 0.0f;
2442         int totalSubsteps;
2443         int substep = 0;
2444         float dtSubdiv;
2445         float gravity[3] = {0.0f, 0.0f, -1.0f};
2446         float gravity_mag;
2447
2448 #if 0  /* UNUSED */
2449            /* get max velocity and lower the dt value if it is too high */
2450         size_t size = sds->res[0] * sds->res[1] * sds->res[2];
2451         float *velX = smoke_get_velocity_x(sds->fluid);
2452         float *velY = smoke_get_velocity_y(sds->fluid);
2453         float *velZ = smoke_get_velocity_z(sds->fluid);
2454         size_t i;
2455 #endif
2456
2457         /* update object state */
2458         invert_m4_m4(sds->imat, ob->obmat);
2459         copy_m4_m4(sds->obmat, ob->obmat);
2460         smoke_set_domain_from_derivedmesh(sds, ob, domain_dm, (sds->flags & MOD_SMOKE_ADAPTIVE_DOMAIN));
2461
2462         /* use global gravity if enabled */
2463         if (scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY) {
2464                 copy_v3_v3(gravity, scene->physics_settings.gravity);
2465                 /* map default value to 1.0 */
2466                 mul_v3_fl(gravity, 1.0f / 9.810f);
2467         }
2468         /* convert gravity to domain space */
2469         gravity_mag = len_v3(gravity);
2470         mul_mat3_m4_v3(sds->imat, gravity);
2471         normalize_v3(gravity);
2472         mul_v3_fl(gravity, gravity_mag);
2473
2474         /* adapt timestep for different framerates, dt = 0.1 is at 25fps */
2475         dt = DT_DEFAULT * (25.0f / fps);
2476         // maximum timestep/"CFL" constraint: dt < 5.0 *dx / maxVel
2477         maxVel = (sds->dx * 5.0f);
2478
2479 #if 0
2480         for (i = 0; i < size; i++) {
2481                 float vtemp = (velX[i] * velX[i] + velY[i] * velY[i] + velZ[i] * velZ[i]);
2482                 if (vtemp > maxVelMag)
2483                         maxVelMag = vtemp;
2484         }
2485 #endif
2486
2487         maxVelMag = sqrtf(maxVelMag) * dt * sds->time_scale;
2488         totalSubsteps = (int)((maxVelMag / maxVel) + 1.0f); /* always round up */
2489         totalSubsteps = (totalSubsteps < 1) ? 1 : totalSubsteps;
2490         totalSubsteps = (totalSubsteps > maxSubSteps) ? maxSubSteps : totalSubsteps;
2491
2492         /* Disable substeps for now, since it results in numerical instability */
2493         totalSubsteps = 1.0f;
2494
2495         dtSubdiv = (float)dt / (float)totalSubsteps;
2496
2497         // printf("totalSubsteps: %d, maxVelMag: %f, dt: %f\n", totalSubsteps, maxVelMag, dt);
2498
2499         for (substep = 0; substep < totalSubsteps; substep++)
2500         {
2501                 // calc animated obstacle velocities
2502                 update_flowsfluids(scene, ob, sds, dtSubdiv, for_render);
2503                 update_obstacles(scene, ob, sds, dtSubdiv, substep, totalSubsteps);
2504
2505                 if (sds->total_cells > 1) {
2506                         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
2507                         smoke_step(sds->fluid, gravity, dtSubdiv);
2508                 }
2509         }
2510 }
2511
2512 static DerivedMesh *createDomainGeometry(SmokeDomainSettings *sds, Object *ob)
2513 {
2514         DerivedMesh *result;
2515         MVert *mverts;
2516         MPoly *mpolys;
2517         MLoop *mloops;
2518         float min[3];
2519         float max[3];
2520         float *co;
2521         MPoly *mp;
2522         MLoop *ml;
2523
2524         int num_verts = 8;
2525         int num_faces = 6;
2526         int i;
2527         float ob_loc[3] = {0};
2528         float ob_cache_loc[3] = {0};
2529
2530         /* dont generate any mesh if there isnt any content */
2531         if (sds->total_cells <= 1) {
2532                 num_verts = 0;
2533                 num_faces = 0;
2534         }
2535
2536         result = CDDM_new(num_verts, 0, 0, num_faces * 4, num_faces);
2537         mverts = CDDM_get_verts(result);
2538         mpolys = CDDM_get_polys(result);
2539         mloops = CDDM_get_loops(result);
2540
2541
2542         if (num_verts) {
2543                 /* volume bounds */
2544                 VECMADD(min, sds->p0, sds->cell_size, sds->res_min);
2545                 VECMADD(max, sds->p0, sds->cell_size, sds->res_max);
2546
2547                 /* set vertices */
2548                 /* top slab */
2549                 co = mverts[0].co; co[0] = min[0]; co[1] = min[1]; co[2] = max[2];
2550                 co = mverts[1].co; co[0] = max[0]; co[1] = min[1]; co[2] = max[2];
2551                 co = mverts[2].co; co[0] = max[0]; co[1] = max[1]; co[2] = max[2];
2552                 co = mverts[3].co; co[0] = min[0]; co[1] = max[1]; co[2] = max[2];
2553                 /* bottom slab */
2554                 co = mverts[4].co; co[0] = min[0]; co[1] = min[1]; co[2] = min[2];
2555                 co = mverts[5].co; co[0] = max[0]; co[1] = min[1]; co[2] = min[2];
2556                 co = mverts[6].co; co[0] = max[0]; co[1] = max[1]; co[2] = min[2];
2557                 co = mverts[7].co; co[0] = min[0]; co[1] = max[1]; co[2] = min[2];
2558
2559                 /* create faces */
2560                 /* top */
2561                 mp = &mpolys[0]; ml = &mloops[0 * 4]; mp->loopstart = 0 * 4; mp->totloop = 4;
2562                 ml[0].v = 0; ml[1].v = 1; ml[2].v = 2; ml[3].v = 3;
2563                 /* right */
2564                 mp = &mpolys[1]; ml = &mloops[1 * 4]; mp->loopstart = 1 * 4; mp->totloop = 4;
2565                 ml[0].v = 2; ml[1].v = 1; ml[2].v = 5; ml[3].v = 6;
2566                 /* bottom */
2567                 mp = &mpolys[2]; ml = &mloops[2 * 4]; mp->loopstart = 2 * 4; mp->totloop = 4;
2568                 ml[0].v = 7; ml[1].v = 6; ml[2].v = 5; ml[3].v = 4;
2569                 /* left */
2570                 mp = &mpolys[3]; ml = &mloops[3 * 4]; mp->loopstart = 3 * 4; mp->totloop = 4;
2571                 ml[0].v = 0; ml[1].v = 3; ml[2].v = 7; ml[3].v = 4;
2572                 /* front */
2573                 mp = &mpolys[4]; ml = &mloops[4 * 4]; mp->loopstart = 4 * 4; mp->totloop = 4;
2574                 ml[0].v = 3; ml[1].v = 2; ml[2].v = 6; ml[3].v = 7;
2575                 /* back */
2576                 mp = &mpolys[5]; ml = &mloops[5 * 4]; mp->loopstart = 5 * 4; mp->totloop = 4;
2577                 ml[0].v = 1; ml[1].v = 0; ml[2].v = 4; ml[3].v = 5;
2578
2579                 /* calculate required shift to match domain's global position
2580                  *  it was originally simulated at (if object moves without smoke step) */
2581                 invert_m4_m4(ob->imat, ob->obmat);
2582                 mul_m4_v3(ob->obmat, ob_loc);
2583                 mul_m4_v3(sds->obmat, ob_cache_loc);
2584                 VECSUB(sds->obj_shift_f, ob_cache_loc, ob_loc);
2585                 /* convert shift to local space and apply to vertices */
2586                 mul_mat3_m4_v3(ob->imat, sds->obj_shift_f);
2587                 /* apply */
2588                 for (i = 0; i < num_verts; i++) {
2589                         add_v3_v3(mverts[i].co, sds->obj_shift_f);
2590                 }
2591         }
2592
2593
2594         CDDM_calc_edges(result);
2595         result->dirty |= DM_DIRTY_NORMALS;
2596         return result;
2597 }
2598
2599 static void smokeModifier_process(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedMesh *dm, bool for_render)
2600 {
2601         if ((smd->type & MOD_SMOKE_TYPE_FLOW))
2602         {
2603                 if (scene->r.cfra >= smd->time)
2604                         smokeModifier_init(smd, ob, scene, dm);
2605
2606                 if (smd->flow->dm) smd->flow->dm->release(smd->flow->dm);
2607                 smd->flow->dm = CDDM_copy(dm);
2608                 DM_ensure_tessface(smd->flow->dm);
2609
2610                 if (scene->r.cfra > smd->time)
2611                 {
2612                         smd->time = scene->r.cfra;
2613                 }
2614                 else if (scene->r.cfra < smd->time)
2615                 {
2616                         smd->time = scene->r.cfra;
2617                         smokeModifier_reset_ex(smd, false);
2618                 }
2619         }
2620         else if (smd->type & MOD_SMOKE_TYPE_COLL)
2621         {
2622                 if (scene->r.cfra >= smd->time)
2623                         smokeModifier_init(smd, ob, scene, dm);
2624
2625                 if (smd->coll)
2626                 {
2627                         if (smd->coll->dm)
2628                                 smd->coll->dm->release(smd->coll->dm);
2629
2630                         smd->coll->dm = CDDM_copy(dm);
2631                         DM_ensure_tessface(smd->coll->dm);
2632                 }
2633
2634                 smd->time = scene->r.cfra;
2635                 if (scene->r.cfra < smd->time)
2636                 {
2637                         smokeModifier_reset_ex(smd, false);
2638                 }
2639         }
2640         else if (smd->type & MOD_SMOKE_TYPE_DOMAIN)
2641         {
2642                 SmokeDomainSettings *sds = smd->domain;
2643                 PointCache *cache = NULL;
2644                 PTCacheID pid;
2645                 int startframe, endframe, framenr;
2646                 float timescale;
2647
2648                 framenr = scene->r.cfra;
2649
2650                 //printf("time: %d\n", scene->r.cfra);
2651
2652                 cache = sds->point_cache[0];
2653                 BKE_ptcache_id_from_smoke(&pid, ob, smd);
2654                 BKE_ptcache_id_time(&pid, scene, framenr, &startframe, &endframe, &timescale);
2655
2656                 if (!smd->domain->fluid || framenr == startframe)
2657                 {
2658                         BKE_ptcache_id_reset(scene, &pid, PTCACHE_RESET_OUTDATED);
2659                         smokeModifier_reset_ex(smd, false);
2660                         BKE_ptcache_validate(cache, framenr);
2661                         cache->flag &= ~PTCACHE_REDO_NEEDED;
2662                 }
2663
2664                 if (!smd->domain->fluid && (framenr != startframe) && (smd->domain->flags & MOD_SMOKE_FILE_LOAD) == 0 && (cache->flag & PTCACHE_BAKED) == 0)
2665                         return;
2666
2667                 smd->domain->flags &= ~MOD_SMOKE_FILE_LOAD;
2668                 CLAMP(framenr, startframe, endframe);
2669
2670                 /* If already viewing a pre/after frame, no need to reload */
2671                 if ((smd->time == framenr) && (framenr != scene->r.cfra))
2672                         return;
2673
2674                 if (smokeModifier_init(smd, ob, scene, dm) == 0)
2675                 {
2676                         printf("bad smokeModifier_init\n");
2677                         return;
2678                 }
2679
2680                 /* try to read from cache */
2681                 if (BKE_ptcache_read(&pid, (float)framenr) == PTCACHE_READ_EXACT) {
2682                         BKE_ptcache_validate(cache, framenr);
2683                         smd->time = framenr;
2684                         return;
2685                 }
2686
2687                 /* only calculate something when we advanced a single frame */
2688                 if (framenr != (int)smd->time + 1)
2689                         return;
2690
2691                 /* don't simulate if viewing start frame, but scene frame is not real start frame */
2692                 if (framenr != scene->r.cfra)
2693                         return;
2694
2695                 tstart();
2696
2697                 /* if on second frame, write cache for first frame */
2698                 if ((int)smd->time == startframe && (cache->flag & PTCACHE_OUTDATED || cache->last_exact == 0)) {
2699                         BKE_ptcache_write(&pid, startframe);
2700                 }
2701
2702                 // set new time
2703                 smd->time = scene->r.cfra;
2704
2705                 /* do simulation */
2706
2707                 // simulate the actual smoke (c++ code in intern/smoke)
2708                 // DG: interesting commenting this line + deactivating loading of noise files
2709                 if (framenr != startframe)
2710                 {
2711                         if (sds->flags & MOD_SMOKE_DISSOLVE) {
2712                                 /* low res dissolve */
2713                                 smoke_dissolve(sds->fluid, sds->diss_speed, sds->flags & MOD_SMOKE_DISSOLVE_LOG);
2714                                 /* high res dissolve */
2715                                 if (sds->wt) {
2716                                         smoke_dissolve_wavelet(sds->wt, sds->diss_speed, sds->flags & MOD_SMOKE_DISSOLVE_LOG);
2717                                 }
2718
2719                         }
2720
2721                         step(scene, ob, smd, dm, scene->r.frs_sec / scene->r.frs_sec_base, for_render);
2722                 }
2723
2724                 // create shadows before writing cache so they get stored
2725                 smoke_calc_transparency(sds, scene);
2726
2727                 if (sds->wt)
2728                 {
2729                         smoke_turbulence_step(sds->wt, sds->fluid);