From cb634b910010c04543cb3361f7a16a261e5b9f89 Mon Sep 17 00:00:00 2001 From: Daniel Genrich Date: Wed, 10 Oct 2012 13:18:07 +0000 Subject: Google Summer of Code project: "Smoke Simulator Improvements & Fire". Documentation & Test blend files: ------------------ http://wiki.blender.org/index.php/User:MiikaH/GSoC-2012-Smoke-Simulator-Improvements Credits: ------------------ Miika Hamalainen (MiikaH): Student / Main programmer Daniel Genrich (Genscher): Mentor / Programmer of merged patches from Smoke2 branch Google: For Google Summer of Code 2012 --- intern/smoke/CMakeLists.txt | 2 + intern/smoke/extern/smoke_API.h | 51 +++- intern/smoke/intern/FLUID_3D.cpp | 394 +++++++++++++++++++++-------- intern/smoke/intern/FLUID_3D.h | 46 +++- intern/smoke/intern/FLUID_3D_SOLVERS.cpp | 1 - intern/smoke/intern/FLUID_3D_STATIC.cpp | 28 --- intern/smoke/intern/WTURBULENCE.cpp | 151 +++++++++++- intern/smoke/intern/WTURBULENCE.h | 19 +- intern/smoke/intern/smoke_API.cpp | 364 ++++++++++++++++++++------- intern/smoke/intern/spectrum.cpp | 411 +++++++++++++++++++++++++++++++ intern/smoke/intern/spectrum.h | 6 + 11 files changed, 1231 insertions(+), 242 deletions(-) create mode 100644 intern/smoke/intern/spectrum.cpp create mode 100644 intern/smoke/intern/spectrum.h (limited to 'intern/smoke') diff --git a/intern/smoke/CMakeLists.txt b/intern/smoke/CMakeLists.txt index 9524f7b2935..11f173a6835 100644 --- a/intern/smoke/CMakeLists.txt +++ b/intern/smoke/CMakeLists.txt @@ -40,6 +40,7 @@ set(SRC intern/FLUID_3D_SOLVERS.cpp intern/FLUID_3D_STATIC.cpp intern/LU_HELPER.cpp + intern/spectrum.cpp intern/SPHERE.cpp intern/WTURBULENCE.cpp intern/smoke_API.cpp @@ -53,6 +54,7 @@ set(SRC intern/LU_HELPER.h intern/MERSENNETWISTER.h intern/OBSTACLE.h + intern/spectrum.h intern/SPHERE.h intern/VEC3.h intern/WAVELET_NOISE.h diff --git a/intern/smoke/extern/smoke_API.h b/intern/smoke/extern/smoke_API.h index a0eb1bf38e0..98c63f08fbd 100644 --- a/intern/smoke/extern/smoke_API.h +++ b/intern/smoke/extern/smoke_API.h @@ -37,17 +37,23 @@ extern "C" { struct FLUID_3D; -// 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 dtdef); +struct FLUID_3D *smoke_init(int *res, float dx, float dtdef, int use_heat, int use_fire, int use_colors); void smoke_free(struct FLUID_3D *fluid); -void smoke_initBlenderRNA(struct FLUID_3D *fluid, float *alpha, float *beta, float *dt_factor, float *vorticity, int *border_colli); -void smoke_step(struct FLUID_3D *fluid, float dtSubdiv); +void smoke_initBlenderRNA(struct FLUID_3D *fluid, float *alpha, float *beta, float *dt_factor, float *vorticity, int *border_colli, float *burning_rate, + float *flame_smoke, float *flame_smoke_color, float *flame_vorticity, float *flame_ignition_temp, float *flame_max_temp); +void smoke_step(struct FLUID_3D *fluid, float gravity[3], float dtSubdiv); float *smoke_get_density(struct FLUID_3D *fluid); +float *smoke_get_flame(struct FLUID_3D *fluid); +float *smoke_get_fuel(struct FLUID_3D *fluid); +float *smoke_get_react(struct FLUID_3D *fluid); +float *smoke_get_color_r(struct FLUID_3D *fluid); +float *smoke_get_color_g(struct FLUID_3D *fluid); +float *smoke_get_color_b(struct FLUID_3D *fluid); +void smoke_get_rgba(struct FLUID_3D *fluid, float *data, int sequential); +void smoke_get_rgba_from_density(struct FLUID_3D *fluid, float color[3], float *data, int sequential); float *smoke_get_heat(struct FLUID_3D *fluid); float *smoke_get_velocity_x(struct FLUID_3D *fluid); float *smoke_get_velocity_y(struct FLUID_3D *fluid); @@ -68,19 +74,44 @@ size_t smoke_get_index2d(int x, int max_x, int y); void smoke_dissolve(struct FLUID_3D *fluid, int speed, int log); // wavelet turbulence functions -struct WTURBULENCE *smoke_turbulence_init(int *res, int amplify, int noisetype); +struct WTURBULENCE *smoke_turbulence_init(int *res, int amplify, int noisetype, int use_fire, int use_colors); 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); +float *smoke_turbulence_get_color_r(struct WTURBULENCE *wt); +float *smoke_turbulence_get_color_g(struct WTURBULENCE *wt); +float *smoke_turbulence_get_color_b(struct WTURBULENCE *wt); +void smoke_turbulence_get_rgba(struct WTURBULENCE *wt, float *data, int sequential); +void smoke_turbulence_get_rgba_from_density(struct WTURBULENCE *wt, float color[3], float *data, int sequential); +float *smoke_turbulence_get_flame(struct WTURBULENCE *wt); +float *smoke_turbulence_get_fuel(struct WTURBULENCE *wt); +float *smoke_turbulence_get_react(struct WTURBULENCE *wt); void smoke_turbulence_get_res(struct WTURBULENCE *wt, int *res); +int smoke_turbulence_get_cells(struct WTURBULENCE *wt); void smoke_turbulence_set_noise(struct WTURBULENCE *wt, int type); void smoke_initWaveletBlenderRNA(struct WTURBULENCE *wt, float *strength); - 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); +/* export */ +void smoke_export(struct FLUID_3D *fluid, float *dt, float *dx, float **dens, float **react, float **flame, float **fuel, float **heat, float **heatold, + float **vx, float **vy, float **vz, float **r, float **g, float **b, unsigned char **obstacles); +void smoke_turbulence_export(struct WTURBULENCE *wt, float **dens, float **react, float **flame, float **fuel, + float **r, float **g, float **b, float **tcu, float **tcv, float **tcw); + +/* flame spectrum */ +void flame_get_spectrum(unsigned char *spec, int width, float t1, float t2); + +/* data fields */ +int smoke_has_heat(struct FLUID_3D *fluid); +int smoke_has_fuel(struct FLUID_3D *fluid); +int smoke_has_colors(struct FLUID_3D *fluid); +int smoke_turbulence_has_fuel(struct WTURBULENCE *wt); +int smoke_turbulence_has_colors(struct WTURBULENCE *wt); + +void smoke_ensure_heat(struct FLUID_3D *fluid); +void smoke_ensure_fire(struct FLUID_3D *fluid, struct WTURBULENCE *wt); +void smoke_ensure_colors(struct FLUID_3D *fluid, struct WTURBULENCE *wt, float init_r, float init_g, float init_b); #ifdef __cplusplus } diff --git a/intern/smoke/intern/FLUID_3D.cpp b/intern/smoke/intern/FLUID_3D.cpp index e006132ea8f..b7b353352e2 100644 --- a/intern/smoke/intern/FLUID_3D.cpp +++ b/intern/smoke/intern/FLUID_3D.cpp @@ -44,16 +44,11 @@ // Construction/Destruction ////////////////////////////////////////////////////////////////////// -FLUID_3D::FLUID_3D(int *res, float *p0, float dtdef) : +FLUID_3D::FLUID_3D(int *res, float dx, float dtdef, int init_heat, int init_fire, int init_colors) : _xRes(res[0]), _yRes(res[1]), _zRes(res[2]), _res(0.0f) { // set simulation consts _dt = dtdef; // just in case. set in step from a RNA factor - - // start point of array - _p0[0] = p0[0]; - _p0[1] = p0[1]; - _p0[2] = p0[2]; _iterations = 100; _tempAmb = 0; @@ -72,7 +67,10 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dtdef) : */ // scale the constants according to the refinement of the grid - _dx = 1.0f / (float)_maxRes; + if (!dx) + _dx = 1.0f / (float)_maxRes; + else + _dx = dx; _constantScaling = 64.0f / _maxRes; _constantScaling = (_constantScaling < 1.0f) ? 1.0f : _constantScaling; _vorticityEps = 2.0f / _constantScaling; // Just in case set a default value @@ -94,8 +92,6 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dtdef) : _zForce = new float[_totalCells]; _density = new float[_totalCells]; _densityOld = new float[_totalCells]; - _heat = new float[_totalCells]; - _heatOld = new float[_totalCells]; _obstacles = new unsigned char[_totalCells]; // set 0 at end of step // For threaded version: @@ -103,7 +99,6 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dtdef) : _yVelocityTemp = new float[_totalCells]; _zVelocityTemp = new float[_totalCells]; _densityTemp = new float[_totalCells]; - _heatTemp = new float[_totalCells]; // DG TODO: check if alloc went fine @@ -111,8 +106,6 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dtdef) : { _density[x] = 0.0f; _densityOld[x] = 0.0f; - _heat[x] = 0.0f; - _heatOld[x] = 0.0f; _xVelocity[x] = 0.0f; _yVelocity[x] = 0.0f; _zVelocity[x] = 0.0f; @@ -128,6 +121,25 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dtdef) : _obstacles[x] = false; } + /* heat */ + _heat = _heatOld = _heatTemp = NULL; + if (init_heat) { + initHeat(); + } + // Fire simulation + _flame = _fuel = _fuelTemp = _fuelOld = NULL; + _react = _reactTemp = _reactOld = NULL; + if (init_fire) { + initFire(); + } + // Smoke color + _color_r = _color_rOld = _color_rTemp = NULL; + _color_g = _color_gOld = _color_gTemp = NULL; + _color_b = _color_bOld = _color_bTemp = NULL; + if (init_colors) { + initColors(0.0f, 0.0f, 0.0f); + } + // boundary conditions of the fluid domain // set default values -> vertically non-colliding _domainBcFront = true; @@ -138,9 +150,70 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dtdef) : _domainBcRight = _domainBcLeft; _colloPrev = 1; // default value +} - setBorderObstacles(); // walls +void FLUID_3D::initHeat() +{ + if (!_heat) { + _heat = new float[_totalCells]; + _heatOld = new float[_totalCells]; + _heatTemp = new float[_totalCells]; + for (int x = 0; x < _totalCells; x++) + { + _heat[x] = 0.0f; + _heatOld[x] = 0.0f; + } + } +} + +void FLUID_3D::initFire() +{ + if (!_flame) { + _flame = new float[_totalCells]; + _fuel = new float[_totalCells]; + _fuelTemp = new float[_totalCells]; + _fuelOld = new float[_totalCells]; + _react = new float[_totalCells]; + _reactTemp = new float[_totalCells]; + _reactOld = new float[_totalCells]; + + for (int x = 0; x < _totalCells; x++) + { + _flame[x] = 0.0f; + _fuel[x] = 0.0f; + _fuelTemp[x] = 0.0f; + _fuelOld[x] = 0.0f; + _react[x] = 0.0f; + _reactTemp[x] = 0.0f; + _reactOld[x] = 0.0f; + } + } +} + +void FLUID_3D::initColors(float init_r, float init_g, float init_b) +{ + if (!_color_r) { + _color_r = new float[_totalCells]; + _color_rOld = new float[_totalCells]; + _color_rTemp = new float[_totalCells]; + _color_g = new float[_totalCells]; + _color_gOld = new float[_totalCells]; + _color_gTemp = new float[_totalCells]; + _color_b = new float[_totalCells]; + _color_bOld = new float[_totalCells]; + _color_bTemp = new float[_totalCells]; + + for (int x = 0; x < _totalCells; x++) + { + _color_r[x] = _density[x] * init_r; + _color_rOld[x] = 0.0f; + _color_g[x] = _density[x] * init_g; + _color_gOld[x] = 0.0f; + _color_b[x] = _density[x] * init_b; + _color_bOld[x] = 0.0f; + } + } } void FLUID_3D::setBorderObstacles() @@ -204,7 +277,6 @@ FLUID_3D::~FLUID_3D() if (_heat) delete[] _heat; if (_heatOld) delete[] _heatOld; if (_obstacles) delete[] _obstacles; - // if (_wTurbulence) delete _wTurbulence; if (_xVelocityTemp) delete[] _xVelocityTemp; if (_yVelocityTemp) delete[] _yVelocityTemp; @@ -212,23 +284,48 @@ FLUID_3D::~FLUID_3D() if (_densityTemp) delete[] _densityTemp; if (_heatTemp) delete[] _heatTemp; + if (_flame) delete[] _flame; + if (_fuel) delete[] _fuel; + if (_fuelTemp) delete[] _fuelTemp; + if (_fuelOld) delete[] _fuelOld; + if (_react) delete[] _react; + if (_reactTemp) delete[] _reactTemp; + if (_reactOld) delete[] _reactOld; + + if (_color_r) delete[] _color_r; + if (_color_rOld) delete[] _color_rOld; + if (_color_rTemp) delete[] _color_rTemp; + if (_color_g) delete[] _color_g; + if (_color_gOld) delete[] _color_gOld; + if (_color_gTemp) delete[] _color_gTemp; + if (_color_b) delete[] _color_b; + if (_color_bOld) delete[] _color_bOld; + if (_color_bTemp) delete[] _color_bTemp; + // printf("deleted fluid\n"); } // init direct access functions from blender -void FLUID_3D::initBlenderRNA(float *alpha, float *beta, float *dt_factor, float *vorticity, int *borderCollision) +void FLUID_3D::initBlenderRNA(float *alpha, float *beta, float *dt_factor, float *vorticity, int *borderCollision, float *burning_rate, + float *flame_smoke, float *flame_smoke_color, float *flame_vorticity, float *flame_ignition_temp, float *flame_max_temp) { _alpha = alpha; _beta = beta; _dtFactor = dt_factor; _vorticityRNA = vorticity; _borderColli = borderCollision; + _burning_rate = burning_rate; + _flame_smoke = flame_smoke; + _flame_smoke_color = flame_smoke_color; + _flame_vorticity = flame_vorticity; + _ignition_temp = flame_ignition_temp; + _max_temp = flame_max_temp; } ////////////////////////////////////////////////////////////////////// // step simulation once ////////////////////////////////////////////////////////////////////// -void FLUID_3D::step(float dt) +void FLUID_3D::step(float dt, float gravity[3]) { #if 0 // If border rules have been changed @@ -281,7 +378,7 @@ void FLUID_3D::step(float dt) wipeBoundariesSL(zBegin, zEnd); addVorticity(zBegin, zEnd); - addBuoyancy(_heat, _density, zBegin, zEnd); + addBuoyancy(_heat, _density, gravity, zBegin, zEnd); addForce(zBegin, zEnd); #if PARALLEL==1 @@ -312,13 +409,15 @@ void FLUID_3D::step(float dt) if (i==0) { #endif - project(); + project(); #if PARALLEL==1 } - else + else if (i==1) { #endif - diffuseHeat(); + if (_heat) { + diffuseHeat(); + } #if PARALLEL==1 } } @@ -338,6 +437,13 @@ void FLUID_3D::step(float dt) SWAP_POINTERS(_density, _densityOld); SWAP_POINTERS(_heat, _heatOld); + SWAP_POINTERS(_fuel, _fuelOld); + SWAP_POINTERS(_react, _reactOld); + + SWAP_POINTERS(_color_r, _color_rOld); + SWAP_POINTERS(_color_g, _color_gOld); + SWAP_POINTERS(_color_b, _color_bOld); + advectMacCormackBegin(0, _zRes); #if PARALLEL==1 @@ -398,11 +504,8 @@ void FLUID_3D::step(float dt) SWAP_POINTERS(_yVelocity, _yForce); SWAP_POINTERS(_zVelocity, _zForce); - - - _totalTime += _dt; - _totalSteps++; + _totalSteps++; for (int i = 0; i < _totalCells; i++) { @@ -643,6 +746,15 @@ void FLUID_3D::wipeBoundaries(int zBegin, int zEnd) setZeroBorder(_yVelocity, _res, zBegin, zEnd); setZeroBorder(_zVelocity, _res, zBegin, zEnd); setZeroBorder(_density, _res, zBegin, zEnd); + if (_fuel) { + setZeroBorder(_fuel, _res, zBegin, zEnd); + setZeroBorder(_react, _res, zBegin, zEnd); + } + if (_color_r) { + setZeroBorder(_color_r, _res, zBegin, zEnd); + setZeroBorder(_color_g, _res, zBegin, zEnd); + setZeroBorder(_color_b, _res, zBegin, zEnd); + } } void FLUID_3D::wipeBoundariesSL(int zBegin, int zEnd) @@ -668,6 +780,15 @@ void FLUID_3D::wipeBoundariesSL(int zBegin, int zEnd) _yVelocity[index] = 0.0f; _zVelocity[index] = 0.0f; _density[index] = 0.0f; + if (_fuel) { + _fuel[index] = 0.0f; + _react[index] = 0.0f; + } + if (_color_r) { + _color_r[index] = 0.0f; + _color_g[index] = 0.0f; + _color_b[index] = 0.0f; + } // right slab index += _xRes - 1; @@ -675,6 +796,15 @@ void FLUID_3D::wipeBoundariesSL(int zBegin, int zEnd) _yVelocity[index] = 0.0f; _zVelocity[index] = 0.0f; _density[index] = 0.0f; + if (_fuel) { + _fuel[index] = 0.0f; + _react[index] = 0.0f; + } + if (_color_r) { + _color_r[index] = 0.0f; + _color_g[index] = 0.0f; + _color_b[index] = 0.0f; + } } ///////////////////////////////////// @@ -690,6 +820,15 @@ void FLUID_3D::wipeBoundariesSL(int zBegin, int zEnd) _yVelocity[index] = 0.0f; _zVelocity[index] = 0.0f; _density[index] = 0.0f; + if (_fuel) { + _fuel[index] = 0.0f; + _react[index] = 0.0f; + } + if (_color_r) { + _color_r[index] = 0.0f; + _color_g[index] = 0.0f; + _color_b[index] = 0.0f; + } // top slab index += slabSize - _xRes; @@ -697,6 +836,15 @@ void FLUID_3D::wipeBoundariesSL(int zBegin, int zEnd) _yVelocity[index] = 0.0f; _zVelocity[index] = 0.0f; _density[index] = 0.0f; + if (_fuel) { + _fuel[index] = 0.0f; + _react[index] = 0.0f; + } + if (_color_r) { + _color_r[index] = 0.0f; + _color_g[index] = 0.0f; + _color_b[index] = 0.0f; + } } @@ -717,6 +865,15 @@ void FLUID_3D::wipeBoundariesSL(int zBegin, int zEnd) _yVelocity[index] = 0.0f; _zVelocity[index] = 0.0f; _density[index] = 0.0f; + if (_fuel) { + _fuel[index] = 0.0f; + _react[index] = 0.0f; + } + if (_color_r) { + _color_r[index] = 0.0f; + _color_g[index] = 0.0f; + _color_b[index] = 0.0f; + } } if (zEnd == _zRes) @@ -735,6 +892,15 @@ void FLUID_3D::wipeBoundariesSL(int zBegin, int zEnd) _yVelocity[indexx] = 0.0f; _zVelocity[indexx] = 0.0f; _density[indexx] = 0.0f; + if (_fuel) { + _fuel[index] = 0.0f; + _react[index] = 0.0f; + } + if (_color_r) { + _color_r[index] = 0.0f; + _color_g[index] = 0.0f; + _color_b[index] = 0.0f; + } } } @@ -781,35 +947,6 @@ void FLUID_3D::project() if(!_domainBcTop) setNeumannZ(_zVelocity, _res, 0, _zRes); else setZeroZ(_zVelocity, _res, 0, _zRes); - /* - { - float maxx = 0, maxy = 0, maxz = 0; - for(unsigned int i = 0; i < _xRes * _yRes * _zRes; i++) - { - if(_xVelocity[i] > maxx) - maxx = _xVelocity[i]; - if(_yVelocity[i] > maxy) - maxy = _yVelocity[i]; - if(_zVelocity[i] > maxz) - maxz = _zVelocity[i]; - } - printf("Max velx: %f, vely: %f, velz: %f\n", maxx, maxy, maxz); - } - */ - - /* - { - float maxvalue = 0; - for(unsigned int i = 0; i < _xRes * _yRes * _zRes; i++) - { - if(_heat[i] > maxvalue) - maxvalue = _heat[i]; - - } - printf("Max heat: %f\n", maxvalue); - } - */ - // calculate divergence index = _slabSize + _xRes + 1; for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes) @@ -1007,7 +1144,6 @@ void FLUID_3D::diffuseHeat() for (int x = 0; x < _totalCells; x++) if (_obstacles[x]) _heat[x] = 0.0f; - } ////////////////////////////////////////////////////////////////////// @@ -1175,7 +1311,7 @@ void FLUID_3D::setObstacleBoundaries(float *_pressure, int zBegin, int zEnd) ////////////////////////////////////////////////////////////////////// // add buoyancy forces ////////////////////////////////////////////////////////////////////// -void FLUID_3D::addBuoyancy(float *heat, float *density, int zBegin, int zEnd) +void FLUID_3D::addBuoyancy(float *heat, float *density, float gravity[3], int zBegin, int zEnd) { int index = zBegin*_slabSize; @@ -1183,7 +1319,10 @@ void FLUID_3D::addBuoyancy(float *heat, float *density, int zBegin, int zEnd) for (int y = 0; y < _yRes; y++) for (int x = 0; x < _xRes; x++, index++) { - _zForce[index] += *_alpha * density[index] + (*_beta * (heat[index] - _tempAmb)); // DG: was _yForce, changed for Blender + float buoyancy = *_alpha * density[index] + (*_beta * (((heat) ? heat[index] : 0.0f) - _tempAmb)); + _xForce[index] -= gravity[0] * buoyancy; + _yForce[index] -= gravity[1] * buoyancy; + _zForce[index] -= gravity[2] * buoyancy; } } @@ -1192,8 +1331,10 @@ void FLUID_3D::addBuoyancy(float *heat, float *density, int zBegin, int zEnd) ////////////////////////////////////////////////////////////////////// void FLUID_3D::addVorticity(int zBegin, int zEnd) { + // set flame vorticity from RNA value + float flame_vorticity = (*_flame_vorticity)/_constantScaling; //int x,y,z,index; - if(_vorticityEps<=0.) return; + if(_vorticityEps+flame_vorticity<=0.) return; int _blockSize=zEnd-zBegin; int _blockTotalCells = _slabSize * (_blockSize+2); @@ -1300,14 +1441,15 @@ void FLUID_3D::addVorticity(int zBegin, int zEnd) float magnitude = sqrtf(N[0] * N[0] + N[1] * N[1] + N[2] * N[2]); if (magnitude > FLT_EPSILON) { + float flame_vort = (_fuel) ? _fuel[index]*flame_vorticity : 0.0f; magnitude = 1.0f / magnitude; N[0] *= magnitude; N[1] *= magnitude; N[2] *= magnitude; - _xForce[index] += (N[1] * _zVorticity[vIndex] - N[2] * _yVorticity[vIndex]) * _dx * eps; - _yForce[index] += (N[2] * _xVorticity[vIndex] - N[0] * _zVorticity[vIndex]) * _dx * eps; - _zForce[index] += (N[0] * _yVorticity[vIndex] - N[1] * _xVorticity[vIndex]) * _dx * eps; + _xForce[index] += (N[1] * _zVorticity[vIndex] - N[2] * _yVorticity[vIndex]) * _dx * (eps + flame_vort); + _yForce[index] += (N[2] * _xVorticity[vIndex] - N[0] * _zVorticity[vIndex]) * _dx * (eps + flame_vort); + _zForce[index] += (N[0] * _yVorticity[vIndex] - N[1] * _xVorticity[vIndex]) * _dx * (eps + flame_vort); } } // if vIndex++; @@ -1328,14 +1470,9 @@ void FLUID_3D::advectMacCormackBegin(int zBegin, int zEnd) { Vec3Int res = Vec3Int(_xRes,_yRes,_zRes); - if(!_domainBcLeft) copyBorderX(_xVelocityOld, res, zBegin, zEnd); - else setZeroX(_xVelocityOld, res, zBegin, zEnd); - - if(!_domainBcFront) copyBorderY(_yVelocityOld, res, zBegin, zEnd); - else setZeroY(_yVelocityOld, res, zBegin, zEnd); - - if(!_domainBcTop) copyBorderZ(_zVelocityOld, res, zBegin, zEnd); - else setZeroZ(_zVelocityOld, res, zBegin, zEnd); + setZeroX(_xVelocityOld, res, zBegin, zEnd); + setZeroY(_yVelocityOld, res, zBegin, zEnd); + setZeroZ(_zVelocityOld, res, zBegin, zEnd); } ////////////////////////////////////////////////////////////////////// @@ -1355,7 +1492,18 @@ void FLUID_3D::advectMacCormackEnd1(int zBegin, int zEnd) // advectFieldMacCormack1(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res) advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _densityOld, _densityTemp, res, zBegin, zEnd); - advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _heatOld, _heatTemp, res, zBegin, zEnd); + if (_heat) { + advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _heatOld, _heatTemp, res, zBegin, zEnd); + } + if (_fuel) { + advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _fuelOld, _fuelTemp, res, zBegin, zEnd); + advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _reactOld, _reactTemp, res, zBegin, zEnd); + } + if (_color_r) { + advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _color_rOld, _color_rTemp, res, zBegin, zEnd); + advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _color_gOld, _color_gTemp, res, zBegin, zEnd); + advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _color_bOld, _color_bTemp, res, zBegin, zEnd); + } advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _xVelocityOld, _xVelocity, res, zBegin, zEnd); advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _yVelocityOld, _yVelocity, res, zBegin, zEnd); advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _zVelocityOld, _zVelocity, res, zBegin, zEnd); @@ -1371,17 +1519,30 @@ void FLUID_3D::advectMacCormackEnd2(int zBegin, int zEnd) const float dt0 = _dt / _dx; Vec3Int res = Vec3Int(_xRes,_yRes,_zRes); - // use force array as temp arrays + // use force array as temp array float* t1 = _xForce; // advectFieldMacCormack2(dt, xVelocity, yVelocity, zVelocity, oldField, newField, tempfield, temp, res, obstacles) + /* finish advection */ advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _densityOld, _density, _densityTemp, t1, res, _obstacles, zBegin, zEnd); - advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _heatOld, _heat, _heatTemp, t1, res, _obstacles, zBegin, zEnd); + if (_heat) { + advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _heatOld, _heat, _heatTemp, t1, res, _obstacles, zBegin, zEnd); + } + if (_fuel) { + advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _fuelOld, _fuel, _fuelTemp, t1, res, _obstacles, zBegin, zEnd); + advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _reactOld, _react, _reactTemp, t1, res, _obstacles, zBegin, zEnd); + } + if (_color_r) { + advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _color_rOld, _color_r, _color_rTemp, t1, res, _obstacles, zBegin, zEnd); + advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _color_gOld, _color_g, _color_gTemp, t1, res, _obstacles, zBegin, zEnd); + advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _color_bOld, _color_b, _color_bTemp, t1, res, _obstacles, zBegin, zEnd); + } advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _xVelocityOld, _xVelocityTemp, _xVelocity, t1, res, _obstacles, zBegin, zEnd); advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _yVelocityOld, _yVelocityTemp, _yVelocity, t1, res, _obstacles, zBegin, zEnd); advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _zVelocityOld, _zVelocityTemp, _zVelocity, t1, res, _obstacles, zBegin, zEnd); + /* set boundary conditions for velocity */ if(!_domainBcLeft) copyBorderX(_xVelocityTemp, res, zBegin, zEnd); else setZeroX(_xVelocityTemp, res, zBegin, zEnd); @@ -1391,40 +1552,71 @@ void FLUID_3D::advectMacCormackEnd2(int zBegin, int zEnd) if(!_domainBcTop) copyBorderZ(_zVelocityTemp, res, zBegin, zEnd); else setZeroZ(_zVelocityTemp, res, zBegin, zEnd); + /* clear data boundaries */ setZeroBorder(_density, res, zBegin, zEnd); - setZeroBorder(_heat, res, zBegin, zEnd); -#if 0 - { - const size_t index_ = _slabSize + _xRes + 1; - int bb=0; - int bt=0; + if (_fuel) { + setZeroBorder(_fuel, res, zBegin, zEnd); + setZeroBorder(_react, res, zBegin, zEnd); + } + if (_color_r) { + setZeroBorder(_color_r, res, zBegin, zEnd); + setZeroBorder(_color_g, res, zBegin, zEnd); + setZeroBorder(_color_b, res, zBegin, zEnd); + } +} - if (zBegin == 0) {bb = 1;} - if (zEnd == _zRes) {bt = 1;} - - for (int z = zBegin + bb; z < zEnd - bt; z++) + +void FLUID_3D::processBurn(float *fuel, float *smoke, float *react, float *flame, float *heat, + float *r, float *g, float *b, int total_cells, float dt) +{ + float burning_rate = *_burning_rate; + float flame_smoke = *_flame_smoke; + float ignition_point = *_ignition_temp; + float temp_max = *_max_temp; + + for (int index = 0; index < total_cells; index++) { - size_t index = index_ +(z-1)*_slabSize; + float orig_fuel = fuel[index]; + float orig_smoke = smoke[index]; + float smoke_emit = 0.0f; + float react_coord = 0.0f; + + /* process fuel */ + fuel[index] -= burning_rate * dt; + if (fuel[index] < 0.0f) fuel[index] = 0.0f; + /* process reaction coordinate */ + if (orig_fuel) { + react[index] *= fuel[index]/orig_fuel; + react_coord = react[index]; + } + else { + react[index] = 0.0f; + } - for (int y = 1; y < _yRes - 1; y++, index += 2) - { - for (int x = 1; x < _xRes - 1; x++, index++) - { - // clean custom velocities from moving obstacles again - if (_obstacles[index]) - { - _xVelocity[index] = - _yVelocity[index] = - _zVelocity[index] = 0.0f; - } - } + /* emit smoke based on fuel burn rate and "flame_smoke" factor */ + smoke_emit = (orig_fuel < 1.0f) ? (1.0f - orig_fuel)*0.5f : 0.0f; + smoke_emit = (smoke_emit + 0.5f) * (orig_fuel-fuel[index]) * 0.1f * flame_smoke; + smoke[index] += smoke_emit; + CLAMP(smoke[index], 0.0f, 1.0f); + + /* model flame temperature curve from the reaction coordinate (fuel) */ + if (react_coord>0.0f) { + /* do a smooth falloff for rest of the values */ + flame[index] = pow(react_coord, 0.5f); + } + else + flame[index] = 0.0f; + + /* set fluid temperature from the flame temperature profile */ + if (heat && flame[index]) + heat[index] = (1.0f-flame[index])*ignition_point + flame[index]*temp_max; + + /* mix new color */ + if (r && smoke_emit) { + float smoke_factor = smoke[index]/(orig_smoke+smoke_emit); + r[index] = (r[index] + _flame_smoke_color[0] * smoke_emit) * smoke_factor; + g[index] = (g[index] + _flame_smoke_color[1] * smoke_emit) * smoke_factor; + b[index] = (b[index] + _flame_smoke_color[2] * smoke_emit) * smoke_factor; } } - } -#endif - - /*int begin=zBegin * _slabSize; - int end=begin + (zEnd - zBegin) * _slabSize; - for (int x = begin; x < end; x++) - _xForce[x] = _yForce[x] = 0.0f;*/ } diff --git a/intern/smoke/intern/FLUID_3D.h b/intern/smoke/intern/FLUID_3D.h index ac77148ce84..8cadf3bc989 100644 --- a/intern/smoke/intern/FLUID_3D.h +++ b/intern/smoke/intern/FLUID_3D.h @@ -46,11 +46,16 @@ class WTURBULENCE; class FLUID_3D { public: - FLUID_3D(int *res, /* int amplify, */ float *p0, float dtdef); + FLUID_3D(int *res, float dx, float dtdef, int init_heat, int init_fire, int init_colors); FLUID_3D() {}; virtual ~FLUID_3D(); - void initBlenderRNA(float *alpha, float *beta, float *dt_factor, float *vorticity, int *border_colli); + void initHeat(); + void initFire(); + void initColors(float init_r, float init_g, float init_b); + + void initBlenderRNA(float *alpha, float *beta, float *dt_factor, float *vorticity, int *border_colli, float *burning_rate, + float *flame_smoke, float *flame_smoke_color, float *flame_vorticity, float *ignition_temp, float *max_temp); // create & allocate vector noise advection void initVectorNoise(int amplify); @@ -58,7 +63,7 @@ class FLUID_3D void addSmokeColumn(); static void addSmokeTestCase(float* field, Vec3Int res); - void step(float dt); + void step(float dt, float gravity[3]); void addObstacle(OBSTACLE* obstacle); const float* xVelocity() { return _xVelocity; }; @@ -115,6 +120,27 @@ class FLUID_3D float* _heatTemp; float* _densityTemp; + // fire simulation + float *_flame; + float *_fuel; + float *_fuelTemp; + float *_fuelOld; + float *_react; + float *_reactTemp; + float *_reactOld; + + // smoke color + float *_color_r; + float *_color_rOld; + float *_color_rTemp; + float *_color_g; + float *_color_gOld; + float *_color_gTemp; + float *_color_b; + float *_color_bOld; + float *_color_bTemp; + + // CG fields int _iterations; @@ -153,14 +179,16 @@ class FLUID_3D void wipeBoundariesSL(int zBegin, int zEnd); void addForce(int zBegin, int zEnd); void addVorticity(int zBegin, int zEnd); - void addBuoyancy(float *heat, float *density, int zBegin, int zEnd); + void addBuoyancy(float *heat, float *density, float gravity[3], int zBegin, int zEnd); // solver stuff void project(); void diffuseHeat(); + void diffuseColor(); void solvePressure(float* field, float* b, unsigned char* skip); void solvePressurePre(float* field, float* b, unsigned char* skip); void solveHeat(float* field, float* b, unsigned char* skip); + void solveDiffusion(float* field, float* b, float* factor); // handle obstacle boundaries @@ -174,6 +202,16 @@ class FLUID_3D void advectMacCormackEnd1(int zBegin, int zEnd); void advectMacCormackEnd2(int zBegin, int zEnd); + /* burning */ + float *_burning_rate; // RNA pointer + float *_flame_smoke; // RNA pointer + float *_flame_smoke_color; // RNA pointer + float *_flame_vorticity; // RNA pointer + float *_ignition_temp; // RNA pointer + float *_max_temp; // RNA pointer + void processBurn(float *fuel, float *smoke, float *react, float *flame, float *heat, + float *r, float *g, float *b, int total_cells, float dt); + // boundary setting functions static void copyBorderX(float* field, Vec3Int res, int zBegin, int zEnd); static void copyBorderY(float* field, Vec3Int res, int zBegin, int zEnd); diff --git a/intern/smoke/intern/FLUID_3D_SOLVERS.cpp b/intern/smoke/intern/FLUID_3D_SOLVERS.cpp index 3cf94eb00a9..42a8b2d6923 100644 --- a/intern/smoke/intern/FLUID_3D_SOLVERS.cpp +++ b/intern/smoke/intern/FLUID_3D_SOLVERS.cpp @@ -165,7 +165,6 @@ void FLUID_3D::solveHeat(float* field, float* b, unsigned char* skip) if (_Acenter) delete[] _Acenter; } - void FLUID_3D::solvePressurePre(float* field, float* b, unsigned char* skip) { int x, y, z; diff --git a/intern/smoke/intern/FLUID_3D_STATIC.cpp b/intern/smoke/intern/FLUID_3D_STATIC.cpp index c7a0ddcd04c..ac485ad983a 100644 --- a/intern/smoke/intern/FLUID_3D_STATIC.cpp +++ b/intern/smoke/intern/FLUID_3D_STATIC.cpp @@ -92,18 +92,10 @@ void FLUID_3D::setNeumannX(float* field, Vec3Int res, int zBegin, int zEnd) // left slab index = y * res[0] + z * slabSize; field[index] = field[index + 2]; - /* only allow outwards flux */ - if(field[index]>0.) field[index] = 0.; - index += 1; - if(field[index]>0.) field[index] = 0.; // right slab index = y * res[0] + z * slabSize + res[0] - 1; field[index] = field[index - 2]; - /* only allow outwards flux */ - if(field[index]<0.) field[index] = 0.; - index -= 1; - if(field[index]<0.) field[index] = 0.; } } @@ -120,18 +112,10 @@ void FLUID_3D::setNeumannY(float* field, Vec3Int res, int zBegin, int zEnd) // front slab index = x + z * slabSize; field[index] = field[index + 2 * res[0]]; - /* only allow outwards flux */ - if(field[index]>0.) field[index] = 0.; - index += res[0]; - if(field[index]>0.) field[index] = 0.; // back slab index = x + z * slabSize + slabSize - res[0]; field[index] = field[index - 2 * res[0]]; - /* only allow outwards flux */ - if(field[index]<0.) field[index] = 0.; - index -= res[0]; - if(field[index]<0.) field[index] = 0.; } } @@ -152,14 +136,6 @@ void FLUID_3D::setNeumannZ(float* field, Vec3Int res, int zBegin, int zEnd) // front slab index = x + y * res[0]; field[index] = field[index + 2 * slabSize]; - /* only allow outwards flux */ - - // DG: Disable this for z-axis. - // The problem is that smoke somehow gets sucked in again - // from the TOP slab when this is enabled - // if(field[index]>0.) field[index] = 0.; - // index += slabSize; - // if(field[index]>0.) field[index] = 0.; } } @@ -170,10 +146,6 @@ void FLUID_3D::setNeumannZ(float* field, Vec3Int res, int zBegin, int zEnd) // back slab index = x + y * res[0] + cellsslab; field[index] = field[index - 2 * slabSize]; - /* only allow outwards flux */ - if(field[index]<0.) field[index] = 0.; - index -= slabSize; - if(field[index]<0.) field[index] = 0.; } } diff --git a/intern/smoke/intern/WTURBULENCE.cpp b/intern/smoke/intern/WTURBULENCE.cpp index 671198065e8..1c89f5d681b 100644 --- a/intern/smoke/intern/WTURBULENCE.cpp +++ b/intern/smoke/intern/WTURBULENCE.cpp @@ -51,7 +51,7 @@ static const float persistence = 0.56123f; ////////////////////////////////////////////////////////////////////// // constructor ////////////////////////////////////////////////////////////////////// -WTURBULENCE::WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify, int noisetype) +WTURBULENCE::WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify, int noisetype, int init_fire, int init_colors) { // if noise magnitude is below this threshold, its contribution // is negilgible, so stop evaluating new octaves @@ -87,12 +87,26 @@ WTURBULENCE::WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify, int no // allocate high resolution density field _totalStepsBig = 0; _densityBig = new float[_totalCellsBig]; - _densityBigOld = new float[_totalCellsBig]; + _densityBigOld = new float[_totalCellsBig]; for(int i = 0; i < _totalCellsBig; i++) { _densityBig[i] = _densityBigOld[i] = 0.; } + + /* fire */ + _flameBig = _fuelBig = _fuelBigOld = NULL; + _reactBig = _reactBigOld = NULL; + if (init_fire) { + initFire(); + } + /* colors */ + _color_rBig = _color_rBigOld = NULL; + _color_gBig = _color_gBigOld = NULL; + _color_bBig = _color_bBigOld = NULL; + if (init_colors) { + initColors(0.0f, 0.0f, 0.0f); + } // allocate & init texture coordinates _tcU = new float[_totalCellsSm]; @@ -128,12 +142,64 @@ WTURBULENCE::WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify, int no */ } +void WTURBULENCE::initFire() +{ + if (!_fuelBig) { + _flameBig = new float[_totalCellsBig]; + _fuelBig = new float[_totalCellsBig]; + _fuelBigOld = new float[_totalCellsBig]; + _reactBig = new float[_totalCellsBig]; + _reactBigOld = new float[_totalCellsBig]; + + for(int i = 0; i < _totalCellsBig; i++) { + _flameBig[i] = + _fuelBig[i] = + _fuelBigOld[i] = 0.; + _reactBig[i] = + _reactBigOld[i] = 0.; + } + } +} + +void WTURBULENCE::initColors(float init_r, float init_g, float init_b) +{ + if (!_color_rBig) { + _color_rBig = new float[_totalCellsBig]; + _color_rBigOld = new float[_totalCellsBig]; + _color_gBig = new float[_totalCellsBig]; + _color_gBigOld = new float[_totalCellsBig]; + _color_bBig = new float[_totalCellsBig]; + _color_bBigOld = new float[_totalCellsBig]; + + for(int i = 0; i < _totalCellsBig; i++) { + _color_rBig[i] = _densityBig[i] * init_r; + _color_rBigOld[i] = 0.0f; + _color_gBig[i] = _densityBig[i] * init_g; + _color_gBigOld[i] = 0.0f; + _color_bBig[i] = _densityBig[i] * init_b; + _color_bBigOld[i] = 0.0f; + } + } +} + ////////////////////////////////////////////////////////////////////// // destructor ////////////////////////////////////////////////////////////////////// WTURBULENCE::~WTURBULENCE() { delete[] _densityBig; delete[] _densityBigOld; + if (_flameBig) delete[] _flameBig; + if (_fuelBig) delete[] _fuelBig; + if (_fuelBigOld) delete[] _fuelBigOld; + if (_reactBig) delete[] _reactBig; + if (_reactBigOld) delete[] _reactBigOld; + + if (_color_rBig) delete[] _color_rBig; + if (_color_rBigOld) delete[] _color_rBigOld; + if (_color_gBig) delete[] _color_gBig; + if (_color_gBigOld) delete[] _color_gBigOld; + if (_color_bBig) delete[] _color_bBig; + if (_color_bBigOld) delete[] _color_bBigOld; delete[] _tcU; delete[] _tcV; @@ -757,8 +823,10 @@ void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, floa // enlarge timestep to match grid const float dt = dtOrg * _amplify; const float invAmp = 1.0f / _amplify; - float *tempBig1 = (float *)calloc(_totalCellsBig, sizeof(float)); - float *tempBig2 = (float *)calloc(_totalCellsBig, sizeof(float)); + float *tempFuelBig = NULL, *tempReactBig = NULL; + float *tempColor_rBig = NULL, *tempColor_gBig = NULL, *tempColor_bBig = NULL; + float *tempDensityBig = (float *)calloc(_totalCellsBig, sizeof(float)); + float *tempBig = (float *)calloc(_totalCellsBig, sizeof(float)); float *bigUx = (float *)calloc(_totalCellsBig, sizeof(float)); float *bigUy = (float *)calloc(_totalCellsBig, sizeof(float)); float *bigUz = (float *)calloc(_totalCellsBig, sizeof(float)); @@ -767,11 +835,21 @@ void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, floa float *eigMin = (float *)calloc(_totalCellsSm, sizeof(float)); float *eigMax = (float *)calloc(_totalCellsSm, sizeof(float)); + if (_fuelBig) { + tempFuelBig = (float *)calloc(_totalCellsBig, sizeof(float)); + tempReactBig = (float *)calloc(_totalCellsBig, sizeof(float)); + } + if (_color_rBig) { + tempColor_rBig = (float *)calloc(_totalCellsBig, sizeof(float)); + tempColor_gBig = (float *)calloc(_totalCellsBig, sizeof(float)); + tempColor_bBig = (float *)calloc(_totalCellsBig, sizeof(float)); + } + memset(_tcTemp, 0, sizeof(float)*_totalCellsSm); // prepare textures - advectTextureCoordinates(dtOrg, xvel,yvel,zvel, tempBig1, tempBig2); + advectTextureCoordinates(dtOrg, xvel,yvel,zvel, tempDensityBig, tempBig); // do wavelet decomposition of energy computeEnergy(_energy, xvel, yvel, zvel, obstacles); @@ -972,6 +1050,11 @@ void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, floa // prepare density for an advection SWAP_POINTERS(_densityBig, _densityBigOld); + SWAP_POINTERS(_fuelBig, _fuelBigOld); + SWAP_POINTERS(_reactBig, _reactBigOld); + SWAP_POINTERS(_color_rBig, _color_rBigOld); + SWAP_POINTERS(_color_gBig, _color_gBigOld); + SWAP_POINTERS(_color_bBig, _color_bBigOld); // based on the maximum velocity present, see if we need to substep, // but cap the maximum number of substeps to 5 @@ -1017,7 +1100,21 @@ void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, floa int zEnd = (int)((float)(i+1)*partSize + 0.5f); #endif FLUID_3D::advectFieldMacCormack1(dtSubdiv, bigUx, bigUy, bigUz, - _densityBigOld, tempBig1, _resBig, zBegin, zEnd); + _densityBigOld, tempDensityBig, _resBig, zBegin, zEnd); + if (_fuelBig) { + FLUID_3D::advectFieldMacCormack1(dtSubdiv, bigUx, bigUy, bigUz, + _fuelBigOld, tempFuelBig, _resBig, zBegin, zEnd); + FLUID_3D::advectFieldMacCormack1(dtSubdiv, bigUx, bigUy, bigUz, + _reactBigOld, tempReactBig, _resBig, zBegin, zEnd); + } + if (_color_rBig) { + FLUID_3D::advectFieldMacCormack1(dtSubdiv, bigUx, bigUy, bigUz, + _color_rBigOld, tempColor_rBig, _resBig, zBegin, zEnd); + FLUID_3D::advectFieldMacCormack1(dtSubdiv, bigUx, bigUy, bigUz, + _color_gBigOld, tempColor_gBig, _resBig, zBegin, zEnd); + FLUID_3D::advectFieldMacCormack1(dtSubdiv, bigUx, bigUy, bigUz, + _color_bBigOld, tempColor_bBig, _resBig, zBegin, zEnd); + } #if PARALLEL==1 } @@ -1030,18 +1127,43 @@ void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, floa int zEnd = (int)((float)(i+1)*partSize + 0.5f); #endif FLUID_3D::advectFieldMacCormack2(dtSubdiv, bigUx, bigUy, bigUz, - _densityBigOld, _densityBig, tempBig1, tempBig2, _resBig, NULL, zBegin, zEnd); + _densityBigOld, _densityBig, tempDensityBig, tempBig, _resBig, NULL, zBegin, zEnd); + if (_fuelBig) { + FLUID_3D::advectFieldMacCormack2(dtSubdiv, bigUx, bigUy, bigUz, + _fuelBigOld, _fuelBig, tempFuelBig, tempBig, _resBig, NULL, zBegin, zEnd); + FLUID_3D::advectFieldMacCormack2(dtSubdiv, bigUx, bigUy, bigUz, + _reactBigOld, _reactBig, tempReactBig, tempBig, _resBig, NULL, zBegin, zEnd); + } + if (_color_rBig) { + FLUID_3D::advectFieldMacCormack2(dtSubdiv, bigUx, bigUy, bigUz, + _color_rBigOld, _color_rBig, tempColor_rBig, tempBig, _resBig, NULL, zBegin, zEnd); + FLUID_3D::advectFieldMacCormack2(dtSubdiv, bigUx, bigUy, bigUz, + _color_gBigOld, _color_gBig, tempColor_gBig, tempBig, _resBig, NULL, zBegin, zEnd); + FLUID_3D::advectFieldMacCormack2(dtSubdiv, bigUx, bigUy, bigUz, + _color_bBigOld, _color_bBig, tempColor_bBig, tempBig, _resBig, NULL, zBegin, zEnd); + } #if PARALLEL==1 } } #endif - if (substep < totalSubsteps - 1) + if (substep < totalSubsteps - 1) { SWAP_POINTERS(_densityBig, _densityBigOld); + SWAP_POINTERS(_fuelBig, _fuelBigOld); + SWAP_POINTERS(_reactBig, _reactBigOld); + SWAP_POINTERS(_color_rBig, _color_rBigOld); + SWAP_POINTERS(_color_gBig, _color_gBigOld); + SWAP_POINTERS(_color_bBig, _color_bBigOld); + } } // substep - free(tempBig1); - free(tempBig2); + free(tempDensityBig); + if (tempFuelBig) free(tempFuelBig); + if (tempReactBig) free(tempReactBig); + if (tempColor_rBig) free(tempColor_rBig); + if (tempColor_gBig) free(tempColor_gBig); + if (tempColor_bBig) free(tempColor_bBig); + free(tempBig); free(bigUx); free(bigUy); free(bigUz); @@ -1050,6 +1172,15 @@ void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, floa // wipe the density borders FLUID_3D::setZeroBorder(_densityBig, _resBig, 0 , _resBig[2]); + if (_fuelBig) { + FLUID_3D::setZeroBorder(_fuelBig, _resBig, 0 , _resBig[2]); + FLUID_3D::setZeroBorder(_reactBig, _resBig, 0 , _resBig[2]); + } + if (_color_rBig) { + FLUID_3D::setZeroBorder(_color_rBig, _resBig, 0 , _resBig[2]); + FLUID_3D::setZeroBorder(_color_gBig, _resBig, 0 , _resBig[2]); + FLUID_3D::setZeroBorder(_color_bBig, _resBig, 0 , _resBig[2]); + } // reset texture coordinates now in preparation for next timestep // Shouldn't do this before generating the noise because then the diff --git a/intern/smoke/intern/WTURBULENCE.h b/intern/smoke/intern/WTURBULENCE.h index f31ca100fdf..1655bd95d32 100644 --- a/intern/smoke/intern/WTURBULENCE.h +++ b/intern/smoke/intern/WTURBULENCE.h @@ -36,10 +36,13 @@ class WTURBULENCE { public: // both config files can be NULL, altCfg might override values from noiseCfg - WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify, int noisetype); + WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify, int noisetype, int init_fire, int init_colors); /// destructor virtual ~WTURBULENCE(); + + void initFire(); + void initColors(float init_r, float init_g, float init_b); void setNoise(int type); void initBlenderRNA(float *strength); @@ -63,6 +66,8 @@ class WTURBULENCE // access functions inline float* getDensityBig() { return _densityBig; } + inline float* getFlameBig() { return _flameBig; } + inline float* getFuelBig() { return _fuelBig; } inline float* getArrayTcU() { return _tcU; } inline float* getArrayTcV() { return _tcV; } inline float* getArrayTcW() { return _tcW; } @@ -111,6 +116,18 @@ class WTURBULENCE float* _densityBig; float* _densityBigOld; + float* _flameBig; + float* _fuelBig; + float* _fuelBigOld; + float* _reactBig; + float* _reactBigOld; + + float* _color_rBig; + float* _color_rBigOld; + float* _color_gBig; + float* _color_gBigOld; + float* _color_bBig; + float* _color_bBigOld; // texture coordinates for noise float* _tcU; diff --git a/intern/smoke/intern/smoke_API.cpp b/intern/smoke/intern/smoke_API.cpp index 4bbf8e0a82b..e51c3176699 100644 --- a/intern/smoke/intern/smoke_API.cpp +++ b/intern/smoke/intern/smoke_API.cpp @@ -30,6 +30,7 @@ #include "FLUID_3D.h" #include "WTURBULENCE.h" +#include "spectrum.h" #include #include @@ -37,22 +38,16 @@ #include "../extern/smoke_API.h" /* to ensure valid prototypes */ -// y in smoke is z in blender -extern "C" FLUID_3D *smoke_init(int *res, float *p0, float dtdef) +extern "C" FLUID_3D *smoke_init(int *res, float dx, float dtdef, int use_heat, int use_fire, int use_colors) { - // smoke lib uses y as top-bottom/vertical axis where blender uses z - FLUID_3D *fluid = new FLUID_3D(res, p0, dtdef); - - // printf("xres: %d, yres: %d, zres: %d\n", res[0], res[1], res[2]); - + FLUID_3D *fluid = new FLUID_3D(res, dx, dtdef, use_heat, use_fire, use_colors); return fluid; } -extern "C" WTURBULENCE *smoke_turbulence_init(int *res, int amplify, int noisetype) +extern "C" WTURBULENCE *smoke_turbulence_init(int *res, int amplify, int noisetype, int use_fire, int use_colors) { - // initialize wavelet turbulence if(amplify) - return new WTURBULENCE(res[0],res[1],res[2], amplify, noisetype); + return new WTURBULENCE(res[0],res[1],res[2], amplify, noisetype, use_fire, use_colors); else return NULL; } @@ -71,7 +66,6 @@ extern "C" void smoke_turbulence_free(WTURBULENCE *wt) extern "C" size_t smoke_get_index(int x, int max_x, int y, int max_y, int z /*, int max_z */) { - // // const int index = x + y * smd->res[0] + z * smd->res[0]*smd->res[1]; return x + y * max_x + z * max_x*max_y; } @@ -80,137 +74,134 @@ extern "C" size_t smoke_get_index2d(int x, int max_x, int y /*, int max_y, int z return x + y * max_x; } -extern "C" void smoke_step(FLUID_3D *fluid, float dtSubdiv) +extern "C" void smoke_step(FLUID_3D *fluid, float gravity[3], float dtSubdiv) { - fluid->step(dtSubdiv); + if (fluid->_fuel) { + fluid->processBurn(fluid->_fuel, fluid->_density, fluid->_react, fluid->_flame, fluid->_heat, + fluid->_color_r, fluid->_color_g, fluid->_color_b, fluid->_totalCells, (*fluid->_dtFactor)*dtSubdiv); + } + fluid->step(dtSubdiv, gravity); } extern "C" void smoke_turbulence_step(WTURBULENCE *wt, FLUID_3D *fluid) { + if (wt->_fuelBig) { + fluid->processBurn(wt->_fuelBig, wt->_densityBig, wt->_reactBig, wt->_flameBig, 0, + wt->_color_rBig, wt->_color_gBig, wt->_color_bBig, wt->_totalCellsBig, fluid->_dt); + } wt->stepTurbulenceFull(fluid->_dt/fluid->_dx, fluid->_xVelocity, fluid->_yVelocity, fluid->_zVelocity, fluid->_obstacles); } -extern "C" void smoke_initBlenderRNA(FLUID_3D *fluid, float *alpha, float *beta, float *dt_factor, float *vorticity, int *border_colli) +extern "C" void smoke_initBlenderRNA(FLUID_3D *fluid, float *alpha, float *beta, float *dt_factor, float *vorticity, int *border_colli, float *burning_rate, + float *flame_smoke, float *flame_smoke_color, float *flame_vorticity, float *flame_ignition_temp, float *flame_max_temp) { - fluid->initBlenderRNA(alpha, beta, dt_factor, vorticity, border_colli); + fluid->initBlenderRNA(alpha, beta, dt_factor, vorticity, border_colli, burning_rate, flame_smoke, flame_smoke_color, flame_vorticity, flame_ignition_temp, flame_max_temp); } -extern "C" void smoke_dissolve(FLUID_3D *fluid, int speed, int log) +extern "C" void smoke_initWaveletBlenderRNA(WTURBULENCE *wt, float *strength) { - float *density = fluid->_density; - //float *densityOld = fluid->_densityOld; - float *heat = fluid->_heat; + wt->initBlenderRNA(strength); +} +static void data_dissolve(float *density, float *heat, float *r, float *g, float *b, int total_cells, int speed, int log) +{ if(log) { /* max density/speed = dydx */ - float dydx = 1.0 / (float)speed; - size_t size= fluid->_xRes * fluid->_yRes * fluid->_zRes; + float fac = 1.0f - (1.0f / (float)speed); - for(size_t i = 0; i < size; i++) + for(size_t i = 0; i < total_cells; i++) { - density[i] *= (1.0 - dydx); - - if(density[i] < 0.0f) - density[i] = 0.0f; - - heat[i] *= (1.0 - dydx); - - /*if(heat[i] < 0.0f) - heat[i] = 0.0f;*/ + /* density */ + density[i] *= fac; + + /* heat */ + if (heat) { + heat[i] *= fac; + } + + /* color */ + if (r) { + r[i] *= fac; + g[i] *= fac; + b[i] *= fac; + } } } else // linear falloff { /* max density/speed = dydx */ - float dydx = 1.0 / (float)speed; - size_t size= fluid->_xRes * fluid->_yRes * fluid->_zRes; + float dydx = 1.0f / (float)speed; - for(size_t i = 0; i < size; i++) + for(size_t i = 0; i < total_cells; i++) { + float d = density[i]; + /* density */ density[i] -= dydx; - if(density[i] < 0.0f) density[i] = 0.0f; - if(abs(heat[i]) < dydx) heat[i] = 0.0f; - else if (heat[i]>0.0f) heat[i] -= dydx; - else if (heat[i]<0.0f) heat[i] += dydx; + /* heat */ + if (heat) { + if(abs(heat[i]) < dydx) heat[i] = 0.0f; + else if (heat[i]>0.0f) heat[i] -= dydx; + else if (heat[i]<0.0f) heat[i] += dydx; + } + + /* color */ + if (r && d) { + r[i] *= (density[i]/d); + g[i] *= (density[i]/d); + b[i] *= (density[i]/d); + } } } } -extern "C" void smoke_dissolve_wavelet(WTURBULENCE *wt, int speed, int log) -{ - float *density = wt->getDensityBig(); - Vec3Int r = wt->getResBig(); - - if(log) - { - /* max density/speed = dydx */ - float dydx = 1.0 / (float)speed; - size_t size= r[0] * r[1] * r[2]; - - for(size_t i = 0; i < size; i++) - { - density[i] *= (1.0 - dydx); - - if(density[i] < 0.0f) - density[i] = 0.0f; - } - } - else // linear falloff - { - /* max density/speed = dydx */ - float dydx = 1.0 / (float)speed; - size_t size= r[0] * r[1] * r[2]; - - for(size_t i = 0; i < size; i++) - { - density[i] -= dydx; - - if(density[i] < 0.0f) - density[i] = 0.0f; - } - } -} - -extern "C" void smoke_initWaveletBlenderRNA(WTURBULENCE *wt, float *strength) +extern "C" void smoke_dissolve(FLUID_3D *fluid, int speed, int log) { - wt->initBlenderRNA(strength); + data_dissolve(fluid->_density, fluid->_heat, fluid->_color_r, fluid->_color_g, fluid->_color_b, fluid->_totalCells, speed, log); } -template < class T > inline T ABS( T a ) +extern "C" void smoke_dissolve_wavelet(WTURBULENCE *wt, int speed, int log) { - return (0 < a) ? a : -a ; + data_dissolve(wt->_densityBig, 0, wt->_color_rBig, wt->_color_gBig, wt->_color_bBig, wt->_totalCellsBig, speed, log); } -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) +extern "C" void smoke_export(FLUID_3D *fluid, float *dt, float *dx, float **dens, float **react, float **flame, float **fuel, float **heat, + float **heatold, float **vx, float **vy, float **vz, float **r, float **g, float **b, unsigned char **obstacles) { *dens = fluid->_density; - *densold = fluid->_densityOld; + *fuel = fluid->_fuel; + *react = fluid->_react; + *flame = fluid->_flame; *heat = fluid->_heat; *heatold = fluid->_heatOld; *vx = fluid->_xVelocity; *vy = fluid->_yVelocity; *vz = fluid->_zVelocity; - *vxold = fluid->_xVelocityOld; - *vyold = fluid->_yVelocityOld; - *vzold = fluid->_zVelocityOld; + *r = fluid->_color_r; + *g = fluid->_color_g; + *b = fluid->_color_b; *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) +extern "C" void smoke_turbulence_export(WTURBULENCE *wt, float **dens, float **react, float **flame, float **fuel, + float **r, float **g, float **b , float **tcu, float **tcv, float **tcw) { if(!wt) return; *dens = wt->_densityBig; - *densold = wt->_densityBigOld; + *fuel = wt->_fuelBig; + *react = wt->_reactBig; + *flame = wt->_flameBig; + *r = wt->_color_rBig; + *g = wt->_color_gBig; + *b = wt->_color_bBig; *tcu = wt->_tcU; *tcv = wt->_tcV; *tcw = wt->_tcW; @@ -221,6 +212,16 @@ extern "C" float *smoke_get_density(FLUID_3D *fluid) return fluid->_density; } +extern "C" float *smoke_get_fuel(FLUID_3D *fluid) +{ + return fluid->_fuel; +} + +extern "C" float *smoke_get_react(FLUID_3D *fluid) +{ + return fluid->_react; +} + extern "C" float *smoke_get_heat(FLUID_3D *fluid) { return fluid->_heat; @@ -256,15 +257,137 @@ extern "C" float *smoke_get_force_z(FLUID_3D *fluid) return fluid->_zForce; } +extern "C" float *smoke_get_flame(FLUID_3D *fluid) +{ + return fluid->_flame; +} + +extern "C" float *smoke_get_color_r(FLUID_3D *fluid) +{ + return fluid->_color_r; +} + +extern "C" float *smoke_get_color_g(FLUID_3D *fluid) +{ + return fluid->_color_g; +} + +extern "C" float *smoke_get_color_b(FLUID_3D *fluid) +{ + return fluid->_color_b; +} + +static void get_rgba(float *r, float *g, float *b, float *a, int total_cells, float *data, int sequential) +{ + int i; + int m = 4, i_g = 1, i_b = 2, i_a = 3; + /* sequential data */ + if (sequential) { + m = 1; + i_g *= total_cells; + i_b *= total_cells; + i_a *= total_cells; + } + + for (i=0; i_color_r, fluid->_color_g, fluid->_color_b, fluid->_density, fluid->_totalCells, data, sequential); +} + +extern "C" void smoke_turbulence_get_rgba(WTURBULENCE *wt, float *data, int sequential) +{ + get_rgba(wt->_color_rBig, wt->_color_gBig, wt->_color_bBig, wt->_densityBig, wt->_totalCellsBig, data, sequential); +} + +/* get a single color premultiplied voxel grid */ +static void get_rgba_from_density(float color[3], float *a, int total_cells, float *data, int sequential) +{ + int i; + int m = 4, i_g = 1, i_b = 2, i_a = 3; + /* sequential data */ + if (sequential) { + m = 1; + i_g *= total_cells; + i_b *= total_cells; + i_a *= total_cells; + } + + for (i=0; i_density, fluid->_totalCells, data, sequential); +} + +extern "C" void smoke_turbulence_get_rgba_from_density(WTURBULENCE *wt, float color[3], float *data, int sequential) +{ + get_rgba_from_density(color, wt->_densityBig, wt->_totalCellsBig, data, sequential); +} + extern "C" float *smoke_turbulence_get_density(WTURBULENCE *wt) { return wt ? wt->getDensityBig() : NULL; } +extern "C" float *smoke_turbulence_get_fuel(WTURBULENCE *wt) +{ + return wt ? wt->getFuelBig() : NULL; +} + +extern "C" float *smoke_turbulence_get_react(WTURBULENCE *wt) +{ + return wt ? wt->_reactBig : NULL; +} + +extern "C" float *smoke_turbulence_get_color_r(WTURBULENCE *wt) +{ + return wt ? wt->_color_rBig : NULL; +} + +extern "C" float *smoke_turbulence_get_color_g(WTURBULENCE *wt) +{ + return wt ? wt->_color_gBig : NULL; +} + +extern "C" float *smoke_turbulence_get_color_b(WTURBULENCE *wt) +{ + return wt ? wt->_color_bBig : NULL; +} + +extern "C" float *smoke_turbulence_get_flame(WTURBULENCE *wt) +{ + return wt ? wt->getFlameBig() : NULL; +} + extern "C" void smoke_turbulence_get_res(WTURBULENCE *wt, int *res) { - if(wt) - { + if(wt) { Vec3Int r = wt->getResBig(); res[0] = r[0]; res[1] = r[1]; @@ -272,6 +395,15 @@ extern "C" void smoke_turbulence_get_res(WTURBULENCE *wt, int *res) } } +extern "C" int smoke_turbulence_get_cells(WTURBULENCE *wt) +{ + if(wt) { + Vec3Int r = wt->getResBig(); + return r[0]*r[1]*r[2]; + } + return 0; +} + extern "C" unsigned char *smoke_get_obstacle(FLUID_3D *fluid) { return fluid->_obstacles; @@ -295,3 +427,61 @@ extern "C" void smoke_turbulence_set_noise(WTURBULENCE *wt, int type) { wt->setNoise(type); } + +extern "C" void flame_get_spectrum(unsigned char *spec, int width, float t1, float t2) +{ + spectrum(t1, t2, width, spec); +} + +extern "C" int smoke_has_heat(FLUID_3D *fluid) +{ + return (fluid->_heat) ? 1 : 0; +} + +extern "C" int smoke_has_fuel(FLUID_3D *fluid) +{ + return (fluid->_fuel) ? 1 : 0; +} + +extern "C" int smoke_has_colors(FLUID_3D *fluid) +{ + return (fluid->_color_r && fluid->_color_g && fluid->_color_b) ? 1 : 0; +} + +extern "C" int smoke_turbulence_has_fuel(WTURBULENCE *wt) +{ + return (wt->_fuelBig) ? 1 : 0; +} + +extern "C" int smoke_turbulence_has_colors(WTURBULENCE *wt) +{ + return (wt->_color_rBig && wt->_color_gBig && wt->_color_bBig) ? 1 : 0; +} + +/* additional field initialization */ +extern "C" void smoke_ensure_heat(FLUID_3D *fluid) +{ + if (fluid) { + fluid->initHeat(); + } +} + +extern "C" void smoke_ensure_fire(FLUID_3D *fluid, WTURBULENCE *wt) +{ + if (fluid) { + fluid->initFire(); + } + if (wt) { + wt->initFire(); + } +} + +extern "C" void smoke_ensure_colors(FLUID_3D *fluid, WTURBULENCE *wt, float init_r, float init_g, float init_b) +{ + if (fluid) { + fluid->initColors(init_r, init_g, init_b); + } + if (wt) { + wt->initColors(init_r, init_g, init_b); + } +} \ No newline at end of file diff --git a/intern/smoke/intern/spectrum.cpp b/intern/smoke/intern/spectrum.cpp new file mode 100644 index 00000000000..34725c12c92 --- /dev/null +++ b/intern/smoke/intern/spectrum.cpp @@ -0,0 +1,411 @@ +/* + Colour Rendering of Spectra + + by John Walker + http://www.fourmilab.ch/ + + Last updated: March 9, 2003 + + This program is in the public domain. + + For complete information about the techniques employed in + this program, see the World-Wide Web document: + + http://www.fourmilab.ch/documents/specrend/ + + The xyz_to_rgb() function, which was wrong in the original + version of this program, was corrected by: + + Andrew J. S. Hamilton 21 May 1999 + Andrew.Hamilton@Colorado.EDU + http://casa.colorado.edu/~ajsh/ + + who also added the gamma correction facilities and + modified constrain_rgb() to work by desaturating the + colour by adding white. + + A program which uses these functions to plot CIE + "tongue" diagrams called "ppmcie" is included in + the Netpbm graphics toolkit: + http://netpbm.sourceforge.net/ + (The program was called cietoppm in earlier + versions of Netpbm.) + +*/ + +#include +#include + +/* A colour system is defined by the CIE x and y coordinates of + its three primary illuminants and the x and y coordinates of + the white point. */ + +struct colourSystem { + char *name; /* Colour system name */ + double xRed, yRed, /* Red x, y */ + xGreen, yGreen, /* Green x, y */ + xBlue, yBlue, /* Blue x, y */ + xWhite, yWhite, /* White point x, y */ + gamma; /* Gamma correction for system */ +}; + +/* White point chromaticities. */ + +#define IlluminantC 0.3101, 0.3162 /* For NTSC television */ +#define IlluminantD65 0.3127, 0.3291 /* For EBU and SMPTE */ +#define IlluminantE 0.33333333, 0.33333333 /* CIE equal-energy illuminant */ + +/* Gamma of nonlinear correction. + + See Charles Poynton's ColorFAQ Item 45 and GammaFAQ Item 6 at: + + http://www.poynton.com/ColorFAQ.html + http://www.poynton.com/GammaFAQ.html + +*/ + +#define GAMMA_REC709 0 /* Rec. 709 */ + +static struct colourSystem + /* Name xRed yRed xGreen yGreen xBlue yBlue White point Gamma */ + NTSCsystem = { "NTSC", 0.67, 0.33, 0.21, 0.71, 0.14, 0.08, IlluminantC, GAMMA_REC709 }, + EBUsystem = { "EBU (PAL/SECAM)", 0.64, 0.33, 0.29, 0.60, 0.15, 0.06, IlluminantD65, GAMMA_REC709 }, + SMPTEsystem = { "SMPTE", 0.630, 0.340, 0.310, 0.595, 0.155, 0.070, IlluminantD65, GAMMA_REC709 }, + HDTVsystem = { "HDTV", 0.670, 0.330, 0.210, 0.710, 0.150, 0.060, IlluminantD65, GAMMA_REC709 }, + CIEsystem = { "CIE", 0.7355, 0.2645, 0.2658, 0.7243, 0.1669, 0.0085, IlluminantE, GAMMA_REC709 }, + Rec709system = { "CIE REC 709", 0.64, 0.33, 0.30, 0.60, 0.15, 0.06, IlluminantD65, GAMMA_REC709 }; + +/* UPVP_TO_XY + + Given 1976 coordinates u', v', determine 1931 chromaticities x, y + +*/ + +void upvp_to_xy(double up, double vp, double *xc, double *yc) +{ + *xc = (9 * up) / ((6 * up) - (16 * vp) + 12); + *yc = (4 * vp) / ((6 * up) - (16 * vp) + 12); +} + +/* XY_TO_UPVP + + Given 1931 chromaticities x, y, determine 1976 coordinates u', v' + +*/ + +void xy_to_upvp(double xc, double yc, double *up, double *vp) +{ + *up = (4 * xc) / ((-2 * xc) + (12 * yc) + 3); + *vp = (9 * yc) / ((-2 * xc) + (12 * yc) + 3); +} + +/* XYZ_TO_RGB + + Given an additive tricolour system CS, defined by the CIE x + and y chromaticities of its three primaries (z is derived + trivially as 1-(x+y)), and a desired chromaticity (XC, YC, + ZC) in CIE space, determine the contribution of each + primary in a linear combination which sums to the desired + chromaticity. If the requested chromaticity falls outside + the Maxwell triangle (colour gamut) formed by the three + primaries, one of the r, g, or b weights will be negative. + + Caller can use constrain_rgb() to desaturate an + outside-gamut colour to the closest representation within + the available gamut and/or norm_rgb to normalise the RGB + components so the largest nonzero component has value 1. + +*/ + +void xyz_to_rgb(struct colourSystem *cs, + double xc, double yc, double zc, + double *r, double *g, double *b) +{ + double xr, yr, zr, xg, yg, zg, xb, yb, zb; + double xw, yw, zw; + double rx, ry, rz, gx, gy, gz, bx, by, bz; + double rw, gw, bw; + + xr = cs->xRed; yr = cs->yRed; zr = 1 - (xr + yr); + xg = cs->xGreen; yg = cs->yGreen; zg = 1 - (xg + yg); + xb = cs->xBlue; yb = cs->yBlue; zb = 1 - (xb + yb); + + xw = cs->xWhite; yw = cs->yWhite; zw = 1 - (xw + yw); + + /* xyz -> rgb matrix, before scaling to white. */ + + rx = (yg * zb) - (yb * zg); ry = (xb * zg) - (xg * zb); rz = (xg * yb) - (xb * yg); + gx = (yb * zr) - (yr * zb); gy = (xr * zb) - (xb * zr); gz = (xb * yr) - (xr * yb); + bx = (yr * zg) - (yg * zr); by = (xg * zr) - (xr * zg); bz = (xr * yg) - (xg * yr); + + /* White scaling factors. + Dividing by yw scales the white luminance to unity, as conventional. */ + + rw = ((rx * xw) + (ry * yw) + (rz * zw)) / yw; + gw = ((gx * xw) + (gy * yw) + (gz * zw)) / yw; + bw = ((bx * xw) + (by * yw) + (bz * zw)) / yw; + + /* xyz -> rgb matrix, correctly scaled to white. */ + + rx = rx / rw; ry = ry / rw; rz = rz / rw; + gx = gx / gw; gy = gy / gw; gz = gz / gw; + bx = bx / bw; by = by / bw; bz = bz / bw; + + /* rgb of the desired point */ + + *r = (rx * xc) + (ry * yc) + (rz * zc); + *g = (gx * xc) + (gy * yc) + (gz * zc); + *b = (bx * xc) + (by * yc) + (bz * zc); +} + +/* INSIDE_GAMUT + + Test whether a requested colour is within the gamut + achievable with the primaries of the current colour + system. This amounts simply to testing whether all the + primary weights are non-negative. */ + +int inside_gamut(double r, double g, double b) +{ + return (r >= 0) && (g >= 0) && (b >= 0); +} + +/* CONSTRAIN_RGB + + If the requested RGB shade contains a negative weight for + one of the primaries, it lies outside the colour gamut + accessible from the given triple of primaries. Desaturate + it by adding white, equal quantities of R, G, and B, enough + to make RGB all positive. The function returns 1 if the + components were modified, zero otherwise. + +*/ + +int constrain_rgb(double *r, double *g, double *b) +{ + double w; + + /* Amount of white needed is w = - min(0, *r, *g, *b) */ + + w = (0 < *r) ? 0 : *r; + w = (w < *g) ? w : *g; + w = (w < *b) ? w : *b; + w = -w; + + /* Add just enough white to make r, g, b all positive. */ + + if (w > 0) { + *r += w; *g += w; *b += w; + return 1; /* Colour modified to fit RGB gamut */ + } + + return 0; /* Colour within RGB gamut */ +} + +/* GAMMA_CORRECT_RGB + + Transform linear RGB values to nonlinear RGB values. Rec. + 709 is ITU-R Recommendation BT. 709 (1990) ``Basic + Parameter Values for the HDTV Standard for the Studio and + for International Programme Exchange'', formerly CCIR Rec. + 709. For details see + + http://www.poynton.com/ColorFAQ.html + http://www.poynton.com/GammaFAQ.html +*/ + +void gamma_correct(const struct colourSystem *cs, double *c) +{ + double gamma; + + gamma = cs->gamma; + + if (gamma == GAMMA_REC709) { + /* Rec. 709 gamma correction. */ + double cc = 0.018; + + if (*c < cc) { + *c *= ((1.099 * pow(cc, 0.45)) - 0.099) / cc; + } else { + *c = (1.099 * pow(*c, 0.45)) - 0.099; + } + } else { + /* Nonlinear colour = (Linear colour)^(1/gamma) */ + *c = pow(*c, 1.0 / gamma); + } +} + +void gamma_correct_rgb(const struct colourSystem *cs, double *r, double *g, double *b) +{ + gamma_correct(cs, r); + gamma_correct(cs, g); + gamma_correct(cs, b); +} + +/* NORM_RGB + + Normalise RGB components so the most intense (unless all + are zero) has a value of 1. + +*/ + +void norm_rgb(double *r, double *g, double *b) +{ +#define Max(a, b) (((a) > (b)) ? (a) : (b)) + double greatest = Max(*r, Max(*g, *b)); + + if (greatest > 0) { + *r /= greatest; + *g /= greatest; + *b /= greatest; + } +#undef Max +} + +/* SPECTRUM_TO_XYZ + + Calculate the CIE X, Y, and Z coordinates corresponding to + a light source with spectral distribution given by the + function SPEC_INTENS, which is called with a series of + wavelengths between 380 and 780 nm (the argument is + expressed in meters), which returns emittance at that + wavelength in arbitrary units. The chromaticity + coordinates of the spectrum are returned in the x, y, and z + arguments which respect the identity: + + x + y + z = 1. +*/ + +void spectrum_to_xyz(double (*spec_intens)(double wavelength), + double *x, double *y, double *z) +{ + int i; + double lambda, X = 0, Y = 0, Z = 0, XYZ; + + /* CIE colour matching functions xBar, yBar, and zBar for + wavelengths from 380 through 780 nanometers, every 5 + nanometers. For a wavelength lambda in this range: + + cie_colour_match[(lambda - 380) / 5][0] = xBar + cie_colour_match[(lambda - 380) / 5][1] = yBar + cie_colour_match[(lambda - 380) / 5][2] = zBar + + To save memory, this table can be declared as floats + rather than doubles; (IEEE) float has enough + significant bits to represent the values. It's declared + as a double here to avoid warnings about "conversion + between floating-point types" from certain persnickety + compilers. */ + + static double cie_colour_match[81][3] = { + {0.0014,0.0000,0.0065}, {0.0022,0.0001,0.0105}, {0.0042,0.0001,0.0201}, + {0.0076,0.0002,0.0362}, {0.0143,0.0004,0.0679}, {0.0232,0.0006,0.1102}, + {0.0435,0.0012,0.2074}, {0.0776,0.0022,0.3713}, {0.1344,0.0040,0.6456}, + {0.2148,0.0073,1.0391}, {0.2839,0.0116,1.3856}, {0.3285,0.0168,1.6230}, + {0.3483,0.0230,1.7471}, {0.3481,0.0298,1.7826}, {0.3362,0.0380,1.7721}, + {0.3187,0.0480,1.7441}, {0.2908,0.0600,1.6692}, {0.2511,0.0739,1.5281}, + {0.1954,0.0910,1.2876}, {0.1421,0.1126,1.0419}, {0.0956,0.1390,0.8130}, + {0.0580,0.1693,0.6162}, {0.0320,0.2080,0.4652}, {0.0147,0.2586,0.3533}, + {0.0049,0.3230,0.2720}, {0.0024,0.4073,0.2123}, {0.0093,0.5030,0.1582}, + {0.0291,0.6082,0.1117}, {0.0633,0.7100,0.0782}, {0.1096,0.7932,0.0573}, + {0.1655,0.8620,0.0422}, {0.2257,0.9149,0.0298}, {0.2904,0.9540,0.0203}, + {0.3597,0.9803,0.0134}, {0.4334,0.9950,0.0087}, {0.5121,1.0000,0.0057}, + {0.5945,0.9950,0.0039}, {0.6784,0.9786,0.0027}, {0.7621,0.9520,0.0021}, + {0.8425,0.9154,0.0018}, {0.9163,0.8700,0.0017}, {0.9786,0.8163,0.0014}, + {1.0263,0.7570,0.0011}, {1.0567,0.6949,0.0010}, {1.0622,0.6310,0.0008}, + {1.0456,0.5668,0.0006}, {1.0026,0.5030,0.0003}, {0.9384,0.4412,0.0002}, + {0.8544,0.3810,0.0002}, {0.7514,0.3210,0.0001}, {0.6424,0.2650,0.0000}, + {0.5419,0.2170,0.0000}, {0.4479,0.1750,0.0000}, {0.3608,0.1382,0.0000}, + {0.2835,0.1070,0.0000}, {0.2187,0.0816,0.0000}, {0.1649,0.0610,0.0000}, + {0.1212,0.0446,0.0000}, {0.0874,0.0320,0.0000}, {0.0636,0.0232,0.0000}, + {0.0468,0.0170,0.0000}, {0.0329,0.0119,0.0000}, {0.0227,0.0082,0.0000}, + {0.0158,0.0057,0.0000}, {0.0114,0.0041,0.0000}, {0.0081,0.0029,0.0000}, + {0.0058,0.0021,0.0000}, {0.0041,0.0015,0.0000}, {0.0029,0.0010,0.0000}, + {0.0020,0.0007,0.0000}, {0.0014,0.0005,0.0000}, {0.0010,0.0004,0.0000}, + {0.0007,0.0002,0.0000}, {0.0005,0.0002,0.0000}, {0.0003,0.0001,0.0000}, + {0.0002,0.0001,0.0000}, {0.0002,0.0001,0.0000}, {0.0001,0.0000,0.0000}, + {0.0001,0.0000,0.0000}, {0.0001,0.0000,0.0000}, {0.0000,0.0000,0.0000} + }; + + for (i = 0, lambda = 380; lambda < 780.1; i++, lambda += 5) { + double Me; + + Me = (*spec_intens)(lambda); + X += Me * cie_colour_match[i][0]; + Y += Me * cie_colour_match[i][1]; + Z += Me * cie_colour_match[i][2]; + } + XYZ = (X + Y + Z); + *x = X / XYZ; + *y = Y / XYZ; + *z = Z / XYZ; +} + +/* BB_SPECTRUM + + Calculate, by Planck's radiation law, the emittance of a black body + of temperature bbTemp at the given wavelength (in metres). */ + +double bbTemp = 5000; /* Hidden temperature argument + to BB_SPECTRUM. */ +double bb_spectrum(double wavelength) +{ + double wlm = wavelength * 1e-9; /* Wavelength in meters */ + + return (3.74183e-16 * pow(wlm, -5.0)) / + (exp(1.4388e-2 / (wlm * bbTemp)) - 1.0); +} + +void xyz_to_lms(double x, double y, double z, double* l, double* m, double* s) +{ + *l = 0.3897*x + 0.6890*y - 0.0787*z; + *m = -0.2298*x + 1.1834*y + 0.0464*z; + *s = z; +} + +void lms_to_xyz(double l, double m, double s, double* x, double *y, double* z) +{ + *x = 1.9102*l - 1.1121*m + 0.2019*s; + *y = 0.3709*l + 0.6290*m + 0.0000*s; + *z = s; +} + +void spectrum(double t1, double t2, int N, unsigned char* d) +{ + int i,j,dj; + double X,Y,Z,R,G,B,L,M,S, Lw, Mw, Sw; + struct colourSystem *cs = &CIEsystem; + + j = 0; dj = 1; + if (t10.1)? B*255 : 0; + j += dj; + } +} diff --git a/intern/smoke/intern/spectrum.h b/intern/smoke/intern/spectrum.h new file mode 100644 index 00000000000..9edd9ad887c --- /dev/null +++ b/intern/smoke/intern/spectrum.h @@ -0,0 +1,6 @@ +#ifndef __SPECTRUM_H +#define __SPECTRUM_H + +void spectrum(double t1, double t2, int n, unsigned char* d); + +#endif -- cgit v1.2.3