Revision: 22287
http://projects.blender.org/plugins/scmsvn/viewcvs.php?view=rev&root=bf-blender&revision=22287
Author: genscher
Date: 2009-08-07 03:07:45 +0200 (Fri, 07 Aug 2009)
Log Message:
-----------
Smoke: test commit of PCG
Modified Paths:
--------------
branches/blender2.5/blender/intern/smoke/intern/FLUID_3D.cpp
branches/blender2.5/blender/intern/smoke/intern/FLUID_3D.h
branches/blender2.5/blender/intern/smoke/intern/FLUID_3D_SOLVERS.cpp
Modified: branches/blender2.5/blender/intern/smoke/intern/FLUID_3D.cpp
===================================================================
--- branches/blender2.5/blender/intern/smoke/intern/FLUID_3D.cpp
2009-08-07 01:05:33 UTC (rev 22286)
+++ branches/blender2.5/blender/intern/smoke/intern/FLUID_3D.cpp
2009-08-07 01:07:45 UTC (rev 22287)
@@ -96,6 +96,8 @@
_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++)
{
@@ -118,6 +120,8 @@
_yVorticity[x] = 0.0f;
_zVorticity[x] = 0.0f;
_residual[x] = 0.0f;
+ _h[x] = 0.0f;
+ _Precond[x] = 0.0f;
_obstacles[x] = false;
}
@@ -189,6 +193,8 @@
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;
@@ -414,7 +420,7 @@
copyBorderAll(_pressure);
// solve Poisson equation
- solvePressure(_pressure, _divergence, _obstacles);
+ solvePressurePre(_pressure, _divergence, _obstacles);
setObstaclePressure();
Modified: branches/blender2.5/blender/intern/smoke/intern/FLUID_3D.h
===================================================================
--- branches/blender2.5/blender/intern/smoke/intern/FLUID_3D.h 2009-08-07
01:05:33 UTC (rev 22286)
+++ branches/blender2.5/blender/intern/smoke/intern/FLUID_3D.h 2009-08-07
01:07:45 UTC (rev 22287)
@@ -96,6 +96,8 @@
float* _yVorticity;
float* _zVorticity;
float* _vorticity;
+ float* _h;
+ float* _Precond;
unsigned char* _obstacles;
// CG fields
@@ -128,6 +130,7 @@
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
Modified: branches/blender2.5/blender/intern/smoke/intern/FLUID_3D_SOLVERS.cpp
===================================================================
--- branches/blender2.5/blender/intern/smoke/intern/FLUID_3D_SOLVERS.cpp
2009-08-07 01:05:33 UTC (rev 22286)
+++ branches/blender2.5/blender/intern/smoke/intern/FLUID_3D_SOLVERS.cpp
2009-08-07 01:07:45 UTC (rev 22287)
@@ -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 > 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 " << maxR << endl;
+}
+
//////////////////////////////////////////////////////////////////////
// solve the poisson equation with CG
//////////////////////////////////////////////////////////////////////
@@ -61,6 +224,7 @@
field[index + _slabSize] * (skip[index + _slabSize] ? 0.0 : -1.0f) );
_residual[index] = (skip[index]) ? 0.0f : _residual[index];
}
+
// d = r
index = _slabSize + _xRes + 1;
@@ -314,6 +478,6 @@
// i = i + 1
i++;
}
- cout << i << " iterations converged to " << maxR << endl;
+ // cout << i << " iterations converged to " << maxR << endl;
}
_______________________________________________
Bf-blender-cvs mailing list
[email protected]
http://lists.blender.org/mailman/listinfo/bf-blender-cvs
