smoke: decoupling high-res, introducing reset for heat+gravity, fixing another obstacle problem

4 files changed, 98 insertions, 49 deletions

diff --git a/intern/smoke/intern/FLUID_3D.cpp b/intern/smoke/intern/FLUID_3D.cpp
index 344883866a8..198d08ee369 100644
--- a/intern/smoke/intern/FLUID_3D.cpp
+++ b/intern/smoke/intern/FLUID_3D.cpp
@@ -59,7 +59,10 @@ FLUID_3D::FLUID_3D(int *res, int amplify, float *p0, float dt) :
 	_maxRes = MAX3(_xRes, _yRes, _zRes);
 	
 	// initialize wavelet turbulence
-	_wTurbulence = new WTURBULENCE(_res[0],_res[1],_res[2], amplify);
+	if(amplify)
+		_wTurbulence = new WTURBULENCE(_res[0],_res[1],_res[2], amplify);
+	else
+		_wTurbulence = NULL;
 	
 	// scale the constants according to the refinement of the grid
 	_dx = 1.0f / (float)_maxRes;
@@ -215,11 +218,11 @@ void FLUID_3D::step()
 	wipeBoundaries();
 
 	// run the solvers
-  addVorticity();
-  addBuoyancy(_heat, _density);
+	addVorticity();
+	addBuoyancy(_heat, _density);
 	addForce();
 	project();
-  diffuseHeat();
+	diffuseHeat();
 
 	// advect everything
 	advectMacCormack();
@@ -249,6 +252,19 @@ void FLUID_3D::step()
   */
 	_totalTime += _dt;
 	_totalSteps++;
+
+	// clear obstacles
+	for (int i = 0; i < _totalCells; i++)
+	{
+		if(_obstacles[i])
+		{
+			_xVelocity[i] =
+			_yVelocity[i] =
+			_zVelocity[i] = 0.0f;
+			_pressure[i] = 0.0f;
+		}
+		_obstacles[i] = 0;
+	}
 }
 
 //////////////////////////////////////////////////////////////////////
@@ -354,6 +370,7 @@ void FLUID_3D::addForce()
 void FLUID_3D::project()
 {
 	int index, x, y, z;
+
 	setObstacleBoundaries();
 
 	// copy out the boundaries
@@ -397,6 +414,8 @@ void FLUID_3D::project()
 	// solve Poisson equation
 	solvePressure(_pressure, _divergence, _obstacles);
 
+	setObstaclePressure();
+
 	// project out solution
 	float invDx = 1.0f / _dx;
 	index = _slabSize + _xRes + 1;
@@ -404,9 +423,12 @@ void FLUID_3D::project()
 		for (y = 1; y < _yRes - 1; y++, index += 2)
 			for (x = 1; x < _xRes - 1; x++, index++)
 			{
-				_xVelocity[index] -= 0.5f * (_pressure[index + 1]     - _pressure[index - 1])     * invDx;
-				_yVelocity[index] -= 0.5f * (_pressure[index + _xRes]  - _pressure[index - _xRes]) * invDx;
-				_zVelocity[index] -= 0.5f * (_pressure[index + _slabSize] - _pressure[index - _slabSize]) * invDx;
+				if(!_obstacles[index])
+				{
+					_xVelocity[index] -= 0.5f * (_pressure[index + 1]     - _pressure[index - 1]) * invDx;
+					_yVelocity[index] -= 0.5f * (_pressure[index + _xRes]  - _pressure[index - _xRes]) * invDx;
+					_zVelocity[index] -= 0.5f * (_pressure[index + _slabSize] - _pressure[index - _slabSize]) * invDx;
+				}
 			}
 }
 
@@ -443,34 +465,8 @@ void FLUID_3D::addObstacle(OBSTACLE* obstacle)
 //////////////////////////////////////////////////////////////////////
 // calculate the obstacle directional types
 //////////////////////////////////////////////////////////////////////
