Sound Branch:

* Last commit before merge to 2.5
* Merge until revision 22324 of 2.5 branch
* Fixes for some compiler warnings
* Updated some build system files

15 files changed, 368 insertions, 106 deletions

diff --git a/intern/audaspace/OpenAL/AUD_OpenALDevice.cpp b/intern/audaspace/OpenAL/AUD_OpenALDevice.cpp
index e4c286c1aa4..f94b11a11b8 100644
--- a/intern/audaspace/OpenAL/AUD_OpenALDevice.cpp
+++ b/intern/audaspace/OpenAL/AUD_OpenALDevice.cpp
@@ -703,7 +703,6 @@ bool AUD_OpenALDevice::pause(AUD_Handle* handle)
 
 bool AUD_OpenALDevice::resume(AUD_Handle* handle)
 {
-	ALint info;
 	lock();
 
 	// only songs that are paused can be resumed
@@ -1196,8 +1195,9 @@ bool AUD_OpenALDevice::setSetting(AUD_3DSetting setting, float value)
 	case AUD_3DS_SPEED_OF_SOUND:
 		alSpeedOfSound(value);
 		return true;
+	default:
+		return false;
 	}
-	return false;
 }
 
 float AUD_OpenALDevice::getSetting(AUD_3DSetting setting)
@@ -1226,8 +1226,9 @@ float AUD_OpenALDevice::getSetting(AUD_3DSetting setting)
 		return alGetFloat(AL_DOPPLER_FACTOR);
 	case AUD_3DS_SPEED_OF_SOUND:
 		return alGetFloat(AL_SPEED_OF_SOUND);
+	default:
+		return std::numeric_limits<float>::quiet_NaN();
 	}
-	return std::numeric_limits<float>::quiet_NaN();
 }
 
 bool AUD_OpenALDevice::updateSource(AUD_Handle* handle, AUD_3DData &data)
@@ -1298,6 +1299,8 @@ bool AUD_OpenALDevice::setSourceSetting(AUD_Handle* handle,
 			alSourcef(source, AL_ROLLOFF_FACTOR, value);
 			result = true;
 			break;
+		default:
+			break;
 		}
 	}
 
@@ -1349,6 +1352,8 @@ float AUD_OpenALDevice::getSourceSetting(AUD_Handle* handle,
 		case AUD_3DSS_ROLLOFF_FACTOR:
 			alGetSourcef(source, AL_ROLLOFF_FACTOR, &result);
 			break;
+		default:
+			break;
 		}
 	}
 
diff --git a/intern/audaspace/OpenAL/AUD_OpenALDevice.h b/intern/audaspace/OpenAL/AUD_OpenALDevice.h
index 8864ccc4e01..074cd3d1924 100644
--- a/intern/audaspace/OpenAL/AUD_OpenALDevice.h
+++ b/intern/audaspace/OpenAL/AUD_OpenALDevice.h
@@ -32,8 +32,8 @@ struct AUD_OpenALHandle;
 struct AUD_OpenALBufferedFactory;
 class AUD_ConverterFactory;
 
-#include <al.h>
-#include <alc.h>
+#include <AL/al.h>
+#include <AL/alc.h>
 #include <list>
 #include <pthread.h>
 
diff --git a/intern/audaspace/intern/AUD_C-API.cpp b/intern/audaspace/intern/AUD_C-API.cpp
index 5039b16de2c..830e1147cbb 100644
--- a/intern/audaspace/intern/AUD_C-API.cpp
+++ b/intern/audaspace/intern/AUD_C-API.cpp
@@ -479,11 +479,9 @@ int AUD_setSoundVolume(AUD_Handle* handle, float volume)
 		{
 			return AUD_device->setCapability(AUD_CAPS_SOURCE_VOLUME, &caps);
 		}
-		catch(AUD_Exception e)
-		{
-			return false;
-		}
+		catch(AUD_Exception e) {}
 	}
+	return false;
 }
 
 int AUD_setSoundPitch(AUD_Handle* handle, float pitch)
@@ -499,11 +497,9 @@ int AUD_setSoundPitch(AUD_Handle* handle, float pitch)
 		{
 			return AUD_device->setCapability(AUD_CAPS_SOURCE_PITCH, &caps);
 		}
