Merge from trunk r49908-r50219

- Add Smoke Density node
- Make density available on GPU
- *TOTALLY* untested! Will most likely crash and you cannot render at the moment since "Volume BSDF" is missing (from blenderartists volume patch)

TODO:
- Fix "svm_node_smoke_density()", atm I have no idea how

CMakeLists.txt

@@ -127,7 +127,7 @@ option(WITH_PYTHON_MODULE "Enable building as a python module which runs without
 option(WITH_BUILDINFO     "Include extra build details (only disable for development & faster builds)" ON)
 option(WITH_IK_ITASC      "Enable ITASC IK solver (only disable for development & for incompatible C++ compilers)" ON)
 option(WITH_IK_SOLVER     "Enable Legacy IK solver (only disable for development)" ON)
-option(WITH_FFTW3         "Enable FFTW3 support (Used for smoke and audio effects)" OFF)
+option(WITH_FFTW3         "Enable FFTW3 support (Used for smoke and audio effects)" ON)
 option(WITH_BULLET        "Enable Bullet (Physics Engine)" ON)
 option(WITH_GAMEENGINE    "Enable Game Engine" ON)
 option(WITH_PLAYER        "Build Player" OFF)
@@ -153,7 +153,7 @@ mark_as_advanced(WITH_AUDASPACE)
 if((UNIX AND NOT APPLE) OR (MINGW))
 	set(PLATFORM_DEFAULT ON)
 else()
-	set(PLATFORM_DEFAULT OFF)
+	set(PLATFORM_DEFAULT ON)
 endif()
 option(WITH_OPENMP        "Enable OpenMP (has to be supported by the compiler)" ${PLATFORM_DEFAULT})
 unset(PLATFORM_DEFAULT)
@@ -1685,9 +1685,9 @@ if(WITH_PYTHON)
 		if(NOT ${PYTHON_NUMPY_PATH} STREQUAL "")
 			if(NOT EXISTS "${PYTHON_NUMPY_PATH}/numpy")
 				message(WARNING "PYTHON_NUMPY_PATH is invalid, numpy not found in '${PYTHON_NUMPY_PATH}' "
-				                "WITH_PYTHON_INSTALL_NUMPY option will be ignored when installing python")
+							"WITH_PYTHON_INSTALL_NUMPY option will be ignored when installing python")
 				set(WITH_PYTHON_INSTALL_NUMPY OFF)
-			endif()
+		endif()
 		# not set, so initialize
 		else()
 			string(REPLACE "." ";" _PY_VER_SPLIT "${PYTHON_VERSION}")
@@ -1717,11 +1717,11 @@ if(WITH_PYTHON)
 				set(WITH_PYTHON_INSTALL_NUMPY OFF)
 			else()
 				message(STATUS "numpy found at '${PYTHON_NUMPY_PATH}'")
-			endif()
+	endif()
 
 			unset(_PY_VER_SPLIT)
 			unset(_PY_VER_MAJOR)
-		endif()
+endif()
 	endif()
 endif()
 

intern/cycles/blender/CMakeLists.txt

@@ -23,6 +23,7 @@ set(SRC
 	blender_mesh.cpp
 	blender_object.cpp
 	blender_particles.cpp
+	blender_smoke.cpp
 	blender_python.cpp
 	blender_session.cpp
 	blender_shader.cpp

intern/cycles/blender/blender_object.cpp

@@ -281,6 +281,9 @@ void BlenderSync::sync_object(BL::Object b_parent, int b_index, BL::Object b_ob,
 		if (need_particle_update)
 			sync_particles(object, b_ob);
 
+		if(object_use_smoke(b_ob))
+			sync_smoke(object, b_ob);
+	
 		object->tag_update(scene);
 	}
 }

intern/cycles/blender/blender_smoke.cpp

@@ -0,0 +1,96 @@
+/*
+ * Copyright 2011, Blender Foundation.
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include "object.h"
+
+#include "mesh.h"
+#include "blender_sync.h"
+#include "blender_util.h"
+
+#include "util_foreach.h"
+
+CCL_NAMESPACE_BEGIN
+
+/* Utilities */
+
+
+/* Smoke Sync */
+
+/* Only looking for Smoke domains */
+// TODO DG: disable rendering of smoke flow??
+bool BlenderSync::object_use_smoke(BL::Object b_ob)
+{
+	BL::Object::modifiers_iterator b_modifiers;
+	for(b_ob.modifiers.begin(b_modifiers); b_modifiers != b_ob.modifiers.end(); ++b_modifiers) {
+		BL::Modifier mod = (*b_modifiers);
+
+		if (mod.is_a(&RNA_SmokeModifier)) {
+			BL::SmokeModifier smd(mod);
+
+			if(smd.smoke_type() == BL::SmokeModifier::smoke_type_DOMAIN) {
+				BL::ID key = (BKE_object_is_modified(b_ob))? b_ob: b_ob.data();
+				Mesh *mesh = mesh_map.find(key);
+				if (mesh) {
+					return mesh->need_attribute(scene, ATTR_STD_SMOKE_DENSITY);
+				}
+			}
+		}
+	}
+
+	return false;
+}
+
+static BL::SmokeModifier *get_smoke(BL::Object b_ob)
+{
+	BL::Object::modifiers_iterator b_modifiers;
+	for(b_ob.modifiers.begin(b_modifiers); b_modifiers != b_ob.modifiers.end(); ++b_modifiers) {
+		BL::Modifier mod = (*b_modifiers);
+
+		if (mod.is_a(&RNA_SmokeModifier)) {
+			BL::SmokeModifier *smd = (BL::SmokeModifier *)(&mod);
+
+			if(smd->smoke_type() == BL::SmokeModifier::smoke_type_DOMAIN) {
+				return smd;
+			}
+		}
+	}
+	
+	return NULL;
+}
+
+void BlenderSync::sync_smoke(Object *ob, BL::Object b_ob)
+{
+	BL::SmokeModifier *smd = get_smoke(b_ob);
+	BL::SmokeDomainSettings sds = smd->domain_settings();
+
+	ob->grid.clear();
+	ob->resolution = get_int3(sds.domain_resolution());
+
+	int length[3];
+	int numcells = rna_SmokeModifier_density_get_length(&sds.ptr, length);
+
+	if(numcells != 0)
+	{
+		vector<float> &grid = ob->grid;
+		grid.reserve(numcells);
+		grid.resize(numcells);
+		rna_SmokeModifier_density_get(&sds.ptr, &grid[0]);
+	}
+}
+
+CCL_NAMESPACE_END

intern/cycles/blender/blender_sync.h

@@ -86,10 +86,12 @@ private:
 	void sync_mesh_motion(BL::Object b_ob, Mesh *mesh, int motion);
 	void sync_camera_motion(BL::Object b_ob, int motion);
 	void sync_particles(Object *ob, BL::Object b_ob);
+	void sync_smoke(Object *ob, BL::Object b_ob);
 
 	/* util */
 	void find_shader(BL::ID id, vector<uint>& used_shaders, int default_shader);
 	bool BKE_object_is_modified(BL::Object b_ob);
+	bool object_use_smoke(BL::Object b_ob);
 	bool object_is_mesh(BL::Object b_ob);
 	bool object_is_light(BL::Object b_ob);
 	bool object_need_particle_update(BL::Object b_ob);

intern/cycles/blender/blender_util.h

@@ -56,6 +56,9 @@ void rna_Scene_frame_set(void *scene, int frame, float subframe);
 void BKE_image_user_frame_calc(void *iuser, int cfra, int fieldnr);
 void BKE_image_user_file_path(void *iuser, void *ima, char *path);
 
+int rna_SmokeModifier_density_get_length(PointerRNA *ptr, int length[3]);
+void rna_SmokeModifier_density_get(PointerRNA *ptr, float *values);
+
 }
 
 CCL_NAMESPACE_BEGIN
@@ -155,6 +158,11 @@ static inline int4 get_int4(BL::Array<int, 4> array)
 	return make_int4(array[0], array[1], array[2], array[3]);
 }
 
