diff options
Diffstat (limited to 'intern/smoke')
-rw-r--r-- | intern/smoke/extern/smoke_API.h | 15 | ||||
-rw-r--r-- | intern/smoke/intern/FLUID_3D.cpp | 79 | ||||
-rw-r--r-- | intern/smoke/intern/FLUID_3D.h | 17 | ||||
-rw-r--r-- | intern/smoke/intern/FLUID_3D_SOLVERS.cpp | 24 | ||||
-rw-r--r-- | intern/smoke/intern/WTURBULENCE.cpp | 499 | ||||
-rw-r--r-- | intern/smoke/intern/WTURBULENCE.h | 21 | ||||
-rw-r--r-- | intern/smoke/intern/smoke_API.cpp | 51 |
7 files changed, 390 insertions, 316 deletions
diff --git a/intern/smoke/extern/smoke_API.h b/intern/smoke/extern/smoke_API.h index f0dba3cc7a4..b21ce473202 100644 --- a/intern/smoke/extern/smoke_API.h +++ b/intern/smoke/extern/smoke_API.h @@ -20,7 +20,7 @@ * The Original Code is Copyright (C) 2009 by Daniel Genrich * All rights reserved. * - * Contributor(s): None + * Contributor(s): Daniel Genrich * * ***** END GPL LICENSE BLOCK ***** */ @@ -32,6 +32,10 @@ extern "C" { #endif +// export +void smoke_export(struct FLUID_3D *fluid, float *dt, float *dx, float **dens, float **densold, float **heat, float **heatold, float **vx, float **vy, float **vz, float **vxold, float **vyold, float **vzold, unsigned char **obstacles); + +// low res struct FLUID_3D *smoke_init(int *res, float *p0, float dt); void smoke_free(struct FLUID_3D *fluid); @@ -57,11 +61,14 @@ 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, int *res); +void smoke_turbulence_get_res(struct WTURBULENCE *wt, unsigned int *res); void smoke_turbulence_set_noise(struct WTURBULENCE *wt, int type); -void smoke_initWaveletBlenderRNA(struct WTURBULENCE *wt, float *strength); +void smoke_turbulence_initBlenderRNA(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); #ifdef __cplusplus } diff --git a/intern/smoke/intern/FLUID_3D.cpp b/intern/smoke/intern/FLUID_3D.cpp index ff66f29143c..89dd893198b 100644 --- a/intern/smoke/intern/FLUID_3D.cpp +++ b/intern/smoke/intern/FLUID_3D.cpp @@ -75,8 +75,6 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dt) : // allocate arrays _totalCells = _xRes * _yRes * _zRes; _slabSize = _xRes * _yRes; - _divergence = new float[_totalCells]; - _pressure = new float[_totalCells]; _xVelocity = new float[_totalCells]; _yVelocity = new float[_totalCells]; _zVelocity = new float[_totalCells]; @@ -86,20 +84,11 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dt) : _xForce = new float[_totalCells]; _yForce = new float[_totalCells]; _zForce = new float[_totalCells]; - _vorticity = new float[_totalCells]; _density = new float[_totalCells]; _densityOld = new float[_totalCells]; _heat = new float[_totalCells]; _heatOld = new float[_totalCells]; - _residual = new float[_totalCells]; - _direction = new float[_totalCells]; - _q = new float[_totalCells]; - _obstacles = new unsigned char[_totalCells]; - _xVorticity = new float[_totalCells]; - _yVorticity = new float[_totalCells]; - _zVorticity = new float[_totalCells]; - _h = new float[_totalCells]; - _Precond = new float[_totalCells]; + _obstacles = new unsigned char[_totalCells]; // set 0 at end of step // DG TODO: check if alloc went fine @@ -109,8 +98,6 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dt) : _densityOld[x] = 0.0f; _heat[x] = 0.0f; _heatOld[x] = 0.0f; - _divergence[x] = 0.0f; - _pressure[x] = 0.0f; _xVelocity[x] = 0.0f; _yVelocity[x] = 0.0f; _zVelocity[x] = 0.0f; @@ -120,19 +107,11 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dt) : _xForce[x] = 0.0f; _yForce[x] = 0.0f; _zForce[x] = 0.0f; - _xVorticity[x] = 0.0f; - _yVorticity[x] = 0.0f; - _zVorticity[x] = 0.0f; - _residual[x] = 0.0f; - _q[x] = 0.0f; - _direction[x] = 0.0f; - _h[x] = 0.0f; - _Precond[x] = 0.0f; _obstacles[x] = false; } // set side obstacles - int index; + size_t index; for (int y = 0; y < _yRes; y++) // z for (int x = 0; x < _xRes; x++) { @@ -177,8 +156,6 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dt) : FLUID_3D::~FLUID_3D() { - if (_divergence) delete[] _divergence; - if (_pressure) delete[] _pressure; if (_xVelocity) delete[] _xVelocity; if (_yVelocity) delete[] _yVelocity; if (_zVelocity) delete[] _zVelocity; @@ -188,23 +165,14 @@ FLUID_3D::~FLUID_3D() if (_xForce) delete[] _xForce; if (_yForce) delete[] _yForce; if (_zForce) delete[] _zForce; - if (_residual) delete[] _residual; - if (_direction) delete[] _direction; - if (_q) delete[] _q; if (_density) delete[] _density; if (_densityOld) delete[] _densityOld; if (_heat) delete[] _heat; if (_heatOld) delete[] _heatOld; - if (_xVorticity) delete[] _xVorticity; - if (_yVorticity) delete[] _yVorticity; - if (_zVorticity) delete[] _zVorticity; - if (_vorticity) delete[] _vorticity; - if (_h) delete[] _h; - if (_Precond) delete[] _Precond; if (_obstacles) delete[] _obstacles; // if (_wTurbulence) delete _wTurbulence; - printf("deleted fluid\n"); + // printf("deleted fluid\n"); } // init direct access functions from blender @@ -263,6 +231,8 @@ void FLUID_3D::step() */ _totalTime += _dt; _totalSteps++; + + memset(_obstacles, 0, sizeof(unsigned char)*_xRes*_yRes*_zRes); } ////////////////////////////////////////////////////////////////////// @@ -300,7 +270,7 @@ void FLUID_3D::artificialDamping(float* field) { ////////////////////////////////////////////////////////////////////// void FLUID_3D::copyBorderAll(float* field) { - int index; + size_t index; for (int y = 0; y < _yRes; y++) for (int x = 0; x < _xRes; x++) { @@ -367,9 +337,16 @@ void FLUID_3D::addForce() ////////////////////////////////////////////////////////////////////// void FLUID_3D::project() { - int index, x, y, z; + int x, y, z; + size_t index; + + float *_pressure = new float[_totalCells]; + float *_divergence = new float[_totalCells]; - setObstacleBoundaries(); + memset(_pressure, 0, sizeof(float)*_totalCells); + memset(_divergence, 0, sizeof(float)*_totalCells); + + setObstacleBoundaries(_pressure); // copy out the boundaries if(DOMAIN_BC_LEFT == 0) setNeumannX(_xVelocity, _res); @@ -414,7 +391,7 @@ void FLUID_3D::project() // solve Poisson equation solvePressurePre(_pressure, _divergence, _obstacles); - setObstaclePressure(); + setObstaclePressure(_pressure); // project out solution float invDx = 1.0f / _dx; @@ -430,6 +407,9 @@ void FLUID_3D::project() _zVelocity[index] -= 0.5f * (_pressure[index + _slabSize] - _pressure[index - _slabSize]) * invDx; } } + + if (_pressure) delete[] _pressure; + if (_divergence) delete[] _divergence; } ////////////////////////////////////////////////////////////////////// @@ -465,7 +445,7 @@ void FLUID_3D::addObstacle(OBSTACLE* obstacle) ////////////////////////////////////////////////////////////////////// // calculate the obstacle directional types ////////////////////////////////////////////////////////////////////// -void FLUID_3D::setObstaclePressure() +void FLUID_3D::setObstaclePressure(float *_pressure) { // tag remaining obstacle blocks for (int z = 1, index = _slabSize + _xRes + 1; @@ -539,7 +519,7 @@ void FLUID_3D::setObstaclePressure() } } -void FLUID_3D::setObstacleBoundaries() +void FLUID_3D::setObstacleBoundaries(float *_pressure) { // cull degenerate obstacles , move to addObstacle? for (int z = 1, index = _slabSize + _xRes + 1; @@ -600,6 +580,18 @@ void FLUID_3D::addVorticity() int x,y,z,index; if(_vorticityEps<=0.) return; + float *_xVorticity, *_yVorticity, *_zVorticity, *_vorticity; + + _xVorticity = new float[_totalCells]; + _yVorticity = new float[_totalCells]; + _zVorticity = new float[_totalCells]; + _vorticity = new float[_totalCells]; + + memset(_xVorticity, 0, sizeof(float)*_totalCells); + memset(_yVorticity, 0, sizeof(float)*_totalCells); + memset(_zVorticity, 0, sizeof(float)*_totalCells); + memset(_vorticity, 0, sizeof(float)*_totalCells); + // calculate vorticity float gridSize = 0.5f / _dx; index = _slabSize + _xRes + 1; @@ -662,6 +654,11 @@ void FLUID_3D::addVorticity() _zForce[index] += (N[0] * _yVorticity[index] - N[1] * _xVorticity[index]) * _dx * eps; } } + + if (_xVorticity) delete[] _xVorticity; + if (_yVorticity) delete[] _yVorticity; + if (_zVorticity) delete[] _zVorticity; + if (_vorticity) delete[] _vorticity; } ////////////////////////////////////////////////////////////////////// diff --git a/intern/smoke/intern/FLUID_3D.h b/intern/smoke/intern/FLUID_3D.h index 78a4cf076e3..9d9e7318204 100644 --- a/intern/smoke/intern/FLUID_3D.h +++ b/intern/smoke/intern/FLUID_3D.h @@ -64,7 +64,7 @@ class FLUID_3D // dimensions int _xRes, _yRes, _zRes, _maxRes; Vec3Int _res; - int _totalCells; + size_t _totalCells; int _slabSize; float _dx; float _p0[3]; @@ -81,7 +81,6 @@ class FLUID_3D float* _densityOld; float* _heat; float* _heatOld; - float* _pressure; float* _xVelocity; float* _yVelocity; float* _zVelocity; @@ -91,19 +90,9 @@ class FLUID_3D float* _xForce; float* _yForce; float* _zForce; - float* _divergence; - float* _xVorticity; - float* _yVorticity; - float* _zVorticity; - float* _vorticity; - float* _h; - float* _Precond; unsigned char* _obstacles; // CG fields - float* _residual; - float* _direction; - float* _q; int _iterations; // simulation constants @@ -134,8 +123,8 @@ class FLUID_3D void solveHeat(float* field, float* b, unsigned char* skip); // handle obstacle boundaries - void setObstacleBoundaries(); - void setObstaclePressure(); + void setObstacleBoundaries(float *_pressure); + void setObstaclePressure(float *_pressure); public: // advection, accessed e.g. by WTURBULENCE class diff --git a/intern/smoke/intern/FLUID_3D_SOLVERS.cpp b/intern/smoke/intern/FLUID_3D_SOLVERS.cpp index a35beaa05d7..7d078d86d61 100644 --- a/intern/smoke/intern/FLUID_3D_SOLVERS.cpp +++ b/intern/smoke/intern/FLUID_3D_SOLVERS.cpp @@ -28,10 +28,17 @@ void FLUID_3D::solvePressurePre(float* field, float* b, unsigned char* skip) { int x, y, z; size_t index; + float *_q, *_Precond, *_h, *_residual, *_direction; // i = 0 int i = 0; + _residual = new float[_totalCells]; // set 0 + _direction = new float[_totalCells]; // set 0 + _q = new float[_totalCells]; // set 0 + _h = new float[_totalCells]; // set 0 + _Precond = new float[_totalCells]; // set 0 + memset(_residual, 0, sizeof(float)*_xRes*_yRes*_zRes); memset(_q, 0, sizeof(float)*_xRes*_yRes*_zRes); memset(_direction, 0, sizeof(float)*_xRes*_yRes*_zRes); @@ -191,11 +198,18 @@ void FLUID_3D::solvePressurePre(float* field, float* b, unsigned char* skip) i++; } // cout << i << " iterations converged to " << sqrt(maxR) << endl; + + if (_h) delete[] _h; + if (_Precond) delete[] _Precond; + if (_residual) delete[] _residual; + if (_direction) delete[] _direction; + if (_q) delete[] _q; } ////////////////////////////////////////////////////////////////////// // solve the poisson equation with CG ////////////////////////////////////////////////////////////////////// +#if 0 void FLUID_3D::solvePressure(float* field, float* b, unsigned char* skip) { int x, y, z; @@ -344,6 +358,7 @@ void FLUID_3D::solvePressure(float* field, float* b, unsigned char* skip) } // cout << i << " iterations converged to " << maxR << endl; } +#endif ////////////////////////////////////////////////////////////////////// // solve the heat equation with CG @@ -353,10 +368,15 @@ void FLUID_3D::solveHeat(float* field, float* b, unsigned char* skip) int x, y, z; size_t index; const float heatConst = _dt * _heatDiffusion / (_dx * _dx); + float *_q, *_residual, *_direction; // i = 0 int i = 0; + _residual = new float[_totalCells]; // set 0 + _direction = new float[_totalCells]; // set 0 + _q = new float[_totalCells]; // set 0 + memset(_residual, 0, sizeof(float)*_xRes*_yRes*_zRes); memset(_q, 0, sizeof(float)*_xRes*_yRes*_zRes); memset(_direction, 0, sizeof(float)*_xRes*_yRes*_zRes); @@ -496,5 +516,9 @@ void FLUID_3D::solveHeat(float* field, float* b, unsigned char* skip) i++; } // cout << i << " iterations converged to " << maxR << endl; + + if (_residual) delete[] _residual; + if (_direction) delete[] _direction; + if (_q) delete[] _q; } diff --git a/intern/smoke/intern/WTURBULENCE.cpp b/intern/smoke/intern/WTURBULENCE.cpp index db7a1b55afa..dd092d4f0cc 100644 --- a/intern/smoke/intern/WTURBULENCE.cpp +++ b/intern/smoke/intern/WTURBULENCE.cpp @@ -81,25 +81,9 @@ WTURBULENCE::WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify, int no _densityBig = new float[_totalCellsBig]; _densityBigOld = new float[_totalCellsBig]; - // 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 = _tempBig2 = - _bigUx = _bigUy = _bigUz = NULL; - _tempBig1 = new float[_totalCellsBig]; - _tempBig2 = new float[_totalCellsBig]; - _bigUx = new float[_totalCellsBig]; - _bigUy = new float[_totalCellsBig]; - _bigUz = new float[_totalCellsBig]; - for(int i = 0; i < _totalCellsBig; i++) { _densityBig[i] = - _densityBigOld[i] = - _bigUx[i] = - _bigUy[i] = - _bigUz[i] = - _tempBig1[i] = - _tempBig2[i] = 0.; + _densityBigOld[i] = 0.; } // allocate & init texture coordinates @@ -154,12 +138,6 @@ WTURBULENCE::~WTURBULENCE() { delete[] _densityBig; delete[] _densityBigOld; - delete[] _bigUx; - delete[] _bigUy; - delete[] _bigUz; - delete[] _tempBig1; - delete[] _tempBig2; - delete[] _tcU; delete[] _tcV; delete[] _tcW; @@ -315,7 +293,7 @@ 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) { +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); @@ -602,12 +580,32 @@ Vec3 WTURBULENCE::WVelocityWithJacobian(Vec3 orgPos, float* xUnwarped, float* yU ////////////////////////////////////////////////////////////////////// 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); + advectTextureCoordinates(dtOrg, xvel,yvel,zvel, _tempBig1, _tempBig2); // compute eigenvalues of the texture coordinates computeEigenvalues(); @@ -744,6 +742,13 @@ void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel, IMAGE::dumpPBRT(_totalStepsBig, pbrtPrefix, _densityBig, _resBig[0],_resBig[1],_resBig[2]); */ _totalStepsBig++; + + delete[] _bigUx; + delete[] _bigUy; + delete[] _bigUz; + + delete[] _tempBig1; + delete[] _tempBig2; } ////////////////////////////////////////////////////////////////////// @@ -752,225 +757,257 @@ void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel, ////////////////////////////////////////////////////////////////////// void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, float* zvel, unsigned char *obstacles) { - // enlarge timestep to match grid - const float dt = dtOrg * _amplify; - const float invAmp = 1.0f / _amplify; + // big velocity macCormack fields + float* _bigUx; + float* _bigUy; + float* _bigUz; - // prepare textures - advectTextureCoordinates(dtOrg, xvel,yvel,zvel); + // temp arrays for BFECC and MacCormack - they have more convenient + // names in the actual implementations + float* _tempBig1; + float* _tempBig2; - // do wavelet decomposition of energy - computeEnergy(xvel, yvel, zvel, obstacles); - decomposeEnergy(); + // 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]; - // zero out coefficients inside of the obstacle - for (int x = 0; x < _totalCellsSm; x++) - if (obstacles[x]) _energy[x] = 0.f; + // enlarge timestep to match grid + const float dt = dtOrg * _amplify; + const float invAmp = 1.0f / _amplify; - // 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 + _bigUx = new float[_totalCellsBig]; + _bigUy = new float[_totalCellsBig]; + _bigUz = new float[_totalCellsBig]; - // 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; + // prepare textures + advectTextureCoordinates(dtOrg, xvel,yvel,zvel, _tempBig1, _tempBig2); + + // do wavelet decomposition of energy + computeEnergy(xvel, yvel, zvel, obstacles); + decomposeEnergy(); + + // 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 - 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; + 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); + } - // scale coefficient for next octave - amplitudeScaled *= persistence; - texCoord *= 2.0f; - } - } + // 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); - // 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 (substep < totalSubsteps - 1) + SWAP_POINTERS(_densityBig, _densityBigOld); + } // substep -#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 + // wipe the density borders + FLUID_3D::setZeroBorder(_densityBig, _resBig); - // prepare density for an advection - SWAP_POINTERS(_densityBig, _densityBigOld); + // 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(); - // 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; + // 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]); - // set boundaries of big velocity grid - FLUID_3D::setZeroX(_bigUx, _resBig); - FLUID_3D::setZeroY(_bigUy, _resBig); - FLUID_3D::setZeroZ(_bigUz, _resBig); + _totalStepsBig++; - // 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); + delete[] _bigUx; + delete[] _bigUy; + delete[] _bigUz; - if (substep < totalSubsteps - 1) - SWAP_POINTERS(_densityBig, _densityBigOld); - } // substep - - // 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(); - - // 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; } diff --git a/intern/smoke/intern/WTURBULENCE.h b/intern/smoke/intern/WTURBULENCE.h index d4e6b0c6a17..c21e002ad48 100644 --- a/intern/smoke/intern/WTURBULENCE.h +++ b/intern/smoke/intern/WTURBULENCE.h @@ -49,7 +49,7 @@ 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); + void advectTextureCoordinates (float dtOrg, float* xvel, float* yvel, float* zvel, float *_tempBig1, float *_tempBig2); void resetTextureCoordinates(); void computeEnergy(float* xvel, float* yvel, float* zvel, unsigned char *obstacles); @@ -73,7 +73,7 @@ class WTURBULENCE // is accessed on through rna gui float *_strength; - protected: + // protected: // enlargement factor from original velocity field / simulation // _Big = _amplify * _Sm int _amplify; @@ -111,26 +111,17 @@ class WTURBULENCE float* _densityBig; float* _densityBigOld; - // 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; - // texture coordinates for noise float* _tcU; float* _tcV; float* _tcW; float* _tcTemp; - float* _eigMin; - float* _eigMax; + float* _eigMin; // no save -dg + float* _eigMax; // no save -dg // wavelet decomposition of velocity energies - float* _energy; + float* _energy; // no save -dg // noise data float* _noiseTile; @@ -140,7 +131,7 @@ class WTURBULENCE int _totalStepsBig; // highest frequency component of wavelet decomposition - float* _highFreqEnergy; + float* _highFreqEnergy; // no save -dg void computeEigenvalues(); void decomposeEnergy(); diff --git a/intern/smoke/intern/smoke_API.cpp b/intern/smoke/intern/smoke_API.cpp index 2e95a576eaf..5a016b51bbe 100644 --- a/intern/smoke/intern/smoke_API.cpp +++ b/intern/smoke/intern/smoke_API.cpp @@ -20,7 +20,7 @@ * The Original Code is Copyright (C) 2009 by Daniel Genrich * All rights reserved. * - * Contributor(s): None + * Contributor(s): Daniel Genrich * * ***** END GPL LICENSE BLOCK ***** */ @@ -137,7 +137,7 @@ extern "C" void smoke_dissolve(FLUID_3D *fluid, int speed, int log) } } -extern "C" void smoke_dissolve_wavelet(WTURBULENCE *wt, int speed, int log) +extern "C" void smoke_turbulence_dissolve(WTURBULENCE *wt, int speed, int log) { float *density = wt->getDensityBig(); Vec3Int r = wt->getResBig(); @@ -172,7 +172,7 @@ extern "C" void smoke_dissolve_wavelet(WTURBULENCE *wt, int speed, int log) } } -extern "C" void smoke_initWaveletBlenderRNA(WTURBULENCE *wt, float *strength) +extern "C" void smoke_turbulence_initBlenderRNA(WTURBULENCE *wt, float *strength) { wt->initBlenderRNA(strength); } @@ -181,6 +181,36 @@ template < class T > inline T ABS( T a ) { return (0 < a) ? a : -a ; } +extern "C" void smoke_export(FLUID_3D *fluid, float *dt, float *dx, float **dens, float **densold, float **heat, float **heatold, float **vx, float **vy, float **vz, float **vxold, float **vyold, float **vzold, unsigned char **obstacles) +{ + *dens = fluid->_density; + *densold = fluid->_densityOld; + *heat = fluid->_heat; + *heatold = fluid->_heatOld; + *vx = fluid->_xVelocity; + *vy = fluid->_yVelocity; + *vz = fluid->_zVelocity; + *vxold = fluid->_xVelocityOld; + *vyold = fluid->_yVelocityOld; + *vzold = fluid->_zVelocityOld; + *obstacles = fluid->_obstacles; + dt = &(fluid->_dt); + dx = &(fluid->_dx); + +} + +extern "C" void smoke_turbulence_export(WTURBULENCE *wt, float **dens, float **densold, float **tcu, float **tcv, float **tcw) +{ + if(!wt) + return; + + *dens = wt->_densityBig; + *densold = wt->_densityBigOld; + *tcu = wt->_tcU; + *tcv = wt->_tcV; + *tcw = wt->_tcW; +} + extern "C" float *smoke_get_density(FLUID_3D *fluid) { return fluid->_density; @@ -193,17 +223,17 @@ extern "C" float *smoke_get_heat(FLUID_3D *fluid) extern "C" float *smoke_get_velocity_x(FLUID_3D *fluid) { - return fluid->_xVorticity; + return fluid->_xVelocity; } extern "C" float *smoke_get_velocity_y(FLUID_3D *fluid) { - return fluid->_yVorticity; + return fluid->_yVelocity; } extern "C" float *smoke_get_velocity_z(FLUID_3D *fluid) { - return fluid->_zVorticity; + return fluid->_zVelocity; } extern "C" float *smoke_turbulence_get_density(WTURBULENCE *wt) @@ -211,14 +241,13 @@ extern "C" float *smoke_turbulence_get_density(WTURBULENCE *wt) return wt ? wt->getDensityBig() : NULL; } -extern "C" void smoke_turbulence_get_res(WTURBULENCE *wt, int *res) +extern "C" void smoke_turbulence_get_res(WTURBULENCE *wt, unsigned int *res) { if(wt) { - Vec3Int r = wt->getResBig(); - res[0] = r[0]; - res[1] = r[1]; - res[2] = r[2]; + res[0] = wt->_resBig[0]; + res[1] = wt->_resBig[1]; + res[2] = wt->_resBig[2]; } } |