-		catch(AUD_Exception e)
-		{
-			return false;
-		}
+		catch(AUD_Exception e) {}
 	}
+	return false;
 }
 
 AUD_Device* AUD_openReadDevice(AUD_Specs specs)
diff --git a/intern/audaspace/intern/AUD_SoftwareDevice.cpp b/intern/audaspace/intern/AUD_SoftwareDevice.cpp
index 1a2a53fed19..174ff8c8979 100644
--- a/intern/audaspace/intern/AUD_SoftwareDevice.cpp
+++ b/intern/audaspace/intern/AUD_SoftwareDevice.cpp
@@ -95,7 +95,7 @@ void AUD_SoftwareDevice::mix(sample_t* buffer, int length)
 	lock();
 
 	AUD_SoftwareHandle* sound;
-	int len, left;
+	int len;
 	sample_t* buf;
 	int sample_size = AUD_SAMPLE_SIZE(m_specs);
 	std::list<AUD_SoftwareHandle*> stopSounds;
diff --git a/intern/smoke/CMakeLists.txt b/intern/smoke/CMakeLists.txt
index 90768e58be0..0db6acb683f 100644
--- a/intern/smoke/CMakeLists.txt
+++ b/intern/smoke/CMakeLists.txt
@@ -25,7 +25,6 @@
 # ***** END GPL LICENSE BLOCK *****
 
 SET(INC ${PNG_INC} ${ZLIB_INC} intern ../../extern/bullet2/src ../memutil ../guardealloc)