+static inline int3 get_int3(BL::Array<int, 3> array)
+{
+	return make_int3(array[0], array[1], array[2]);
+}
+
 static inline uint get_layer(BL::Array<int, 20> array)
 {
 	uint layer = 0;

intern/cycles/kernel/kernel_textures.h

@@ -55,6 +55,9 @@ KERNEL_TEX(float2, texture_float2, __light_background_conditional_cdf)
 /* particles */
 KERNEL_TEX(float4, texture_float4, __particles)
 
+/* Smoke */
+KERNEL_TEX(float, texture_float, __smoke_density)
+
 /* shaders */
 KERNEL_TEX(uint4, texture_uint4, __svm_nodes)
 KERNEL_TEX(uint, texture_uint, __shader_flag)

intern/cycles/kernel/kernel_types.h

@@ -363,6 +363,7 @@ typedef enum AttributeStandard {
 	ATTR_STD_MOTION_PRE,
 	ATTR_STD_MOTION_POST,
 	ATTR_STD_PARTICLE,
+	ATTR_STD_SMOKE_DENSITY,
 	ATTR_STD_NUM,
 
 	ATTR_STD_NOT_FOUND = ~0

intern/cycles/kernel/svm/svm.h

@@ -272,6 +272,9 @@ __device_noinline void svm_eval_nodes(KernelGlobals *kg, ShaderData *sd, ShaderT
 			case NODE_PARTICLE_INFO:
 				svm_node_particle_info(kg, sd, stack, node.y, node.z);
 				break;
+			case NODE_SMOKE_DENSITY:
+				svm_node_smoke_density(kg, sd, stack, node.y, node.z);
+				break;
 #endif
 			case NODE_CONVERT:
 				svm_node_convert(sd, stack, node.y, node.z, node.w);

intern/cycles/kernel/svm/svm_geometry.h

@@ -122,5 +122,22 @@ __device void svm_node_particle_info(KernelGlobals *kg, ShaderData *sd, float *s
 	}
 }
 
+
+/* Smoke Density */
+
+__device void svm_node_smoke_density(KernelGlobals *kg, ShaderData *sd, float *stack, uint type, uint out_offset)
+{
+	float data;
+
+	switch(type) {
+		case NODE_INFO_SMO_DEN: {
+			// DG TODO uint particle_id = object_particle_id(kg, sd->object);
+			// DG TODO data = smoke_density(kg, WHICH CELL IN GRID?);
+			// DG TODO stack_store_float(stack, out_offset, data);
+			break;
+		}
+	}
+}
+
 CCL_NAMESPACE_END
 

intern/cycles/kernel/svm/svm_types.h

@@ -89,7 +89,8 @@ typedef enum NodeType {
 	NODE_MIN_MAX,
 	NODE_LIGHT_FALLOFF,
 	NODE_OBJECT_INFO,
-	NODE_PARTICLE_INFO
+	NODE_PARTICLE_INFO,
+	NODE_SMOKE_DENSITY
 } NodeType;
 
 typedef enum NodeAttributeType {
@@ -119,6 +120,10 @@ typedef enum NodeParticleInfo {
 	NODE_INFO_PAR_LIFETIME
 } NodeParticleInfo;
 
+typedef enum NodeSmokeDensity {
+	NODE_INFO_SMO_DEN
+} NodeSmokeDensity;
+
 typedef enum NodeLightPath {
 	NODE_LP_camera = 0,
 	NODE_LP_shadow,

intern/cycles/render/nodes.cpp

@@ -1843,6 +1843,38 @@ void ParticleInfoNode::compile(OSLCompiler& compiler)
 	compiler.add(this, "node_particle_info");
 }
 
+/* Smoke Density */
+
+SmokeDensityNode::SmokeDensityNode()
+: ShaderNode("particle_info")
+{
+	add_output("Density", SHADER_SOCKET_FLOAT);
+}
+
+void SmokeDensityNode::attributes(AttributeRequestSet *attributes)
+{
+	if(!output("Density")->links.empty())
+		attributes->add(ATTR_STD_SMOKE_DENSITY);
+
+	ShaderNode::attributes(attributes);
+}
+
+void SmokeDensityNode::compile(SVMCompiler& compiler)
+{
+	ShaderOutput *out;
+	
+	out = output("Density");
+	if(!out->links.empty()) {
+		compiler.stack_assign(out);
+		compiler.add_node(NODE_SMOKE_DENSITY, NODE_INFO_SMO_DEN, out->stack_offset);
+	}
+}
+
+void SmokeDensityNode::compile(OSLCompiler& compiler)
+{
+	compiler.add(this, "node_smoke_density");
+}
+
 /* Value */
 
 ValueNode::ValueNode()

intern/cycles/render/nodes.h

@@ -296,6 +296,12 @@ public:
 	void attributes(AttributeRequestSet *attributes);
 };
 
+class SmokeDensityNode : public ShaderNode {
+public:
+	SHADER_NODE_CLASS(SmokeDensityNode)
+	void attributes(AttributeRequestSet *attributes);
+};
+
 class ValueNode : public ShaderNode {
 public:
 	SHADER_NODE_CLASS(ValueNode)

intern/cycles/render/object.cpp

@@ -167,7 +167,7 @@ void ObjectManager::device_update_transforms(Device *device, DeviceScene *dscene
 		float uniform_scale;
 		float surface_area = 0.0f;
 		float pass_id = ob->pass_id;
-		float random_number = (float)ob->random_id * (1.0f/(float)0xFFFFFFFF);
+		float random_number = (float)ob->random_id * (1.0f/(float)0xFFFFFFFF)+0.5f;
 		
 		if(transform_uniform_scale(tfm, uniform_scale)) {
 			map<Mesh*, float>::iterator it = surface_area_map.find(mesh);
@@ -281,6 +281,37 @@ void ObjectManager::device_update_particles(Device *device, DeviceScene *dscene,
 	device->tex_alloc("__particles", dscene->particles);
 }
 
+void ObjectManager::device_update_smoke(Device *device, DeviceScene *dscene, Scene *scene, Progress& progress)
+{
+	/* count smoke cells.
+	 * adds one dummy particle at the beginning to avoid invalid lookups,
+	 * in case a shader uses particle info without actual particle data.
+	 */
+	int num_cells = 1;
+	foreach(Object *ob, scene->objects)
+		num_cells += ob->grid.size();
+	
+	float *density = dscene->smoke_density.resize(num_cells);
+	
+	/* dummy particle */
+	// DG TODO density[0] = 0.0f;
+	
+	int i = 0;
+	foreach(Object *ob, scene->objects) {
+		/* pack in texture */
+		for(i = 0; i < ob->grid.size(); i++) {
+
+			// DG TODO: use "*PARTICLE_SIZE"?
+			density[i] = ob->grid[i];
+			
+			if(progress.get_cancel()) return;
+
+		}
+	}
+	
+	device->tex_alloc("__smoke_density", dscene->smoke_density);
+}
+
 void ObjectManager::device_update(Device *device, DeviceScene *dscene, Scene *scene, Progress& progress)
 {
 	if(!need_update)
@@ -311,6 +342,11 @@ void ObjectManager::device_update(Device *device, DeviceScene *dscene, Scene *sc
 	
 	if(progress.get_cancel()) return;
 	
+	progress.set_status("Updating Objects", "Copying Smoke Density to device");
+	device_update_smoke(device, dscene, scene, progress);
+	
+	if(progress.get_cancel()) return;
+
 	need_update = false;
 }
 

intern/cycles/render/object.h

@@ -58,6 +58,10 @@ public:
 	int particle_id;
 	vector<Particle> particles;
 
+	/* Voxel / 3D volume data */
+	int3 resolution;
+	vector<float> grid;
+
 	Object();
 	~Object();
 
@@ -79,6 +83,7 @@ public:
 	void device_update(Device *device, DeviceScene *dscene, Scene *scene, Progress& progress);
 	void device_update_transforms(Device *device, DeviceScene *dscene, Scene *scene, Progress& progress);
 	void device_update_particles(Device *device, DeviceScene *dscene, Scene *scene, Progress& progress);
+	void device_update_smoke(Device *device, DeviceScene *dscene, Scene *scene, Progress& progress);
 	void device_free(Device *device, DeviceScene *dscene);
 
 	void tag_update(Scene *scene);

intern/cycles/render/scene.h

@@ -85,6 +85,9 @@ public:
 	/* particles */
 	device_vector<float4> particles;
 
+	/* smoke */
+	device_vector<float> smoke_density;
+
 	/* shaders */
 	device_vector<uint4> svm_nodes;
 	device_vector<uint> shader_flag;

intern/smoke/CMakeLists.txt

@@ -25,8 +25,7 @@
 
 set(INC
 	intern
-	../memutil
+	../../extern/Eigen3
-	../../extern/bullet2/src
 )
 
 set(INC_SYS

intern/smoke/SConscript

@@ -9,11 +9,10 @@ defs = ''
 if env['WITH_BF_OPENMP']:
     if env['OURPLATFORM'] == 'linuxcross':
         incs += ' ' + env['BF_OPENMP_INC']
-
+        defs += ' PARALLEL=1'
-    defs += ' PARALLEL=1'
 
 incs += ' ' + env['BF_PNG_INC'] + ' ' + env['BF_ZLIB_INC'] 
-incs += ' intern ../../extern/bullet2/src ../memutil ../guardealloc '
+incs += ' intern ../memutil ../guardealloc ../../extern/Eigen3 '
 
 if env['WITH_BF_FFTW3']:
     defs += ' WITH_FFTW3'

intern/smoke/extern/smoke_API.h

@@ -53,6 +53,8 @@ float *smoke_get_velocity_x(struct FLUID_3D *fluid);
 float *smoke_get_velocity_y(struct FLUID_3D *fluid);
 float *smoke_get_velocity_z(struct FLUID_3D *fluid);
 
+float *smoke_get_pressure(struct FLUID_3D *fluid);
+
 /* Moving obstacle velocity provided by blender */
 void smoke_get_ob_velocity(struct FLUID_3D *fluid, float **x, float **y, float **z);
 

intern/smoke/intern/FLUID_3D.cpp

@@ -36,10 +36,16 @@
 
 #include "float.h"
 
+#include <iostream>
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+
 #if PARALLEL==1
 #include <omp.h>
 #endif // PARALLEL 
 
+
 //////////////////////////////////////////////////////////////////////
 // Construction/Destruction
 //////////////////////////////////////////////////////////////////////
@@ -98,6 +104,9 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dtdef) :
 	_heatOld      = new float[_totalCells];
 	_obstacles    = new unsigned char[_totalCells]; // set 0 at end of step
 
+	// For debugging purposes
+	_pressure	  = new float[_totalCells];
+
 	// For threaded version:
 	_xVelocityTemp = new float[_totalCells];
 	_yVelocityTemp = new float[_totalCells];
@@ -126,13 +135,14 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dtdef) :
 		_yForce[x]       = 0.0f;
 		_zForce[x]       = 0.0f;
 		_obstacles[x]    = false;
+		_pressure[x]	 = 0.0f;
 	}
 
 	// boundary conditions of the fluid domain
 	// set default values -> vertically non-colliding
-	_domainBcFront = true;
+	_domainBcFront = false;
 	_domainBcTop = false;
-	_domainBcLeft = true;
+	_domainBcLeft = false;
 	_domainBcBack = _domainBcFront;
 	_domainBcBottom = _domainBcTop;
 	_domainBcRight	= _domainBcLeft;
@@ -206,6 +216,10 @@ FLUID_3D::~FLUID_3D()
 	if (_obstacles) delete[] _obstacles;
     // if (_wTurbulence) delete _wTurbulence;
 
+	// for debugging purposes
+	if(_pressure)
+		delete[] _pressure;
+
 	if (_xVelocityTemp) delete[] _xVelocityTemp;
 	if (_yVelocityTemp) delete[] _yVelocityTemp;
 	if (_zVelocityTemp) delete[] _zVelocityTemp;
@@ -230,16 +244,6 @@ void FLUID_3D::initBlenderRNA(float *alpha, float *beta, float *dt_factor, float
 //////////////////////////////////////////////////////////////////////
 void FLUID_3D::step(float dt)
 {
-#if 0
-	// If border rules have been changed
-	if (_colloPrev != *_borderColli) {
-		printf("Border collisions changed\n");
-		
-		// DG TODO: Need to check that no animated obstacle flags are overwritten
-		setBorderCollisions();
-	}
-#endif
-
 	// DG: TODO for the moment redo border for every timestep since it's been deleted every time by moving obstacles
 	setBorderCollisions();
 
@@ -407,6 +411,17 @@ void FLUID_3D::step(float dt)
 	for (int i = 0; i < _totalCells; i++)
 	{
 		_xForce[i] = _yForce[i] = _zForce[i] = 0.0f;
+
+		/*
+		if(_density[i] < FLT_EPSILON)
+		{
+			// _density[i] = 0.0f;
+
+			_xVelocity[i] =
+			_yVelocity[i] = 
+			_zVelocity[i] = 0.0f;
+		}
+		*/
 	}
 
 }
@@ -739,6 +754,7 @@ void FLUID_3D::wipeBoundariesSL(int zBegin, int zEnd)
 	}
 
 }
+
 //////////////////////////////////////////////////////////////////////
 // add forces to velocity field
 //////////////////////////////////////////////////////////////////////
@@ -754,6 +770,7 @@ void FLUID_3D::addForce(int zBegin, int zEnd)
 		_zVelocityTemp[i] = _zVelocity[i] + _dt * _zForce[i];
 	}
 }
+
 //////////////////////////////////////////////////////////////////////
 // project into divergence free field
 //////////////////////////////////////////////////////////////////////
@@ -761,14 +778,46 @@ void FLUID_3D::project()
 {
 	int x, y, z;
 	size_t index;
+#if USE_NEW_CG == 1
+	VectorXf p(_totalCells);
 
-	float *_pressure = new float[_totalCells];
+	ArrayXd A0(_totalCells);
+	ArrayXd Ai(_totalCells);
+	ArrayXd Aj(_totalCells);
+	ArrayXd Ak(_totalCells);
+
+	A0.setZero(_totalCells);
+	Ai.setZero(_totalCells);
+	Aj.setZero(_totalCells);
+	Ak.setZero(_totalCells);
+	p.setZero(_totalCells);
+
+
+	ArrayXd gti(_totalCells); // gridToIndex
+	unsigned int linearIndex = 0;
+
+	for (z = 0; z < _zRes; z++)
+		for (y = 0; y < _yRes; y++)
+			for (x = 0; x < _xRes; x++)
+				if(!_obstacles[FINDEX(x, y, z)])
+				{
+					gti[FINDEX(x, y, z)] = linearIndex;
+					linearIndex++;
+				}
+
+
+	VectorXf b(linearIndex);
+	b.setZero(linearIndex);
+	SparseMatrix<float,RowMajor> A(linearIndex, linearIndex);
+	A.reserve(VectorXi::Constant(linearIndex, 7));
+#endif
+	// float *_pressure = new float[_totalCells];
 	float *_divergence   = new float[_totalCells];
 
 	memset(_pressure, 0, sizeof(float)*_totalCells);
 	memset(_divergence, 0, sizeof(float)*_totalCells);
 	
-	// set velocity and pressure inside of obstacles to zero
+	// set VELOCITY and PRESSURE inside of obstacles cells to zero
 	setObstacleBoundaries(_pressure, 0, _zRes);
 	
 	// copy out the boundaries
@@ -781,35 +830,6 @@ void FLUID_3D::project()
 	if(!_domainBcTop) setNeumannZ(_zVelocity, _res, 0, _zRes);
 	else setZeroZ(_zVelocity, _res, 0, _zRes);
 
-	/*
-	{
-		float maxx = 0, maxy = 0, maxz = 0;
-		for(unsigned int i = 0; i < _xRes * _yRes * _zRes; i++)
-		{
-			if(_xVelocity[i] > maxx)
-				maxx = _xVelocity[i];
-			if(_yVelocity[i] > maxy)
-				maxy = _yVelocity[i];
-			if(_zVelocity[i] > maxz)
-				maxz = _zVelocity[i];
-		}
-		printf("Max velx: %f, vely: %f, velz: %f\n", maxx, maxy, maxz);
-	}
-	*/
-
-	/*
-	{
-		float maxvalue = 0;
-		for(unsigned int i = 0; i < _xRes * _yRes * _zRes; i++)
-		{
-			if(_heat[i] > maxvalue)
-				maxvalue = _heat[i];
-
-		}
-		printf("Max heat: %f\n", maxvalue);
-	}
-	*/
-
 	// calculate divergence
 	index = _slabSize + _xRes + 1;
 	for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
@@ -823,7 +843,6 @@ void FLUID_3D::project()
 					continue;
 				}
 				
-
 				float xright = _xVelocity[index + 1];
 				float xleft  = _xVelocity[index - 1];
 				float yup    = _yVelocity[index + _xRes];
@@ -838,12 +857,12 @@ void FLUID_3D::project()
 				if(_obstacles[index+_slabSize]) ztop    = - _zVelocity[index];
 				if(_obstacles[index-_slabSize]) zbottom = - _zVelocity[index];
 
-				if(_obstacles[index+1] & 8)			xright	+= _xVelocityOb[index + 1];
+				if(_obstacles[index+1])			xright	+= _xVelocityOb[index + 1];
-				if(_obstacles[index-1] & 8)			xleft	+= _xVelocityOb[index - 1];
+				if(_obstacles[index-1])			xleft	+= _xVelocityOb[index - 1];
-				if(_obstacles[index+_xRes] & 8)		yup		+= _yVelocityOb[index + _xRes];
+				if(_obstacles[index+_xRes])		yup		+= _yVelocityOb[index + _xRes];
-				if(_obstacles[index-_xRes] & 8)		ydown	+= _yVelocityOb[index - _xRes];
+				if(_obstacles[index-_xRes])		ydown	+= _yVelocityOb[index - _xRes];
-				if(_obstacles[index+_slabSize] & 8) ztop    += _zVelocityOb[index + _slabSize];
+				if(_obstacles[index+_slabSize]) ztop    += _zVelocityOb[index + _slabSize];
-				if(_obstacles[index-_slabSize] & 8) zbottom += _zVelocityOb[index - _slabSize];
+				if(_obstacles[index-_slabSize]) zbottom += _zVelocityOb[index - _slabSize];
 
 				_divergence[index] = -_dx * 0.5f * (
 						xright - xleft +
@@ -854,37 +873,162 @@ void FLUID_3D::project()
 				_pressure[index] = 0.0f;
 			}
 
-
+	float scale = 1.0; // DG TODO: make this global and incooperate this into other functions
-	/*
+#if USE_NEW_CG == 1
+	for (z = 0; z < _zRes; z++)
+		for (y = 0; y < _yRes; y++)
+			for (x = 0; x < _xRes; x++)
 	{
-		float maxvalue = 0;
+				if(!_obstacles[FINDEX(x, y, z)])
-		for(unsigned int i = 0; i < _xRes * _yRes * _zRes; i++)
+				{
-		{
+					// +x
-			if(_divergence[i] > maxvalue)
+					if(x < _xRes -1)
-				maxvalue = _divergence[i];
+						{ if(!_obstacles[FINDEX(x + 1, y, z)]) A0[FINDEX(x, y, z)] += scale; }
+					else
+						A0[FINDEX(x, y, z)] += scale;
+
+					// -x
+					if(x > 0)
+						{ if(!_obstacles[FINDEX(x - 1, y, z)]) A0[FINDEX(x, y, z)] += scale; }
+					else
+						A0[FINDEX(x, y, z)] += scale;
+
+					// +y
+					if(y < _yRes - 1)
+						{ if(!_obstacles[FINDEX(x, y + 1, z)]) A0[FINDEX(x, y, z)] += scale; }
+					else
+						A0[FINDEX(x, y, z)] += scale;
+
+					// -y
+					if(y > 0)
+						{ if(!_obstacles[FINDEX(x, y - 1, z)]) A0[FINDEX(x, y, z)] += scale; }
+					else
+						A0[FINDEX(x, y, z)] += scale;
+
+					// +z
+					if(z < _zRes - 1)
+						{ if(!_obstacles[FINDEX(x, y, z + 1)]) A0[FINDEX(x, y, z)] += scale; }
+					else
+						A0[FINDEX(x, y, z)] += scale;
+
+					// -z
+					if(z > 0)
+						{ if(!_obstacles[FINDEX(x, y, z - 1)]) A0[FINDEX(x, y, z)] += scale; }
+					else
+						A0[FINDEX(x, y, z)] += scale;
+				}
 
-		}
-		printf("Max divergence: %f\n", maxvalue);
 	}
-	*/
 
-	copyBorderAll(_pressure, 0, _zRes);
+	for (z = 0; z < _zRes; z++)
-
+		for (y = 0; y < _yRes; y++)
-	/*
+			for (x = 0; x < _xRes; x++)
 	{
-		float maxvalue = 0;
+				if(!_obstacles[FINDEX(x, y, z)])
-		for(unsigned int i = 0; i < _xRes * _yRes * _zRes; i++)
+				{
-		{
+					// x
-			if(_pressure[i] > maxvalue)
+					if(x < _xRes - 1)
-				maxvalue = _pressure[i];
+						{ if(!_obstacles[FINDEX(x + 1, y, z)]) Ai[FINDEX(x, y, z)] = -scale; }
-		}
+					else
-		printf("Max pressure BEFORE: %f\n", maxvalue);
+						Ai[FINDEX(x, y, z)] = -scale;
+
+					// y
+					if(y < _yRes - 1)
+						{ if(!_obstacles[FINDEX(x, y + 1, z)]) Aj[FINDEX(x, y, z)] = -scale; }
+					else
+						Aj[FINDEX(x, y, z)] = -scale;
+
+					// z
+					if(z < _zRes - 1)
+						{ if(!_obstacles[FINDEX(x, y, z + 1)]) Ak[FINDEX(x, y, z)] = -scale; }
+					else
+						Ak[FINDEX(x, y, z)] = -scale;
+				}
+
 	}
-	*/
+
+	unsigned int rowCount = 0;
+	for (z = 0; z < _zRes; z++)
+		for (y = 0; y < _yRes; y++)
+			for (x = 0; x < _xRes; x++)
+		if (!_obstacles[FINDEX(x,y,z)])
+		{
+			rowCount = gti(FINDEX(x, y, z));
+
+			/*
+			if((x > 0) && (Ai[FINDEX(x - 1, y, z)] < 0))
+				A.insert(rowCount, gti(FINDEX(x - 1, y, z))) = Ai[FINDEX(x - 1, y, z)];
+			if((y > 0) && (Aj[FINDEX(x, y - 1, z)] < 0))
+				A.insert(rowCount, gti(FINDEX(x, y - 1, z))) = Aj[FINDEX(x, y - 1, z)];
+			if((z > 0) && (Ak[FINDEX(x, y, z - 1)] < 0))
+				A.insert(rowCount, gti(FINDEX(x, y, z - 1))) = Ak[FINDEX(x, y, z - 1)];
+			*/
+		
+			if(A0[FINDEX(x, y, z)] > 0)
+				A.insert(rowCount, gti(FINDEX(x, y, z))) = A0[FINDEX(x, y, z)];
+
+			if((x + 1 < _xRes) && (Ai[FINDEX(x, y, z)] < 0))
+				A.insert(rowCount, gti(FINDEX(x + 1, y, z))) = Ai[FINDEX(x, y, z)];
+			if((y + 1 < _yRes) && (Aj[FINDEX(x, y, z)] < 0))
+				A.insert(rowCount, gti(FINDEX(x, y + 1, z))) = Aj[FINDEX(x, y, z)];
+			if((z + 1 < _zRes) && (Ak[FINDEX(x, y, z)] < 0))
+				A.insert(rowCount, gti(FINDEX(x, y, z + 1))) = Ak[FINDEX(x, y, z)];
+
+			// rowCount++;
+		}
+
+	for (z = 0; z < _zRes; z++)
+		for (y = 0; y < _yRes; y++)
+			for (x = 0; x < _xRes; x++)
+		if (!_obstacles[FINDEX(x,y,z)])
+			b[gti(FINDEX(x, y, z))] = _divergence[FINDEX(x,y,z)];
+#endif
+	// copyBorderAll(_pressure, 0, _zRes);
 
 	// solve Poisson equation
-	solvePressurePre(_pressure, _divergence, _obstacles);
+#if USE_NEW_CG == 1
+	A.makeCompressed();
 
+	VectorXf result(linearIndex);
+	result.setZero(linearIndex);
+	ConjugateGradient< SparseMatrix<float, RowMajor>, Upper > solver;
+	solver.setMaxIterations(300);
+	solver.setTolerance(1e-06);
+
+	solver.compute(A);
+	if(solver.info() != Success) 
+	{
+	  // decomposition failed
+		printf("Solver error: Cannot create decomposition: %d\n", solver.info());
+	}
+	else
+	{
+		result = solver.solve(b);
+		if(solver.info() != Success) 
+		{
+		  // solving failed
+			printf("Solver error: Cannot solve equation: %d\n", solver.info());
+		}
+		else
+			printf("Solver finished.\n");
+
+		std::cout << "#iterations:     " << solver.iterations() << std::endl;
+		std::cout << "estimated error: " << solver.error()      << std::endl;
+	}
+#else
+	solvePressurePre(_pressure, _divergence, _obstacles);
+#endif
+
+#if 0
+	for(unsigned int i = 0; i < _xRes * _yRes * _zRes; i++)
+	{
+		float value = (_pressure[i] - result[i]);
+		if(value > 0.01)
+			printf("error: p: %f, b: %f\n", _pressure[i], result[i]);
+	}
+#endif
+
+#if 1
 	{
 		float maxvalue = 0;
 		for(unsigned int i = 0; i < _xRes * _yRes * _zRes; i++)
@@ -906,8 +1050,9 @@ void FLUID_3D::project()
 		}
 		// printf("Max pressure: %f, dx: %f\n", maxvalue, _dx);
 	}
