diff options
21 files changed, 1184 insertions, 1157 deletions
diff --git a/intern/smoke/extern/smoke_API.h b/intern/smoke/extern/smoke_API.h index 1a3edce2344..5607df70cf3 100644 --- a/intern/smoke/extern/smoke_API.h +++ b/intern/smoke/extern/smoke_API.h @@ -63,11 +63,11 @@ void smoke_turbulence_free(struct WTURBULENCE *wt); void smoke_turbulence_step(struct WTURBULENCE *wt, struct FLUID_3D *fluid); float *smoke_turbulence_get_density(struct WTURBULENCE *wt); -void smoke_turbulence_get_res(struct WTURBULENCE *wt, unsigned int *res); +void smoke_turbulence_get_res(struct WTURBULENCE *wt, int *res); void smoke_turbulence_set_noise(struct WTURBULENCE *wt, int type); -void smoke_turbulence_initBlenderRNA(struct WTURBULENCE *wt, float *strength); +void smoke_initWaveletBlenderRNA(struct WTURBULENCE *wt, float *strength); -void smoke_turbulence_dissolve(struct WTURBULENCE *wt, int speed, int log); +void smoke_dissolve_wavelet(struct WTURBULENCE *wt, int speed, int log); // export void smoke_turbulence_export(struct WTURBULENCE *wt, float **dens, float **densold, float **tcu, float **tcv, float **tcw); diff --git a/intern/smoke/intern/FLUID_3D.cpp b/intern/smoke/intern/FLUID_3D.cpp index 89dd893198b..7574a0ec16e 100644 --- a/intern/smoke/intern/FLUID_3D.cpp +++ b/intern/smoke/intern/FLUID_3D.cpp @@ -27,9 +27,9 @@ #include <zlib.h> // boundary conditions of the fluid domain -#define DOMAIN_BC_FRONT 1 -#define DOMAIN_BC_TOP 0 -#define DOMAIN_BC_LEFT 1 +#define DOMAIN_BC_FRONT 0 // z +#define DOMAIN_BC_TOP 1 // y +#define DOMAIN_BC_LEFT 1 // x #define DOMAIN_BC_BACK DOMAIN_BC_FRONT #define DOMAIN_BC_BOTTOM DOMAIN_BC_TOP #define DOMAIN_BC_RIGHT DOMAIN_BC_LEFT @@ -111,47 +111,42 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dt) : } // set side obstacles - size_t index; - for (int y = 0; y < _yRes; y++) // z - for (int x = 0; x < _xRes; x++) - { - // front slab - index = x + y * _xRes; - if(DOMAIN_BC_BOTTOM==1) _obstacles[index] = 1; + int index; + for (int y = 0; y < _yRes; y++) + for (int x = 0; x < _xRes; x++) + { + // front slab + index = x + y * _xRes; + if(DOMAIN_BC_FRONT==1) _obstacles[index] = 1; - // back slab - index += _totalCells - _slabSize; - if(DOMAIN_BC_TOP==1) _obstacles[index] = 1; - } - for (int z = 0; z < _zRes; z++) // y - for (int x = 0; x < _xRes; x++) - { - // bottom slab - index = x + z * _slabSize; - if(DOMAIN_BC_FRONT==1) _obstacles[index] = 1; + // back slab + index += _totalCells - _slabSize; + if(DOMAIN_BC_BACK==1) _obstacles[index] = 1; + } - // top slab - index += _slabSize - _xRes; - if(DOMAIN_BC_BACK==1) _obstacles[index] = 1; - } - for (int z = 0; z < _zRes; z++) // x - for (int y = 0; y < _yRes; y++) - { - // left slab - index = y * _xRes + z * _slabSize; - if(DOMAIN_BC_LEFT==1) _obstacles[index] = 1; + for (int z = 0; z < _zRes; z++) + for (int x = 0; x < _xRes; x++) + { + // bottom slab + index = x + z * _slabSize; + if(DOMAIN_BC_BOTTOM==1) _obstacles[index] = 1; - // right slab - index += _xRes - 1; - if(DOMAIN_BC_RIGHT==1) _obstacles[index] = 1; - } + // top slab + index += _slabSize - _xRes; + if(DOMAIN_BC_TOP==1) _obstacles[index] = 1; + } - /* - SPHERE *obsSphere = NULL; - obsSphere = new SPHERE(0.375,0.5,0.375, 0.1); // for 4 to 3 domain - addObstacle(obsSphere); - delete obsSphere; - */ + for (int z = 0; z < _zRes; z++) + for (int y = 0; y < _yRes; y++) + { + // left slab + index = y * _xRes + z * _slabSize; + if(DOMAIN_BC_LEFT==1) _obstacles[index] = 1; + + // right slab + index += _xRes - 1; + if(DOMAIN_BC_RIGHT==1) _obstacles[index] = 1; + } } FLUID_3D::~FLUID_3D() @@ -191,7 +186,7 @@ void FLUID_3D::step() for (int i = 0; i < _totalCells; i++) { _xForce[i] = _yForce[i] = _zForce[i] = 0.0f; - _obstacles[i] &= ~2; + // _obstacles[i] &= ~2; } wipeBoundaries(); @@ -232,7 +227,8 @@ void FLUID_3D::step() _totalTime += _dt; _totalSteps++; - memset(_obstacles, 0, sizeof(unsigned char)*_xRes*_yRes*_zRes); + // todo xxx dg: only clear obstacles, not boundaries + // memset(_obstacles, 0, sizeof(unsigned char)*_xRes*_yRes*_zRes); } ////////////////////////////////////////////////////////////////////// @@ -270,7 +266,7 @@ void FLUID_3D::artificialDamping(float* field) { ////////////////////////////////////////////////////////////////////// void FLUID_3D::copyBorderAll(float* field) { - size_t index; + int index; for (int y = 0; y < _yRes; y++) for (int x = 0; x < _xRes; x++) { @@ -350,13 +346,13 @@ void FLUID_3D::project() // copy out the boundaries if(DOMAIN_BC_LEFT == 0) setNeumannX(_xVelocity, _res); - else setZeroX(_xVelocity, _res); + else setZeroX(_xVelocity, _res); - if(DOMAIN_BC_TOP == 0) setNeumannZ(_zVelocity, _res); - else setZeroZ(_zVelocity, _res); + if(DOMAIN_BC_TOP == 0) setNeumannY(_yVelocity, _res); + else setZeroY(_yVelocity, _res); - if(DOMAIN_BC_FRONT == 0) setNeumannY(_yVelocity, _res); - else setZeroY(_yVelocity, _res); + if(DOMAIN_BC_FRONT == 0) setNeumannZ(_zVelocity, _res); + else setZeroZ(_zVelocity, _res); // calculate divergence index = _slabSize + _xRes + 1; @@ -400,12 +396,16 @@ void FLUID_3D::project() for (y = 1; y < _yRes - 1; y++, index += 2) for (x = 1; x < _xRes - 1; x++, index++) { - if(!_obstacles[index]) + // 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; - } + }/* + else + { + _xVelocity[index] = _yVelocity[index] = _zVelocity[index] = 0.0f; + }*/ } if (_pressure) delete[] _pressure; @@ -669,13 +669,13 @@ void FLUID_3D::advectMacCormack() Vec3Int res = Vec3Int(_xRes,_yRes,_zRes); if(DOMAIN_BC_LEFT == 0) copyBorderX(_xVelocity, res); - else setZeroX(_xVelocity, res); + else setZeroX(_xVelocity, res); - if(DOMAIN_BC_TOP == 0) copyBorderZ(_zVelocity, res); - else setZeroZ(_zVelocity, res); + if(DOMAIN_BC_TOP == 0) copyBorderY(_yVelocity, res); + else setZeroY(_yVelocity, res); - if(DOMAIN_BC_FRONT == 0) copyBorderY(_yVelocity, res); - else setZeroY(_yVelocity, res); + if(DOMAIN_BC_FRONT == 0) copyBorderZ(_zVelocity, res); + else setZeroZ(_zVelocity, res); SWAP_POINTERS(_xVelocity, _xVelocityOld); SWAP_POINTERS(_yVelocity, _yVelocityOld); @@ -698,13 +698,13 @@ void FLUID_3D::advectMacCormack() advectFieldMacCormack(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _zVelocityOld, _zVelocity, t1,t2, res, _obstacles); if(DOMAIN_BC_LEFT == 0) copyBorderX(_xVelocity, res); - else setZeroX(_xVelocity, res); + else setZeroX(_xVelocity, res); - if(DOMAIN_BC_TOP == 0) copyBorderZ(_zVelocity, res); - else setZeroZ(_zVelocity, res); + if(DOMAIN_BC_TOP == 0) copyBorderY(_yVelocity, res); + else setZeroY(_yVelocity, res); - if(DOMAIN_BC_FRONT == 0) copyBorderY(_yVelocity, res); - else setZeroY(_yVelocity, res); + if(DOMAIN_BC_FRONT == 0) copyBorderZ(_zVelocity, res); + else setZeroZ(_zVelocity, res); setZeroBorder(_density, res); setZeroBorder(_heat, res); diff --git a/intern/smoke/intern/FLUID_3D_STATIC.cpp b/intern/smoke/intern/FLUID_3D_STATIC.cpp index 4474129beea..4909c071c3d 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]; - size_t index; + int index; for (int z = 0; z < res[2]; z++) for (int y = 0; y < res[1]; y++) { @@ -92,6 +92,18 @@ void FLUID_3D::setNeumannX(float* field, Vec3Int res) index += res[0] - 1; field[index] = field[index - 2]; } + + // fix, force top slab to only allow outwards flux + for (int y = 0; y < res[1]; y++) + for (int z = 0; z < res[2]; z++) + { + // top slab + int index = y * res[0] + z * slabSize; + index += res[0] - 1; + if(field[index]<0.) field[index] = 0.; + index -= 1; + if(field[index]<0.) field[index] = 0.; + } } ////////////////////////////////////////////////////////////////////// @@ -100,18 +112,31 @@ void FLUID_3D::setNeumannX(float* field, Vec3Int res) void FLUID_3D::setNeumannY(float* field, Vec3Int res) { const int slabSize = res[0] * res[1]; - size_t index; + int index; for (int z = 0; z < res[2]; z++) for (int x = 0; x < res[0]; x++) { - // front slab + // bottom slab index = x + z * slabSize; field[index] = field[index + 2 * res[0]]; - // back slab + // top slab index += slabSize - res[0]; field[index] = field[index - 2 * res[0]]; } + + // fix, force top slab to only allow outwards flux + for (int z = 0; z < res[2]; z++) + for (int x = 0; x < res[0]; x++) + { + // top slab + int index = x + z * slabSize; + index += slabSize - res[0]; + if(field[index]<0.) field[index] = 0.; + index -= res[0]; + if(field[index]<0.) field[index] = 0.; + } + } ////////////////////////////////////////////////////////////////////// @@ -121,15 +146,15 @@ void FLUID_3D::setNeumannZ(float* field, Vec3Int res) { const int slabSize = res[0] * res[1]; const int totalCells = res[0] * res[1] * res[2]; - size_t index; + int index; for (int y = 0; y < res[1]; y++) for (int x = 0; x < res[0]; x++) { - // bottom slab + // front slab index = x + y * res[0]; field[index] = field[index + 2 * slabSize]; - // top slab + // back slab index += totalCells - slabSize; field[index] = field[index - 2 * slabSize]; } @@ -139,11 +164,11 @@ void FLUID_3D::setNeumannZ(float* field, Vec3Int res) for (int x = 0; x < res[0]; x++) { // top slab - index = x + y * res[0]; + int index = x + y * res[0]; index += totalCells - slabSize; - if(field[index]<0.) field[index] = 0.0f; + if(field[index]<0.) field[index] = 0.; index -= slabSize; - if(field[index]<0.) field[index] = 0.0f; + if(field[index]<0.) field[index] = 0.; } } @@ -231,13 +256,14 @@ void FLUID_3D::copyBorderX(float* field, Vec3Int res) void FLUID_3D::copyBorderY(float* field, Vec3Int res) { const int slabSize = res[0] * res[1]; + const int totalCells = res[0] * res[1] * res[2]; int index; for (int z = 0; z < res[2]; z++) for (int x = 0; x < res[0]; x++) { // bottom slab index = x + z * slabSize; - field[index] = field[index + res[0]]; + field[index] = field[index + res[0]]; // top slab index += slabSize - res[0]; field[index] = field[index - res[0]]; @@ -253,7 +279,7 @@ void FLUID_3D::copyBorderZ(float* field, Vec3Int res) { // front slab index = x + y * res[0]; - field[index] = field[index + slabSize]; + field[index] = field[index + slabSize]; // back slab index += totalCells - slabSize; field[index] = field[index - slabSize]; @@ -269,18 +295,21 @@ void FLUID_3D::advectFieldSemiLagrange(const float dt, const float* velx, const const int xres = res[0]; const int yres = res[1]; const int zres = res[2]; + static int hits = 0; + static int total = 0; const int slabSize = res[0] * res[1]; // scale dt up to grid resolution -#if PARALLEL==1 -#pragma omp parallel for schedule(static) +#if PARALLEL==1 && !_WIN32 +#pragma omp parallel +#pragma omp for schedule(static) #endif for (int z = 0; z < zres; z++) for (int y = 0; y < yres; y++) for (int x = 0; x < xres; x++) { const int index = x + y * xres + z * xres*yres; - + // backtrace float xTrace = x - dt * velx[index]; float yTrace = y - dt * vely[index]; @@ -337,8 +366,8 @@ void FLUID_3D::advectFieldSemiLagrange(const float dt, const float* velx, const // // comments are the pseudocode from selle's paper ////////////////////////////////////////////////////////////////////// -void FLUID_3D::advectFieldMacCormack(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity, - float* oldField, float* newField, float* temp1, float* temp2, Vec3Int res, const unsigned char* obstacles) +void FLUID_3D::advectFieldMacCormack(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity, + float* oldField, float* newField, float* temp1, float* temp2, Vec3Int res, const unsigned char* obstacles) { float* phiHatN = temp1; float* phiHatN1 = temp2; @@ -359,7 +388,7 @@ void FLUID_3D::advectFieldMacCormack(const float dt, const float* xVelocity, con advectFieldSemiLagrange( -1.0*dt, xVelocity, yVelocity, zVelocity, phiHatN1, phiHatN, res); // phiN1 = phiHatN1 + (phiN - phiHatN) / 2 - const int border = 0; + const int border = 0; for (int z = border; z < sz-border; z++) for (int y = border; y < sy-border; y++) for (int x = border; x < sx-border; x++) { @@ -376,7 +405,7 @@ void FLUID_3D::advectFieldMacCormack(const float dt, const float* xVelocity, con // if the error estimate was bad, revert to first order clampOutsideRays(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res, obstacles, phiHatN1); -} +} ////////////////////////////////////////////////////////////////////// @@ -454,17 +483,18 @@ void FLUID_3D::clampExtrema(const float dt, const float* velx, const float* vely } ////////////////////////////////////////////////////////////////////// -// Reverts any backtraces that go into boundaries back to first +// Reverts any backtraces that go into boundaries back to first // order -- in this case the error correction term was totally // incorrect ////////////////////////////////////////////////////////////////////// -void FLUID_3D::clampOutsideRays(const float dt, const float* velx, const float* vely, const float* velz, - float* oldField, float* newField, Vec3Int res, const unsigned char* obstacles, const float *oldAdvection) +void FLUID_3D::clampOutsideRays(const float dt, const float* velx, const float* vely, const float* velz, + float* oldField, float* newField, Vec3Int res, const unsigned char* obstacles, const float *oldAdvection) { const int sx= res[0]; const int sy= res[1]; const int sz= res[2]; const int slabSize = res[0] * res[1]; + for (int z = 1; z < sz-1; z++) for (int y = 1; y < sy-1; y++) for (int x = 1; x < sx-1; x++) @@ -479,7 +509,7 @@ void FLUID_3D::clampOutsideRays(const float dt, const float* velx, const float* float zTrace = z - dt * velz[index]; // see if it goes outside the boundaries - bool hasObstacle = + bool hasObstacle = (zTrace < 1.0f) || (zTrace > sz - 2.0f) || (yTrace < 1.0f) || (yTrace > sy - 2.0f) || (xTrace < 1.0f) || (xTrace > sx - 2.0f) || @@ -515,7 +545,7 @@ void FLUID_3D::clampOutsideRays(const float dt, const float* velx, const float* int z0 = (int)zBackward; int z1 = z0 + 1; if(obstacles && !hasObstacle) { - hasObstacle = hasObstacle || + hasObstacle = hasObstacle || obstacles[x0 + y0 * sx + z0*slabSize] || obstacles[x0 + y1 * sx + z0*slabSize] || obstacles[x1 + y0 * sx + z0*slabSize] || @@ -535,7 +565,7 @@ void FLUID_3D::clampOutsideRays(const float dt, const float* velx, const float* z0 = (int)zTrace; z1 = z0 + 1; if(obstacles && !hasObstacle) { - hasObstacle = hasObstacle || + hasObstacle = hasObstacle || obstacles[x0 + y0 * sx + z0*slabSize] || obstacles[x0 + y1 * sx + z0*slabSize] || obstacles[x1 + y0 * sx + z0*slabSize] || @@ -577,84 +607,8 @@ void FLUID_3D::clampOutsideRays(const float dt, const float* velx, const float* u1 * (s0 * (t0 * oldField[i001] + t1 * oldField[i011]) + s1 * (t0 * oldField[i101] + - t1 * oldField[i111])); + t1 * oldField[i111])); } } // xyz } -////////////////////////////////////////////////////////////////////// -// image output -////////////////////////////////////////////////////////////////////// -/* -void FLUID_3D::writeImageSliceXY(const float *field, Vec3Int res, int slice, string prefix, int picCnt, float scale) { - writeProjectedIntern(field, res, 0,1, prefix, picCnt, scale); -} -void FLUID_3D::writeImageSliceYZ(const float *field, Vec3Int res, int slice, string prefix, int picCnt, float scale) { - writeProjectedIntern(field, res, 1,2, prefix, picCnt, scale); -} -void FLUID_3D::writeImageSliceXZ(const float *field, Vec3Int res, int slice, string prefix, int picCnt, float scale) { - writeProjectedIntern(field, res, 0,2, prefix, picCnt, scale); -} -*/ - -////////////////////////////////////////////////////////////////////// -// Helper function for projecting densities along a dimension -////////////////////////////////////////////////////////////////////// -/* -static int getOtherDir(int dir1, int dir2) { - switch(dir1) { - case 0: - switch(dir2) { - case 1: return 2; - case 2: return 1; } - break; - case 1: - switch(dir2) { - case 0: return 2; - case 2: return 0; } - break; - case 2: - switch(dir2) { - case 0: return 1; - case 1: return 0; } - break; - default: - return 0; - } - return 0; -} -*/ - -////////////////////////////////////////////////////////////////////// -// average densities along third spatial direction -////////////////////////////////////////////////////////////////////// -/* -void FLUID_3D::writeProjectedIntern(const float *field, Vec3Int res, - int dir1, int dir2, string prefix, int picCnt, float scale) { - const int nitems = res[dir1]*res[dir2]; - const int otherDir = getOtherDir(dir1,dir2); - float *buf = new float[nitems]; - Vec3Int min = Vec3Int(0); - Vec3Int max = res; - - min[otherDir] = 0; - max[otherDir] = res[otherDir]; - float div = 1./(float)MIN3V(res); // normalize for shorter sides, old: res[otherDir]; - div *= 4.; //slightly increase contrast - for(int i=0; i<nitems; i++) buf[i] = 0.; - - Vec3Int cnt = 0; - for (cnt[2] = min[2]; cnt[2] < max[2]; cnt[2]++) { - for (cnt[1] = min[1]; cnt[1] < max[1]; cnt[1]++) - for (cnt[0] = min[0]; cnt[0] < max[0]; cnt[0]++) - { - const int index = cnt[0] + cnt[1] * res[0] + cnt[2] * res[0]*res[1]; - const int bufindex = cnt[dir1] + cnt[dir2] * res[dir1]; - buf[bufindex] += field[index] * scale *div; - } - } - // IMAGE::dumpNumberedPNG(picCnt, prefix, buf, res[dir1], res[dir2]); - delete[] buf; -} -*/ - diff --git a/intern/smoke/intern/WTURBULENCE.cpp b/intern/smoke/intern/WTURBULENCE.cpp index dd092d4f0cc..a1b2aaf30f2 100644 --- a/intern/smoke/intern/WTURBULENCE.cpp +++ b/intern/smoke/intern/WTURBULENCE.cpp @@ -92,10 +92,6 @@ WTURBULENCE::WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify, int no _tcW = new float[_totalCellsSm]; _tcTemp = new float[_totalCellsSm]; - // allocate & init energy terms - _energy = new float[_totalCellsSm]; - _highFreqEnergy = new float[_totalCellsSm]; - // map all const float dx = 1./(float)(_resSm[0]); const float dy = 1./(float)(_resSm[1]); @@ -109,15 +105,8 @@ WTURBULENCE::WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify, int no _tcV[index] = y*dy; _tcW[index] = z*dz; _tcTemp[index] = 0.; - _energy[index] = 0.; } - // allocate eigenvalue arrays - _eigMin = new float[_totalCellsSm]; - _eigMax = new float[_totalCellsSm]; - for(int i=0; i < _totalCellsSm; i++) - _eigMin[i] = _eigMax[i] = 0.; - // noise tiles _noiseTile = new float[noiseTileSize * noiseTileSize * noiseTileSize]; /* @@ -143,12 +132,7 @@ WTURBULENCE::~WTURBULENCE() { delete[] _tcW; delete[] _tcTemp; - delete[] _eigMin; - delete[] _eigMax; delete[] _noiseTile; - - delete[] _energy; - delete[] _highFreqEnergy; } ////////////////////////////////////////////////////////////////////// @@ -293,34 +277,34 @@ static float minDz(int x, int y, int z, float* input, Vec3Int res) // handle texture coordinates (advection, reset, eigenvalues), // Beware -- uses big density maccormack as temporary arrays ////////////////////////////////////////////////////////////////////// -void WTURBULENCE::advectTextureCoordinates (float dtOrg, float* xvel, float* yvel, float* zvel, float *_tempBig1, float *_tempBig2) { +void WTURBULENCE::advectTextureCoordinates (float dtOrg, float* xvel, float* yvel, float* zvel, float *tempBig1, float *tempBig2) { // advection SWAP_POINTERS(_tcTemp, _tcU); FLUID_3D::copyBorderX(_tcTemp, _resSm); FLUID_3D::copyBorderY(_tcTemp, _resSm); FLUID_3D::copyBorderZ(_tcTemp, _resSm); FLUID_3D::advectFieldMacCormack(dtOrg, xvel, yvel, zvel, - _tcTemp, _tcU, _tempBig1, _tempBig2, _resSm, NULL); + _tcTemp, _tcU, tempBig1, tempBig2, _resSm, NULL); SWAP_POINTERS(_tcTemp, _tcV); FLUID_3D::copyBorderX(_tcTemp, _resSm); FLUID_3D::copyBorderY(_tcTemp, _resSm); FLUID_3D::copyBorderZ(_tcTemp, _resSm); FLUID_3D::advectFieldMacCormack(dtOrg, xvel, yvel, zvel, - _tcTemp, _tcV, _tempBig1, _tempBig2, _resSm, NULL); + _tcTemp, _tcV, tempBig1, tempBig2, _resSm, NULL); SWAP_POINTERS(_tcTemp, _tcW); FLUID_3D::copyBorderX(_tcTemp, _resSm); FLUID_3D::copyBorderY(_tcTemp, _resSm); FLUID_3D::copyBorderZ(_tcTemp, _resSm); FLUID_3D::advectFieldMacCormack(dtOrg, xvel, yvel, zvel, - _tcTemp, _tcW, _tempBig1, _tempBig2, _resSm, NULL); + _tcTemp, _tcW, tempBig1, tempBig2, _resSm, NULL); } ////////////////////////////////////////////////////////////////////// // Compute the eigenvalues of the advected texture ////////////////////////////////////////////////////////////////////// -void WTURBULENCE::computeEigenvalues() { +void WTURBULENCE::computeEigenvalues(float *_eigMin, float *_eigMax) { // stats float maxeig = -1.; float mineig = 10.; @@ -365,7 +349,7 @@ void WTURBULENCE::computeEigenvalues() { ////////////////////////////////////////////////////////////////////// // advect & reset texture coordinates based on eigenvalues ////////////////////////////////////////////////////////////////////// -void WTURBULENCE::resetTextureCoordinates() +void WTURBULENCE::resetTextureCoordinates(float *_eigMin, float *_eigMax) { // allowed deformation of the textures const float limit = 2.f; @@ -396,7 +380,7 @@ void WTURBULENCE::resetTextureCoordinates() // Compute the highest frequency component of the wavelet // decomposition ////////////////////////////////////////////////////////////////////// -void WTURBULENCE::decomposeEnergy() +void WTURBULENCE::decomposeEnergy(float *_energy, float *_highFreqEnergy) { // do the decomposition -- the goal here is to have // the energy with the high frequency component stomped out @@ -432,7 +416,7 @@ void WTURBULENCE::decomposeEnergy() // compute velocity from energies and march into obstacles // for wavelet decomposition ////////////////////////////////////////////////////////////////////// -void WTURBULENCE::computeEnergy(float* xvel, float* yvel, float* zvel, unsigned char *obstacles) +void WTURBULENCE::computeEnergy(float *_energy, float* xvel, float* yvel, float* zvel, unsigned char *obstacles) { // compute everywhere for (int x = 0; x < _totalCellsSm; x++) @@ -491,7 +475,7 @@ void WTURBULENCE::computeEnergy(float* xvel, float* yvel, float* zvel, unsigned } if (valid > 0) { - _energy[index] = sum / valid; + _energy[index] = sum / (float)valid; obstacles[index] = MARCHED; } } @@ -516,9 +500,9 @@ void WTURBULENCE::computeEnergy(float* xvel, float* yvel, float* zvel, unsigned Vec3 WTURBULENCE::WVelocity(Vec3 orgPos) { // arbitrarily offset evaluation points - const Vec3 p1 = orgPos + Vec3(NOISE_TILE_SIZE/2,0,0); - const Vec3 p2 = orgPos + Vec3(0,NOISE_TILE_SIZE/2,0); - const Vec3 p3 = orgPos + Vec3(0,0,NOISE_TILE_SIZE/2); + const Vec3 p1 = orgPos + Vec3(NOISE_TILE_SIZE/2.0,0,0); + const Vec3 p2 = orgPos + Vec3(0,NOISE_TILE_SIZE/2.0,0); + const Vec3 p3 = orgPos + Vec3(0,0,NOISE_TILE_SIZE/2.0); const float f1y = WNoiseDy(p1, _noiseTile); const float f1z = WNoiseDz(p1, _noiseTile); @@ -542,9 +526,9 @@ Vec3 WTURBULENCE::WVelocity(Vec3 orgPos) Vec3 WTURBULENCE::WVelocityWithJacobian(Vec3 orgPos, float* xUnwarped, float* yUnwarped, float* zUnwarped) { // arbitrarily offset evaluation points - const Vec3 p1 = orgPos + Vec3(NOISE_TILE_SIZE/2,0,0); - const Vec3 p2 = orgPos + Vec3(0,NOISE_TILE_SIZE/2,0); - const Vec3 p3 = orgPos + Vec3(0,0,NOISE_TILE_SIZE/2); + const Vec3 p1 = orgPos + Vec3(NOISE_TILE_SIZE/2.0,0,0); + const Vec3 p2 = orgPos + Vec3(0,NOISE_TILE_SIZE/2.0,0); + const Vec3 p3 = orgPos + Vec3(0,0,NOISE_TILE_SIZE/2.0); Vec3 final; final[0] = WNoiseDx(p1, _noiseTile); @@ -575,44 +559,40 @@ Vec3 WTURBULENCE::WVelocityWithJacobian(Vec3 orgPos, float* xUnwarped, float* yU return ret; } + ////////////////////////////////////////////////////////////////////// // perform an actual noise advection step ////////////////////////////////////////////////////////////////////// void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel, float* zvel, unsigned char *obstacles) { - // big velocity macCormack fields - float* _bigUx; - float* _bigUy; - float* _bigUz; - - // temp arrays for BFECC and MacCormack - they have more convenient - // names in the actual implementations - float* _tempBig1; - float* _tempBig2; - - // allocate high resolution velocity field. Note that this is only - // necessary because we use MacCormack advection. For semi-Lagrangian - // advection, these arrays are not necessary. - _tempBig1 = new float[_totalCellsBig]; - _tempBig2 = new float[_totalCellsBig]; - - // enlarge timestep to match grid - const float dt = dtOrg * _amplify; - const float invAmp = 1.0f / _amplify; - - _bigUx = new float[_totalCellsBig]; - _bigUy = new float[_totalCellsBig]; - _bigUz = new float[_totalCellsBig]; - - // prepare textures - advectTextureCoordinates(dtOrg, xvel,yvel,zvel, _tempBig1, _tempBig2); + // enlarge timestep to match grid + const float dt = dtOrg * _amplify; + const float invAmp = 1.0f / _amplify; + float *tempBig1 = new float[_totalCellsBig]; + float *tempBig2 = new float[_totalCellsBig]; + float *bigUx = new float[_totalCellsBig]; + float *bigUy = new float[_totalCellsBig]; + float *bigUz = new float[_totalCellsBig]; + float *_energy = new float[_totalCellsSm]; + float *highFreqEnergy = new float[_totalCellsSm]; + float *eigMin = new float[_totalCellsSm]; + float *eigMax = new float[_totalCellsSm]; + + memset(tempBig1, 0, sizeof(float)*_totalCellsBig); + memset(tempBig2, 0, sizeof(float)*_totalCellsBig); + memset(highFreqEnergy, 0, sizeof(float)*_totalCellsSm); + memset(eigMin, 0, sizeof(float)*_totalCellsSm); + memset(eigMax, 0, sizeof(float)*_totalCellsSm); + + // prepare textures + advectTextureCoordinates(dtOrg, xvel,yvel,zvel, tempBig1, tempBig2); // compute eigenvalues of the texture coordinates - computeEigenvalues(); + computeEigenvalues(eigMin, eigMax); // do wavelet decomposition of energy - computeEnergy(xvel, yvel, zvel, obstacles); - decomposeEnergy(); + computeEnergy(_energy, xvel, yvel, zvel, obstacles); + decomposeEnergy(_energy, highFreqEnergy); // zero out coefficients inside of the obstacle for (int x = 0; x < _totalCellsSm; x++) @@ -647,7 +627,7 @@ void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel, // retrieve wavelet energy at highest frequency float energy = INTERPOLATE::lerp3d( - _highFreqEnergy, posSm[0],posSm[1],posSm[2], _xResSm, _yResSm, _zResSm); + highFreqEnergy, posSm[0],posSm[1],posSm[2], _xResSm, _yResSm, _zResSm); // base amplitude for octave 0 float coefficient = sqrtf(2.0f * fabs(energy)); @@ -656,8 +636,8 @@ void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel, // add noise to velocity, but only if the turbulence is // sufficiently undeformed, and the energy is large enough // to make a difference - const bool addNoise = _eigMax[indexSmall] < 2. && - _eigMin[indexSmall] > 0.5; + const bool addNoise = eigMax[indexSmall] < 2. && + eigMin[indexSmall] > 0.5; if (addNoise && amplitude > _cullingThreshold) { // base amplitude for octave 0 float amplitudeScaled = amplitude; @@ -679,21 +659,21 @@ void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel, // If you wanted to save memory, you would instead perform a // semi-Lagrangian backtrace for the current grid cell here. Then // you could just throw the velocity away. - _bigUx[index] = vel[0]; - _bigUy[index] = vel[1]; - _bigUz[index] = vel[2]; + bigUx[index] = vel[0]; + bigUy[index] = vel[1]; + bigUz[index] = vel[2]; // compute the velocity magnitude for substepping later - const float velMag = _bigUx[index] * _bigUx[index] + - _bigUy[index] * _bigUy[index] + - _bigUz[index] * _bigUz[index]; + const float velMag = bigUx[index] * bigUx[index] + + bigUy[index] * bigUy[index] + + bigUz[index] * bigUz[index]; if (velMag > maxVelocity) maxVelocity = velMag; // zero out velocity inside obstacles float obsCheck = INTERPOLATE::lerp3dToFloat( obstacles, posSm[0], posSm[1], posSm[2], _xResSm, _yResSm, _zResSm); if (obsCheck > 0.95) - _bigUx[index] = _bigUy[index] = _bigUz[index] = 0.; + bigUx[index] = bigUy[index] = bigUz[index] = 0.; } // prepare density for an advection @@ -701,24 +681,23 @@ void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel, // based on the maximum velocity present, see if we need to substep, // but cap the maximum number of substeps to 5 - const int maxSubSteps = 25; - const int maxVel = 5; + const int maxSubSteps = 5; maxVelocity = sqrt(maxVelocity) * dt; - int totalSubsteps = (int)(maxVelocity / (float)maxVel); + int totalSubsteps = (int)(maxVelocity / (float)maxSubSteps); totalSubsteps = (totalSubsteps < 1) ? 1 : totalSubsteps; totalSubsteps = (totalSubsteps > maxSubSteps) ? maxSubSteps : totalSubsteps; const float dtSubdiv = dt / (float)totalSubsteps; // set boundaries of big velocity grid - FLUID_3D::setZeroX(_bigUx, _resBig); - FLUID_3D::setZeroY(_bigUy, _resBig); - FLUID_3D::setZeroZ(_bigUz, _resBig); + FLUID_3D::setZeroX(bigUx, _resBig); + FLUID_3D::setZeroY(bigUy, _resBig); + FLUID_3D::setZeroZ(bigUz, _resBig); // do the MacCormack advection, with substepping if necessary for(int substep = 0; substep < totalSubsteps; substep++) { - FLUID_3D::advectFieldMacCormack(dtSubdiv, _bigUx, _bigUy, _bigUz, - _densityBigOld, _densityBig, _tempBig1, _tempBig2, _resBig, NULL); + FLUID_3D::advectFieldMacCormack(dtSubdiv, bigUx, bigUy, bigUz, + _densityBigOld, _densityBig, tempBig1, tempBig2, _resBig, NULL); if (substep < totalSubsteps - 1) SWAP_POINTERS(_densityBig, _densityBigOld); @@ -730,25 +709,21 @@ void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel, // reset texture coordinates now in preparation for next timestep // Shouldn't do this before generating the noise because then the // eigenvalues stored do not reflect the underlying texture coordinates - resetTextureCoordinates(); + resetTextureCoordinates(eigMin, eigMax); - // output files - /* - string prefix = string("./amplified.preview/density_bigxy_"); - FLUID_3D::writeImageSliceXY(_densityBig, _resBig, _resBig[2]/2, prefix, _totalStepsBig, 1.0f); - //string df3Prefix = string("./df3/density_big_"); - //IMAGE::dumpDF3(_totalStepsBig, df3Prefix, _densityBig, _resBig[0],_resBig[1],_resBig[2]); - string pbrtPrefix = string("./pbrt/density_big_"); - IMAGE::dumpPBRT(_totalStepsBig, pbrtPrefix, _densityBig, _resBig[0],_resBig[1],_resBig[2]); - */ - _totalStepsBig++; + delete[] tempBig1; + delete[] tempBig2; + delete[] bigUx; + delete[] bigUy; + delete[] bigUz; + delete[] _energy; + delete[] highFreqEnergy; - delete[] _bigUx; - delete[] _bigUy; - delete[] _bigUz; + delete[] eigMin; + delete[] eigMax; + - delete[] _tempBig1; - delete[] _tempBig2; + _totalStepsBig++; } ////////////////////////////////////////////////////////////////////// @@ -757,257 +732,258 @@ void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel, ////////////////////////////////////////////////////////////////////// void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, float* zvel, unsigned char *obstacles) { - // big velocity macCormack fields - float* _bigUx; - float* _bigUy; - float* _bigUz; - - // temp arrays for BFECC and MacCormack - they have more convenient - // names in the actual implementations - float* _tempBig1; - float* _tempBig2; - - // allocate high resolution velocity field. Note that this is only - // necessary because we use MacCormack advection. For semi-Lagrangian - // advection, these arrays are not necessary. - _tempBig1 = new float[_totalCellsBig]; - _tempBig2 = new float[_totalCellsBig]; - // enlarge timestep to match grid const float dt = dtOrg * _amplify; const float invAmp = 1.0f / _amplify; - - _bigUx = new float[_totalCellsBig]; - _bigUy = new float[_totalCellsBig]; - _bigUz = new float[_totalCellsBig]; + float *tempBig1 = new float[_totalCellsBig]; + float *tempBig2 = new float[_totalCellsBig]; + float *bigUx = new float[_totalCellsBig]; + float *bigUy = new float[_totalCellsBig]; + float *bigUz = new float[_totalCellsBig]; + float *_energy = new float[_totalCellsSm]; + float *highFreqEnergy = new float[_totalCellsSm]; + float *eigMin = new float[_totalCellsSm]; + float *eigMax = new float[_totalCellsSm]; + + memset(highFreqEnergy, 0, sizeof(float)*_totalCellsSm); + memset(eigMin, 0, sizeof(float)*_totalCellsSm); + memset(eigMax, 0, sizeof(float)*_totalCellsSm); // prepare textures - advectTextureCoordinates(dtOrg, xvel,yvel,zvel, _tempBig1, _tempBig2); + advectTextureCoordinates(dtOrg, xvel,yvel,zvel, tempBig1, tempBig2); // do wavelet decomposition of energy - computeEnergy(xvel, yvel, zvel, obstacles); - decomposeEnergy(); + computeEnergy(_energy, xvel, yvel, zvel, obstacles); - // zero out coefficients inside of the obstacle for (int x = 0; x < _totalCellsSm; x++) if (obstacles[x]) _energy[x] = 0.f; - // parallel region setup - float maxVelMagThreads[8] = { -1., -1., -1., -1., -1., -1., -1., -1. }; - - #if PARALLEL==1 - #pragma omp parallel - #endif - { float maxVelMag1 = 0.; - #if PARALLEL==1 - const int id = omp_get_thread_num(); /*, num = omp_get_num_threads(); */ - #endif - - // vector noise main loop - #if PARALLEL==1 - #pragma omp for schedule(static) - #endif - for (int zSmall = 0; zSmall < _zResSm; zSmall++) - for (int ySmall = 0; ySmall < _yResSm; ySmall++) - for (int xSmall = 0; xSmall < _xResSm; xSmall++) - { - const int indexSmall = xSmall + ySmall * _xResSm + zSmall * _slabSizeSm; - - // compute jacobian - float jacobian[3][3] = { - { minDx(xSmall, ySmall, zSmall, _tcU, _resSm), minDx(xSmall, ySmall, zSmall, _tcV, _resSm), minDx(xSmall, ySmall, zSmall, _tcW, _resSm) } , - { minDy(xSmall, ySmall, zSmall, _tcU, _resSm), minDy(xSmall, ySmall, zSmall, _tcV, _resSm), minDy(xSmall, ySmall, zSmall, _tcW, _resSm) } , - { minDz(xSmall, ySmall, zSmall, _tcU, _resSm), minDz(xSmall, ySmall, zSmall, _tcV, _resSm), minDz(xSmall, ySmall, zSmall, _tcW, _resSm) } - }; - - // get LU factorization of texture jacobian and apply - // it to unit vectors - JAMA::LU<float> LU = computeLU3x3(jacobian); - float xUnwarped[] = {1.0f, 0.0f, 0.0f}; - float yUnwarped[] = {0.0f, 1.0f, 0.0f}; - float zUnwarped[] = {0.0f, 0.0f, 1.0f}; - float xWarped[] = {1.0f, 0.0f, 0.0f}; - float yWarped[] = {0.0f, 1.0f, 0.0f}; - float zWarped[] = {0.0f, 0.0f, 1.0f}; - bool nonSingular = LU.isNonsingular(); - - #if 0 - // UNUSED - float eigMax = 10.0f; - float eigMin = 0.1f; - #endif + decomposeEnergy(_energy, highFreqEnergy); - if (nonSingular) - { - solveLU3x3(LU, xUnwarped, xWarped); - solveLU3x3(LU, yUnwarped, yWarped); - solveLU3x3(LU, zUnwarped, zWarped); - - // compute the eigenvalues while we have the Jacobian available - Vec3 eigenvalues = Vec3(1.); - computeEigenvalues3x3( &eigenvalues[0], jacobian); - _eigMax[indexSmall] = MAX3V(eigenvalues); - _eigMin[indexSmall] = MIN3V(eigenvalues); - } + // zero out coefficients inside of the obstacle + for (int x = 0; x < _totalCellsSm; x++) + if (obstacles[x]) highFreqEnergy[x] = 0.f; - // make sure to skip one on the beginning and end - int xStart = (xSmall == 0) ? 1 : 0; - int xEnd = (xSmall == _xResSm - 1) ? _amplify - 1 : _amplify; - int yStart = (ySmall == 0) ? 1 : 0; - int yEnd = (ySmall == _yResSm - 1) ? _amplify - 1 : _amplify; - int zStart = (zSmall == 0) ? 1 : 0; - int zEnd = (zSmall == _zResSm - 1) ? _amplify - 1 : _amplify; - - for (int zBig = zStart; zBig < zEnd; zBig++) - for (int yBig = yStart; yBig < yEnd; yBig++) - for (int xBig = xStart; xBig < xEnd; xBig++) - { - const int x = xSmall * _amplify + xBig; - const int y = ySmall * _amplify + yBig; - const int z = zSmall * _amplify + zBig; - - // get unit position for both fine and coarse grid - const Vec3 pos = Vec3(x,y,z); - const Vec3 posSm = pos * invAmp; - - // get grid index for both fine and coarse grid - const int index = x + y *_xResBig + z *_slabSizeBig; - - // get a linearly interpolated velocity and texcoords - // from the coarse grid - Vec3 vel = INTERPOLATE::lerp3dVec( xvel,yvel,zvel, - posSm[0], posSm[1], posSm[2], _xResSm,_yResSm,_zResSm); - Vec3 uvw = INTERPOLATE::lerp3dVec( _tcU,_tcV,_tcW, - posSm[0], posSm[1], posSm[2], _xResSm,_yResSm,_zResSm); - - // multiply the texture coordinate by _resSm so that turbulence - // synthesis begins at the first octave that the coarse grid - // cannot capture - Vec3 texCoord = Vec3(uvw[0] * _resSm[0], - uvw[1] * _resSm[1], - uvw[2] * _resSm[2]); - - // retrieve wavelet energy at highest frequency - float energy = INTERPOLATE::lerp3d( - _highFreqEnergy, posSm[0],posSm[1],posSm[2], _xResSm, _yResSm, _zResSm); - - // base amplitude for octave 0 - float coefficient = sqrtf(2.0f * fabs(energy)); - 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 - // to make a difference - const bool addNoise = _eigMax[indexSmall] < 2. && - _eigMin[indexSmall] > 0.5; - - if (addNoise && amplitude > _cullingThreshold) { - // base amplitude for octave 0 - float amplitudeScaled = amplitude; - - for (int octave = 0; octave < _octaves; octave++) - { - // multiply the vector noise times the maximum allowed - // noise amplitude at this octave, and add it to the total - vel += WVelocityWithJacobian(texCoord, &xUnwarped[0], &yUnwarped[0], &zUnwarped[0]) * amplitudeScaled; - - // scale coefficient for next octave - amplitudeScaled *= persistence; - texCoord *= 2.0f; - } - } - - // Store velocity + turbulence in big grid for maccormack step - // - // If you wanted to save memory, you would instead perform a - // semi-Lagrangian backtrace for the current grid cell here. Then - // you could just throw the velocity away. - _bigUx[index] = vel[0]; - _bigUy[index] = vel[1]; - _bigUz[index] = vel[2]; - - // compute the velocity magnitude for substepping later - const float velMag = _bigUx[index] * _bigUx[index] + - _bigUy[index] * _bigUy[index] + - _bigUz[index] * _bigUz[index]; - if (velMag > maxVelMag1) maxVelMag1 = velMag; - - // zero out velocity inside obstacles - float obsCheck = INTERPOLATE::lerp3dToFloat( - obstacles, posSm[0], posSm[1], posSm[2], _xResSm, _yResSm, _zResSm); - - if (obsCheck > 0.95) - _bigUx[index] = _bigUy[index] = _bigUz[index] = 0.; - } // xyz - - #if PARALLEL==1 - maxVelMagThreads[id] = maxVelMag1; - #else - maxVelMagThreads[0] = maxVelMag1; - #endif - } - } // omp - - // compute maximum over threads - float maxVelMag = maxVelMagThreads[0]; - #if PARALLEL==1 - for (int i = 1; i < 8; i++) - if (maxVelMag < maxVelMagThreads[i]) - maxVelMag = maxVelMagThreads[i]; - #endif - - // prepare density for an advection - SWAP_POINTERS(_densityBig, _densityBigOld); - - // based on the maximum velocity present, see if we need to substep, - // but cap the maximum number of substeps to 5 - const int maxSubSteps = 25; - const int maxVel = 5; - 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; - - // set boundaries of big velocity grid - FLUID_3D::setZeroX(_bigUx, _resBig); - FLUID_3D::setZeroY(_bigUy, _resBig); - FLUID_3D::setZeroZ(_bigUz, _resBig); - - // do the MacCormack advection, with substepping if necessary - for(int substep = 0; substep < totalSubsteps; substep++) - { - FLUID_3D::advectFieldMacCormack(dtSubdiv, _bigUx, _bigUy, _bigUz, - _densityBigOld, _densityBig, _tempBig1, _tempBig2, _resBig, NULL); + Vec3Int ressm(_xResSm, _yResSm, _zResSm); + FLUID_3D::setNeumannX(highFreqEnergy, ressm); + FLUID_3D::setNeumannY(highFreqEnergy, ressm); + FLUID_3D::setNeumannZ(highFreqEnergy, ressm); + + // parallel region setup + float maxVelMagThreads[8] = { -1., -1., -1., -1., -1., -1., -1., -1. }; +#if PARALLEL==1 +#pragma omp parallel +#endif + { float maxVelMag1 = 0.; +#if PARALLEL==1 + const int id = omp_get_thread_num(); /*, num = omp_get_num_threads(); */ +#endif + + // vector noise main loop +#if PARALLEL==1 +#pragma omp for schedule(static) +#endif + for (int zSmall = 0; zSmall < _zResSm; zSmall++) + for (int ySmall = 0; ySmall < _yResSm; ySmall++) + for (int xSmall = 0; xSmall < _xResSm; xSmall++) + { + const int indexSmall = xSmall + ySmall * _xResSm + zSmall * _slabSizeSm; + + // compute jacobian + float jacobian[3][3] = { + { minDx(xSmall, ySmall, zSmall, _tcU, _resSm), minDx(xSmall, ySmall, zSmall, _tcV, _resSm), minDx(xSmall, ySmall, zSmall, _tcW, _resSm) } , + { minDy(xSmall, ySmall, zSmall, _tcU, _resSm), minDy(xSmall, ySmall, zSmall, _tcV, _resSm), minDy(xSmall, ySmall, zSmall, _tcW, _resSm) } , + { minDz(xSmall, ySmall, zSmall, _tcU, _resSm), minDz(xSmall, ySmall, zSmall, _tcV, _resSm), minDz(xSmall, ySmall, zSmall, _tcW, _resSm) } + }; + + // get LU factorization of texture jacobian and apply + // it to unit vectors + JAMA::LU<float> LU = computeLU3x3(jacobian); + float xUnwarped[] = {1.0f, 0.0f, 0.0f}; + float yUnwarped[] = {0.0f, 1.0f, 0.0f}; + float zUnwarped[] = {0.0f, 0.0f, 1.0f}; + float xWarped[] = {1.0f, 0.0f, 0.0f}; + float yWarped[] = {0.0f, 1.0f, 0.0f}; + float zWarped[] = {0.0f, 0.0f, 1.0f}; + bool nonSingular = LU.isNonsingular(); +#if 0 + // UNUSED + float eigMax = 10.0f; + float eigMin = 0.1f; +#endif + if (nonSingular) + { + solveLU3x3(LU, xUnwarped, xWarped); + solveLU3x3(LU, yUnwarped, yWarped); + solveLU3x3(LU, zUnwarped, zWarped); + + // compute the eigenvalues while we have the Jacobian available + Vec3 eigenvalues = Vec3(1.); + computeEigenvalues3x3( &eigenvalues[0], jacobian); + eigMax[indexSmall] = MAX3V(eigenvalues); + eigMin[indexSmall] = MIN3V(eigenvalues); + } + + // make sure to skip one on the beginning and end + int xStart = (xSmall == 0) ? 1 : 0; + int xEnd = (xSmall == _xResSm - 1) ? _amplify - 1 : _amplify; + int yStart = (ySmall == 0) ? 1 : 0; + int yEnd = (ySmall == _yResSm - 1) ? _amplify - 1 : _amplify; + int zStart = (zSmall == 0) ? 1 : 0; + int zEnd = (zSmall == _zResSm - 1) ? _amplify - 1 : _amplify; + + for (int zBig = zStart; zBig < zEnd; zBig++) + for (int yBig = yStart; yBig < yEnd; yBig++) + for (int xBig = xStart; xBig < xEnd; xBig++) + { + const int x = xSmall * _amplify + xBig; + const int y = ySmall * _amplify + yBig; + const int z = zSmall * _amplify + zBig; + + // get unit position for both fine and coarse grid + const Vec3 pos = Vec3(x,y,z); + const Vec3 posSm = pos * invAmp; + + // get grid index for both fine and coarse grid + const int index = x + y *_xResBig + z *_slabSizeBig; + + // get a linearly interpolated velocity and texcoords + // from the coarse grid + Vec3 vel = INTERPOLATE::lerp3dVec( xvel,yvel,zvel, + posSm[0], posSm[1], posSm[2], _xResSm,_yResSm,_zResSm); + Vec3 uvw = INTERPOLATE::lerp3dVec( _tcU,_tcV,_tcW, + posSm[0], posSm[1], posSm[2], _xResSm,_yResSm,_zResSm); + + // multiply the texture coordinate by _resSm so that turbulence + // synthesis begins at the first octave that the coarse grid + // cannot capture + Vec3 texCoord = Vec3(uvw[0] * _resSm[0], + uvw[1] * _resSm[1], + uvw[2] * _resSm[2]); + + // retrieve wavelet energy at highest frequency + float energy = INTERPOLATE::lerp3d( + highFreqEnergy, posSm[0],posSm[1],posSm[2], _xResSm, _yResSm, _zResSm); + + // base amplitude for octave 0 + float coefficient = sqrtf(2.0f * fabs(energy)); + 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 + // to make a difference + const bool addNoise = eigMax[indexSmall] < 2. && + eigMin[indexSmall] > 0.5; + if (addNoise && amplitude > _cullingThreshold) { + // base amplitude for octave 0 + float amplitudeScaled = amplitude; - if (substep < totalSubsteps - 1) - SWAP_POINTERS(_densityBig, _densityBigOld); - } // substep + for (int octave = 0; octave < _octaves; octave++) + { + // multiply the vector noise times the maximum allowed + // noise amplitude at this octave, and add it to the total + vel += WVelocityWithJacobian(texCoord, &xUnwarped[0], &yUnwarped[0], &zUnwarped[0]) * amplitudeScaled; - // wipe the density borders - FLUID_3D::setZeroBorder(_densityBig, _resBig); + // scale coefficient for next octave + amplitudeScaled *= persistence; + texCoord *= 2.0f; + } + } - // reset texture coordinates now in preparation for next timestep - // Shouldn't do this before generating the noise because then the - // eigenvalues stored do not reflect the underlying texture coordinates - resetTextureCoordinates(); + // Store velocity + turbulence in big grid for maccormack step + // + // If you wanted to save memory, you would instead perform a + // semi-Lagrangian backtrace for the current grid cell here. Then + // you could just throw the velocity away. + bigUx[index] = vel[0]; + bigUy[index] = vel[1]; + bigUz[index] = vel[2]; + + // compute the velocity magnitude for substepping later + const float velMag = bigUx[index] * bigUx[index] + + bigUy[index] * bigUy[index] + + bigUz[index] * bigUz[index]; + if (velMag > maxVelMag1) maxVelMag1 = velMag; + + // zero out velocity inside obstacles + float obsCheck = INTERPOLATE::lerp3dToFloat( + obstacles, posSm[0], posSm[1], posSm[2], _xResSm, _yResSm, _zResSm); + if (obsCheck > 0.95) + bigUx[index] = bigUy[index] = bigUz[index] = 0.; + } // xyz - // output files - // string prefix = string("./amplified.preview/density_bigxy_"); - // FLUID_3D::writeImageSliceXY(_densityBig, _resBig, _resBig[2]/2, prefix, _totalStepsBig, 1.0f); - //string df3prefix = string("./df3/density_big_"); - //IMAGE::dumpDF3(_totalStepsBig, df3prefix, _densityBig, _resBig[0],_resBig[1],_resBig[2]); - // string pbrtPrefix = string("./pbrt/density_big_"); - // IMAGE::dumpPBRT(_totalStepsBig, pbrtPrefix, _densityBig, _resBig[0],_resBig[1],_resBig[2]); +#if PARALLEL==1 + maxVelMagThreads[id] = maxVelMag1; +#else + maxVelMagThreads[0] = maxVelMag1; +#endif + } + } // omp + + // compute maximum over threads + float maxVelMag = maxVelMagThreads[0]; +#if PARALLEL==1 + for (int i = 1; i < 8; i++) + if (maxVelMag < maxVelMagThreads[i]) + maxVelMag = maxVelMagThreads[i]; +#endif - _totalStepsBig++; + // prepare density for an advection + SWAP_POINTERS(_densityBig, _densityBigOld); - delete[] _bigUx; - delete[] _bigUy; - delete[] _bigUz; + // based on the maximum velocity present, see if we need to substep, + // but cap the maximum number of substeps to 5 + const int maxSubSteps = 25; + const int maxVel = 5; + 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; - delete[] _tempBig1; - delete[] _tempBig2; + // set boundaries of big velocity grid + FLUID_3D::setZeroX(bigUx, _resBig); + FLUID_3D::setZeroY(bigUy, _resBig); + FLUID_3D::setZeroZ(bigUz, _resBig); + + // do the MacCormack advection, with substepping if necessary + for(int substep = 0; substep < totalSubsteps; substep++) + { + FLUID_3D::advectFieldMacCormack(dtSubdiv, bigUx, bigUy, bigUz, + _densityBigOld, _densityBig, tempBig1, tempBig2, _resBig, NULL); + + if (substep < totalSubsteps - 1) + SWAP_POINTERS(_densityBig, _densityBigOld); + } // substep + + delete[] tempBig1; + delete[] tempBig2; + delete[] bigUx; + delete[] bigUy; + delete[] bigUz; + delete[] _energy; + delete[] highFreqEnergy; + + // wipe the density borders + FLUID_3D::setZeroBorder(_densityBig, _resBig); + + // reset texture coordinates now in preparation for next timestep + // Shouldn't do this before generating the noise because then the + // eigenvalues stored do not reflect the underlying texture coordinates + resetTextureCoordinates(eigMin, eigMax); + + delete[] eigMin; + delete[] eigMax; + + // output files + // string prefix = string("./amplified.preview/density_bigxy_"); + // FLUID_3D::writeImageSliceXY(_densityBig, _resBig, _resBig[2]/2, prefix, _totalStepsBig, 1.0f); + //string df3prefix = string("./df3/density_big_"); + //IMAGE::dumpDF3(_totalStepsBig, df3prefix, _densityBig, _resBig[0],_resBig[1],_resBig[2]); + // string pbrtPrefix = string("./pbrt/density_big_"); + // IMAGE::dumpPBRT(_totalStepsBig, pbrtPrefix, _densityBig, _resBig[0],_resBig[1],_resBig[2]); + + _totalStepsBig++; } diff --git a/intern/smoke/intern/WTURBULENCE.h b/intern/smoke/intern/WTURBULENCE.h index c21e002ad48..0aa978e9e52 100644 --- a/intern/smoke/intern/WTURBULENCE.h +++ b/intern/smoke/intern/WTURBULENCE.h @@ -49,10 +49,10 @@ class WTURBULENCE void stepTurbulenceFull(float dt, float* xvel, float* yvel, float* zvel, unsigned char *obstacles); // texcoord functions - void advectTextureCoordinates (float dtOrg, float* xvel, float* yvel, float* zvel, float *_tempBig1, float *_tempBig2); - void resetTextureCoordinates(); + void advectTextureCoordinates(float dtOrg, float* xvel, float* yvel, float* zvel, float *tempBig1, float *tempBig2); + void resetTextureCoordinates(float *_eigMin, float *_eigMax); - void computeEnergy(float* xvel, float* yvel, float* zvel, unsigned char *obstacles); + void computeEnergy(float *energy, float* xvel, float* yvel, float* zvel, unsigned char *obstacles); // evaluate wavelet noise function Vec3 WVelocity(Vec3 p); @@ -63,8 +63,6 @@ class WTURBULENCE inline float* getArrayTcU() { return _tcU; } inline float* getArrayTcV() { return _tcV; } inline float* getArrayTcW() { return _tcW; } - inline float* getArrayEigMin() { return _eigMin; } - inline float* getArrayEigMax() { return _eigMax; } inline Vec3Int getResSm() { return _resSm; } // small resolution inline Vec3Int getResBig() { return _resBig; } @@ -81,15 +79,15 @@ class WTURBULENCE // noise settings float _cullingThreshold; - float _noiseStrength; - float _noiseSizeScale; - bool _uvwAdvection; - bool _uvwReset; - float _noiseTimeanimSpeed; - int _dumpInterval; - int _noiseControlType; + // float _noiseStrength; + // float _noiseSizeScale; + // bool _uvwAdvection; + // bool _uvwReset; + // float _noiseTimeanimSpeed; + // int _dumpInterval; + // nt _noiseControlType; // debug, scale density for projections output images - float _outputScale; + // float _outputScale; // noise resolution int _xResBig; @@ -117,24 +115,15 @@ class WTURBULENCE float* _tcW; float* _tcTemp; - float* _eigMin; // no save -dg - float* _eigMax; // no save -dg - - // wavelet decomposition of velocity energies - float* _energy; // no save -dg - // noise data float* _noiseTile; //float* _noiseTileExt; // step counter int _totalStepsBig; - - // highest frequency component of wavelet decomposition - float* _highFreqEnergy; // no save -dg - void computeEigenvalues(); - void decomposeEnergy(); + void computeEigenvalues(float *_eigMin, float *_eigMax); + void decomposeEnergy(float *energy, float *_highFreqEnergy); }; #endif // WTURBULENCE_H diff --git a/intern/smoke/intern/smoke_API.cpp b/intern/smoke/intern/smoke_API.cpp index 5a016b51bbe..058831dddbb 100644 --- a/intern/smoke/intern/smoke_API.cpp +++ b/intern/smoke/intern/smoke_API.cpp @@ -137,7 +137,7 @@ extern "C" void smoke_dissolve(FLUID_3D *fluid, int speed, int log) } } -extern "C" void smoke_turbulence_dissolve(WTURBULENCE *wt, int speed, int log) +extern "C" void smoke_dissolve_wavelet(WTURBULENCE *wt, int speed, int log) { float *density = wt->getDensityBig(); Vec3Int r = wt->getResBig(); @@ -172,7 +172,7 @@ extern "C" void smoke_turbulence_dissolve(WTURBULENCE *wt, int speed, int log) } } -extern "C" void smoke_turbulence_initBlenderRNA(WTURBULENCE *wt, float *strength) +extern "C" void smoke_initWaveletBlenderRNA(WTURBULENCE *wt, float *strength) { wt->initBlenderRNA(strength); } @@ -241,13 +241,14 @@ extern "C" float *smoke_turbulence_get_density(WTURBULENCE *wt) return wt ? wt->getDensityBig() : NULL; } -extern "C" void smoke_turbulence_get_res(WTURBULENCE *wt, unsigned int *res) +extern "C" void smoke_turbulence_get_res(WTURBULENCE *wt, int *res) { if(wt) { - res[0] = wt->_resBig[0]; - res[1] = wt->_resBig[1]; - res[2] = wt->_resBig[2]; + Vec3Int r = wt->getResBig(); + res[0] = r[0]; + res[1] = r[1]; + res[2] = r[2]; } } diff --git a/release/ui/buttons_physics_smoke.py b/release/ui/buttons_physics_smoke.py index c87f71bff42..83c1ffc2e9a 100644 --- a/release/ui/buttons_physics_smoke.py +++ b/release/ui/buttons_physics_smoke.py @@ -3,11 +3,6 @@ import bpy from buttons_particle import point_cache_ui -def smoke_panel_enabled_low(smd): - if smd.smoke_type == 'TYPE_DOMAIN': - return smd.domain.point_cache.baked==False - return True - class PhysicButtonsPanel(bpy.types.Panel): __space_type__ = 'PROPERTIES' __region_type__ = 'WINDOW' @@ -45,8 +40,6 @@ class PHYSICS_PT_smoke(PhysicButtonsPanel): split.itemL() if md: - - # layout.enabled = smoke_panel_enabled(md) layout.itemR(md, "smoke_type", expand=True) if md.smoke_type == 'TYPE_DOMAIN': @@ -66,7 +59,7 @@ class PHYSICS_PT_smoke(PhysicButtonsPanel): col.itemR(domain, "dissolve_smoke", text="Dissolve") sub = col.column() sub.active = domain.dissolve_smoke - sub.itemR(domain, "dissolve_speed", text="Speed") + sub.itemR(domain, "dissolve_speed", text="Time") sub.itemR(domain, "dissolve_smoke_log", text="Slow") elif md.smoke_type == 'TYPE_FLOW': @@ -90,14 +83,17 @@ class PHYSICS_PT_smoke(PhysicButtonsPanel): #elif md.smoke_type == 'TYPE_COLL': # layout.itemS() - + class PHYSICS_PT_smoke_groups(PhysicButtonsPanel): __label__ = "Smoke Groups" __default_closed__ = True def poll(self, context): md = context.smoke - return md and (md.smoke_type == 'TYPE_DOMAIN') + if md: + return (md.smoke_type == 'TYPE_DOMAIN') + + return False def draw(self, context): layout = self.layout @@ -116,7 +112,7 @@ class PHYSICS_PT_smoke_groups(PhysicButtonsPanel): col = split.column() col.itemL(text="Collision Group:") col.itemR(group, "coll_group", text="") - + class PHYSICS_PT_smoke_cache(PhysicButtonsPanel): __label__ = "Smoke Cache" __default_closed__ = True @@ -129,7 +125,7 @@ class PHYSICS_PT_smoke_cache(PhysicButtonsPanel): layout = self.layout md = context.smoke.domain_settings - cache = md.point_cache + cache = md.point_cache_low point_cache_ui(self, cache, cache.baked==False, 0, 1) @@ -152,7 +148,7 @@ class PHYSICS_PT_smoke_highres(PhysicButtonsPanel): md = context.smoke.domain_settings split = layout.split() - + col = split.column() col.itemL(text="Resolution:") col.itemR(md, "amplify", text="Divisions") @@ -161,66 +157,26 @@ class PHYSICS_PT_smoke_highres(PhysicButtonsPanel): col.itemL(text="Noise Method:") col.row().itemR(md, "noise_type", text="") col.itemR(md, "strength") - col.itemR(md, "show_highres") + sub.itemR(md, "viewhighres") class PHYSICS_PT_smoke_cache_highres(PhysicButtonsPanel): - __label__ = "Smoke Cache" + __label__ = "Smoke High Resolution Cache" __default_closed__ = True def poll(self, context): - return (context.smoke) + md = context.smoke + return md and (md.smoke_type == 'TYPE_DOMAIN') and md.domain_settings.highres def draw(self, context): layout = self.layout - md = context.smoke - - cache = md.point_cache - - layout.set_context_pointer("PointCache", cache) - - row = layout.row() - row.template_list(cache, "point_cache_list", cache, "active_point_cache_index") - col = row.column(align=True) - col.itemO("ptcache.add_new", icon='ICON_ZOOMIN', text="") - col.itemO("ptcache.remove", icon='ICON_ZOOMOUT', text="") - - row = layout.row() - row.itemR(cache, "name") - - row = layout.row() - row.itemR(cache, "start_frame") - row.itemR(cache, "end_frame") - - row = layout.row() - - if cache.baked == True: - row.itemO("ptcache.free_bake", text="Free Bake") - else: - row.item_booleanO("ptcache.bake", "bake", True, text="Bake") - - subrow = row.row() - subrow.enabled = cache.frames_skipped or cache.outdated - subrow.itemO("ptcache.bake", "bake", False, text="Calculate to Current Frame") - - row = layout.row() - #row.enabled = smoke_panel_enabled(psys) - row.itemO("ptcache.bake_from_cache", text="Current Cache to Bake") - - row = layout.row() - #row.enabled = smoke_panel_enabled(psys) - - layout.itemL(text=cache.info) - - layout.itemS() + md = context.smoke.domain_settings + cache = md.point_cache_high - row = layout.row() - row.itemO("ptcache.bake_all", "bake", True, text="Bake All Dynamics") - row.itemO("ptcache.free_bake_all", text="Free All Bakes") - layout.itemO("ptcache.bake_all", "bake", False, text="Update All Dynamics to current frame") - + point_cache_ui(self, cache, cache.baked==False, 0, 1) + bpy.types.register(PHYSICS_PT_smoke) bpy.types.register(PHYSICS_PT_smoke_cache) +bpy.types.register(PHYSICS_PT_smoke_highres) bpy.types.register(PHYSICS_PT_smoke_groups) -#bpy.types.register(PHYSICS_PT_smoke_highres) -#bpy.types.register(PHYSICS_PT_smoke_cache_highres) +bpy.types.register(PHYSICS_PT_smoke_cache_highres) diff --git a/source/blender/blenkernel/BKE_pointcache.h b/source/blender/blenkernel/BKE_pointcache.h index b7ab07b0f91..5ae10d736fd 100644 --- a/source/blender/blenkernel/BKE_pointcache.h +++ b/source/blender/blenkernel/BKE_pointcache.h @@ -232,6 +232,7 @@ void BKE_ptcache_id_from_softbody(PTCacheID *pid, struct Object *ob, struct Soft void BKE_ptcache_id_from_particles(PTCacheID *pid, struct Object *ob, struct ParticleSystem *psys); void BKE_ptcache_id_from_cloth(PTCacheID *pid, struct Object *ob, struct ClothModifierData *clmd); void BKE_ptcache_id_from_smoke(PTCacheID *pid, struct Object *ob, struct SmokeModifierData *smd); +void BKE_ptcache_id_from_smoke_turbulence(PTCacheID *pid, struct Object *ob, struct SmokeModifierData *smd); void BKE_ptcache_ids_from_object(struct ListBase *lb, struct Object *ob); diff --git a/source/blender/blenkernel/BKE_smoke.h b/source/blender/blenkernel/BKE_smoke.h index fddcf0fea83..0f8e9c5edf5 100644 --- a/source/blender/blenkernel/BKE_smoke.h +++ b/source/blender/blenkernel/BKE_smoke.h @@ -32,13 +32,15 @@ #ifndef BKE_SMOKE_H_ #define BKE_SMOKE_H_ -typedef int (*bresenham_callback) (float *input, int res[3], int *pixel, float *tRay); +typedef float (*bresenham_callback) (float *result, float *input, int res[3], int *pixel, float *tRay, float correct); void smokeModifier_do(struct SmokeModifierData *smd, struct Scene *scene, struct Object *ob, struct DerivedMesh *dm, int useRenderParams, int isFinalCalc); void smokeModifier_free (struct SmokeModifierData *smd); void smokeModifier_reset(struct SmokeModifierData *smd); +void smokeModifier_reset_turbulence(struct SmokeModifierData *smd); void smokeModifier_createType(struct SmokeModifierData *smd); +long long smoke_get_mem_req(int xres, int yres, int zres, int amplify); #endif /* BKE_SMOKE_H_ */ diff --git a/source/blender/blenkernel/intern/pointcache.c b/source/blender/blenkernel/intern/pointcache.c index 912acc4786f..1ab1daf1782 100644 --- a/source/blender/blenkernel/intern/pointcache.c +++ b/source/blender/blenkernel/intern/pointcache.c @@ -471,6 +471,19 @@ static int ptcache_totpoint_smoke(void *smoke_v) return 0; } +/* Smoke functions */ +static int ptcache_totpoint_smoke_turbulence(void *smoke_v) +{ + SmokeModifierData *smd= (SmokeModifierData *)smoke_v; + SmokeDomainSettings *sds = smd->domain; + + if(sds->wt) { + return sds->res_wt[0]*sds->res_wt[1]*sds->res_wt[2]; + } + else + return 0; +} + // forward decleration static int ptcache_file_write(PTCacheFile *pf, void *f, size_t tot, int size); @@ -521,7 +534,7 @@ static int ptcache_compress_write(PTCacheFile *pf, unsigned char *in, unsigned i } static int ptcache_write_smoke(PTCacheFile *pf, void *smoke_v) -{ +{ SmokeModifierData *smd= (SmokeModifierData *)smoke_v; SmokeDomainSettings *sds = smd->domain; @@ -535,7 +548,7 @@ static int ptcache_write_smoke(PTCacheFile *pf, void *smoke_v) smoke_export(sds->fluid, &dt, &dx, &dens, &densold, &heat, &heatold, &vx, &vy, &vz, &vxold, &vyold, &vzold, &obstacles); - ptcache_compress_write(pf, (unsigned char *)sds->view3d, in_len*4, out, mode); + ptcache_compress_write(pf, (unsigned char *)sds->shadow, in_len, out, mode); ptcache_compress_write(pf, (unsigned char *)dens, in_len, out, mode); ptcache_compress_write(pf, (unsigned char *)densold, in_len, out, mode); ptcache_compress_write(pf, (unsigned char *)heat, in_len, out, mode); @@ -554,36 +567,36 @@ static int ptcache_write_smoke(PTCacheFile *pf, void *smoke_v) return 1; } - return 0; } -/* static int ptcache_write_smoke_turbulence(PTCacheFile *pf, void *smoke_v) { SmokeModifierData *smd= (SmokeModifierData *)smoke_v; SmokeDomainSettings *sds = smd->domain; if(sds->wt) { - unsigned int res_big[3]; - size_t res = sds->res[0]*sds->res[1]*sds->res[2]; + unsigned int res_big_array[3]; + unsigned int res_big; + unsigned int res = sds->res[0]*sds->res[1]*sds->res[2]; float *dens, *densold, *tcu, *tcv, *tcw; unsigned int in_len = sizeof(float)*(unsigned int)res; - unsigned int in_len_big = sizeof(float) * (unsigned int)res_big; + unsigned int in_len_big; unsigned char *out; int mode; - smoke_turbulence_get_res(sds->wt, res_big); - mode = res_big[0]*res_big[1]*res_big[2] >= 1000000 ? 2 : 1; + smoke_turbulence_get_res(sds->wt, res_big_array); + res_big = res_big_array[0]*res_big_array[1]*res_big_array[2]; + mode = res_big >= 1000000 ? 2 : 1; + in_len_big = sizeof(float) * (unsigned int)res_big; smoke_turbulence_export(sds->wt, &dens, &densold, &tcu, &tcv, &tcw); out = (unsigned char *)MEM_callocN(LZO_OUT_LEN(in_len_big), "pointcache_lzo_buffer"); - ptcache_compress_write(pf, (unsigned char *)dens, in_len_big, out, mode); ptcache_compress_write(pf, (unsigned char *)densold, in_len_big, out, mode); - MEM_freeN(out); + out = (unsigned char *)MEM_callocN(LZO_OUT_LEN(in_len), "pointcache_lzo_buffer"); ptcache_compress_write(pf, (unsigned char *)tcu, in_len, out, mode); ptcache_compress_write(pf, (unsigned char *)tcv, in_len, out, mode); @@ -594,7 +607,6 @@ static int ptcache_write_smoke_turbulence(PTCacheFile *pf, void *smoke_v) } return 0; } -*/ // forward decleration static int ptcache_file_read(PTCacheFile *pf, void *f, size_t tot, int size); @@ -649,7 +661,7 @@ static void ptcache_read_smoke(PTCacheFile *pf, void *smoke_v) smoke_export(sds->fluid, &dt, &dx, &dens, &densold, &heat, &heatold, &vx, &vy, &vz, &vxold, &vyold, &vzold, &obstacles); - ptcache_compress_read(pf, (unsigned char *)sds->view3d, out_len*4); + ptcache_compress_read(pf, (unsigned char *)sds->shadow, out_len); ptcache_compress_read(pf, (unsigned char*)dens, out_len); ptcache_compress_read(pf, (unsigned char*)densold, out_len); ptcache_compress_read(pf, (unsigned char*)heat, out_len); @@ -666,26 +678,32 @@ static void ptcache_read_smoke(PTCacheFile *pf, void *smoke_v) } } -/* static void ptcache_read_smoke_turbulence(PTCacheFile *pf, void *smoke_v) { SmokeModifierData *smd= (SmokeModifierData *)smoke_v; SmokeDomainSettings *sds = smd->domain; if(sds->fluid) { - unsigned int res[3]; + unsigned int res = sds->res[0]*sds->res[1]*sds->res[2]; + unsigned int res_big, res_big_array[3]; float *dens, *densold, *tcu, *tcv, *tcw; unsigned int out_len = sizeof(float)*(unsigned int)res; + unsigned int out_len_big; - smoke_turbulence_get_res(sds->wt, res); + smoke_turbulence_get_res(sds->wt, res_big_array); + res_big = res_big_array[0]*res_big_array[1]*res_big_array[2]; + out_len_big = sizeof(float) * (unsigned int)res_big; smoke_turbulence_export(sds->wt, &dens, &densold, &tcu, &tcv, &tcw); - ptcache_compress_read(pf, (unsigned char*)dens, out_len); - + ptcache_compress_read(pf, (unsigned char*)dens, out_len_big); + ptcache_compress_read(pf, (unsigned char*)densold, out_len_big); + + ptcache_compress_read(pf, (unsigned char*)tcu, out_len); + ptcache_compress_read(pf, (unsigned char*)tcv, out_len); + ptcache_compress_read(pf, (unsigned char*)tcw, out_len); } } -*/ void BKE_ptcache_id_from_smoke(PTCacheID *pid, struct Object *ob, struct SmokeModifierData *smd) { @@ -716,7 +734,7 @@ void BKE_ptcache_id_from_smoke(PTCacheID *pid, struct Object *ob, struct SmokeMo pid->write_header= ptcache_write_basic_header; pid->read_header= ptcache_read_basic_header; - pid->data_types= (1<<BPHYS_DATA_LOCATION); // bogus values tot make pointcache happy + pid->data_types= (1<<BPHYS_DATA_LOCATION); // bogus values to make pointcache happy pid->info_types= 0; } @@ -736,13 +754,13 @@ void BKE_ptcache_id_from_smoke_turbulence(PTCacheID *pid, struct Object *ob, str pid->cache_ptr= &sds->point_cache[1]; pid->ptcaches= &sds->ptcaches[1]; - pid->totpoint= pid->totwrite= ptcache_totpoint_smoke; + pid->totpoint= pid->totwrite= ptcache_totpoint_smoke_turbulence; pid->write_elem= NULL; pid->read_elem= NULL; - pid->read_stream = ptcache_read_smoke; - pid->write_stream = ptcache_write_smoke; + pid->read_stream = ptcache_read_smoke_turbulence; + pid->write_stream = ptcache_write_smoke_turbulence; pid->interpolate_elem= NULL; @@ -820,6 +838,10 @@ void BKE_ptcache_ids_from_object(ListBase *lb, Object *ob) pid= MEM_callocN(sizeof(PTCacheID), "PTCacheID"); BKE_ptcache_id_from_smoke(pid, ob, (SmokeModifierData*)md); BLI_addtail(lb, pid); + + pid= MEM_callocN(sizeof(PTCacheID), "PTCacheID"); + BKE_ptcache_id_from_smoke_turbulence(pid, ob, (SmokeModifierData*)md); + BLI_addtail(lb, pid); } } } @@ -1824,6 +1846,8 @@ int BKE_ptcache_id_reset(Scene *scene, PTCacheID *pid, int mode) psys_reset(pid->calldata, PSYS_RESET_DEPSGRAPH); else if(pid->type == PTCACHE_TYPE_SMOKE_DOMAIN) smokeModifier_reset(pid->calldata); + else if(pid->type == PTCACHE_TYPE_SMOKE_HIGHRES) + smokeModifier_reset(pid->calldata); } if(clear) BKE_ptcache_id_clear(pid, PTCACHE_CLEAR_ALL, 0); @@ -1878,6 +1902,9 @@ int BKE_ptcache_object_reset(Scene *scene, Object *ob, int mode) { BKE_ptcache_id_from_smoke(&pid, ob, (SmokeModifierData*)md); reset |= BKE_ptcache_id_reset(scene, &pid, mode); + + BKE_ptcache_id_from_smoke_turbulence(&pid, ob, (SmokeModifierData*)md); + reset |= BKE_ptcache_id_reset(scene, &pid, mode); } } } diff --git a/source/blender/blenkernel/intern/smoke.c b/source/blender/blenkernel/intern/smoke.c index 223d48012df..026ccf4fc7d 100644 --- a/source/blender/blenkernel/intern/smoke.c +++ b/source/blender/blenkernel/intern/smoke.c @@ -92,10 +92,10 @@ static void tend ( void ) { QueryPerformanceCounter ( &liCurrentTime ); } -//static double tval() -//{ -// return ((double)( (liCurrentTime.QuadPart - liStartTime.QuadPart)* (double)1000.0/(double)liFrequency.QuadPart )); -//} +static double tval() +{ + return ((double)( (liCurrentTime.QuadPart - liStartTime.QuadPart)* (double)1000.0/(double)liFrequency.QuadPart )); +} #else #include <sys/time.h> static struct timeval _tstart, _tend; @@ -125,7 +125,6 @@ struct SmokeModifierData; // forward declerations static void get_cell(float *p0, int res[3], float dx, float *pos, int *cell, int correct); void calcTriangleDivs(Object *ob, MVert *verts, int numverts, MFace *tris, int numfaces, int numtris, int **tridivs, float cell_len); -void smoke_prepare_View(SmokeModifierData *smd, float framenr, float *light, int have_light); #define TRI_UVOFFSET (1./4.) @@ -167,6 +166,7 @@ int smokeModifier_init (SmokeModifierData *smd, Object *ob, Scene *scene, Derive // calc other res with max_res provided VECSUB(size, max, min); + // prevent crash when initializing a plane as domain if((size[0] < FLT_EPSILON) || (size[1] < FLT_EPSILON) || (size[2] < FLT_EPSILON)) return 0; @@ -215,22 +215,28 @@ int smokeModifier_init (SmokeModifierData *smd, Object *ob, Scene *scene, Derive // dt max is 0.1 smd->domain->fluid = smoke_init(smd->domain->res, smd->domain->p0, 0.1); smd->time = scene->r.cfra; - smd->domain->firstframe = smd->time; - if(!smd->domain->wt && (smd->domain->flags & MOD_SMOKE_HIGHRES)) + if(smd->domain->flags & MOD_SMOKE_HIGHRES) { - smd->domain->wt = smoke_turbulence_init(smd->domain->res, smd->domain->amplify + 1, smd->domain->noise); - smoke_turbulence_initBlenderRNA(smd->domain->wt, &smd->domain->strength); + smd->domain->wt = smoke_turbulence_init(smd->domain->res, smd->domain->amplify + 1, smd->domain->noise); + smd->domain->res_wt[0] = smd->domain->res[0] * (smd->domain->amplify + 1); + smd->domain->res_wt[1] = smd->domain->res[1] * (smd->domain->amplify + 1); + smd->domain->res_wt[2] = smd->domain->res[2] * (smd->domain->amplify + 1); + smd->domain->dx_wt = smd->domain->dx / (smd->domain->amplify + 1); + printf("smd->domain->amplify: %d\n", smd->domain->amplify); + printf("(smd->domain->flags & MOD_SMOKE_HIGHRES)\n"); } - if(!smd->domain->view3d) - { - // TVox is for transparency - smd->domain->view3d = MEM_callocN(sizeof(float)*smd->domain->res[0]*smd->domain->res[1]*smd->domain->res[2]*4, "Smoke_tVox"); - } + if(!smd->domain->shadow) + smd->domain->shadow = MEM_callocN(sizeof(float) * smd->domain->res[0] * smd->domain->res[1] * smd->domain->res[2], "SmokeDomainShadow"); smoke_initBlenderRNA(smd->domain->fluid, &(smd->domain->alpha), &(smd->domain->beta)); + if(smd->domain->wt) + { + smoke_initWaveletBlenderRNA(smd->domain->wt, &(smd->domain->strength)); + printf("smoke_initWaveletBlenderRNA\n"); + } return 1; } else if((smd->type & MOD_SMOKE_TYPE_FLOW) && smd->flow) @@ -270,13 +276,11 @@ int smokeModifier_init (SmokeModifierData *smd, Object *ob, Scene *scene, Derive SmokeCollSettings *scs = smd->coll; MVert *mvert = dm->getVertArray(dm); MFace *mface = dm->getFaceArray(dm); - size_t i = 0; - int divs = 0; + int i = 0, divs = 0; int *tridivs = NULL; float cell_len = 1.0 / 50.0; // for res = 50 - size_t newdivs = 0; - //size_t max_points = 0; - size_t quads = 0, facecounter = 0; + int newdivs = 0; + int quads = 0, facecounter = 0; // copy obmat Mat4CpyMat4(scs->mat, ob->obmat); @@ -314,7 +318,7 @@ int smokeModifier_init (SmokeModifierData *smd, Object *ob, Scene *scene, Derive int again = 0; do { - size_t j, k; + int j, k; int divs1 = tridivs[3 * facecounter + 0]; int divs2 = tridivs[3 * facecounter + 1]; //int divs3 = tridivs[3 * facecounter + 2]; @@ -521,10 +525,10 @@ void smokeModifier_freeDomain(SmokeModifierData *smd) { if(smd->domain) { - // free visualisation buffers - if(smd->domain->view3d) - MEM_freeN(smd->domain->view3d); - + if(smd->domain->shadow) + MEM_freeN(smd->domain->shadow); + smd->domain->shadow = NULL; + if(smd->domain->fluid) smoke_free(smd->domain->fluid); @@ -583,41 +587,46 @@ void smokeModifier_freeCollision(SmokeModifierData *smd) } } +void smokeModifier_reset_turbulence(struct SmokeModifierData *smd) +{ + if(smd && smd->domain && smd->domain->wt) + { + smoke_turbulence_free(smd->domain->wt); + smd->domain->wt = NULL; + } + + smd->domain->point_cache[1]->flag &= ~PTCACHE_SIMULATION_VALID; + smd->domain->point_cache[1]->flag |= PTCACHE_OUTDATED; + smd->domain->point_cache[1]->simframe= 0; + smd->domain->point_cache[1]->last_exact= 0; +} + void smokeModifier_reset(struct SmokeModifierData *smd) { if(smd) { if(smd->domain) { - if(smd->domain->view3d) - MEM_freeN(smd->domain->view3d); - smd->domain->view3d = NULL; - - smd->domain->tex = NULL; + if(smd->domain->shadow) + MEM_freeN(smd->domain->shadow); + smd->domain->shadow = NULL; if(smd->domain->fluid) { smoke_free(smd->domain->fluid); smd->domain->fluid = NULL; } - - if(smd->domain->wt) - { - smoke_turbulence_free(smd->domain->wt); - smd->domain->wt = NULL; - } - + smd->domain->point_cache[0]->flag &= ~PTCACHE_SIMULATION_VALID; smd->domain->point_cache[0]->flag |= PTCACHE_OUTDATED; smd->domain->point_cache[0]->simframe= 0; smd->domain->point_cache[0]->last_exact= 0; - smd->domain->point_cache[1]->flag &= ~PTCACHE_SIMULATION_VALID; - smd->domain->point_cache[1]->flag |= PTCACHE_OUTDATED; - smd->domain->point_cache[1]->simframe= 0; - smd->domain->point_cache[1]->last_exact= 0; + smokeModifier_reset_turbulence(smd); - // printf("reset_domain\n"); + smd->time = -1; + + // printf("reset domain end\n"); } else if(smd->flow) { @@ -685,22 +694,21 @@ void smokeModifier_createType(struct SmokeModifierData *smd) /* set some standard values */ smd->domain->fluid = NULL; + smd->domain->wt = NULL; smd->domain->eff_group = NULL; smd->domain->fluid_group = NULL; smd->domain->coll_group = NULL; smd->domain->maxres = 32; + smd->domain->amplify = 1; + smd->domain->omega = 1.0; smd->domain->alpha = -0.001; smd->domain->beta = 0.1; smd->domain->flags = MOD_SMOKE_DISSOLVE_LOG; - smd->domain->diss_speed = 5; - smd->domain->strength = 2.0f; - smd->domain->amplify = 1; + smd->domain->strength = 2.0; smd->domain->noise = MOD_SMOKE_NOISEWAVE; - smd->domain->wt = NULL; - + smd->domain->diss_speed = 5; // init 3dview buffer - smd->domain->view3d = NULL; - smd->domain->tex = NULL; + smd->domain->viewsettings = 0; } else if(smd->type & MOD_SMOKE_TYPE_FLOW) { @@ -734,9 +742,314 @@ void smokeModifier_createType(struct SmokeModifierData *smd) } } -// forward declaration -void smoke_simulate_domain(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedMesh *dm); +// forward decleration +void smoke_calc_transparency(float *result, float *input, float *p0, float *p1, int res[3], float dx, float *light, bresenham_callback cb, float correct); +static float calc_voxel_transp(float *result, float *input, int res[3], int *pixel, float *tRay, float correct); +static int get_lamp(Scene *scene, float *light) +{ + Base *base_tmp = NULL; + for(base_tmp = scene->base.first; base_tmp; base_tmp= base_tmp->next) + { + if(base_tmp->object->type == OB_LAMP) + { + Lamp *la = (Lamp *)base_tmp->object->data; + + if(la->type == LA_LOCAL) + { + VECCOPY(light, base_tmp->object->obmat[3]); + return 1; + } + } + } + return 0; +} + +static void smoke_calc_domain(Scene *scene, Object *ob, SmokeModifierData *smd) +{ + SmokeDomainSettings *sds = smd->domain; + GroupObject *go = NULL; + Base *base = NULL; + + // do flows and fluids + if(1) + { + Object *otherobj = NULL; + ModifierData *md = NULL; + if(sds->fluid_group) // we use groups since we have 2 domains + go = sds->fluid_group->gobject.first; + else + base = scene->base.first; + while(base || go) + { + otherobj = NULL; + if(sds->fluid_group) + { + if(go->ob) + otherobj = go->ob; + } + else + otherobj = base->object; + if(!otherobj) + { + if(sds->fluid_group) + go = go->next; + else + base= base->next; + + continue; + } + + md = modifiers_findByType(otherobj, eModifierType_Smoke); + + // check for active smoke modifier + if(md && md->mode & (eModifierMode_Realtime | eModifierMode_Render)) + { + SmokeModifierData *smd2 = (SmokeModifierData *)md; + + // check for initialized smoke object + if((smd2->type & MOD_SMOKE_TYPE_FLOW) && smd2->flow) + { + // we got nice flow object + SmokeFlowSettings *sfs = smd2->flow; + + if(sfs->psys && sfs->psys->part && sfs->psys->part->type==PART_EMITTER) // is particle system selected + { + ParticleSystem *psys = sfs->psys; + ParticleSettings *part=psys->part; + ParticleData *pa = NULL; + int p = 0; + float *density = smoke_get_density(sds->fluid); + float *bigdensity = smoke_turbulence_get_density(sds->wt); + float *heat = smoke_get_heat(sds->fluid); + float *velocity_x = smoke_get_velocity_x(sds->fluid); + float *velocity_y = smoke_get_velocity_y(sds->fluid); + float *velocity_z = smoke_get_velocity_z(sds->fluid); + unsigned char *obstacle = smoke_get_obstacle(sds->fluid); + int bigres[3]; + + // mostly copied from particle code + for(p=0, pa=psys->particles; p<psys->totpart; p++, pa++) + { + int cell[3]; + size_t i = 0; + size_t index = 0; + int badcell = 0; + if(pa->alive == PARS_KILLED) continue; + else if(pa->alive == PARS_UNBORN && (part->flag & PART_UNBORN)==0) continue; + else if(pa->alive == PARS_DEAD && (part->flag & PART_DIED)==0) continue; + else if(pa->flag & (PARS_UNEXIST+PARS_NO_DISP)) continue; + // VECCOPY(pos, pa->state.co); + // Mat4MulVecfl (ob->imat, pos); + // 1. get corresponding cell + get_cell(smd->domain->p0, smd->domain->res, smd->domain->dx, pa->state.co, cell, 0); + // check if cell is valid (in the domain boundary) + for(i = 0; i < 3; i++) + { + if((cell[i] > sds->res[i] - 1) || (cell[i] < 0)) + { + badcell = 1; + break; + } + } + if(badcell) + continue; + // 2. set cell values (heat, density and velocity) + index = smoke_get_index(cell[0], sds->res[0], cell[1], sds->res[1], cell[2]); + if(!(sfs->type & MOD_SMOKE_FLOW_TYPE_OUTFLOW) && !(obstacle[index] & 2)) // this is inflow + { + // heat[index] += sfs->temp * 0.1; + // density[index] += sfs->density * 0.1; + heat[index] = sfs->temp; + density[index] = sfs->density; + + /* + velocity_x[index] = pa->state.vel[0]; + velocity_y[index] = pa->state.vel[1]; + velocity_z[index] = pa->state.vel[2]; + */ + + // obstacle[index] |= 2; + // we need different handling for the high-res feature + if(bigdensity) + { + // init all surrounding cells according to amplification, too + int i, j, k; + + smoke_turbulence_get_res(smd->domain->wt, bigres); + + for(i = 0; i < smd->domain->amplify + 1; i++) + for(j = 0; j < smd->domain->amplify + 1; j++) + for(k = 0; k < smd->domain->amplify + 1; k++) + { + index = smoke_get_index((smd->domain->amplify + 1)* cell[0] + i, bigres[0], (smd->domain->amplify + 1)* cell[1] + j, bigres[1], (smd->domain->amplify + 1)* cell[2] + k); + bigdensity[index] = sfs->density; + } + } + } + else if(sfs->type & MOD_SMOKE_FLOW_TYPE_OUTFLOW) // outflow + { + heat[index] = 0.f; + density[index] = 0.f; + velocity_x[index] = 0.f; + velocity_y[index] = 0.f; + velocity_z[index] = 0.f; + // we need different handling for the high-res feature + if(bigdensity) + { + // init all surrounding cells according to amplification, too + int i, j, k; + smoke_turbulence_get_res(smd->domain->wt, bigres); + + for(i = 0; i < smd->domain->amplify + 1; i++) + for(j = 0; j < smd->domain->amplify + 1; j++) + for(k = 0; k < smd->domain->amplify + 1; k++) + { + index = smoke_get_index((smd->domain->amplify + 1)* cell[0] + i, bigres[0], (smd->domain->amplify + 1)* cell[1] + j, bigres[1], (smd->domain->amplify + 1)* cell[2] + k); + bigdensity[index] = 0.f; + } + } + } // particles loop + } + } + else + { + /* + for() + { + // no psys + BVHTreeNearest nearest; + nearest.index = -1; + nearest.dist = FLT_MAX; + + BLI_bvhtree_find_nearest(sfs->bvh->tree, pco, &nearest, sfs->bvh->nearest_callback, sfs->bvh); + }*/ + } + } + } + if(sds->fluid_group) + go = go->next; + else + base= base->next; + } + } + + // do effectors + /* + if(sds->eff_group) + { + for(go = sds->eff_group->gobject.first; go; go = go->next) + { + if(go->ob) + { + if(ob->pd) + { + + } + } + } + } + */ + + // do collisions + if(1) + { + Object *otherobj = NULL; + ModifierData *md = NULL; + + if(sds->coll_group) // we use groups since we have 2 domains + go = sds->coll_group->gobject.first; + else + base = scene->base.first; + + while(base || go) + { + otherobj = NULL; + if(sds->coll_group) + { + if(go->ob) + otherobj = go->ob; + } + else + otherobj = base->object; + if(!otherobj) + { + if(sds->coll_group) + go = go->next; + else + base= base->next; + continue; + } + md = modifiers_findByType(otherobj, eModifierType_Smoke); + + // check for active smoke modifier + if(md && md->mode & (eModifierMode_Realtime | eModifierMode_Render)) + { + SmokeModifierData *smd2 = (SmokeModifierData *)md; + + if((smd2->type & MOD_SMOKE_TYPE_COLL) && smd2->coll) + { + // we got nice collision object + SmokeCollSettings *scs = smd2->coll; + size_t i, j; + unsigned char *obstacles = smoke_get_obstacle(smd->domain->fluid); + + for(i = 0; i < scs->numpoints; i++) + { + int badcell = 0; + size_t index = 0; + int cell[3]; + + // 1. get corresponding cell + get_cell(smd->domain->p0, smd->domain->res, smd->domain->dx, &scs->points[3 * i], cell, 0); + + // check if cell is valid (in the domain boundary) + for(j = 0; j < 3; j++) + if((cell[j] > sds->res[j] - 1) || (cell[j] < 0)) + { + badcell = 1; + break; + } + + if(badcell) + continue; + // 2. set cell values (heat, density and velocity) + index = smoke_get_index(cell[0], sds->res[0], cell[1], sds->res[1], cell[2]); + + // printf("cell[0]: %d, cell[1]: %d, cell[2]: %d\n", cell[0], cell[1], cell[2]); + // printf("res[0]: %d, res[1]: %d, res[2]: %d, index: %d\n\n", sds->res[0], sds->res[1], sds->res[2], index); + obstacles[index] = 1; + // for moving gobstacles + /* + const LbmFloat maxVelVal = 0.1666; + const LbmFloat maxusqr = maxVelVal*maxVelVal*3. *1.5; + + LbmVec objvel = vec2L((mMOIVertices[n]-mMOIVerticesOld[n]) /dvec); + { + const LbmFloat usqr = (objvel[0]*objvel[0]+objvel[1]*objvel[1]+objvel[2]*objvel[2])*1.5; + USQRMAXCHECK(usqr, objvel[0],objvel[1],objvel[2], mMaxVlen, mMxvx,mMxvy,mMxvz); + if(usqr>maxusqr) { + // cutoff at maxVelVal + for(int jj=0; jj<3; jj++) { + if(objvel[jj]>0.) objvel[jj] = maxVelVal; + if(objvel[jj]<0.) objvel[jj] = -maxVelVal; + } + } + } + const LbmFloat dp=dot(objvel, vec2L((*pNormals)[n]) ); + const LbmVec oldov=objvel; // debug + objvel = vec2L((*pNormals)[n]) *dp; + */ + } + } + } + if(sds->coll_group) + go = go->next; + else + base= base->next; + } + } +} void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedMesh *dm, int useRenderParams, int isFinalCalc) { if((smd->type & MOD_SMOKE_TYPE_FLOW)) @@ -788,16 +1101,15 @@ void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedM } else if(smd->type & MOD_SMOKE_TYPE_DOMAIN) { - PointCache *cache; + PointCache *cache = NULL; PTCacheID pid; + PointCache *cache_wt = NULL; + PTCacheID pid_wt; float timescale; - int cache_result = 0; - int startframe, endframe, framenr; + int cache_result = 0, cache_result_wt = 0; + int startframe, endframe, framenr, badloading = 0; SmokeDomainSettings *sds = smd->domain; - float light[3] = {0.0,0.0,0.0}; - int have_lamp = 0; - - // printf("smd->type & MOD_SMOKE_TYPE_DOMAIN\n"); + float light[3]; framenr = scene->r.cfra; @@ -806,51 +1118,26 @@ void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedM BKE_ptcache_id_from_smoke(&pid, ob, smd); BKE_ptcache_id_time(&pid, scene, framenr, &startframe, &endframe, ×cale); + cache_wt = sds->point_cache[1]; + BKE_ptcache_id_from_smoke_turbulence(&pid_wt, ob, smd); + /* handle continuous simulation with the play button */ if(BKE_ptcache_get_continue_physics()) { - cache->flag &= ~PTCACHE_SIMULATION_VALID; - cache->simframe= 0; - cache->last_exact= 0; - - if(!smokeModifier_init(smd, ob, scene, dm)) - return; - - if(!smd->domain->fluid) - return; - - smoke_simulate_domain(smd, scene, ob, dm); - - { - Base *base_tmp = NULL; - - for(base_tmp = scene->base.first; base_tmp; base_tmp= base_tmp->next) - { - if(base_tmp->object->type == OB_LAMP) - { - Lamp *la = (Lamp *)base_tmp->object->data; - - if(la->type == LA_LOCAL) - { - VECCOPY(light, base_tmp->object->obmat[3]); - have_lamp = 1; - break; - } - } - } - } - - smoke_prepare_View(smd, (float)framenr, light, have_lamp); - + // TODO return; } - + if(framenr < startframe) { cache->flag &= ~PTCACHE_SIMULATION_VALID; cache->simframe= 0; cache->last_exact= 0; + cache_wt->flag &= ~PTCACHE_SIMULATION_VALID; + cache_wt->simframe= 0; + cache_wt->last_exact= 0; + // we got back in time, reset smoke in this case (TODO: use cache later) // smd->time = scene->r.cfra; // smokeModifier_reset(smd); @@ -860,14 +1147,25 @@ void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedM else if(framenr > endframe) { framenr = endframe; + + // we load last valid frame here + // and don't update the smd->time variable later + badloading = 1; } if(!(cache->flag & PTCACHE_SIMULATION_VALID)) { - // printf("reseting\n"); BKE_ptcache_id_reset(scene, &pid, PTCACHE_RESET_OUTDATED); } - + if(!(cache_wt->flag & PTCACHE_SIMULATION_VALID)) + { + BKE_ptcache_id_reset(scene, &pid_wt, PTCACHE_RESET_OUTDATED); + } + + if(smd->time == -1 && framenr!= startframe) + return; + + if(!smokeModifier_init(smd, ob, scene, dm)) return; @@ -875,7 +1173,7 @@ void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedM return; /* try to read from cache */ - cache_result = BKE_ptcache_read_cache(&pid, (float)framenr, scene->r.frs_sec); + cache_result = BKE_ptcache_read_cache(&pid, (float)framenr, scene->r.frs_sec); // printf("cache_result: %d\n", cache_result); if(cache_result == PTCACHE_READ_EXACT) @@ -886,417 +1184,148 @@ void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedM cache->simframe= framenr; sds->v3dnum = framenr; + if(!badloading) + smd->time = scene->r.cfra; + + // check for wt cache + if(sds->wt) + { + cache_result_wt = BKE_ptcache_read_cache(&pid_wt, (float)framenr, scene->r.frs_sec); + // printf("cache_result_wt: %d\n", cache_result_wt); + + // error handling + if(cache_result_wt == PTCACHE_READ_EXACT) + { + cache_wt->flag |= PTCACHE_SIMULATION_VALID; + cache_wt->simframe= framenr; + } + else if(cache_result_wt==PTCACHE_READ_OLD) + { + BKE_ptcache_id_reset(scene, &pid_wt, PTCACHE_RESET_FREE); + cache_wt->flag |= PTCACHE_SIMULATION_VALID; + } + else if(ob->id.lib || (cache_wt->flag & PTCACHE_BAKED)) + { + // if baked and nothing in cache, do nothing + cache_wt->flag &= ~PTCACHE_SIMULATION_VALID; + cache_wt->simframe= 0; + cache_wt->last_exact= 0; + } + } + // printf("PTCACHE_READ_EXACT\n"); return; } else if(cache_result==PTCACHE_READ_OLD) { BKE_ptcache_id_reset(scene, &pid, PTCACHE_RESET_FREE); - - // printf("PTCACHE_READ_OLD\n"); - cache->flag |= PTCACHE_SIMULATION_VALID; + + BKE_ptcache_id_reset(scene, &pid_wt, PTCACHE_RESET_FREE); + cache_wt->flag |= PTCACHE_SIMULATION_VALID; } else if(ob->id.lib || (cache->flag & PTCACHE_BAKED)) { - /* if baked and nothing in cache, do nothing */ + // if baked and nothing in cache, do nothing cache->flag &= ~PTCACHE_SIMULATION_VALID; cache->simframe= 0; cache->last_exact= 0; + cache_wt->flag &= ~PTCACHE_SIMULATION_VALID; + cache_wt->simframe= 0; + cache_wt->last_exact= 0; + // printf("PTCACHE_BAKED\n"); return; } - else if((cache_result==0) && (startframe!=framenr) && !(cache->flag & PTCACHE_SIMULATION_VALID)) + /* + else if((cache_result==0) && ((startframe!=framenr) && !(cache->flag & PTCACHE_SIMULATION_VALID || (framenr == smd->time)))) { cache->flag &= ~PTCACHE_SIMULATION_VALID; cache->simframe= 0; cache->last_exact= 0; return; - } - + }*/ + + // printf("framenr: %d, time: %f\n", framenr, smd->time); + /* do simulation */ // low res cache->flag |= PTCACHE_SIMULATION_VALID; cache->simframe= framenr; - smoke_simulate_domain(smd, scene, ob, dm); - if(sds->wt) - smoke_turbulence_step(sds->wt, sds->fluid); - { - Base *base_tmp = NULL; - - for(base_tmp = scene->base.first; base_tmp; base_tmp= base_tmp->next) - { - if(base_tmp->object->type == OB_LAMP) - { - Lamp *la = (Lamp *)base_tmp->object->data; - - if(la->type == LA_LOCAL) - { - VECCOPY(light, base_tmp->object->obmat[3]); - have_lamp = 1; - break; - } - } - } + cache_wt->flag |= PTCACHE_SIMULATION_VALID; + cache_wt->simframe= framenr; } - - smoke_prepare_View(smd, (float)framenr, light, have_lamp); - - BKE_ptcache_write_cache(&pid, framenr); - - // printf("Writing cache_low, %d\n", framenr); - - - tend(); - // printf ( "Frame: %d, Time: %f\n", (int)smd->time, ( float ) tval() ); - } -} - -void smoke_simulate_domain(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedMesh *dm) -{ - GroupObject *go = NULL; - Base *base = NULL; - SmokeDomainSettings *sds = smd->domain; - - tstart(); - - if(sds->flags & MOD_SMOKE_DISSOLVE) - smoke_dissolve(sds->fluid, sds->diss_speed, sds->flags & MOD_SMOKE_DISSOLVE_LOG); - - // do flows and fluids - if(1) - { - Object *otherobj = NULL; - ModifierData *md = NULL; - - if(sds->fluid_group) // we use groups since we have 2 domains - go = sds->fluid_group->gobject.first; - else - base = scene->base.first; - - while(base || go) - { - otherobj = NULL; - - if(sds->fluid_group) - { - if(go->ob) - otherobj = go->ob; - } - else - otherobj = base->object; - - if(!otherobj) - { - if(sds->fluid_group) - go = go->next; - else - base= base->next; - - continue; - } - - md = modifiers_findByType(otherobj, eModifierType_Smoke); - - // check for active smoke modifier - if(md && md->mode & (eModifierMode_Realtime | eModifierMode_Render)) - { - SmokeModifierData *smd2 = (SmokeModifierData *)md; - // check for initialized smoke object - if((smd2->type & MOD_SMOKE_TYPE_FLOW) && smd2->flow) - { - // we got nice flow object - SmokeFlowSettings *sfs = smd2->flow; - - if(sfs->psys && sfs->psys->part && sfs->psys->part->type==PART_EMITTER) // is particle system selected - { - ParticleSystem *psys = sfs->psys; - ParticleSettings *part=psys->part; - ParticleData *pa = NULL; - int p = 0; - float *density = smoke_get_density(sds->fluid); - // float *bigdensity = smoke_turbulence_get_density(sds->wt); - float *heat = smoke_get_heat(sds->fluid); - float *velocity_x = smoke_get_velocity_x(sds->fluid); - float *velocity_y = smoke_get_velocity_y(sds->fluid); - float *velocity_z = smoke_get_velocity_z(sds->fluid); - unsigned char *obstacle = smoke_get_obstacle(sds->fluid); - - // debug printf("found flow psys\n"); - - // mostly copied from particle code - for(p=0, pa=psys->particles; p<psys->totpart; p++, pa++) - { - int cell[3]; - size_t i = 0; - size_t index = 0; - int badcell = 0; - - if(pa->alive == PARS_KILLED) continue; - else if(pa->alive == PARS_UNBORN && (part->flag & PART_UNBORN)==0) continue; - else if(pa->alive == PARS_DEAD && (part->flag & PART_DIED)==0) continue; - else if(pa->flag & (PARS_UNEXIST+PARS_NO_DISP)) continue; - - // VECCOPY(pos, pa->state.co); - // Mat4MulVecfl (ob->imat, pos); - - // 1. get corresponding cell - get_cell(sds->p0, sds->res, sds->dx, pa->state.co, cell, 0); - - // check if cell is valid (in the domain boundary) - for(i = 0; i < 3; i++) - { - if((cell[i] > sds->res[i] - 1) || (cell[i] < 0)) - { - badcell = 1; - break; - } - } - - if(badcell) - continue; - - // 2. set cell values (heat, density and velocity) - index = smoke_get_index(cell[0], sds->res[0], cell[1], sds->res[1], cell[2]); - - if(!(sfs->type & MOD_SMOKE_FLOW_TYPE_OUTFLOW) && !(obstacle[index] & 2)) // this is inflow - { - // heat[index] += sfs->temp * 0.1; - // density[index] += sfs->density * 0.1; - - heat[index] = sfs->temp; - density[index] = sfs->density; - - /* - velocity_x[index] = pa->state.vel[0]; - velocity_y[index] = pa->state.vel[1]; - velocity_z[index] = pa->state.vel[2]; - */ - obstacle[index] |= 2; - } - else if(sfs->type & MOD_SMOKE_FLOW_TYPE_OUTFLOW) // outflow - { - heat[index] = 0.f; - density[index] = 0.f; - velocity_x[index] = 0.f; - velocity_y[index] = 0.f; - velocity_z[index] = 0.f; - } - } - } - else - { - /* - for() - { - // no psys - BVHTreeNearest nearest; - - nearest.index = -1; - nearest.dist = FLT_MAX; - - BLI_bvhtree_find_nearest(sfs->bvh->tree, pco, &nearest, sfs->bvh->nearest_callback, sfs->bvh); - }*/ - } - } - } + tstart(); - if(sds->fluid_group) - go = go->next; - else - base= base->next; - } - } - - // do effectors - /* - if(sds->eff_group) - { - for(go = sds->eff_group->gobject.first; go; go = go->next) - { - if(go->ob) - { - if(ob->pd) - { - - } - } + if(sds->flags & MOD_SMOKE_DISSOLVE) + { + smoke_dissolve(sds->fluid, sds->diss_speed, sds->flags & MOD_SMOKE_DISSOLVE_LOG); } - } - */ - - // do collisions - if(1) - { - Object *otherobj = NULL; - ModifierData *md = NULL; - if(sds->coll_group) // we use groups since we have 2 domains - go = sds->coll_group->gobject.first; - else - base = scene->base.first; + smoke_calc_domain(scene, ob, smd); + + // set new time + smd->time = scene->r.cfra; - while(base || go) + // frame 1 is start, don't simulate anything + if(smd->time == 1) { - otherobj = NULL; - - if(sds->coll_group) - { - if(go->ob) - otherobj = go->ob; - } - else - otherobj = base->object; + // set new time + smd->time = scene->r.cfra; - if(!otherobj) - { - if(sds->coll_group) - go = go->next; - else - base= base->next; + BKE_ptcache_write_cache(&pid, framenr); + if(sds->wt) + BKE_ptcache_write_cache(&pid_wt, framenr); - continue; - } - - md = modifiers_findByType(otherobj, eModifierType_Smoke); - - // check for active smoke modifier - if(md && md->mode & (eModifierMode_Realtime | eModifierMode_Render)) - { - SmokeModifierData *smd2 = (SmokeModifierData *)md; + if(get_lamp(scene, light)) + smoke_calc_transparency(sds->shadow, smoke_get_density(sds->fluid), sds->p0, sds->p1, sds->res, sds->dx, light, calc_voxel_transp, -7.0*sds->dx); - if((smd2->type & MOD_SMOKE_TYPE_COLL) && smd2->coll) - { - // we got nice collision object - SmokeCollSettings *scs = smd2->coll; - size_t i, j; - unsigned char *obstacles = smoke_get_obstacle(smd->domain->fluid); + return; + } - for(i = 0; i < scs->numpoints; i++) - { - int badcell = 0; - size_t index = 0; - int cell[3]; + // simulate the actual smoke (c++ code in intern/smoke) + smoke_step(sds->fluid, smd->time); + BKE_ptcache_write_cache(&pid, framenr); - // 1. get corresponding cell - get_cell(sds->p0, sds->res, sds->dx, &scs->points[3 * i], cell, 0); - - // check if cell is valid (in the domain boundary) - for(j = 0; j < 3; j++) - if((cell[j] > sds->res[j] - 1) || (cell[j] < 0)) - { - badcell = 1; - break; - } - - if(badcell) - continue; + if(sds->wt) + { - // 2. set cell values (heat, density and velocity) - index = smoke_get_index(cell[0], sds->res[0], cell[1], sds->res[1], cell[2]); - - // printf("cell[0]: %d, cell[1]: %d, cell[2]: %d\n", cell[0], cell[1], cell[2]); - // printf("res[0]: %d, res[1]: %d, res[2]: %d, index: %d\n\n", sds->res[0], sds->res[1], sds->res[2], index); - - obstacles[index] = 1; + if(sds->flags & MOD_SMOKE_DISSOLVE) + smoke_dissolve_wavelet(sds->wt, sds->diss_speed, sds->flags & MOD_SMOKE_DISSOLVE_LOG); - // for moving gobstacles - /* - const LbmFloat maxVelVal = 0.1666; - const LbmFloat maxusqr = maxVelVal*maxVelVal*3. *1.5; + smoke_turbulence_step(sds->wt, sds->fluid); + BKE_ptcache_write_cache(&pid_wt, framenr); + } - LbmVec objvel = vec2L((mMOIVertices[n]-mMOIVerticesOld[n]) /dvec); { - const LbmFloat usqr = (objvel[0]*objvel[0]+objvel[1]*objvel[1]+objvel[2]*objvel[2])*1.5; - USQRMAXCHECK(usqr, objvel[0],objvel[1],objvel[2], mMaxVlen, mMxvx,mMxvy,mMxvz); - if(usqr>maxusqr) { - // cutoff at maxVelVal - for(int jj=0; jj<3; jj++) { - if(objvel[jj]>0.) objvel[jj] = maxVelVal; - if(objvel[jj]<0.) objvel[jj] = -maxVelVal; - } - } } - - const LbmFloat dp=dot(objvel, vec2L((*pNormals)[n]) ); - const LbmVec oldov=objvel; // debug - objvel = vec2L((*pNormals)[n]) *dp; - */ - } - } - } + if(get_lamp(scene, light)) + smoke_calc_transparency(sds->shadow, smoke_get_density(sds->fluid), sds->p0, sds->p1, sds->res, sds->dx, light, calc_voxel_transp, -7.0*sds->dx); - if(sds->coll_group) - go = go->next; - else - base= base->next; - } + tend(); + // printf ( "Frame: %d, Time: %f\n", (int)smd->time, ( float ) tval() ); } - - // set new time - smd->time = scene->r.cfra; - - // simulate the actual smoke (c++ code in intern/smoke) - smoke_step(sds->fluid, smd->time); } -static int calc_voxel_transp(float *input, int res[3], int *pixel, float *tRay) +static float calc_voxel_transp(float *result, float *input, int res[3], int *pixel, float *tRay, float correct) { const size_t index = smoke_get_index(pixel[0], res[0], pixel[1], res[1], pixel[2]); // T_ray *= T_vox - *tRay *= input[index*4]; - - return *tRay; -} - -// forward decleration -void smoke_calc_transparency(float *result, float *p0, float *p1, int res[3], float dx, float *light, bresenham_callback cb); - -// update necessary information for 3dview -void smoke_prepare_View(SmokeModifierData *smd, float framenr, float *light, int have_light) -{ - float *density = NULL; - int x, y, z; - size_t cells, i; - SmokeDomainSettings *sds = smd->domain; - - // update 3dview - density = smoke_get_density(smd->domain->fluid); - for(x = 0; x < smd->domain->res[0]; x++) - for(y = 0; y < smd->domain->res[1]; y++) - for(z = 0; z < smd->domain->res[2]; z++) - { - size_t index; - - index = smoke_get_index(x, smd->domain->res[0], y, smd->domain->res[1], z); - // Transparency computation - // formula taken from "Visual Simulation of Smoke" / Fedkiw et al. pg. 4 - // T_vox = exp(-C_ext * h) - // C_ext/sigma_t = density * C_ext - smd->domain->view3d[index * 4] = smd->domain->view3d[index * 4 + 1] = - smd->domain->view3d[index * 4 + 2] = exp(-density[index] * 7.0 * smd->domain->dx); - smd->domain->view3d[index * 4 + 3] = 1.0 - smd->domain->view3d[index * 4]; - - } - - if(have_light) + *tRay *= exp(input[index]*correct); + + if(result[index] < 0.0f) { - smoke_calc_transparency(sds->view3d, sds->p0, sds->p1, sds->res, sds->dx, light, calc_voxel_transp); +#pragma omp critical + result[index] = *tRay; + } - cells = smd->domain->res[0]*smd->domain->res[1]*smd->domain->res[2]; - for(i = 0; i < cells; i++) - { - smd->domain->view3d[i * 4] = smd->domain->view3d[i * 4 + 1] = - smd->domain->view3d[i * 4 + 2] = smd->domain->view3d[i * 4 + 1] * smd->domain->view3d[i * 4 + 0]; - } - } - smd->domain->v3dnum = framenr; + return *tRay; } long long smoke_get_mem_req(int xres, int yres, int zres, int amplify) @@ -1317,7 +1346,7 @@ long long smoke_get_mem_req(int xres, int yres, int zres, int amplify) return totalMB; } -static void bresenham_linie_3D(int x1, int y1, int z1, int x2, int y2, int z2, float *tRay, bresenham_callback cb, float *input, int res[3]) +static void bresenham_linie_3D(int x1, int y1, int z1, int x2, int y2, int z2, float *tRay, bresenham_callback cb, float *result, float *input, int res[3], float correct) { int dx, dy, dz, i, l, m, n, x_inc, y_inc, z_inc, err_1, err_2, dx2, dy2, dz2; int pixel[3]; @@ -1344,8 +1373,8 @@ static void bresenham_linie_3D(int x1, int y1, int z1, int x2, int y2, int z2, f err_1 = dy2 - l; err_2 = dz2 - l; for (i = 0; i < l; i++) { - if(cb(input, res, pixel, tRay) < 0.0) - return; + if(cb(result, input, res, pixel, tRay, correct) <= FLT_EPSILON) + break; if (err_1 > 0) { pixel[1] += y_inc; err_1 -= dx2; @@ -1362,8 +1391,8 @@ static void bresenham_linie_3D(int x1, int y1, int z1, int x2, int y2, int z2, f err_1 = dx2 - m; err_2 = dz2 - m; for (i = 0; i < m; i++) { - if(cb(input, res, pixel, tRay) < 0.0f) - return; + if(cb(result, input, res, pixel, tRay, correct) <= FLT_EPSILON) + break; if (err_1 > 0) { pixel[0] += x_inc; err_1 -= dy2; @@ -1380,8 +1409,8 @@ static void bresenham_linie_3D(int x1, int y1, int z1, int x2, int y2, int z2, f err_1 = dy2 - n; err_2 = dx2 - n; for (i = 0; i < n; i++) { - if(cb(input, res, pixel, tRay) < 0.0f) - return; + if(cb(result, input, res, pixel, tRay, correct) <= FLT_EPSILON) + break; if (err_1 > 0) { pixel[1] += y_inc; err_1 -= dz2; @@ -1395,7 +1424,7 @@ static void bresenham_linie_3D(int x1, int y1, int z1, int x2, int y2, int z2, f pixel[2] += z_inc; } } - cb(input, res, pixel, tRay); + cb(result, input, res, pixel, tRay, correct); } static void get_cell(float *p0, int res[3], float dx, float *pos, int *cell, int correct) @@ -1419,12 +1448,12 @@ static void get_cell(float *p0, int res[3], float dx, float *pos, int *cell, int } } -void smoke_calc_transparency(float *result, float *p0, float *p1, int res[3], float dx, float *light, bresenham_callback cb) +void smoke_calc_transparency(float *result, float *input, float *p0, float *p1, int res[3], float dx, float *light, bresenham_callback cb, float correct) { int x, y, z; float bv[6]; - // x + memset(result, -1, sizeof(float)*res[0]*res[1]*res[2]); // x bv[0] = p0[0]; bv[1] = p1[0]; // y @@ -1434,7 +1463,7 @@ void smoke_calc_transparency(float *result, float *p0, float *p1, int res[3], fl bv[4] = p0[2]; bv[5] = p1[2]; -// #pragma omp parallel for schedule(static) private(y, z) +#pragma omp parallel for schedule(static) private(y, z) for(x = 0; x < res[0]; x++) for(y = 0; y < res[1]; y++) for(z = 0; z < res[2]; z++) @@ -1447,6 +1476,8 @@ void smoke_calc_transparency(float *result, float *p0, float *p1, int res[3], fl index = smoke_get_index(x, res[0], y, res[1], z); + if(result[index] >= 0.0f) + continue; voxelCenter[0] = p0[0] + dx * x + dx * 0.5; voxelCenter[1] = p0[1] + dx * y + dx * 0.5; voxelCenter[2] = p0[2] + dx * z + dx * 0.5; @@ -1463,10 +1494,11 @@ void smoke_calc_transparency(float *result, float *p0, float *p1, int res[3], fl get_cell(p0, res, dx, light, cell, 1); } - bresenham_linie_3D(cell[0], cell[1], cell[2], x, y, z, &tRay, cb, result, res); + bresenham_linie_3D(cell[0], cell[1], cell[2], x, y, z, &tRay, cb, result, input, res, correct); // convention -> from a RGBA float array, use G value for tRay - result[index*4 + 1] = tRay; +// #pragma omp critical + result[index] = tRay; } } diff --git a/source/blender/blenloader/intern/readfile.c b/source/blender/blenloader/intern/readfile.c index 618e587d674..b6fc9f7af12 100644 --- a/source/blender/blenloader/intern/readfile.c +++ b/source/blender/blenloader/intern/readfile.c @@ -3685,8 +3685,8 @@ static void direct_link_modifiers(FileData *fd, ListBase *lb) smd->domain->smd = smd; smd->domain->fluid = NULL; - smd->domain->view3d = NULL; - smd->domain->tex = NULL; + // smd->domain->view3d = NULL; + // smd->domain->tex = NULL; direct_link_pointcache_list(fd, &(smd->domain->ptcaches[0]), &(smd->domain->point_cache[0])); direct_link_pointcache_list(fd, &(smd->domain->ptcaches[1]), &(smd->domain->point_cache[1])); diff --git a/source/blender/editors/space_view3d/drawobject.c b/source/blender/editors/space_view3d/drawobject.c index 20cab3b8aeb..7ed029f3eaf 100644 --- a/source/blender/editors/space_view3d/drawobject.c +++ b/source/blender/editors/space_view3d/drawobject.c @@ -5310,13 +5310,60 @@ void draw_object(Scene *scene, ARegion *ar, View3D *v3d, Base *base, int flag) } /* draw code for smoke */ + if((md = modifiers_findByType(ob, eModifierType_Smoke))) { - md = modifiers_findByType(ob, eModifierType_Smoke); - if (md) { - SmokeModifierData *smd = (SmokeModifierData *)md; - if(smd->type & MOD_SMOKE_TYPE_DOMAIN && smd->domain && smd->domain->fluid) { - GPU_create_smoke(smd); - draw_volume(scene, ar, v3d, base, smd->domain->tex, smd->domain->res); + SmokeModifierData *smd = (SmokeModifierData *)md; + + // draw collision objects + if((smd->type & MOD_SMOKE_TYPE_COLL) && smd->coll) + { + /*SmokeCollSettings *scs = smd->coll; + if(scs->points) + { + size_t i; + + wmLoadMatrix(rv3d->viewmat); + + if(col || (ob->flag & SELECT)) cpack(0xFFFFFF); + glDepthMask(GL_FALSE); + glEnable(GL_BLEND); + + + // glPointSize(3.0); + bglBegin(GL_POINTS); + + for(i = 0; i < scs->numpoints; i++) + { + bglVertex3fv(&scs->points[3*i]); + } + + bglEnd(); + glPointSize(1.0); + + wmMultMatrix(ob->obmat); + glDisable(GL_BLEND); + glDepthMask(GL_TRUE); + if(col) cpack(col); + + } + */ + } + + // only draw domains + if(smd->domain && smd->domain->fluid) + { + if(!smd->domain->wt || !(smd->domain->viewsettings & MOD_SMOKE_VIEW_SHOWBIG)) + { + smd->domain->tex = NULL; + GPU_create_smoke(smd, 0); + draw_volume(scene, ar, v3d, base, smd->domain->tex, smd->domain->res, smd->domain->dx, smd->domain->tex_shadow); + GPU_free_smoke(smd); + } + else if(smd->domain->wt || (smd->domain->viewsettings & MOD_SMOKE_VIEW_SHOWBIG)) + { + smd->domain->tex = NULL; + GPU_create_smoke(smd, 1); + draw_volume(scene, ar, v3d, base, smd->domain->tex, smd->domain->res_wt, smd->domain->dx_wt, smd->domain->tex_shadow); GPU_free_smoke(smd); } } diff --git a/source/blender/editors/space_view3d/drawvolume.c b/source/blender/editors/space_view3d/drawvolume.c index 342acfe92b3..c8eda10566c 100644 --- a/source/blender/editors/space_view3d/drawvolume.c +++ b/source/blender/editors/space_view3d/drawvolume.c @@ -191,7 +191,7 @@ static int larger_pow2(int n) return n*2; } -void draw_volume(Scene *scene, ARegion *ar, View3D *v3d, Base *base, GPUTexture *tex, int res[3]) +void draw_volume(Scene *scene, ARegion *ar, View3D *v3d, Base *base, GPUTexture *tex, int res[3], float dx, GPUTexture *tex_shadow) { Object *ob = base->object; RegionView3D *rv3d= ar->regiondata; @@ -204,6 +204,27 @@ void draw_volume(Scene *scene, ARegion *ar, View3D *v3d, Base *base, GPUTexture float cor[3] = {1.,1.,1.}; int gl_depth = 0, gl_blend = 0; + /* Fragment program to calculate the 3dview of smoke */ + /* using 2 textures, density and shadow */ + const char *text = "!!ARBfp1.0\n" + "PARAM dx = program.local[0];\n" + "PARAM darkness = program.local[1];\n" + "PARAM f = {1.442695041, 1.442695041, 1.442695041, 0.01};\n" + "TEMP temp, shadow, value;\n" + "TEX temp, fragment.texcoord[0], texture[0], 3D;\n" + "TEX shadow, fragment.texcoord[0], texture[1], 3D;\n" + "MUL value, temp, darkness;\n" + "MUL value, value, dx;\n" + "MUL value, value, f;\n" + "EX2 temp, -value.r;\n" + "SUB temp.a, 1.0, temp.r;\n" + "MUL temp.r, temp.r, shadow.r;\n" + "MUL temp.g, temp.g, shadow.r;\n" + "MUL temp.b, temp.b, shadow.r;\n" + "MOV result.color, temp;\n" + "END\n"; + unsigned int prog; + glGetBooleanv(GL_BLEND, (GLboolean *)&gl_blend); glGetBooleanv(GL_DEPTH_TEST, (GLboolean *)&gl_depth); @@ -234,7 +255,23 @@ void draw_volume(Scene *scene, ARegion *ar, View3D *v3d, Base *base, GPUTexture } } + if(GLEW_ARB_fragment_program) + { + glGenProgramsARB(1, &prog); + glEnable(GL_FRAGMENT_PROGRAM_ARB); + + glBindProgramARB(GL_FRAGMENT_PROGRAM_ARB, prog); + glProgramStringARB(GL_FRAGMENT_PROGRAM_ARB, GL_PROGRAM_FORMAT_ASCII_ARB, (GLsizei)strlen(text), text); + + // cell spacing + glProgramLocalParameter4fARB (GL_FRAGMENT_PROGRAM_ARB, 0, dx, dx, dx, 1.0); + // custom parameter for smoke style (higher = thicker) + glProgramLocalParameter4fARB (GL_FRAGMENT_PROGRAM_ARB, 1, 7.0, 7.0, 7.0, 1.0); + } + GPU_texture_bind(tex, 0); + if(tex_shadow) + GPU_texture_bind(tex_shadow, 1); if (!GLEW_ARB_texture_non_power_of_two) { cor[0] = (float)res[0]/(float)larger_pow2(res[0]); @@ -289,8 +326,17 @@ void draw_volume(Scene *scene, ARegion *ar, View3D *v3d, Base *base, GPUTexture n++; } + if(tex_shadow) + GPU_texture_unbind(tex_shadow); GPU_texture_unbind(tex); + if(GLEW_ARB_fragment_program) + { + glDisable(GL_FRAGMENT_PROGRAM_ARB); + glDeleteProgramsARB(1, &prog); + } + + MEM_freeN(points); if(!gl_blend) diff --git a/source/blender/editors/space_view3d/view3d_intern.h b/source/blender/editors/space_view3d/view3d_intern.h index e5e85cf9d16..00b0b5c4fd1 100644 --- a/source/blender/editors/space_view3d/view3d_intern.h +++ b/source/blender/editors/space_view3d/view3d_intern.h @@ -157,7 +157,7 @@ ARegion *view3d_has_buttons_region(ScrArea *sa); ARegion *view3d_has_tools_region(ScrArea *sa); /* draw_volume.c */ -void draw_volume(struct Scene *scene, struct ARegion *ar, struct View3D *v3d, struct Base *base, struct GPUTexture *tex, int res[3]); +void draw_volume(struct Scene *scene, struct ARegion *ar, struct View3D *v3d, struct Base *base, struct GPUTexture *tex, int res[3], float dx, struct GPUTexture *tex_shadow); #endif /* ED_VIEW3D_INTERN_H */ diff --git a/source/blender/gpu/CMakeLists.txt b/source/blender/gpu/CMakeLists.txt index 85f4233a251..279596e5ad7 100644 --- a/source/blender/gpu/CMakeLists.txt +++ b/source/blender/gpu/CMakeLists.txt @@ -28,7 +28,7 @@ FILE(GLOB SRC intern/*.c) SET(INC . ../blenlib ../blenkernel ../makesdna ../include - ../../../extern/glew/include ../../../intern/guardedalloc ../imbuf) + ../../../extern/glew/include ../../../intern/guardedalloc ../../../intern/smoke/extern ../imbuf) BLENDERLIB(bf_gpu "${SRC}" "${INC}") diff --git a/source/blender/gpu/GPU_draw.h b/source/blender/gpu/GPU_draw.h index 00d0e3131e5..fabe1420e83 100644 --- a/source/blender/gpu/GPU_draw.h +++ b/source/blender/gpu/GPU_draw.h @@ -115,7 +115,7 @@ void GPU_free_images(void); /* smoke drawing functions */ void GPU_free_smoke(struct SmokeModifierData *smd); -void GPU_create_smoke(struct SmokeModifierData *smd); +void GPU_create_smoke(struct SmokeModifierData *smd, int highres); #ifdef __cplusplus } diff --git a/source/blender/gpu/intern/gpu_draw.c b/source/blender/gpu/intern/gpu_draw.c index a81c7e03455..75e8073aafd 100644 --- a/source/blender/gpu/intern/gpu_draw.c +++ b/source/blender/gpu/intern/gpu_draw.c @@ -64,6 +64,8 @@ #include "GPU_material.h" #include "GPU_draw.h" +#include "smoke_API.h" + /* These are some obscure rendering functions shared between the * game engine and the blender, in this module to avoid duplicaten * and abstract them away from the rest a bit */ @@ -754,13 +756,21 @@ void GPU_free_smoke(SmokeModifierData *smd) if(smd->domain->tex) GPU_texture_free(smd->domain->tex); smd->domain->tex = NULL; + + if(smd->domain->tex_shadow) + GPU_texture_free(smd->domain->tex_shadow); + smd->domain->tex_shadow = NULL; } } -void GPU_create_smoke(SmokeModifierData *smd) +void GPU_create_smoke(SmokeModifierData *smd, int highres) { - if(smd->type & MOD_SMOKE_TYPE_DOMAIN && smd->domain && !smd->domain->tex) - smd->domain->tex = GPU_texture_create_3D(smd->domain->res[0], smd->domain->res[1], smd->domain->res[2], smd->domain->view3d); + if(smd->type & MOD_SMOKE_TYPE_DOMAIN && smd->domain && !smd->domain->tex && !highres) + smd->domain->tex = GPU_texture_create_3D(smd->domain->res[0], smd->domain->res[1], smd->domain->res[2], smoke_get_density(smd->domain->fluid)); + else if(smd->type & MOD_SMOKE_TYPE_DOMAIN && smd->domain && !smd->domain->tex && highres) + smd->domain->tex = GPU_texture_create_3D(smd->domain->res_wt[0], smd->domain->res_wt[1], smd->domain->res_wt[2], smoke_turbulence_get_density(smd->domain->wt)); + + smd->domain->tex_shadow = GPU_texture_create_3D(smd->domain->res[0], smd->domain->res[1], smd->domain->res[2], smd->domain->shadow); } void GPU_free_image(Image *ima) diff --git a/source/blender/gpu/intern/gpu_extensions.c b/source/blender/gpu/intern/gpu_extensions.c index 850b46dc28c..d7b54e425fd 100644 --- a/source/blender/gpu/intern/gpu_extensions.c +++ b/source/blender/gpu/intern/gpu_extensions.c @@ -346,17 +346,17 @@ GPUTexture *GPU_texture_create_3D(int w, int h, int depth, float *fpixels) tex->number = 0; glBindTexture(tex->target, tex->bindcode); - type = GL_UNSIGNED_BYTE; - format = GL_RGBA; - internalformat = GL_RGBA8; + type = GL_FLOAT; // GL_UNSIGNED_BYTE + format = GL_RED; + internalformat = GL_RED; - if (fpixels) - pixels = GPU_texture_convert_pixels(w*h*depth, fpixels); + //if (fpixels) + // pixels = GPU_texture_convert_pixels(w*h*depth, fpixels); glTexImage3D(tex->target, 0, internalformat, tex->w, tex->h, tex->depth, 0, format, type, 0); if (fpixels) { - glTexSubImage3D(tex->target, 0, 0, 0, 0, w, h, depth, format, type, pixels); + glTexSubImage3D(tex->target, 0, 0, 0, 0, w, h, depth, format, type, fpixels); } glTexParameterfv(GL_TEXTURE_3D, GL_TEXTURE_BORDER_COLOR, vfBorderColor); diff --git a/source/blender/makesdna/DNA_smoke_types.h b/source/blender/makesdna/DNA_smoke_types.h index 7c6c7fab9e4..4e4714cdaa1 100644 --- a/source/blender/makesdna/DNA_smoke_types.h +++ b/source/blender/makesdna/DNA_smoke_types.h @@ -38,9 +38,8 @@ #define MOD_SMOKE_NOISEWAVE (1<<0) #define MOD_SMOKE_NOISEFFT (1<<1) #define MOD_SMOKE_NOISECURL (1<<2) - /* viewsettings */ -#define MOD_SMOKE_SHOWHIGHRES (1<<0) /* show high resolution */ +#define MOD_SMOKE_VIEW_SHOWBIG (1<<0) typedef struct SmokeDomainSettings { struct SmokeModifierData *smd; /* for fast RNA access */ @@ -48,35 +47,34 @@ typedef struct SmokeDomainSettings { struct Group *fluid_group; struct Group *eff_group; // effector group for e.g. wind force struct Group *coll_group; // collision objects group + struct WTURBULENCE *wt; // WTURBULENCE object, if active struct GPUTexture *tex; - float *view3d; /* voxel data for display */ - unsigned int v3dnum; /* number of frame in view3d buffer */ + struct GPUTexture *tex_wt; + struct GPUTexture *tex_shadow; + float *shadow; float p0[3]; /* start point of BB */ float p1[3]; /* end point of BB */ float dx; /* edge length of one cell */ - float firstframe; - float lastframe; + float omega; /* smoke color - from 0 to 1 */ float temp; /* fluid temperature */ float tempAmb; /* ambient temperature */ float alpha; float beta; int res[3]; /* domain resolution */ + int amplify; /* wavelet amplification */ int maxres; /* longest axis on the BB gets this resolution assigned */ int flags; /* show up-res or low res, etc */ + int pad; int viewsettings; + short noise; /* noise type: wave, curl, anisotropic */ short diss_percent; - short pad; int diss_speed;/* in frames */ - struct PointCache *point_cache[2]; /* definition is in DNA_object_force.h */ - struct ListBase ptcaches[2]; - struct WTURBULENCE *wt; // WTURBULENCE object, if active - int pad3; float strength; int res_wt[3]; - int maxres_wt; - short noise; /* noise type: wave, curl, anisotropic */ - short pad2; - int amplify; + float dx_wt; + int v3dnum; + struct PointCache *point_cache[2]; /* definition is in DNA_object_force.h */ + struct ListBase ptcaches[2]; } SmokeDomainSettings; diff --git a/source/blender/makesrna/intern/rna_smoke.c b/source/blender/makesrna/intern/rna_smoke.c index b3192b110f4..943129c7169 100644 --- a/source/blender/makesrna/intern/rna_smoke.c +++ b/source/blender/makesrna/intern/rna_smoke.c @@ -46,7 +46,6 @@ #include "BKE_context.h" #include "BKE_depsgraph.h" #include "BKE_particle.h" -#include "BKE_pointcache.h" #include "ED_object.h" @@ -79,15 +78,6 @@ static void rna_Smoke_reset_dependancy(bContext *C, PointerRNA *ptr) rna_Smoke_dependency_update(C, ptr); } -#if 0 -static void rna_Smoke_redraw(bContext *C, PointerRNA *ptr) -{ - SmokeDomainSettings *settings = (SmokeDomainSettings*)ptr->data; - - settings->flags |= MOD_SMOKE_VIEW_REDRAWNICE; -} -#endif - static char *rna_SmokeDomainSettings_path(PointerRNA *ptr) { SmokeDomainSettings *settings = (SmokeDomainSettings*)ptr->data; @@ -139,6 +129,29 @@ static void rna_def_smoke_domain_settings(BlenderRNA *brna) RNA_def_property_ui_text(prop, "Max Res", "Maximal resolution used in the fluid domain."); RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset"); + prop= RNA_def_property(srna, "amplify", PROP_INT, PROP_NONE); + RNA_def_property_int_sdna(prop, NULL, "amplify"); + RNA_def_property_range(prop, 1, 10); + RNA_def_property_ui_range(prop, 1, 10, 1, 0); + RNA_def_property_ui_text(prop, "Amplification", "Enhance the resolution of smoke by this factor using noise."); + RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset"); + + prop= RNA_def_property(srna, "highres", PROP_BOOLEAN, PROP_NONE); + RNA_def_property_boolean_sdna(prop, NULL, "flags", MOD_SMOKE_HIGHRES); + RNA_def_property_ui_text(prop, "High res", "Enable high resolution (using amplification)."); + RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset"); + + prop= RNA_def_property(srna, "viewhighres", PROP_BOOLEAN, PROP_NONE); + RNA_def_property_boolean_sdna(prop, NULL, "viewsettings", MOD_SMOKE_VIEW_SHOWBIG); + RNA_def_property_ui_text(prop, "Show High Resolution", "Show high resolution (using amplification)."); + RNA_def_property_update(prop, NC_OBJECT|ND_DRAW, NULL); + + prop= RNA_def_property(srna, "noise_type", PROP_ENUM, PROP_NONE); + RNA_def_property_enum_sdna(prop, NULL, "noise"); + RNA_def_property_enum_items(prop, prop_noise_type_items); + RNA_def_property_ui_text(prop, "Noise Method", "Noise method which is used for creating the high resolution"); + RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset"); + prop= RNA_def_property(srna, "alpha", PROP_FLOAT, PROP_NONE); RNA_def_property_float_sdna(prop, NULL, "alpha"); RNA_def_property_range(prop, -5.0, 5.0); @@ -174,61 +187,36 @@ static void rna_def_smoke_domain_settings(BlenderRNA *brna) RNA_def_property_ui_text(prop, "Effector Group", "Limit effectors to this group."); RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset_dependancy"); + prop= RNA_def_property(srna, "strength", PROP_FLOAT, PROP_NONE); + RNA_def_property_float_sdna(prop, NULL, "strength"); + RNA_def_property_range(prop, 1.0, 10.0); + RNA_def_property_ui_range(prop, 1.0, 10.0, 1, 2); + RNA_def_property_ui_text(prop, "Strength", "Strength of wavelet noise"); + RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset"); + prop= RNA_def_property(srna, "dissolve_speed", PROP_INT, PROP_NONE); RNA_def_property_int_sdna(prop, NULL, "diss_speed"); RNA_def_property_range(prop, 1.0, 100.0); RNA_def_property_ui_range(prop, 1.0, 1000.0, 1, 0); RNA_def_property_ui_text(prop, "Dissolve Speed", "Dissolve Speed"); - RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset"); - - prop= RNA_def_property(srna, "highres", PROP_BOOLEAN, PROP_NONE); - RNA_def_property_boolean_sdna(prop, NULL, "flags", MOD_SMOKE_HIGHRES); - RNA_def_property_ui_text(prop, "High Resolution Smoke", "Enable high resolution smoke"); - RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, NULL); + RNA_def_property_update(prop, 0, NULL); prop= RNA_def_property(srna, "dissolve_smoke", PROP_BOOLEAN, PROP_NONE); RNA_def_property_boolean_sdna(prop, NULL, "flags", MOD_SMOKE_DISSOLVE); RNA_def_property_ui_text(prop, "Dissolve Smoke", "Enable smoke to disappear over time."); - RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset"); + RNA_def_property_update(prop, 0, NULL); prop= RNA_def_property(srna, "dissolve_smoke_log", PROP_BOOLEAN, PROP_NONE); RNA_def_property_boolean_sdna(prop, NULL, "flags", MOD_SMOKE_DISSOLVE_LOG); RNA_def_property_ui_text(prop, "Logarithmic dissolve", "Using 1/x "); - RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset"); + RNA_def_property_update(prop, 0, NULL); - prop= RNA_def_property(srna, "point_cache", PROP_POINTER, PROP_NEVER_NULL); + prop= RNA_def_property(srna, "point_cache_low", PROP_POINTER, PROP_NEVER_NULL); RNA_def_property_pointer_sdna(prop, NULL, "point_cache[0]"); - RNA_def_property_struct_type(prop, "PointCache"); RNA_def_property_ui_text(prop, "Point Cache", ""); - prop= RNA_def_property(srna, "show_highres", PROP_BOOLEAN, PROP_NONE); - RNA_def_property_boolean_sdna(prop, NULL, "viewsettings", MOD_SMOKE_SHOWHIGHRES); - RNA_def_property_ui_text(prop, "High res", "Show high resolution (using amplification)."); - RNA_def_property_update(prop, NC_OBJECT|ND_DRAW, NULL); - - prop= RNA_def_property(srna, "noise_type", PROP_ENUM, PROP_NONE); - RNA_def_property_enum_sdna(prop, NULL, "noise"); - RNA_def_property_enum_items(prop, prop_noise_type_items); - RNA_def_property_ui_text(prop, "Noise Method", "Noise method which is used for creating the high resolution"); - RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset"); - - prop= RNA_def_property(srna, "amplify", PROP_INT, PROP_NONE); - RNA_def_property_int_sdna(prop, NULL, "amplify"); - RNA_def_property_range(prop, 1, 10); - RNA_def_property_ui_range(prop, 1, 10, 1, 0); - RNA_def_property_ui_text(prop, "Amplification", "Enhance the resolution of smoke by this factor using noise."); - RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset"); - - prop= RNA_def_property(srna, "strength", PROP_FLOAT, PROP_NONE); - RNA_def_property_float_sdna(prop, NULL, "strength"); - RNA_def_property_range(prop, 1.0, 10.0); - RNA_def_property_ui_range(prop, 1.0, 10.0, 1, 2); - RNA_def_property_ui_text(prop, "Strength", "Strength of wavelet noise"); - RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset"); - - prop= RNA_def_property(srna, "point_cache_turbulence", PROP_POINTER, PROP_NEVER_NULL); + prop= RNA_def_property(srna, "point_cache_high", PROP_POINTER, PROP_NEVER_NULL); RNA_def_property_pointer_sdna(prop, NULL, "point_cache[1]"); - RNA_def_property_struct_type(prop, "PointCache"); RNA_def_property_ui_text(prop, "Point Cache", ""); } @@ -247,26 +235,26 @@ static void rna_def_smoke_flow_settings(BlenderRNA *brna) RNA_def_property_range(prop, 0.001, 1); RNA_def_property_ui_range(prop, 0.001, 1.0, 1.0, 4); RNA_def_property_ui_text(prop, "Density", ""); - RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, NULL); + RNA_def_property_update(prop, 0, NULL); // NC_OBJECT|ND_MODIFIER prop= RNA_def_property(srna, "temperature", PROP_FLOAT, PROP_NONE); RNA_def_property_float_sdna(prop, NULL, "temp"); RNA_def_property_range(prop, -10, 10); RNA_def_property_ui_range(prop, -10, 10, 1, 1); RNA_def_property_ui_text(prop, "Temp. Diff.", "Temperature difference to ambientt temperature."); - RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, NULL); + RNA_def_property_update(prop, 0, NULL); prop= RNA_def_property(srna, "psys", PROP_POINTER, PROP_NONE); RNA_def_property_pointer_sdna(prop, NULL, "psys"); RNA_def_property_struct_type(prop, "ParticleSystem"); RNA_def_property_flag(prop, PROP_EDITABLE); RNA_def_property_ui_text(prop, "Particle Systems", "Particle systems emitted from the object."); - RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset_dependancy"); + RNA_def_property_update(prop, 0, "rna_Smoke_reset_dependancy"); prop= RNA_def_property(srna, "outflow", PROP_BOOLEAN, PROP_NONE); RNA_def_property_boolean_sdna(prop, NULL, "type", MOD_SMOKE_FLOW_TYPE_OUTFLOW); RNA_def_property_ui_text(prop, "Outflow", "Deletes smoke from simulation"); - RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, NULL); + RNA_def_property_update(prop, 0, NULL); } static void rna_def_smoke_coll_settings(BlenderRNA *brna) |