Smoke simulator: Add flow subframes and ability to set custom particle size.

Previously it was nearly impossible to have fast moving objects emitting smoke or they would just leave behind a row of smoke poofs instead of continious stream of smoke. Now it's possible to set number of subframes for each smoke flow.

Another new thing is ability to set size of smoke flow particles instead of using closest smoke cell. This also works with my earlier "full sample" commit, so no more blocky particles either. :)

For more info check my blog post: http://www.miikahweb.com/en/blog/2013/05/17/blender-smoke-subframes

This commit also includes couple of fixes I spotted while testing:
* Fix: dissolve was applied at different time for low res and high res simulations.
* Fix: full sample setting didn't get copied with domain.

1 files changed, 424 insertions, 58 deletions

diff --git a/source/blender/blenkernel/intern/smoke.c b/source/blender/blenkernel/intern/smoke.c
index 1c66c5ccf77..6ccbfe9378b 100644
--- a/source/blender/blenkernel/intern/smoke.c
+++ b/source/blender/blenkernel/intern/smoke.c
@@ -54,6 +54,9 @@
 #include "BLI_utildefines.h"
 #include "BLI_voxel.h"
 
+#include "DNA_anim_types.h"
+#include "DNA_armature_types.h"
+#include "DNA_constraint_types.h"
 #include "DNA_customdata_types.h"
 #include "DNA_group_types.h"
 #include "DNA_lamp_types.h"
@@ -65,16 +68,21 @@
 #include "DNA_scene_types.h"
 #include "DNA_smoke_types.h"
 
+#include "BKE_animsys.h"
+#include "BKE_armature.h"
 #include "BKE_bvhutils.h"
 #include "BKE_cdderivedmesh.h"
 #include "BKE_collision.h"
+#include "BKE_constraint.h"
 #include "BKE_customdata.h"
 #include "BKE_deform.h"
 #include "BKE_DerivedMesh.h"
 #include "BKE_effect.h"
 #include "BKE_modifier.h"
+#include "BKE_object.h"
 #include "BKE_particle.h"
 #include "BKE_pointcache.h"
+#include "BKE_scene.h"
 #include "BKE_smoke.h"
 
 #include "RE_shader_ext.h"
@@ -561,11 +569,14 @@ void smokeModifier_createType(struct SmokeModifierData *smd)
 			smd->flow->density = 1.0f;
 			smd->flow->fuel_amount = 1.0f;
 			smd->flow->temp = 1.0f;
-			smd->flow->flags = MOD_SMOKE_FLOW_ABSOLUTE;
+			smd->flow->flags = MOD_SMOKE_FLOW_ABSOLUTE | MOD_SMOKE_FLOW_USE_PART_SIZE;
 			smd->flow->vel_multi = 1.0f;
+			smd->flow->volume_density = 0.0f;
 			smd->flow->surface_distance = 1.5f;
 			smd->flow->source = MOD_SMOKE_FLOW_SOURCE_MESH;
 			smd->flow->texture_size = 1.0f;
+			smd->flow->particle_size = 1.0f;
+			smd->flow->subframes = 0;
 
 			smd->flow->color[0] = 0.7f;
 			smd->flow->color[1] = 0.7f;