-
+#endif
-	setObstaclePressure(_pressure, 0, _zRes);
+	// DG TODO: check this function, for now this is done in the next function
+	// setObstaclePressure(_pressure, 0, _zRes);
 
 	// project out solution
 	float invDx = 1.0f / _dx;
@@ -916,17 +1061,55 @@ void FLUID_3D::project()
 		for (y = 1; y < _yRes - 1; y++, index += 2)
 			for (x = 1; x < _xRes - 1; x++, index++)
 			{
+				float vMask[3] = {1.0f, 1.0f, 1.0f}, vObst[3] = {0, 0, 0};
+				float vR = 0.0f, vL = 0.0f, vT = 0.0f, vB = 0.0f, vD = 0.0f, vU = 0.0f;
+#if USE_NEW_CG == 1
+				float pC = result[gti(FINDEX(x, y, z))]; // center
+				float pR = result[gti(FINDEX(x + 1, y, z))]; // right
+				float pL = result[gti(FINDEX(x - 1, y, z))]; // left
+				float pU = result[gti(FINDEX(x, y + 1, z))]; // top
+				float pD = result[gti(FINDEX(x, y - 1, z))]; // bottom
+				float pT = result[gti(FINDEX(x, y, z + 1))]; // Up
+				float pB = result[gti(FINDEX(x, y, z - 1))]; // Down
+#else
+
+				float pC = _pressure[index]; // center
+				float pR = _pressure[index + 1]; // right
+				float pL = _pressure[index - 1]; // left
+				float pU = _pressure[index + _xRes]; // Up
+				float pD = _pressure[index - _xRes]; // Down
+				float pT = _pressure[index + _slabSize]; // top
+				float pB = _pressure[index - _slabSize]; // bottom
+#endif
 				if(!_obstacles[index])
 				{
-					_xVelocity[index] -= 0.5f * (_pressure[index + 1]     - _pressure[index - 1]) * invDx;
+					// DG TODO: What if obstacle is left + right and one is moving?
-					_yVelocity[index] -= 0.5f * (_pressure[index + _xRes]  - _pressure[index - _xRes]) * invDx;
+					if(_obstacles[index+1])			{ pR = pC; vObst[0] = _xVelocityOb[index + 1];			vMask[0] = 0; }
-					_zVelocity[index] -= 0.5f * (_pressure[index + _slabSize] - _pressure[index - _slabSize]) * invDx;
+					if(_obstacles[index-1])			{ pL = pC; vObst[0]	= _xVelocityOb[index - 1];			vMask[0] = 0; }
+					if(_obstacles[index+_xRes])		{ pU = pC; vObst[1]	= _yVelocityOb[index + _xRes];		vMask[1] = 0; }
+					if(_obstacles[index-_xRes])		{ pD = pC; vObst[1]	= _yVelocityOb[index - _xRes];		vMask[1] = 0; }
+					if(_obstacles[index+_slabSize]) { pT = pC; vObst[2] = _zVelocityOb[index + _slabSize];	vMask[2] = 0; }
+					if(_obstacles[index-_slabSize]) { pB = pC; vObst[2]	= _zVelocityOb[index - _slabSize];	vMask[2] = 0; }
+
+					_xVelocity[index] -= 0.5f * (pR - pL) * invDx;
+					_yVelocity[index] -= 0.5f * (pU - pD) * invDx;
+					_zVelocity[index] -= 0.5f * (pT - pB) * invDx;
+
+					_xVelocity[index] = (vMask[0] * _xVelocity[index]) + vObst[0];
+					_yVelocity[index] = (vMask[1] * _yVelocity[index]) + vObst[1];
+					_zVelocity[index] = (vMask[2] * _zVelocity[index]) + vObst[2];
+				}
+				else
+				{
+					_xVelocity[index] = _xVelocityOb[index];
+					_yVelocity[index] = _yVelocityOb[index];
+					_zVelocity[index] = _zVelocityOb[index];
 				}
 			}
 
-	setObstacleVelocity(0, _zRes);
+	// setObstacleVelocity(0, _zRes);
 
-	if (_pressure) delete[] _pressure;
+	// if (_pressure) delete[] _pressure;
 	if (_divergence) delete[] _divergence;
 }
 
@@ -1154,8 +1337,9 @@ void FLUID_3D::setObstacleBoundaries(float *_pressure, int zBegin, int zEnd)
 					if (top)  counter++;
 					if (bottom) counter++;
 
-					if (counter < 3)
+					// DG TODO: Do not shrink obstacles for now
-						_obstacles[index] = EMPTY;
+					// if (counter < 3)
+					//	_obstacles[index] = EMPTY;
 				}
 				if (_obstacles[index])
 				{
@@ -1190,6 +1374,9 @@ void FLUID_3D::addBuoyancy(float *heat, float *density, int zBegin, int zEnd)
 //////////////////////////////////////////////////////////////////////
 // add vorticity to the force field
 //////////////////////////////////////////////////////////////////////
+#define VORT_VEL(i, j) \
+	((_obstacles[obpos[(i)]] & 8) ? ((abs(objvelocity[(j)][obpos[(i)]]) > FLT_EPSILON) ? objvelocity[(j)][obpos[(i)]] : velocity[(j)][index]) : velocity[(j)][obpos[(i)]])
+
 void FLUID_3D::addVorticity(int zBegin, int zEnd)
 {
 	//int x,y,z,index;
@@ -1225,9 +1412,18 @@ void FLUID_3D::addVorticity(int zBegin, int zEnd)
 	float gridSize = 0.5f / _dx;
 	//index = _slabSize + _xRes + 1;
 
+	float *velocity[3];
+	float *objvelocity[3];
+
+	velocity[0] = _xVelocity;
+	velocity[1] = _yVelocity;
+	velocity[2] = _zVelocity;
+
+	objvelocity[0] = _xVelocityOb;
+	objvelocity[1] = _yVelocityOb;
+	objvelocity[2] = _zVelocityOb;
 
 	size_t vIndex=_xRes + 1;
-
 	for (int z = zBegin + bb1; z < (zEnd - bt1); z++)
 	{
 		size_t index = index_ +(z-1)*_slabSize;
@@ -1237,25 +1433,47 @@ void FLUID_3D::addVorticity(int zBegin, int zEnd)
 		{
 			for (int x = 1; x < _xRes - 1; x++, index++)
 			{
+				if (!_obstacles[index])
+				{
+					int obpos[6];
 
-				int up    = _obstacles[index + _xRes] ? index : index + _xRes;
+					obpos[0] = (_obstacles[index + _xRes] == 1) ? index : index + _xRes; // up
-				int down  = _obstacles[index - _xRes] ? index : index - _xRes;
+					obpos[1] = (_obstacles[index - _xRes] == 1) ? index : index - _xRes; // down
-				float dy  = (up == index || down == index) ? 1.0f / _dx : gridSize;
+					float dy = (obpos[0] == index || obpos[1] == index) ? 1.0f / _dx : gridSize;
-				int out   = _obstacles[index + _slabSize] ? index : index + _slabSize;
-				int in    = _obstacles[index - _slabSize] ? index : index - _slabSize;
-				float dz  = (out == index || in == index) ? 1.0f / _dx : gridSize;
-				int right = _obstacles[index + 1] ? index : index + 1;
-				int left  = _obstacles[index - 1] ? index : index - 1;
-				float dx  = (right == index || left == index) ? 1.0f / _dx : gridSize;
 
-				_xVorticity[vIndex] = (_zVelocity[up] - _zVelocity[down]) * dy + (-_yVelocity[out] + _yVelocity[in]) * dz;
+					obpos[2]  = (_obstacles[index + _slabSize] == 1) ? index : index + _slabSize; // out
-				_yVorticity[vIndex] = (_xVelocity[out] - _xVelocity[in]) * dz + (-_zVelocity[right] + _zVelocity[left]) * dx;
+					obpos[3]  = (_obstacles[index - _slabSize] == 1) ? index : index - _slabSize; // in
-				_zVorticity[vIndex] = (_yVelocity[right] - _yVelocity[left]) * dx + (-_xVelocity[up] + _xVelocity[down])* dy;
+					float dz  = (obpos[2] == index || obpos[3] == index) ? 1.0f / _dx : gridSize;
 
-				_vorticity[vIndex] = sqrtf(_xVorticity[vIndex] * _xVorticity[vIndex] +
+					obpos[4] = (_obstacles[index + 1] == 1) ? index : index + 1; // right
+					obpos[5] = (_obstacles[index - 1] == 1) ? index : index - 1; // left
+					float dx = (obpos[4] == index || obpos[5] == index) ? 1.0f / _dx : gridSize;
+
+					float xV[2], yV[2], zV[2];
+
+					zV[1] = VORT_VEL(0, 2);
+					zV[0] = VORT_VEL(1, 2);
+					yV[1] = VORT_VEL(2, 1);
+					yV[0] = VORT_VEL(3, 1);
+					_xVorticity[vIndex] = (zV[1] - zV[0]) * dy + (-yV[1] + yV[0]) * dz;
+
+					xV[1] = VORT_VEL(2, 0);
+					xV[0] = VORT_VEL(3, 0);
+					zV[1] = VORT_VEL(4, 2);
+					zV[0] = VORT_VEL(5, 2);
+					_yVorticity[vIndex] = (xV[1] - xV[0]) * dz + (-zV[1] + zV[0]) * dx;
+
+					yV[1] = VORT_VEL(4, 1);
+					yV[0] = VORT_VEL(5, 1);
+					xV[1] = VORT_VEL(0, 0);
+					xV[0] = VORT_VEL(1, 0);
+					_zVorticity[vIndex] = (yV[1] - yV[0]) * dx + (-xV[1] + xV[0])* dy;
+
+					_vorticity[vIndex] = sqrtf(_xVorticity[vIndex] * _xVorticity[vIndex] +
 						_yVorticity[vIndex] * _yVorticity[vIndex] +
 						_zVorticity[vIndex] * _zVorticity[vIndex]);
 
+				}
 				vIndex++;
 			}
 			vIndex+=2;
@@ -1284,15 +1502,18 @@ void FLUID_3D::addVorticity(int zBegin, int zEnd)
 				{
 					float N[3];
 
-					int up    = _obstacles[index + _xRes] ? vIndex : vIndex + _xRes;
+					int up    = (_obstacles[index + _xRes] == 1) ? vIndex : vIndex + _xRes;
-					int down  = _obstacles[index - _xRes] ? vIndex : vIndex - _xRes;
+					int down  = (_obstacles[index - _xRes] == 1) ? vIndex : vIndex - _xRes;
 					float dy  = (up == vIndex || down == vIndex) ? 1.0f / _dx : gridSize;
-					int out   = _obstacles[index + _slabSize] ? vIndex : vIndex + _slabSize;
+
-					int in    = _obstacles[index - _slabSize] ? vIndex : vIndex - _slabSize;
+					int out   = (_obstacles[index + _slabSize] == 1) ? vIndex : vIndex + _slabSize;
+					int in    = (_obstacles[index - _slabSize] == 1) ? vIndex : vIndex - _slabSize;
 					float dz  = (out == vIndex || in == vIndex) ? 1.0f / _dx : gridSize;
-					int right = _obstacles[index + 1] ? vIndex : vIndex + 1;
+
-					int left  = _obstacles[index - 1] ? vIndex : vIndex - 1;
+					int right = (_obstacles[index + 1] == 1) ? vIndex : vIndex + 1;
+					int left  = (_obstacles[index - 1] == 1) ? vIndex : vIndex - 1;
 					float dx  = (right == vIndex || left == vIndex) ? 1.0f / _dx : gridSize;
+
 					N[0] = (_vorticity[right] - _vorticity[left]) * dx;
 					N[1] = (_vorticity[up] - _vorticity[down]) * dy;
 					N[2] = (_vorticity[out] - _vorticity[in]) * dz;
@@ -1309,20 +1530,20 @@ void FLUID_3D::addVorticity(int zBegin, int zEnd)
 						_yForce[index] += (N[2] * _xVorticity[vIndex] - N[0] * _zVorticity[vIndex]) * _dx * eps;
 						_zForce[index] += (N[0] * _yVorticity[vIndex] - N[1] * _xVorticity[vIndex]) * _dx * eps;
 					}
-					}	// if
+				}	// if
-					vIndex++;
+				vIndex++;
-					}	// x loop
+			}	// x loop
-				vIndex+=2;
+			vIndex+=2;
-				}		// y loop
+		}		// y loop
-			//vIndex+=2*_xRes;
+		//vIndex+=2*_xRes;
-		}				// z loop
+	}				// z loop
 
 	if (_xVorticity) delete[] _xVorticity;
 	if (_yVorticity) delete[] _yVorticity;
 	if (_zVorticity) delete[] _zVorticity;
 	if (_vorticity) delete[] _vorticity;
 }
-
+#undef CURL_VEL
 
 void FLUID_3D::advectMacCormackBegin(int zBegin, int zEnd)
 {
@@ -1393,38 +1614,10 @@ void FLUID_3D::advectMacCormackEnd2(int zBegin, int zEnd)
 
 	setZeroBorder(_density, res, zBegin, zEnd);
 	setZeroBorder(_heat, res, zBegin, zEnd);
-#if 0
-	{
-	const size_t index_ = _slabSize + _xRes + 1;
-	int bb=0;
-	int bt=0;
-
-	if (zBegin == 0) {bb = 1;}
-	if (zEnd == _zRes) {bt = 1;}
-	
-	for (int z = zBegin + bb; z < zEnd - bt; z++)
-	{
-		size_t index = index_ +(z-1)*_slabSize;
-
-		for (int y = 1; y < _yRes - 1; y++, index += 2)
-		{
-			for (int x = 1; x < _xRes - 1; x++, index++)
-			{
-				// clean custom velocities from moving obstacles again
-				if (_obstacles[index])
-				{
-					_xVelocity[index] =
-					_yVelocity[index] =
-					_zVelocity[index] = 0.0f;
-				}
-			}
-		}
-	}
-	}
-#endif
 
 	/*int begin=zBegin * _slabSize;
 	int end=begin + (zEnd - zBegin) * _slabSize;
   for (int x = begin; x < end; x++)
     _xForce[x] = _yForce[x] = 0.0f;*/
 }
+

