[Bf-blender-cvs] SVN commit: /data/svn/bf-blender [47269] branches/smoke2/intern/smoke/ intern: Eigen3 "Kinda" working - still weird splitting up going on in 3dview
Daniel Genrich
daniel.genrich at gmx.net
Thu May 31 12:54:28 CEST 2012
Revision: 47269
http://projects.blender.org/scm/viewvc.php?view=rev&root=bf-blender&revision=47269
Author: genscher
Date: 2012-05-31 10:54:27 +0000 (Thu, 31 May 2012)
Log Message:
-----------
Eigen3 "Kinda" working - still weird splitting up going on in 3dview
Modified Paths:
--------------
branches/smoke2/intern/smoke/intern/FLUID_3D.cpp
branches/smoke2/intern/smoke/intern/FLUID_3D.h
branches/smoke2/intern/smoke/intern/FLUID_3D_SOLVERS.cpp
branches/smoke2/intern/smoke/intern/FLUID_3D_STATIC.cpp
Modified: branches/smoke2/intern/smoke/intern/FLUID_3D.cpp
===================================================================
--- branches/smoke2/intern/smoke/intern/FLUID_3D.cpp 2012-05-31 10:38:11 UTC (rev 47268)
+++ branches/smoke2/intern/smoke/intern/FLUID_3D.cpp 2012-05-31 10:54:27 UTC (rev 47269)
@@ -45,8 +45,6 @@
#include <omp.h>
#endif // PARALLEL
-using namespace std;
-using namespace Eigen;
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
@@ -762,8 +760,6 @@
}
}
-#define INDEX(indexX, indexY, indexZ) ((indexZ) * _xRes * _yRes + (indexY) * _xRes + (indexX))
-
//////////////////////////////////////////////////////////////////////
// project into divergence free field
//////////////////////////////////////////////////////////////////////
@@ -777,8 +773,6 @@
SparseMatrix<float,RowMajor> A(_totalCells, _totalCells);
A.reserve(VectorXi::Constant(_totalCells, 7));
- ArrayXd gridToIndex(_totalCells);
-
ArrayXd A0(_totalCells);
ArrayXd Ai(_totalCells);
ArrayXd Aj(_totalCells);
@@ -787,14 +781,6 @@
b.setZero(_totalCells);
p.setZero(_totalCells);
- unsigned int linearIndex = 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++)
- if (!_obstacles[INDEX(x,y,z)])
- gridToIndex(INDEX(x,y,z)) = linearIndex++;
-
// float *_pressure = new float[_totalCells];
float *_divergence = new float[_totalCells];
@@ -856,7 +842,7 @@
// Pressure is zero anyway since now a local array is used
_pressure[index] = 0.0f;
}
-#if 0
+
float scale = 1.0; // DG TODO: make this global and incooperate this into other functions
index = _slabSize + _xRes + 1;
@@ -864,61 +850,71 @@
for (y = 1; y < _yRes - 1; y++, index += 2)
for (x = 1; x < _xRes - 1; x++, index++)
{
- if(!_obstacles[INDEX(x, y, z)])
- {
- if(!_obstacles[INDEX(x + 1, y, z)]) A0[INDEX(x, y, z)] += scale;
- if(!_obstacles[INDEX(x - 1, y, z)]) A0[INDEX(x, y, z)] += scale;
- if(!_obstacles[INDEX(x, y + 1, z)]) A0[INDEX(x, y, z)] += scale;
- if(!_obstacles[INDEX(x, y - 1, z)]) A0[INDEX(x, y, z)] += scale;
- if(!_obstacles[INDEX(x, y, z + 1)]) A0[INDEX(x, y, z)] += scale;
- if(!_obstacles[INDEX(x, y, z - 1)]) A0[INDEX(x, y, z)] += scale;
+ if(!_obstacles[FINDEX(x, y, z)])
+ {
+ if(!_obstacles[FINDEX(x + 1, y, z)]) A0[FINDEX(x, y, z)] += scale;
+ if(!_obstacles[FINDEX(x - 1, y, z)]) A0[FINDEX(x, y, z)] += scale;
+ if(!_obstacles[FINDEX(x, y + 1, z)]) A0[FINDEX(x, y, z)] += scale;
+ if(!_obstacles[FINDEX(x, y - 1, z)]) A0[FINDEX(x, y, z)] += scale;
+ if(!_obstacles[FINDEX(x, y, z + 1)]) A0[FINDEX(x, y, z)] += scale;
+ if(!_obstacles[FINDEX(x, y, z - 1)]) A0[FINDEX(x, y, z)] += scale;
+ }
- if(!_obstacles[INDEX(x + 1, y, z)]) Ai[INDEX(x, y, z)] = -scale;
- if(!_obstacles[INDEX(x, y + 1, z)]) Aj[INDEX(x, y, z)] = -scale;
- if(!_obstacles[INDEX(x, y, z + 1)]) Ak[INDEX(x, y, z)] = -scale;
- }
+ }
+ for (z = 0; z < _zRes - 1; z++)
+ for (y = 0; y < _yRes - 1; y++)
+ for (x = 0; x < _xRes - 1; x++)
+ {
+ if(!_obstacles[FINDEX(x, y, z)])
+ {
+ if(!_obstacles[FINDEX(x + 1, y, z)]) Ai[FINDEX(x, y, z)] = -scale;
+ if(!_obstacles[FINDEX(x, y + 1, z)]) Aj[FINDEX(x, y, z)] = -scale;
+ if(!_obstacles[FINDEX(x, y, z + 1)]) Ak[FINDEX(x, y, z)] = -scale;
+ }
+
}
- unsigned int valCount = 0;
unsigned int rowCount = 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++)
- if (!_obstacles[INDEX(x,y,z)])
+ for (z = 1; z < _zRes - 1; z++)
+ for (y = 1; y < _yRes - 1; y++)
+ for (x = 1; x < _xRes - 1; x++)
+ if (!_obstacles[FINDEX(x,y,z)])
{
- if(Ai[INDEX(x - 1, y, z)] < 0)
- A.insert(rowCount, gridToIndex(INDEX(x - 1, y, z))) = Ai[INDEX(x - 1, y, z)];
- if(Aj[INDEX(x, y - 1, z)] < 0)
- A.insert(rowCount, gridToIndex(INDEX(x, y - 1, z))) = Aj[INDEX(x, y - 1, z)];
- if(Ak[INDEX(x, y, z - 1)] < 0)
- A.insert(rowCount, gridToIndex(INDEX(x, y, z - 1))) = Ak[INDEX(x, y, z - 1)];
+ rowCount = FINDEX(x, y, z);
- A.insert(rowCount, gridToIndex(INDEX(x, y, z))) = A0[INDEX(x, y, z)];
+ if(Ai[FINDEX(x - 1, y, z)] < 0)
+ A.insert(rowCount, FINDEX(x - 1, y, z)) = Ai[FINDEX(x - 1, y, z)];
+ if(Aj[FINDEX(x, y - 1, z)] < 0)
+ A.insert(rowCount, FINDEX(x, y - 1, z)) = Aj[FINDEX(x, y - 1, z)];
+ if(Ak[FINDEX(x, y, z - 1)] < 0)
+ A.insert(rowCount, FINDEX(x, y, z - 1)) = Ak[FINDEX(x, y, z - 1)];
- if(Ai[INDEX(x, y, z)] < 0)
- A.insert(rowCount, gridToIndex(INDEX(x + 1, y, z))) = Ai[INDEX(x, y, z)];
- if(Aj[INDEX(x, y, z)] < 0)
- A.insert(rowCount, gridToIndex(INDEX(x, y + 1, z))) = Aj[INDEX(x, y, z)];
- if(Ak[INDEX(x, y, z)] < 0)
- A.insert(rowCount, gridToIndex(INDEX(x, y, z + 1))) = Ak[INDEX(x, y, z)];
+ if(A0[FINDEX(x, y, z)] > 0)
+ A.insert(rowCount, FINDEX(x, y, z)) = A0[FINDEX(x, y, z)];
- rowCount++;
+ if(Ai[FINDEX(x, y, z)] < 0)
+ A.insert(rowCount, FINDEX(x + 1, y, z)) = Ai[FINDEX(x, y, z)];
+ if(Aj[FINDEX(x, y, z)] < 0)
+ A.insert(rowCount, FINDEX(x, y + 1, z)) = Aj[FINDEX(x, y, z)];
+ if(Ak[FINDEX(x, y, z)] < 0)
+ A.insert(rowCount, FINDEX(x, y, z + 1)) = Ak[FINDEX(x, y, z)];
}
+ // DG TODO: for schleifen fixen
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 (!_obstacles[INDEX(x,y,z)])
- b[gridToIndex(INDEX(x, y, z))] = _divergence[INDEX(x,y,z)];
+ if (!_obstacles[FINDEX(x,y,z)])
+ b[FINDEX(x, y, z)] = _divergence[FINDEX(x,y,z)];
A.makeCompressed();
- ConjugateGradient<SparseMatrix<float>, Eigen::Symmetric > solver;
+
+ ConjugateGradient< SparseMatrix<float, RowMajor> > solver;
solver.setMaxIterations(100);
- solver.setTolerance(1e-06);
+ solver.setTolerance(1e-05);
solver.compute(A);
if(solver.info() != Success)
@@ -936,14 +932,24 @@
}
else
printf("Solver finished.\n");
+
+ std::cout << "#iterations: " << solver.iterations() << std::endl;
+ std::cout << "estimated error: " << solver.error() << std::endl;
}
-#endif
+
// copyBorderAll(_pressure, 0, _zRes);
// solve Poisson equation
- solvePressurePre(_pressure, _divergence, _obstacles);
+ solvePressurePre(_pressure, _divergence, _obstacles, b, A);
+ for(unsigned int i = 0; i < _xRes * _yRes * _zRes; i++)
{
+ float value = (_pressure[i] - b[i]);
+ if(value > 0.01)
+ printf("error: p: %f, b: %f\n", _pressure[i], b[i]);
+ }
+#if 0
+ {
float maxvalue = 0;
for(unsigned int i = 0; i < _xRes * _yRes * _zRes; i++)
{
@@ -964,7 +970,7 @@
}
// printf("Max pressure: %f, dx: %f\n", maxvalue, _dx);
}
-
+#endif
// DG TODO: check this function, for now this is done in the next function
// setObstaclePressure(_pressure, 0, _zRes);
@@ -978,13 +984,13 @@
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;
- float pC = _pressure[index]; // center
- float pR = _pressure[index + 1]; // right
- float pL = _pressure[index - 1]; // left
- float pT = _pressure[index + _slabSize]; // top
- float pB = _pressure[index - _slabSize]; // bottom
- float pU = _pressure[index + _xRes]; // Up
- float pD = _pressure[index - _xRes]; // Down
+ float pC = b[index]; // center
+ float pR = b[index + 1]; // right
+ float pL = b[index - 1]; // left
+ float pT = b[index + _slabSize]; // top
+ float pB = b[index - _slabSize]; // bottom
+ float pU = b[index + _xRes]; // Up
+ float pD = b[index - _xRes]; // Down
if(!_obstacles[index])
{
Modified: branches/smoke2/intern/smoke/intern/FLUID_3D.h
===================================================================
--- branches/smoke2/intern/smoke/intern/FLUID_3D.h 2012-05-31 10:38:11 UTC (rev 47268)
+++ branches/smoke2/intern/smoke/intern/FLUID_3D.h 2012-05-31 10:54:27 UTC (rev 47269)
@@ -43,6 +43,13 @@
using namespace BasicVector;
class WTURBULENCE;
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+using namespace Eigen;
+
+// Fluid index
+#define FINDEX(indexX, indexY, indexZ) ((indexZ) * _xRes * _yRes + (indexY) * _xRes + (indexX))
+
class FLUID_3D
{
public:
@@ -163,7 +170,8 @@
void diffuseHeat();
void solvePressure(float* field, float* b, unsigned char* skip);
- void solvePressurePre(float* field, float* b, unsigned char* skip);
+ void solvePressurePre(float* field, float* b, unsigned char* skip, VectorXf &myb, SparseMatrix<float,RowMajor> &A);
+
void solvePressureJacobian(float* p, float* d, unsigned char* ob);
void solveHeat(float* field, float* b, unsigned char* skip);
Modified: branches/smoke2/intern/smoke/intern/FLUID_3D_SOLVERS.cpp
===================================================================
--- branches/smoke2/intern/smoke/intern/FLUID_3D_SOLVERS.cpp 2012-05-31 10:38:11 UTC (rev 47268)
+++ branches/smoke2/intern/smoke/intern/FLUID_3D_SOLVERS.cpp 2012-05-31 10:54:27 UTC (rev 47269)
@@ -1,6 +1,6 @@
/** \file smoke/intern/FLUID_3D_SOLVERS.cpp
- * \ingroup smoke
- */
+* \ingroup smoke
+*/
//////////////////////////////////////////////////////////////////////
// This file is part of Wavelet Turbulence.
//
@@ -29,8 +29,11 @@
#include "FLUID_3D.h"
#include <cstring>
-#define SOLVER_ACCURACY 1e-06
+#define SOLVER_ACCURACY 1e-07
+#include <iostream>
+#include <fstream>
+
//////////////////////////////////////////////////////////////////////
// solve the heat equation with CG
//////////////////////////////////////////////////////////////////////
@@ -50,123 +53,123 @@
_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
- 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++)
- {
- // if the cell is a variable
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list