@@ -615,6 +626,7 @@ void smokeModifier_copy(struct SmokeModifierData *smd, struct SmokeModifierData
 		tsmd->domain->amplify = smd->domain->amplify;
 		tsmd->domain->maxres = smd->domain->maxres;
 		tsmd->domain->flags = smd->domain->flags;
+		tsmd->domain->highres_sampling = smd->domain->highres_sampling;
 		tsmd->domain->viewsettings = smd->domain->viewsettings;
 		tsmd->domain->noise = smd->domain->noise;
 		tsmd->domain->diss_speed = smd->domain->diss_speed;
@@ -648,6 +660,8 @@ void smokeModifier_copy(struct SmokeModifierData *smd, struct SmokeModifierData
 		tsmd->flow->temp = smd->flow->temp;
 		tsmd->flow->volume_density = smd->flow->volume_density;
 		tsmd->flow->surface_distance = smd->flow->surface_distance;
+		tsmd->flow->particle_size = smd->flow->particle_size;
+		tsmd->flow->subframes = smd->flow->subframes;
 
 		tsmd->flow->texture_size = smd->flow->texture_size;
 		tsmd->flow->texture_offset = smd->flow->texture_offset;
@@ -694,6 +708,10 @@ static int get_lamp(Scene *scene, float *light)
 	return found_lamp;
 }
 
+/**********************************************************
+ *	Obstacles
+ **********************************************************/
+
 static void obstacles_from_derivedmesh(Object *coll_ob, SmokeDomainSettings *sds, SmokeCollSettings *scs, unsigned char *obstacle_map, float *velocityX, float *velocityY, float *velocityZ, float dt)
 {
 	if (!scs->dm) return;
@@ -898,11 +916,107 @@ static void update_obstacles(Scene *scene, Object *ob, SmokeDomainSettings *sds,
 }
 
 
+/**********************************************************
+ *	Object subframe update method from dynamicpaint.c
+ **********************************************************/
+
+/* set "ignore cache" flag for all caches on this object */
+static void object_cacheIgnoreClear(Object *ob, int state)
+{
+	ListBase pidlist;
+	PTCacheID *pid;
+	BKE_ptcache_ids_from_object(&pidlist, ob, NULL, 0);
+
+	for (pid = pidlist.first; pid; pid = pid->next) {
+		if (pid->cache) {
+			if (state)
+				pid->cache->flag |= PTCACHE_IGNORE_CLEAR;
+			else
+				pid->cache->flag &= ~PTCACHE_IGNORE_CLEAR;
+		}
+	}
+
+	BLI_freelistN(&pidlist);
+}
+
+static int subframe_updateObject(Scene *scene, Object *ob, int update_mesh, int parent_recursion, float frame)
+{
+	SmokeModifierData *smd = (SmokeModifierData *)modifiers_findByType(ob, eModifierType_Smoke);
+	bConstraint *con;
+
+	/* if other is dynamic paint canvas, don't update */
+	if (smd && (smd->type & MOD_SMOKE_TYPE_DOMAIN))
+		return 1;
+
+	/* if object has parents, update them too */
+	if (parent_recursion) {
+		int recursion = parent_recursion - 1;
+		int is_domain = 0;
+		if (ob->parent) is_domain += subframe_updateObject(scene, ob->parent, 0, recursion, frame);
+		if (ob->track) is_domain += subframe_updateObject(scene, ob->track, 0, recursion, frame);
+
+		/* skip subframe if object is parented
+		 *  to vertex of a dynamic paint canvas */
+		if (is_domain && (ob->partype == PARVERT1 || ob->partype == PARVERT3))
+			return 0;
+
+		/* also update constraint targets */
+		for (con = ob->constraints.first; con; con = con->next) {
+			bConstraintTypeInfo *cti = BKE_constraint_get_typeinfo(con);
+			ListBase targets = {NULL, NULL};
+
+			if (cti && cti->get_constraint_targets) {
+				bConstraintTarget *ct;
+				cti->get_constraint_targets(con, &targets);
+				for (ct = targets.first; ct; ct = ct->next) {
+					if (ct->tar)
+						subframe_updateObject(scene, ct->tar, 0, recursion, frame);
+				}
+				/* free temp targets */
+				if (cti->flush_constraint_targets)
+					cti->flush_constraint_targets(con, &targets, 0);
+			}
+		}
+	}
+
+	/* was originally OB_RECALC_ALL - TODO - which flags are really needed??? */
+	ob->recalc |= OB_RECALC_OB | OB_RECALC_DATA | OB_RECALC_TIME;
+	BKE_animsys_evaluate_animdata(scene, &ob->id, ob->adt, frame, ADT_RECALC_ANIM);
+	if (update_mesh) {
+		/* ignore cache clear during subframe updates
+		 *  to not mess up cache validity */
+		object_cacheIgnoreClear(ob, 1);
+		BKE_object_handle_update(scene, ob);
+		object_cacheIgnoreClear(ob, 0);
+	}
+	else
+		BKE_object_where_is_calc_time(scene, ob, frame);
+
+	/* for curve following objects, parented curve has to be updated too */
+	if (ob->type == OB_CURVE) {
+		Curve *cu = ob->data;
+		BKE_animsys_evaluate_animdata(scene, &cu->id, cu->adt, frame, ADT_RECALC_ANIM);
+	}
+	/* and armatures... */
+	if (ob->type == OB_ARMATURE) {
+		bArmature *arm = ob->data;
+		BKE_animsys_evaluate_animdata(scene, &arm->id, arm->adt, frame, ADT_RECALC_ANIM);
+		BKE_pose_where_is(scene, ob);
+	}
+
+	return 0;
+}
+
+/**********************************************************
+ *	Flow emission code
+ **********************************************************/
+
 typedef struct EmissionMap {
 	float *influence;
 	float *influence_high;
 	float *velocity;
 	int min[3], max[3], res[3];
+	int hmin[3], hmax[3], hres[3];
 	int total_cells, valid;
 } EmissionMap;
 
@@ -935,7 +1049,7 @@ static void clampBoundsInDomain(SmokeDomainSettings *sds, int min[3], int max[3]
 
 		/* adapt to velocity */
 		if (min_vel && min_vel[i] < 0.0f) {
-			min[i] += (int)ceil(min_vel[i] * dt);
+			min[i] += (int)floor(min_vel[i] * dt);
 		}
 		if (max_vel && max_vel[i] > 0.0f) {
 			max[i] += (int)ceil(max_vel[i] * dt);
@@ -967,8 +1081,16 @@ static void em_allocateData(EmissionMap *em, int use_velocity, int hires_mul)
 	/* allocate high resolution map if required */
 	if (hires_mul > 1) {
 		int total_cells_high = em->total_cells * (hires_mul * hires_mul * hires_mul);
+
+		for (i = 0; i < 3; i++) {
+			em->hmin[i] = em->min[i] * hires_mul;
+			em->hmax[i] = em->max[i] * hires_mul;
+			em->hres[i] = em->res[i] * hires_mul;
+		}
+
 		em->influence_high = MEM_callocN(sizeof(float) * total_cells_high, "smoke_flow_influence_high");
 	}
+	em->valid = 1;
 }
 
 static void em_freeData(EmissionMap *em)
@@ -981,8 +1103,108 @@ static void em_freeData(EmissionMap *em)
 		MEM_freeN(em->velocity);
 }
 
+static void em_combineMaps(EmissionMap *output, EmissionMap *em2, int hires_multiplier, int additive, float sample_size)
+{
+	int i, x,y,z;
+
+	/* copyfill input 1 struct and clear output for new allocation */
+	EmissionMap em1;
+	memcpy(&em1, output, sizeof(EmissionMap));
+	memset(output, 0, sizeof(EmissionMap));
+
+	for (i = 0; i < 3; i++) {
+		if (em1.valid) {
+			output->min[i] = MIN2(em1.min[i], em2->min[i]);
+			output->max[i] = MAX2(em1.max[i], em2->max[i]);
+		}
+		else {
+			output->min[i] = em2->min[i];
+			output->max[i] = em2->max[i];
+		}
+	}
+	/* allocate output map */
+	em_allocateData(output, (em1.velocity || em2->velocity), hires_multiplier);
+
+	/* base resolution inputs */
+	for (x = output->min[0]; x < output->max[0]; x++)
+		for (y = output->min[1]; y < output->max[1]; y++)
+			for (z = output->min[2]; z < output->max[2]; z++) {
+				int index_out = smoke_get_index(x - output->min[0], output->res[0], y - output->min[1], output->res[1], z - output->min[2]);
+
+				/* initialize with first input if in range */
+				if (x >= em1.min[0] && x < em1.max[0] &&
+					y >= em1.min[1] && y < em1.max[1] &&
+					z >= em1.min[2] && z < em1.max[2]) {
+					int index_in = smoke_get_index(x - em1.min[0], em1.res[0], y - em1.min[1], em1.res[1], z - em1.min[2]);
+
+					/* values */
+					output->influence[index_out] = em1.influence[index_in];
+					if (output->velocity) {
+						output->velocity[index_out] = em1.velocity[index_in];
+					}
+				}
+
+				/* apply second input if in range */
+				if (x >= em2->min[0] && x < em2->max[0] &&
+					y >= em2->min[1] && y < em2->max[1] &&
+					z >= em2->min[2] && z < em2->max[2]) {
+					int index_in = smoke_get_index(x - em2->min[0], em2->res[0], y - em2->min[1], em2->res[1], z - em2->min[2]);
+
+					/* values */
+					if (additive) {
+						output->influence[index_out] += em2->influence[index_in] * sample_size;
+					}
+					else {
+						output->influence[index_out] = MAX2(em2->influence[index_in], output->influence[index_out]);
+					}
+					if (output->velocity) {
+						/* last sample replaces the velocity */
+						output->velocity[index_out] = ADD_IF_LOWER(output->velocity[index_out], em2->velocity[index_in]);
+					}
+				}
+	} // low res loop
+
+
+
+	/* initialize high resolution input if available */
+	if (output->influence_high) {
+		for (x = output->hmin[0]; x < output->hmax[0]; x++)
+			for (y = output->hmin[1]; y < output->hmax[1]; y++)
+				for (z = output->hmin[2]; z < output->hmax[2]; z++) {
+					int index_out = smoke_get_index(x - output->hmin[0], output->hres[0], y - output->hmin[1], output->hres[1], z - output->hmin[2]);
+
+					/* initialize with first input if in range */
+					if (x >= em1.hmin[0] && x < em1.hmax[0] &&
+						y >= em1.hmin[1] && y < em1.hmax[1] &&
+						z >= em1.hmin[2] && z < em1.hmax[2]) {
+						int index_in = smoke_get_index(x - em1.hmin[0], em1.hres[0], y - em1.hmin[1], em1.hres[1], z - em1.hmin[2]);
+						/* values */
+						output->influence_high[index_out] = em1.influence_high[index_in];
+					}
+
+					/* apply second input if in range */
+					if (x >= em2->hmin[0] && x < em2->hmax[0] &&
+						y >= em2->hmin[1] && y < em2->hmax[1] &&
+						z >= em2->hmin[2] && z < em2->hmax[2]) {
+						int index_in = smoke_get_index(x - em2->hmin[0], em2->hres[0], y - em2->hmin[1], em2->hres[1], z - em2->hmin[2]);
+
+						/* values */
+						if (additive) {
+							output->influence_high[index_out] += em2->influence_high[index_in] * sample_size;
+						}
+						else {
+							output->influence_high[index_out] = MAX2(em2->influence_high[index_in], output->influence_high[index_out]);
+						}
+					}
+		} // high res loop
+	}
+
+	/* free original data */
+	em_freeData(&em1);
+}
+
 
-static void emit_from_particles(Object *flow_ob, SmokeDomainSettings *sds, SmokeFlowSettings *sfs, EmissionMap *em, Scene *scene, float time, float dt)
+static void emit_from_particles(Object *flow_ob, SmokeDomainSettings *sds, SmokeFlowSettings *sfs, EmissionMap *em, Scene *scene, float dt)
 {
 	if (sfs && sfs->psys && sfs->psys->part && ELEM(sfs->psys->part->type, PART_EMITTER, PART_FLUID)) // is particle system selected
 	{
@@ -993,23 +1215,44 @@ static void emit_from_particles(Object *flow_ob, SmokeDomainSettings *sds, Smoke
 		int totpart = psys->totpart, totchild;
 		int p = 0;
 		int valid_particles = 0;
+		int bounds_margin = 1;
+
+		/* radius based flow */
+		float solid = sfs->particle_size * 0.5f;
+		float smooth = 0.5f; /* add 0.5 cells of linear falloff to reduce aliasing */
+		int hires_multiplier = 1;
+		int i,z;
+		KDTree *tree;
 
 		sim.scene = scene;
 		sim.ob = flow_ob;
 		sim.psys = psys;
 		sim.rng = BLI_rng_new(psys->seed);
 
-		if (psys->part->type == PART_HAIR)
-		{
+		/* initialize particle cache */
+		if (psys->part->type == PART_HAIR) {
 			// TODO: PART_HAIR not supported whatsoever
 			totchild = 0;
 		}
-		else
+		else {
 			totchild = psys->totchild * psys->part->disp / 100;
+		}
 
 		particle_pos = MEM_callocN(sizeof(float) * (totpart + totchild) * 3, "smoke_flow_particles");
 		particle_vel = MEM_callocN(sizeof(float) * (totpart + totchild) * 3, "smoke_flow_particles");
 
+		/* setup particle radius emission if enabled */
+		if (sfs->flags & MOD_SMOKE_FLOW_USE_PART_SIZE) {
+			tree = BLI_kdtree_new(psys->totpart);
+
+			/* check need for high resolution map */
+			if ((sds->flags & MOD_SMOKE_HIGHRES) && (sds->highres_sampling == SM_HRES_FULLSAMPLE)) {
+				hires_multiplier = sds->amplify + 1;
+			}
+
+			bounds_margin = (int)ceil(solid + smooth);
+		}
+
 		/* calculate local position for each particle */
 		for (p = 0; p < totpart + totchild; p++)
 		{
@@ -1026,7 +1269,7 @@ static void emit_from_particles(Object *flow_ob, SmokeDomainSettings *sds, Smoke
 					continue;
 			}
 
-			state.time = time;
+			state.time = BKE_scene_frame_get(scene); /* use scene time */
 			if (psys_get_particle_state(&sim, p, &state, 0) == 0)
 				continue;
 
@@ -1039,45 +1282,121 @@ static void emit_from_particles(Object *flow_ob, SmokeDomainSettings *sds, Smoke
 			copy_v3_v3(&particle_vel[valid_particles * 3], state.vel);
 			mul_mat3_m4_v3(sds->imat, &particle_vel[valid_particles * 3]);
 
+			if (sfs->flags & MOD_SMOKE_FLOW_USE_PART_SIZE) {
+				BLI_kdtree_insert(tree, p, pos, NULL);
+			}
+
 			/* calculate emission map bounds */
 			em_boundInsert(em, pos);
 			valid_particles++;
 		}
 
 		/* set emission map */
-		clampBoundsInDomain(sds, em->min, em->max, NULL, NULL, 1, dt);
-		em_allocateData(em, sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY, 0);
+		clampBoundsInDomain(sds, em->min, em->max, NULL, NULL, bounds_margin, dt);
+		em_allocateData(em, sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY, hires_multiplier);
 
-		for (p = 0; p < valid_particles; p++)
+		if (!(sfs->flags & MOD_SMOKE_FLOW_USE_PART_SIZE)) {
+			for (p = 0; p < valid_particles; p++)
+			{
+				int cell[3];
+				size_t i = 0;
+				size_t index = 0;
+				int badcell = 0;
+
+				/* 1. get corresponding cell */
+				cell[0] = floor(particle_pos[p * 3]) - em->min[0];
+				cell[1] = floor(particle_pos[p * 3 + 1]) - em->min[1];
+				cell[2] = floor(particle_pos[p * 3 + 2]) - em->min[2];
+				/* check if cell is valid (in the domain boundary) */
+				for (i = 0; i < 3; i++) {
+					if ((cell[i] > em->res[i] - 1) || (cell[i] < 0)) {
+						badcell = 1;
+						break;
+					}
+				}
+				if (badcell)
+					continue;
+				/* get cell index */
+				index = smoke_get_index(cell[0], em->res[0], cell[1], em->res[1], cell[2]);
+				/* Add influence to emission map */
+				em->influence[index] = 1.0f;
+				/* Uses particle velocity as initial velocity for smoke */
+				if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY && (psys->part->phystype != PART_PHYS_NO))
+				{
+					VECADDFAC(&em->velocity[index * 3], &em->velocity[index * 3], &particle_vel[p * 3], sfs->vel_multi);
+				}
+			}   // particles loop
+		}
+		else // MOD_SMOKE_FLOW_USE_PART_SIZE
 		{
-			int cell[3];
-			size_t i = 0;
-			size_t index = 0;
-			int badcell = 0;
-
-			/* 1. get corresponding cell */
-			cell[0] = floor(particle_pos[p * 3]) - em->min[0];
-			cell[1] = floor(particle_pos[p * 3 + 1]) - em->min[1];
-			cell[2] = floor(particle_pos[p * 3 + 2]) - em->min[2];
-			/* check if cell is valid (in the domain boundary) */
+			int min[3], max[3], res[3];
+			float hr = 1.0f / ((float)hires_multiplier);
+			/* slightly adjust high res antialias smoothness based on number of divisions
+			 * to allow smaller details but yet not differing too much from the low res size */
+			float hr_smooth = smooth * pow(hr, 1.0f/3.0f);
+
+			/* setup loop bounds */
 			for (i = 0; i < 3; i++) {
-				if ((cell[i] > em->res[i] - 1) || (cell[i] < 0)) {
-					badcell = 1;
-					break;
-				}
+				min[i] = em->min[i] * hires_multiplier;
+				max[i] = em->max[i] * hires_multiplier;
+				res[i] = em->res[i] * hires_multiplier;
 			}
-			if (badcell)
-				continue;
-			/* get cell index */
-			index = smoke_get_index(cell[0], em->res[0], cell[1], em->res[1], cell[2]);
-			/* Add influence to emission map */
-			em->influence[index] = 1.0f;
-			/* Uses particle velocity as initial velocity for smoke */
-			if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY && (psys->part->phystype != PART_PHYS_NO))
-			{
-				VECADDFAC(&em->velocity[index * 3], &em->velocity[index * 3], &particle_vel[p * 3], sfs->vel_multi);
+
+			BLI_kdtree_balance(tree);
+
+			/* begin thread safe malloc */
+			BLI_begin_threaded_malloc();
+
+			#pragma omp parallel for schedule(static)
+			for (z = min[2]; z < max[2]; z++) {
+				int x, y;
+				for (x = min[0]; x < max[0]; x++)
+					for (y = min[1]; y < max[1]; y++) {
+						/* take low res samples where possible */
+						if (hires_multiplier <= 1 || !(x % hires_multiplier || y % hires_multiplier || z % hires_multiplier)) {
+							/* get low res space coordinates */
+							int lx = x / hires_multiplier;
+							int ly = y / hires_multiplier;
+							int lz = z / hires_multiplier;
+
+							int index = smoke_get_index(lx - em->min[0], em->res[0], ly - em->min[1], em->res[1], lz - em->min[2]);
+							float ray_start[3] = {((float)lx) + 0.5f, ((float)ly) + 0.5f, ((float)lz) + 0.5f};
+
+							/* find particle distance from the kdtree */
+							KDTreeNearest nearest;
+							float range = solid + smooth;
+							BLI_kdtree_find_nearest(tree, ray_start, NULL, &nearest);
+
+							if (nearest.dist < range) {
+								em->influence[index] = (nearest.dist < solid) ? 1.0f : (1.0f - (nearest.dist-solid) / smooth);
+							}
+						}
+
+						/* take high res samples if required */
+						if (hires_multiplier > 1) {
+							/* get low res space coordinates */
+							float lx = ((float)x) * hr;
+							float ly = ((float)y) * hr;
+							float lz = ((float)z) * hr;
+
+							int index = smoke_get_index(x - min[0], res[0], y - min[1], res[1], z - min[2]);
+							float ray_start[3] = {lx + 0.5f*hr, ly + 0.5f*hr, lz + 0.5f*hr};
+
+							/* find particle distance from the kdtree */
+							KDTreeNearest nearest;
+							float range = solid + hr_smooth;
+							BLI_kdtree_find_nearest(tree, ray_start, NULL, &nearest);
+
+							if (nearest.dist < range) {
+								em->influence_high[index] = (nearest.dist < solid) ? 1.0f : (1.0f - (nearest.dist-solid) / smooth);
+							}
+						}
+
+					}
 			}
-		}   // particles loop
+			BLI_end_threaded_malloc();
+			BLI_kdtree_free(tree);
+		}
 
 		/* free data */
 		if (particle_pos)
@@ -1108,7 +1427,7 @@ static void get_texture_value(Tex *texture, float tex_co[3], TexResult *texres)
 	}
 }
 
-static void sample_derived_mesh(SmokeFlowSettings *sfs, MVert *mvert, MTFace *tface, MFace *mface, float *influence_map, float *velocity_map, int index, int base_res[3], float flow_center[3], BVHTreeFromMesh *treeData, float ray_start[3],
+static void sample_derivedmesh(SmokeFlowSettings *sfs, MVert *mvert, MTFace *tface, MFace *mface, float *influence_map, float *velocity_map, int index, int base_res[3], float flow_center[3], BVHTreeFromMesh *treeData, float ray_start[3],
 								float *vert_vel, int has_velocity, int defgrp_index, MDeformVert *dvert, float x, float y, float z)
 {
 	float ray_dir[3] = {1.0f, 0.0f, 0.0f};
@@ -1336,7 +1655,7 @@ static void emit_from_derivedmesh(Object *flow_ob, SmokeDomainSettings *sds, Smo
 							int index = smoke_get_index(lx - em->min[0], em->res[0], ly - em->min[1], em->res[1], lz - em->min[2]);
 							float ray_start[3] = {((float)lx) + 0.5f, ((float)ly) + 0.5f, ((float)lz) + 0.5f};
 
-							sample_derived_mesh(sfs, mvert, tface, mface, em->influence, em->velocity, index, sds->base_res, flow_center, &treeData, ray_start,
+							sample_derivedmesh(sfs, mvert, tface, mface, em->influence, em->velocity, index, sds->base_res, flow_center, &treeData, ray_start,
 												vert_vel, has_velocity, defgrp_index, dvert, (float)lx, (float)ly, (float)lz);
 						}
 
@@ -1351,7 +1670,7 @@ static void emit_from_derivedmesh(Object *flow_ob, SmokeDomainSettings *sds, Smo
 							int index = smoke_get_index(x - min[0], res[0], y - min[1], res[1], z - min[2]);
 							float ray_start[3] = {lx + 0.5f*hr, ly + 0.5f*hr, lz + 0.5f*hr};
 
-							sample_derived_mesh(sfs, mvert, tface, mface, em->influence_high, NULL, index, sds->base_res, flow_center, &treeData, ray_start,
+							sample_derivedmesh(sfs, mvert, tface, mface, em->influence_high, NULL, index, sds->base_res, flow_center, &treeData, ray_start,
 												vert_vel, has_velocity, defgrp_index, dvert, lx, ly, lz); /* x,y,z needs to be always lowres */
 						}
 
@@ -1369,6 +1688,10 @@ static void emit_from_derivedmesh(Object *flow_ob, SmokeDomainSettings *sds, Smo
 	}
 }
 
+/**********************************************************
+ *	Smoke step
+ **********************************************************/
+
 static void adjustDomainResolution(SmokeDomainSettings *sds, int new_shift[3], EmissionMap *emaps, unsigned int numflowobj, float dt)
 {
 	int min[3] = {32767, 32767, 32767}, max[3] = {-32767, -32767, -32767}, res[3];
@@ -1658,8 +1981,8 @@ BLI_INLINE void apply_inflow_fields(SmokeFlowSettings *sfs, float emission_value
 	float dens_flow = (sfs->type == MOD_SMOKE_FLOW_TYPE_FIRE) ? 0.0f : emission_value * sfs->density;
 	float fuel_flow = emission_value * sfs->fuel_amount;
 	/* add heat */
-	if (heat) {
-		heat[index] = ADD_IF_LOWER(heat[index], emission_value * sfs->temp);
+	if (heat && emission_value > 0.0f) {
+		heat[index] = ADD_IF_LOWER(heat[index], sfs->temp);
 	}
 	/* absolute */
 	if (absolute_flow) {
@@ -1706,7 +2029,7 @@ BLI_INLINE void apply_inflow_fields(SmokeFlowSettings *sfs, float emission_value
 	}
 }
 
-static void update_flowsfluids(Scene *scene, Object *ob, SmokeDomainSettings *sds, float time, float dt)
+static void update_flowsfluids(Scene *scene, Object *ob, SmokeDomainSettings *sds, float dt)
 {
 	Object **flowobjs = NULL;
 	EmissionMap *emaps = NULL;
@@ -1764,13 +2087,61 @@ static void update_flowsfluids(Scene *scene, Object *ob, SmokeDomainSettings *sd
 		{
 			// we got nice flow object
 			SmokeFlowSettings *sfs = smd2->flow;
+			int subframes = sfs->subframes;
 			EmissionMap *em = &emaps[flowIndex];
 
-			if (sfs->source == MOD_SMOKE_FLOW_SOURCE_PARTICLES) {
-				emit_from_particles(collob, sds, sfs, em, scene, time, dt);
+			/* just sample flow directly to emission map if no subframes */
+			if (!subframes) {
+				if (sfs->source == MOD_SMOKE_FLOW_SOURCE_PARTICLES) {
+					emit_from_particles(collob, sds, sfs, em, scene, dt);
+				}
+				else {
+					emit_from_derivedmesh(collob, sds, sfs, em, dt);
+				}
 			}
+			/* sample subframes */
 			else {
-				emit_from_derivedmesh(collob, sds, sfs, em, dt);
+				int scene_frame = scene->r.cfra;
+				float scene_subframe = scene->r.subframe;
+				int subframe;
+				for (subframe = 0; subframe <= subframes; subframe++) {
+					EmissionMap em_temp = {0};
+					float sample_size = 1.0f / (float)(subframes+1);
+					float prev_frame_pos = sample_size * (float)(subframe+1);
+					float sdt = dt * sample_size;
+					int hires_multiplier = 1;
+
+					if ((sds->flags & MOD_SMOKE_HIGHRES) && (sds->highres_sampling == SM_HRES_FULLSAMPLE)) {
+						hires_multiplier = sds->amplify + 1;
+					}
+
+					/* set scene frame to match previous frame + subframe
+					 * or use current frame for last sample */
+					if (subframe < subframes) {
+						scene->r.cfra = scene_frame - 1;
+						scene->r.subframe = prev_frame_pos;
+					}
+					else {
+						scene->r.cfra = scene_frame;
+						scene->r.subframe = 0.0f;
+					}
+
+					if (sfs->source == MOD_SMOKE_FLOW_SOURCE_PARTICLES) {
+						/* emit_from_particles() updates timestep internally */
+						emit_from_particles(collob, sds, sfs, &em_temp, scene, sdt);
+					}
+					else { /* MOD_SMOKE_FLOW_SOURCE_MESH */
+						/* update flow object frame */
+						subframe_updateObject(scene, collob, 1, 5, BKE_scene_frame_get(scene));
+
+						/* apply flow */
+						emit_from_derivedmesh(collob, sds, sfs, &em_temp, sdt);
+					}
+
+					/* combine emission maps */
+					em_combineMaps(em, &em_temp, hires_multiplier, !(sfs->flags & MOD_SMOKE_FLOW_ABSOLUTE), sample_size);
+					em_freeData(&em_temp);
+				}
 			}
 
 			/* update required data fields */
@@ -2140,7 +2511,7 @@ static void step(Scene *scene, Object *ob, SmokeModifierData *smd, DerivedMesh *
 	for (substep = 0; substep < totalSubsteps; substep++)
 	{
 		// calc animated obstacle velocities
-		update_flowsfluids(scene, ob, sds, smd->time, dtSubdiv);
+		update_flowsfluids(scene, ob, sds, dtSubdiv);
 		update_obstacles(scene, ob, sds, dtSubdiv, substep, totalSubsteps);
 
 		if (sds->total_cells > 1) {
@@ -2336,16 +2707,6 @@ static void smokeModifier_process(SmokeModifierData *smd, Scene *scene, Object *
 
 		/* if on second frame, write cache for first frame */
 		if ((int)smd->time == startframe && (cache->flag & PTCACHE_OUTDATED || cache->last_exact == 0)) {
-			// create shadows straight after domain initialization so we get nice shadows for startframe, too
-			smoke_calc_transparency(sds, scene);
-
-			if (sds->wt && sds->total_cells > 1)
-			{
-				if (sds->flags & MOD_SMOKE_DISSOLVE)
-					smoke_dissolve_wavelet(sds->wt, sds->diss_speed, sds->flags & MOD_SMOKE_DISSOLVE_LOG);
-				smoke_turbulence_step(sds->wt, sds->fluid);
-			}
-
 			BKE_ptcache_write(&pid, startframe);
 		}
 
@@ -2358,8 +2719,15 @@ static void smokeModifier_process(SmokeModifierData *smd, Scene *scene, Object *
 		// DG: interesting commenting this line + deactivating loading of noise files
 		if (framenr != startframe)
 		{
-			if (sds->flags & MOD_SMOKE_DISSOLVE)
+			if (sds->flags & MOD_SMOKE_DISSOLVE) {
+				/* low res dissolve */
 				smoke_dissolve(sds->fluid, sds->diss_speed, sds->flags & MOD_SMOKE_DISSOLVE_LOG);
+				/* high res dissolve */
+				if (sds->wt) {
+					smoke_dissolve_wavelet(sds->wt, sds->diss_speed, sds->flags & MOD_SMOKE_DISSOLVE_LOG);
+				}
+
+			}
 
 			step(scene, ob, smd, dm, scene->r.frs_sec / scene->r.frs_sec_base);
 		}
@@ -2369,8 +2737,6 @@ static void smokeModifier_process(SmokeModifierData *smd, Scene *scene, Object *
 
 		if (sds->wt)
 		{
-			if (sds->flags & MOD_SMOKE_DISSOLVE)
-				smoke_dissolve_wavelet(sds->wt, sds->diss_speed, sds->flags & MOD_SMOKE_DISSOLVE_LOG);
 			smoke_turbulence_step(sds->wt, sds->fluid);
 		}