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