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