Cleanup: misc spelling fixes
[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
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 TaskParallelTLS *__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
918     // but is declared as "does not move".
919
920     {
921       vert_vel = MEM_callocN(sizeof(float) * numverts * 3, "smoke_obs_velocity");
922
923       if (scs->numverts != numverts || !scs->verts_old) {
924         if (scs->verts_old) {
925           MEM_freeN(scs->verts_old);
926         }
927
928         scs->verts_old = MEM_callocN(sizeof(float) * numverts * 3, "smoke_obs_verts_old");
929         scs->numverts = numverts;
930       }
931       else {
932         has_velocity = true;
933       }
934     }
935
936     /*  Transform collider vertices to
937      *   domain grid space for fast lookups */
938     for (i = 0; i < numverts; i++) {
939       float n[3];
940       float co[3];
941
942       /* vert pos */
943       mul_m4_v3(coll_ob->obmat, mvert[i].co);
944       smoke_pos_to_cell(sds, mvert[i].co);
945
946       /* vert normal */
947       normal_short_to_float_v3(n, mvert[i].no);
948       mul_mat3_m4_v3(coll_ob->obmat, n);
949       mul_mat3_m4_v3(sds->imat, n);
950       normalize_v3(n);
951       normal_float_to_short_v3(mvert[i].no, n);
952
953       /* vert velocity */
954       add_v3fl_v3fl_v3i(co, mvert[i].co, sds->shift);
955       if (has_velocity) {
956         sub_v3_v3v3(&vert_vel[i * 3], co, &scs->verts_old[i * 3]);
957         mul_v3_fl(&vert_vel[i * 3], sds->dx / dt);
958       }
959       copy_v3_v3(&scs->verts_old[i * 3], co);
960     }
961
962     if (BKE_bvhtree_from_mesh_get(&treeData, me, BVHTREE_FROM_LOOPTRI, 4)) {
963       ObstaclesFromDMData data = {
964           .sds = sds,
965           .mvert = mvert,
966           .mloop = mloop,
967           .looptri = looptri,
968           .tree = &treeData,
969           .obstacle_map = obstacle_map,
970           .has_velocity = has_velocity,
971           .vert_vel = vert_vel,
972           .velocityX = velocityX,
973           .velocityY = velocityY,
974           .velocityZ = velocityZ,
975           .num_obstacles = num_obstacles,
976       };
977       TaskParallelSettings settings;
978       BLI_parallel_range_settings_defaults(&settings);
979       settings.scheduling_mode = TASK_SCHEDULING_DYNAMIC;
980       BLI_task_parallel_range(
981           sds->res_min[2], sds->res_max[2], &data, obstacles_from_mesh_task_cb, &settings);
982     }
983     /* free bvh tree */
984     free_bvhtree_from_mesh(&treeData);
985     BKE_id_free(NULL, me);
986
987     if (vert_vel) {
988       MEM_freeN(vert_vel);
989     }
990   }
991 }
992
993 /* Animated obstacles: dx_step = ((x_new - x_old) / totalsteps) * substep */
994 static void update_obstacles(Depsgraph *depsgraph,
995                              Object *ob,
996                              SmokeDomainSettings *sds,
997                              float dt,
998                              int UNUSED(substep),
999                              int UNUSED(totalsteps))
1000 {
1001   Object **collobjs = NULL;
1002   unsigned int numcollobj = 0;
1003
1004   unsigned int collIndex;
1005   unsigned char *obstacles = smoke_get_obstacle(sds->fluid);
1006   float *velx = NULL;
1007   float *vely = NULL;
1008   float *velz = NULL;
1009   float *velxOrig = smoke_get_velocity_x(sds->fluid);
1010   float *velyOrig = smoke_get_velocity_y(sds->fluid);
1011   float *velzOrig = smoke_get_velocity_z(sds->fluid);
1012   float *density = smoke_get_density(sds->fluid);
1013   float *fuel = smoke_get_fuel(sds->fluid);
1014   float *flame = smoke_get_flame(sds->fluid);
1015   float *r = smoke_get_color_r(sds->fluid);
1016   float *g = smoke_get_color_g(sds->fluid);
1017   float *b = smoke_get_color_b(sds->fluid);
1018   unsigned int z;
1019
1020   int *num_obstacles = MEM_callocN(sizeof(int) * sds->res[0] * sds->res[1] * sds->res[2],
1021                                    "smoke_num_obstacles");
1022
1023   smoke_get_ob_velocity(sds->fluid, &velx, &vely, &velz);
1024
1025   // TODO: delete old obstacle flags
1026   for (z = 0; z < sds->res[0] * sds->res[1] * sds->res[2]; z++) {
1027     if (obstacles[z] & 8)  // Do not delete static obstacles
1028     {
1029       obstacles[z] = 0;
1030     }
1031
1032     velx[z] = 0;
1033     vely[z] = 0;
1034     velz[z] = 0;
1035   }
1036
1037   collobjs = BKE_collision_objects_create(
1038       depsgraph, ob, sds->coll_group, &numcollobj, eModifierType_Smoke);
1039
1040   // update obstacle tags in cells
1041   for (collIndex = 0; collIndex < numcollobj; collIndex++) {
1042     Object *collob = collobjs[collIndex];
1043     SmokeModifierData *smd2 = (SmokeModifierData *)modifiers_findByType(collob,
1044                                                                         eModifierType_Smoke);
1045
1046     // DG TODO: check if modifier is active?
1047
1048     if ((smd2->type & MOD_SMOKE_TYPE_COLL) && smd2->coll) {
1049       SmokeCollSettings *scs = smd2->coll;
1050       obstacles_from_mesh(collob, sds, scs, obstacles, velx, vely, velz, num_obstacles, dt);
1051     }
1052   }
1053
1054   BKE_collision_objects_free(collobjs);
1055
1056   /* obstacle cells should not contain any velocity from the smoke simulation */
1057   for (z = 0; z < sds->res[0] * sds->res[1] * sds->res[2]; z++) {
1058     if (obstacles[z]) {
1059       velxOrig[z] = 0;
1060       velyOrig[z] = 0;
1061       velzOrig[z] = 0;
1062       density[z] = 0;
1063       if (fuel) {
1064         fuel[z] = 0;
1065         flame[z] = 0;
1066       }
1067       if (r) {
1068         r[z] = 0;
1069         g[z] = 0;
1070         b[z] = 0;
1071       }
1072     }
1073     /* average velocities from multiple obstacles in one cell */
1074     if (num_obstacles[z]) {
1075       velx[z] /= num_obstacles[z];
1076       vely[z] /= num_obstacles[z];
1077       velz[z] /= num_obstacles[z];
1078     }
1079   }
1080
1081   MEM_freeN(num_obstacles);
1082 }
1083
1084 /**********************************************************
1085  * Flow emission code
1086  **********************************************************/
1087
1088 typedef struct EmissionMap {
1089   float *influence;
1090   float *influence_high;
1091   float *velocity;
1092   int min[3], max[3], res[3];
1093   int hmin[3], hmax[3], hres[3];
1094   int total_cells, valid;
1095 } EmissionMap;
1096
1097 static void em_boundInsert(EmissionMap *em, float point[3])
1098 {
1099   int i = 0;
1100   if (!em->valid) {
1101     for (; i < 3; i++) {
1102       em->min[i] = (int)floor(point[i]);
1103       em->max[i] = (int)ceil(point[i]);
1104     }
1105     em->valid = 1;
1106   }
1107   else {
1108     for (; i < 3; i++) {
1109       if (point[i] < em->min[i]) {
1110         em->min[i] = (int)floor(point[i]);
1111       }
1112       if (point[i] > em->max[i]) {
1113         em->max[i] = (int)ceil(point[i]);
1114       }
1115     }
1116   }
1117 }
1118
1119 static void clampBoundsInDomain(SmokeDomainSettings *sds,
1120                                 int min[3],
1121                                 int max[3],
1122                                 float *min_vel,
1123                                 float *max_vel,
1124                                 int margin,
1125                                 float dt)
1126 {
1127   int i;
1128   for (i = 0; i < 3; i++) {
1129     int adapt = (sds->flags & MOD_SMOKE_ADAPTIVE_DOMAIN) ? sds->adapt_res : 0;
1130     /* add margin */
1131     min[i] -= margin;
1132     max[i] += margin;
1133
1134     /* adapt to velocity */
1135     if (min_vel && min_vel[i] < 0.0f) {
1136       min[i] += (int)floor(min_vel[i] * dt);
1137     }
1138     if (max_vel && max_vel[i] > 0.0f) {
1139       max[i] += (int)ceil(max_vel[i] * dt);
1140     }
1141
1142     /* clamp within domain max size */
1143     CLAMP(min[i], -adapt, sds->base_res[i] + adapt);
1144     CLAMP(max[i], -adapt, sds->base_res[i] + adapt);
1145   }
1146 }
1147
1148 static void em_allocateData(EmissionMap *em, bool use_velocity, int hires_mul)
1149 {
1150   int i, res[3];
1151
1152   for (i = 0; i < 3; i++) {
1153     res[i] = em->max[i] - em->min[i];
1154     if (res[i] <= 0) {
1155       return;
1156     }
1157   }
1158   em->total_cells = res[0] * res[1] * res[2];
1159   copy_v3_v3_int(em->res, res);
1160
1161   em->influence = MEM_callocN(sizeof(float) * em->total_cells, "smoke_flow_influence");
1162   if (use_velocity) {
1163     em->velocity = MEM_callocN(sizeof(float) * em->total_cells * 3, "smoke_flow_velocity");
1164   }
1165
1166   /* allocate high resolution map if required */
1167   if (hires_mul > 1) {
1168     int total_cells_high = em->total_cells * (hires_mul * hires_mul * hires_mul);
1169
1170     for (i = 0; i < 3; i++) {
1171       em->hmin[i] = em->min[i] * hires_mul;
1172       em->hmax[i] = em->max[i] * hires_mul;
1173       em->hres[i] = em->res[i] * hires_mul;
1174     }
1175
1176     em->influence_high = MEM_callocN(sizeof(float) * total_cells_high,
1177                                      "smoke_flow_influence_high");
1178   }
1179   em->valid = 1;
1180 }
1181
1182 static void em_freeData(EmissionMap *em)
1183 {
1184   if (em->influence) {
1185     MEM_freeN(em->influence);
1186   }
1187   if (em->influence_high) {
1188     MEM_freeN(em->influence_high);
1189   }
1190   if (em->velocity) {
1191     MEM_freeN(em->velocity);
1192   }
1193 }
1194
1195 static void em_combineMaps(
1196     EmissionMap *output, EmissionMap *em2, int hires_multiplier, int additive, float sample_size)
1197 {
1198   int i, x, y, z;
1199
1200   /* copyfill input 1 struct and clear output for new allocation */
1201   EmissionMap em1;
1202   memcpy(&em1, output, sizeof(EmissionMap));
1203   memset(output, 0, sizeof(EmissionMap));
1204
1205   for (i = 0; i < 3; i++) {
1206     if (em1.valid) {
1207       output->min[i] = MIN2(em1.min[i], em2->min[i]);
1208       output->max[i] = MAX2(em1.max[i], em2->max[i]);
1209     }
1210     else {
1211       output->min[i] = em2->min[i];
1212       output->max[i] = em2->max[i];
1213     }
1214   }
1215   /* allocate output map */
1216   em_allocateData(output, (em1.velocity || em2->velocity), hires_multiplier);
1217
1218   /* base resolution inputs */
1219   for (x = output->min[0]; x < output->max[0]; x++) {
1220     for (y = output->min[1]; y < output->max[1]; y++) {
1221       for (z = output->min[2]; z < output->max[2]; z++) {
1222         int index_out = smoke_get_index(x - output->min[0],
1223                                         output->res[0],
1224                                         y - output->min[1],
1225                                         output->res[1],
1226                                         z - output->min[2]);
1227
1228         /* initialize with first input if in range */
1229         if (x >= em1.min[0] && x < em1.max[0] && y >= em1.min[1] && y < em1.max[1] &&
1230             z >= em1.min[2] && z < em1.max[2]) {
1231           int index_in = smoke_get_index(
1232               x - em1.min[0], em1.res[0], y - em1.min[1], em1.res[1], z - em1.min[2]);
1233
1234           /* values */
1235           output->influence[index_out] = em1.influence[index_in];
1236           if (output->velocity && em1.velocity) {
1237             copy_v3_v3(&output->velocity[index_out * 3], &em1.velocity[index_in * 3]);
1238           }
1239         }
1240
1241         /* apply second input if in range */
1242         if (x >= em2->min[0] && x < em2->max[0] && y >= em2->min[1] && y < em2->max[1] &&
1243             z >= em2->min[2] && z < em2->max[2]) {
1244           int index_in = smoke_get_index(
1245               x - em2->min[0], em2->res[0], y - em2->min[1], em2->res[1], z - em2->min[2]);
1246
1247           /* values */
1248           if (additive) {
1249             output->influence[index_out] += em2->influence[index_in] * sample_size;
1250           }
1251           else {
1252             output->influence[index_out] = MAX2(em2->influence[index_in],
1253                                                 output->influence[index_out]);
1254           }
1255           if (output->velocity && em2->velocity) {
1256             /* last sample replaces the velocity */
1257             output->velocity[index_out * 3] = ADD_IF_LOWER(output->velocity[index_out * 3],
1258                                                            em2->velocity[index_in * 3]);
1259             output->velocity[index_out * 3 + 1] = ADD_IF_LOWER(output->velocity[index_out * 3 + 1],
1260                                                                em2->velocity[index_in * 3 + 1]);
1261             output->velocity[index_out * 3 + 2] = ADD_IF_LOWER(output->velocity[index_out * 3 + 2],
1262                                                                em2->velocity[index_in * 3 + 2]);
1263           }
1264         }
1265       }  // low res loop
1266     }
1267   }
1268
1269   /* initialize high resolution input if available */
1270   if (output->influence_high) {
1271     for (x = output->hmin[0]; x < output->hmax[0]; x++) {
1272       for (y = output->hmin[1]; y < output->hmax[1]; y++) {
1273         for (z = output->hmin[2]; z < output->hmax[2]; z++) {
1274           int index_out = smoke_get_index(x - output->hmin[0],
1275                                           output->hres[0],
1276                                           y - output->hmin[1],
1277                                           output->hres[1],
1278                                           z - output->hmin[2]);
1279
1280           /* initialize with first input if in range */
1281           if (x >= em1.hmin[0] && x < em1.hmax[0] && y >= em1.hmin[1] && y < em1.hmax[1] &&
1282               z >= em1.hmin[2] && z < em1.hmax[2]) {
1283             int index_in = smoke_get_index(
1284                 x - em1.hmin[0], em1.hres[0], y - em1.hmin[1], em1.hres[1], z - em1.hmin[2]);
1285             /* values */
1286             output->influence_high[index_out] = em1.influence_high[index_in];
1287           }
1288
1289           /* apply second input if in range */
1290           if (x >= em2->hmin[0] && x < em2->hmax[0] && y >= em2->hmin[1] && y < em2->hmax[1] &&
1291               z >= em2->hmin[2] && z < em2->hmax[2]) {
1292             int index_in = smoke_get_index(
1293                 x - em2->hmin[0], em2->hres[0], y - em2->hmin[1], em2->hres[1], z - em2->hmin[2]);
1294
1295             /* values */
1296             if (additive) {
1297               output->influence_high[index_out] += em2->influence_high[index_in] * sample_size;
1298             }
1299             else {
1300               output->influence_high[index_out] = MAX2(em2->influence_high[index_in],
1301                                                        output->influence_high[index_out]);
1302             }
1303           }
1304         }  // high res loop
1305       }
1306     }
1307   }
1308
1309   /* free original data */
1310   em_freeData(&em1);
1311 }
1312
1313 typedef struct EmitFromParticlesData {
1314   SmokeFlowSettings *sfs;
1315   KDTree_3d *tree;
1316   int hires_multiplier;
1317
1318   EmissionMap *em;
1319   float *particle_vel;
1320   float hr;
1321
1322   int *min, *max, *res;
1323
1324   float solid;
1325   float smooth;
1326   float hr_smooth;
1327 } EmitFromParticlesData;
1328
1329 static void emit_from_particles_task_cb(void *__restrict userdata,
1330                                         const int z,
1331                                         const TaskParallelTLS *__restrict UNUSED(tls))
1332 {
1333   EmitFromParticlesData *data = userdata;
1334   SmokeFlowSettings *sfs = data->sfs;
1335   EmissionMap *em = data->em;
1336   const int hires_multiplier = data->hires_multiplier;
1337
1338   for (int x = data->min[0]; x < data->max[0]; x++) {
1339     for (int y = data->min[1]; y < data->max[1]; y++) {
1340       /* take low res samples where possible */
1341       if (hires_multiplier <= 1 ||
1342           !(x % hires_multiplier || y % hires_multiplier || z % hires_multiplier)) {
1343         /* get low res space coordinates */
1344         const int lx = x / hires_multiplier;
1345         const int ly = y / hires_multiplier;
1346         const int lz = z / hires_multiplier;
1347
1348         const int index = smoke_get_index(
1349             lx - em->min[0], em->res[0], ly - em->min[1], em->res[1], lz - em->min[2]);
1350         const float ray_start[3] = {((float)lx) + 0.5f, ((float)ly) + 0.5f, ((float)lz) + 0.5f};
1351
1352         /* find particle distance from the kdtree */
1353         KDTreeNearest_3d nearest;
1354         const float range = data->solid + data->smooth;
1355         BLI_kdtree_3d_find_nearest(data->tree, ray_start, &nearest);
1356
1357         if (nearest.dist < range) {
1358           em->influence[index] = (nearest.dist < data->solid) ?
1359                                      1.0f :
1360                                      (1.0f - (nearest.dist - data->solid) / data->smooth);
1361           /* Uses particle velocity as initial velocity for smoke */
1362           if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY &&
1363               (sfs->psys->part->phystype != PART_PHYS_NO)) {
1364             madd_v3_v3fl(
1365                 &em->velocity[index * 3], &data->particle_vel[nearest.index * 3], sfs->vel_multi);
1366           }
1367         }
1368       }
1369
1370       /* take high res samples if required */
1371       if (hires_multiplier > 1) {
1372         /* get low res space coordinates */
1373         const float lx = ((float)x) * data->hr;
1374         const float ly = ((float)y) * data->hr;
1375         const float lz = ((float)z) * data->hr;
1376
1377         const int index = smoke_get_index(
1378             x - data->min[0], data->res[0], y - data->min[1], data->res[1], z - data->min[2]);
1379         const float ray_start[3] = {
1380             lx + 0.5f * data->hr, ly + 0.5f * data->hr, lz + 0.5f * data->hr};
1381
1382         /* find particle distance from the kdtree */
1383         KDTreeNearest_3d nearest;
1384         const float range = data->solid + data->hr_smooth;
1385         BLI_kdtree_3d_find_nearest(data->tree, ray_start, &nearest);
1386
1387         if (nearest.dist < range) {
1388           em->influence_high[index] = (nearest.dist < data->solid) ?
1389                                           1.0f :
1390                                           (1.0f - (nearest.dist - data->solid) / data->smooth);
1391         }
1392       }
1393     }
1394   }
1395 }
1396
1397 static void emit_from_particles(Object *flow_ob,
1398                                 SmokeDomainSettings *sds,
1399                                 SmokeFlowSettings *sfs,
1400                                 EmissionMap *em,
1401                                 Depsgraph *depsgraph,
1402                                 Scene *scene,
1403                                 float dt)
1404 {
1405   /* Is particle system selected. */
1406   if (sfs && sfs->psys && sfs->psys->part &&
1407       ELEM(sfs->psys->part->type, PART_EMITTER, PART_FLUID)) {
1408     ParticleSimulationData sim;
1409     ParticleSystem *psys = sfs->psys;
1410     float *particle_pos;
1411     float *particle_vel;
1412     int totpart = psys->totpart, totchild;
1413     int p = 0;
1414     int valid_particles = 0;
1415     int bounds_margin = 1;
1416
1417     /* radius based flow */
1418     const float solid = sfs->particle_size * 0.5f;
1419     const float smooth = 0.5f; /* add 0.5 cells of linear falloff to reduce aliasing */
1420     int hires_multiplier = 1;
1421     KDTree_3d *tree = NULL;
1422
1423     sim.depsgraph = depsgraph;
1424     sim.scene = scene;
1425     sim.ob = flow_ob;
1426     sim.psys = psys;
1427     sim.psys->lattice_deform_data = psys_create_lattice_deform_data(&sim);
1428
1429     /* prepare curvemapping tables */
1430     if ((psys->part->child_flag & PART_CHILD_USE_CLUMP_CURVE) && psys->part->clumpcurve) {
1431       curvemapping_changed_all(psys->part->clumpcurve);
1432     }
1433     if ((psys->part->child_flag & PART_CHILD_USE_ROUGH_CURVE) && psys->part->roughcurve) {
1434       curvemapping_changed_all(psys->part->roughcurve);
1435     }
1436     if ((psys->part->child_flag & PART_CHILD_USE_TWIST_CURVE) && psys->part->twistcurve) {
1437       curvemapping_changed_all(psys->part->twistcurve);
1438     }
1439
1440     /* initialize particle cache */
1441     if (psys->part->type == PART_HAIR) {
1442       // TODO: PART_HAIR not supported whatsoever
1443       totchild = 0;
1444     }
1445     else {
1446       totchild = psys->totchild * psys->part->disp / 100;
1447     }
1448
1449     particle_pos = MEM_callocN(sizeof(float) * (totpart + totchild) * 3, "smoke_flow_particles");
1450     particle_vel = MEM_callocN(sizeof(float) * (totpart + totchild) * 3, "smoke_flow_particles");
1451
1452     /* setup particle radius emission if enabled */
1453     if (sfs->flags & MOD_SMOKE_FLOW_USE_PART_SIZE) {
1454       tree = BLI_kdtree_3d_new(psys->totpart + psys->totchild);
1455
1456       /* check need for high resolution map */
1457       if ((sds->flags & MOD_SMOKE_HIGHRES) && (sds->highres_sampling == SM_HRES_FULLSAMPLE)) {
1458         hires_multiplier = sds->amplify + 1;
1459       }
1460
1461       bounds_margin = (int)ceil(solid + smooth);
1462     }
1463
1464     /* calculate local position for each particle */
1465     for (p = 0; p < totpart + totchild; p++) {
1466       ParticleKey state;
1467       float *pos;
1468       if (p < totpart) {
1469         if (psys->particles[p].flag & (PARS_NO_DISP | PARS_UNEXIST)) {
1470           continue;
1471         }
1472       }
1473       else {
1474         /* handle child particle */
1475         ChildParticle *cpa = &psys->child[p - totpart];
1476         if (psys->particles[cpa->parent].flag & (PARS_NO_DISP | PARS_UNEXIST)) {
1477           continue;
1478         }
1479       }
1480
1481       state.time = DEG_get_ctime(depsgraph); /* use depsgraph time */
1482       if (psys_get_particle_state(&sim, p, &state, 0) == 0) {
1483         continue;
1484       }
1485
1486       /* location */
1487       pos = &particle_pos[valid_particles * 3];
1488       copy_v3_v3(pos, state.co);
1489       smoke_pos_to_cell(sds, pos);
1490
1491       /* velocity */
1492       copy_v3_v3(&particle_vel[valid_particles * 3], state.vel);
1493       mul_mat3_m4_v3(sds->imat, &particle_vel[valid_particles * 3]);
1494
1495       if (sfs->flags & MOD_SMOKE_FLOW_USE_PART_SIZE) {
1496         BLI_kdtree_3d_insert(tree, valid_particles, pos);
1497       }
1498
1499       /* calculate emission map bounds */
1500       em_boundInsert(em, pos);
1501       valid_particles++;
1502     }
1503
1504     /* set emission map */
1505     clampBoundsInDomain(sds, em->min, em->max, NULL, NULL, bounds_margin, dt);
1506     em_allocateData(em, sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY, hires_multiplier);
1507
1508     if (!(sfs->flags & MOD_SMOKE_FLOW_USE_PART_SIZE)) {
1509       for (p = 0; p < valid_particles; p++) {
1510         int cell[3];
1511         size_t i = 0;
1512         size_t index = 0;
1513         int badcell = 0;
1514
1515         /* 1. get corresponding cell */
1516         cell[0] = floor(particle_pos[p * 3]) - em->min[0];
1517         cell[1] = floor(particle_pos[p * 3 + 1]) - em->min[1];
1518         cell[2] = floor(particle_pos[p * 3 + 2]) - em->min[2];
1519         /* check if cell is valid (in the domain boundary) */
1520         for (i = 0; i < 3; i++) {
1521           if ((cell[i] > em->res[i] - 1) || (cell[i] < 0)) {
1522             badcell = 1;
1523             break;
1524           }
1525         }
1526         if (badcell) {
1527           continue;
1528         }
1529         /* get cell index */
1530         index = smoke_get_index(cell[0], em->res[0], cell[1], em->res[1], cell[2]);
1531         /* Add influence to emission map */
1532         em->influence[index] = 1.0f;
1533         /* Uses particle velocity as initial velocity for smoke */
1534         if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY && (psys->part->phystype != PART_PHYS_NO)) {
1535           madd_v3_v3fl(&em->velocity[index * 3], &particle_vel[p * 3], sfs->vel_multi);
1536         }
1537       }  // particles loop
1538     }
1539     else if (valid_particles > 0) {  // MOD_SMOKE_FLOW_USE_PART_SIZE
1540       int min[3], max[3], res[3];
1541       const float hr = 1.0f / ((float)hires_multiplier);
1542       /* slightly adjust high res antialias smoothness based on number of divisions
1543        * to allow smaller details but yet not differing too much from the low res size */
1544       const float hr_smooth = smooth * powf(hr, 1.0f / 3.0f);
1545
1546       /* setup loop bounds */
1547       for (int i = 0; i < 3; i++) {
1548         min[i] = em->min[i] * hires_multiplier;
1549         max[i] = em->max[i] * hires_multiplier;
1550         res[i] = em->res[i] * hires_multiplier;
1551       }
1552
1553       BLI_kdtree_3d_balance(tree);
1554
1555       EmitFromParticlesData data = {
1556           .sfs = sfs,
1557           .tree = tree,
1558           .hires_multiplier = hires_multiplier,
1559           .hr = hr,
1560           .em = em,
1561           .particle_vel = particle_vel,
1562           .min = min,
1563           .max = max,
1564           .res = res,
1565           .solid = solid,
1566           .smooth = smooth,
1567           .hr_smooth = hr_smooth,
1568       };
1569
1570       TaskParallelSettings settings;
1571       BLI_parallel_range_settings_defaults(&settings);
1572       settings.scheduling_mode = TASK_SCHEDULING_DYNAMIC;
1573       BLI_task_parallel_range(min[2], max[2], &data, emit_from_particles_task_cb, &settings);
1574     }
1575
1576     if (sfs->flags & MOD_SMOKE_FLOW_USE_PART_SIZE) {
1577       BLI_kdtree_3d_free(tree);
1578     }
1579
1580     /* free data */
1581     if (particle_pos) {
1582       MEM_freeN(particle_pos);
1583     }
1584     if (particle_vel) {
1585       MEM_freeN(particle_vel);
1586     }
1587   }
1588 }
1589
1590 static void sample_mesh(SmokeFlowSettings *sfs,
1591                         const MVert *mvert,
1592                         const MLoop *mloop,
1593                         const MLoopTri *mlooptri,
1594                         const MLoopUV *mloopuv,
1595                         float *influence_map,
1596                         float *velocity_map,
1597                         int index,
1598                         const int base_res[3],
1599                         float flow_center[3],
1600                         BVHTreeFromMesh *treeData,
1601                         const float ray_start[3],
1602                         const float *vert_vel,
1603                         bool has_velocity,
1604                         int defgrp_index,
1605                         MDeformVert *dvert,
1606                         float x,
1607                         float y,
1608                         float z)
1609 {
1610   float ray_dir[3] = {1.0f, 0.0f, 0.0f};
1611   BVHTreeRayHit hit = {0};
1612   BVHTreeNearest nearest = {0};
1613
1614   float volume_factor = 0.0f;
1615   float sample_str = 0.0f;
1616
1617   hit.index = -1;
1618   hit.dist = 9999;
1619   nearest.index = -1;
1620   nearest.dist_sq = sfs->surface_distance *
1621                     sfs->surface_distance; /* find_nearest uses squared distance */
1622
1623   /* Check volume collision */
1624   if (sfs->volume_density) {
1625     if (BLI_bvhtree_ray_cast(treeData->tree,
1626                              ray_start,
1627                              ray_dir,
1628                              0.0f,
1629                              &hit,
1630                              treeData->raycast_callback,
1631                              treeData) != -1) {
1632       float dot = ray_dir[0] * hit.no[0] + ray_dir[1] * hit.no[1] + ray_dir[2] * hit.no[2];
1633       /* If ray and hit face normal are facing same direction
1634        * hit point is inside a closed mesh. */
1635       if (dot >= 0) {
1636         /* Also cast a ray in opposite direction to make sure
1637          * point is at least surrounded by two faces */
1638         negate_v3(ray_dir);
1639         hit.index = -1;
1640         hit.dist = 9999;
1641
1642         BLI_bvhtree_ray_cast(
1643             treeData->tree, ray_start, ray_dir, 0.0f, &hit, treeData->raycast_callback, treeData);
1644         if (hit.index != -1) {
1645           volume_factor = sfs->volume_density;
1646         }
1647       }
1648     }
1649   }
1650
1651   /* find the nearest point on the mesh */
1652   if (BLI_bvhtree_find_nearest(
1653           treeData->tree, ray_start, &nearest, treeData->nearest_callback, treeData) != -1) {
1654     float weights[3];
1655     int v1, v2, v3, f_index = nearest.index;
1656     float n1[3], n2[3], n3[3], hit_normal[3];
1657
1658     /* emit from surface based on distance */
1659     if (sfs->surface_distance) {
1660       sample_str = sqrtf(nearest.dist_sq) / sfs->surface_distance;
1661       CLAMP(sample_str, 0.0f, 1.0f);
1662       sample_str = pow(1.0f - sample_str, 0.5f);
1663     }
1664     else {
1665       sample_str = 0.0f;
1666     }
1667
1668     /* calculate barycentric weights for nearest point */
1669     v1 = mloop[mlooptri[f_index].tri[0]].v;
1670     v2 = mloop[mlooptri[f_index].tri[1]].v;
1671     v3 = mloop[mlooptri[f_index].tri[2]].v;
1672     interp_weights_tri_v3(weights, mvert[v1].co, mvert[v2].co, mvert[v3].co, nearest.co);
1673
1674     if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY && velocity_map) {
1675       /* apply normal directional velocity */
1676       if (sfs->vel_normal) {
1677         /* interpolate vertex normal vectors to get nearest point normal */
1678         normal_short_to_float_v3(n1, mvert[v1].no);
1679         normal_short_to_float_v3(n2, mvert[v2].no);
1680         normal_short_to_float_v3(n3, mvert[v3].no);
1681         interp_v3_v3v3v3(hit_normal, n1, n2, n3, weights);
1682         normalize_v3(hit_normal);
1683         /* apply normal directional and random velocity
1684          * - TODO: random disabled for now since it doesn't really work well
1685          *   as pressure calc smoothens it out. */
1686         velocity_map[index * 3] += hit_normal[0] * sfs->vel_normal * 0.25f;
1687         velocity_map[index * 3 + 1] += hit_normal[1] * sfs->vel_normal * 0.25f;
1688         velocity_map[index * 3 + 2] += hit_normal[2] * sfs->vel_normal * 0.25f;
1689         /* TODO: for fire emitted from mesh surface we can use
1690          * Vf = Vs + (Ps/Pf - 1)*S to model gaseous expansion from solid to fuel */
1691       }
1692       /* apply object velocity */
1693       if (has_velocity && sfs->vel_multi) {
1694         float hit_vel[3];
1695         interp_v3_v3v3v3(
1696             hit_vel, &vert_vel[v1 * 3], &vert_vel[v2 * 3], &vert_vel[v3 * 3], weights);
1697         velocity_map[index * 3] += hit_vel[0] * sfs->vel_multi;
1698         velocity_map[index * 3 + 1] += hit_vel[1] * sfs->vel_multi;
1699         velocity_map[index * 3 + 2] += hit_vel[2] * sfs->vel_multi;
1700       }
1701     }
1702
1703     /* apply vertex group influence if used */
1704     if (defgrp_index != -1 && dvert) {
1705       float weight_mask = defvert_find_weight(&dvert[v1], defgrp_index) * weights[0] +
1706                           defvert_find_weight(&dvert[v2], defgrp_index) * weights[1] +
1707                           defvert_find_weight(&dvert[v3], defgrp_index) * weights[2];
1708       sample_str *= weight_mask;
1709     }
1710
1711     /* apply emission texture */
1712     if ((sfs->flags & MOD_SMOKE_FLOW_TEXTUREEMIT) && sfs->noise_texture) {
1713       float tex_co[3] = {0};
1714       TexResult texres;
1715
1716       if (sfs->texture_type == MOD_SMOKE_FLOW_TEXTURE_MAP_AUTO) {
1717         tex_co[0] = ((x - flow_center[0]) / base_res[0]) / sfs->texture_size;
1718         tex_co[1] = ((y - flow_center[1]) / base_res[1]) / sfs->texture_size;
1719         tex_co[2] = ((z - flow_center[2]) / base_res[2] - sfs->texture_offset) / sfs->texture_size;
1720       }
1721       else if (mloopuv) {
1722         const float *uv[3];
1723         uv[0] = mloopuv[mlooptri[f_index].tri[0]].uv;
1724         uv[1] = mloopuv[mlooptri[f_index].tri[1]].uv;
1725         uv[2] = mloopuv[mlooptri[f_index].tri[2]].uv;
1726
1727         interp_v2_v2v2v2(tex_co, UNPACK3(uv), weights);
1728
1729         /* map between -1.0f and 1.0f */
1730         tex_co[0] = tex_co[0] * 2.0f - 1.0f;
1731         tex_co[1] = tex_co[1] * 2.0f - 1.0f;
1732         tex_co[2] = sfs->texture_offset;
1733       }
1734       texres.nor = NULL;
1735       BKE_texture_get_value(NULL, sfs->noise_texture, tex_co, &texres, false);
1736       sample_str *= texres.tin;
1737     }
1738   }
1739
1740   /* multiply initial velocity by emitter influence */
1741   if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY && velocity_map) {
1742     mul_v3_fl(&velocity_map[index * 3], sample_str);
1743   }
1744
1745   /* apply final influence based on volume factor */
1746   influence_map[index] = MAX2(volume_factor, sample_str);
1747 }
1748
1749 typedef struct EmitFromDMData {
1750   SmokeDomainSettings *sds;
1751   SmokeFlowSettings *sfs;
1752   const MVert *mvert;
1753   const MLoop *mloop;
1754   const MLoopTri *mlooptri;
1755   const MLoopUV *mloopuv;
1756   MDeformVert *dvert;
1757   int defgrp_index;
1758
1759   BVHTreeFromMesh *tree;
1760   int hires_multiplier;
1761   float hr;
1762
1763   EmissionMap *em;
1764   bool has_velocity;
1765   float *vert_vel;
1766
1767   float *flow_center;
1768   int *min, *max, *res;
1769 } EmitFromDMData;
1770
1771 static void emit_from_mesh_task_cb(void *__restrict userdata,
1772                                    const int z,
1773                                    const TaskParallelTLS *__restrict UNUSED(tls))
1774 {
1775   EmitFromDMData *data = userdata;
1776   EmissionMap *em = data->em;
1777   const int hires_multiplier = data->hires_multiplier;
1778
1779   for (int x = data->min[0]; x < data->max[0]; x++) {
1780     for (int y = data->min[1]; y < data->max[1]; y++) {
1781       /* take low res samples where possible */
1782       if (hires_multiplier <= 1 ||
1783           !(x % hires_multiplier || y % hires_multiplier || z % hires_multiplier)) {
1784         /* get low res space coordinates */
1785         const int lx = x / hires_multiplier;
1786         const int ly = y / hires_multiplier;
1787         const int lz = z / hires_multiplier;
1788
1789         const int index = smoke_get_index(
1790             lx - em->min[0], em->res[0], ly - em->min[1], em->res[1], lz - em->min[2]);
1791         const float ray_start[3] = {((float)lx) + 0.5f, ((float)ly) + 0.5f, ((float)lz) + 0.5f};
1792
1793         sample_mesh(data->sfs,
1794                     data->mvert,
1795                     data->mloop,
1796                     data->mlooptri,
1797                     data->mloopuv,
1798                     em->influence,
1799                     em->velocity,
1800                     index,
1801                     data->sds->base_res,
1802                     data->flow_center,
1803                     data->tree,
1804                     ray_start,
1805                     data->vert_vel,
1806                     data->has_velocity,
1807                     data->defgrp_index,
1808                     data->dvert,
1809                     (float)lx,
1810                     (float)ly,
1811                     (float)lz);
1812       }
1813
1814       /* take high res samples if required */
1815       if (hires_multiplier > 1) {
1816         /* get low res space coordinates */
1817         const float lx = ((float)x) * data->hr;
1818         const float ly = ((float)y) * data->hr;
1819         const float lz = ((float)z) * data->hr;
1820
1821         const int index = smoke_get_index(
1822             x - data->min[0], data->res[0], y - data->min[1], data->res[1], z - data->min[2]);
1823         const float ray_start[3] = {
1824             lx + 0.5f * data->hr,
1825             ly + 0.5f * data->hr,
1826             lz + 0.5f * data->hr,
1827         };
1828
1829         sample_mesh(data->sfs,
1830                     data->mvert,
1831                     data->mloop,
1832                     data->mlooptri,
1833                     data->mloopuv,
1834                     em->influence_high,
1835                     NULL,
1836                     index,
1837                     data->sds->base_res,
1838                     data->flow_center,
1839                     data->tree,
1840                     ray_start,
1841                     data->vert_vel,
1842                     data->has_velocity,
1843                     data->defgrp_index,
1844                     data->dvert,
1845                     /* x,y,z needs to be always lowres */
1846                     lx,
1847                     ly,
1848                     lz);
1849       }
1850     }
1851   }
1852 }
1853
1854 static void emit_from_mesh(
1855     Object *flow_ob, SmokeDomainSettings *sds, SmokeFlowSettings *sfs, EmissionMap *em, float dt)
1856 {
1857   if (sfs->mesh) {
1858     Mesh *me;
1859     int defgrp_index = sfs->vgroup_density - 1;
1860     MDeformVert *dvert = NULL;
1861     MVert *mvert = NULL;
1862     const MLoopTri *mlooptri = NULL;
1863     const MLoopUV *mloopuv = NULL;
1864     const MLoop *mloop = NULL;
1865     BVHTreeFromMesh treeData = {NULL};
1866     int numOfVerts, i;
1867     float flow_center[3] = {0};
1868
1869     float *vert_vel = NULL;
1870     int has_velocity = 0;
1871     int min[3], max[3], res[3];
1872     int hires_multiplier = 1;
1873
1874     /* copy mesh for thread safety because we modify it,
1875      * main issue is its VertArray being modified, then replaced and freed
1876      */
1877     me = BKE_mesh_copy_for_eval(sfs->mesh, true);
1878
1879     /* Duplicate vertices to modify. */
1880     if (me->mvert) {
1881       me->mvert = MEM_dupallocN(me->mvert);
1882       CustomData_set_layer(&me->vdata, CD_MVERT, me->mvert);
1883     }
1884
1885     BKE_mesh_ensure_normals(me);
1886     mvert = me->mvert;
1887     numOfVerts = me->totvert;
1888     dvert = CustomData_get_layer(&me->vdata, CD_MDEFORMVERT);
1889     mloopuv = CustomData_get_layer_named(&me->ldata, CD_MLOOPUV, sfs->uvlayer_name);
1890     mloop = me->mloop;
1891     mlooptri = BKE_mesh_runtime_looptri_ensure(me);
1892
1893     if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY) {
1894       vert_vel = MEM_callocN(sizeof(float) * numOfVerts * 3, "smoke_flow_velocity");
1895
1896       if (sfs->numverts != numOfVerts || !sfs->verts_old) {
1897         if (sfs->verts_old) {
1898           MEM_freeN(sfs->verts_old);
1899         }
1900         sfs->verts_old = MEM_callocN(sizeof(float) * numOfVerts * 3, "smoke_flow_verts_old");
1901         sfs->numverts = numOfVerts;
1902       }
1903       else {
1904         has_velocity = 1;
1905       }
1906     }
1907
1908     /*  Transform mesh vertices to
1909      *   domain grid space for fast lookups */
1910     for (i = 0; i < numOfVerts; i++) {
1911       float n[3];
1912       /* vert pos */
1913       mul_m4_v3(flow_ob->obmat, mvert[i].co);
1914       smoke_pos_to_cell(sds, mvert[i].co);
1915       /* vert normal */
1916       normal_short_to_float_v3(n, mvert[i].no);
1917       mul_mat3_m4_v3(flow_ob->obmat, n);
1918       mul_mat3_m4_v3(sds->imat, n);
1919       normalize_v3(n);
1920       normal_float_to_short_v3(mvert[i].no, n);
1921       /* vert velocity */
1922       if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY) {
1923         float co[3];
1924         add_v3fl_v3fl_v3i(co, mvert[i].co, sds->shift);
1925         if (has_velocity) {
1926           sub_v3_v3v3(&vert_vel[i * 3], co, &sfs->verts_old[i * 3]);
1927           mul_v3_fl(&vert_vel[i * 3], sds->dx / dt);
1928         }
1929         copy_v3_v3(&sfs->verts_old[i * 3], co);
1930       }
1931
1932       /* calculate emission map bounds */
1933       em_boundInsert(em, mvert[i].co);
1934     }
1935     mul_m4_v3(flow_ob->obmat, flow_center);
1936     smoke_pos_to_cell(sds, flow_center);
1937
1938     /* check need for high resolution map */
1939     if ((sds->flags & MOD_SMOKE_HIGHRES) && (sds->highres_sampling == SM_HRES_FULLSAMPLE)) {
1940       hires_multiplier = sds->amplify + 1;
1941     }
1942
1943     /* set emission map */
1944     clampBoundsInDomain(sds, em->min, em->max, NULL, NULL, (int)ceil(sfs->surface_distance), dt);
1945     em_allocateData(em, sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY, hires_multiplier);
1946
1947     /* setup loop bounds */
1948     for (i = 0; i < 3; i++) {
1949       min[i] = em->min[i] * hires_multiplier;
1950       max[i] = em->max[i] * hires_multiplier;
1951       res[i] = em->res[i] * hires_multiplier;
1952     }
1953
1954     if (BKE_bvhtree_from_mesh_get(&treeData, me, BVHTREE_FROM_LOOPTRI, 4)) {
1955       const float hr = 1.0f / ((float)hires_multiplier);
1956
1957       EmitFromDMData data = {
1958           .sds = sds,
1959           .sfs = sfs,
1960           .mvert = mvert,
1961           .mloop = mloop,
1962           .mlooptri = mlooptri,
1963           .mloopuv = mloopuv,
1964           .dvert = dvert,
1965           .defgrp_index = defgrp_index,
1966           .tree = &treeData,
1967           .hires_multiplier = hires_multiplier,
1968           .hr = hr,
1969           .em = em,
1970           .has_velocity = has_velocity,
1971           .vert_vel = vert_vel,
1972           .flow_center = flow_center,
1973           .min = min,
1974           .max = max,
1975           .res = res,
1976       };
1977
1978       TaskParallelSettings settings;
1979       BLI_parallel_range_settings_defaults(&settings);
1980       settings.scheduling_mode = TASK_SCHEDULING_DYNAMIC;
1981       BLI_task_parallel_range(min[2], max[2], &data, emit_from_mesh_task_cb, &settings);
1982     }
1983     /* free bvh tree */
1984     free_bvhtree_from_mesh(&treeData);
1985
1986     if (vert_vel) {
1987       MEM_freeN(vert_vel);
1988     }
1989
1990     if (me->mvert) {
1991       MEM_freeN(me->mvert);
1992     }
1993     BKE_id_free(NULL, me);
1994   }
1995 }
1996
1997 /**********************************************************
1998  *  Smoke step
1999  **********************************************************/
2000
2001 static void adjustDomainResolution(SmokeDomainSettings *sds,
2002                                    int new_shift[3],
2003                                    EmissionMap *emaps,
2004                                    unsigned int numflowobj,
2005                                    float dt)
2006 {
2007   const int block_size = sds->amplify + 1;
2008   int min[3] = {32767, 32767, 32767}, max[3] = {-32767, -32767, -32767}, res[3];
2009   int total_cells = 1, res_changed = 0, shift_changed = 0;
2010   float min_vel[3], max_vel[3];
2011   int x, y, z;
2012   float *density = smoke_get_density(sds->fluid);
2013   float *fuel = smoke_get_fuel(sds->fluid);
2014   float *bigdensity = smoke_turbulence_get_density(sds->wt);
2015   float *bigfuel = smoke_turbulence_get_fuel(sds->wt);
2016   float *vx = smoke_get_velocity_x(sds->fluid);
2017   float *vy = smoke_get_velocity_y(sds->fluid);
2018   float *vz = smoke_get_velocity_z(sds->fluid);
2019   int wt_res[3];
2020
2021   if (sds->flags & MOD_SMOKE_HIGHRES && sds->wt) {
2022     smoke_turbulence_get_res(sds->wt, wt_res);
2023   }
2024
2025   INIT_MINMAX(min_vel, max_vel);
2026
2027   /* Calculate bounds for current domain content */
2028   for (x = sds->res_min[0]; x < sds->res_max[0]; x++) {
2029     for (y = sds->res_min[1]; y < sds->res_max[1]; y++) {
2030       for (z = sds->res_min[2]; z < sds->res_max[2]; z++) {
2031         int xn = x - new_shift[0];
2032         int yn = y - new_shift[1];
2033         int zn = z - new_shift[2];
2034         int index;
2035         float max_den;
2036
2037         /* skip if cell already belongs to new area */
2038         if (xn >= min[0] && xn <= max[0] && yn >= min[1] && yn <= max[1] && zn >= min[2] &&
2039             zn <= max[2]) {
2040           continue;
2041         }
2042
2043         index = smoke_get_index(x - sds->res_min[0],
2044                                 sds->res[0],
2045                                 y - sds->res_min[1],
2046                                 sds->res[1],
2047                                 z - sds->res_min[2]);
2048         max_den = (fuel) ? MAX2(density[index], fuel[index]) : density[index];
2049
2050         /* check high resolution bounds if max density isnt already high enough */
2051         if (max_den < sds->adapt_threshold && sds->flags & MOD_SMOKE_HIGHRES && sds->wt) {
2052           int i, j, k;
2053           /* high res grid index */
2054           int xx = (x - sds->res_min[0]) * block_size;
2055           int yy = (y - sds->res_min[1]) * block_size;
2056           int zz = (z - sds->res_min[2]) * block_size;
2057
2058           for (i = 0; i < block_size; i++) {
2059             for (j = 0; j < block_size; j++) {
2060               for (k = 0; k < block_size; k++) {
2061                 int big_index = smoke_get_index(xx + i, wt_res[0], yy + j, wt_res[1], zz + k);
2062                 float den = (bigfuel) ? MAX2(bigdensity[big_index], bigfuel[big_index]) :
2063                                         bigdensity[big_index];
2064                 if (den > max_den) {
2065                   max_den = den;
2066                 }
2067               }
2068             }
2069           }
2070         }
2071
2072         /* content bounds (use shifted coordinates) */
2073         if (max_den >= sds->adapt_threshold) {
2074           if (min[0] > xn) {
2075             min[0] = xn;
2076           }
2077           if (min[1] > yn) {
2078             min[1] = yn;
2079           }
2080           if (min[2] > zn) {
2081             min[2] = zn;
2082           }
2083           if (max[0] < xn) {
2084             max[0] = xn;
2085           }
2086           if (max[1] < yn) {
2087             max[1] = yn;
2088           }
2089           if (max[2] < zn) {
2090             max[2] = zn;
2091           }
2092         }
2093
2094         /* velocity bounds */
2095         if (min_vel[0] > vx[index]) {
2096           min_vel[0] = vx[index];
2097         }
2098         if (min_vel[1] > vy[index]) {
2099           min_vel[1] = vy[index];
2100         }
2101         if (min_vel[2] > vz[index]) {
2102           min_vel[2] = vz[index];
2103         }
2104         if (max_vel[0] < vx[index]) {
2105           max_vel[0] = vx[index];
2106         }
2107         if (max_vel[1] < vy[index]) {
2108           max_vel[1] = vy[index];
2109         }
2110         if (max_vel[2] < vz[index]) {
2111           max_vel[2] = vz[index];
2112         }
2113       }
2114     }
2115   }
2116
2117   /* also apply emission maps */
2118   for (int i = 0; i < numflowobj; i++) {
2119     EmissionMap *em = &emaps[i];
2120
2121     for (x = em->min[0]; x < em->max[0]; x++) {
2122       for (y = em->min[1]; y < em->max[1]; y++) {
2123         for (z = em->min[2]; z < em->max[2]; z++) {
2124           int index = smoke_get_index(
2125               x - em->min[0], em->res[0], y - em->min[1], em->res[1], z - em->min[2]);
2126           float max_den = em->influence[index];
2127
2128           /* density bounds */
2129           if (max_den >= sds->adapt_threshold) {
2130             if (min[0] > x) {
2131               min[0] = x;
2132             }
2133             if (min[1] > y) {
2134               min[1] = y;
2135             }
2136             if (min[2] > z) {
2137               min[2] = z;
2138             }
2139             if (max[0] < x) {
2140               max[0] = x;
2141             }
2142             if (max[1] < y) {
2143               max[1] = y;
2144             }
2145             if (max[2] < z) {
2146               max[2] = z;
2147             }
2148           }
2149         }
2150       }
2151     }
2152   }
2153
2154   /* calculate new bounds based on these values */
2155   mul_v3_fl(min_vel, 1.0f / sds->dx);
2156   mul_v3_fl(max_vel, 1.0f / sds->dx);
2157   clampBoundsInDomain(sds, min, max, min_vel, max_vel, sds->adapt_margin + 1, dt);
2158
2159   for (int i = 0; i < 3; i++) {
2160     /* calculate new resolution */
2161     res[i] = max[i] - min[i];
2162     total_cells *= res[i];
2163
2164     if (new_shift[i]) {
2165       shift_changed = 1;
2166     }
2167
2168     /* if no content set minimum dimensions */
2169     if (res[i] <= 0) {
2170       int j;
2171       for (j = 0; j < 3; j++) {
2172         min[j] = 0;
2173         max[j] = 1;
2174         res[j] = 1;
2175       }
2176       res_changed = 1;
2177       total_cells = 1;
2178       break;
2179     }
2180     if (min[i] != sds->res_min[i] || max[i] != sds->res_max[i]) {
2181       res_changed = 1;
2182     }
2183   }
2184
2185   if (res_changed || shift_changed) {
2186     struct FLUID_3D *fluid_old = sds->fluid;
2187     struct WTURBULENCE *turb_old = sds->wt;
2188     /* allocate new fluid data */
2189     BKE_smoke_reallocate_fluid(sds, sds->dx, res, 0);
2190     if (sds->flags & MOD_SMOKE_HIGHRES) {
2191       BKE_smoke_reallocate_highres_fluid(sds, sds->dx, res, 0);
2192     }
2193
2194     /* copy values from old fluid to new */
2195     if (sds->total_cells > 1 && total_cells > 1) {
2196       /* low res smoke */
2197       float *o_dens, *o_react, *o_flame, *o_fuel, *o_heat, *o_heatold, *o_vx, *o_vy, *o_vz, *o_r,
2198           *o_g, *o_b;
2199       float *n_dens, *n_react, *n_flame, *n_fuel, *n_heat, *n_heatold, *n_vx, *n_vy, *n_vz, *n_r,
2200           *n_g, *n_b;
2201       float dummy;
2202       unsigned char *dummy_p;
2203       /* high res smoke */
2204       int wt_res_old[3];
2205       float *o_wt_dens, *o_wt_react, *o_wt_flame, *o_wt_fuel, *o_wt_tcu, *o_wt_tcv, *o_wt_tcw,
2206           *o_wt_r, *o_wt_g, *o_wt_b;
2207       float *n_wt_dens, *n_wt_react, *n_wt_flame, *n_wt_fuel, *n_wt_tcu, *n_wt_tcv, *n_wt_tcw,
2208           *n_wt_r, *n_wt_g, *n_wt_b;
2209
2210       smoke_export(fluid_old,
2211                    &dummy,
2212                    &dummy,
2213                    &o_dens,
2214                    &o_react,
2215                    &o_flame,
2216                    &o_fuel,
2217                    &o_heat,
2218                    &o_heatold,
2219                    &o_vx,
2220                    &o_vy,
2221                    &o_vz,
2222                    &o_r,
2223                    &o_g,
2224                    &o_b,
2225                    &dummy_p);
2226       smoke_export(sds->fluid,
2227                    &dummy,
2228                    &dummy,
2229                    &n_dens,
2230                    &n_react,
2231                    &n_flame,
2232                    &n_fuel,
2233                    &n_heat,
2234                    &n_heatold,
2235                    &n_vx,
2236                    &n_vy,
2237                    &n_vz,
2238                    &n_r,
2239                    &n_g,
2240                    &n_b,
2241                    &dummy_p);
2242
2243       if (sds->flags & MOD_SMOKE_HIGHRES) {
2244         smoke_turbulence_export(turb_old,
2245                                 &o_wt_dens,
2246                                 &o_wt_react,
2247                                 &o_wt_flame,
2248                                 &o_wt_fuel,
2249                                 &o_wt_r,
2250                                 &o_wt_g,
2251                                 &o_wt_b,
2252                                 &o_wt_tcu,
2253                                 &o_wt_tcv,
2254                                 &o_wt_tcw);
2255         smoke_turbulence_get_res(turb_old, wt_res_old);
2256         smoke_turbulence_export(sds->wt,
2257                                 &n_wt_dens,
2258                                 &n_wt_react,
2259                                 &n_wt_flame,
2260                                 &n_wt_fuel,
2261                                 &n_wt_r,
2262                                 &n_wt_g,
2263                                 &n_wt_b,
2264                                 &n_wt_tcu,
2265                                 &n_wt_tcv,
2266                                 &n_wt_tcw);
2267       }
2268
2269       for (x = sds->res_min[0]; x < sds->res_max[0]; x++) {
2270         for (y = sds->res_min[1]; y < sds->res_max[1]; y++) {
2271           for (z = sds->res_min[2]; z < sds->res_max[2]; z++) {
2272             /* old grid index */
2273             int xo = x - sds->res_min[0];
2274             int yo = y - sds->res_min[1];
2275             int zo = z - sds->res_min[2];
2276             int index_old = smoke_get_index(xo, sds->res[0], yo, sds->res[1], zo);
2277             /* new grid index */
2278             int xn = x - min[0] - new_shift[0];
2279             int yn = y - min[1] - new_shift[1];
2280             int zn = z - min[2] - new_shift[2];
2281             int index_new = smoke_get_index(xn, res[0], yn, res[1], zn);
2282
2283             /* skip if outside new domain */
2284             if (xn < 0 || xn >= res[0] || yn < 0 || yn >= res[1] || zn < 0 || zn >= res[2]) {
2285               continue;
2286             }
2287
2288             /* copy data */
2289             n_dens[index_new] = o_dens[index_old];
2290             /* heat */
2291             if (n_heat && o_heat) {
2292               n_heat[index_new] = o_heat[index_old];
2293               n_heatold[index_new] = o_heatold[index_old];
2294             }
2295             /* fuel */
2296             if (n_fuel && o_fuel) {
2297               n_flame[index_new] = o_flame[index_old];
2298               n_fuel[index_new] = o_fuel[index_old];
2299               n_react[index_new] = o_react[index_old];
2300             }
2301             /* color */
2302             if (o_r && n_r) {
2303               n_r[index_new] = o_r[index_old];
2304               n_g[index_new] = o_g[index_old];
2305               n_b[index_new] = o_b[index_old];
2306             }
2307             n_vx[index_new] = o_vx[index_old];
2308             n_vy[index_new] = o_vy[index_old];
2309             n_vz[index_new] = o_vz[index_old];
2310
2311             if (sds->flags & MOD_SMOKE_HIGHRES && turb_old) {
2312               int i, j, k;
2313               /* old grid index */
2314               int xx_o = xo * block_size;
2315               int yy_o = yo * block_size;
2316               int zz_o = zo * block_size;
2317               /* new grid index */
2318               int xx_n = xn * block_size;
2319               int yy_n = yn * block_size;
2320               int zz_n = zn * block_size;
2321
2322               n_wt_tcu[index_new] = o_wt_tcu[index_old];
2323               n_wt_tcv[index_new] = o_wt_tcv[index_old];
2324               n_wt_tcw[index_new] = o_wt_tcw[index_old];
2325
2326               for (i = 0; i < block_size; i++) {
2327                 for (j = 0; j < block_size; j++) {
2328                   for (k = 0; k < block_size; k++) {
2329                     int big_index_old = smoke_get_index(
2330                         xx_o + i, wt_res_old[0], yy_o + j, wt_res_old[1], zz_o + k);
2331                     int big_index_new = smoke_get_index(
2332                         xx_n + i, sds->res_wt[0], yy_n + j, sds->res_wt[1], zz_n + k);
2333                     /* copy data */
2334                     n_wt_dens[big_index_new] = o_wt_dens[big_index_old];
2335                     if (n_wt_flame && o_wt_flame) {
2336                       n_wt_flame[big_index_new] = o_wt_flame[big_index_old];
2337                       n_wt_fuel[big_index_new] = o_wt_fuel[big_index_old];
2338                       n_wt_react[big_index_new] = o_wt_react[big_index_old];
2339                     }
2340                     if (n_wt_r && o_wt_r) {
2341                       n_wt_r[big_index_new] = o_wt_r[big_index_old];
2342                       n_wt_g[big_index_new] = o_wt_g[big_index_old];
2343                       n_wt_b[big_index_new] = o_wt_b[big_index_old];
2344                     }
2345                   }
2346                 }
2347               }
2348             }
2349           }
2350         }
2351       }
2352     }
2353     smoke_free(fluid_old);
2354     if (turb_old) {
2355       smoke_turbulence_free(turb_old);
2356     }
2357
2358     /* set new domain dimensions */
2359     copy_v3_v3_int(sds->res_min, min);
2360     copy_v3_v3_int(sds->res_max, max);
2361     copy_v3_v3_int(sds->res, res);
2362     sds->total_cells = total_cells;
2363   }
2364 }
2365
2366 BLI_INLINE void apply_outflow_fields(int index,
2367                                      float *density,
2368                                      float *heat,
2369                                      float *fuel,
2370                                      float *react,
2371                                      float *color_r,
2372                                      float *color_g,
2373                                      float *color_b)
2374 {
2375   density[index] = 0.f;
2376   if (heat) {
2377     heat[index] = 0.f;
2378   }
2379   if (fuel) {
2380     fuel[index] = 0.f;
2381     react[index] = 0.f;
2382   }
2383   if (color_r) {
2384     color_r[index] = 0.f;
2385     color_g[index] = 0.f;
2386     color_b[index] = 0.f;
2387   }
2388 }
2389
2390 BLI_INLINE void apply_inflow_fields(SmokeFlowSettings *sfs,
2391                                     float emission_value,
2392                                     int index,
2393                                     float *density,
2394                                     float *heat,
2395                                     float *fuel,
2396                                     float *react,
2397                                     float *color_r,
2398                                     float *color_g,
2399                                     float *color_b)
2400 {
2401   int absolute_flow = (sfs->flags & MOD_SMOKE_FLOW_ABSOLUTE);
2402   float dens_old = density[index];
2403   // float fuel_old = (fuel) ? fuel[index] : 0.0f;  /* UNUSED */
2404   float dens_flow = (sfs->type == MOD_SMOKE_FLOW_TYPE_FIRE) ? 0.0f : emission_value * sfs->density;
2405   float fuel_flow = emission_value * sfs->fuel_amount;
2406   /* add heat */
2407   if (heat && emission_value > 0.0f) {
2408     heat[index] = ADD_IF_LOWER(heat[index], sfs->temp);
2409   }
2410   /* absolute */
2411   if (absolute_flow) {
2412     if (sfs->type != MOD_SMOKE_FLOW_TYPE_FIRE) {
2413       if (dens_flow > density[index]) {
2414         density[index] = dens_flow;
2415       }
2416     }
2417     if (sfs->type != MOD_SMOKE_FLOW_TYPE_SMOKE && fuel && fuel_flow) {
2418       if (fuel_flow > fuel[index]) {
2419         fuel[index] = fuel_flow;
2420       }
2421     }
2422   }
2423   /* additive */
2424   else {
2425     if (sfs->type != MOD_SMOKE_FLOW_TYPE_FIRE) {
2426       density[index] += dens_flow;
2427       CLAMP(density[index], 0.0f, 1.0f);
2428     }
2429     if (sfs->type != MOD_SMOKE_FLOW_TYPE_SMOKE && fuel && sfs->fuel_amount) {
2430       fuel[index] += fuel_flow;
2431       CLAMP(fuel[index], 0.0f, 10.0f);
2432     }
2433   }
2434
2435   /* set color */
2436   if (color_r && dens_flow) {
2437     float total_dens = density[index] / (dens_old + dens_flow);
2438     color_r[index] = (color_r[index] + sfs->color[0] * dens_flow) * total_dens;
2439     color_g[index] = (color_g[index] + sfs->color[1] * dens_flow) * total_dens;
2440     color_b[index] = (color_b[index] + sfs->color[2] * dens_flow) * total_dens;
2441   }
2442
2443   /* set fire reaction coordinate */
2444   if (fuel && fuel[index] > FLT_EPSILON) {
2445     /* instead of using 1.0 for all new fuel add slight falloff
2446      * to reduce flow blockiness */
2447     float value = 1.0f - pow2f(1.0f - emission_value);
2448
2449     if (value > react[index]) {
2450       float f = fuel_flow / fuel[index];
2451       react[index] = value * f + (1.0f - f) * react[index];
2452       CLAMP(react[index], 0.0f, value);
2453     }
2454   }
2455 }
2456
2457 static void update_flowsfluids(
2458     Depsgraph *depsgraph, Scene *scene, Object *ob, SmokeDomainSettings *sds, float dt)
2459 {
2460   Object **flowobjs = NULL;
2461   EmissionMap *emaps = NULL;
2462   unsigned int numflowobj = 0;
2463   unsigned int flowIndex;
2464   int new_shift[3] = {0};
2465   int active_fields = sds->active_fields;
2466
2467   /* calculate domain shift for current frame if using adaptive domain */
2468   if (sds->flags & MOD_SMOKE_ADAPTIVE_DOMAIN) {
2469     int total_shift[3];
2470     float frame_shift_f[3];
2471     float ob_loc[3] = {0};
2472
2473     mul_m4_v3(ob->obmat, ob_loc);
2474
2475     sub_v3_v3v3(frame_shift_f, ob_loc, sds->prev_loc);
2476     copy_v3_v3(sds->prev_loc, ob_loc);
2477     /* convert global space shift to local "cell" space */
2478     mul_mat3_m4_v3(sds->imat, frame_shift_f);
2479     frame_shift_f[0] = frame_shift_f[0] / sds->cell_size[0];
2480     frame_shift_f[1] = frame_shift_f[1] / sds->cell_size[1];
2481     frame_shift_f[2] = frame_shift_f[2] / sds->cell_size[2];
2482     /* add to total shift */
2483     add_v3_v3(sds->shift_f, frame_shift_f);
2484     /* convert to integer */
2485     total_shift[0] = (int)(floorf(sds->shift_f[0]));
2486     total_shift[1] = (int)(floorf(sds->shift_f[1]));
2487     total_shift[2] = (int)(floorf(sds->shift_f[2]));
2488     sub_v3_v3v3_int(new_shift, total_shift, sds->shift);
2489     copy_v3_v3_int(sds->shift, total_shift);
2490
2491     /* calculate new domain boundary points so that smoke doesn't slide on sub-cell movement */
2492     sds->p0[0] = sds->dp0[0] - sds->cell_size[0] * (sds->shift_f[0] - total_shift[0] - 0.5f);
2493     sds->p0[1] = sds->dp0[1] - sds->cell_size[1] * (sds->shift_f[1] - total_shift[1] - 0.5f);
2494     sds->p0[2] = sds->dp0[2] - sds->cell_size[2] * (sds->shift_f[2] - total_shift[2] - 0.5f);
2495     sds->p1[0] = sds->p0[0] + sds->cell_size[0] * sds->base_res[0];
2496     sds->p1[1] = sds->p0[1] + sds->cell_size[1] * sds->base_res[1];
2497     sds->p1[2] = sds->p0[2] + sds->cell_size[2] * sds->base_res[2];
2498   }
2499
2500   flowobjs = BKE_collision_objects_create(
2501       depsgraph, ob, sds->fluid_group, &numflowobj, eModifierType_Smoke);
2502
2503   /* init emission maps for each flow */
2504   emaps = MEM_callocN(sizeof(struct EmissionMap) * numflowobj, "smoke_flow_maps");
2505
2506   /* Prepare flow emission maps */
2507   for (flowIndex = 0; flowIndex < numflowobj; flowIndex++) {
2508     Object *collob = flowobjs[flowIndex];
2509     SmokeModifierData *smd2 = (SmokeModifierData *)modifiers_findByType(collob,
2510                                                                         eModifierType_Smoke);
2511
2512     // check for initialized smoke object
2513     if ((smd2->type & MOD_SMOKE_TYPE_FLOW) && smd2->flow) {
2514       // we got nice flow object
2515       SmokeFlowSettings *sfs = smd2->flow;
2516       int subframes = sfs->subframes;
2517       EmissionMap *em = &emaps[flowIndex];
2518
2519       /* just sample flow directly to emission map if no subframes */
2520       if (!subframes) {
2521         if (sfs->source == MOD_SMOKE_FLOW_SOURCE_PARTICLES) {
2522           emit_from_particles(collob, sds, sfs, em, depsgraph, scene, dt);
2523         }
2524         else {
2525           emit_from_mesh(collob, sds, sfs, em, dt);
2526         }
2527       }
2528       /* sample subframes */
2529       else {
2530 #  if 0
2531         int scene_frame = (int)DEG_get_ctime(depsgraph);
2532 #  endif
2533         // float scene_subframe = scene->r.subframe;  // UNUSED
2534         int subframe;
2535         for (subframe = 0; subframe <= subframes; subframe++) {
2536           EmissionMap em_temp = {NULL};
2537           float sample_size = 1.0f / (float)(subframes + 1);
2538 #  if 0
2539           float prev_frame_pos = sample_size * (float)(subframe + 1);
2540 #  endif
2541           float sdt = dt * sample_size;
2542           int hires_multiplier = 1;
2543
2544           if ((sds->flags & MOD_SMOKE_HIGHRES) && (sds->highres_sampling == SM_HRES_FULLSAMPLE)) {
2545             hires_multiplier = sds->amplify + 1;
2546           }
2547
2548           /* TODO: setting the scene frame no longer works with the new depsgraph. */
2549 #  if 0
2550           /* set scene frame to match previous frame + subframe
2551            * or use current frame for last sample */
2552           if (subframe < subframes) {
2553             scene->r.cfra = scene_frame - 1;
2554             scene->r.subframe = prev_frame_pos;
2555           }
2556           else {
2557             scene->r.cfra = scene_frame;
2558             scene->r.subframe = 0.0f;
2559           }
2560 #  endif
2561
2562           if (sfs->source == MOD_SMOKE_FLOW_SOURCE_PARTICLES) {
2563             /* emit_from_particles() updates timestep internally */
2564             emit_from_particles(collob, sds, sfs, &em_temp, depsgraph, scene, sdt);
2565             if (!(sfs->flags & MOD_SMOKE_FLOW_USE_PART_SIZE)) {
2566               hires_multiplier = 1;
2567             }
2568           }
2569           else { /* MOD_SMOKE_FLOW_SOURCE_MESH */
2570             /* update flow object frame */
2571             BLI_mutex_lock(&object_update_lock);
2572             BKE_object_modifier_update_subframe(
2573                 depsgraph, scene, collob, true, 5, DEG_get_ctime(depsgraph), eModifierType_Smoke);
2574             BLI_mutex_unlock(&object_update_lock);
2575
2576             /* apply flow */
2577             emit_from_mesh(collob, sds, sfs, &em_temp, sdt);
2578           }
2579
2580           /* combine emission maps */
2581           em_combineMaps(em,
2582                          &em_temp,
2583                          hires_multiplier,
2584                          !(sfs->flags & MOD_SMOKE_FLOW_ABSOLUTE),
2585                          sample_size);
2586           em_freeData(&em_temp);
2587         }
2588       }
2589
2590       /* update required data fields */
2591       if (em->total_cells && sfs->type != MOD_SMOKE_FLOW_TYPE_OUTFLOW) {
2592         /* activate heat field if flow produces any heat */
2593         if (sfs->temp) {
2594           active_fields |= SM_ACTIVE_HEAT;
2595         }
2596         /* activate fuel field if flow adds any fuel */
2597         if (sfs->type != MOD_SMOKE_FLOW_TYPE_SMOKE && sfs->fuel_amount) {
2598           active_fields |= SM_ACTIVE_FIRE;
2599         }
2600         /* activate color field if flows add smoke with varying colors */
2601         if (sfs->type != MOD_SMOKE_FLOW_TYPE_FIRE && sfs->density) {
2602           if (!(active_fields & SM_ACTIVE_COLOR_SET)) {
2603             copy_v3_v3(sds->active_color, sfs->color);
2604             active_fields |= SM_ACTIVE_COLOR_SET;
2605           }
2606           else if (!equals_v3v3(sds->active_color, sfs->color)) {
2607             copy_v3_v3(sds->active_color, sfs->color);
2608             active_fields |= SM_ACTIVE_COLORS;
2609           }
2610         }
2611       }
2612     }
2613   }
2614
2615   /* monitor active fields based on domain settings */
2616   /* if domain has fire, activate new fields if required */
2617   if (active_fields & SM_ACTIVE_FIRE) {
2618     /* heat is always needed for fire */
2619     active_fields |= SM_ACTIVE_HEAT;
2620     /* also activate colors if domain smoke color differs from active color */
2621     if (!(active_fields & SM_ACTIVE_COLOR_SET)) {
2622       copy_v3_v3(sds->active_color, sds->flame_smoke_color);
2623       active_fields |= SM_ACTIVE_COLOR_SET;
2624     }
2625     else if (!equals_v3v3(sds->active_color, sds->flame_smoke_color)) {
2626       copy_v3_v3(sds->active_color, sds->flame_smoke_color);
2627       active_fields |= SM_ACTIVE_COLORS;
2628     }
2629   }
2630
2631   /* Adjust domain size if needed */
2632   if (sds->flags & MOD_SMOKE_ADAPTIVE_DOMAIN) {
2633     adjustDomainResolution(sds, new_shift, emaps, numflowobj, dt);
2634   }
2635
2636   /* Initialize new data fields if any */
2637   if (active_fields & SM_ACTIVE_HEAT) {
2638     smoke_ensure_heat(sds->fluid);
2639   }
2640   if (active_fields & SM_ACTIVE_FIRE) {
2641     smoke_ensure_fire(sds->fluid, sds->wt);
2642   }
2643   if (active_fields & SM_ACTIVE_COLORS) {
2644     /* initialize all smoke with "active_color" */
2645     smoke_ensure_colors(
2646         sds->fluid, sds->wt, sds->active_color[0], sds->active_color[1], sds->active_color[2]);
2647   }
2648   sds->active_fields = active_fields;
2649
2650   /* Apply emission data */
2651   if (sds->fluid) {
2652     for (flowIndex = 0; flowIndex < numflowobj; flowIndex++) {
2653       Object *collob = flowobjs[flowIndex];
2654       SmokeModifierData *smd2 = (SmokeModifierData *)modifiers_findByType(collob,
2655                                                                           eModifierType_Smoke);
2656
2657       // check for initialized smoke object
2658       if ((smd2->type & MOD_SMOKE_TYPE_FLOW) && smd2->flow) {
2659         // we got nice flow object
2660         SmokeFlowSettings *sfs = smd2->flow;
2661         EmissionMap *em = &emaps[flowIndex];
2662
2663         float *density = smoke_get_density(sds->fluid);
2664         float *color_r = smoke_get_color_r(sds->fluid);
2665         float *color_g = smoke_get_color_g(sds->fluid);
2666         float *color_b = smoke_get_color_b(sds->fluid);
2667         float *fuel = smoke_get_fuel(sds->fluid);
2668         float *react = smoke_get_react(sds->fluid);
2669         float *bigdensity = smoke_turbulence_get_density(sds->wt);
2670         float *bigfuel = smoke_turbulence_get_fuel(sds->wt);
2671         float *bigreact = smoke_turbulence_get_react(sds->wt);
2672         float *bigcolor_r = smoke_turbulence_get_color_r(sds->wt);
2673         float *bigcolor_g = smoke_turbulence_get_color_g(sds->wt);
2674         float *bigcolor_b = smoke_turbulence_get_color_b(sds->wt);
2675         float *heat = smoke_get_heat(sds->fluid);
2676         float *velocity_x = smoke_get_velocity_x(sds->fluid);
2677         float *velocity_y = smoke_get_velocity_y(sds->fluid);
2678         float *velocity_z = smoke_get_velocity_z(sds->fluid);
2679         // unsigned char *obstacle = smoke_get_obstacle(sds->fluid);
2680         // DG TODO UNUSED unsigned char *obstacleAnim = smoke_get_obstacle_anim(sds->fluid);
2681         int bigres[3];
2682         float *velocity_map = em->velocity;
2683         float *emission_map = em->influence;
2684         float *emission_map_high = em->influence_high;
2685
2686         int ii, jj, kk, gx, gy, gz, ex, ey, ez, dx, dy, dz, block_size;
2687         size_t e_index, d_index, index_big;
2688
2689         // loop through every emission map cell
2690         for (gx = em->min[0]; gx < em->max[0]; gx++) {
2691           for (gy = em->min[1]; gy < em->max[1]; gy++) {
2692             for (gz = em->min[2]; gz < em->max[2]; gz++) {
2693               /* get emission map index */
2694               ex = gx - em->min[0];
2695               ey = gy - em->min[1];
2696               ez = gz - em->min[2];
2697               e_index = smoke_get_index(ex, em->res[0], ey, em->res[1], ez);
2698
2699               /* get domain index */
2700               dx = gx - sds->res_min[0];
2701               dy = gy - sds->res_min[1];
2702               dz = gz - sds->res_min[2];
2703               d_index = smoke_get_index(dx, sds->res[0], dy, sds->res[1], dz);
2704               /* make sure emission cell is inside the new domain boundary */
2705               if (dx < 0 || dy < 0 || dz < 0 || dx >= sds->res[0] || dy >= sds->res[1] ||
2706                   dz >= sds->res[2]) {
2707                 continue;
2708               }
2709
2710               if (sfs->type == MOD_SMOKE_FLOW_TYPE_OUTFLOW) {  // outflow
2711                 apply_outflow_fields(
2712                     d_index, density, heat, fuel, react, color_r, color_g, color_b);
2713               }
2714               else {  // inflow
2715                 apply_inflow_fields(sfs,
2716                                     emission_map[e_index],
2717                                     d_index,
2718                                     density,
2719                                     heat,
2720                                     fuel,
2721                                     react,
2722                                     color_r,
2723                                     color_g,
2724                                     color_b);
2725
2726                 /* initial velocity */
2727                 if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY) {
2728                   velocity_x[d_index] = ADD_IF_LOWER(velocity_x[d_index],
2729                                                      velocity_map[e_index * 3]);
2730                   velocity_y[d_index] = ADD_IF_LOWER(velocity_y[d_index],
2731                                                      velocity_map[e_index * 3 + 1]);
2732                   velocity_z[d_index] = ADD_IF_LOWER(velocity_z[d_index],
2733                                                      velocity_map[e_index * 3 + 2]);
2734                 }
2735               }
2736
2737               /* loop through high res blocks if high res enabled */
2738               if (bigdensity) {
2739                 /* Neighbor cell emission densities
2740                  * (for high resolution smoke smooth interpolation). */
2741                 float c000, c001, c010, c011, c100, c101, c110, c111;
2742
2743                 smoke_turbulence_get_res(sds->wt, bigres);
2744                 block_size = sds->amplify + 1;  // high res block size
2745
2746                 c000 = (ex > 0 && ey > 0 && ez > 0) ?
2747                            emission_map[smoke_get_index(
2748                                ex - 1, em->res[0], ey - 1, em->res[1], ez - 1)] :
2749                            0;
2750                 c001 =
2751                     (ex > 0 && ey > 0) ?
2752                         emission_map[smoke_get_index(ex - 1, em->res[0], ey - 1, em->res[1], ez)] :
2753                         0;
2754                 c010 =
2755                     (ex > 0 && ez > 0) ?
2756                         emission_map[smoke_get_index(ex - 1, em->res[0], ey, em->res[1], ez - 1)] :
2757                         0;
2758                 c011 = (ex > 0) ?
2759                            emission_map[smoke_get_index(ex - 1, em->res[0], ey, em->res[1], ez)] :
2760                            0;
2761
2762                 c100 =
2763                     (ey > 0 && ez > 0) ?
2764                         emission_map[smoke_get_index(ex, em->res[0], ey - 1, em->res[1], ez - 1)] :
2765                         0;
2766                 c101 = (ey > 0) ?
2767                            emission_map[smoke_get_index(ex, em->res[0], ey - 1, em->res[1], ez)] :
2768                            0;
2769                 c110 = (ez > 0) ?
2770                            emission_map[smoke_get_index(ex, em->res[0], ey, em->res[1], ez - 1)] :
2771                            0;
2772                 c111 = emission_map[smoke_get_index(
2773                     ex, em->res[0], ey, em->res[1], ez)];  // this cell
2774
2775                 for (ii = 0; ii < block_size; ii++) {
2776                   for (jj = 0; jj < block_size; jj++) {
2777                     for (kk = 0; kk < block_size; kk++) {
2778
2779                       float fx, fy, fz, interpolated_value;
2780                       int shift_x = 0, shift_y = 0, shift_z = 0;
2781
2782                       /* Use full sample emission map if enabled and available */
2783                       if ((sds->highres_sampling == SM_HRES_FULLSAMPLE) && emission_map_high) {
2784                         interpolated_value =
2785                             emission_map_high[smoke_get_index(ex * block_size + ii,
2786                                                               em->res[0] * block_size,
2787                                                               ey * block_size + jj,
2788                                                               em->res[1] * block_size,
2789                                                               ez * block_size + kk)];  // this cell
2790                       }
2791                       else if (sds->highres_sampling == SM_HRES_NEAREST) {
2792                         /* without interpolation use same low resolution
2793                          * block value for all hi-res blocks */
2794                         interpolated_value = c111;
2795                       }
2796                       /* Fall back to interpolated */
2797                       else {
2798                         /* get relative block position
2799                          * for interpolation smoothing */
2800                         fx = (float)ii / block_size + 0.5f / block_size;
2801                         fy = (float)jj / block_size + 0.5f / block_size;
2802                         fz = (float)kk / block_size + 0.5f / block_size;
2803
2804                         /* calculate trilinear interpolation */
2805                         interpolated_value = c000 * (1 - fx) * (1 - fy) * (1 - fz) +
2806                                              c100 * fx * (1 - fy) * (1 - fz) +
2807                                              c010 * (1 - fx) * fy * (1 - fz) +
2808                                              c001 * (1 - fx) * (1 - fy) * fz +
2809                                              c101 * fx * (1 - fy) * fz +
2810                                              c011 * (1 - fx) * fy * fz +
2811                                              c110 * fx * fy * (1 - fz) + c111 * fx * fy * fz;
2812
2813                         /* add some contrast / sharpness
2814                          * depending on hi-res block size */
2815                         interpolated_value = (interpolated_value - 0.4f) * (block_size / 2) + 0.4f;
2816                         CLAMP(interpolated_value, 0.0f, 1.0f);
2817
2818                         /* shift smoke block index
2819                          * (because pixel center is actually
2820                          * in halfway of the low res block) */
2821                         shift_x = (dx < 1) ? 0 : block_size / 2;
2822                         shift_y = (dy < 1) ? 0 : block_size / 2;
2823                         shift_z = (dz < 1) ? 0 : block_size / 2;
2824                       }
2825
2826                       /* get shifted index for current high resolution block */
2827                       index_big = smoke_get_index(block_size * dx + ii - shift_x,
2828                                                   bigres[0],
2829                                                   block_size * dy + jj - shift_y,
2830                                                   bigres[1],
2831                                                   block_size * dz + kk - shift_z);
2832
2833                       if (sfs->type == MOD_SMOKE_FLOW_TYPE_OUTFLOW) {  // outflow
2834                         if (interpolated_value) {
2835                           apply_outflow_fields(index_big,
2836                                                bigdensity,
2837                                                NULL,
2838                                                bigfuel,
2839                                                bigreact,
2840                                                bigcolor_r,
2841                                                bigcolor_g,
2842                                                bigcolor_b);
2843                         }
2844                       }
2845                       else {  // inflow
2846                         apply_inflow_fields(sfs,
2847                                             interpolated_value,
2848                                             index_big,
2849                                             bigdensity,
2850                                             NULL,
2851                                             bigfuel,
2852                                             bigreact,
2853                                             bigcolor_r,
2854                                             bigcolor_g,
2855                                             bigcolor_b);
2856                       }
2857                     }  // hires loop
2858                   }
2859                 }
2860               }  // bigdensity
2861             }    // low res loop
2862           }
2863         }
2864
2865         // free emission maps
2866         em_freeData(em);
2867
2868       }  // end emission
2869     }
2870   }
2871
2872   BKE_collision_objects_free(flowobjs);
2873   if (emaps) {
2874     MEM_freeN(emaps);
2875   }
2876 }
2877
2878 typedef struct UpdateEffectorsData {
2879   Scene *scene;
2880   SmokeDomainSettings *sds;
2881   ListBase *effectors;
2882
2883   float *density;
2884   float *fuel;
2885   float *force_x;
2886   float *force_y;
2887   float *force_z;
2888   float *velocity_x;
2889   float *velocity_y;
2890   float *velocity_z;
2891   unsigned char *obstacle;
2892 } UpdateEffectorsData;
2893
2894 static void update_effectors_task_cb(void *__restrict userdata,
2895                                      const int x,
2896                                      const TaskParallelTLS *__restrict UNUSED(tls))
2897 {
2898   UpdateEffectorsData *data = userdata;
2899   SmokeDomainSettings *sds = data->sds;
2900
2901   for (int y = 0; y < sds->res[1]; y++) {
2902     for (int z = 0; z < sds->res[2]; z++) {
2903       EffectedPoint epoint;
2904       float mag;
2905       float voxelCenter[3] = {0, 0, 0}, vel[3] = {0, 0, 0}, retvel[3] = {0, 0, 0};
2906       const unsigned int index = smoke_get_index(x, sds->res[0], y, sds->res[1], z);
2907
2908       if (((data->fuel ? MAX2(data->density[index], data->fuel[index]) : data->density[index]) <
2909            FLT_EPSILON) ||
2910           data->obstacle[index]) {
2911         continue;
2912       }
2913
2914       vel[0] = data->velocity_x[index];
2915       vel[1] = data->velocity_y[index];
2916       vel[2] = data->velocity_z[index];
2917
2918       /* convert vel to global space */
2919       mag = len_v3(vel);
2920       mul_mat3_m4_v3(sds->obmat, vel);
2921       normalize_v3(vel);
2922       mul_v3_fl(vel, mag);
2923
2924       voxelCenter[0] = sds->p0[0] + sds->cell_size[0] * ((float)(x + sds->res_min[0]) + 0.5f);
2925       voxelCenter[1] = sds->p0[1] + sds->cell_size[1] * ((float)(y + sds->res_min[1]) + 0.5f);
2926       voxelCenter[2] = sds->p0[2] + sds->cell_size[2] * ((float)(z + sds->res_min[2]) + 0.5f);
2927       mul_m4_v3(sds->obmat, voxelCenter);
2928
2929       pd_point_from_loc(data->scene, voxelCenter, vel, index, &epoint);
2930       BKE_effectors_apply(data->effectors, NULL, sds->effector_weights, &epoint, retvel, NULL);
2931
2932       /* convert retvel to local space */
2933       mag = len_v3(retvel);
2934       mul_mat3_m4_v3(sds->imat, retvel);
2935       normalize_v3(retvel);
2936       mul_v3_fl(retvel, mag);
2937
2938       // TODO dg - do in force!
2939       data->force_x[index] = min_ff(max_ff(-1.0f, retvel[0] * 0.2f), 1.0f);
2940       data->force_y[index] = min_ff(max_ff(-1.0f, retvel[1] * 0.2f), 1.0f);
2941       data->force_z[index] = min_ff(max_ff(-1.0f, retvel[2] * 0.2f), 1.0f);
2942     }
2943   }
2944 }
2945
2946 static void update_effectors(
2947     Depsgraph *depsgraph, Scene *scene, Object *ob, SmokeDomainSettings *sds, float UNUSED(dt))
2948 {
2949   ListBase *effectors;
2950   /* make sure smoke flow influence is 0.0f */
2951   sds->effector_weights->weight[PFIELD_SMOKEFLOW] = 0.0f;
2952   effectors = BKE_effectors_create(depsgraph, ob, NULL, sds->effector_weights);
2953
2954   if (effectors) {
2955     // precalculate wind forces
2956     UpdateEffectorsData data;
2957     data.scene = scene;
2958     data.sds = sds;
2959     data.effectors = effectors;
2960     data.density = smoke_get_density(sds->fluid);
2961     data.fuel = smoke_get_fuel(sds->fluid);
2962     data.force_x = smoke_get_force_x(sds->fluid);
2963     data.force_y = smoke_get_force_y(sds->fluid);
2964     data.force_z = smoke_get_force_z(sds->fluid);
2965     data.velocity_x = smoke_get_velocity_x(sds->fluid);
2966     data.velocity_y = smoke_get_velocity_y(sds->fluid);
2967     data.velocity_z = smoke_get_velocity_z(sds->fluid);
2968     data.obstacle = smoke_get_obstacle(sds->fluid);
2969
2970     TaskParallelSettings settings;
2971     BLI_parallel_range_settings_defaults(&settings);
2972     settings.scheduling_mode = TASK_SCHEDULING_DYNAMIC;
2973     BLI_task_parallel_range(0, sds->res[0], &data, update_effectors_task_cb, &settings);
2974   }
2975
2976   BKE_effectors_free(effectors);
2977 }
2978
2979 static void step(Depsgraph *depsgraph,
2980                  Scene *scene,
2981                  Object *ob,
2982                  SmokeModifierData *smd,
2983                  Mesh *domain_me,
2984                  float fps)
2985 {
2986   SmokeDomainSettings *sds = smd->domain;
2987   /* stability values copied from wturbulence.cpp */
2988   const int maxSubSteps = 25;
2989   float maxVel;
2990   // maxVel should be 1.5 (1.5 cell max movement) * dx (cell size)
2991
2992   float dt;
2993   float maxVelMag = 0.0f;
2994   int totalSubsteps;
2995   int substep = 0;
2996   float dtSubdiv;
2997   float gravity[3] = {0.0f, 0.0f, -1.0f};
2998   float gravity_mag;
2999
3000   /* update object state */
3001   invert_m4_m4(sds->imat, ob->obmat);
3002   copy_m4_m4(sds->obmat, ob->obmat);
3003   smoke_set_domain_from_mesh(sds, ob, domain_me, (sds->flags & MOD_SMOKE_ADAPTIVE_DOMAIN) != 0);
3004
3005   /* use global gravity if enabled */
3006   if (scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY) {
3007     copy_v3_v3(gravity, scene->physics_settings.gravity);
3008     /* map default value to 1.0 */
3009     mul_v3_fl(gravity, 1.0f / 9.810f);
3010   }
3011   /* convert gravity to domain space */
3012   gravity_mag = len_v3(gravity);
3013   mul_mat3_m4_v3(sds->imat, gravity);
3014   normalize_v3(gravity);
3015   mul_v3_fl(gravity, gravity_mag);
3016
3017   /* adapt timestep for different framerates, dt = 0.1 is at 25fps */
3018   dt = DT_DEFAULT * (25.0f / fps);
3019   // maximum timestep/"CFL" constraint: dt < 5.0 *dx / maxVel
3020   maxVel = (sds->dx * 5.0f);
3021
3022   maxVelMag = sqrtf(maxVelMag) * dt * sds->time_scale;
3023   totalSubsteps = (int)((maxVelMag / maxVel) + 1.0f); /* always round up */
3024   totalSubsteps = (totalSubsteps < 1) ? 1 : totalSubsteps;
3025   totalSubsteps = (totalSubsteps > maxSubSteps) ? maxSubSteps : totalSubsteps;
3026
3027   /* Disable substeps for now, since it results in numerical instability */
3028   totalSubsteps = 1.0f;
3029
3030   dtSubdiv = (float)dt / (float)totalSubsteps;
3031
3032   // printf("totalSubsteps: %d, maxVelMag: %f, dt: %f\n", totalSubsteps, maxVelMag, dt);
3033
3034   for (substep = 0; substep < totalSubsteps; substep++) {
3035     // calc animated obstacle velocities
3036     update_flowsfluids(depsgraph, scene, ob, sds, dtSubdiv);
3037     update_obstacles(depsgraph, ob, sds, dtSubdiv, substep, totalSubsteps);
3038
3039     if (sds->total_cells > 1) {
3040       // DG TODO? problem --> uses forces instead of velocity,
3041       // need to check how they need to be changed with variable dt.
3042       update_effectors(depsgraph, scene, ob, sds, dtSubdiv);
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
3490   /* convert light pos to sim cell space */