intern/smoke/intern/FLUID_3D.h

@@ -43,6 +43,16 @@ using namespace std;
 using namespace BasicVector;
 class WTURBULENCE;
 
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+using namespace Eigen;
+
+#define USE_NEW_CG 0
+
+// Fluid index
+#define FINDEX(indexX, indexY, indexZ) ((indexZ) * _xRes * _yRes + (indexY) * _xRes + (indexX))
+
 class FLUID_3D  
 {
 	public:
@@ -105,9 +115,12 @@ class FLUID_3D
 		float* _xForce;
 		float* _yForce;
 		float* _zForce;
-		unsigned char*  _obstacles; /* only used (useful) for static obstacles like domain boundaries */
+		unsigned char*  _obstacles; /* only used (usefull) for static obstacles like domain boundaries */
 		unsigned char*  _obstaclesAnim;
 
+		// for debug purposes
+		float *_pressure;
+
 		// Required for proper threading:
 		float* _xVelocityTemp;
 		float* _yVelocityTemp;
@@ -159,7 +172,23 @@ class FLUID_3D
 		void project();
 		void diffuseHeat();
 		void solvePressure(float* field, float* b, unsigned char* skip);
+		
+#if USE_NEW_CG == 1
+		void solvePressurePre(VectorXf &b, SparseMatrix<float,RowMajor> &A, ArrayXd &gti, VectorXf &result);
+#else
 		void solvePressurePre(float* field, float* b, unsigned char* skip);
+
+		// diagonal preconditioner
+		void precond_init_diag(float *precond, float* field, unsigned char* skip, size_t index);
+		void precond_apply_diag(float *dest, float* source, float *precond, unsigned char* skip, size_t index);
+
+		// modified incomplete cholesky preconditioner
+		void precond_init_mic(float *precond, float* field, unsigned char* skip, size_t index);
+		void precond_apply_mic(float *dest, float* source, float *precond, unsigned char* skip, size_t index);
+#endif
+
+		void solvePressureJacobian(float* p, float* d, unsigned char* ob);
+
 		void solveHeat(float* field, float* b, unsigned char* skip);
 
 

intern/smoke/intern/FLUID_3D_SOLVERS.cpp

@@ -1,6 +1,6 @@
 /** \file smoke/intern/FLUID_3D_SOLVERS.cpp
- *  \ingroup smoke
+*  \ingroup smoke
- */
+*/
 //////////////////////////////////////////////////////////////////////
 // This file is part of Wavelet Turbulence.
 // 
@@ -50,121 +50,167 @@ void FLUID_3D::solveHeat(float* field, float* b, unsigned char* skip)
 	_q            = new float[_totalCells]; // set 0
 	_Acenter       = new float[_totalCells]; // set 0
 
-  	memset(_residual, 0, sizeof(float)*_totalCells);
+	memset(_residual, 0, sizeof(float)*_totalCells);
 	memset(_q, 0, sizeof(float)*_totalCells);
 	memset(_direction, 0, sizeof(float)*_totalCells);
 	memset(_Acenter, 0, sizeof(float)*_totalCells);
 
 	float deltaNew = 0.0f;
 
-  // r = b - Ax
+	// r = b - Ax
-  index = _slabSize + _xRes + 1;
+	index = _slabSize + _xRes + 1;
-  for (z = 1; z < _zRes - 1; z++, index += twoxr)
+	for (z = 1; z < _zRes - 1; z++, index += twoxr)
-    for (y = 1; y < _yRes - 1; y++, index += 2)
+		for (y = 1; y < _yRes - 1; y++, index += 2)
-      for (x = 1; x < _xRes - 1; x++, index++)
+			for (x = 1; x < _xRes - 1; x++, index++)
-      {
+			{
-        // if the cell is a variable
+				// if the cell is a variable
-        _Acenter[index] = 1.0f;
+				_Acenter[index] = 1.0f;
-        if (!skip[index])
+				if (!skip[index])
-        {
+				{
-          // set the matrix to the Poisson stencil in order
+					// set the matrix to the Poisson stencil in order
-          if (!skip[index + 1]) _Acenter[index] += heatConst;
+					if (!skip[index + 1]) _Acenter[index] += heatConst;
-          if (!skip[index - 1]) _Acenter[index] += heatConst;
+					if (!skip[index - 1]) _Acenter[index] += heatConst;
-          if (!skip[index + _xRes]) _Acenter[index] += heatConst;
+					if (!skip[index + _xRes]) _Acenter[index] += heatConst;
-          if (!skip[index - _xRes]) _Acenter[index] += heatConst;
+					if (!skip[index - _xRes]) _Acenter[index] += heatConst;
-          if (!skip[index + _slabSize]) _Acenter[index] += heatConst;
+					if (!skip[index + _slabSize]) _Acenter[index] += heatConst;
-          if (!skip[index - _slabSize]) _Acenter[index] += heatConst;
+					if (!skip[index - _slabSize]) _Acenter[index] += heatConst;
 
-		  _residual[index] = b[index] - (_Acenter[index] * field[index] + 
+					_residual[index] = b[index] - (_Acenter[index] * field[index] + 
-          field[index - 1] * (skip[index - 1] ? 0.0 : -heatConst) + 
+						field[index - 1] * (skip[index - 1] ? 0.0 : -heatConst) + 
-          field[index + 1] * (skip[index + 1] ? 0.0 : -heatConst) +
+						field[index + 1] * (skip[index + 1] ? 0.0 : -heatConst) +
-          field[index - _xRes] * (skip[index - _xRes] ? 0.0 : -heatConst) + 
+						field[index - _xRes] * (skip[index - _xRes] ? 0.0 : -heatConst) + 
-          field[index + _xRes] * (skip[index + _xRes] ? 0.0 : -heatConst) +
+						field[index + _xRes] * (skip[index + _xRes] ? 0.0 : -heatConst) +
-          field[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -heatConst) + 
+						field[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -heatConst) + 
-          field[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -heatConst));
+						field[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -heatConst));
-        }
+				}
-		else
+				else
-		{
+				{
-        _residual[index] = 0.0f;
+					_residual[index] = 0.0f;
-		}
+				}
 
-		_direction[index] = _residual[index];
+				_direction[index] = _residual[index];
-		deltaNew += _residual[index] * _residual[index];
+				deltaNew += _residual[index] * _residual[index];
-      }
+			}
 
 
-  // While deltaNew > (eps^2) * delta0
+			// While deltaNew > (eps^2) * delta0
-  const float eps  = SOLVER_ACCURACY;
+			const float eps  = SOLVER_ACCURACY;
-  float maxR = 2.0f * eps;
+			float maxR = 2.0f * eps;
-  while ((i < _iterations) && (maxR > eps))
+			while ((i < _iterations) && (maxR > eps))
-  {
+			{
-    // q = Ad
+				// q = Ad
-	float alpha = 0.0f;
+				float alpha = 0.0f;
 
-    index = _slabSize + _xRes + 1;
+				index = _slabSize + _xRes + 1;
-    for (z = 1; z < _zRes - 1; z++, index += twoxr)
+				for (z = 1; z < _zRes - 1; z++, index += twoxr)
-      for (y = 1; y < _yRes - 1; y++, index += 2)
+					for (y = 1; y < _yRes - 1; y++, index += 2)
-        for (x = 1; x < _xRes - 1; x++, index++)
+						for (x = 1; x < _xRes - 1; x++, index++)
-        {
+						{
-          // if the cell is a variable
+							// if the cell is a variable
-          if (!skip[index])
+							if (!skip[index])
-          {
+							{
 
-			_q[index] = (_Acenter[index] * _direction[index] + 
+								_q[index] = (_Acenter[index] * _direction[index] + 
-            _direction[index - 1] * (skip[index - 1] ? 0.0 : -heatConst) + 
+									_direction[index - 1] * (skip[index - 1] ? 0.0 : -heatConst) + 
-            _direction[index + 1] * (skip[index + 1] ? 0.0 : -heatConst) +
+									_direction[index + 1] * (skip[index + 1] ? 0.0 : -heatConst) +
-            _direction[index - _xRes] * (skip[index - _xRes] ? 0.0 : -heatConst) + 
+									_direction[index - _xRes] * (skip[index - _xRes] ? 0.0 : -heatConst) + 
-            _direction[index + _xRes] * (skip[index + _xRes] ? 0.0 : -heatConst) +
+									_direction[index + _xRes] * (skip[index + _xRes] ? 0.0 : -heatConst) +
-            _direction[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -heatConst) + 
+									_direction[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -heatConst) + 
-            _direction[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -heatConst));
+									_direction[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -heatConst));
-          }
+							}
-		  else
+							else
-		  {
+							{
-          _q[index] = 0.0f;
+								_q[index] = 0.0f;
-		  }
+							}
-		  alpha += _direction[index] * _q[index];
+							alpha += _direction[index] * _q[index];
-        }
+						}
 
-    if (fabs(alpha) > 0.0f)
+						if (fabs(alpha) > 0.0f)
-      alpha = deltaNew / alpha;
+							alpha = deltaNew / alpha;
 
-	float deltaOld = deltaNew;
+						float deltaOld = deltaNew;
-	deltaNew = 0.0f;
+						deltaNew = 0.0f;
 
-	maxR = 0.0f;
+						maxR = 0.0f;
 
-    index = _slabSize + _xRes + 1;
+						index = _slabSize + _xRes + 1;
-    for (z = 1; z < _zRes - 1; z++, index += twoxr)
+						for (z = 1; z < _zRes - 1; z++, index += twoxr)
-      for (y = 1; y < _yRes - 1; y++, index += 2)
+							for (y = 1; y < _yRes - 1; y++, index += 2)
-        for (x = 1; x < _xRes - 1; x++, index++)
+								for (x = 1; x < _xRes - 1; x++, index++)
-		{
+								{
-          field[index] += alpha * _direction[index];
+									field[index] += alpha * _direction[index];
 
-		  _residual[index] -= alpha * _q[index];
+									_residual[index] -= alpha * _q[index];
-          maxR = (_residual[index] > maxR) ? _residual[index] : maxR;
+									maxR = (_residual[index] > maxR) ? _residual[index] : maxR;
 
-		  deltaNew += _residual[index] * _residual[index];
+									deltaNew += _residual[index] * _residual[index];
-		}
+								}
 
-    float beta = deltaNew / deltaOld;
+								float beta = deltaNew / deltaOld;
 
-    index = _slabSize + _xRes + 1;
+								index = _slabSize + _xRes + 1;
-    for (z = 1; z < _zRes - 1; z++, index += twoxr)
+								for (z = 1; z < _zRes - 1; z++, index += twoxr)
-      for (y = 1; y < _yRes - 1; y++, index += 2)
+									for (y = 1; y < _yRes - 1; y++, index += 2)
-        for (x = 1; x < _xRes - 1; x++, index++)
+										for (x = 1; x < _xRes - 1; x++, index++)
-         _direction[index] = _residual[index] + beta * _direction[index];
+											_direction[index] = _residual[index] + beta * _direction[index];
 
 
-    i++;
+								i++;
-  }
+			}
-  // cout << i << " iterations converged to " << maxR << endl;
+			// cout << i << " iterations converged to " << maxR << endl;
 
-	if (_residual) delete[] _residual;
+			if (_residual) delete[] _residual;
-	if (_direction) delete[] _direction;
+			if (_direction) delete[] _direction;
-	if (_q)       delete[] _q;
+			if (_q)       delete[] _q;
-	if (_Acenter)  delete[] _Acenter;
+			if (_Acenter)  delete[] _Acenter;
 }
 
