diff options
author | Daniel Genrich <daniel.genrich@gmx.net> | 2012-07-02 04:40:42 +0400 |
---|---|---|
committer | Daniel Genrich <daniel.genrich@gmx.net> | 2012-07-02 04:40:42 +0400 |
commit | ec66cd4910d24f31ac46ce9d9ae3cf4755c515bb (patch) | |
tree | a67ef7c567d40b50ef95c9a1a6bf1bef13d15ecb | |
parent | bb5e072320a8d4d519c0803dc6243aeff7f4b20c (diff) |
- Add CG helper for difficult setups
- Remove unsued code
-rw-r--r-- | intern/smoke/intern/FLUID_3D_SOLVERS.cpp | 102 |
1 files changed, 40 insertions, 62 deletions
diff --git a/intern/smoke/intern/FLUID_3D_SOLVERS.cpp b/intern/smoke/intern/FLUID_3D_SOLVERS.cpp index 10ee31de49d..98fc2c99c3e 100644 --- a/intern/smoke/intern/FLUID_3D_SOLVERS.cpp +++ b/intern/smoke/intern/FLUID_3D_SOLVERS.cpp @@ -307,9 +307,48 @@ void FLUID_3D::solvePressurePre(float* field, float* b, unsigned char* skip) SMOKE_LOOP { field[index] += alpha * _direction[index]; + } - _residual[index] -= alpha * _q[index]; + if(i % 50) + { + SMOKE_LOOP + { + // 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) ); + } + else + { + _residual[index] = 0.0f; + } + } + } + else + { + SMOKE_LOOP + { + _residual[index] -= alpha * _q[index]; + } + } + SMOKE_LOOP + { // Apply preconditioner precond_apply_diag(_h, _residual, _Precond, skip, index); // _h[index] = _Precond[index] * _residual[index]; @@ -317,7 +356,6 @@ void FLUID_3D::solvePressurePre(float* field, float* b, unsigned char* skip) tmp = _residual[index] * _h[index]; deltaNew += tmp; maxR = (tmp > maxR) ? tmp : maxR; - } @@ -341,63 +379,3 @@ void FLUID_3D::solvePressurePre(float* field, float* b, unsigned char* skip) if (_direction) delete[] _direction; if (_q) delete[] _q; } -#if 0 -void FLUID_3D::solvePressureJacobian(float* p, float* d, unsigned char* ob) -{ - float *_tmp = new float[_totalCells]; - unsigned int x, y, z; - - float maxdp = 0.0f; - - do - { - maxdp = 0.0f; - - size_t 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++) - { - float dC = d[index]; - float pC = p[index]; // center - - float pR = p[index + 1]; // right - float pL = p[index - 1]; // left - float pT = p[index + _slabSize]; // top - float pB = p[index - _slabSize]; // bottom - float pU = p[index + _xRes]; // Up - float pD = p[index - _xRes]; // Down - - if(ob[index + 1]) pR = pC; - if(ob[index - 1]) pL = pC; - if(ob[index + _slabSize]) pT = pC; - if(ob[index - _slabSize]) pB = pC; - if(ob[index + _xRes]) pU = pC; - if(ob[index - _xRes]) pD = pC; - - _tmp[index] = (pR + pL + pT + pB + pU + pD + dC) / 6.0f; - - if(ob[index]) - _tmp[index] = 0; - } - - for(unsigned int i = 0; i < _totalCells; i++) - { - float dp = _tmp[i] - p[i]; - - if(dp < 0) - dp *= -1.0f; - - if(dp > maxdp) - maxdp = dp; - - p[i] = _tmp[i]; - } - printf("maxdp: %f\n", maxdp); - } while(maxdp > SOLVER_ACCURACY); - - printf("\n"); - - if (_tmp) delete[] _tmp; -} -#endif |