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