+#define SMOKE_LOOP \
+	index = _slabSize + _xRes + 1; \
+	for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes) \
+	for (y = 1; y < _yRes - 1; y++, index += 2) \
+	for (x = 1; x < _xRes - 1; x++, index++) 
+
+void FLUID_3D::precond_init_diag(float *precond, float* source, unsigned char* skip, size_t index)
+{
+	// if the cell is a variable
+	float Acenter = 0.0f;
+	if (!skip[index])
+	{
+		// set the matrix to the Poisson stencil in order
+		if (!skip[index + 1]) Acenter += 1.;
+		if (!skip[index - 1]) Acenter += 1.;
+		if (!skip[index + _xRes]) Acenter += 1.;
+		if (!skip[index - _xRes]) Acenter += 1.;
+		if (!skip[index + _slabSize]) Acenter += 1.;
+		if (!skip[index - _slabSize]) Acenter += 1.;
+	}
+
+	// P^-1
+	if(Acenter < 1.0)
+		precond[index] = 0.0;
+	else
+		precond[index] = 1.0 / Acenter;
+}
+
+void FLUID_3D::precond_apply_diag(float *dest, float* source, float *precond, unsigned char* skip, size_t index)
+{
+	dest[index] = source[index] * precond[index];
+}
+
+#define SQUARE(a) ((a)*(a))
+void FLUID_3D::precond_init_mic(float *precond, float* source, unsigned char* skip, size_t index)
+{
+	const Real tau = 0.97;
+	const Real sigma = 0.25;
+
+	// TODO
+}
+
+void FLUID_3D::precond_apply_mic(float *dest, float* source, float *precond, unsigned char* skip, size_t index)
+{
+	// TODO
+}
 
 void FLUID_3D::solvePressurePre(float* field, float* b, unsigned char* skip)
 {
@@ -190,136 +236,142 @@ void FLUID_3D::solvePressurePre(float* field, float* b, unsigned char* skip)
 	float deltaNew = 0.0f;
 
 	// r = b - Ax
-	index = _slabSize + _xRes + 1;
+	SMOKE_LOOP
-	for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
+	{
-		for (y = 1; y < _yRes - 1; y++, index += 2)
+		// r = b - Ax, with x = 0
-		  for (x = 1; x < _xRes - 1; x++, index++)
+		_residual[index] = b[index];
-		  {
+
+		// init diagonal preconditioner
+		precond_init_diag(_Precond, field, skip, index);
+
+		// p = P^-1 * r
+		precond_apply_diag(_direction, _residual, _Precond, skip, index);
+		// _direction[index] = _residual[index] * _Precond[index];
+
+		deltaNew += _residual[index] * _direction[index];
+	}
+
+
+	// While deltaNew > (eps^2) * delta0
+	const float eps  = SOLVER_ACCURACY;
+	//while ((i < _iterations) && (deltaNew > eps*delta0))
+	float maxR = 2.0f * eps;
+	// while (i < _iterations)
+	while ((i < _iterations) && (maxR > 0.001*eps))
+	{
+
+		float alpha = 0.0f;
+
+		SMOKE_LOOP
+		{
 			// if the cell is a variable
 			float Acenter = 0.0f;
 			if (!skip[index])
 			{
-			  // set the matrix to the Poisson stencil in order
+				// set the matrix to the Poisson stencil in order
-			  if (!skip[index + 1]) Acenter += 1.;
+				if (!skip[index + 1]) Acenter += 1.;
-			  if (!skip[index - 1]) Acenter += 1.;
+				if (!skip[index - 1]) Acenter += 1.;
-			  if (!skip[index + _xRes]) Acenter += 1.;
+				if (!skip[index + _xRes]) Acenter += 1.;
-			  if (!skip[index - _xRes]) Acenter += 1.;
+				if (!skip[index - _xRes]) Acenter += 1.;
-			  if (!skip[index + _slabSize]) Acenter += 1.;
+				if (!skip[index + _slabSize]) Acenter += 1.;
-			  if (!skip[index - _slabSize]) Acenter += 1.;
+				if (!skip[index - _slabSize]) Acenter += 1.;
 
-			  _residual[index] = b[index] - (Acenter * field[index] +  
+				_q[index] = Acenter * _direction[index] +  
-			  field[index - 1] * (skip[index - 1] ? 0.0 : -1.0f)+ 
+					_direction[index - 1] * (skip[index - 1] ? 0.0 : -1.0f) + 
-			  field[index + 1] * (skip[index + 1] ? 0.0 : -1.0f)+
+					_direction[index + 1] * (skip[index + 1] ? 0.0 : -1.0f) +
-			  field[index - _xRes] * (skip[index - _xRes] ? 0.0 : -1.0f)+ 
+					_direction[index - _xRes] * (skip[index - _xRes] ? 0.0 : -1.0f) + 
-			  field[index + _xRes] * (skip[index + _xRes] ? 0.0 : -1.0f)+
+					_direction[index + _xRes] * (skip[index + _xRes] ? 0.0 : -1.0f)+
-			  field[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -1.0f)+ 
+					_direction[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -1.0f) + 
-			  field[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -1.0f) );
+					_direction[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -1.0f);
 			}
 			else
 			{
-			_residual[index] = 0.0f;
+				_q[index] = 0.0f;
 			}
 
-			// P^-1
+			alpha += _direction[index] * _q[index];
-			if(Acenter < 1.0)
-				_Precond[index] = 0.0;
-			else
-				_Precond[index] = 1.0 / Acenter;
-
-			// p = P^-1 * r
-			_direction[index] = _residual[index] * _Precond[index];
-
-			deltaNew += _residual[index] * _direction[index];
-		  }
-
-
-  // While deltaNew > (eps^2) * delta0
-  const float eps  = SOLVER_ACCURACY;
-  //while ((i < _iterations) && (deltaNew > eps*delta0))
-  float maxR = 2.0f * eps;
-  // while (i < _iterations)
-  while ((i < _iterations) && (maxR > 0.001*eps))
-  {
-
-	float alpha = 0.0f;
-
-    index = _slabSize + _xRes + 1;
-    for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-      for (y = 1; y < _yRes - 1; y++, index += 2)
-        for (x = 1; x < _xRes - 1; x++, index++)
-        {
-          // if the cell is a variable
-          float Acenter = 0.0f;
-          if (!skip[index])
-          {
-            // set the matrix to the Poisson stencil in order
-            if (!skip[index + 1]) Acenter += 1.;
-            if (!skip[index - 1]) Acenter += 1.;
-            if (!skip[index + _xRes]) Acenter += 1.;
-            if (!skip[index - _xRes]) Acenter += 1.;
-            if (!skip[index + _slabSize]) Acenter += 1.;
-            if (!skip[index - _slabSize]) Acenter += 1.;
-
-			_q[index] = Acenter * _direction[index] +  
-            _direction[index - 1] * (skip[index - 1] ? 0.0 : -1.0f) + 
-            _direction[index + 1] * (skip[index + 1] ? 0.0 : -1.0f) +
-            _direction[index - _xRes] * (skip[index - _xRes] ? 0.0 : -1.0f) + 
-            _direction[index + _xRes] * (skip[index + _xRes] ? 0.0 : -1.0f)+
-            _direction[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -1.0f) + 
-            _direction[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -1.0f);
-          }
-		  else
-		  {
-          _q[index] = 0.0f;
-		  }
-
-		  alpha += _direction[index] * _q[index];
-        }
-
-
-    if (fabs(alpha) > 0.0f)
-      alpha = deltaNew / alpha;
-
-	float deltaOld = deltaNew;
-	deltaNew = 0.0f;
-
-	maxR = 0.0;
-
-	float tmp;
-
-    // x = x + alpha * d
-    index = _slabSize + _xRes + 1;
-    for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-      for (y = 1; y < _yRes - 1; y++, index += 2)
-        for (x = 1; x < _xRes - 1; x++, index++)
-		{
-          field[index] += alpha * _direction[index];
-
-		  _residual[index] -= alpha * _q[index];
-
-		  _h[index] = _Precond[index] * _residual[index];
-
-		  tmp = _residual[index] * _h[index];
-		  deltaNew += tmp;
-		  maxR = (tmp > maxR) ? tmp : maxR;
-
 		}
 
 
-    // beta = deltaNew / deltaOld
+		if (fabs(alpha) > 0.0f)
-    float beta = deltaNew / deltaOld;
+			alpha = deltaNew / alpha;
 
-    // d = h + beta * d
+		float deltaOld = deltaNew;
-    index = _slabSize + _xRes + 1;
+		deltaNew = 0.0f;
-    for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
-      for (y = 1; y < _yRes - 1; y++, index += 2)
-        for (x = 1; x < _xRes - 1; x++, index++)
-          _direction[index] = _h[index] + beta * _direction[index];
 
-    // i = i + 1
+		maxR = 0.0;
-    i++;
+
-  }
+		float tmp;
-  // cout << i << " iterations converged to " << sqrt(maxR) << endl;
+
+		// x = x + alpha * d
+		SMOKE_LOOP
+		{
+			field[index] += alpha * _direction[index];
+		}
+
+		if(i % 50)
+		{
+			SMOKE_LOOP
+			{
+				// if the cell is a variable
+				float Acenter = 0.0f;
+				if (!skip[index])
+				{
+					// set the matrix to the Poisson stencil in order
+					if (!skip[index + 1]) Acenter += 1.;
+					if (!skip[index - 1]) Acenter += 1.;
+					if (!skip[index + _xRes]) Acenter += 1.;
+					if (!skip[index - _xRes]) Acenter += 1.;
+					if (!skip[index + _slabSize]) Acenter += 1.;
+					if (!skip[index - _slabSize]) Acenter += 1.;
+
+					_residual[index] = b[index] - (Acenter * field[index] +  
+						field[index - 1] * (skip[index - 1] ? 0.0 : -1.0f) + 
+						field[index + 1] * (skip[index + 1] ? 0.0 : -1.0f) +
+						field[index - _xRes] * (skip[index - _xRes] ? 0.0 : -1.0f) + 
+						field[index + _xRes] * (skip[index + _xRes] ? 0.0 : -1.0f)+
+						field[index - _slabSize] * (skip[index - _slabSize] ? 0.0 : -1.0f) + 
+						field[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -1.0f) );
+				}
+				else
+				{
+					_residual[index] = 0.0f;
+				}
+			}
+		}
+		else
+		{
+			SMOKE_LOOP
+			{
+				_residual[index] -= alpha * _q[index];
+			}
+		}
+
+		SMOKE_LOOP
+		{
+			// Apply preconditioner
+			precond_apply_diag(_h, _residual, _Precond, skip, index);
+			// _h[index] = _Precond[index] * _residual[index];
+
+			tmp = _residual[index] * _h[index];
+			deltaNew += tmp;
+			maxR = (tmp > maxR) ? tmp : maxR;
+		}
+
+
+		// beta = deltaNew / deltaOld
+		float beta = deltaNew / deltaOld;
+
+		// d = h + beta * d
+		SMOKE_LOOP
+		{
+			_direction[index] = _h[index] + beta * _direction[index];
+		}
+
+		// i = i + 1
+		i++;
+	}
+	// cout << i << " iterations converged to " << sqrt(maxR) << endl;
 
 	if (_h) delete[] _h;
 	if (_Precond) delete[] _Precond;

intern/smoke/intern/FLUID_3D_STATIC.cpp

@@ -1,6 +1,6 @@
 /** \file smoke/intern/FLUID_3D_STATIC.cpp
- *  \ingroup smoke
+*  \ingroup smoke
- */
+*/
 //////////////////////////////////////////////////////////////////////
 // This file is part of Wavelet Turbulence.
 //
@@ -39,12 +39,12 @@
 //////////////////////////////////////////////////////////////////////
 /*
 void FLUID_3D::addSmokeColumn() {
-	addSmokeTestCase(_density, _res, 1.0);
+addSmokeTestCase(_density, _res, 1.0);
-	// addSmokeTestCase(_zVelocity, _res, 1.0);
+// addSmokeTestCase(_zVelocity, _res, 1.0);
-	addSmokeTestCase(_heat, _res, 1.0);
+addSmokeTestCase(_heat, _res, 1.0);
-	if (_wTurbulence) {
+if (_wTurbulence) {
-		addSmokeTestCase(_wTurbulence->getDensityBig(), _wTurbulence->getResBig(), 1.0);
+addSmokeTestCase(_wTurbulence->getDensityBig(), _wTurbulence->getResBig(), 1.0);
-	}
+}
 }
 */
 
@@ -105,7 +105,7 @@ void FLUID_3D::setNeumannX(float* field, Vec3Int res, int zBegin, int zEnd)
 			index -= 1;
 			if(field[index]<0.) field[index] = 0.;
 		}
- }
+}
 
 //////////////////////////////////////////////////////////////////////
 // set y direction to Neumann boundary conditions
@@ -229,29 +229,29 @@ void FLUID_3D::setZeroZ(float* field, Vec3Int res, int zBegin, int zEnd)
 
 	int index = 0;
 	if ((zBegin == 0))
-	for (int y = 0; y < res[1]; y++)
-		for (int x = 0; x < res[0]; x++, index++)
-		{
-			// front slab
-			field[index] = 0.0f;
-    }
-
-	if (zEnd == res[2])
-	{
-		index=0;
-		int indexx=0;
-		const int cellsslab = totalCells - slabSize;
-
 		for (int y = 0; y < res[1]; y++)
 			for (int x = 0; x < res[0]; x++, index++)
 			{
-
+				// front slab
-				// back slab
+				field[index] = 0.0f;
-				indexx = index + cellsslab;
-				field[indexx] = 0.0f;
 			}
-	}
+
- }
+			if (zEnd == res[2])
+			{
+				index=0;
+				int indexx=0;
+				const int cellsslab = totalCells - slabSize;
+
+				for (int y = 0; y < res[1]; y++)
+					for (int x = 0; x < res[0]; x++, index++)
+					{
+
+						// back slab
+						indexx = index + cellsslab;
+						field[indexx] = 0.0f;
+					}
+			}
+}
 //////////////////////////////////////////////////////////////////////
 // copy grid boundary
 //////////////////////////////////////////////////////////////////////
@@ -294,34 +294,34 @@ void FLUID_3D::copyBorderZ(float* field, Vec3Int res, int zBegin, int zEnd)
 	int index=0;
 
 	if ((zBegin == 0))
-	for (int y = 0; y < res[1]; y++)
+		for (int y = 0; y < res[1]; y++)
-		for (int x = 0; x < res[0]; x++, index++)
+			for (int x = 0; x < res[0]; x++, index++)
-		{
+			{
-			field[index] = field[index + slabSize]; 
+				field[index] = field[index + slabSize]; 
-		}
+			}
 
-	if ((zEnd == res[2]))
+			if ((zEnd == res[2]))
-	{
+			{
 
-	index=0;
+				index=0;
-	int indexx=0;
+				int indexx=0;
-	const int cellsslab = totalCells - slabSize;
+				const int cellsslab = totalCells - slabSize;
 
-	for (int y = 0; y < res[1]; y++)
+				for (int y = 0; y < res[1]; y++)
-		for (int x = 0; x < res[0]; x++, index++)
+					for (int x = 0; x < res[0]; x++, index++)
-		{
+					{
-			// back slab
+						// back slab
-			indexx = index + cellsslab;
+						indexx = index + cellsslab;
-			field[indexx] = field[indexx - slabSize];
+						field[indexx] = field[indexx - slabSize];
-		}
+					}
-	}
+			}
 }
 
 /////////////////////////////////////////////////////////////////////
 // advect field with the semi lagrangian method
 //////////////////////////////////////////////////////////////////////
 void FLUID_3D::advectFieldSemiLagrange(const float dt, const float* velx, const float* vely,  const float* velz,
-		float* oldField, float* newField, Vec3Int res, int zBegin, int zEnd)
+									   float* oldField, float* newField, Vec3Int res, int zBegin, int zEnd)
 {
 	const int xres = res[0];
 	const int yres = res[1];
@@ -335,7 +335,7 @@ void FLUID_3D::advectFieldSemiLagrange(const float dt, const float* velx, const
 			{
 				const int index = x + y * xres + z * xres*yres;
 
-        // backtrace
+				// backtrace
 				float xTrace = x - dt * velx[index];
 				float yTrace = y - dt * vely[index];
 				float zTrace = z - dt * velz[index];
@@ -376,13 +376,13 @@ void FLUID_3D::advectFieldSemiLagrange(const float dt, const float* velx, const
 				// interpolate
 				// (indices could be computed once)
 				newField[index] = u0 * (s0 * (t0 * oldField[i000] +
-							t1 * oldField[i010]) +
+					t1 * oldField[i010]) +
-						s1 * (t0 * oldField[i100] +
+					s1 * (t0 * oldField[i100] +
-							t1 * oldField[i110])) +
+					t1 * oldField[i110])) +
 					u1 * (s0 * (t0 * oldField[i001] +
-								t1 * oldField[i011]) +
+					t1 * oldField[i011]) +
-							s1 * (t0 * oldField[i101] +
+					s1 * (t0 * oldField[i101] +
-								t1 * oldField[i111]));
+					t1 * oldField[i111]));
 			}
 }
 
@@ -393,14 +393,14 @@ void FLUID_3D::advectFieldSemiLagrange(const float dt, const float* velx, const
 // comments are the pseudocode from selle's paper
 //////////////////////////////////////////////////////////////////////
 void FLUID_3D::advectFieldMacCormack1(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity, 
-				float* oldField, float* tempResult, Vec3Int res, int zBegin, int zEnd)
+									  float* oldField, float* tempResult, Vec3Int res, int zBegin, int zEnd)
 {
 	/*const int sx= res[0];
 	const int sy= res[1];
 	const int sz= res[2];
 
 	for (int x = 0; x < sx * sy * sz; x++)
-		phiHatN[x] = phiHatN1[x] = oldField[x];*/	// not needed as all the values are written first
+	phiHatN[x] = phiHatN1[x] = oldField[x];*/	// not needed as all the values are written first
 
 	float*& phiN    = oldField;
 	float*& phiN1   = tempResult;
@@ -414,7 +414,7 @@ void FLUID_3D::advectFieldMacCormack1(const float dt, const float* xVelocity, co
 
 
 void FLUID_3D::advectFieldMacCormack2(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity, 
-				float* oldField, float* newField, float* tempResult, float* temp1, Vec3Int res, const unsigned char* obstacles, int zBegin, int zEnd)
+									  float* oldField, float* newField, float* tempResult, float* temp1, Vec3Int res, const unsigned char* obstacles, int zBegin, int zEnd)
 {
 	float* phiHatN  = tempResult;
 	float* t1  = temp1;
@@ -438,15 +438,15 @@ void FLUID_3D::advectFieldMacCormack2(const float dt, const float* xVelocity, co
 				phiN1[index] = phiHatN[index] + (phiN[index] - t1[index]) * 0.50f;
 				//phiN1[index] = phiHatN1[index]; // debug, correction off
 			}
-	copyBorderX(phiN1, res, zBegin, zEnd);
+			copyBorderX(phiN1, res, zBegin, zEnd);
-	copyBorderY(phiN1, res, zBegin, zEnd);
+			copyBorderY(phiN1, res, zBegin, zEnd);
-	copyBorderZ(phiN1, res, zBegin, zEnd);
+			copyBorderZ(phiN1, res, zBegin, zEnd);
 
-	// clamp any newly created extrema
+			// clamp any newly created extrema
-	clampExtrema(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res, zBegin, zEnd);		// uses wide data from old field and velocities (both are whole)
+			clampExtrema(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res, zBegin, zEnd);		// uses wide data from old field and velocities (both are whole)
 
-	// if the error estimate was bad, revert to first order
+			// if the error estimate was bad, revert to first order
-	clampOutsideRays(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res, obstacles, phiHatN, zBegin, zEnd);	// phiHatN is only used at cells within thread range, so its ok
+			clampOutsideRays(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res, obstacles, phiHatN, zBegin, zEnd);	// phiHatN is only used at cells within thread range, so its ok
 
 } 
 
@@ -455,7 +455,7 @@ void FLUID_3D::advectFieldMacCormack2(const float dt, const float* xVelocity, co
 // Clamp the extrema generated by the BFECC error correction
 //////////////////////////////////////////////////////////////////////
 void FLUID_3D::clampExtrema(const float dt, const float* velx, const float* vely,  const float* velz,
-		float* oldField, float* newField, Vec3Int res, int zBegin, int zEnd)
+							float* oldField, float* newField, Vec3Int res, int zBegin, int zEnd)
 {
 	const int xres= res[0];
 	const int yres= res[1];
@@ -539,7 +539,7 @@ void FLUID_3D::clampExtrema(const float dt, const float* velx, const float* vely
 // incorrect
 //////////////////////////////////////////////////////////////////////
 void FLUID_3D::clampOutsideRays(const float dt, const float* velx, const float* vely,  const float* velz,
-				float* oldField, float* newField, Vec3Int res, const unsigned char* obstacles, const float *oldAdvection, int zBegin, int zEnd)
+								float* oldField, float* newField, Vec3Int res, const unsigned char* obstacles, const float *oldAdvection, int zBegin, int zEnd)
 {
 	const int sx= res[0];
 	const int sy= res[1];
@@ -657,14 +657,14 @@ void FLUID_3D::clampOutsideRays(const float dt, const float* velx, const float*
 
 					// interpolate, (indices could be computed once)
 					newField[index] = u0 * (s0 * (
-								t0 * oldField[i000] +
+						t0 * oldField[i000] +
-								t1 * oldField[i010]) +
+						t1 * oldField[i010]) +
-							s1 * (t0 * oldField[i100] +
+						s1 * (t0 * oldField[i100] +
-								t1 * oldField[i110])) +
+						t1 * oldField[i110])) +
 						u1 * (s0 * (t0 * oldField[i001] +
-									t1 * oldField[i011]) +
+						t1 * oldField[i011]) +
-								s1 * (t0 * oldField[i101] +
+						s1 * (t0 * oldField[i101] +
-									t1 * oldField[i111])); 
+						t1 * oldField[i111])); 
 				}
 			} // xyz
 }