-# ${FFTW3_INC}
 
 FILE(GLOB SRC intern/*.cpp)
 
@@ -33,6 +32,11 @@ IF(WITH_OPENMP)
     ADD_DEFINITIONS(-DPARALLEL=1)
 ENDIF(WITH_OPENMP)
 
+IF(WITH_FFTW3)
+  ADD_DEFINITIONS(-DFFTW3=1)
+  SET(INC ${INC} ${FFTW3_INC})
+ENDIF(WITH_FFTW3)
+
 
 BLENDERLIB(bf_smoke "${SRC}" "${INC}")
 #, libtype='blender', priority = 0 )
diff --git a/intern/smoke/SConscript b/intern/smoke/SConscript
index c427764e19e..af5bf1aeb20 100644
--- a/intern/smoke/SConscript
+++ b/intern/smoke/SConscript
@@ -10,6 +10,9 @@ if env['WITH_BF_OPENMP']:
 
 incs = env['BF_PNG_INC'] + ' ' + env['BF_ZLIB_INC'] 
 incs += ' intern ../../extern/bullet2/src ../memutil ../guardealloc '
-# incs += env['BF_FFTW3_INC'] 
+
+if env['WITH_BF_FFTW3']:
+    defs += ' FFTW3=1'
+    incs += env['BF_FFTW3_INC'] 
 
 env.BlenderLib ('bf_smoke', sources, Split(incs), Split(defs), libtype=['intern'], priority=[40] )
diff --git a/intern/smoke/extern/smoke_API.h b/intern/smoke/extern/smoke_API.h
index b7819d06edc..ba14a247bf1 100644
--- a/intern/smoke/extern/smoke_API.h
+++ b/intern/smoke/extern/smoke_API.h
@@ -36,11 +36,9 @@ struct FLUID_3D *smoke_init(int *res, int amplify, float *p0, float *p1, float d
 void smoke_free(struct FLUID_3D *fluid);
 
 void smoke_initBlenderRNA(struct FLUID_3D *fluid, float *alpha, float *beta);
-
 void smoke_step(struct FLUID_3D *fluid);
 
 float *smoke_get_density(struct FLUID_3D *fluid);
-float *smoke_get_bigdensity(struct FLUID_3D *fluid);
 float *smoke_get_heat(struct FLUID_3D *fluid);
 float *smoke_get_velocity_x(struct FLUID_3D *fluid);
 float *smoke_get_velocity_y(struct FLUID_3D *fluid);
@@ -51,9 +49,15 @@ unsigned char *smoke_get_obstacle(struct FLUID_3D *fluid);
 size_t smoke_get_index(int x, int max_x, int y, int max_y, int z);
 size_t smoke_get_index2d(int x, int max_x, int y);
 
-void smoke_set_noise(struct FLUID_3D *fluid, int type);
+// wavelet turbulence functions
+struct WTURBULENCE *smoke_turbulence_init(int *res, int amplify, int noisetype);
+void smoke_turbulence_free(struct WTURBULENCE *wt);
+void smoke_turbulence_step(struct WTURBULENCE *wt, struct FLUID_3D *fluid);
 
-void smoke_get_bigres(struct FLUID_3D *fluid, int *res);
+float *smoke_turbulence_get_density(struct WTURBULENCE *wt);
+void smoke_turbulence_get_res(struct WTURBULENCE *wt, int *res);
+void smoke_turbulence_set_noise(struct WTURBULENCE *wt, int type);
+void smoke_initWaveletBlenderRNA(struct WTURBULENCE *wt, float *strength);
 
 #ifdef __cplusplus
 }
diff --git a/intern/smoke/intern/FFT_NOISE.h b/intern/smoke/intern/FFT_NOISE.h
index 9ae9682e2a0..e1000f72fb3 100644
--- a/intern/smoke/intern/FFT_NOISE.h
+++ b/intern/smoke/intern/FFT_NOISE.h
@@ -22,7 +22,7 @@
 #ifndef FFT_NOISE_H_
 #define FFT_NOISE_H_
 
-#if 0
+#if FFTW3==1
 #include <iostream>
 #include <fftw3.h>
 #include <MERSENNETWISTER.h>
diff --git a/intern/smoke/intern/FLUID_3D.cpp b/intern/smoke/intern/FLUID_3D.cpp
index 344883866a8..7e3a00417c1 100644
--- a/intern/smoke/intern/FLUID_3D.cpp
+++ b/intern/smoke/intern/FLUID_3D.cpp
@@ -59,7 +59,12 @@ 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, noisetype);
+	else
+		_wTurbulence = NULL;
+	*/
 	
 	// scale the constants according to the refinement of the grid
 	_dx = 1.0f / (float)_maxRes;
@@ -93,6 +98,8 @@ FLUID_3D::FLUID_3D(int *res, int amplify, float *p0, float dt) :
 	_xVorticity   = new float[_totalCells];
 	_yVorticity   = new float[_totalCells];
 	_zVorticity   = new float[_totalCells];
+	_h			  = new float[_totalCells];
+	_Precond	  = new float[_totalCells];
 
 	for (int x = 0; x < _totalCells; x++)
 	{
@@ -115,6 +122,8 @@ FLUID_3D::FLUID_3D(int *res, int amplify, float *p0, float dt) :
 		_yVorticity[x]   = 0.0f;
 		_zVorticity[x]   = 0.0f;
 		_residual[x]     = 0.0f;
+		_h[x]			 = 0.0f;
+		_Precond[x]		 = 0.0f;
 		_obstacles[x]    = false;
 	}
 
@@ -186,8 +195,10 @@ FLUID_3D::~FLUID_3D()
 	if (_yVorticity) delete[] _yVorticity;
 	if (_zVorticity) delete[] _zVorticity;
 	if (_vorticity) delete[] _vorticity;
+	if (_h) delete[] _h;
+	if (_Precond) delete[] _Precond;
 	if (_obstacles) delete[] _obstacles;
-    if (_wTurbulence) delete _wTurbulence;
+    // if (_wTurbulence) delete _wTurbulence;
 
     printf("deleted fluid\n");
 }
@@ -197,10 +208,6 @@ void FLUID_3D::initBlenderRNA(float *alpha, float *beta)
 {
 	_alpha = alpha;
 	_beta = beta;
-	
-	// XXX TODO DEBUG
-	// *_alpha = 0;
-	// *_beta = 0;
 }
 
 //////////////////////////////////////////////////////////////////////
@@ -215,21 +222,21 @@ void FLUID_3D::step()
 	wipeBoundaries();
 
 	// run the solvers
-  addVorticity();
-  addBuoyancy(_heat, _density);
+	addVorticity();
+	addBuoyancy(_heat, _density);
 	addForce();
 	project();
-  diffuseHeat();
+	diffuseHeat();
 
 	// advect everything
 	advectMacCormack();
 