-void FLUID_3D::setObstacleBoundaries()
+void FLUID_3D::setObstaclePressure()
 {
-	// cull degenerate obstacles , move to addObstacle?
-	for (int z = 1, index = _slabSize + _xRes + 1;
-			z < _zRes - 1; z++, index += 2 * _xRes)
-		for (int y = 1; y < _yRes - 1; y++, index += 2)
-			for (int x = 1; x < _xRes - 1; x++, index++)
-				if (_obstacles[index] != EMPTY)
-				{
-					const int top   = _obstacles[index + _slabSize];
-					const int bottom= _obstacles[index - _slabSize];
-					const int up    = _obstacles[index + _xRes];
-					const int down  = _obstacles[index - _xRes];
-					const int left  = _obstacles[index - 1];
-					const int right = _obstacles[index + 1];
-
-					int counter = 0;
-					if (up)    counter++;
-					if (down)  counter++;
-					if (left)  counter++;
-					if (right) counter++;
-					if (top)  counter++;
-					if (bottom) counter++;
-
-					if (counter < 3)
-						_obstacles[index] = EMPTY;
-				}
-
 	// tag remaining obstacle blocks
 	for (int z = 1, index = _slabSize + _xRes + 1;
 			z < _zRes - 1; z++, index += 2 * _xRes)
@@ -478,7 +474,7 @@ void FLUID_3D::setObstacleBoundaries()
 			for (int x = 1; x < _xRes - 1; x++, index++)
 		{
 			// could do cascade of ifs, but they are a pain
-			if (_obstacles[index] != EMPTY)
+			if (_obstacles[index])
 			{
 				const int top   = _obstacles[index + _slabSize];
 				const int bottom= _obstacles[index - _slabSize];
@@ -498,7 +494,7 @@ void FLUID_3D::setObstacleBoundaries()
 				_pressure[index] = 0.0f;
 
 				// average pressure neighbors
-				float pcnt = 0.;
+				float pcnt = 0., vp = 0.;
 				if (left && !right) {
 					_pressure[index] += _pressure[index + 1];
 					pcnt += 1.;
@@ -516,14 +512,24 @@ void FLUID_3D::setObstacleBoundaries()
 					pcnt += 1.;
 				}
 				if (top && !bottom) {
-					_pressure[index] += _pressure[index - _xRes];
+					_pressure[index] += _pressure[index - _slabSize];
 					pcnt += 1.;
+					// _zVelocity[index] +=  - _zVelocity[index - _slabSize];
+					// vp += 1.0;
 				}
 				if (!top && bottom) {
-					_pressure[index] += _pressure[index + _xRes];
+					_pressure[index] += _pressure[index + _slabSize];
 					pcnt += 1.;
+					// _zVelocity[index] +=  - _zVelocity[index + _slabSize];
+					// vp += 1.0;
 				}
-				_pressure[index] /= pcnt;
+				
+				if(pcnt > 0.000001f)
+				 	_pressure[index] /= pcnt;
+
+				// test - dg
+				// if(vp > 0.000001f)
+				//  	_zVelocity[index] /= vp;
 
 				// TODO? set correct velocity bc's
 				// velocities are only set to zero right now
@@ -533,6 +539,44 @@ void FLUID_3D::setObstacleBoundaries()
 		}
 }
 
+void FLUID_3D::setObstacleBoundaries()
+{
+	// cull degenerate obstacles , move to addObstacle?
+	for (int z = 1, index = _slabSize + _xRes + 1;
+			z < _zRes - 1; z++, index += 2 * _xRes)
+		for (int y = 1; y < _yRes - 1; y++, index += 2)
+			for (int x = 1; x < _xRes - 1; x++, index++)
+			{
+				if (_obstacles[index] != EMPTY)
+				{
+					const int top   = _obstacles[index + _slabSize];
+					const int bottom= _obstacles[index - _slabSize];
+					const int up    = _obstacles[index + _xRes];
+					const int down  = _obstacles[index - _xRes];
+					const int left  = _obstacles[index - 1];
+					const int right = _obstacles[index + 1];
+
+					int counter = 0;
+					if (up)    counter++;
+					if (down)  counter++;
+					if (left)  counter++;
+					if (right) counter++;
+					if (top)  counter++;
+					if (bottom) counter++;
+
+					if (counter < 3)
+						_obstacles[index] = EMPTY;
+				}
+				if (_obstacles[index])
+				{
+					_xVelocity[index] =
+					_yVelocity[index] =
+					_zVelocity[index] = 0.0f;
+					_pressure[index] = 0.0f;
+				}
+			}
+}
+
 //////////////////////////////////////////////////////////////////////
 // add buoyancy forces
 //////////////////////////////////////////////////////////////////////
diff --git a/intern/smoke/intern/FLUID_3D.h b/intern/smoke/intern/FLUID_3D.h
index 2d212caa6d3..7548f48250d 100644
--- a/intern/smoke/intern/FLUID_3D.h
+++ b/intern/smoke/intern/FLUID_3D.h
@@ -132,6 +132,7 @@ class FLUID_3D
 
 		// handle obstacle boundaries
 		void setObstacleBoundaries();
+		void setObstaclePressure();
 
 	public:
 		// advection, accessed e.g. by WTURBULENCE class
diff --git a/intern/smoke/intern/FLUID_3D_STATIC.cpp b/intern/smoke/intern/FLUID_3D_STATIC.cpp
index f2bbcf33075..7d3fe0c1c8d 100644
--- a/intern/smoke/intern/FLUID_3D_STATIC.cpp
+++ b/intern/smoke/intern/FLUID_3D_STATIC.cpp
@@ -80,7 +80,7 @@ void FLUID_3D::addSmokeTestCase(float* field, Vec3Int res, float value)
 void FLUID_3D::setNeumannX(float* field, Vec3Int res)
 {
 	const int slabSize = res[0] * res[1];
-	int index;
+	size_t index;
 	for (int z = 0; z < res[2]; z++)
 		for (int y = 0; y < res[1]; y++)
 		{
@@ -100,7 +100,7 @@ void FLUID_3D::setNeumannX(float* field, Vec3Int res)
 void FLUID_3D::setNeumannY(float* field, Vec3Int res)
 {
 	const int slabSize = res[0] * res[1];
-	int index;
+	size_t index;
 	for (int z = 0; z < res[2]; z++)
 		for (int x = 0; x < res[0]; x++)
 		{
@@ -121,7 +121,7 @@ void FLUID_3D::setNeumannZ(float* field, Vec3Int res)
 {
 	const int slabSize = res[0] * res[1];
 	const int totalCells = res[0] * res[1] * res[2];
-	int index;
+	size_t index;
 	for (int y = 0; y < res[1]; y++)
 		for (int x = 0; x < res[0]; x++)
 		{
@@ -141,10 +141,11 @@ void FLUID_3D::setNeumannZ(float* field, Vec3Int res)
 			// top slab
 			index = x + y * res[0];
 			index += totalCells - slabSize;
-			if(field[index]<0.) field[index] = 0.;
+			if(field[index]<0.) field[index] = 0.0f;
 			index -= slabSize;
-			if(field[index]<0.) field[index] = 0.;
+			if(field[index]<0.) field[index] = 0.0f;
 		}
+		
 }
 
 //////////////////////////////////////////////////////////////////////
diff --git a/intern/smoke/intern/smoke_API.cpp b/intern/smoke/intern/smoke_API.cpp
index 9c835cf87ee..281761bd387 100644
--- a/intern/smoke/intern/smoke_API.cpp
+++ b/intern/smoke/intern/smoke_API.cpp
@@ -88,15 +88,18 @@ extern "C" float *smoke_get_velocity_z(FLUID_3D *fluid)
 
 extern "C" float *smoke_get_bigdensity(FLUID_3D *fluid)
 {
-	return fluid->_wTurbulence->getDensityBig();
+	return fluid->_wTurbulence ? fluid->_wTurbulence->getDensityBig() : NULL;
 }
 
 extern "C" void smoke_get_bigres(FLUID_3D *fluid, int *res)
 {
-	Vec3Int r = fluid->_wTurbulence->getResBig();
-	res[0] = r[0];
-	res[1] = r[1];
-	res[2] = r[2];
+	if(fluid->_wTurbulence)
+	{
+		Vec3Int r = fluid->_wTurbulence->getResBig();
+		res[0] = r[0];
+		res[1] = r[1];
+		res[2] = r[2];
+	}
 }
 
 extern "C" unsigned char *smoke_get_obstacle(FLUID_3D *fluid)