intern/smoke/intern/smoke_API.cpp

@@ -215,6 +215,11 @@ extern "C" void smoke_turbulence_export(WTURBULENCE *wt, float **dens, float **d
 	*tcw = wt->_tcW;
 }
 
+extern "C" float *smoke_get_pressure(FLUID_3D *fluid)
+{
+	return fluid->_pressure;
+}
+
 extern "C" float *smoke_get_density(FLUID_3D *fluid)
 {
 	return fluid->_density;

source/blender/blenkernel/intern/smoke.c

@@ -103,7 +103,7 @@ static void tend ( void )
 {
 	QueryPerformanceCounter ( &liCurrentTime );
 }
-static double UNUSED_FUNCTION(tval)( void )
+static double tval( void )
 {
 	return ((double)( (liCurrentTime.QuadPart - liStartTime.QuadPart)* (double)1000.0/(double)liFrequency.QuadPart ));
 }
@@ -120,7 +120,7 @@ static void tend ( void )
 	gettimeofday ( &_tend,&tz );
 }
 
-static double UNUSED_FUNCTION(tval)( void )
+static double tval( void )
 {
 	double t1, t2;
 	t1 = ( double ) _tstart.tv_sec*1000 + ( double ) _tstart.tv_usec/ ( 1000 );
@@ -140,9 +140,7 @@ struct SmokeModifierData;
 #define DT_DEFAULT 0.1f
 
 /* forward declerations */
-static void calcTriangleDivs(Object *ob, MVert *verts, int numverts, MFace *tris, int numfaces, int numtris, int **tridivs, float cell_len);
 static void get_cell(const float p0[3], const int res[3], float dx, const float pos[3], int cell[3], int correct);
-static void fill_scs_points(Object *ob, DerivedMesh *dm, SmokeCollSettings *scs);
 
 #else /* WITH_SMOKE */
 
@@ -280,444 +278,19 @@ static int smokeModifier_init (SmokeModifierData *smd, Object *ob, Scene *scene,
 	}
 	else if((smd->type & MOD_SMOKE_TYPE_COLL))
 	{
-		// todo: delete this when loading colls work -dg
-
 		if(!smd->coll)
 		{
 			smokeModifier_createType(smd);
 		}
 
-		if(!smd->coll->points)
+		smd->time = scene->r.cfra;
-		{
-			// init collision points
-			SmokeCollSettings *scs = smd->coll;
 
-			smd->time = scene->r.cfra;
-
-			// copy obmat
-			copy_m4_m4(scs->mat, ob->obmat);
-			copy_m4_m4(scs->mat_old, ob->obmat);
-
-			DM_ensure_tessface(dm);
-			fill_scs_points(ob, dm, scs);
-		}
-
-		if(!smd->coll->bvhtree)
-		{
-			smd->coll->bvhtree = NULL; // bvhtree_build_from_smoke ( ob->obmat, dm->getTessFaceArray(dm), dm->getNumTessFaces(dm), dm->getVertArray(dm), dm->getNumVerts(dm), 0.0 );
-		}
 		return 1;
 	}
 
 	return 2;
 }
 
-static void fill_scs_points(Object *ob, DerivedMesh *dm, SmokeCollSettings *scs)
-{
-	MVert *mvert = dm->getVertArray(dm);
-	MFace *mface = dm->getTessFaceArray(dm);
-	int i = 0, divs = 0;
-
-	// DG TODO: need to do this dynamically according to the domain object!
-	float cell_len = scs->dx;
-	int newdivs = 0;
-	int quads = 0, facecounter = 0;
-
-	// count quads
-	for(i = 0; i < dm->getNumTessFaces(dm); i++)
-	{
-		if(mface[i].v4)
-			quads++;
-	}
-
-	scs->numtris = dm->getNumTessFaces(dm) + quads;
-	scs->tridivs = NULL;
-	calcTriangleDivs(ob, mvert, dm->getNumVerts(dm), mface,  dm->getNumTessFaces(dm), scs->numtris, &(scs->tridivs), cell_len);
-
-	// count triangle divisions
-	for(i = 0; i < dm->getNumTessFaces(dm) + quads; i++)
-	{
-		divs += (scs->tridivs[3 * i] + 1) * (scs->tridivs[3 * i + 1] + 1) * (scs->tridivs[3 * i + 2] + 1);
-	}
-
-	scs->points = MEM_callocN(sizeof(float) * (dm->getNumVerts(dm) + divs) * 3, "SmokeCollPoints");
-	scs->points_old = MEM_callocN(sizeof(float) * (dm->getNumVerts(dm) + divs) * 3, "SmokeCollPointsOld");
-
-	for(i = 0; i < dm->getNumVerts(dm); i++)
-	{
-		float tmpvec[3];
-		copy_v3_v3(tmpvec, mvert[i].co);
-		// mul_m4_v3(ob->obmat, tmpvec); // DG: use local coordinates, we save MAT anyway
-		copy_v3_v3(&scs->points[i * 3], tmpvec);
-	}
-	
-	for(i = 0, facecounter = 0; i < dm->getNumTessFaces(dm); i++)
-	{
-		int again = 0;
-		do
-		{
-			int j, k;
-			int divs1 = scs->tridivs[3 * facecounter + 0];
-			int divs2 = scs->tridivs[3 * facecounter + 1];
-			//int divs3 = scs->tridivs[3 * facecounter + 2];
-			float side1[3], side2[3], trinormorg[3], trinorm[3];
-			
-			if(again == 1 && mface[i].v4)
-			{
-				sub_v3_v3v3(side1,  mvert[ mface[i].v3 ].co, mvert[ mface[i].v1 ].co);
-				sub_v3_v3v3(side2,  mvert[ mface[i].v4 ].co, mvert[ mface[i].v1 ].co);
-			}
-			else {
-				sub_v3_v3v3(side1,  mvert[ mface[i].v2 ].co, mvert[ mface[i].v1 ].co);
-				sub_v3_v3v3(side2,  mvert[ mface[i].v3 ].co, mvert[ mface[i].v1 ].co);
-			}
-
-			cross_v3_v3v3(trinormorg, side1, side2);
-			normalize_v3(trinormorg);
-			copy_v3_v3(trinorm, trinormorg);
-			mul_v3_fl(trinorm, 0.25 * cell_len);
-
-			for(j = 0; j <= divs1; j++)
-			{
-				for(k = 0; k <= divs2; k++)
-				{
-					float p1[3], p2[3], p3[3], p[3]={0,0,0}; 
-					const float uf = (float)(j + TRI_UVOFFSET) / (float)(divs1 + 0.0);
-					const float vf = (float)(k + TRI_UVOFFSET) / (float)(divs2 + 0.0);
-					float tmpvec[3];
-					
-					if(uf+vf > 1.0) 
-					{
-						// printf("bigger - divs1: %d, divs2: %d\n", divs1, divs2);
-						continue;
-					}
-
-					copy_v3_v3(p1, mvert[ mface[i].v1 ].co);
-					if(again == 1 && mface[i].v4)
-					{
-						copy_v3_v3(p2, mvert[ mface[i].v3 ].co);
-						copy_v3_v3(p3, mvert[ mface[i].v4 ].co);
-					}
-					else {
-						copy_v3_v3(p2, mvert[ mface[i].v2 ].co);
-						copy_v3_v3(p3, mvert[ mface[i].v3 ].co);
-					}
-
-					mul_v3_fl(p1, (1.0-uf-vf));
-					mul_v3_fl(p2, uf);
-					mul_v3_fl(p3, vf);
-					
-					add_v3_v3v3(p, p1, p2);
-					add_v3_v3(p, p3);
-
-					if(newdivs > divs)
-						printf("mem problem\n");
-
-					// mMovPoints.push_back(p + trinorm);
-					add_v3_v3v3(tmpvec, p, trinorm);
-					// mul_m4_v3(ob->obmat, tmpvec); // DG: use local coordinates, we save MAT anyway
-					copy_v3_v3(&scs->points[3 * (dm->getNumVerts(dm) + newdivs)], tmpvec);
-					newdivs++;
-
-					if(newdivs > divs)
-						printf("mem problem\n");
-
-					// mMovPoints.push_back(p - trinorm);
-					copy_v3_v3(tmpvec, p);
-					sub_v3_v3(tmpvec, trinorm);
-					// mul_m4_v3(ob->obmat, tmpvec); // DG: use local coordinates, we save MAT anyway
-					copy_v3_v3(&scs->points[3 * (dm->getNumVerts(dm) + newdivs)], tmpvec);
-					newdivs++;
-				}
-			}
-
-			if(again == 0 && mface[i].v4)
-				again++;
-			else
-				again = 0;
-
-			facecounter++;
-
-		} while(again!=0);
-	}
-
-	scs->numverts = dm->getNumVerts(dm);
-	// DG TODO: also save triangle count?
-
-	scs->numpoints = dm->getNumVerts(dm) + newdivs;
-
-	for(i = 0; i < scs->numpoints * 3; i++)
-	{
-		scs->points_old[i] = scs->points[i];
-	}
-}
-
-
-static void fill_scs_points_anim(Object *UNUSED(ob), DerivedMesh *dm, SmokeCollSettings *scs)
-{
-	MVert *mvert = dm->getVertArray(dm);
-	MFace *mface = dm->getTessFaceArray(dm);
-	int quads = 0, numtris = 0, facecounter = 0;
-	unsigned int i = 0;
-	int divs = 0, newdivs = 0;
-	
-	// DG TODO: need to do this dynamically according to the domain object!
-	float cell_len = scs->dx;
-
-	// count quads
-	for(i = 0; i < dm->getNumTessFaces(dm); i++)
-	{
-		if(mface[i].v4)
-			quads++;
-	}
-
-	numtris = dm->getNumTessFaces(dm) + quads;
-
-	// check if mesh changed topology
-	if(scs->numtris != numtris)
-		return;
-	if(scs->numverts != dm->getNumVerts(dm))
-		return;
-
-	// update new positions
-	for(i = 0; i < dm->getNumVerts(dm); i++)
-	{
-		float tmpvec[3];
-		copy_v3_v3(tmpvec, mvert[i].co);
-		copy_v3_v3(&scs->points[i * 3], tmpvec);
-	}
-
-	// for every triangle // update div points
-	for(i = 0, facecounter = 0; i < dm->getNumTessFaces(dm); i++)
-	{
-		int again = 0;
-		do
-		{
-			int j, k;
-			int divs1 = scs->tridivs[3 * facecounter + 0];
-			int divs2 = scs->tridivs[3 * facecounter + 1];
-			float srcside1[3], srcside2[3], destside1[3], destside2[3], src_trinormorg[3], dest_trinormorg[3], src_trinorm[3], dest_trinorm[3];
-			
-			if(again == 1 && mface[i].v4)
-			{
-				sub_v3_v3v3(srcside1,  &scs->points_old[mface[i].v3 * 3], &scs->points_old[mface[i].v1 * 3]);
-				sub_v3_v3v3(destside1,  &scs->points[mface[i].v3 * 3], &scs->points[mface[i].v1 * 3]);
-
-				sub_v3_v3v3(srcside2,  &scs->points_old[mface[i].v4 * 3], &scs->points_old[mface[i].v1 * 3]);
-				sub_v3_v3v3(destside2,  &scs->points[mface[i].v4 * 3], &scs->points[mface[i].v1 * 3]);
-			}
-			else {
-				sub_v3_v3v3(srcside1,  &scs->points_old[mface[i].v2 * 3], &scs->points_old[mface[i].v1 * 3]);
-				sub_v3_v3v3(destside1,  &scs->points[mface[i].v2 * 3], &scs->points[mface[i].v1 * 3]);
-
-				sub_v3_v3v3(srcside2,  &scs->points_old[mface[i].v3 * 3], &scs->points_old[mface[i].v1 * 3]);
-				sub_v3_v3v3(destside2,  &scs->points[mface[i].v3 * 3], &scs->points[mface[i].v1 * 3]);
-			}
-
-			cross_v3_v3v3(src_trinormorg, srcside1, srcside2);
-			cross_v3_v3v3(dest_trinormorg, destside1, destside2);
-
-			normalize_v3(src_trinormorg);
-			normalize_v3(dest_trinormorg);
-
-			copy_v3_v3(src_trinorm, src_trinormorg);
-			copy_v3_v3(dest_trinorm, dest_trinormorg);
-
-			mul_v3_fl(src_trinorm, 0.25 * cell_len);
-			mul_v3_fl(dest_trinorm, 0.25 * cell_len);
-
-			for(j = 0; j <= divs1; j++)
-			{
-				for(k = 0; k <= divs2; k++)
-				{
-					float src_p1[3], src_p2[3], src_p3[3], src_p[3]={0,0,0};
-					float dest_p1[3], dest_p2[3], dest_p3[3], dest_p[3]={0,0,0};
-					const float uf = (float)(j + TRI_UVOFFSET) / (float)(divs1 + 0.0);
-					const float vf = (float)(k + TRI_UVOFFSET) / (float)(divs2 + 0.0);
-					float src_tmpvec[3], dest_tmpvec[3];
-					
-					if(uf+vf > 1.0) 
-					{
-						// printf("bigger - divs1: %d, divs2: %d\n", divs1, divs2);
-						continue;
-					}
-
-					copy_v3_v3(src_p1, &scs->points_old[mface[i].v1 * 3]);
-					copy_v3_v3(dest_p1, &scs->points[mface[i].v1 * 3]);
-					if(again == 1 && mface[i].v4)
-					{
-						copy_v3_v3(src_p2, &scs->points_old[mface[i].v3 * 3]);
-						copy_v3_v3(dest_p2, &scs->points[mface[i].v3 * 3]);
-
-						copy_v3_v3(src_p3,&scs->points_old[mface[i].v4 * 3]);
-						copy_v3_v3(dest_p3, &scs->points[mface[i].v4 * 3]);
-					}
-					else {
-						copy_v3_v3(src_p2, &scs->points_old[mface[i].v2 * 3]);
-						copy_v3_v3(dest_p2, &scs->points[mface[i].v2 * 3]);
-						copy_v3_v3(src_p3, &scs->points_old[mface[i].v3 * 3]);
-						copy_v3_v3(dest_p3, &scs->points[mface[i].v3 * 3]);
-					}
-
-					mul_v3_fl(src_p1, (1.0-uf-vf));
-					mul_v3_fl(dest_p1, (1.0-uf-vf));
-
-					mul_v3_fl(src_p2, uf);
-					mul_v3_fl(dest_p2, uf);
-
-					mul_v3_fl(src_p3, vf);
-					mul_v3_fl(dest_p3, vf);
-
-					add_v3_v3v3(src_p, src_p1, src_p2);
-					add_v3_v3v3(dest_p, dest_p1, dest_p2);
-
-					add_v3_v3(src_p, src_p3);
-					add_v3_v3(dest_p, dest_p3);
-
-					if(newdivs > divs)
-						printf("mem problem\n");
-
-					// mMovPoints.push_back(p + trinorm);
-					add_v3_v3v3(src_tmpvec, src_p, src_trinorm);
-					add_v3_v3v3(dest_tmpvec, dest_p, dest_trinorm);
-
-					// mul_m4_v3(ob->obmat, tmpvec); // DG: use local coordinates, we save MAT anyway
-					copy_v3_v3(&scs->points_old[3 * (dm->getNumVerts(dm) + newdivs)], src_tmpvec);
-					copy_v3_v3(&scs->points[3 * (dm->getNumVerts(dm) + newdivs)], dest_tmpvec);
-					newdivs++;
-
-					if(newdivs > divs)
-						printf("mem problem\n");
-
-					// mMovPoints.push_back(p - trinorm);
-					copy_v3_v3(src_tmpvec, src_p);
-					copy_v3_v3(dest_tmpvec, dest_p);
-
-					sub_v3_v3(src_tmpvec, src_trinorm);
-					sub_v3_v3(dest_tmpvec, dest_trinorm);
-
-					// mul_m4_v3(ob->obmat, tmpvec); // DG: use local coordinates, we save MAT anyway
-					copy_v3_v3(&scs->points_old[3 * (dm->getNumVerts(dm) + newdivs)], src_tmpvec);
-					copy_v3_v3(&scs->points[3 * (dm->getNumVerts(dm) + newdivs)], dest_tmpvec);
-					newdivs++;
-				}
-			}
-
-			if(again == 0 && mface[i].v4)
-				again++;
-			else
-				again = 0;
-
-			facecounter++;
-
-		} while(again!=0);
-	}
-
-	// scs->numpoints = dm->getNumVerts(dm) + newdivs;
-
-}
-
-/*! init triangle divisions */
-static void calcTriangleDivs(Object *ob, MVert *verts, int UNUSED(numverts), MFace *faces, int numfaces, int numtris, int **tridivs, float cell_len)
-{
-	// mTriangleDivs1.resize( faces.size() );
-	// mTriangleDivs2.resize( faces.size() );
-	// mTriangleDivs3.resize( faces.size() );
-
-	size_t i = 0, facecounter = 0;
-	float maxscale[3] = {1,1,1}; // = channelFindMaxVf(mcScale); get max scale value
-	float maxpart = ABS(maxscale[0]);
-	float scaleFac = 0;
-	float fsTri = 0;
-	if(ABS(maxscale[1])>maxpart) maxpart = ABS(maxscale[1]);
-	if(ABS(maxscale[2])>maxpart) maxpart = ABS(maxscale[2]);
-	scaleFac = 1.0 / maxpart;
-	// featureSize = mLevel[mMaxRefine].nodeSize
-	fsTri = cell_len * 0.75 * scaleFac; // fsTri = cell_len * 0.9;
-
-	if(*tridivs)
-		MEM_freeN(*tridivs);
-
-	*tridivs = MEM_callocN(sizeof(int) * numtris * 3, "Smoke_Tridivs");
-
-	for(i = 0, facecounter = 0; i < numfaces; i++) 
-	{
-		float p0[3], p1[3], p2[3];
-		float side1[3];
-		float side2[3];
-		float side3[3];
-		int divs1=0, divs2=0, divs3=0;
-
-		copy_v3_v3(p0, verts[faces[i].v1].co);
-		mul_m4_v3(ob->obmat, p0);
-		copy_v3_v3(p1, verts[faces[i].v2].co);
-		mul_m4_v3(ob->obmat, p1);
-		copy_v3_v3(p2, verts[faces[i].v3].co);
-		mul_m4_v3(ob->obmat, p2);
-
-		sub_v3_v3v3(side1, p1, p0);
-		sub_v3_v3v3(side2, p2, p0);
-		sub_v3_v3v3(side3, p1, p2);
-
-		if(dot_v3v3(side1, side1) > fsTri*fsTri)
-		{ 
-			float tmp = normalize_v3(side1);
-			divs1 = (int)ceil(tmp/fsTri); 
-		}
-		if(dot_v3v3(side2, side2) > fsTri*fsTri)
-		{ 
-			float tmp = normalize_v3(side2);
-			divs2 = (int)ceil(tmp/fsTri); 
-			
-			/*		
-			// debug
-			if(i==0)
-				printf("b tmp: %f, fsTri: %f, divs2: %d\n", tmp, fsTri, divs2);
-			*/
-			
-		}
-
-		(*tridivs)[3 * facecounter + 0] = divs1;
-		(*tridivs)[3 * facecounter + 1] = divs2;
-		(*tridivs)[3 * facecounter + 2] = divs3;
-
-		// TODO quad case
-		if(faces[i].v4)
-		{
-			divs1=0, divs2=0, divs3=0;
-
-			facecounter++;
-			
-			copy_v3_v3(p0, verts[faces[i].v3].co);
-			mul_m4_v3(ob->obmat, p0);
-			copy_v3_v3(p1, verts[faces[i].v4].co);
-			mul_m4_v3(ob->obmat, p1);
-			copy_v3_v3(p2, verts[faces[i].v1].co);
-			mul_m4_v3(ob->obmat, p2);
-
-			sub_v3_v3v3(side1, p1, p0);
-			sub_v3_v3v3(side2, p2, p0);
-			sub_v3_v3v3(side3, p1, p2);
-
-			if(dot_v3v3(side1, side1) > fsTri*fsTri)
-			{ 
-				float tmp = normalize_v3(side1);
-				divs1 = (int)ceil(tmp/fsTri); 
-			}
-			if(dot_v3v3(side2, side2) > fsTri*fsTri)
-			{ 
-				float tmp = normalize_v3(side2);
-				divs2 = (int)ceil(tmp/fsTri); 
-			}
-
-			(*tridivs)[3 * facecounter + 0] = divs1;
-			(*tridivs)[3 * facecounter + 1] = divs2;
-			(*tridivs)[3 * facecounter + 2] = divs3;
-		}
-		facecounter++;
-	}
-}
-
 #endif /* WITH_SMOKE */
 
 static void smokeModifier_freeDomain(SmokeModifierData *smd)