-	if(_wTurbulence) {
-		_wTurbulence->stepTurbulenceFull(_dt/_dx,
-				_xVelocity, _yVelocity, _zVelocity, _obstacles);
+	// if(_wTurbulence) {
+	// 	_wTurbulence->stepTurbulenceFull(_dt/_dx,
+	//			_xVelocity, _yVelocity, _zVelocity, _obstacles);
 		// _wTurbulence->stepTurbulenceReadable(_dt/_dx,
 		//  _xVelocity, _yVelocity, _zVelocity, _obstacles);
-	}
+	// }
 /*
  // no file output
   float *src = _density;
@@ -248,7 +255,7 @@ void FLUID_3D::step()
   IMAGE::dumpPBRT(_totalSteps, pbrtPrefix, _density, _res[0],_res[1],_res[2]);
   */
 	_totalTime += _dt;
-	_totalSteps++;
+	_totalSteps++;	
 }
 
 //////////////////////////////////////////////////////////////////////
@@ -354,6 +361,7 @@ void FLUID_3D::addForce()
 void FLUID_3D::project()
 {
 	int index, x, y, z;
+
 	setObstacleBoundaries();
 
 	// copy out the boundaries
@@ -390,12 +398,16 @@ void FLUID_3D::project()
 						xright - xleft +
 						yup - ydown +
 						ztop - zbottom );
-				_pressure[index] = 0.0f;
+
+				// DG: commenting this helps CG to get a better start, 10-20% speed improvement
+				// _pressure[index] = 0.0f;
 			}
 	copyBorderAll(_pressure);
 
 	// solve Poisson equation
-	solvePressure(_pressure, _divergence, _obstacles);
+	solvePressurePre(_pressure, _divergence, _obstacles);
+
+	setObstaclePressure();
 
 	// project out solution
 	float invDx = 1.0f / _dx;
@@ -404,9 +416,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 +458,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 +467,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 +487,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 +505,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 +532,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..ffda34021c3 100644
--- a/intern/smoke/intern/FLUID_3D.h
+++ b/intern/smoke/intern/FLUID_3D.h
@@ -27,7 +27,7 @@
 #include <cmath>
 #include <iostream>
 #include "OBSTACLE.h"
-#include "WTURBULENCE.h"
+// #include "WTURBULENCE.h"
 #include "VEC3.h"
 
 using namespace std;
@@ -96,6 +96,8 @@ class FLUID_3D
 		float* _yVorticity;
 		float* _zVorticity;
 		float* _vorticity;
+		float* _h;
+		float* _Precond;
 		unsigned char*  _obstacles;
 
 		// CG fields
@@ -113,7 +115,7 @@ class FLUID_3D
 		float _tempAmb; /* ambient temperature */
 
 		// WTURBULENCE object, if active
-		WTURBULENCE* _wTurbulence;
+		// WTURBULENCE* _wTurbulence;
 
 		// boundary setting functions
 		void copyBorderAll(float* field);
@@ -128,10 +130,12 @@ class FLUID_3D
 		void project();
 		void diffuseHeat();
 		void solvePressure(float* field, float* b, unsigned char* skip);
+		void solvePressurePre(float* field, float* b, unsigned char* skip);
 		void solveHeat(float* field, float* b, unsigned char* skip);
 
 		// handle obstacle boundaries
 		void setObstacleBoundaries();
+		void setObstaclePressure();
 
 	public:
 		// advection, accessed e.g. by WTURBULENCE class
diff --git a/intern/smoke/intern/FLUID_3D_SOLVERS.cpp b/intern/smoke/intern/FLUID_3D_SOLVERS.cpp
index 5fd8f72d79d..51929e1194b 100644
--- a/intern/smoke/intern/FLUID_3D_SOLVERS.cpp
+++ b/intern/smoke/intern/FLUID_3D_SOLVERS.cpp
@@ -21,8 +21,171 @@
 //////////////////////////////////////////////////////////////////////
 
 #include "FLUID_3D.h"
+
 #define SOLVER_ACCURACY 1e-06
 
