[Bf-blender-cvs] SVN commit: /data/svn/bf-blender [48055] branches/smoke2/intern/smoke/ intern/FLUID_3D_SOLVERS.cpp: - Add old CG for reference

Daniel Genrich daniel.genrich at gmx.net
Mon Jun 18 21:49:47 CEST 2012


Revision: 48055
          http://projects.blender.org/scm/viewvc.php?view=rev&root=bf-blender&revision=48055
Author:   genscher
Date:     2012-06-18 19:49:47 +0000 (Mon, 18 Jun 2012)
Log Message:
-----------
- Add old CG for reference

Modified Paths:
--------------
    branches/smoke2/intern/smoke/intern/FLUID_3D_SOLVERS.cpp

Modified: branches/smoke2/intern/smoke/intern/FLUID_3D_SOLVERS.cpp
===================================================================
--- branches/smoke2/intern/smoke/intern/FLUID_3D_SOLVERS.cpp	2012-06-18 19:37:01 UTC (rev 48054)
+++ branches/smoke2/intern/smoke/intern/FLUID_3D_SOLVERS.cpp	2012-06-18 19:49:47 UTC (rev 48055)
@@ -128,37 +128,37 @@
 							alpha += _direction[index] * _q[index];
 						}
 
-						if (fabs(alpha) > 0.0f)
-							alpha = deltaNew / alpha;
+				if (fabs(alpha) > 0.0f)
+					alpha = deltaNew / alpha;
 
-						float deltaOld = deltaNew;
-						deltaNew = 0.0f;
+				float deltaOld = deltaNew;
+				deltaNew = 0.0f;
 
-						maxR = 0.0f;
+				maxR = 0.0f;
 
-						index = _slabSize + _xRes + 1;
-						for (z = 1; z < _zRes - 1; z++, index += twoxr)
-							for (y = 1; y < _yRes - 1; y++, index += 2)
-								for (x = 1; x < _xRes - 1; x++, index++)
-								{
-									field[index] += alpha * _direction[index];
+				index = _slabSize + _xRes + 1;
+				for (z = 1; z < _zRes - 1; z++, index += twoxr)
+					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];
-									maxR = (_residual[index] > maxR) ? _residual[index] : maxR;
+							_residual[index] -= alpha * _q[index];
+							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;
-								for (z = 1; z < _zRes - 1; z++, index += twoxr)
-									for (y = 1; y < _yRes - 1; y++, index += 2)
-										for (x = 1; x < _xRes - 1; x++, index++)
-											_direction[index] = _residual[index] + beta * _direction[index];
+				index = _slabSize + _xRes + 1;
+				for (z = 1; z < _zRes - 1; z++, index += twoxr)
+					for (y = 1; y < _yRes - 1; y++, index += 2)
+						for (x = 1; x < _xRes - 1; x++, index++)
+							_direction[index] = _residual[index] + beta * _direction[index];
 
 
-								i++;
+				i++;
 			}
 			// cout << i << " iterations converged to " << maxR << endl;
 
@@ -169,6 +169,161 @@
 }
 
 
+#if 1
+void FLUID_3D::solvePressurePre(float* field, float* b, unsigned char* skip)
+{
+  int x, y, z, index;
+  float *_q, *_h, *_residual, *_direction;
+
+  // i = 0
+  int i = 0;
+
+ _residual  = new float[_totalCells]; // set 0
+ _h  = new float[_totalCells]; // set 0
+ _q  = new float[_totalCells]; // set 0
+ _direction  = new float[_totalCells]; // set 0
+
+ memset(_residual, 0, sizeof(float)*_xRes*_yRes*_zRes);
+ memset(_h, 0, sizeof(float)*_xRes*_yRes*_zRes);
+ memset(_q, 0, sizeof(float)*_xRes*_yRes*_zRes);
+ memset(_direction, 0, sizeof(float)*_xRes*_yRes*_zRes);
+
+  // 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];
+      }
+
+  // d = 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++)
+        _direction[index] = _residual[index];
+
+  // deltaNew = transpose(r) * r
+  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] * _residual[index];
+
+  // delta0 = deltaNew
+  float delta0 = deltaNew;
+
+  // While deltaNew > (eps^2) * delta0
+  const float eps  = SOLVER_ACCURACY;
+  float maxR = 2.0f * eps;
+  while ((i < _iterations) && (maxR > eps))
+  {
+    // q = Ad
+    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.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++)
+        {
+          _residual[index] -= alpha * _q[index];
+          maxR = (_residual[index] > maxR) ? _residual[index] : maxR;
+        }
+
+    // deltaOld = deltaNew
+    float deltaOld = deltaNew;
+
+    // deltaNew = transpose(r) * r
+    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] * _residual[index];
+
+    // beta = deltaNew / deltaOld
+    float beta = deltaNew / deltaOld;
+
+    // d = r + 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] = _residual[index] + beta * _direction[index];
+
+    // i = i + 1
+    i++;
+  }
+  cout << i << " iterations converged to " << maxR << endl;
+}
+#else
 #define SMOKE_LOOP \
 		index = _slabSize + _xRes + 1; \
 		for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes) \
@@ -344,3 +499,4 @@
 	if (_z) delete[] _z;
 	if (_q) delete[] _q;
 }
+#endif
\ No newline at end of file




More information about the Bf-blender-cvs mailing list