Cleanup: style, use braces for blenkernel
[blender.git] / source / blender / blenkernel / intern / smoke.c
1 /*
2  * This program is free software; you can redistribute it and/or
3  * modify it under the terms of the GNU General Public License
4  * as published by the Free Software Foundation; either version 2
5  * of the License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program; if not, write to the Free Software Foundation,
14  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
15  *
16  * The Original Code is Copyright (C) Blender Foundation.
17  * All rights reserved.
18  */
19
20 /** \file
21  * \ingroup bke
22  */
23
24 /* Part of the code copied from elbeem fluid library, copyright by Nils Thuerey */
25
26 #include "MEM_guardedalloc.h"
27
28 #include <float.h>
29 #include <math.h>
30 #include <stdio.h>
31 #include <string.h> /* memset */
32
33 #include "BLI_blenlib.h"
34 #include "BLI_math.h"
35 #include "BLI_kdopbvh.h"
36 #include "BLI_threads.h"
37 #include "BLI_utildefines.h"
38
39 #include "DNA_anim_types.h"
40 #include "DNA_armature_types.h"
41 #include "DNA_constraint_types.h"
42 #include "DNA_customdata_types.h"
43 #include "DNA_light_types.h"
44 #include "DNA_mesh_types.h"
45 #include "DNA_meshdata_types.h"
46 #include "DNA_modifier_types.h"
47 #include "DNA_object_types.h"
48 #include "DNA_particle_types.h"
49 #include "DNA_scene_types.h"
50 #include "DNA_smoke_types.h"
51
52 #include "BKE_appdir.h"
53 #include "BKE_animsys.h"
54 #include "BKE_armature.h"
55 #include "BKE_bvhutils.h"
56 #include "BKE_collision.h"
57 #include "BKE_colortools.h"
58 #include "BKE_constraint.h"
59 #include "BKE_customdata.h"
60 #include "BKE_deform.h"
61 #include "BKE_effect.h"
62 #include "BKE_library.h"
63 #include "BKE_mesh.h"
64 #include "BKE_mesh_runtime.h"
65 #include "BKE_modifier.h"
66 #include "BKE_object.h"
67 #include "BKE_particle.h"
68 #include "BKE_pointcache.h"
69 #include "BKE_scene.h"
70 #include "BKE_smoke.h"
71 #include "BKE_texture.h"
72
73 #include "DEG_depsgraph.h"
74 #include "DEG_depsgraph_query.h"
75
76 #include "RE_shader_ext.h"
77
78 #include "GPU_glew.h"
79
80 /* UNUSED so far, may be enabled later */
81 /* #define USE_SMOKE_COLLISION_DM */
82
83 //#define DEBUG_TIME
84
85 #include "smoke_API.h"
86
87 #ifdef DEBUG_TIME
88 #  include "PIL_time.h"
89 #endif
90
91 #ifdef WITH_SMOKE
92 #  include "BLI_task.h"
93 #  include "BLI_kdtree.h"
94 #  include "BLI_voxel.h"
95
96 static ThreadMutex object_update_lock = BLI_MUTEX_INITIALIZER;
97
98 struct Mesh;
99 struct Object;
100 struct Scene;
101 struct SmokeModifierData;
102
103 // timestep default value for nice appearance 0.1f
104 #  define DT_DEFAULT 0.1f
105
106 #  define ADD_IF_LOWER_POS(a, b) (min_ff((a) + (b), max_ff((a), (b))))
107 #  define ADD_IF_LOWER_NEG(a, b) (max_ff((a) + (b), min_ff((a), (b))))
108 #  define ADD_IF_LOWER(a, b) (((b) > 0) ? ADD_IF_LOWER_POS((a), (b)) : ADD_IF_LOWER_NEG((a), (b)))
109
110 #else /* WITH_SMOKE */
111
112 /* Stubs to use when smoke is disabled */
113 struct WTURBULENCE *smoke_turbulence_init(int *UNUSED(res),
114                                           int UNUSED(amplify),
115                                           int UNUSED(noisetype),
116                                           const char *UNUSED(noisefile_path),
117                                           int UNUSED(use_fire),
118                                           int UNUSED(use_colors))
119 {
120   return NULL;
121 }
122 //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; }
123 void smoke_free(struct FLUID_3D *UNUSED(fluid))
124 {
125 }
126 float *smoke_get_density(struct FLUID_3D *UNUSED(fluid))
127 {
128   return NULL;
129 }
130 void smoke_turbulence_free(struct WTURBULENCE *UNUSED(wt))
131 {
132 }
133 void smoke_initWaveletBlenderRNA(struct WTURBULENCE *UNUSED(wt), float *UNUSED(strength))
134 {
135 }
136 void smoke_initBlenderRNA(struct FLUID_3D *UNUSED(fluid),
137                           float *UNUSED(alpha),
138                           float *UNUSED(beta),
139                           float *UNUSED(dt_factor),
140                           float *UNUSED(vorticity),
141                           int *UNUSED(border_colli),
142                           float *UNUSED(burning_rate),
143                           float *UNUSED(flame_smoke),
144                           float *UNUSED(flame_smoke_color),
145                           float *UNUSED(flame_vorticity),
146                           float *UNUSED(flame_ignition_temp),
147                           float *UNUSED(flame_max_temp))
148 {
149 }
150 struct Mesh *smokeModifier_do(SmokeModifierData *UNUSED(smd),
151                               Depsgraph *UNUSED(depsgraph),
152                               Scene *UNUSED(scene),
153                               Object *UNUSED(ob),
154                               Mesh *UNUSED(me))
155 {
156   return NULL;
157 }
158 float BKE_smoke_get_velocity_at(struct Object *UNUSED(ob),
159                                 float UNUSED(position[3]),
160                                 float UNUSED(velocity[3]))
161 {
162   return 0.0f;
163 }
164
165 #endif /* WITH_SMOKE */
166
167 #ifdef WITH_SMOKE
168
169 void BKE_smoke_reallocate_fluid(SmokeDomainSettings *sds, float dx, int res[3], int free_old)
170 {
171   int use_heat = (sds->active_fields & SM_ACTIVE_HEAT);
172   int use_fire = (sds->active_fields & SM_ACTIVE_FIRE);
173   int use_colors = (sds->active_fields & SM_ACTIVE_COLORS);
174
175   if (free_old && sds->fluid) {
176     smoke_free(sds->fluid);
177   }
178   if (!min_iii(res[0], res[1], res[2])) {
179     sds->fluid = NULL;
180     return;
181   }
182   sds->fluid = smoke_init(res, dx, DT_DEFAULT, use_heat, use_fire, use_colors);
183   smoke_initBlenderRNA(sds->fluid,
184                        &(sds->alpha),
185                        &(sds->beta),
186                        &(sds->time_scale),
187                        &(sds->vorticity),
188                        &(sds->border_collisions),
189                        &(sds->burning_rate),
190                        &(sds->flame_smoke),
191                        sds->flame_smoke_color,
192                        &(sds->flame_vorticity),
193                        &(sds->flame_ignition),
194                        &(sds->flame_max_temp));
195
196   /* reallocate shadow buffer */
197   if (sds->shadow) {
198     MEM_freeN(sds->shadow);
199   }
200   sds->shadow = MEM_callocN(sizeof(float) * res[0] * res[1] * res[2], "SmokeDomainShadow");
201 }
202
203 void BKE_smoke_reallocate_highres_fluid(SmokeDomainSettings *sds,
204                                         float dx,
205                                         int res[3],
206                                         int free_old)
207 {
208   int use_fire = (sds->active_fields & (SM_ACTIVE_HEAT | SM_ACTIVE_FIRE));
209   int use_colors = (sds->active_fields & SM_ACTIVE_COLORS);
210
211   if (free_old && sds->wt) {
212     smoke_turbulence_free(sds->wt);
213   }
214   if (!min_iii(res[0], res[1], res[2])) {
215     sds->wt = NULL;
216     return;
217   }
218
219   /* smoke_turbulence_init uses non-threadsafe functions from fftw3 lib (like fftw_plan & co). */
220   BLI_thread_lock(LOCK_FFTW);
221
222   sds->wt = smoke_turbulence_init(
223       res, sds->amplify + 1, sds->noise, BKE_tempdir_session(), use_fire, use_colors);
224
225   BLI_thread_unlock(LOCK_FFTW);
226
227   sds->res_wt[0] = res[0] * (sds->amplify + 1);
228   sds->res_wt[1] = res[1] * (sds->amplify + 1);
229   sds->res_wt[2] = res[2] * (sds->amplify + 1);
230   sds->dx_wt = dx / (sds->amplify + 1);
231   smoke_initWaveletBlenderRNA(sds->wt, &(sds->strength));
232 }
233
234 /* convert global position to domain cell space */
235 static void smoke_pos_to_cell(SmokeDomainSettings *sds, float pos[3])
236 {
237   mul_m4_v3(sds->imat, pos);
238   sub_v3_v3(pos, sds->p0);
239   pos[0] *= 1.0f / sds->cell_size[0];
240   pos[1] *= 1.0f / sds->cell_size[1];
241   pos[2] *= 1.0f / sds->cell_size[2];
242 }
243
244 /* set domain transformations and base resolution from object mesh */
245 static void smoke_set_domain_from_mesh(SmokeDomainSettings *sds,
246                                        Object *ob,
247                                        Mesh *me,
248                                        bool init_resolution)
249 {
250   size_t i;
251   float min[3] = {FLT_MAX, FLT_MAX, FLT_MAX}, max[3] = {-FLT_MAX, -FLT_MAX, -FLT_MAX};
252   float size[3];
253   MVert *verts = me->mvert;
254   float scale = 0.0;
255   int res;
256
257   res = sds->maxres;
258
259   // get BB of domain
260   for (i = 0; i < me->totvert; i++) {
261     // min BB
262     min[0] = MIN2(min[0], verts[i].co[0]);
263     min[1] = MIN2(min[1], verts[i].co[1]);
264     min[2] = MIN2(min[2], verts[i].co[2]);
265
266     // max BB
267     max[0] = MAX2(max[0], verts[i].co[0]);
268     max[1] = MAX2(max[1], verts[i].co[1]);
269     max[2] = MAX2(max[2], verts[i].co[2]);
270   }
271
272   /* set domain bounds */
273   copy_v3_v3(sds->p0, min);
274   copy_v3_v3(sds->p1, max);
275   sds->dx = 1.0f / res;
276
277   /* calculate domain dimensions */
278   sub_v3_v3v3(size, max, min);
279   if (init_resolution) {
280     zero_v3_int(sds->base_res);
281     copy_v3_v3(sds->cell_size, size);
282   }
283   /* apply object scale */
284   for (i = 0; i < 3; i++) {
285     size[i] = fabsf(size[i] * ob->scale[i]);
286   }
287   copy_v3_v3(sds->global_size, size);
288   copy_v3_v3(sds->dp0, min);
289
290   invert_m4_m4(sds->imat, ob->obmat);
291
292   // prevent crash when initializing a plane as domain
293   if (!init_resolution || (size[0] < FLT_EPSILON) || (size[1] < FLT_EPSILON) ||
294       (size[2] < FLT_EPSILON)) {
295     return;
296   }
297
298   /* define grid resolutions from longest domain side */
299   if (size[0] >= MAX2(size[1], size[2])) {
300     scale = res / size[0];
301     sds->scale = size[0] / fabsf(ob->scale[0]);
302     sds->base_res[0] = res;
303     sds->base_res[1] = max_ii((int)(size[1] * scale + 0.5f), 4);
304     sds->base_res[2] = max_ii((int)(size[2] * scale + 0.5f), 4);
305   }
306   else if (size[1] >= MAX2(size[0], size[2])) {
307     scale = res / size[1];
308     sds->scale = size[1] / fabsf(ob->scale[1]);
309     sds->base_res[0] = max_ii((int)(size[0] * scale + 0.5f), 4);
310     sds->base_res[1] = res;
311     sds->base_res[2] = max_ii((int)(size[2] * scale + 0.5f), 4);
312   }
313   else {
314     scale = res / size[2];
315     sds->scale = size[2] / fabsf(ob->scale[2]);
316     sds->base_res[0] = max_ii((int)(size[0] * scale + 0.5f), 4);
317     sds->base_res[1] = max_ii((int)(size[1] * scale + 0.5f), 4);
318     sds->base_res[2] = res;
319   }
320
321   /* set cell size */
322   sds->cell_size[0] /= (float)sds->base_res[0];
323   sds->cell_size[1] /= (float)sds->base_res[1];
324   sds->cell_size[2] /= (float)sds->base_res[2];
325 }
326
327 static int smokeModifier_init(SmokeModifierData *smd, Object *ob, int scene_framenr, Mesh *me)
328 {
329   if ((smd->type & MOD_SMOKE_TYPE_DOMAIN) && smd->domain && !smd->domain->fluid) {
330     SmokeDomainSettings *sds = smd->domain;
331     int res[3];
332     /* set domain dimensions from mesh */
333     smoke_set_domain_from_mesh(sds, ob, me, true);
334     /* reset domain values */
335     zero_v3_int(sds->shift);
336     zero_v3(sds->shift_f);
337     add_v3_fl(sds->shift_f, 0.5f);
338     zero_v3(sds->prev_loc);
339     mul_m4_v3(ob->obmat, sds->prev_loc);
340     copy_m4_m4(sds->obmat, ob->obmat);
341
342     /* set resolutions */
343     if (smd->domain->flags & MOD_SMOKE_ADAPTIVE_DOMAIN) {
344       res[0] = res[1] = res[2] = 1; /* use minimum res for adaptive init */
345     }
346     else {
347       copy_v3_v3_int(res, sds->base_res);
348     }
349     copy_v3_v3_int(sds->res, res);
350     sds->total_cells = sds->res[0] * sds->res[1] * sds->res[2];
351     sds->res_min[0] = sds->res_min[1] = sds->res_min[2] = 0;
352     copy_v3_v3_int(sds->res_max, res);
353
354     /* allocate fluid */
355     BKE_smoke_reallocate_fluid(sds, sds->dx, sds->res, 0);
356
357     smd->time = scene_framenr;
358
359     /* allocate highres fluid */
360     if (sds->flags & MOD_SMOKE_HIGHRES) {
361       BKE_smoke_reallocate_highres_fluid(sds, sds->dx, sds->res, 0);
362     }
363     /* allocate shadow buffer */
364     if (!sds->shadow) {
365       sds->shadow = MEM_callocN(sizeof(float) * sds->res[0] * sds->res[1] * sds->res[2],
366                                 "SmokeDomainShadow");
367     }
368     return 1;
369   }
370   else if ((smd->type & MOD_SMOKE_TYPE_FLOW) && smd->flow) {
371     smd->time = scene_framenr;
372
373     return 1;
374   }
375   else if ((smd->type & MOD_SMOKE_TYPE_COLL)) {
376     if (!smd->coll) {
377       smokeModifier_createType(smd);
378     }
379
380     smd->time = scene_framenr;
381
382     return 1;
383   }
384
385   return 2;
386 }
387
388 #endif /* WITH_SMOKE */
389
390 static void smokeModifier_freeDomain(SmokeModifierData *smd)
391 {
392   if (smd->domain) {
393     if (smd->domain->shadow) {
394       MEM_freeN(smd->domain->shadow);
395     }
396     smd->domain->shadow = NULL;
397
398     if (smd->domain->fluid) {
399       smoke_free(smd->domain->fluid);
400     }
401
402     if (smd->domain->fluid_mutex) {
403       BLI_rw_mutex_free(smd->domain->fluid_mutex);
404     }
405
406     if (smd->domain->wt) {
407       smoke_turbulence_free(smd->domain->wt);
408     }
409
410     if (smd->domain->effector_weights) {
411       MEM_freeN(smd->domain->effector_weights);
412     }
413     smd->domain->effector_weights = NULL;
414
415     if (!(smd->modifier.flag & eModifierFlag_SharedCaches)) {
416       BKE_ptcache_free_list(&(smd->domain->ptcaches[0]));
417       smd->domain->point_cache[0] = NULL;
418     }
419
420     if (smd->domain->coba) {
421       MEM_freeN(smd->domain->coba);
422     }
423
424     MEM_freeN(smd->domain);
425     smd->domain = NULL;
426   }
427 }
428
429 static void smokeModifier_freeFlow(SmokeModifierData *smd)
430 {
431   if (smd->flow) {
432     if (smd->flow->mesh) {
433       BKE_id_free(NULL, smd->flow->mesh);
434     }
435     if (smd->flow->verts_old) {
436       MEM_freeN(smd->flow->verts_old);
437     }
438     MEM_freeN(smd->flow);
439     smd->flow = NULL;
440   }
441 }
442
443 static void smokeModifier_freeCollision(SmokeModifierData *smd)
444 {
445   if (smd->coll) {
446     SmokeCollSettings *scs = smd->coll;
447
448     if (scs->numverts) {
449       if (scs->verts_old) {
450         MEM_freeN(scs->verts_old);
451         scs->verts_old = NULL;
452       }
453     }
454
455     if (smd->coll->mesh) {
456       BKE_id_free(NULL, smd->coll->mesh);
457     }
458     smd->coll->mesh = NULL;
459
460     MEM_freeN(smd->coll);
461     smd->coll = NULL;
462   }
463 }
464
465 void smokeModifier_reset_turbulence(struct SmokeModifierData *smd)
466 {
467   if (smd && smd->domain && smd->domain->wt) {
468     smoke_turbulence_free(smd->domain->wt);
469     smd->domain->wt = NULL;
470   }
471 }
472
473 static void smokeModifier_reset_ex(struct SmokeModifierData *smd, bool need_lock)
474 {
475   if (smd) {
476     if (smd->domain) {
477       if (smd->domain->shadow) {
478         MEM_freeN(smd->domain->shadow);
479       }
480       smd->domain->shadow = NULL;
481
482       if (smd->domain->fluid) {
483         if (need_lock) {
484           BLI_rw_mutex_lock(smd->domain->fluid_mutex, THREAD_LOCK_WRITE);
485         }
486
487         smoke_free(smd->domain->fluid);
488         smd->domain->fluid = NULL;
489
490         if (need_lock) {
491           BLI_rw_mutex_unlock(smd->domain->fluid_mutex);
492         }
493       }
494
495       smokeModifier_reset_turbulence(smd);
496
497       smd->time = -1;
498       smd->domain->total_cells = 0;
499       smd->domain->active_fields = 0;
500     }
501     else if (smd->flow) {
502       if (smd->flow->verts_old) {
503         MEM_freeN(smd->flow->verts_old);
504       }
505       smd->flow->verts_old = NULL;
506       smd->flow->numverts = 0;
507     }
508     else if (smd->coll) {
509       SmokeCollSettings *scs = smd->coll;
510
511       if (scs->numverts && scs->verts_old) {
512         MEM_freeN(scs->verts_old);
513         scs->verts_old = NULL;
514       }
515     }
516   }
517 }
518
519 void smokeModifier_reset(struct SmokeModifierData *smd)
520 {
521   smokeModifier_reset_ex(smd, true);
522 }
523
524 void smokeModifier_free(SmokeModifierData *smd)
525 {
526   if (smd) {
527     smokeModifier_freeDomain(smd);
528     smokeModifier_freeFlow(smd);
529     smokeModifier_freeCollision(smd);
530   }
531 }
532
533 void smokeModifier_createType(struct SmokeModifierData *smd)
534 {
535   if (smd) {
536     if (smd->type & MOD_SMOKE_TYPE_DOMAIN) {
537       if (smd->domain) {
538         smokeModifier_freeDomain(smd);
539       }
540
541       smd->domain = MEM_callocN(sizeof(SmokeDomainSettings), "SmokeDomain");
542
543       smd->domain->smd = smd;
544
545       smd->domain->point_cache[0] = BKE_ptcache_add(&(smd->domain->ptcaches[0]));
546       smd->domain->point_cache[0]->flag |= PTCACHE_DISK_CACHE;
547       smd->domain->point_cache[0]->step = 1;
548
549       /* Deprecated */
550       smd->domain->point_cache[1] = NULL;
551       BLI_listbase_clear(&smd->domain->ptcaches[1]);
552       /* set some standard values */
553       smd->domain->fluid = NULL;
554       smd->domain->fluid_mutex = BLI_rw_mutex_alloc();
555       smd->domain->wt = NULL;
556       smd->domain->eff_group = NULL;
557       smd->domain->fluid_group = NULL;
558       smd->domain->coll_group = NULL;
559       smd->domain->maxres = 32;
560       smd->domain->amplify = 1;
561       smd->domain->alpha = -0.001;
562       smd->domain->beta = 0.1;
563       smd->domain->time_scale = 1.0;
564       smd->domain->vorticity = 2.0;
565       smd->domain->border_collisions = SM_BORDER_OPEN;  // open domain
566       smd->domain->flags = MOD_SMOKE_DISSOLVE_LOG;
567       smd->domain->highres_sampling = SM_HRES_FULLSAMPLE;
568       smd->domain->strength = 2.0;
569       smd->domain->noise = MOD_SMOKE_NOISEWAVE;
570       smd->domain->diss_speed = 5;
571       smd->domain->active_fields = 0;
572
573       smd->domain->adapt_margin = 4;
574       smd->domain->adapt_res = 0;
575       smd->domain->adapt_threshold = 0.02f;
576
577       smd->domain->burning_rate = 0.75f;
578       smd->domain->flame_smoke = 1.0f;
579       smd->domain->flame_vorticity = 0.5f;
580       smd->domain->flame_ignition = 1.5f;
581       smd->domain->flame_max_temp = 3.0f;
582       /* color */
583       smd->domain->flame_smoke_color[0] = 0.7f;
584       smd->domain->flame_smoke_color[1] = 0.7f;
585       smd->domain->flame_smoke_color[2] = 0.7f;
586
587       smd->domain->viewsettings = MOD_SMOKE_VIEW_SHOW_HIGHRES;
588       smd->domain->effector_weights = BKE_effector_add_weights(NULL);
589
590 #ifdef WITH_OPENVDB_BLOSC
591       smd->domain->openvdb_comp = VDB_COMPRESSION_BLOSC;
592 #else
593       smd->domain->openvdb_comp = VDB_COMPRESSION_ZIP;
594 #endif
595       smd->domain->data_depth = 0;
596       smd->domain->cache_file_format = PTCACHE_FILE_PTCACHE;
597
598       smd->domain->display_thickness = 1.0f;
599       smd->domain->slice_method = MOD_SMOKE_SLICE_VIEW_ALIGNED;
600       smd->domain->axis_slice_method = AXIS_SLICE_FULL;
601       smd->domain->slice_per_voxel = 5.0f;
602       smd->domain->slice_depth = 0.5f;
603       smd->domain->slice_axis = 0;
604       smd->domain->vector_scale = 1.0f;
605
606       smd->domain->coba = NULL;
607       smd->domain->coba_field = FLUID_FIELD_DENSITY;
608
609       smd->domain->clipping = 1e-3f;
610     }
611     else if (smd->type & MOD_SMOKE_TYPE_FLOW) {
612       if (smd->flow) {
613         smokeModifier_freeFlow(smd);
614       }
615
616       smd->flow = MEM_callocN(sizeof(SmokeFlowSettings), "SmokeFlow");
617
618       smd->flow->smd = smd;
619
620       /* set some standard values */
621       smd->flow->density = 1.0f;
622       smd->flow->fuel_amount = 1.0f;
623       smd->flow->temp = 1.0f;
624       smd->flow->flags = MOD_SMOKE_FLOW_ABSOLUTE | MOD_SMOKE_FLOW_USE_PART_SIZE;
625       smd->flow->vel_multi = 1.0f;
626       smd->flow->volume_density = 0.0f;
627       smd->flow->surface_distance = 1.5f;
628       smd->flow->source = MOD_SMOKE_FLOW_SOURCE_MESH;
629       smd->flow->texture_size = 1.0f;
630       smd->flow->particle_size = 1.0f;
631       smd->flow->subframes = 0;
632
633       smd->flow->color[0] = 0.7f;
634       smd->flow->color[1] = 0.7f;
635       smd->flow->color[2] = 0.7f;
636
637       smd->flow->mesh = NULL;
638       smd->flow->psys = NULL;
639     }
640     else if (smd->type & MOD_SMOKE_TYPE_COLL) {
641       if (smd->coll) {
642         smokeModifier_freeCollision(smd);
643       }
644
645       smd->coll = MEM_callocN(sizeof(SmokeCollSettings), "SmokeColl");
646
647       smd->coll->smd = smd;
648       smd->coll->verts_old = NULL;
649       smd->coll->numverts = 0;
650       smd->coll->type = 0;  // static obstacle
651       smd->coll->mesh = NULL;
652     }
653   }
654 }
655
656 void smokeModifier_copy(const struct SmokeModifierData *smd,
657                         struct SmokeModifierData *tsmd,
658                         const int flag)
659 {
660   tsmd->type = smd->type;
661   tsmd->time = smd->time;
662
663   smokeModifier_createType(tsmd);
664
665   if (tsmd->domain) {
666     SmokeDomainSettings *tsds = tsmd->domain;
667     SmokeDomainSettings *sds = smd->domain;
668
669     BKE_ptcache_free_list(&(tsds->ptcaches[0]));
670
671     if (flag & LIB_ID_CREATE_NO_MAIN) {
672       /* Share the cache with the original object's modifier. */
673       tsmd->modifier.flag |= eModifierFlag_SharedCaches;
674       tsds->point_cache[0] = sds->point_cache[0];
675       tsds->ptcaches[0] = sds->ptcaches[0];
676     }
677     else {
678       tsds->point_cache[0] = BKE_ptcache_copy_list(
679           &(tsds->ptcaches[0]), &(sds->ptcaches[0]), flag);
680     }
681
682     tsds->fluid_group = sds->fluid_group;
683     tsds->coll_group = sds->coll_group;
684
685     tsds->adapt_margin = sds->adapt_margin;
686     tsds->adapt_res = sds->adapt_res;
687     tsds->adapt_threshold = sds->adapt_threshold;
688
689     tsds->alpha = sds->alpha;
690     tsds->beta = sds->beta;
691     tsds->amplify = sds->amplify;
692     tsds->maxres = sds->maxres;
693     tsds->flags = sds->flags;
694     tsds->highres_sampling = sds->highres_sampling;
695     tsds->viewsettings = sds->viewsettings;
696     tsds->noise = sds->noise;
697     tsds->diss_speed = sds->diss_speed;
698     tsds->strength = sds->strength;
699
700     tsds->border_collisions = sds->border_collisions;
701     tsds->vorticity = sds->vorticity;
702     tsds->time_scale = sds->time_scale;
703
704     tsds->burning_rate = sds->burning_rate;
705     tsds->flame_smoke = sds->flame_smoke;
706     tsds->flame_vorticity = sds->flame_vorticity;
707     tsds->flame_ignition = sds->flame_ignition;
708     tsds->flame_max_temp = sds->flame_max_temp;
709     copy_v3_v3(tsds->flame_smoke_color, sds->flame_smoke_color);
710
711     MEM_freeN(tsds->effector_weights);
712     tsds->effector_weights = MEM_dupallocN(sds->effector_weights);
713     tsds->openvdb_comp = sds->openvdb_comp;
714     tsds->data_depth = sds->data_depth;
715     tsds->cache_file_format = sds->cache_file_format;
716
717     tsds->display_thickness = sds->display_thickness;
718     tsds->slice_method = sds->slice_method;
719     tsds->axis_slice_method = sds->axis_slice_method;
720     tsds->slice_per_voxel = sds->slice_per_voxel;
721     tsds->slice_depth = sds->slice_depth;
722     tsds->slice_axis = sds->slice_axis;
723     tsds->interp_method = sds->interp_method;
724     tsds->draw_velocity = sds->draw_velocity;
725     tsds->vector_draw_type = sds->vector_draw_type;
726     tsds->vector_scale = sds->vector_scale;
727
728     tsds->use_coba = sds->use_coba;
729     tsds->coba_field = sds->coba_field;
730     if (sds->coba) {
731       tsds->coba = MEM_dupallocN(sds->coba);
732     }
733
734     tsds->clipping = sds->clipping;
735   }
736   else if (tsmd->flow) {
737     SmokeFlowSettings *tsfs = tsmd->flow;
738     SmokeFlowSettings *sfs = smd->flow;
739
740     tsfs->psys = sfs->psys;
741     tsfs->noise_texture = sfs->noise_texture;
742
743     tsfs->vel_multi = sfs->vel_multi;
744     tsfs->vel_normal = sfs->vel_normal;
745     tsfs->vel_random = sfs->vel_random;
746
747     tsfs->density = sfs->density;
748     copy_v3_v3(tsfs->color, sfs->color);
749     tsfs->fuel_amount = sfs->fuel_amount;
750     tsfs->temp = sfs->temp;
751     tsfs->volume_density = sfs->volume_density;
752     tsfs->surface_distance = sfs->surface_distance;
753     tsfs->particle_size = sfs->particle_size;
754     tsfs->subframes = sfs->subframes;
755
756     tsfs->texture_size = sfs->texture_size;
757     tsfs->texture_offset = sfs->texture_offset;
758     BLI_strncpy(tsfs->uvlayer_name, sfs->uvlayer_name, sizeof(tsfs->uvlayer_name));
759     tsfs->vgroup_density = sfs->vgroup_density;
760
761     tsfs->type = sfs->type;
762     tsfs->source = sfs->source;
763     tsfs->texture_type = sfs->texture_type;
764     tsfs->flags = sfs->flags;
765   }
766   else if (tsmd->coll) {
767     /* leave it as initialized, collision settings is mostly caches */
768   }
769 }
770
771 #ifdef WITH_SMOKE
772
773 // forward declaration
774 static void smoke_calc_transparency(SmokeDomainSettings *sds, ViewLayer *view_layer);
775 static float calc_voxel_transp(
776     float *result, float *input, int res[3], int *pixel, float *tRay, float correct);
777
778 static int get_light(ViewLayer *view_layer, float *light)
779 {
780   Base *base_tmp = NULL;
781   int found_light = 0;
782
783   // try to find a lamp, preferably local
784   for (base_tmp = FIRSTBASE(view_layer); base_tmp; base_tmp = base_tmp->next) {
785     if (base_tmp->object->type == OB_LAMP) {
786       Light *la = base_tmp->object->data;
787
788       if (la->type == LA_LOCAL) {
789         copy_v3_v3(light, base_tmp->object->obmat[3]);
790         return 1;
791       }
792       else if (!found_light) {
793         copy_v3_v3(light, base_tmp->object->obmat[3]);
794         found_light = 1;
795       }
796     }
797   }
798
799   return found_light;
800 }
801
802 /**********************************************************
803  * Obstacles
804  **********************************************************/
805
806 typedef struct ObstaclesFromDMData {
807   SmokeDomainSettings *sds;
808   const MVert *mvert;
809   const MLoop *mloop;
810   const MLoopTri *looptri;
811   BVHTreeFromMesh *tree;
812   unsigned char *obstacle_map;
813
814   bool has_velocity;
815   float *vert_vel;
816   float *velocityX, *velocityY, *velocityZ;
817   int *num_obstacles;
818 } ObstaclesFromDMData;
819
820 static void obstacles_from_mesh_task_cb(void *__restrict userdata,
821                                         const int z,
822                                         const ParallelRangeTLS *__restrict UNUSED(tls))
823 {
824   ObstaclesFromDMData *data = userdata;
825   SmokeDomainSettings *sds = data->sds;
826
827   /* slightly rounded-up sqrt(3 * (0.5)^2) == max. distance of cell boundary along the diagonal */
828   const float surface_distance = 0.867f;
829
830   for (int x = sds->res_min[0]; x < sds->res_max[0]; x++) {
831     for (int y = sds->res_min[1]; y < sds->res_max[1]; y++) {
832       const int index = smoke_get_index(
833           x - sds->res_min[0], sds->res[0], y - sds->res_min[1], sds->res[1], z - sds->res_min[2]);
834
835       float ray_start[3] = {(float)x + 0.5f, (float)y + 0.5f, (float)z + 0.5f};
836       BVHTreeNearest nearest = {0};
837       nearest.index = -1;
838       nearest.dist_sq = surface_distance *
839                         surface_distance; /* find_nearest uses squared distance */
840
841       /* find the nearest point on the mesh */
842       if (BLI_bvhtree_find_nearest(
843               data->tree->tree, ray_start, &nearest, data->tree->nearest_callback, data->tree) !=
844           -1) {
845         const MLoopTri *lt = &data->looptri[nearest.index];
846         float weights[3];
847         int v1, v2, v3;
848
849         /* calculate barycentric weights for nearest point */
850         v1 = data->mloop[lt->tri[0]].v;
851         v2 = data->mloop[lt->tri[1]].v;
852         v3 = data->mloop[lt->tri[2]].v;
853         interp_weights_tri_v3(
854             weights, data->mvert[v1].co, data->mvert[v2].co, data->mvert[v3].co, nearest.co);
855
856         // DG TODO
857         if (data->has_velocity) {
858           /* apply object velocity */
859           {
860             float hit_vel[3];
861             interp_v3_v3v3v3(hit_vel,
862                              &data->vert_vel[v1 * 3],
863                              &data->vert_vel[v2 * 3],
864                              &data->vert_vel[v3 * 3],
865                              weights);
866             data->velocityX[index] += hit_vel[0];
867             data->velocityY[index] += hit_vel[1];
868             data->velocityZ[index] += hit_vel[2];
869           }
870         }
871
872         /* tag obstacle cells */
873         data->obstacle_map[index] = 1;
874
875         if (data->has_velocity) {
876           data->obstacle_map[index] |= 8;
877           data->num_obstacles[index]++;
878         }
879       }
880     }
881   }
882 }
883
884 static void obstacles_from_mesh(Object *coll_ob,
885                                 SmokeDomainSettings *sds,
886                                 SmokeCollSettings *scs,
887                                 unsigned char *obstacle_map,
888                                 float *velocityX,
889                                 float *velocityY,
890                                 float *velocityZ,
891                                 int *num_obstacles,
892                                 float dt)
893 {
894   if (!scs->mesh) {
895     return;
896   }
897   {
898     Mesh *me = NULL;
899     MVert *mvert = NULL;
900     const MLoopTri *looptri;
901     const MLoop *mloop;
902     BVHTreeFromMesh treeData = {NULL};
903     int numverts, i;
904
905     float *vert_vel = NULL;
906     bool has_velocity = false;
907
908     me = BKE_mesh_copy_for_eval(scs->mesh, true);
909     BKE_mesh_ensure_normals(me);
910     mvert = me->mvert;
911     mloop = me->mloop;
912     looptri = BKE_mesh_runtime_looptri_ensure(me);
913     numverts = me->totvert;
914
915     // DG TODO
916     // if (scs->type > SM_COLL_STATIC)
917     // if line above is used, the code is in trouble if the object moves but is declared as "does not move"
918
919     {
920       vert_vel = MEM_callocN(sizeof(float) * numverts * 3, "smoke_obs_velocity");
921
922       if (scs->numverts != numverts || !scs->verts_old) {
923         if (scs->verts_old) {
924           MEM_freeN(scs->verts_old);
925         }
926
927         scs->verts_old = MEM_callocN(sizeof(float) * numverts * 3, "smoke_obs_verts_old");
928         scs->numverts = numverts;
929       }
930       else {
931         has_velocity = true;
932       }
933     }
934
935     /*  Transform collider vertices to
936      *   domain grid space for fast lookups */
937     for (i = 0; i < numverts; i++) {
938       float n[3];
939       float co[3];
940
941       /* vert pos */
942       mul_m4_v3(coll_ob->obmat, mvert[i].co);
943       smoke_pos_to_cell(sds, mvert[i].co);
944
945       /* vert normal */
946       normal_short_to_float_v3(n, mvert[i].no);
947       mul_mat3_m4_v3(coll_ob->obmat, n);
948       mul_mat3_m4_v3(sds->imat, n);
949       normalize_v3(n);
950       normal_float_to_short_v3(mvert[i].no, n);
951
952       /* vert velocity */
953       add_v3fl_v3fl_v3i(co, mvert[i].co, sds->shift);
954       if (has_velocity) {
955         sub_v3_v3v3(&vert_vel[i * 3], co, &scs->verts_old[i * 3]);
956         mul_v3_fl(&vert_vel[i * 3], sds->dx / dt);
957       }
958       copy_v3_v3(&scs->verts_old[i * 3], co);
959     }
960
961     if (BKE_bvhtree_from_mesh_get(&treeData, me, BVHTREE_FROM_LOOPTRI, 4)) {
962       ObstaclesFromDMData data = {
963           .sds = sds,
964           .mvert = mvert,
965           .mloop = mloop,
966           .looptri = looptri,
967           .tree = &treeData,
968           .obstacle_map = obstacle_map,
969           .has_velocity = has_velocity,
970           .vert_vel = vert_vel,
971           .velocityX = velocityX,
972           .velocityY = velocityY,
973           .velocityZ = velocityZ,
974           .num_obstacles = num_obstacles,
975       };
976       ParallelRangeSettings settings;
977       BLI_parallel_range_settings_defaults(&settings);
978       settings.scheduling_mode = TASK_SCHEDULING_DYNAMIC;
979       BLI_task_parallel_range(
980           sds->res_min[2], sds->res_max[2], &data, obstacles_from_mesh_task_cb, &settings);
981     }
982     /* free bvh tree */
983     free_bvhtree_from_mesh(&treeData);
984     BKE_id_free(NULL, me);
985
986     if (vert_vel) {
987       MEM_freeN(vert_vel);
988     }
989   }
990 }
991
992 /* Animated obstacles: dx_step = ((x_new - x_old) / totalsteps) * substep */
993 static void update_obstacles(Depsgraph *depsgraph,
994                              Object *ob,
995                              SmokeDomainSettings *sds,
996                              float dt,
997                              int UNUSED(substep),
998                              int UNUSED(totalsteps))
999 {
1000   Object **collobjs = NULL;
1001   unsigned int numcollobj = 0;
1002
1003   unsigned int collIndex;
1004   unsigned char *obstacles = smoke_get_obstacle(sds->fluid);
1005   float *velx = NULL;
1006   float *vely = NULL;
1007   float *velz = NULL;
1008   float *velxOrig = smoke_get_velocity_x(sds->fluid);
1009   float *velyOrig = smoke_get_velocity_y(sds->fluid);
1010   float *velzOrig = smoke_get_velocity_z(sds->fluid);
1011   float *density = smoke_get_density(sds->fluid);
1012   float *fuel = smoke_get_fuel(sds->fluid);
1013   float *flame = smoke_get_flame(sds->fluid);
1014   float *r = smoke_get_color_r(sds->fluid);
1015   float *g = smoke_get_color_g(sds->fluid);
1016   float *b = smoke_get_color_b(sds->fluid);
1017   unsigned int z;
1018
1019   int *num_obstacles = MEM_callocN(sizeof(int) * sds->res[0] * sds->res[1] * sds->res[2],
1020                                    "smoke_num_obstacles");
1021
1022   smoke_get_ob_velocity(sds->fluid, &velx, &vely, &velz);
1023
1024   // TODO: delete old obstacle flags
1025   for (z = 0; z < sds->res[0] * sds->res[1] * sds->res[2]; z++) {
1026     if (obstacles[z] & 8)  // Do not delete static obstacles
1027     {
1028       obstacles[z] = 0;
1029     }
1030
1031     velx[z] = 0;
1032     vely[z] = 0;
1033     velz[z] = 0;
1034   }
1035
1036   collobjs = BKE_collision_objects_create(
1037       depsgraph, ob, sds->coll_group, &numcollobj, eModifierType_Smoke);
1038
1039   // update obstacle tags in cells
1040   for (collIndex = 0; collIndex < numcollobj; collIndex++) {
1041     Object *collob = collobjs[collIndex];
1042     SmokeModifierData *smd2 = (SmokeModifierData *)modifiers_findByType(collob,
1043                                                                         eModifierType_Smoke);
1044
1045     // DG TODO: check if modifier is active?
1046
1047     if ((smd2->type & MOD_SMOKE_TYPE_COLL) && smd2->coll) {
1048       SmokeCollSettings *scs = smd2->coll;
1049       obstacles_from_mesh(collob, sds, scs, obstacles, velx, vely, velz, num_obstacles, dt);
1050     }
1051   }
1052
1053   BKE_collision_objects_free(collobjs);
1054
1055   /* obstacle cells should not contain any velocity from the smoke simulation */
1056   for (z = 0; z < sds->res[0] * sds->res[1] * sds->res[2]; z++) {
1057     if (obstacles[z]) {
1058       velxOrig[z] = 0;
1059       velyOrig[z] = 0;
1060       velzOrig[z] = 0;
1061       density[z] = 0;
1062       if (fuel) {
1063         fuel[z] = 0;
1064         flame[z] = 0;
1065       }
1066       if (r) {
1067         r[z] = 0;
1068         g[z] = 0;
1069         b[z] = 0;
1070       }
1071     }
1072     /* average velocities from multiple obstacles in one cell */
1073     if (num_obstacles[z]) {
1074       velx[z] /= num_obstacles[z];
1075       vely[z] /= num_obstacles[z];
1076       velz[z] /= num_obstacles[z];
1077     }
1078   }
1079
1080   MEM_freeN(num_obstacles);
1081 }
1082
1083 /**********************************************************
1084  * Flow emission code
1085  **********************************************************/
1086
1087 typedef struct EmissionMap {
1088   float *influence;
1089   float *influence_high;
1090   float *velocity;
1091   int min[3], max[3], res[3];
1092   int hmin[3], hmax[3], hres[3];
1093   int total_cells, valid;
1094 } EmissionMap;
1095
1096 static void em_boundInsert(EmissionMap *em, float point[3])
1097 {
1098   int i = 0;
1099   if (!em->valid) {
1100     for (; i < 3; i++) {
1101       em->min[i] = (int)floor(point[i]);
1102       em->max[i] = (int)ceil(point[i]);
1103     }
1104     em->valid = 1;
1105   }
1106   else {
1107     for (; i < 3; i++) {
1108       if (point[i] < em->min[i]) {
1109         em->min[i] = (int)floor(point[i]);
1110       }
1111       if (point[i] > em->max[i]) {
1112         em->max[i] = (int)ceil(point[i]);
1113       }
1114     }
1115   }
1116 }
1117
1118 static void clampBoundsInDomain(SmokeDomainSettings *sds,
1119                                 int min[3],
1120                                 int max[3],
1121                                 float *min_vel,
1122                                 float *max_vel,
1123                                 int margin,
1124                                 float dt)
1125 {
1126   int i;
1127   for (i = 0; i < 3; i++) {
1128     int adapt = (sds->flags & MOD_SMOKE_ADAPTIVE_DOMAIN) ? sds->adapt_res : 0;
1129     /* add margin */
1130     min[i] -= margin;
1131     max[i] += margin;
1132
1133     /* adapt to velocity */
1134     if (min_vel && min_vel[i] < 0.0f) {
1135       min[i] += (int)floor(min_vel[i] * dt);
1136     }
1137     if (max_vel && max_vel[i] > 0.0f) {
1138       max[i] += (int)ceil(max_vel[i] * dt);
1139     }
1140
1141     /* clamp within domain max size */
1142     CLAMP(min[i], -adapt, sds->base_res[i] + adapt);
1143     CLAMP(max[i], -adapt, sds->base_res[i] + adapt);
1144   }
1145 }
1146
1147 static void em_allocateData(EmissionMap *em, bool use_velocity, int hires_mul)
1148 {
1149   int i, res[3];
1150
1151   for (i = 0; i < 3; i++) {
1152     res[i] = em->max[i] - em->min[i];
1153     if (res[i] <= 0) {
1154       return;
1155     }
1156   }
1157   em->total_cells = res[0] * res[1] * res[2];
1158   copy_v3_v3_int(em->res, res);
1159
1160   em->influence = MEM_callocN(sizeof(float) * em->total_cells, "smoke_flow_influence");
1161   if (use_velocity) {
1162     em->velocity = MEM_callocN(sizeof(float) * em->total_cells * 3, "smoke_flow_velocity");
1163   }
1164
1165   /* allocate high resolution map if required */
1166   if (hires_mul > 1) {
1167     int total_cells_high = em->total_cells * (hires_mul * hires_mul * hires_mul);
1168
1169     for (i = 0; i < 3; i++) {
1170       em->hmin[i] = em->min[i] * hires_mul;
1171       em->hmax[i] = em->max[i] * hires_mul;
1172       em->hres[i] = em->res[i] * hires_mul;
1173     }
1174
1175     em->influence_high = MEM_callocN(sizeof(float) * total_cells_high,
1176                                      "smoke_flow_influence_high");
1177   }
1178   em->valid = 1;
1179 }
1180
1181 static void em_freeData(EmissionMap *em)
1182 {
1183   if (em->influence) {
1184     MEM_freeN(em->influence);
1185   }
1186   if (em->influence_high) {
1187     MEM_freeN(em->influence_high);
1188   }
1189   if (em->velocity) {
1190     MEM_freeN(em->velocity);
1191   }
1192 }
1193
1194 static void em_combineMaps(
1195     EmissionMap *output, EmissionMap *em2, int hires_multiplier, int additive, float sample_size)
1196 {
1197   int i, x, y, z;
1198
1199   /* copyfill input 1 struct and clear output for new allocation */
1200   EmissionMap em1;
1201   memcpy(&em1, output, sizeof(EmissionMap));
1202   memset(output, 0, sizeof(EmissionMap));
1203
1204   for (i = 0; i < 3; i++) {
1205     if (em1.valid) {
1206       output->min[i] = MIN2(em1.min[i], em2->min[i]);
1207       output->max[i] = MAX2(em1.max[i], em2->max[i]);
1208     }
1209     else {
1210       output->min[i] = em2->min[i];
1211       output->max[i] = em2->max[i];
1212     }
1213   }
1214   /* allocate output map */
1215   em_allocateData(output, (em1.velocity || em2->velocity), hires_multiplier);
1216
1217   /* base resolution inputs */
1218   for (x = output->min[0]; x < output->max[0]; x++) {
1219     for (y = output->min[1]; y < output->max[1]; y++) {
1220       for (z = output->min[2]; z < output->max[2]; z++) {
1221         int index_out = smoke_get_index(x - output->min[0],
1222                                         output->res[0],
1223                                         y - output->min[1],
1224                                         output->res[1],
1225                                         z - output->min[2]);
1226
1227         /* initialize with first input if in range */
1228         if (x >= em1.min[0] && x < em1.max[0] && y >= em1.min[1] && y < em1.max[1] &&
1229             z >= em1.min[2] && z < em1.max[2]) {
1230           int index_in = smoke_get_index(
1231               x - em1.min[0], em1.res[0], y - em1.min[1], em1.res[1], z - em1.min[2]);
1232
1233           /* values */
1234           output->influence[index_out] = em1.influence[index_in];
1235           if (output->velocity && em1.velocity) {
1236             copy_v3_v3(&output->velocity[index_out * 3], &em1.velocity[index_in * 3]);
1237           }
1238         }
1239
1240         /* apply second input if in range */
1241         if (x >= em2->min[0] && x < em2->max[0] && y >= em2->min[1] && y < em2->max[1] &&
1242             z >= em2->min[2] && z < em2->max[2]) {
1243           int index_in = smoke_get_index(
1244               x - em2->min[0], em2->res[0], y - em2->min[1], em2->res[1], z - em2->min[2]);
1245
1246           /* values */
1247           if (additive) {
1248             output->influence[index_out] += em2->influence[index_in] * sample_size;
1249           }
1250           else {
1251             output->influence[index_out] = MAX2(em2->influence[index_in],
1252                                                 output->influence[index_out]);
1253           }
1254           if (output->velocity && em2->velocity) {
1255             /* last sample replaces the velocity */
1256             output->velocity[index_out * 3] = ADD_IF_LOWER(output->velocity[index_out * 3],
1257                                                            em2->velocity[index_in * 3]);
1258             output->velocity[index_out * 3 + 1] = ADD_IF_LOWER(output->velocity[index_out * 3 + 1],
1259                                                                em2->velocity[index_in * 3 + 1]);
1260             output->velocity[index_out * 3 + 2] = ADD_IF_LOWER(output->velocity[index_out * 3 + 2],
1261                                                                em2->velocity[index_in * 3 + 2]);
1262           }
1263         }
1264       }  // low res loop
1265     }
1266   }
1267
1268   /* initialize high resolution input if available */
1269   if (output->influence_high) {
1270     for (x = output->hmin[0]; x < output->hmax[0]; x++) {
1271       for (y = output->hmin[1]; y < output->hmax[1]; y++) {
1272         for (z = output->hmin[2]; z < output->hmax[2]; z++) {
1273           int index_out = smoke_get_index(x - output->hmin[0],
1274                                           output->hres[0],
1275                                           y - output->hmin[1],
1276                                           output->hres[1],
1277                                           z - output->hmin[2]);
1278
1279           /* initialize with first input if in range */
1280           if (x >= em1.hmin[0] && x < em1.hmax[0] && y >= em1.hmin[1] && y < em1.hmax[1] &&
1281               z >= em1.hmin[2] && z < em1.hmax[2]) {
1282             int index_in = smoke_get_index(
1283                 x - em1.hmin[0], em1.hres[0], y - em1.hmin[1], em1.hres[1], z - em1.hmin[2]);
1284             /* values */
1285             output->influence_high[index_out] = em1.influence_high[index_in];
1286           }
1287
1288           /* apply second input if in range */
1289           if (x >= em2->hmin[0] && x < em2->hmax[0] && y >= em2->hmin[1] && y < em2->hmax[1] &&
1290               z >= em2->hmin[2] && z < em2->hmax[2]) {
1291             int index_in = smoke_get_index(
1292                 x - em2->hmin[0], em2->hres[0], y - em2->hmin[1], em2->hres[1], z - em2->hmin[2]);
1293
1294             /* values */
1295             if (additive) {
1296               output->influence_high[index_out] += em2->influence_high[index_in] * sample_size;
1297             }
1298             else {
1299               output->influence_high[index_out] = MAX2(em2->influence_high[index_in],
1300                                                        output->influence_high[index_out]);
1301             }
1302           }
1303         }  // high res loop
1304       }
1305     }
1306   }
1307
1308   /* free original data */
1309   em_freeData(&em1);
1310 }
1311
1312 typedef struct EmitFromParticlesData {
1313   SmokeFlowSettings *sfs;
1314   KDTree_3d *tree;
1315   int hires_multiplier;
1316
1317   EmissionMap *em;
1318   float *particle_vel;
1319   float hr;
1320
1321   int *min, *max, *res;
1322
1323   float solid;
1324   float smooth;
1325   float hr_smooth;
1326 } EmitFromParticlesData;
1327
1328 static void emit_from_particles_task_cb(void *__restrict userdata,
1329                                         const int z,
1330                                         const ParallelRangeTLS *__restrict UNUSED(tls))
1331 {
1332   EmitFromParticlesData *data = userdata;
1333   SmokeFlowSettings *sfs = data->sfs;
1334   EmissionMap *em = data->em;
1335   const int hires_multiplier = data->hires_multiplier;
1336
1337   for (int x = data->min[0]; x < data->max[0]; x++) {
1338     for (int y = data->min[1]; y < data->max[1]; y++) {
1339       /* take low res samples where possible */
1340       if (hires_multiplier <= 1 ||
1341           !(x % hires_multiplier || y % hires_multiplier || z % hires_multiplier)) {
1342         /* get low res space coordinates */
1343         const int lx = x / hires_multiplier;
1344         const int ly = y / hires_multiplier;
1345         const int lz = z / hires_multiplier;
1346
1347         const int index = smoke_get_index(
1348             lx - em->min[0], em->res[0], ly - em->min[1], em->res[1], lz - em->min[2]);
1349         const float ray_start[3] = {((float)lx) + 0.5f, ((float)ly) + 0.5f, ((float)lz) + 0.5f};
1350
1351         /* find particle distance from the kdtree */
1352         KDTreeNearest_3d nearest;
1353         const float range = data->solid + data->smooth;
1354         BLI_kdtree_3d_find_nearest(data->tree, ray_start, &nearest);
1355
1356         if (nearest.dist < range) {
1357           em->influence[index] = (nearest.dist < data->solid) ?
1358                                      1.0f :
1359                                      (1.0f - (nearest.dist - data->solid) / data->smooth);
1360           /* Uses particle velocity as initial velocity for smoke */
1361           if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY &&
1362               (sfs->psys->part->phystype != PART_PHYS_NO)) {
1363             madd_v3_v3fl(
1364                 &em->velocity[index * 3], &data->particle_vel[nearest.index * 3], sfs->vel_multi);
1365           }
1366         }
1367       }
1368
1369       /* take high res samples if required */
1370       if (hires_multiplier > 1) {
1371         /* get low res space coordinates */
1372         const float lx = ((float)x) * data->hr;
1373         const float ly = ((float)y) * data->hr;
1374         const float lz = ((float)z) * data->hr;
1375
1376         const int index = smoke_get_index(
1377             x - data->min[0], data->res[0], y - data->min[1], data->res[1], z - data->min[2]);
1378         const float ray_start[3] = {
1379             lx + 0.5f * data->hr, ly + 0.5f * data->hr, lz + 0.5f * data->hr};
1380
1381         /* find particle distance from the kdtree */
1382         KDTreeNearest_3d nearest;
1383         const float range = data->solid + data->hr_smooth;
1384         BLI_kdtree_3d_find_nearest(data->tree, ray_start, &nearest);
1385
1386         if (nearest.dist < range) {
1387           em->influence_high[index] = (nearest.dist < data->solid) ?
1388                                           1.0f :
1389                                           (1.0f - (nearest.dist - data->solid) / data->smooth);
1390         }
1391       }
1392     }
1393   }
1394 }
1395
1396 static void emit_from_particles(Object *flow_ob,
1397                                 SmokeDomainSettings *sds,
1398                                 SmokeFlowSettings *sfs,
1399                                 EmissionMap *em,
1400                                 Depsgraph *depsgraph,
1401                                 Scene *scene,
1402                                 float dt)
1403 {
1404   if (sfs && sfs->psys && sfs->psys->part &&
1405       ELEM(sfs->psys->part->type, PART_EMITTER, PART_FLUID))  // is particle system selected
1406   {
1407     ParticleSimulationData sim;
1408     ParticleSystem *psys = sfs->psys;
1409     float *particle_pos;
1410     float *particle_vel;
1411     int totpart = psys->totpart, totchild;
1412     int p = 0;
1413     int valid_particles = 0;
1414     int bounds_margin = 1;
1415
1416     /* radius based flow */
1417     const float solid = sfs->particle_size * 0.5f;
1418     const float smooth = 0.5f; /* add 0.5 cells of linear falloff to reduce aliasing */
1419     int hires_multiplier = 1;
1420     KDTree_3d *tree = NULL;
1421
1422     sim.depsgraph = depsgraph;
1423     sim.scene = scene;
1424     sim.ob = flow_ob;
1425     sim.psys = psys;
1426     sim.psys->lattice_deform_data = psys_create_lattice_deform_data(&sim);
1427
1428     /* prepare curvemapping tables */
1429     if ((psys->part->child_flag & PART_CHILD_USE_CLUMP_CURVE) && psys->part->clumpcurve) {
1430       curvemapping_changed_all(psys->part->clumpcurve);
1431     }
1432     if ((psys->part->child_flag & PART_CHILD_USE_ROUGH_CURVE) && psys->part->roughcurve) {
1433       curvemapping_changed_all(psys->part->roughcurve);
1434     }
1435     if ((psys->part->child_flag & PART_CHILD_USE_TWIST_CURVE) && psys->part->twistcurve) {
1436       curvemapping_changed_all(psys->part->twistcurve);
1437     }
1438
1439     /* initialize particle cache */
1440     if (psys->part->type == PART_HAIR) {
1441       // TODO: PART_HAIR not supported whatsoever
1442       totchild = 0;
1443     }
1444     else {
1445       totchild = psys->totchild * psys->part->disp / 100;
1446     }
1447
1448     particle_pos = MEM_callocN(sizeof(float) * (totpart + totchild) * 3, "smoke_flow_particles");
1449     particle_vel = MEM_callocN(sizeof(float) * (totpart + totchild) * 3, "smoke_flow_particles");
1450
1451     /* setup particle radius emission if enabled */
1452     if (sfs->flags & MOD_SMOKE_FLOW_USE_PART_SIZE) {
1453       tree = BLI_kdtree_3d_new(psys->totpart + psys->totchild);
1454
1455       /* check need for high resolution map */
1456       if ((sds->flags & MOD_SMOKE_HIGHRES) && (sds->highres_sampling == SM_HRES_FULLSAMPLE)) {
1457         hires_multiplier = sds->amplify + 1;
1458       }
1459
1460       bounds_margin = (int)ceil(solid + smooth);
1461     }
1462
1463     /* calculate local position for each particle */
1464     for (p = 0; p < totpart + totchild; p++) {
1465       ParticleKey state;
1466       float *pos;
1467       if (p < totpart) {
1468         if (psys->particles[p].flag & (PARS_NO_DISP | PARS_UNEXIST)) {
1469           continue;
1470         }
1471       }
1472       else {
1473         /* handle child particle */
1474         ChildParticle *cpa = &psys->child[p - totpart];
1475         if (psys->particles[cpa->parent].flag & (PARS_NO_DISP | PARS_UNEXIST)) {
1476           continue;
1477         }
1478       }
1479
1480       state.time = DEG_get_ctime(depsgraph); /* use depsgraph time */
1481       if (psys_get_particle_state(&sim, p, &state, 0) == 0) {
1482         continue;
1483       }
1484
1485       /* location */
1486       pos = &particle_pos[valid_particles * 3];
1487       copy_v3_v3(pos, state.co);
1488       smoke_pos_to_cell(sds, pos);
1489
1490       /* velocity */
1491       copy_v3_v3(&particle_vel[valid_particles * 3], state.vel);
1492       mul_mat3_m4_v3(sds->imat, &particle_vel[valid_particles * 3]);
1493
1494       if (sfs->flags & MOD_SMOKE_FLOW_USE_PART_SIZE) {
1495         BLI_kdtree_3d_insert(tree, valid_particles, pos);
1496       }
1497
1498       /* calculate emission map bounds */
1499       em_boundInsert(em, pos);
1500       valid_particles++;
1501     }
1502
1503     /* set emission map */
1504     clampBoundsInDomain(sds, em->min, em->max, NULL, NULL, bounds_margin, dt);
1505     em_allocateData(em, sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY, hires_multiplier);
1506
1507     if (!(sfs->flags & MOD_SMOKE_FLOW_USE_PART_SIZE)) {
1508       for (p = 0; p < valid_particles; p++) {
1509         int cell[3];
1510         size_t i = 0;
1511         size_t index = 0;
1512         int badcell = 0;
1513
1514         /* 1. get corresponding cell */
1515         cell[0] = floor(particle_pos[p * 3]) - em->min[0];
1516         cell[1] = floor(particle_pos[p * 3 + 1]) - em->min[1];
1517         cell[2] = floor(particle_pos[p * 3 + 2]) - em->min[2];
1518         /* check if cell is valid (in the domain boundary) */
1519         for (i = 0; i < 3; i++) {
1520           if ((cell[i] > em->res[i] - 1) || (cell[i] < 0)) {
1521             badcell = 1;
1522             break;
1523           }
1524         }
1525         if (badcell) {
1526           continue;
1527         }
1528         /* get cell index */
1529         index = smoke_get_index(cell[0], em->res[0], cell[1], em->res[1], cell[2]);
1530         /* Add influence to emission map */
1531         em->influence[index] = 1.0f;
1532         /* Uses particle velocity as initial velocity for smoke */
1533         if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY && (psys->part->phystype != PART_PHYS_NO)) {
1534           madd_v3_v3fl(&em->velocity[index * 3], &particle_vel[p * 3], sfs->vel_multi);
1535         }
1536       }  // particles loop
1537     }
1538     else if (valid_particles > 0) {  // MOD_SMOKE_FLOW_USE_PART_SIZE
1539       int min[3], max[3], res[3];
1540       const float hr = 1.0f / ((float)hires_multiplier);
1541       /* slightly adjust high res antialias smoothness based on number of divisions
1542        * to allow smaller details but yet not differing too much from the low res size */
1543       const float hr_smooth = smooth * powf(hr, 1.0f / 3.0f);
1544
1545       /* setup loop bounds */
1546       for (int i = 0; i < 3; i++) {
1547         min[i] = em->min[i] * hires_multiplier;
1548         max[i] = em->max[i] * hires_multiplier;
1549         res[i] = em->res[i] * hires_multiplier;
1550       }
1551
1552       BLI_kdtree_3d_balance(tree);
1553
1554       EmitFromParticlesData data = {
1555           .sfs = sfs,
1556           .tree = tree,
1557           .hires_multiplier = hires_multiplier,
1558           .hr = hr,
1559           .em = em,
1560           .particle_vel = particle_vel,
1561           .min = min,
1562           .max = max,
1563           .res = res,
1564           .solid = solid,
1565           .smooth = smooth,
1566           .hr_smooth = hr_smooth,
1567       };
1568
1569       ParallelRangeSettings settings;
1570       BLI_parallel_range_settings_defaults(&settings);
1571       settings.scheduling_mode = TASK_SCHEDULING_DYNAMIC;
1572       BLI_task_parallel_range(min[2], max[2], &data, emit_from_particles_task_cb, &settings);
1573     }
1574
1575     if (sfs->flags & MOD_SMOKE_FLOW_USE_PART_SIZE) {
1576       BLI_kdtree_3d_free(tree);
1577     }
1578
1579     /* free data */
1580     if (particle_pos) {
1581       MEM_freeN(particle_pos);
1582     }
1583     if (particle_vel) {
1584       MEM_freeN(particle_vel);
1585     }
1586   }
1587 }
1588
1589 static void sample_mesh(SmokeFlowSettings *sfs,
1590                         const MVert *mvert,
1591                         const MLoop *mloop,
1592                         const MLoopTri *mlooptri,
1593                         const MLoopUV *mloopuv,
1594                         float *influence_map,
1595                         float *velocity_map,
1596                         int index,
1597                         const int base_res[3],
1598                         float flow_center[3],
1599                         BVHTreeFromMesh *treeData,
1600                         const float ray_start[3],
1601                         const float *vert_vel,
1602                         bool has_velocity,
1603                         int defgrp_index,
1604                         MDeformVert *dvert,
1605                         float x,
1606                         float y,
1607                         float z)
1608 {
1609   float ray_dir[3] = {1.0f, 0.0f, 0.0f};
1610   BVHTreeRayHit hit = {0};
1611   BVHTreeNearest nearest = {0};
1612
1613   float volume_factor = 0.0f;
1614   float sample_str = 0.0f;
1615
1616   hit.index = -1;
1617   hit.dist = 9999;
1618   nearest.index = -1;
1619   nearest.dist_sq = sfs->surface_distance *
1620                     sfs->surface_distance; /* find_nearest uses squared distance */
1621
1622   /* Check volume collision */
1623   if (sfs->volume_density) {
1624     if (BLI_bvhtree_ray_cast(treeData->tree,
1625                              ray_start,
1626                              ray_dir,
1627                              0.0f,
1628                              &hit,
1629                              treeData->raycast_callback,
1630                              treeData) != -1) {
1631       float dot = ray_dir[0] * hit.no[0] + ray_dir[1] * hit.no[1] + ray_dir[2] * hit.no[2];
1632       /* If ray and hit face normal are facing same direction
1633        * hit point is inside a closed mesh. */
1634       if (dot >= 0) {
1635         /* Also cast a ray in opposite direction to make sure
1636          * point is at least surrounded by two faces */
1637         negate_v3(ray_dir);
1638         hit.index = -1;
1639         hit.dist = 9999;
1640
1641         BLI_bvhtree_ray_cast(
1642             treeData->tree, ray_start, ray_dir, 0.0f, &hit, treeData->raycast_callback, treeData);
1643         if (hit.index != -1) {
1644           volume_factor = sfs->volume_density;
1645         }
1646       }
1647     }
1648   }
1649
1650   /* find the nearest point on the mesh */
1651   if (BLI_bvhtree_find_nearest(
1652           treeData->tree, ray_start, &nearest, treeData->nearest_callback, treeData) != -1) {
1653     float weights[3];
1654     int v1, v2, v3, f_index = nearest.index;
1655     float n1[3], n2[3], n3[3], hit_normal[3];
1656
1657     /* emit from surface based on distance */
1658     if (sfs->surface_distance) {
1659       sample_str = sqrtf(nearest.dist_sq) / sfs->surface_distance;
1660       CLAMP(sample_str, 0.0f, 1.0f);
1661       sample_str = pow(1.0f - sample_str, 0.5f);
1662     }
1663     else {
1664       sample_str = 0.0f;
1665     }
1666
1667     /* calculate barycentric weights for nearest point */
1668     v1 = mloop[mlooptri[f_index].tri[0]].v;
1669     v2 = mloop[mlooptri[f_index].tri[1]].v;
1670     v3 = mloop[mlooptri[f_index].tri[2]].v;
1671     interp_weights_tri_v3(weights, mvert[v1].co, mvert[v2].co, mvert[v3].co, nearest.co);
1672
1673     if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY && velocity_map) {
1674       /* apply normal directional velocity */
1675       if (sfs->vel_normal) {
1676         /* interpolate vertex normal vectors to get nearest point normal */
1677         normal_short_to_float_v3(n1, mvert[v1].no);
1678         normal_short_to_float_v3(n2, mvert[v2].no);
1679         normal_short_to_float_v3(n3, mvert[v3].no);
1680         interp_v3_v3v3v3(hit_normal, n1, n2, n3, weights);
1681         normalize_v3(hit_normal);
1682         /* apply normal directional and random velocity
1683          * - TODO: random disabled for now since it doesn't really work well as pressure calc smoothens it out... */
1684         velocity_map[index * 3] += hit_normal[0] * sfs->vel_normal * 0.25f;
1685         velocity_map[index * 3 + 1] += hit_normal[1] * sfs->vel_normal * 0.25f;
1686         velocity_map[index * 3 + 2] += hit_normal[2] * sfs->vel_normal * 0.25f;
1687         /* TODO: for fire emitted from mesh surface we can use
1688          * Vf = Vs + (Ps/Pf - 1)*S to model gaseous expansion from solid to fuel */
1689       }
1690       /* apply object velocity */
1691       if (has_velocity && sfs->vel_multi) {
1692         float hit_vel[3];
1693         interp_v3_v3v3v3(
1694             hit_vel, &vert_vel[v1 * 3], &vert_vel[v2 * 3], &vert_vel[v3 * 3], weights);
1695         velocity_map[index * 3] += hit_vel[0] * sfs->vel_multi;
1696         velocity_map[index * 3 + 1] += hit_vel[1] * sfs->vel_multi;
1697         velocity_map[index * 3 + 2] += hit_vel[2] * sfs->vel_multi;
1698       }
1699     }
1700
1701     /* apply vertex group influence if used */
1702     if (defgrp_index != -1 && dvert) {
1703       float weight_mask = defvert_find_weight(&dvert[v1], defgrp_index) * weights[0] +
1704                           defvert_find_weight(&dvert[v2], defgrp_index) * weights[1] +
1705                           defvert_find_weight(&dvert[v3], defgrp_index) * weights[2];
1706       sample_str *= weight_mask;
1707     }
1708
1709     /* apply emission texture */
1710     if ((sfs->flags & MOD_SMOKE_FLOW_TEXTUREEMIT) && sfs->noise_texture) {
1711       float tex_co[3] = {0};
1712       TexResult texres;
1713
1714       if (sfs->texture_type == MOD_SMOKE_FLOW_TEXTURE_MAP_AUTO) {
1715         tex_co[0] = ((x - flow_center[0]) / base_res[0]) / sfs->texture_size;
1716         tex_co[1] = ((y - flow_center[1]) / base_res[1]) / sfs->texture_size;
1717         tex_co[2] = ((z - flow_center[2]) / base_res[2] - sfs->texture_offset) / sfs->texture_size;
1718       }
1719       else if (mloopuv) {
1720         const float *uv[3];
1721         uv[0] = mloopuv[mlooptri[f_index].tri[0]].uv;
1722         uv[1] = mloopuv[mlooptri[f_index].tri[1]].uv;
1723         uv[2] = mloopuv[mlooptri[f_index].tri[2]].uv;
1724
1725         interp_v2_v2v2v2(tex_co, UNPACK3(uv), weights);
1726
1727         /* map between -1.0f and 1.0f */
1728         tex_co[0] = tex_co[0] * 2.0f - 1.0f;
1729         tex_co[1] = tex_co[1] * 2.0f - 1.0f;
1730         tex_co[2] = sfs->texture_offset;
1731       }
1732       texres.nor = NULL;
1733       BKE_texture_get_value(NULL, sfs->noise_texture, tex_co, &texres, false);
1734       sample_str *= texres.tin;
1735     }
1736   }
1737
1738   /* multiply initial velocity by emitter influence */
1739   if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY && velocity_map) {
1740     mul_v3_fl(&velocity_map[index * 3], sample_str);
1741   }
1742
1743   /* apply final influence based on volume factor */
1744   influence_map[index] = MAX2(volume_factor, sample_str);
1745 }
1746
1747 typedef struct EmitFromDMData {
1748   SmokeDomainSettings *sds;
1749   SmokeFlowSettings *sfs;
1750   const MVert *mvert;
1751   const MLoop *mloop;
1752   const MLoopTri *mlooptri;
1753   const MLoopUV *mloopuv;
1754   MDeformVert *dvert;
1755   int defgrp_index;
1756
1757   BVHTreeFromMesh *tree;
1758   int hires_multiplier;
1759   float hr;
1760
1761   EmissionMap *em;
1762   bool has_velocity;
1763   float *vert_vel;
1764
1765   float *flow_center;
1766   int *min, *max, *res;
1767 } EmitFromDMData;
1768
1769 static void emit_from_mesh_task_cb(void *__restrict userdata,
1770                                    const int z,
1771                                    const ParallelRangeTLS *__restrict UNUSED(tls))
1772 {
1773   EmitFromDMData *data = userdata;
1774   EmissionMap *em = data->em;
1775   const int hires_multiplier = data->hires_multiplier;
1776
1777   for (int x = data->min[0]; x < data->max[0]; x++) {
1778     for (int y = data->min[1]; y < data->max[1]; y++) {
1779       /* take low res samples where possible */
1780       if (hires_multiplier <= 1 ||
1781           !(x % hires_multiplier || y % hires_multiplier || z % hires_multiplier)) {
1782         /* get low res space coordinates */
1783         const int lx = x / hires_multiplier;
1784         const int ly = y / hires_multiplier;
1785         const int lz = z / hires_multiplier;
1786
1787         const int index = smoke_get_index(
1788             lx - em->min[0], em->res[0], ly - em->min[1], em->res[1], lz - em->min[2]);
1789         const float ray_start[3] = {((float)lx) + 0.5f, ((float)ly) + 0.5f, ((float)lz) + 0.5f};
1790
1791         sample_mesh(data->sfs,
1792                     data->mvert,
1793                     data->mloop,
1794                     data->mlooptri,
1795                     data->mloopuv,
1796                     em->influence,
1797                     em->velocity,
1798                     index,
1799                     data->sds->base_res,
1800                     data->flow_center,
1801                     data->tree,
1802                     ray_start,
1803                     data->vert_vel,
1804                     data->has_velocity,
1805                     data->defgrp_index,
1806                     data->dvert,
1807                     (float)lx,
1808                     (float)ly,
1809                     (float)lz);
1810       }
1811
1812       /* take high res samples if required */
1813       if (hires_multiplier > 1) {
1814         /* get low res space coordinates */
1815         const float lx = ((float)x) * data->hr;
1816         const float ly = ((float)y) * data->hr;
1817         const float lz = ((float)z) * data->hr;
1818
1819         const int index = smoke_get_index(
1820             x - data->min[0], data->res[0], y - data->min[1], data->res[1], z - data->min[2]);
1821         const float ray_start[3] = {
1822             lx + 0.5f * data->hr,
1823             ly + 0.5f * data->hr,
1824             lz + 0.5f * data->hr,
1825         };
1826
1827         sample_mesh(data->sfs,
1828                     data->mvert,
1829                     data->mloop,
1830                     data->mlooptri,
1831                     data->mloopuv,
1832                     em->influence_high,
1833                     NULL,
1834                     index,
1835                     data->sds->base_res,
1836                     data->flow_center,
1837                     data->tree,
1838                     ray_start,
1839                     data->vert_vel,
1840                     data->has_velocity,
1841                     data->defgrp_index,
1842                     data->dvert,
1843                     /* x,y,z needs to be always lowres */
1844                     lx,
1845                     ly,
1846                     lz);
1847       }
1848     }
1849   }
1850 }
1851
1852 static void emit_from_mesh(
1853     Object *flow_ob, SmokeDomainSettings *sds, SmokeFlowSettings *sfs, EmissionMap *em, float dt)
1854 {
1855   if (sfs->mesh) {
1856     Mesh *me;
1857     int defgrp_index = sfs->vgroup_density - 1;
1858     MDeformVert *dvert = NULL;
1859     MVert *mvert = NULL;
1860     const MLoopTri *mlooptri = NULL;
1861     const MLoopUV *mloopuv = NULL;
1862     const MLoop *mloop = NULL;
1863     BVHTreeFromMesh treeData = {NULL};
1864     int numOfVerts, i;
1865     float flow_center[3] = {0};
1866
1867     float *vert_vel = NULL;
1868     int has_velocity = 0;
1869     int min[3], max[3], res[3];
1870     int hires_multiplier = 1;
1871
1872     /* copy mesh for thread safety because we modify it,
1873      * main issue is its VertArray being modified, then replaced and freed
1874      */
1875     me = BKE_mesh_copy_for_eval(sfs->mesh, true);
1876
1877     /* Duplicate vertices to modify. */
1878     if (me->mvert) {
1879       me->mvert = MEM_dupallocN(me->mvert);
1880       CustomData_set_layer(&me->vdata, CD_MVERT, me->mvert);
1881     }
1882
1883     BKE_mesh_ensure_normals(me);
1884     mvert = me->mvert;
1885     numOfVerts = me->totvert;
1886     dvert = CustomData_get_layer(&me->vdata, CD_MDEFORMVERT);
1887     mloopuv = CustomData_get_layer_named(&me->ldata, CD_MLOOPUV, sfs->uvlayer_name);
1888     mloop = me->mloop;
1889     mlooptri = BKE_mesh_runtime_looptri_ensure(me);
1890
1891     if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY) {
1892       vert_vel = MEM_callocN(sizeof(float) * numOfVerts * 3, "smoke_flow_velocity");
1893
1894       if (sfs->numverts != numOfVerts || !sfs->verts_old) {
1895         if (sfs->verts_old) {
1896           MEM_freeN(sfs->verts_old);
1897         }
1898         sfs->verts_old = MEM_callocN(sizeof(float) * numOfVerts * 3, "smoke_flow_verts_old");
1899         sfs->numverts = numOfVerts;
1900       }
1901       else {
1902         has_velocity = 1;
1903       }
1904     }
1905
1906     /*  Transform mesh vertices to
1907      *   domain grid space for fast lookups */
1908     for (i = 0; i < numOfVerts; i++) {
1909       float n[3];
1910       /* vert pos */
1911       mul_m4_v3(flow_ob->obmat, mvert[i].co);
1912       smoke_pos_to_cell(sds, mvert[i].co);
1913       /* vert normal */
1914       normal_short_to_float_v3(n, mvert[i].no);
1915       mul_mat3_m4_v3(flow_ob->obmat, n);
1916       mul_mat3_m4_v3(sds->imat, n);
1917       normalize_v3(n);
1918       normal_float_to_short_v3(mvert[i].no, n);
1919       /* vert velocity */
1920       if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY) {
1921         float co[3];
1922         add_v3fl_v3fl_v3i(co, mvert[i].co, sds->shift);
1923         if (has_velocity) {
1924           sub_v3_v3v3(&vert_vel[i * 3], co, &sfs->verts_old[i * 3]);
1925           mul_v3_fl(&vert_vel[i * 3], sds->dx / dt);
1926         }
1927         copy_v3_v3(&sfs->verts_old[i * 3], co);
1928       }
1929
1930       /* calculate emission map bounds */
1931       em_boundInsert(em, mvert[i].co);
1932     }
1933     mul_m4_v3(flow_ob->obmat, flow_center);
1934     smoke_pos_to_cell(sds, flow_center);
1935
1936     /* check need for high resolution map */
1937     if ((sds->flags & MOD_SMOKE_HIGHRES) && (sds->highres_sampling == SM_HRES_FULLSAMPLE)) {
1938       hires_multiplier = sds->amplify + 1;
1939     }
1940
1941     /* set emission map */
1942     clampBoundsInDomain(sds, em->min, em->max, NULL, NULL, (int)ceil(sfs->surface_distance), dt);
1943     em_allocateData(em, sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY, hires_multiplier);
1944
1945     /* setup loop bounds */
1946     for (i = 0; i < 3; i++) {
1947       min[i] = em->min[i] * hires_multiplier;
1948       max[i] = em->max[i] * hires_multiplier;
1949       res[i] = em->res[i] * hires_multiplier;
1950     }
1951
1952     if (BKE_bvhtree_from_mesh_get(&treeData, me, BVHTREE_FROM_LOOPTRI, 4)) {
1953       const float hr = 1.0f / ((float)hires_multiplier);
1954
1955       EmitFromDMData data = {
1956           .sds = sds,
1957           .sfs = sfs,
1958           .mvert = mvert,
1959           .mloop = mloop,
1960           .mlooptri = mlooptri,
1961           .mloopuv = mloopuv,
1962           .dvert = dvert,
1963           .defgrp_index = defgrp_index,
1964           .tree = &treeData,
1965           .hires_multiplier = hires_multiplier,
1966           .hr = hr,
1967           .em = em,
1968           .has_velocity = has_velocity,
1969           .vert_vel = vert_vel,
1970           .flow_center = flow_center,
1971           .min = min,
1972           .max = max,
1973           .res = res,
1974       };
1975
1976       ParallelRangeSettings settings;
1977       BLI_parallel_range_settings_defaults(&settings);
1978       settings.scheduling_mode = TASK_SCHEDULING_DYNAMIC;
1979       BLI_task_parallel_range(min[2], max[2], &data, emit_from_mesh_task_cb, &settings);
1980     }
1981     /* free bvh tree */
1982     free_bvhtree_from_mesh(&treeData);
1983
1984     if (vert_vel) {
1985       MEM_freeN(vert_vel);
1986     }
1987
1988     if (me->mvert) {
1989       MEM_freeN(me->mvert);
1990     }
1991     BKE_id_free(NULL, me);
1992   }
1993 }
1994
1995 /**********************************************************
1996  *  Smoke step
1997  **********************************************************/
1998
1999 static void adjustDomainResolution(SmokeDomainSettings *sds,
2000                                    int new_shift[3],
2001                                    EmissionMap *emaps,
2002                                    unsigned int numflowobj,
2003                                    float dt)
2004 {
2005   const int block_size = sds->amplify + 1;
2006   int min[3] = {32767, 32767, 32767}, max[3] = {-32767, -32767, -32767}, res[3];
2007   int total_cells = 1, res_changed = 0, shift_changed = 0;
2008   float min_vel[3], max_vel[3];
2009   int x, y, z;
2010   float *density = smoke_get_density(sds->fluid);
2011   float *fuel = smoke_get_fuel(sds->fluid);
2012   float *bigdensity = smoke_turbulence_get_density(sds->wt);
2013   float *bigfuel = smoke_turbulence_get_fuel(sds->wt);
2014   float *vx = smoke_get_velocity_x(sds->fluid);
2015   float *vy = smoke_get_velocity_y(sds->fluid);
2016   float *vz = smoke_get_velocity_z(sds->fluid);
2017   int wt_res[3];
2018
2019   if (sds->flags & MOD_SMOKE_HIGHRES && sds->wt) {
2020     smoke_turbulence_get_res(sds->wt, wt_res);
2021   }
2022
2023   INIT_MINMAX(min_vel, max_vel);
2024
2025   /* Calculate bounds for current domain content */
2026   for (x = sds->res_min[0]; x < sds->res_max[0]; x++) {
2027     for (y = sds->res_min[1]; y < sds->res_max[1]; y++) {
2028       for (z = sds->res_min[2]; z < sds->res_max[2]; z++) {
2029         int xn = x - new_shift[0];
2030         int yn = y - new_shift[1];
2031         int zn = z - new_shift[2];
2032         int index;
2033         float max_den;
2034
2035         /* skip if cell already belongs to new area */
2036         if (xn >= min[0] && xn <= max[0] && yn >= min[1] && yn <= max[1] && zn >= min[2] &&
2037             zn <= max[2]) {
2038           continue;
2039         }
2040
2041         index = smoke_get_index(x - sds->res_min[0],
2042                                 sds->res[0],
2043                                 y - sds->res_min[1],
2044                                 sds->res[1],
2045                                 z - sds->res_min[2]);
2046         max_den = (fuel) ? MAX2(density[index], fuel[index]) : density[index];
2047
2048         /* check high resolution bounds if max density isnt already high enough */
2049         if (max_den < sds->adapt_threshold && sds->flags & MOD_SMOKE_HIGHRES && sds->wt) {
2050           int i, j, k;
2051           /* high res grid index */
2052           int xx = (x - sds->res_min[0]) * block_size;
2053           int yy = (y - sds->res_min[1]) * block_size;
2054           int zz = (z - sds->res_min[2]) * block_size;
2055
2056           for (i = 0; i < block_size; i++) {
2057             for (j = 0; j < block_size; j++) {
2058               for (k = 0; k < block_size; k++) {
2059                 int big_index = smoke_get_index(xx + i, wt_res[0], yy + j, wt_res[1], zz + k);
2060                 float den = (bigfuel) ? MAX2(bigdensity[big_index], bigfuel[big_index]) :
2061                                         bigdensity[big_index];
2062                 if (den > max_den) {
2063                   max_den = den;
2064                 }
2065               }
2066             }
2067           }
2068         }
2069
2070         /* content bounds (use shifted coordinates) */
2071         if (max_den >= sds->adapt_threshold) {
2072           if (min[0] > xn) {
2073             min[0] = xn;
2074           }
2075           if (min[1] > yn) {
2076             min[1] = yn;
2077           }
2078           if (min[2] > zn) {
2079             min[2] = zn;
2080           }
2081           if (max[0] < xn) {
2082             max[0] = xn;
2083           }
2084           if (max[1] < yn) {
2085             max[1] = yn;
2086           }
2087           if (max[2] < zn) {
2088             max[2] = zn;
2089           }
2090         }
2091
2092         /* velocity bounds */
2093         if (min_vel[0] > vx[index]) {
2094           min_vel[0] = vx[index];
2095         }
2096         if (min_vel[1] > vy[index]) {
2097           min_vel[1] = vy[index];
2098         }
2099         if (min_vel[2] > vz[index]) {
2100           min_vel[2] = vz[index];
2101         }
2102         if (max_vel[0] < vx[index]) {
2103           max_vel[0] = vx[index];
2104         }
2105         if (max_vel[1] < vy[index]) {
2106           max_vel[1] = vy[index];
2107         }
2108         if (max_vel[2] < vz[index]) {
2109           max_vel[2] = vz[index];
2110         }
2111       }
2112     }
2113   }
2114
2115   /* also apply emission maps */
2116   for (int i = 0; i < numflowobj; i++) {
2117     EmissionMap *em = &emaps[i];
2118
2119     for (x = em->min[0]; x < em->max[0]; x++) {
2120       for (y = em->min[1]; y < em->max[1]; y++) {
2121         for (z = em->min[2]; z < em->max[2]; z++) {
2122           int index = smoke_get_index(
2123               x - em->min[0], em->res[0], y - em->min[1], em->res[1], z - em->min[2]);
2124           float max_den = em->influence[index];
2125
2126           /* density bounds */
2127           if (max_den >= sds->adapt_threshold) {
2128             if (min[0] > x) {
2129               min[0] = x;
2130             }
2131             if (min[1] > y) {
2132               min[1] = y;
2133             }
2134             if (min[2] > z) {
2135               min[2] = z;
2136             }
2137             if (max[0] < x) {
2138               max[0] = x;
2139             }
2140             if (max[1] < y) {
2141               max[1] = y;
2142             }
2143             if (max[2] < z) {
2144               max[2] = z;
2145             }
2146           }
2147         }
2148       }
2149     }
2150   }
2151
2152   /* calculate new bounds based on these values */
2153   mul_v3_fl(min_vel, 1.0f / sds->dx);
2154   mul_v3_fl(max_vel, 1.0f / sds->dx);
2155   clampBoundsInDomain(sds, min, max, min_vel, max_vel, sds->adapt_margin + 1, dt);
2156
2157   for (int i = 0; i < 3; i++) {
2158     /* calculate new resolution */
2159     res[i] = max[i] - min[i];
2160     total_cells *= res[i];
2161
2162     if (new_shift[i]) {
2163       shift_changed = 1;
2164     }
2165
2166     /* if no content set minimum dimensions */
2167     if (res[i] <= 0) {
2168       int j;
2169       for (j = 0; j < 3; j++) {
2170         min[j] = 0;
2171         max[j] = 1;
2172         res[j] = 1;
2173       }
2174       res_changed = 1;
2175       total_cells = 1;
2176       break;
2177     }
2178     if (min[i] != sds->res_min[i] || max[i] != sds->res_max[i]) {
2179       res_changed = 1;
2180     }
2181   }
2182
2183   if (res_changed || shift_changed) {
2184     struct FLUID_3D *fluid_old = sds->fluid;
2185     struct WTURBULENCE *turb_old = sds->wt;
2186     /* allocate new fluid data */
2187     BKE_smoke_reallocate_fluid(sds, sds->dx, res, 0);
2188     if (sds->flags & MOD_SMOKE_HIGHRES) {
2189       BKE_smoke_reallocate_highres_fluid(sds, sds->dx, res, 0);
2190     }
2191
2192     /* copy values from old fluid to new */
2193     if (sds->total_cells > 1 && total_cells > 1) {
2194       /* low res smoke */
2195       float *o_dens, *o_react, *o_flame, *o_fuel, *o_heat, *o_heatold, *o_vx, *o_vy, *o_vz, *o_r,
2196           *o_g, *o_b;
2197       float *n_dens, *n_react, *n_flame, *n_fuel, *n_heat, *n_heatold, *n_vx, *n_vy, *n_vz, *n_r,
2198           *n_g, *n_b;
2199       float dummy;
2200       unsigned char *dummy_p;
2201       /* high res smoke */
2202       int wt_res_old[3];
2203       float *o_wt_dens, *o_wt_react, *o_wt_flame, *o_wt_fuel, *o_wt_tcu, *o_wt_tcv, *o_wt_tcw,
2204           *o_wt_r, *o_wt_g, *o_wt_b;
2205       float *n_wt_dens, *n_wt_react, *n_wt_flame, *n_wt_fuel, *n_wt_tcu, *n_wt_tcv, *n_wt_tcw,
2206           *n_wt_r, *n_wt_g, *n_wt_b;
2207
2208       smoke_export(fluid_old,
2209                    &dummy,
2210                    &dummy,
2211                    &o_dens,
2212                    &o_react,
2213                    &o_flame,
2214                    &o_fuel,
2215                    &o_heat,
2216                    &o_heatold,
2217                    &o_vx,
2218                    &o_vy,
2219                    &o_vz,
2220                    &o_r,
2221                    &o_g,
2222                    &o_b,
2223                    &dummy_p);
2224       smoke_export(sds->fluid,
2225                    &dummy,
2226                    &dummy,
2227                    &n_dens,
2228                    &n_react,
2229                    &n_flame,
2230                    &n_fuel,
2231                    &n_heat,
2232                    &n_heatold,
2233                    &n_vx,
2234                    &n_vy,
2235                    &n_vz,
2236                    &n_r,
2237                    &n_g,
2238                    &n_b,
2239                    &dummy_p);
2240
2241       if (sds->flags & MOD_SMOKE_HIGHRES) {
2242         smoke_turbulence_export(turb_old,
2243                                 &o_wt_dens,
2244                                 &o_wt_react,
2245                                 &o_wt_flame,
2246                                 &o_wt_fuel,
2247                                 &o_wt_r,
2248                                 &o_wt_g,
2249                                 &o_wt_b,
2250                                 &o_wt_tcu,
2251                                 &o_wt_tcv,
2252                                 &o_wt_tcw);
2253         smoke_turbulence_get_res(turb_old, wt_res_old);
2254         smoke_turbulence_export(sds->wt,
2255                                 &n_wt_dens,
2256                                 &n_wt_react,
2257                                 &n_wt_flame,
2258                                 &n_wt_fuel,
2259                                 &n_wt_r,
2260                                 &n_wt_g,
2261                                 &n_wt_b,
2262                                 &n_wt_tcu,
2263                                 &n_wt_tcv,
2264                                 &n_wt_tcw);
2265       }
2266
2267       for (x = sds->res_min[0]; x < sds->res_max[0]; x++) {
2268         for (y = sds->res_min[1]; y < sds->res_max[1]; y++) {
2269           for (z = sds->res_min[2]; z < sds->res_max[2]; z++) {
2270             /* old grid index */
2271             int xo = x - sds->res_min[0];
2272             int yo = y - sds->res_min[1];
2273             int zo = z - sds->res_min[2];
2274             int index_old = smoke_get_index(xo, sds->res[0], yo, sds->res[1], zo);
2275             /* new grid index */
2276             int xn = x - min[0] - new_shift[0];
2277             int yn = y - min[1] - new_shift[1];
2278             int zn = z - min[2] - new_shift[2];
2279             int index_new = smoke_get_index(xn, res[0], yn, res[1], zn);
2280
2281             /* skip if outside new domain */
2282             if (xn < 0 || xn >= res[0] || yn < 0 || yn >= res[1] || zn < 0 || zn >= res[2]) {
2283               continue;
2284             }
2285
2286             /* copy data */
2287             n_dens[index_new] = o_dens[index_old];
2288             /* heat */
2289             if (n_heat && o_heat) {
2290               n_heat[index_new] = o_heat[index_old];
2291               n_heatold[index_new] = o_heatold[index_old];
2292             }
2293             /* fuel */
2294             if (n_fuel && o_fuel) {
2295               n_flame[index_new] = o_flame[index_old];
2296               n_fuel[index_new] = o_fuel[index_old];
2297               n_react[index_new] = o_react[index_old];
2298             }
2299             /* color */
2300             if (o_r && n_r) {
2301               n_r[index_new] = o_r[index_old];
2302               n_g[index_new] = o_g[index_old];
2303               n_b[index_new] = o_b[index_old];
2304             }
2305             n_vx[index_new] = o_vx[index_old];
2306             n_vy[index_new] = o_vy[index_old];
2307             n_vz[index_new] = o_vz[index_old];
2308
2309             if (sds->flags & MOD_SMOKE_HIGHRES && turb_old) {
2310               int i, j, k;
2311               /* old grid index */
2312               int xx_o = xo * block_size;
2313               int yy_o = yo * block_size;
2314               int zz_o = zo * block_size;
2315               /* new grid index */
2316               int xx_n = xn * block_size;
2317               int yy_n = yn * block_size;
2318               int zz_n = zn * block_size;
2319
2320               n_wt_tcu[index_new] = o_wt_tcu[index_old];
2321               n_wt_tcv[index_new] = o_wt_tcv[index_old];
2322               n_wt_tcw[index_new] = o_wt_tcw[index_old];
2323
2324               for (i = 0; i < block_size; i++) {
2325                 for (j = 0; j < block_size; j++) {
2326                   for (k = 0; k < block_size; k++) {
2327                     int big_index_old = smoke_get_index(
2328                         xx_o + i, wt_res_old[0], yy_o + j, wt_res_old[1], zz_o + k);
2329                     int big_index_new = smoke_get_index(
2330                         xx_n + i, sds->res_wt[0], yy_n + j, sds->res_wt[1], zz_n + k);
2331                     /* copy data */
2332                     n_wt_dens[big_index_new] = o_wt_dens[big_index_old];
2333                     if (n_wt_flame && o_wt_flame) {
2334                       n_wt_flame[big_index_new] = o_wt_flame[big_index_old];
2335                       n_wt_fuel[big_index_new] = o_wt_fuel[big_index_old];
2336                       n_wt_react[big_index_new] = o_wt_react[big_index_old];
2337                     }
2338                     if (n_wt_r && o_wt_r) {
2339                       n_wt_r[big_index_new] = o_wt_r[big_index_old];
2340                       n_wt_g[big_index_new] = o_wt_g[big_index_old];
2341                       n_wt_b[big_index_new] = o_wt_b[big_index_old];
2342                     }
2343                   }
2344                 }
2345               }
2346             }
2347           }
2348         }
2349       }
2350     }
2351     smoke_free(fluid_old);
2352     if (turb_old) {
2353       smoke_turbulence_free(turb_old);
2354     }
2355
2356     /* set new domain dimensions */
2357     copy_v3_v3_int(sds->res_min, min);
2358     copy_v3_v3_int(sds->res_max, max);
2359     copy_v3_v3_int(sds->res, res);
2360     sds->total_cells = total_cells;
2361   }
2362 }
2363
2364 BLI_INLINE void apply_outflow_fields(int index,
2365                                      float *density,
2366                                      float *heat,
2367                                      float *fuel,
2368                                      float *react,
2369                                      float *color_r,
2370                                      float *color_g,
2371                                      float *color_b)
2372 {
2373   density[index] = 0.f;
2374   if (heat) {
2375     heat[index] = 0.f;
2376   }
2377   if (fuel) {
2378     fuel[index] = 0.f;
2379     react[index] = 0.f;
2380   }
2381   if (color_r) {
2382     color_r[index] = 0.f;
2383     color_g[index] = 0.f;
2384     color_b[index] = 0.f;
2385   }
2386 }
2387
2388 BLI_INLINE void apply_inflow_fields(SmokeFlowSettings *sfs,
2389                                     float emission_value,
2390                                     int index,
2391                                     float *density,
2392                                     float *heat,
2393                                     float *fuel,
2394                                     float *react,
2395                                     float *color_r,
2396                                     float *color_g,
2397                                     float *color_b)
2398 {
2399   int absolute_flow = (sfs->flags & MOD_SMOKE_FLOW_ABSOLUTE);
2400   float dens_old = density[index];
2401   // float fuel_old = (fuel) ? fuel[index] : 0.0f;  /* UNUSED */
2402   float dens_flow = (sfs->type == MOD_SMOKE_FLOW_TYPE_FIRE) ? 0.0f : emission_value * sfs->density;
2403   float fuel_flow = emission_value * sfs->fuel_amount;
2404   /* add heat */
2405   if (heat && emission_value > 0.0f) {
2406     heat[index] = ADD_IF_LOWER(heat[index], sfs->temp);
2407   }
2408   /* absolute */
2409   if (absolute_flow) {
2410     if (sfs->type != MOD_SMOKE_FLOW_TYPE_FIRE) {
2411       if (dens_flow > density[index]) {
2412         density[index] = dens_flow;
2413       }
2414     }
2415     if (sfs->type != MOD_SMOKE_FLOW_TYPE_SMOKE && fuel && fuel_flow) {
2416       if (fuel_flow > fuel[index]) {
2417         fuel[index] = fuel_flow;
2418       }
2419     }
2420   }
2421   /* additive */
2422   else {
2423     if (sfs->type != MOD_SMOKE_FLOW_TYPE_FIRE) {
2424       density[index] += dens_flow;
2425       CLAMP(density[index], 0.0f, 1.0f);
2426     }
2427     if (sfs->type != MOD_SMOKE_FLOW_TYPE_SMOKE && fuel && sfs->fuel_amount) {
2428       fuel[index] += fuel_flow;
2429       CLAMP(fuel[index], 0.0f, 10.0f);
2430     }
2431   }
2432
2433   /* set color */
2434   if (color_r && dens_flow) {
2435     float total_dens = density[index] / (dens_old + dens_flow);
2436     color_r[index] = (color_r[index] + sfs->color[0] * dens_flow) * total_dens;
2437     color_g[index] = (color_g[index] + sfs->color[1] * dens_flow) * total_dens;
2438     color_b[index] = (color_b[index] + sfs->color[2] * dens_flow) * total_dens;
2439   }
2440
2441   /* set fire reaction coordinate */
2442   if (fuel && fuel[index] > FLT_EPSILON) {
2443     /* instead of using 1.0 for all new fuel add slight falloff
2444      * to reduce flow blockiness */
2445     float value = 1.0f - pow2f(1.0f - emission_value);
2446
2447     if (value > react[index]) {
2448       float f = fuel_flow / fuel[index];
2449       react[index] = value * f + (1.0f - f) * react[index];
2450       CLAMP(react[index], 0.0f, value);
2451     }
2452   }
2453 }
2454
2455 static void update_flowsfluids(
2456     Depsgraph *depsgraph, Scene *scene, Object *ob, SmokeDomainSettings *sds, float dt)
2457 {
2458   Object **flowobjs = NULL;
2459   EmissionMap *emaps = NULL;
2460   unsigned int numflowobj = 0;
2461   unsigned int flowIndex;
2462   int new_shift[3] = {0};
2463   int active_fields = sds->active_fields;
2464
2465   /* calculate domain shift for current frame if using adaptive domain */
2466   if (sds->flags & MOD_SMOKE_ADAPTIVE_DOMAIN) {
2467     int total_shift[3];
2468     float frame_shift_f[3];
2469     float ob_loc[3] = {0};
2470
2471     mul_m4_v3(ob->obmat, ob_loc);
2472
2473     sub_v3_v3v3(frame_shift_f, ob_loc, sds->prev_loc);
2474     copy_v3_v3(sds->prev_loc, ob_loc);
2475     /* convert global space shift to local "cell" space */
2476     mul_mat3_m4_v3(sds->imat, frame_shift_f);
2477     frame_shift_f[0] = frame_shift_f[0] / sds->cell_size[0];
2478     frame_shift_f[1] = frame_shift_f[1] / sds->cell_size[1];
2479     frame_shift_f[2] = frame_shift_f[2] / sds->cell_size[2];
2480     /* add to total shift */
2481     add_v3_v3(sds->shift_f, frame_shift_f);
2482     /* convert to integer */
2483     total_shift[0] = (int)(floorf(sds->shift_f[0]));
2484     total_shift[1] = (int)(floorf(sds->shift_f[1]));
2485     total_shift[2] = (int)(floorf(sds->shift_f[2]));
2486     sub_v3_v3v3_int(new_shift, total_shift, sds->shift);
2487     copy_v3_v3_int(sds->shift, total_shift);
2488
2489     /* calculate new domain boundary points so that smoke doesn't slide on sub-cell movement */
2490     sds->p0[0] = sds->dp0[0] - sds->cell_size[0] * (sds->shift_f[0] - total_shift[0] - 0.5f);
2491     sds->p0[1] = sds->dp0[1] - sds->cell_size[1] * (sds->shift_f[1] - total_shift[1] - 0.5f);
2492     sds->p0[2] = sds->dp0[2] - sds->cell_size[2] * (sds->shift_f[2] - total_shift[2] - 0.5f);
2493     sds->p1[0] = sds->p0[0] + sds->cell_size[0] * sds->base_res[0];
2494     sds->p1[1] = sds->p0[1] + sds->cell_size[1] * sds->base_res[1];
2495     sds->p1[2] = sds->p0[2] + sds->cell_size[2] * sds->base_res[2];
2496   }
2497
2498   flowobjs = BKE_collision_objects_create(
2499       depsgraph, ob, sds->fluid_group, &numflowobj, eModifierType_Smoke);
2500
2501   /* init emission maps for each flow */
2502   emaps = MEM_callocN(sizeof(struct EmissionMap) * numflowobj, "smoke_flow_maps");
2503
2504   /* Prepare flow emission maps */
2505   for (flowIndex = 0; flowIndex < numflowobj; flowIndex++) {
2506     Object *collob = flowobjs[flowIndex];
2507     SmokeModifierData *smd2 = (SmokeModifierData *)modifiers_findByType(collob,
2508                                                                         eModifierType_Smoke);
2509
2510     // check for initialized smoke object
2511     if ((smd2->type & MOD_SMOKE_TYPE_FLOW) && smd2->flow) {
2512       // we got nice flow object
2513       SmokeFlowSettings *sfs = smd2->flow;
2514       int subframes = sfs->subframes;
2515       EmissionMap *em = &emaps[flowIndex];
2516
2517       /* just sample flow directly to emission map if no subframes */
2518       if (!subframes) {
2519         if (sfs->source == MOD_SMOKE_FLOW_SOURCE_PARTICLES) {
2520           emit_from_particles(collob, sds, sfs, em, depsgraph, scene, dt);
2521         }
2522         else {
2523           emit_from_mesh(collob, sds, sfs, em, dt);
2524         }
2525       }
2526       /* sample subframes */
2527       else {
2528 #  if 0
2529         int scene_frame = (int)DEG_get_ctime(depsgraph);
2530 #  endif
2531         // float scene_subframe = scene->r.subframe;  // UNUSED
2532         int subframe;
2533         for (subframe = 0; subframe <= subframes; subframe++) {
2534           EmissionMap em_temp = {NULL};
2535           float sample_size = 1.0f / (float)(subframes + 1);
2536 #  if 0
2537           float prev_frame_pos = sample_size * (float)(subframe + 1);
2538 #  endif
2539           float sdt = dt * sample_size;
2540           int hires_multiplier = 1;
2541
2542           if ((sds->flags & MOD_SMOKE_HIGHRES) && (sds->highres_sampling == SM_HRES_FULLSAMPLE)) {
2543             hires_multiplier = sds->amplify + 1;
2544           }
2545
2546           /* TODO: setting the scene frame no longer works with the new depsgraph. */
2547 #  if 0
2548           /* set scene frame to match previous frame + subframe
2549            * or use current frame for last sample */
2550           if (subframe < subframes) {
2551             scene->r.cfra = scene_frame - 1;
2552             scene->r.subframe = prev_frame_pos;
2553           }
2554           else {
2555             scene->r.cfra = scene_frame;
2556             scene->r.subframe = 0.0f;
2557           }
2558 #  endif
2559
2560           if (sfs->source == MOD_SMOKE_FLOW_SOURCE_PARTICLES) {
2561             /* emit_from_particles() updates timestep internally */
2562             emit_from_particles(collob, sds, sfs, &em_temp, depsgraph, scene, sdt);
2563             if (!(sfs->flags & MOD_SMOKE_FLOW_USE_PART_SIZE)) {
2564               hires_multiplier = 1;
2565             }
2566           }
2567           else { /* MOD_SMOKE_FLOW_SOURCE_MESH */
2568             /* update flow object frame */
2569             BLI_mutex_lock(&object_update_lock);
2570             BKE_object_modifier_update_subframe(
2571                 depsgraph, scene, collob, true, 5, DEG_get_ctime(depsgraph), eModifierType_Smoke);
2572             BLI_mutex_unlock(&object_update_lock);
2573
2574             /* apply flow */
2575             emit_from_mesh(collob, sds, sfs, &em_temp, sdt);
2576           }
2577
2578           /* combine emission maps */
2579           em_combineMaps(em,
2580                          &em_temp,
2581                          hires_multiplier,
2582                          !(sfs->flags & MOD_SMOKE_FLOW_ABSOLUTE),
2583                          sample_size);
2584           em_freeData(&em_temp);
2585         }
2586       }
2587
2588       /* update required data fields */
2589       if (em->total_cells && sfs->type != MOD_SMOKE_FLOW_TYPE_OUTFLOW) {
2590         /* activate heat field if flow produces any heat */
2591         if (sfs->temp) {
2592           active_fields |= SM_ACTIVE_HEAT;
2593         }
2594         /* activate fuel field if flow adds any fuel */
2595         if (sfs->type != MOD_SMOKE_FLOW_TYPE_SMOKE && sfs->fuel_amount) {
2596           active_fields |= SM_ACTIVE_FIRE;
2597         }
2598         /* activate color field if flows add smoke with varying colors */
2599         if (sfs->type != MOD_SMOKE_FLOW_TYPE_FIRE && sfs->density) {
2600           if (!(active_fields & SM_ACTIVE_COLOR_SET)) {
2601             copy_v3_v3(sds->active_color, sfs->color);
2602             active_fields |= SM_ACTIVE_COLOR_SET;
2603           }
2604           else if (!equals_v3v3(sds->active_color, sfs->color)) {
2605             copy_v3_v3(sds->active_color, sfs->color);
2606             active_fields |= SM_ACTIVE_COLORS;
2607           }
2608         }
2609       }
2610     }
2611   }
2612
2613   /* monitor active fields based on domain settings */
2614   /* if domain has fire, activate new fields if required */
2615   if (active_fields & SM_ACTIVE_FIRE) {
2616     /* heat is always needed for fire */
2617     active_fields |= SM_ACTIVE_HEAT;
2618     /* also activate colors if domain smoke color differs from active color */
2619     if (!(active_fields & SM_ACTIVE_COLOR_SET)) {
2620       copy_v3_v3(sds->active_color, sds->flame_smoke_color);
2621       active_fields |= SM_ACTIVE_COLOR_SET;
2622     }
2623     else if (!equals_v3v3(sds->active_color, sds->flame_smoke_color)) {
2624       copy_v3_v3(sds->active_color, sds->flame_smoke_color);
2625       active_fields |= SM_ACTIVE_COLORS;
2626     }
2627   }
2628
2629   /* Adjust domain size if needed */
2630   if (sds->flags & MOD_SMOKE_ADAPTIVE_DOMAIN) {
2631     adjustDomainResolution(sds, new_shift, emaps, numflowobj, dt);
2632   }
2633
2634   /* Initialize new data fields if any */
2635   if (active_fields & SM_ACTIVE_HEAT) {
2636     smoke_ensure_heat(sds->fluid);
2637   }
2638   if (active_fields & SM_ACTIVE_FIRE) {
2639     smoke_ensure_fire(sds->fluid, sds->wt);
2640   }
2641   if (active_fields & SM_ACTIVE_COLORS) {
2642     /* initialize all smoke with "active_color" */
2643     smoke_ensure_colors(
2644         sds->fluid, sds->wt, sds->active_color[0], sds->active_color[1], sds->active_color[2]);
2645   }
2646   sds->active_fields = active_fields;
2647
2648   /* Apply emission data */
2649   if (sds->fluid) {
2650     for (flowIndex = 0; flowIndex < numflowobj; flowIndex++) {
2651       Object *collob = flowobjs[flowIndex];
2652       SmokeModifierData *smd2 = (SmokeModifierData *)modifiers_findByType(collob,
2653                                                                           eModifierType_Smoke);
2654
2655       // check for initialized smoke object
2656       if ((smd2->type & MOD_SMOKE_TYPE_FLOW) && smd2->flow) {
2657         // we got nice flow object
2658         SmokeFlowSettings *sfs = smd2->flow;
2659         EmissionMap *em = &emaps[flowIndex];
2660
2661         float *density = smoke_get_density(sds->fluid);
2662         float *color_r = smoke_get_color_r(sds->fluid);
2663         float *color_g = smoke_get_color_g(sds->fluid);
2664         float *color_b = smoke_get_color_b(sds->fluid);
2665         float *fuel = smoke_get_fuel(sds->fluid);
2666         float *react = smoke_get_react(sds->fluid);
2667         float *bigdensity = smoke_turbulence_get_density(sds->wt);
2668         float *bigfuel = smoke_turbulence_get_fuel(sds->wt);
2669         float *bigreact = smoke_turbulence_get_react(sds->wt);
2670         float *bigcolor_r = smoke_turbulence_get_color_r(sds->wt);
2671         float *bigcolor_g = smoke_turbulence_get_color_g(sds->wt);
2672         float *bigcolor_b = smoke_turbulence_get_color_b(sds->wt);
2673         float *heat = smoke_get_heat(sds->fluid);
2674         float *velocity_x = smoke_get_velocity_x(sds->fluid);
2675         float *velocity_y = smoke_get_velocity_y(sds->fluid);
2676         float *velocity_z = smoke_get_velocity_z(sds->fluid);
2677         //unsigned char *obstacle = smoke_get_obstacle(sds->fluid);
2678         // DG TODO UNUSED unsigned char *obstacleAnim = smoke_get_obstacle_anim(sds->fluid);
2679         int bigres[3];
2680         float *velocity_map = em->velocity;
2681         float *emission_map = em->influence;
2682         float *emission_map_high = em->influence_high;
2683
2684         int ii, jj, kk, gx, gy, gz, ex, ey, ez, dx, dy, dz, block_size;
2685         size_t e_index, d_index, index_big;
2686
2687         // loop through every emission map cell
2688         for (gx = em->min[0]; gx < em->max[0]; gx++) {
2689           for (gy = em->min[1]; gy < em->max[1]; gy++) {
2690             for (gz = em->min[2]; gz < em->max[2]; gz++) {
2691               /* get emission map index */
2692               ex = gx - em->min[0];
2693               ey = gy - em->min[1];
2694               ez = gz - em->min[2];
2695               e_index = smoke_get_index(ex, em->res[0], ey, em->res[1], ez);
2696
2697               /* get domain index */
2698               dx = gx - sds->res_min[0];
2699               dy = gy - sds->res_min[1];
2700               dz = gz - sds->res_min[2];
2701               d_index = smoke_get_index(dx, sds->res[0], dy, sds->res[1], dz);
2702               /* make sure emission cell is inside the new domain boundary */
2703               if (dx < 0 || dy < 0 || dz < 0 || dx >= sds->res[0] || dy >= sds->res[1] ||
2704                   dz >= sds->res[2]) {
2705                 continue;
2706               }
2707
2708               if (sfs->type == MOD_SMOKE_FLOW_TYPE_OUTFLOW) {  // outflow
2709                 apply_outflow_fields(
2710                     d_index, density, heat, fuel, react, color_r, color_g, color_b);
2711               }
2712               else {  // inflow
2713                 apply_inflow_fields(sfs,
2714                                     emission_map[e_index],
2715                                     d_index,
2716                                     density,
2717                                     heat,
2718                                     fuel,
2719                                     react,
2720                                     color_r,
2721                                     color_g,
2722                                     color_b);
2723
2724                 /* initial velocity */
2725                 if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY) {
2726                   velocity_x[d_index] = ADD_IF_LOWER(velocity_x[d_index],
2727                                                      velocity_map[e_index * 3]);
2728                   velocity_y[d_index] = ADD_IF_LOWER(velocity_y[d_index],
2729                                                      velocity_map[e_index * 3 + 1]);
2730                   velocity_z[d_index] = ADD_IF_LOWER(velocity_z[d_index],
2731                                                      velocity_map[e_index * 3 + 2]);
2732                 }
2733               }
2734
2735               /* loop through high res blocks if high res enabled */
2736               if (bigdensity) {
2737                 // neighbor cell emission densities (for high resolution smoke smooth interpolation)
2738                 float c000, c001, c010, c011, c100, c101, c110, c111;
2739
2740                 smoke_turbulence_get_res(sds->wt, bigres);
2741                 block_size = sds->amplify + 1;  // high res block size
2742
2743                 c000 = (ex > 0 && ey > 0 && ez > 0) ?
2744                            emission_map[smoke_get_index(
2745                                ex - 1, em->res[0], ey - 1, em->res[1], ez - 1)] :
2746                            0;
2747                 c001 =
2748                     (ex > 0 && ey > 0) ?
2749                         emission_map[smoke_get_index(ex - 1, em->res[0], ey - 1, em->res[1], ez)] :
2750                         0;
2751                 c010 =
2752                     (ex > 0 && ez > 0) ?
2753                         emission_map[smoke_get_index(ex - 1, em->res[0], ey, em->res[1], ez - 1)] :
2754                         0;
2755                 c011 = (ex > 0) ?
2756                            emission_map[smoke_get_index(ex - 1, em->res[0], ey, em->res[1], ez)] :
2757                            0;
2758
2759                 c100 =
2760                     (ey > 0 && ez > 0) ?
2761                         emission_map[smoke_get_index(ex, em->res[0], ey - 1, em->res[1], ez - 1)] :
2762                         0;
2763                 c101 = (ey > 0) ?
2764                            emission_map[smoke_get_index(ex, em->res[0], ey - 1, em->res[1], ez)] :
2765                            0;
2766                 c110 = (ez > 0) ?
2767                            emission_map[smoke_get_index(ex, em->res[0], ey, em->res[1], ez - 1)] :
2768                            0;
2769                 c111 = emission_map[smoke_get_index(
2770                     ex, em->res[0], ey, em->res[1], ez)];  // this cell
2771
2772                 for (ii = 0; ii < block_size; ii++) {
2773                   for (jj = 0; jj < block_size; jj++) {
2774                     for (kk = 0; kk < block_size; kk++) {
2775
2776                       float fx, fy, fz, interpolated_value;
2777                       int shift_x = 0, shift_y = 0, shift_z = 0;
2778
2779                       /* Use full sample emission map if enabled and available */
2780                       if ((sds->highres_sampling == SM_HRES_FULLSAMPLE) && emission_map_high) {
2781                         interpolated_value =
2782                             emission_map_high[smoke_get_index(ex * block_size + ii,
2783                                                               em->res[0] * block_size,
2784                                                               ey * block_size + jj,
2785                                                               em->res[1] * block_size,
2786                                                               ez * block_size + kk)];  // this cell
2787                       }
2788                       else if (sds->highres_sampling == SM_HRES_NEAREST) {
2789                         /* without interpolation use same low resolution
2790                          * block value for all hi-res blocks */
2791                         interpolated_value = c111;
2792                       }
2793                       /* Fall back to interpolated */
2794                       else {
2795                         /* get relative block position
2796                          * for interpolation smoothing */
2797                         fx = (float)ii / block_size + 0.5f / block_size;
2798                         fy = (float)jj / block_size + 0.5f / block_size;
2799                         fz = (float)kk / block_size + 0.5f / block_size;
2800
2801                         /* calculate trilinear interpolation */
2802                         interpolated_value = c000 * (1 - fx) * (1 - fy) * (1 - fz) +
2803                                              c100 * fx * (1 - fy) * (1 - fz) +
2804                                              c010 * (1 - fx) * fy * (1 - fz) +
2805                                              c001 * (1 - fx) * (1 - fy) * fz +
2806                                              c101 * fx * (1 - fy) * fz +
2807                                              c011 * (1 - fx) * fy * fz +
2808                                              c110 * fx * fy * (1 - fz) + c111 * fx * fy * fz;
2809
2810                         /* add some contrast / sharpness
2811                          * depending on hi-res block size */
2812                         interpolated_value = (interpolated_value - 0.4f) * (block_size / 2) + 0.4f;
2813                         CLAMP(interpolated_value, 0.0f, 1.0f);
2814
2815                         /* shift smoke block index
2816                          * (because pixel center is actually
2817                          * in halfway of the low res block) */
2818                         shift_x = (dx < 1) ? 0 : block_size / 2;
2819                         shift_y = (dy < 1) ? 0 : block_size / 2;
2820                         shift_z = (dz < 1) ? 0 : block_size / 2;
2821                       }
2822
2823                       /* get shifted index for current high resolution block */
2824                       index_big = smoke_get_index(block_size * dx + ii - shift_x,
2825                                                   bigres[0],
2826                                                   block_size * dy + jj - shift_y,
2827                                                   bigres[1],
2828                                                   block_size * dz + kk - shift_z);
2829
2830                       if (sfs->type == MOD_SMOKE_FLOW_TYPE_OUTFLOW) {  // outflow
2831                         if (interpolated_value) {
2832                           apply_outflow_fields(index_big,
2833                                                bigdensity,
2834                                                NULL,
2835                                                bigfuel,
2836                                                bigreact,
2837                                                bigcolor_r,
2838                                                bigcolor_g,
2839                                                bigcolor_b);
2840                         }
2841                       }
2842                       else {  // inflow
2843                         apply_inflow_fields(sfs,
2844                                             interpolated_value,
2845                                             index_big,
2846                                             bigdensity,
2847                                             NULL,
2848                                             bigfuel,
2849                                             bigreact,
2850                                             bigcolor_r,
2851                                             bigcolor_g,
2852                                             bigcolor_b);
2853                       }
2854                     }  // hires loop
2855                   }
2856                 }
2857               }  // bigdensity
2858             }    // low res loop
2859           }
2860         }
2861
2862         // free emission maps
2863         em_freeData(em);
2864
2865       }  // end emission
2866     }
2867   }
2868
2869   BKE_collision_objects_free(flowobjs);
2870   if (emaps) {
2871     MEM_freeN(emaps);
2872   }
2873 }
2874
2875 typedef struct UpdateEffectorsData {
2876   Scene *scene;
2877   SmokeDomainSettings *sds;
2878   ListBase *effectors;
2879
2880   float *density;
2881   float *fuel;
2882   float *force_x;
2883   float *force_y;
2884   float *force_z;
2885   float *velocity_x;
2886   float *velocity_y;
2887   float *velocity_z;
2888   unsigned char *obstacle;
2889 } UpdateEffectorsData;
2890
2891 static void update_effectors_task_cb(void *__restrict userdata,
2892                                      const int x,
2893                                      const ParallelRangeTLS *__restrict UNUSED(tls))
2894 {
2895   UpdateEffectorsData *data = userdata;
2896   SmokeDomainSettings *sds = data->sds;
2897
2898   for (int y = 0; y < sds->res[1]; y++) {
2899     for (int z = 0; z < sds->res[2]; z++) {
2900       EffectedPoint epoint;
2901       float mag;
2902       float voxelCenter[3] = {0, 0, 0}, vel[3] = {0, 0, 0}, retvel[3] = {0, 0, 0};
2903       const unsigned int index = smoke_get_index(x, sds->res[0], y, sds->res[1], z);
2904
2905       if (((data->fuel ? MAX2(data->density[index], data->fuel[index]) : data->density[index]) <
2906            FLT_EPSILON) ||
2907           data->obstacle[index]) {
2908         continue;
2909       }
2910
2911       vel[0] = data->velocity_x[index];
2912       vel[1] = data->velocity_y[index];
2913       vel[2] = data->velocity_z[index];
2914
2915       /* convert vel to global space */
2916       mag = len_v3(vel);
2917       mul_mat3_m4_v3(sds->obmat, vel);
2918       normalize_v3(vel);
2919       mul_v3_fl(vel, mag);
2920
2921       voxelCenter[0] = sds->p0[0] + sds->cell_size[0] * ((float)(x + sds->res_min[0]) + 0.5f);
2922       voxelCenter[1] = sds->p0[1] + sds->cell_size[1] * ((float)(y + sds->res_min[1]) + 0.5f);
2923       voxelCenter[2] = sds->p0[2] + sds->cell_size[2] * ((float)(z + sds->res_min[2]) + 0.5f);
2924       mul_m4_v3(sds->obmat, voxelCenter);
2925
2926       pd_point_from_loc(data->scene, voxelCenter, vel, index, &epoint);
2927       BKE_effectors_apply(data->effectors, NULL, sds->effector_weights, &epoint, retvel, NULL);
2928
2929       /* convert retvel to local space */
2930       mag = len_v3(retvel);
2931       mul_mat3_m4_v3(sds->imat, retvel);
2932       normalize_v3(retvel);
2933       mul_v3_fl(retvel, mag);
2934
2935       // TODO dg - do in force!
2936       data->force_x[index] = min_ff(max_ff(-1.0f, retvel[0] * 0.2f), 1.0f);
2937       data->force_y[index] = min_ff(max_ff(-1.0f, retvel[1] * 0.2f), 1.0f);
2938       data->force_z[index] = min_ff(max_ff(-1.0f, retvel[2] * 0.2f), 1.0f);
2939     }
2940   }
2941 }
2942
2943 static void update_effectors(
2944     Depsgraph *depsgraph, Scene *scene, Object *ob, SmokeDomainSettings *sds, float UNUSED(dt))
2945 {
2946   ListBase *effectors;
2947   /* make sure smoke flow influence is 0.0f */
2948   sds->effector_weights->weight[PFIELD_SMOKEFLOW] = 0.0f;
2949   effectors = BKE_effectors_create(depsgraph, ob, NULL, sds->effector_weights);
2950
2951   if (effectors) {
2952     // precalculate wind forces
2953     UpdateEffectorsData data;
2954     data.scene = scene;
2955     data.sds = sds;
2956     data.effectors = effectors;
2957     data.density = smoke_get_density(sds->fluid);
2958     data.fuel = smoke_get_fuel(sds->fluid);
2959     data.force_x = smoke_get_force_x(sds->fluid);
2960     data.force_y = smoke_get_force_y(sds->fluid);
2961     data.force_z = smoke_get_force_z(sds->fluid);
2962     data.velocity_x = smoke_get_velocity_x(sds->fluid);
2963     data.velocity_y = smoke_get_velocity_y(sds->fluid);
2964     data.velocity_z = smoke_get_velocity_z(sds->fluid);
2965     data.obstacle = smoke_get_obstacle(sds->fluid);
2966
2967     ParallelRangeSettings settings;
2968     BLI_parallel_range_settings_defaults(&settings);
2969     settings.scheduling_mode = TASK_SCHEDULING_DYNAMIC;
2970     BLI_task_parallel_range(0, sds->res[0], &data, update_effectors_task_cb, &settings);
2971   }
2972
2973   BKE_effectors_free(effectors);
2974 }
2975
2976 static void step(Depsgraph *depsgraph,
2977                  Scene *scene,
2978                  Object *ob,
2979                  SmokeModifierData *smd,
2980                  Mesh *domain_me,
2981                  float fps)
2982 {
2983   SmokeDomainSettings *sds = smd->domain;
2984   /* stability values copied from wturbulence.cpp */
2985   const int maxSubSteps = 25;
2986   float maxVel;
2987   // maxVel should be 1.5 (1.5 cell max movement) * dx (cell size)
2988
2989   float dt;
2990   float maxVelMag = 0.0f;
2991   int totalSubsteps;
2992   int substep = 0;
2993   float dtSubdiv;
2994   float gravity[3] = {0.0f, 0.0f, -1.0f};
2995   float gravity_mag;
2996
2997   /* update object state */
2998   invert_m4_m4(sds->imat, ob->obmat);
2999   copy_m4_m4(sds->obmat, ob->obmat);
3000   smoke_set_domain_from_mesh(sds, ob, domain_me, (sds->flags & MOD_SMOKE_ADAPTIVE_DOMAIN) != 0);
3001
3002   /* use global gravity if enabled */
3003   if (scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY) {
3004     copy_v3_v3(gravity, scene->physics_settings.gravity);
3005     /* map default value to 1.0 */
3006     mul_v3_fl(gravity, 1.0f / 9.810f);
3007   }
3008   /* convert gravity to domain space */
3009   gravity_mag = len_v3(gravity);
3010   mul_mat3_m4_v3(sds->imat, gravity);
3011   normalize_v3(gravity);
3012   mul_v3_fl(gravity, gravity_mag);
3013
3014   /* adapt timestep for different framerates, dt = 0.1 is at 25fps */
3015   dt = DT_DEFAULT * (25.0f / fps);
3016   // maximum timestep/"CFL" constraint: dt < 5.0 *dx / maxVel
3017   maxVel = (sds->dx * 5.0f);
3018
3019   maxVelMag = sqrtf(maxVelMag) * dt * sds->time_scale;
3020   totalSubsteps = (int)((maxVelMag / maxVel) + 1.0f); /* always round up */
3021   totalSubsteps = (totalSubsteps < 1) ? 1 : totalSubsteps;
3022   totalSubsteps = (totalSubsteps > maxSubSteps) ? maxSubSteps : totalSubsteps;
3023
3024   /* Disable substeps for now, since it results in numerical instability */
3025   totalSubsteps = 1.0f;
3026
3027   dtSubdiv = (float)dt / (float)totalSubsteps;
3028
3029   // printf("totalSubsteps: %d, maxVelMag: %f, dt: %f\n", totalSubsteps, maxVelMag, dt);
3030
3031   for (substep = 0; substep < totalSubsteps; substep++) {
3032     // calc animated obstacle velocities
3033     update_flowsfluids(depsgraph, scene, ob, sds, dtSubdiv);
3034     update_obstacles(depsgraph, ob, sds, dtSubdiv, substep, totalSubsteps);
3035
3036     if (sds->total_cells > 1) {
3037       update_effectors(
3038           depsgraph,
3039           scene,
3040           ob,
3041           sds,
3042           dtSubdiv);  // DG TODO? problem --> uses forces instead of velocity, need to check how they need to be changed with variable dt
3043       smoke_step(sds->fluid, gravity, dtSubdiv);
3044     }
3045   }
3046 }
3047
3048 static Mesh *createDomainGeometry(SmokeDomainSettings *sds, Object *ob)
3049 {
3050   Mesh *result;
3051   MVert *mverts;
3052   MPoly *mpolys;
3053   MLoop *mloops;
3054   float min[3];
3055   float max[3];
3056   float *co;
3057   MPoly *mp;
3058   MLoop *ml;
3059
3060   int num_verts = 8;
3061   int num_faces = 6;
3062   int i;
3063   float ob_loc[3] = {0};
3064   float ob_cache_loc[3] = {0};
3065
3066   /* dont generate any mesh if there isnt any content */
3067   if (sds->total_cells <= 1) {
3068     num_verts = 0;
3069     num_faces = 0;
3070   }
3071
3072   result = BKE_mesh_new_nomain(num_verts, 0, 0, num_faces * 4, num_faces);
3073   mverts = result->mvert;
3074   mpolys = result->mpoly;
3075   mloops = result->mloop;
3076
3077   if (num_verts) {
3078     /* volume bounds */
3079     madd_v3fl_v3fl_v3fl_v3i(min, sds->p0, sds->cell_size, sds->res_min);
3080     madd_v3fl_v3fl_v3fl_v3i(max, sds->p0, sds->cell_size, sds->res_max);
3081
3082     /* set vertices */
3083     /* top slab */
3084     co = mverts[0].co;
3085     co[0] = min[0];
3086     co[1] = min[1];
3087     co[2] = max[2];
3088     co = mverts[1].co;
3089     co[0] = max[0];
3090     co[1] = min[1];
3091     co[2] = max[2];
3092     co = mverts[2].co;
3093     co[0] = max[0];
3094     co[1] = max[1];
3095     co[2] = max[2];
3096     co = mverts[3].co;
3097     co[0] = min[0];
3098     co[1] = max[1];
3099     co[2] = max[2];
3100     /* bottom slab */
3101     co = mverts[4].co;
3102     co[0] = min[0];
3103     co[1] = min[1];
3104     co[2] = min[2];
3105     co = mverts[5].co;
3106     co[0] = max[0];
3107     co[1] = min[1];
3108     co[2] = min[2];
3109     co = mverts[6].co;
3110     co[0] = max[0];
3111     co[1] = max[1];
3112     co[2] = min[2];
3113     co = mverts[7].co;
3114     co[0] = min[0];
3115     co[1] = max[1];
3116     co[2] = min[2];
3117
3118     /* create faces */
3119     /* top */
3120     mp = &mpolys[0];
3121     ml = &mloops[0 * 4];
3122     mp->loopstart = 0 * 4;
3123     mp->totloop = 4;
3124     ml[0].v = 0;
3125     ml[1].v = 1;
3126     ml[2].v = 2;
3127     ml[3].v = 3;
3128     /* right */
3129     mp = &mpolys[1];
3130     ml = &mloops[1 * 4];
3131     mp->loopstart = 1 * 4;
3132     mp->totloop = 4;
3133     ml[0].v = 2;
3134     ml[1].v = 1;
3135     ml[2].v = 5;
3136     ml[3].v = 6;
3137     /* bottom */
3138     mp = &mpolys[2];
3139     ml = &mloops[2 * 4];
3140     mp->loopstart = 2 * 4;
3141     mp->totloop = 4;
3142     ml[0].v = 7;
3143     ml[1].v = 6;
3144     ml[2].v = 5;
3145     ml[3].v = 4;
3146     /* left */
3147     mp = &mpolys[3];
3148     ml = &mloops[3 * 4];
3149     mp->loopstart = 3 * 4;
3150     mp->totloop = 4;
3151     ml[0].v = 0;
3152     ml[1].v = 3;
3153     ml[2].v = 7;
3154     ml[3].v = 4;
3155     /* front */
3156     mp = &mpolys[4];
3157     ml = &mloops[4 * 4];
3158     mp->loopstart = 4 * 4;
3159     mp->totloop = 4;
3160     ml[0].v = 3;
3161     ml[1].v = 2;
3162     ml[2].v = 6;
3163     ml[3].v = 7;
3164     /* back */
3165     mp = &mpolys[5];
3166     ml = &mloops[5 * 4];
3167     mp->loopstart = 5 * 4;
3168     mp->totloop = 4;
3169     ml[0].v = 1;
3170     ml[1].v = 0;
3171     ml[2].v = 4;
3172     ml[3].v = 5;
3173
3174     /* calculate required shift to match domain's global position
3175      * it was originally simulated at (if object moves without smoke step) */
3176     invert_m4_m4(ob->imat, ob->obmat);
3177     mul_m4_v3(ob->obmat, ob_loc);
3178     mul_m4_v3(sds->obmat, ob_cache_loc);
3179     sub_v3_v3v3(sds->obj_shift_f, ob_cache_loc, ob_loc);
3180     /* convert shift to local space and apply to vertices */
3181     mul_mat3_m4_v3(ob->imat, sds->obj_shift_f);
3182     /* apply */
3183     for (i = 0; i < num_verts; i++) {
3184       add_v3_v3(mverts[i].co, sds->obj_shift_f);
3185     }
3186   }
3187
3188   BKE_mesh_calc_edges(result, false, false);
3189   result->runtime.cd_dirty_vert |= CD_MASK_NORMAL;
3190   return result;
3191 }
3192
3193 static void smokeModifier_process(
3194     SmokeModifierData *smd, Depsgraph *depsgraph, Scene *scene, Object *ob, Mesh *me)
3195 {
3196   const int scene_framenr = (int)DEG_get_ctime(depsgraph);
3197
3198   if ((smd->type & MOD_SMOKE_TYPE_FLOW)) {
3199     if (scene_framenr >= smd->time) {
3200       smokeModifier_init(smd, ob, scene_framenr, me);
3201     }
3202
3203     if (smd->flow->mesh) {
3204       BKE_id_free(NULL, smd->flow->mesh);
3205     }
3206     smd->flow->mesh = BKE_mesh_copy_for_eval(me, false);
3207
3208     if (scene_framenr > smd->time) {
3209       smd->time = scene_framenr;
3210     }
3211     else if (scene_framenr < smd->time) {
3212       smd->time = scene_framenr;
3213       smokeModifier_reset_ex(smd, false);
3214     }
3215   }
3216   else if (smd->type & MOD_SMOKE_TYPE_COLL) {
3217     if (scene_framenr >= smd->time) {
3218       smokeModifier_init(smd, ob, scene_framenr, me);
3219     }
3220
3221     if (smd->coll) {
3222       if (smd->coll->mesh) {
3223         BKE_id_free(NULL, smd->coll->mesh);
3224       }
3225
3226       smd->coll->mesh = BKE_mesh_copy_for_eval(me, false);
3227     }
3228
3229     smd->time = scene_framenr;
3230     if (scene_framenr < smd->time) {
3231       smokeModifier_reset_ex(smd, false);
3232     }
3233   }
3234   else if (smd->type & MOD_SMOKE_TYPE_DOMAIN) {
3235     SmokeDomainSettings *sds = smd->domain;
3236     PointCache *cache = NULL;
3237     PTCacheID pid;
3238     int startframe, endframe, framenr;
3239     float timescale;
3240
3241     framenr = scene_framenr;
3242
3243     cache = sds->point_cache[0];
3244     BKE_ptcache_id_from_smoke(&pid, ob, smd);
3245     BKE_ptcache_id_time(&pid, scene, framenr, &startframe, &endframe, &timescale);
3246
3247     if (!smd->domain->fluid || framenr == startframe) {
3248       BKE_ptcache_id_reset(scene, &pid, PTCACHE_RESET_OUTDATED);
3249       smokeModifier_reset_ex(smd, false);
3250       BKE_ptcache_validate(cache, framenr);
3251       cache->flag &= ~PTCACHE_REDO_NEEDED;
3252     }
3253
3254     if (!smd->domain->fluid && (framenr != startframe) &&
3255         (smd->domain->flags & MOD_SMOKE_FILE_LOAD) == 0 && (cache->flag & PTCACHE_BAKED) == 0) {
3256       return;
3257     }
3258
3259     smd->domain->flags &= ~MOD_SMOKE_FILE_LOAD;
3260     CLAMP(framenr, startframe, endframe);
3261
3262     /* If already viewing a pre/after frame, no need to reload */
3263     if ((smd->time == framenr) && (framenr != scene_framenr)) {
3264       return;
3265     }
3266
3267     if (smokeModifier_init(smd, ob, scene_framenr, me) == 0) {
3268       printf("bad smokeModifier_init\n");
3269       return;
3270     }
3271
3272     /* only calculate something when we advanced a single frame */
3273     /* don't simulate if viewing start frame, but scene frame is not real start frame */
3274     bool can_simulate = (framenr == (int)smd->time + 1) && (framenr == scene_framenr);
3275
3276     /* try to read from cache */
3277     if (BKE_ptcache_read(&pid, (float)framenr, can_simulate) == PTCACHE_READ_EXACT) {
3278       BKE_ptcache_validate(cache, framenr);
3279       smd->time = framenr;
3280       return;
3281     }
3282
3283     if (!can_simulate) {
3284       return;
3285     }
3286
3287 #  ifdef DEBUG_TIME
3288     double start = PIL_check_seconds_timer();
3289 #  endif
3290
3291     /* if on second frame, write cache for first frame */
3292     if ((int)smd->time == startframe &&
3293         (cache->flag & PTCACHE_OUTDATED || cache->last_exact == 0)) {
3294       BKE_ptcache_write(&pid, startframe);
3295     }
3296
3297     // set new time
3298     smd->time = scene_framenr;
3299
3300     /* do simulation */
3301
3302     // simulate the actual smoke (c++ code in intern/smoke)
3303     // DG: interesting commenting this line + deactivating loading of noise files
3304     if (framenr != startframe) {
3305       if (sds->flags & MOD_SMOKE_DISSOLVE) {
3306         /* low res dissolve */
3307         smoke_dissolve(sds->fluid, sds->diss_speed, sds->flags & MOD_SMOKE_DISSOLVE_LOG);
3308         /* high res dissolve */
3309         if (sds->wt) {
3310           smoke_dissolve_wavelet(sds->wt, sds->diss_speed, sds->flags & MOD_SMOKE_DISSOLVE_LOG);
3311         }
3312       }
3313
3314       step(depsgraph, scene, ob, smd, me, scene->r.frs_sec / scene->r.frs_sec_base);
3315     }
3316
3317     // create shadows before writing cache so they get stored
3318     smoke_calc_transparency(sds, DEG_get_evaluated_view_layer(depsgraph));
3319
3320     if (sds->wt && sds->total_cells > 1) {
3321       smoke_turbulence_step(sds->wt, sds->fluid);
3322     }
3323
3324     BKE_ptcache_validate(cache, framenr);
3325     if (framenr != startframe) {
3326       BKE_ptcache_write(&pid, framenr);
3327     }
3328
3329 #  ifdef DEBUG_TIME
3330     double end = PIL_check_seconds_timer();
3331     printf("Frame: %d, Time: %f\n\n", (int)smd->time, (float)(end - start));
3332 #  endif
3333   }
3334 }
3335
3336 struct Mesh *smokeModifier_do(
3337     SmokeModifierData *smd, Depsgraph *depsgraph, Scene *scene, Object *ob, Mesh *me)
3338 {
3339   /* lock so preview render does not read smoke data while it gets modified */
3340   if ((smd->type & MOD_SMOKE_TYPE_DOMAIN) && smd->domain) {
3341     BLI_rw_mutex_lock(smd->domain->fluid_mutex, THREAD_LOCK_WRITE);
3342   }
3343
3344   smokeModifier_process(smd, depsgraph, scene, ob, me);
3345
3346   if ((smd->type & MOD_SMOKE_TYPE_DOMAIN) && smd->domain) {
3347     BLI_rw_mutex_unlock(smd->domain->fluid_mutex);
3348   }
3349
3350   /* return generated geometry for adaptive domain */
3351   Mesh *result;
3352   if (smd->type & MOD_SMOKE_TYPE_DOMAIN && smd->domain &&
3353       smd->domain->flags & MOD_SMOKE_ADAPTIVE_DOMAIN && smd->domain->base_res[0]) {
3354     result = createDomainGeometry(smd->domain, ob);
3355   }
3356   else {
3357     result = BKE_mesh_copy_for_eval(me, false);
3358   }
3359   /* XXX This is really not a nice hack, but until root of the problem is understood,
3360    * this should be an acceptable workaround I think.
3361    * See T58492 for details on the issue. */
3362   result->texflag |= ME_AUTOSPACE;
3363   return result;
3364 }
3365
3366 static float calc_voxel_transp(
3367     float *result, float *input, int res[3], int *pixel, float *tRay, float correct)
3368 {
3369   const size_t index = smoke_get_index(pixel[0], res[0], pixel[1], res[1], pixel[2]);
3370
3371   // T_ray *= T_vox
3372   *tRay *= expf(input[index] * correct);
3373
3374   if (result[index] < 0.0f) {
3375     result[index] = *tRay;
3376   }
3377
3378   return *tRay;
3379 }
3380
3381 static void bresenham_linie_3D(int x1,
3382                                int y1,
3383                                int z1,
3384                                int x2,
3385                                int y2,
3386                                int z2,
3387                                float *tRay,
3388                                bresenham_callback cb,
3389                                float *result,
3390                                float *input,
3391                                int res[3],
3392                                float correct)
3393 {
3394   int dx, dy, dz, i, l, m, n, x_inc, y_inc, z_inc, err_1, err_2, dx2, dy2, dz2;
3395   int pixel[3];
3396
3397   pixel[0] = x1;
3398   pixel[1] = y1;
3399   pixel[2] = z1;
3400
3401   dx = x2 - x1;
3402   dy = y2 - y1;
3403   dz = z2 - z1;
3404
3405   x_inc = (dx < 0) ? -1 : 1;
3406   l = abs(dx);
3407   y_inc = (dy < 0) ? -1 : 1;
3408   m = abs(dy);
3409   z_inc = (dz < 0) ? -1 : 1;
3410   n = abs(dz);
3411   dx2 = l << 1;
3412   dy2 = m << 1;
3413   dz2 = n << 1;
3414
3415   if ((l >= m) && (l >= n)) {
3416     err_1 = dy2 - l;
3417     err_2 = dz2 - l;
3418     for (i = 0; i < l; i++) {
3419       if (cb(result, input, res, pixel, tRay, correct) <= FLT_EPSILON) {
3420         break;
3421       }
3422       if (err_1 > 0) {
3423         pixel[1] += y_inc;
3424         err_1 -= dx2;
3425       }
3426       if (err_2 > 0) {
3427         pixel[2] += z_inc;
3428         err_2 -= dx2;
3429       }
3430       err_1 += dy2;
3431       err_2 += dz2;
3432       pixel[0] += x_inc;
3433     }
3434   }
3435   else if ((m >= l) && (m >= n)) {
3436     err_1 = dx2 - m;
3437     err_2 = dz2 - m;
3438     for (i = 0; i < m; i++) {
3439       if (cb(result, input, res, pixel, tRay, correct) <= FLT_EPSILON) {
3440         break;
3441       }
3442       if (err_1 > 0) {
3443         pixel[0] += x_inc;
3444         err_1 -= dy2;
3445       }
3446       if (err_2 > 0) {
3447         pixel[2] += z_inc;
3448         err_2 -= dy2;
3449       }
3450       err_1 += dx2;
3451       err_2 += dz2;
3452       pixel[1] += y_inc;
3453     }
3454   }
3455   else {
3456     err_1 = dy2 - n;
3457     err_2 = dx2 - n;
3458     for (i = 0; i < n; i++) {
3459       if (cb(result, input, res, pixel, tRay, correct) <= FLT_EPSILON) {
3460         break;
3461       }
3462       if (err_1 > 0) {
3463         pixel[1] += y_inc;
3464         err_1 -= dz2;
3465       }
3466       if (err_2 > 0) {
3467         pixel[0] += x_inc;
3468         err_2 -= dz2;
3469       }
3470       err_1 += dy2;
3471       err_2 += dx2;
3472       pixel[2] += z_inc;
3473     }
3474   }
3475   cb(result, input, res, pixel, tRay, correct);
3476 }
3477
3478 static void smoke_calc_transparency(SmokeDomainSettings *sds, ViewLayer *view_layer)
3479 {
3480   float bv[6] = {0};
3481   float light[3];
3482   int a, z, slabsize = sds->res[0] * sds->res[1], size = sds->res[0] * sds->res[1] * sds->res[2];
3483   float *density = smoke_get_density(sds->fluid);
3484   float correct = -7.0f * sds->dx;
3485
3486   if (!get_light(view_layer, light)) {
3487     return;
3488   }
3489