+void FLUID_3D::solvePressurePre(float* field, float* b, unsigned char* skip)
+{
+	int x, y, z, index;
+
+	// i = 0
+	int i = 0;
+
+	// r = b - Ax
+	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.;
+			}
+		    
+			_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) );
+			_residual[index] = (skip[index]) ? 0.0f : _residual[index];
+
+			// P^-1
+			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 = transpose(r) * p
+	float deltaNew = 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++)
+			deltaNew += _residual[index] * _direction[index];
+
+	// delta0 = deltaNew
+	float delta0 = deltaNew;
+
+  // 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))
+  {
+    // (s) q = Ad (p)
+    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);
+          _q[index] = (skip[index]) ? 0.0f : _q[index];
+        }
+
+    // alpha = deltaNew / (transpose(d) * q)
+    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++)
+          alpha += _direction[index] * _q[index];
+    if (fabs(alpha) > 0.0f)
+      alpha = deltaNew / alpha;
+
+    // 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];
+
+    // r = r - alpha * q
+	maxR = 0.0;
+    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++)
+        {
+          _residual[index] -= alpha * _q[index];
+		  // maxR = (_residual[index] > maxR) ? _residual[index] : maxR;
+        }
+
+	// if(maxR <= eps)
+	// 	break;
+
+	// h = P^-1 * r
+	 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++)
+		{
+			_h[index] = _Precond[index] * _residual[index];
+		}
+
+    // deltaOld = deltaNew
+    float deltaOld = deltaNew;
+
+    // deltaNew = transpose(r) * h
+    deltaNew = 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++)
+		{
+          deltaNew += _residual[index] * _h[index];
+		  maxR = (_residual[index]* _h[index] > maxR) ? _residual[index]* _h[index] : maxR;
+		}
+
+    // beta = deltaNew / deltaOld
+    float beta = deltaNew / deltaOld;
+
+    // d = h + beta * 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++)
+          _direction[index] = _h[index] + beta * _direction[index];
+
+    // i = i + 1
+    i++;
+  }
+  // cout << i << " iterations converged to " << sqrt(maxR) << endl;
+}
+
 //////////////////////////////////////////////////////////////////////
 // solve the poisson equation with CG
 //////////////////////////////////////////////////////////////////////
@@ -61,6 +224,7 @@ void FLUID_3D::solvePressure(float* field, float* b, unsigned char* skip)
           field[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -1.0f) );
         _residual[index] = (skip[index]) ? 0.0f : _residual[index];
       }
+	
 
   // d = r
   index = _slabSize + _xRes + 1;
@@ -166,7 +330,7 @@ void FLUID_3D::solvePressure(float* field, float* b, unsigned char* skip)
     // i = i + 1
     i++;
   }
-  cout << i << " iterations converged to " << maxR << endl;
+  // cout << i << " iterations converged to " << maxR << endl;
 }
 
 //////////////////////////////////////////////////////////////////////
@@ -314,6 +478,6 @@ void FLUID_3D::solveHeat(float* field, float* b, unsigned char* skip)
     // i = i + 1
     i++;
   }
-  cout << i << " iterations converged to " << maxR << endl;
+  // cout << i << " iterations converged to " << maxR << endl;
 }
 
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/WTURBULENCE.cpp b/intern/smoke/intern/WTURBULENCE.cpp
index 022c8f28252..06d8d84e821 100644
--- a/intern/smoke/intern/WTURBULENCE.cpp
+++ b/intern/smoke/intern/WTURBULENCE.cpp
@@ -43,7 +43,7 @@ static const float persistence = 0.56123f;
 //////////////////////////////////////////////////////////////////////
 // constructor
 //////////////////////////////////////////////////////////////////////
-WTURBULENCE::WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify)
+WTURBULENCE::WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify, int noisetype)
 {
 	// if noise magnitude is below this threshold, its contribution
 	// is negilgible, so stop evaluating new octaves
@@ -53,10 +53,10 @@ WTURBULENCE::WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify)
 	_amplify = amplify;
 	
 	// manually adjust the overall amount of turbulence
-	_strength = 2.;
+	// DG - RNA-fied _strength = 2.;
 	
 	// add the corresponding octaves of noise
-	_octaves = (int)log((float)_amplify) / log(2.0f); // XXX DEBUG/ TODO: int casting correct? - dg
+	_octaves = (int)(log((float)_amplify) / log(2.0f) + 0.5); // XXX DEBUG/ TODO: int casting correct? - dg
 	
 	// noise resolution
 	_xResBig = _amplify * xResSm;