@@ -769,36 +342,18 @@ static void smokeModifier_freeCollision(SmokeModifierData *smd)
 	{
 		SmokeCollSettings *scs = smd->coll;
 
-		if(scs->numpoints)
+		if(scs->numverts)
 		{
-			if(scs->points)
+			if(scs->verts_old)
 			{
-				MEM_freeN(scs->points);
+				MEM_freeN(scs->verts_old);
-				scs->points = NULL;
+				scs->verts_old = NULL;
-			}
-			if(scs->points_old)
-			{
-				MEM_freeN(scs->points_old);
-				scs->points_old = NULL;
-			}
-			if(scs->tridivs)
-			{
-				MEM_freeN(scs->tridivs);
-				scs->tridivs = NULL;
 			}
 		}
 
-		if(scs->bvhtree)
-		{
-			BLI_bvhtree_free(scs->bvhtree);
-			scs->bvhtree = NULL;
-		}
-
-#ifdef USE_SMOKE_COLLISION_DM
 		if(smd->coll->dm)
 			smd->coll->dm->release(smd->coll->dm);
 		smd->coll->dm = NULL;
-#endif
 
 		MEM_freeN(smd->coll);
 		smd->coll = NULL;
@@ -851,21 +406,10 @@ void smokeModifier_reset(struct SmokeModifierData *smd)
 		{
 			SmokeCollSettings *scs = smd->coll;
 
-			if(scs->numpoints && scs->points)
+			if(scs->numverts && scs->verts_old)
 			{
-				MEM_freeN(scs->points);
+				MEM_freeN(scs->verts_old);
-				scs->points = NULL;
+				scs->verts_old = NULL;
-			
-				if(scs->points_old)
-				{
-					MEM_freeN(scs->points_old);
-					scs->points_old = NULL;
-				}
-				if(scs->tridivs)
-				{
-					MEM_freeN(scs->tridivs);
-					scs->tridivs = NULL;
-				}
 			}
 		}
 	}
@@ -950,15 +494,10 @@ void smokeModifier_createType(struct SmokeModifierData *smd)
 			smd->coll = MEM_callocN(sizeof(SmokeCollSettings), "SmokeColl");
 
 			smd->coll->smd = smd;
-			smd->coll->points = NULL;
+			smd->coll->verts_old = NULL;
-			smd->coll->points_old = NULL;
+			smd->coll->numverts = 0;
-			smd->coll->tridivs = NULL;
-			smd->coll->vel = NULL;
-			smd->coll->numpoints = 0;
-			smd->coll->numtris = 0;
-			smd->coll->bvhtree = NULL;
 			smd->coll->type = 0; // static obstacle
-			smd->coll->dx = 1.0f / 50.0f;
+			smd->coll->dm = NULL;
 
 #ifdef USE_SMOKE_COLLISION_DM
 			smd->coll->dm = NULL;
@@ -1078,6 +617,153 @@ static void smoke_calc_domain(Scene *UNUSED(scene), Object *UNUSED(ob), SmokeMod
 #endif
 }
 