@@ -136,8 +136,11 @@ WTURBULENCE::WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify)
 	
 	// noise tiles
 	_noiseTile = new float[noiseTileSize * noiseTileSize * noiseTileSize];
+	/*
 	std::string noiseTileFilename = std::string("noise.wavelets");
 	generateTile_WAVELET(_noiseTile, noiseTileFilename);
+	*/
+	setNoise(noisetype);
 	/*
 	std::string noiseTileFilename = std::string("noise.fft");
 	generatTile_FFT(_noiseTile, noiseTileFilename);
@@ -173,19 +176,21 @@ WTURBULENCE::~WTURBULENCE() {
 //////////////////////////////////////////////////////////////////////
 // Change noise type
 //
-// type (1<<1) = wavelet / 2
-// type (1<<2) = FFT / 4
-// type (1<<3) = curl / 8
+// type (1<<0) = wavelet / 2
+// type (1<<1) = FFT / 4
+// type (1<<2) = curl / 8
 //////////////////////////////////////////////////////////////////////
 void WTURBULENCE::setNoise(int type)
 {
-	if(type == 4) // FFT
+	if(type == (1<<1)) // FFT
 	{
 		// needs fft
-		// std::string noiseTileFilename = std::string("noise.fft");
-		// generatTile_FFT(_noiseTile, noiseTileFilename);
+		#if FFTW3==1
+		std::string noiseTileFilename = std::string("noise.fft");
+		generatTile_FFT(_noiseTile, noiseTileFilename);
+		#endif
 	}
-	else if(type == 8) // curl
+	else if(type == (1<<2)) // curl
 	{
 		// TODO: not supported yet
 	}
@@ -196,6 +201,12 @@ void WTURBULENCE::setNoise(int type)
 	}
 }
 
+// init direct access functions from blender
+void WTURBULENCE::initBlenderRNA(float *strength)
+{
+	_strength = strength;
+}
+
 //////////////////////////////////////////////////////////////////////
 // Get the smallest valid x derivative
 //
@@ -642,7 +653,7 @@ void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel,
 
         // base amplitude for octave 0
         float coefficient = sqrtf(2.0f * fabs(energy));
-        const float amplitude = _strength * fabs(0.5 * coefficient) * persistence;
+        const float amplitude = *_strength * fabs(0.5 * coefficient) * persistence;
 
         // add noise to velocity, but only if the turbulence is
         // sufficiently undeformed, and the energy is large enough
@@ -854,7 +865,7 @@ void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, floa
 
       // base amplitude for octave 0
       float coefficient = sqrtf(2.0f * fabs(energy));
-      const float amplitude = _strength * fabs(0.5 * coefficient) * persistence;
+      const float amplitude = *_strength * fabs(0.5 * coefficient) * persistence;
 
       // add noise to velocity, but only if the turbulence is
       // sufficiently undeformed, and the energy is large enough
@@ -925,7 +936,7 @@ void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, floa
   maxVelMag = sqrt(maxVelMag) * dt;
   int totalSubsteps = (int)(maxVelMag / (float)maxVel);
   totalSubsteps = (totalSubsteps < 1) ? 1 : totalSubsteps;
-  
+  // printf("totalSubsteps: %d\n", totalSubsteps);
   totalSubsteps = (totalSubsteps > maxSubSteps) ? maxSubSteps : totalSubsteps;
   const float dtSubdiv = dt / (float)totalSubsteps;
 
diff --git a/intern/smoke/intern/WTURBULENCE.h b/intern/smoke/intern/WTURBULENCE.h
index 858a47b7dd1..d4e6b0c6a17 100644
--- a/intern/smoke/intern/WTURBULENCE.h
+++ b/intern/smoke/intern/WTURBULENCE.h
@@ -33,12 +33,13 @@ class WTURBULENCE
 {
 	public:
 		// both config files can be NULL, altCfg might override values from noiseCfg
-		WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify);
+		WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify, int noisetype);
 
 		/// destructor
 		virtual ~WTURBULENCE();
 		
 		void setNoise(int type);
+		void initBlenderRNA(float *strength);
 
 		// step more readable version -- no rotation correction
 		void stepTurbulenceReadable(float dt, float* xvel, float* yvel, float* zvel, unsigned char *obstacles);
@@ -69,13 +70,15 @@ class WTURBULENCE
 		inline Vec3Int getResBig() { return _resBig; }
 		inline int getOctaves() { return _octaves; }
 
+		// is accessed on through rna gui
+		float *_strength;
+
 	protected:
 		// enlargement factor from original velocity field / simulation
 		// _Big = _amplify * _Sm
 		int _amplify;
 		int _octaves;
-		float _strength;
-
+		
 		// noise settings
 		float _cullingThreshold;
 		float _noiseStrength;
diff --git a/intern/smoke/intern/smoke_API.cpp b/intern/smoke/intern/smoke_API.cpp
index 9c835cf87ee..7d476191c90 100644
--- a/intern/smoke/intern/smoke_API.cpp
+++ b/intern/smoke/intern/smoke_API.cpp
@@ -26,6 +26,7 @@
  */
 
 #include "FLUID_3D.h"
+#include "WTURBULENCE.h"
 
 #include <stdio.h>
 #include <stdlib.h>
@@ -41,22 +42,48 @@ extern "C" FLUID_3D *smoke_init(int *res, int amplify, float *p0, float *p1, flo
 	return fluid;
 }
 
+extern "C" WTURBULENCE *smoke_turbulence_init(int *res, int amplify, int noisetype)
+{
+	// initialize wavelet turbulence
+	if(amplify)
+		return new WTURBULENCE(res[0],res[1],res[2], amplify, noisetype);
+	else 
+		return NULL;
+}
+
 extern "C" void smoke_free(FLUID_3D *fluid)
 {
 	delete fluid;
 	fluid = NULL;
 }
 
+extern "C" void smoke_turbulence_free(WTURBULENCE *wt)
+{
+	 delete wt;
+	 wt = NULL;
+}
+
 extern "C" void smoke_step(FLUID_3D *fluid)
 {
 	fluid->step();
 }
 
+extern "C" void smoke_turbulence_step(WTURBULENCE *wt, FLUID_3D *fluid)
+{
+	if(wt) 
+		wt->stepTurbulenceFull(fluid->_dt/fluid->_dx, fluid->_xVelocity, fluid->_yVelocity, fluid->_zVelocity, fluid->_obstacles); 
+}
+
 extern "C" void smoke_initBlenderRNA(FLUID_3D *fluid, float *alpha, float *beta)
 {
 	fluid->initBlenderRNA(alpha, beta);
 }
 
+extern "C" void smoke_initWaveletBlenderRNA(WTURBULENCE *wt, float *strength)
+{
+	wt->initBlenderRNA(strength);
+}
+
 template < class T > inline T ABS( T a ) {
 	return (0 < a) ? a : -a ;
 }
@@ -86,17 +113,20 @@ extern "C" float *smoke_get_velocity_z(FLUID_3D *fluid)
 	return fluid->_zVorticity;
 }
 
-extern "C" float *smoke_get_bigdensity(FLUID_3D *fluid)
+extern "C" float *smoke_turbulence_get_density(WTURBULENCE *wt)
 {
-	return fluid->_wTurbulence->getDensityBig();
+	return wt ? wt->getDensityBig() : NULL;
 }
 
-extern "C" void smoke_get_bigres(FLUID_3D *fluid, int *res)
+extern "C" void smoke_turbulence_get_res(WTURBULENCE *wt, int *res)
 {
-	Vec3Int r = fluid->_wTurbulence->getResBig();
-	res[0] = r[0];
-	res[1] = r[1];
-	res[2] = r[2];
+	if(wt)
+	{
+		Vec3Int r = wt->getResBig();
+		res[0] = r[0];
+		res[1] = r[1];
+		res[2] = r[2];
+	}
 }
 
 extern "C" unsigned char *smoke_get_obstacle(FLUID_3D *fluid)
@@ -115,7 +145,7 @@ extern "C" size_t smoke_get_index2d(int x, int max_x, int y /*, int max_y, int z
 	return x + y * max_x;
 }
 
-extern "C" void smoke_set_noise(FLUID_3D *fluid, int type)
+extern "C" void smoke_turbulence_set_noise(WTURBULENCE *wt, int type)
 {
-	fluid->_wTurbulence->setNoise(type);
+	wt->setNoise(type);
 }