+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;
+	{
+		DerivedMesh *dm = NULL;
+		MDeformVert *dvert = NULL;
+		MVert *mvert = NULL;
+		MFace *mface = NULL;
+		BVHTreeFromMesh treeData = {0};
+		int numverts, i, z;
+		int *res = sds->res;
+
+		float surface_distance = 0.6;
+
+		float *vert_vel = NULL;
+		int has_velocity = 0;
+
+		tstart();
+
+		dm = CDDM_copy(scs->dm);
+		CDDM_calc_normals(dm);
+		mvert = dm->getVertArray(dm);
+		mface = dm->getTessFaceArray(dm);
+		numverts = dm->getNumVerts(dm);
+		dvert = dm->getVertDataArray(dm, CD_MDEFORMVERT);
+
+		// DG TODO
+		// if(scs->type > SM_COLL_STATIC)
+		// if line above is used, the code is in trouble if the object moves but is declared as "does not move"
+
+		// if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY) 
+		{
+			vert_vel = MEM_callocN(sizeof(float) * numverts * 3, "smoke_obs_velocity");
+
+			if (scs->numverts != numverts || !scs->verts_old) {
+				if (scs->verts_old) MEM_freeN(scs->verts_old);
+
+				scs->verts_old = MEM_callocN(sizeof(float) * numverts * 3, "smoke_obs_verts_old");
+				scs->numverts = numverts;
+			}
+			else {
+				has_velocity = 1;
+			}
+		}
+
+		/*	Transform collider vertices to
+		*   domain grid space for fast lookups */
+		for (i = 0; i < numverts; i++) {
+			float n[3];
+
+			/* vert pos */
+			mul_m4_v3(coll_ob->obmat, mvert[i].co);
+			sub_v3_v3(mvert[i].co, sds->p0);
+			mul_v3_fl(mvert[i].co, (1.0f / sds->dx) / sds->scale);
+
+			/* vert normal */
+			normal_short_to_float_v3(n, mvert[i].no);
+			mul_mat3_m4_v3(coll_ob->obmat, n);
+			normalize_v3(n);
+			normal_float_to_short_v3(mvert[i].no, n);
+
+			/* vert velocity */
+			// DG TODO 
+			// if (sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY) 
+			{
+				if (has_velocity) 
+				{
+					sub_v3_v3v3(&vert_vel[i*3], mvert[i].co, &scs->verts_old[i*3]);
+					mul_v3_fl(&vert_vel[i*3], sds->dx/dt);
+				}
+				copy_v3_v3(&scs->verts_old[i*3], mvert[i].co);
+			}
+		}
+
+		if (bvhtree_from_mesh_faces(&treeData, dm, 0.0f, 4, 6)) {
+			#pragma omp parallel for schedule(static)
+			for (z = 0; z < res[2]; z++) {
+				int x,y;
+				for (x = 0; x < res[0]; x++)
+				for (y = 0; y < res[1]; y++) {
+					int index = x + y*res[0] + z*res[0]*res[1];
+
+					float ray_start[3] = {(float)x + 0.5f, (float)y + 0.5f, (float)z + 0.5f};
+					float ray_dir[3] = {1.0f, 0.0f, 0.0f};
+
+					BVHTreeRayHit hit = {0};
+					BVHTreeNearest nearest = {0};
+
+					hit.index = -1;
+					hit.dist = 9999;
+					nearest.index = -1;
+					nearest.dist = surface_distance * surface_distance; /* find_nearest uses squared distance */
+
+					/* find the nearest point on the mesh */
+					if (BLI_bvhtree_find_nearest(treeData.tree, ray_start, &nearest, treeData.nearest_callback, &treeData) != -1) {
+						float weights[4];
+						int v1, v2, v3, f_index = nearest.index;
+						float n1[3], n2[3], n3[3], hit_normal[3];
+
+						/* calculate barycentric weights for nearest point */
+						v1 = mface[f_index].v1;
+						v2 = (nearest.flags & BVH_ONQUAD) ? mface[f_index].v3 : mface[f_index].v2;
+						v3 = (nearest.flags & BVH_ONQUAD) ? mface[f_index].v4 : mface[f_index].v3;
+						interp_weights_face_v3(weights, mvert[v1].co, mvert[v2].co, mvert[v3].co, NULL, nearest.co);
+
+						// DG TODO
+						// if(sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY) 
+						if(has_velocity)
+						{
+							/* interpolate vertex normal vectors to get nearest point normal */
+							normal_short_to_float_v3(n1, mvert[v1].no);
+							normal_short_to_float_v3(n2, mvert[v2].no);
+							normal_short_to_float_v3(n3, mvert[v3].no);
+							interp_v3_v3v3v3(hit_normal, n1, n2, n3, weights);
+							normalize_v3(hit_normal);
+							
+							/* apply object velocity */
+							{
+								float hit_vel[3];
+								interp_v3_v3v3v3(hit_vel, &vert_vel[v1*3], &vert_vel[v2*3], &vert_vel[v3*3], weights);
+								velocityX[index] += hit_vel[0];
+								velocityY[index] += hit_vel[1];
+								velocityZ[index] += hit_vel[2];
+							}
+						}
+
+						/* tag obstacle cells */
+						obstacle_map[index] = 1;
+
+						if(has_velocity)
+							obstacle_map[index] |= 8;
+					}
+				}
+			}
+		}
+		/* free bvh tree */
+		free_bvhtree_from_mesh(&treeData);
+		dm->release(dm);
+
+		if (vert_vel) MEM_freeN(vert_vel);
+
+		tend();
+		printf ( "Time: %f\n\n", ( float ) tval() );
+
+	}
+}
+
 /* Animated obstacles: dx_step = ((x_new - x_old) / totalsteps) * substep */
 static void update_obstacles(Scene *scene, Object *ob, SmokeDomainSettings *sds, float dt, int substep, int totalsteps)
 {
@@ -1092,7 +778,7 @@ static void update_obstacles(Scene *scene, Object *ob, SmokeDomainSettings *sds,
 	float *velxOrig = smoke_get_velocity_x(sds->fluid);
 	float *velyOrig = smoke_get_velocity_y(sds->fluid);
 	float *velzOrig = smoke_get_velocity_z(sds->fluid);
-	// float *density = smoke_get_density(sds->fluid);
+	float *density = smoke_get_density(sds->fluid);
 	unsigned int z;
 
 	smoke_get_ob_velocity(sds->fluid, &velx, &vely, &velz);
@@ -1100,15 +786,6 @@ static void update_obstacles(Scene *scene, Object *ob, SmokeDomainSettings *sds,
 	// TODO: delete old obstacle flags
 	for(z = 0; z < sds->res[0] * sds->res[1] * sds->res[2]; z++)
 	{
-		if(obstacles[z])
-		{
-			// density[z] = 0;
-
-			velxOrig[z] = 0;
-			velyOrig[z] = 0;
-			velzOrig[z] = 0;
-		}
-
 		if(obstacles[z] & 8) // Do not delete static obstacles
 		{
 			obstacles[z] = 0;
@@ -1119,6 +796,7 @@ static void update_obstacles(Scene *scene, Object *ob, SmokeDomainSettings *sds,
 		velz[z] = 0;
 	}
 
+
 	collobjs = get_collisionobjects(scene, ob, sds->coll_group, &numcollobj, eModifierType_Smoke);
 
 	// update obstacle tags in cells
@@ -1129,94 +807,29 @@ static void update_obstacles(Scene *scene, Object *ob, SmokeDomainSettings *sds,
 
 		// DG TODO: check if modifier is active?
 		
-		if((smd2->type & MOD_SMOKE_TYPE_COLL) && smd2->coll && smd2->coll->points && smd2->coll->points_old)
+		if((smd2->type & MOD_SMOKE_TYPE_COLL) && smd2->coll)
 		{
 			SmokeCollSettings *scs = smd2->coll;
-			unsigned int i;
 			
-			/*
+			obstacles_from_derivedmesh(collob, sds, scs, obstacles, velx, vely, velz, dt);
-			// DG TODO: support static cobstacles, but basicly we could even support static + rigid with one set of code
-			if(scs->type > SM_COLL_STATIC)
-			*/
-
-			/* Handle collisions */
-			for(i = 0; i < scs->numpoints; i++)
-			{
-				// 1. get corresponding cell
-				int cell[3];
-				float pos[3], oldpos[3], vel[3];
-				float cPos[3], cOldpos[3]; /* current position in substeps */
-				int badcell = 0;
-				size_t index;
-				int j;
-
-				// translate local points into global positions
-				copy_v3_v3(cPos, &scs->points[3 * i]);
-				mul_m4_v3(scs->mat, cPos);
-				copy_v3_v3(pos, cPos);
-
-				copy_v3_v3(cOldpos, &scs->points_old[3 * i]);
-				mul_m4_v3(scs->mat_old, cOldpos);
-				copy_v3_v3(oldpos, cOldpos);
-
-				/* support for rigid bodies, armatures etc */
-				{
-					float tmp[3];
-
-					/* x_current = x_old + (x_new - x_old) * step_current / steps_total */
-					float mulStep = (float)(((float)substep) / ((float)totalsteps));
-
-					sub_v3_v3v3(tmp, cPos, cOldpos);
-					mul_v3_fl(tmp, mulStep);
-					add_v3_v3(cOldpos, tmp);
-				}
-
-				sub_v3_v3v3(vel, pos, oldpos);
-				/* Scale velocity to incorperate the object movement during this step */
-				mul_v3_fl(vel, 1.0 / (totalsteps * dt * sds->scale));
-				// mul_v3_fl(vel, 1.0 / dt);
-
-				// DG TODO: cap velocity to maxVelMag (or maxvel)
-
-				// oldpos + velocity * dt = newpos
-				get_cell(sds->p0, sds->res, sds->dx*sds->scale, cOldpos /* use current position here instead of "pos" */, cell, 0);
-
-				// check if cell is valid (in the domain boundary)
-				for(j = 0; j < 3; j++)
-					if((cell[j] > sds->res[j] - 1) || (cell[j] < 0))
-					{
-						badcell = 1;
-						break;
-					}
-														
-				if(badcell)
-					continue;
-
-				// 2. set cell values (heat, density and velocity)
-				index = smoke_get_index(cell[0], sds->res[0], cell[1], sds->res[1], cell[2]);
-
-				// Don't overwrite existing obstacles
-				if(obstacles[index])
-					continue;
-												
-				// printf("cell[0]: %d, cell[1]: %d, cell[2]: %d\n", cell[0], cell[1], cell[2]);								
-				// printf("res[0]: %d, res[1]: %d, res[2]: %d, index: %d\n\n", sds->res[0], sds->res[1], sds->res[2], index);																	
-				obstacles[index] = 1 | 8 /* ANIMATED */;
-
-				if(len_v3(vel) > FLT_EPSILON)
-				{
-					// Collision object is moving
-
-					velx[index] = vel[0]; // use "+="?
-					vely[index] = vel[1];
-					velz[index] = vel[2];
-				}
-			}
 		}
 	}
 
 	if(collobjs)
 		MEM_freeN(collobjs);
+
+	/* obstacle cells should not contain any velocity from the smoke simulation */
+	for(z = 0; z < sds->res[0] * sds->res[1] * sds->res[2]; z++)
+	{
+		if(obstacles[z])
+		{
+			velxOrig[z] = 0;
+			velyOrig[z] = 0;
+			velzOrig[z] = 0;
+
+			density[z] = 0;
+		}
+	}
 }
 
 static void update_flowsfluids(Scene *scene, Object *ob, SmokeDomainSettings *sds, float time)
@@ -1656,88 +1269,22 @@ void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedM
 	}
 	else if(smd->type & MOD_SMOKE_TYPE_COLL)
 	{
-		/* Check if domain resolution changed */
-		/* DG TODO: can this be solved more elegant using dependancy graph? */
-		{
-			SmokeCollSettings *scs = smd->coll;
-			Base *base = scene->base.first;
-			int changed = 0;
-			float dx = FLT_MAX;
-			float scale = 1.0f;
-			int haveDomain = 0;
-
-			for ( ; base; base = base->next) 
-			{
-				SmokeModifierData *smd2 = (SmokeModifierData *)modifiers_findByType(base->object, eModifierType_Smoke);
-
-				if (smd2 && (smd2->type & MOD_SMOKE_TYPE_DOMAIN) && smd2->domain)
-				{
-					SmokeDomainSettings *sds = smd2->domain;
-
-					if(sds->dx * sds->scale < dx)
-					{
-						dx = sds->dx;
-						scale = sds->scale;
-						changed = 1;
-					}
-
-					haveDomain = 1;
-				}
-			}
-
-			if(!haveDomain)
-				return;
-			
-			if(changed)
-			{
-				if(dx*scale != scs->dx)
-				{
-					scs->dx = dx*scale;
-					smokeModifier_reset(smd);
-				}
-			}
-		}
-
 		if(scene->r.cfra >= smd->time)
 			smokeModifier_init(smd, ob, scene, dm);
 
+		if(smd->coll)
+		{
+			if (smd->coll->dm) 
+				smd->coll->dm->release(smd->coll->dm);
+
+			smd->coll->dm = CDDM_copy(dm);
+			DM_ensure_tessface(smd->coll->dm);
+		}
+
 		if(scene->r.cfra > smd->time)
 		{
-			unsigned int i;
 			SmokeCollSettings *scs = smd->coll;
-			float *points_old = scs->points_old;
-			float *points = scs->points;
-			unsigned int numpoints = scs->numpoints; 
-
-			// XXX TODO <-- DG: what is TODO here?
 			smd->time = scene->r.cfra;
-
-			// rigid movement support
-			copy_m4_m4(scs->mat_old, scs->mat);
-			copy_m4_m4(scs->mat, ob->obmat);
-
-			if(scs->type != SM_COLL_ANIMATED) // if(not_animated)
-			{
-				// nothing to do, "mat" is already up to date
-			}
-			else
-			{
-				// XXX TODO: need to update positions + divs
-
-				if(scs->numverts != dm->getNumVerts(dm))
-				{
-					// DG TODO: reset modifier?
-					return;
-				}
-
-				for(i = 0; i < numpoints * 3; i++)
-				{
-					points_old[i] = points[i];
-				}
-
-				DM_ensure_tessface(dm);
-				fill_scs_points_anim(ob, dm, scs);
-			}
 		}
 		else if(scene->r.cfra < smd->time)
 		{
@@ -1803,7 +1350,7 @@ void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedM
 		if (framenr != scene->r.cfra)
 			return;
 
-		tstart();
+		// tstart();
 
 		smoke_calc_domain(scene, ob, smd);
 
@@ -1855,7 +1402,7 @@ void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedM
 		if(framenr != startframe)
 			BKE_ptcache_write(&pid, framenr);
 
-		tend();
+		//tend();
 		// printf ( "Frame: %d, Time: %f\n\n", (int)smd->time, (float) tval() );
 	}
 }

source/blender/blenloader/intern/readfile.c

@@ -4402,8 +4402,9 @@ static void direct_link_modifiers(FileData *fd, ListBase *lb)
 				smd->coll = newdataadr(fd, smd->coll);
 				smd->coll->smd = smd;
 				if (smd->coll) {
-					smd->coll->points = NULL;
+					smd->coll->verts_old = NULL;
-					smd->coll->numpoints = 0;
+					smd->coll->numverts = 0;
+					smd->coll->dm = NULL;
 				}
 				else
 					smd->type = 0;

source/blender/editors/space_view3d/drawobject.c

@@ -6996,6 +6996,8 @@ void draw_object(Scene *scene, ARegion *ar, View3D *v3d, Base *base, const short
 				            smd->domain->res, smd->domain->dx,
 				            smd->domain->tex_shadow);
 				GPU_free_smoke(smd);
+				// draw_smoke_heat(smd->domain);
+				// draw_smoke_velocity(smd->domain);
 // #endif
 #if 0
 				int x, y, z;

source/blender/editors/space_view3d/drawvolume.c

@@ -36,6 +36,7 @@
 
 #include "DNA_scene_types.h"
 #include "DNA_screen_types.h"
+#include "DNA_smoke_types.h"
 #include "DNA_view3d_types.h"
 
 #include "BLI_utildefines.h"
@@ -454,3 +455,83 @@ void draw_volume(ARegion *ar, GPUTexture *tex, float min[3], float max[3], int r
 		glDepthMask(GL_TRUE);	
 	}
 }
+
+void draw_smoke_velocity(SmokeDomainSettings *domain)
+{
+	float x,y,z;
+	int *res = domain->res;
+	float *vel_x = smoke_get_velocity_x(domain->fluid);
+	float *vel_y = smoke_get_velocity_y(domain->fluid);
+	float *vel_z = smoke_get_velocity_z(domain->fluid);
+	float *density = smoke_get_density(domain->fluid);
+
+	float *min = domain->p0;
+	float cell_size = domain->dx * domain->scale;
+	float step_size = ((float)MAX3(res[0], res[1], res[2]))/16.f;
+	float vf = domain->scale / 16.f * 2.f;
+
+	glColor3f(1.0f, 0.5f, 0.0f);
+	glLineWidth(1.0f);
+
+	for (x=0; x<res[0]; x+=step_size)
+		for (y=0; y<res[1]; y+=step_size)
+			for (z=0; z<res[2]; z+=step_size) {
+				int index = floor(x) + floor(y)*res[0] + floor(z)*res[0]*res[1];
+
+				float pos[3] = {min[0]+((float)x + 0.5f)*cell_size, min[1]+((float)y + 0.5f)*cell_size, min[2]+((float)z + 0.5f)*cell_size};
+				float vel = sqrtf(vel_x[index]*vel_x[index] + vel_y[index]*vel_y[index] + vel_z[index]*vel_z[index]);
+
+				if (vel > 0.01f) {
+					float col_g = 1.0f - vel;
+					CLAMP(col_g, 0.0f, 1.0f);
+					glColor3f(1.0f, col_g, 0.0f);
+					glPointSize(10.0f * vel);
+
+					glBegin(GL_LINES);
+					glVertex3f(pos[0], pos[1], pos[2]);
+					glVertex3f(pos[0]+vel_x[index]*vf, pos[1]+vel_y[index]*vf, pos[2]+vel_z[index]*vf);
+					glEnd();
+					glBegin(GL_POINTS);
+					glVertex3f(pos[0]+vel_x[index]*vf, pos[1]+vel_y[index]*vf, pos[2]+vel_z[index]*vf);
+					glEnd();
+				}
+	}
+}
+
+void draw_smoke_heat(SmokeDomainSettings *domain)
+{
+	float x,y,z;
+	int *res = domain->res;
+	// float *heat = smoke_get_heat(domain->fluid);
+	unsigned char *heat = smoke_get_obstacle(domain->fluid);
+
+	float *min = domain->p0;
+	float cell_size = domain->dx * domain->scale;
+	//float step_size = ((float)MAX3(res[0], res[1], res[2]))/16.f;
+	float step_size = 1.0f;
+
+	for (x=0; x<res[0]; x+=step_size)
+		for (y=0; y<res[1]; y+=step_size)
+			for (z=0; z<res[2]; z+=step_size) {
+				int index = floor(x) + floor(y)*res[0] + floor(z)*res[0]*res[1];
+
+				float pos[3] = {min[0]+((float)x + 0.5f)*cell_size, min[1]+((float)y + 0.5f)*cell_size, min[2]+((float)z + 0.5f)*cell_size};
+
+				if (heat[index] > 0) {
+					GLUquadricObj *qobj = gluNewQuadric();
+					gluQuadricDrawStyle(qobj, GLU_FILL);
+
+					glColor3f(0, 0, 1.0f);
+					
+					glPushMatrix();
+
+					glTranslatef(pos[0], pos[1], pos[2]);
+					glScalef(domain->dx * domain->scale * 0.1, domain->dx * domain->scale * 0.1, domain->dx * domain->scale * 0.1);
+					gluSphere(qobj, 1.0, 8, 5);	
+
+					glPopMatrix();
+	
+					gluDeleteQuadric(qobj);
+				}
+	}
+}

source/blender/editors/space_view3d/view3d_intern.h

@@ -213,6 +213,8 @@ extern const char *view3d_context_dir[]; /* doc access */
 
 /* draw_volume.c */
 void draw_volume(struct ARegion *ar, struct GPUTexture *tex, float min[3], float max[3], int res[3], float dx, struct GPUTexture *tex_shadow);
+void draw_smoke_velocity(struct SmokeDomainSettings *domain);
+void draw_smoke_heat(struct SmokeDomainSettings *domain);
 
 /* workaround for trivial but noticeable camera bug caused by imprecision
  * between view border calculation in 2D/3D space, workaround for bug [#28037].

source/blender/makesdna/DNA_smoke_types.h

@@ -131,28 +131,14 @@ typedef struct SmokeFlowSettings {
 	int flags; /* absolute emission etc*/
 } SmokeFlowSettings;
 
-
-// struct BVHTreeFromMesh *bvh;
-// float mat[4][4];
-// float mat_old[4][4];
-
 /* collision objects (filled with smoke) */
 typedef struct SmokeCollSettings {
 	struct SmokeModifierData *smd; /* for fast RNA access */
-	struct BVHTree *bvhtree; /* bounding volume hierarchy for this cloth object */
+	struct DerivedMesh *dm;
-	float *points;
+	float *verts_old;
-	float *points_old;
+	int numverts;
-	float *vel; // UNUSED
-	int *tridivs;
-	float mat[4][4];
-	float mat_old[4][4];
-	int numpoints;
-	int numverts; // check if mesh changed
-	int numtris;
-	float dx; /* global domain cell length taken from (scale / resolution) */
 	short type; // static = 0, rigid = 1, dynamic = 2
 	short pad;
-	int pad2;
 } SmokeCollSettings;
 
 #endif

source/blender/makesrna/intern/rna_smoke.c

@@ -113,7 +113,7 @@ static char *rna_SmokeCollSettings_path(PointerRNA *ptr)
 	return BLI_sprintfN("modifiers[\"%s\"].coll_settings", md->name);
 }
 
-static int rna_SmokeModifier_density_get_length(PointerRNA *ptr, int length[RNA_MAX_ARRAY_DIMENSION])
+int rna_SmokeModifier_density_get_length(PointerRNA *ptr, int length[3])
 {
 	SmokeDomainSettings *settings = (SmokeDomainSettings *)ptr->data;
 
@@ -133,7 +133,7 @@ static int rna_SmokeModifier_density_get_length(PointerRNA *ptr, int length[RNA_
 	return length[0];
 }
 
-static void rna_SmokeModifier_density_get(PointerRNA *ptr, float *values)
+void rna_SmokeModifier_density_get(PointerRNA *ptr, float *values)
 {
 	SmokeDomainSettings *settings = (SmokeDomainSettings *)ptr->data;
 	float *density = smoke_get_density(settings->fluid);