diff options
Diffstat (limited to 'intern/smoke/intern')
| -rw-r--r-- | intern/smoke/intern/FLUID_3D.cpp | 228 | ||||
| -rw-r--r-- | intern/smoke/intern/FLUID_3D.h | 17 | ||||
| -rw-r--r-- | intern/smoke/intern/OBSTACLE.h | 6 | ||||
| -rw-r--r-- | intern/smoke/intern/WTURBULENCE.cpp | 9 | ||||
| -rw-r--r-- | intern/smoke/intern/smoke_API.cpp | 53 |
5 files changed, 258 insertions, 55 deletions
diff --git a/intern/smoke/intern/FLUID_3D.cpp b/intern/smoke/intern/FLUID_3D.cpp index 9f036cc6d2f..04971f898e9 100644 --- a/intern/smoke/intern/FLUID_3D.cpp +++ b/intern/smoke/intern/FLUID_3D.cpp @@ -34,6 +34,8 @@ #include "SPHERE.h" #include <zlib.h> +#include "float.h" + #if PARALLEL==1 #include <omp.h> #endif // PARALLEL @@ -42,11 +44,11 @@ // Construction/Destruction ////////////////////////////////////////////////////////////////////// -FLUID_3D::FLUID_3D(int *res, float *p0) : +FLUID_3D::FLUID_3D(int *res, float *p0, float dtdef) : _xRes(res[0]), _yRes(res[1]), _zRes(res[2]), _res(0.0f) { // set simulation consts - _dt = DT_DEFAULT; // just in case. set in step from a RNA factor + _dt = dtdef; // just in case. set in step from a RNA factor // start point of array _p0[0] = p0[0]; @@ -81,6 +83,9 @@ FLUID_3D::FLUID_3D(int *res, float *p0) : _xVelocity = new float[_totalCells]; _yVelocity = new float[_totalCells]; _zVelocity = new float[_totalCells]; + _xVelocityOb = new float[_totalCells]; + _yVelocityOb = new float[_totalCells]; + _zVelocityOb = new float[_totalCells]; _xVelocityOld = new float[_totalCells]; _yVelocityOld = new float[_totalCells]; _zVelocityOld = new float[_totalCells]; @@ -111,6 +116,9 @@ FLUID_3D::FLUID_3D(int *res, float *p0) : _xVelocity[x] = 0.0f; _yVelocity[x] = 0.0f; _zVelocity[x] = 0.0f; + _xVelocityOb[x] = 0.0f; + _yVelocityOb[x] = 0.0f; + _zVelocityOb[x] = 0.0f; _xVelocityOld[x] = 0.0f; _yVelocityOld[x] = 0.0f; _zVelocityOld[x] = 0.0f; @@ -131,9 +139,15 @@ FLUID_3D::FLUID_3D(int *res, float *p0) : _colloPrev = 1; // default value + setBorderObstacles(); // walls + +} +void FLUID_3D::setBorderObstacles() +{ + // set side obstacles - int index; + unsigned int index; for (int y = 0; y < _yRes; y++) for (int x = 0; x < _xRes; x++) { @@ -169,7 +183,6 @@ FLUID_3D::FLUID_3D(int *res, float *p0) : index += _xRes - 1; if(_domainBcRight==1) _obstacles[index] = 1; } - } FLUID_3D::~FLUID_3D() @@ -177,6 +190,9 @@ FLUID_3D::~FLUID_3D() if (_xVelocity) delete[] _xVelocity; if (_yVelocity) delete[] _yVelocity; if (_zVelocity) delete[] _zVelocity; + if (_xVelocityOb) delete[] _xVelocityOb; + if (_yVelocityOb) delete[] _yVelocityOb; + if (_zVelocityOb) delete[] _zVelocityOb; if (_xVelocityOld) delete[] _xVelocityOld; if (_yVelocityOld) delete[] _yVelocityOld; if (_zVelocityOld) delete[] _zVelocityOld; @@ -214,10 +230,18 @@ void FLUID_3D::initBlenderRNA(float *alpha, float *beta, float *dt_factor, float ////////////////////////////////////////////////////////////////////// void FLUID_3D::step(float dt) { +#if 0 // If border rules have been changed if (_colloPrev != *_borderColli) { + printf("Border collisions changed\n"); + + // DG TODO: Need to check that no animated obstacle flags are overwritten setBorderCollisions(); } +#endif + + // DG: TODO for the moment redo border for every timestep since it's been deleted every time by moving obstacles + setBorderCollisions(); // set delta time by dt_factor @@ -786,6 +810,7 @@ void FLUID_3D::project() memset(_pressure, 0, sizeof(float)*_totalCells); memset(_divergence, 0, sizeof(float)*_totalCells); + // set velocity and pressure inside of obstacles to zero setObstacleBoundaries(_pressure, 0, _zRes); // copy out the boundaries @@ -798,12 +823,49 @@ void FLUID_3D::project() if(_domainBcTop == 0) 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) for (y = 1; y < _yRes - 1; y++, index += 2) for (x = 1; x < _xRes - 1; x++, index++) { + + if(_obstacles[index]) + { + _divergence[index] = 0.0f; + continue; + } + + float xright = _xVelocity[index + 1]; float xleft = _xVelocity[index - 1]; float yup = _yVelocity[index + _xRes]; @@ -811,26 +873,82 @@ void FLUID_3D::project() float ztop = _zVelocity[index + _slabSize]; float zbottom = _zVelocity[index - _slabSize]; - if(_obstacles[index+1]) xright = - _xVelocity[index]; + if(_obstacles[index+1]) xright = - _xVelocity[index]; // DG: += if(_obstacles[index-1]) xleft = - _xVelocity[index]; if(_obstacles[index+_xRes]) yup = - _yVelocity[index]; if(_obstacles[index-_xRes]) ydown = - _yVelocity[index]; if(_obstacles[index+_slabSize]) ztop = - _zVelocity[index]; if(_obstacles[index-_slabSize]) zbottom = - _zVelocity[index]; + if(_obstacles[index+1] & 8) xright += _xVelocityOb[index + 1]; + if(_obstacles[index-1] & 8) xleft += _xVelocityOb[index - 1]; + if(_obstacles[index+_xRes] & 8) yup += _yVelocityOb[index + _xRes]; + if(_obstacles[index-_xRes] & 8) ydown += _yVelocityOb[index - _xRes]; + if(_obstacles[index+_slabSize] & 8) ztop += _zVelocityOb[index + _slabSize]; + if(_obstacles[index-_slabSize] & 8) zbottom += _zVelocityOb[index - _slabSize]; + _divergence[index] = -_dx * 0.5f * ( xright - xleft + yup - ydown + ztop - zbottom ); - // DG: commenting this helps CG to get a better start, 10-20% speed improvement - // _pressure[index] = 0.0f; + // Pressure is zero anyway since now a local array is used + _pressure[index] = 0.0f; } + + + /* + { + float maxvalue = 0; + for(unsigned int i = 0; i < _xRes * _yRes * _zRes; i++) + { + if(_divergence[i] > maxvalue) + maxvalue = _divergence[i]; + + } + printf("Max divergence: %f\n", maxvalue); + } + */ + copyBorderAll(_pressure, 0, _zRes); + /* + { + float maxvalue = 0; + for(unsigned int i = 0; i < _xRes * _yRes * _zRes; i++) + { + if(_pressure[i] > maxvalue) + maxvalue = _pressure[i]; + } + printf("Max pressure BEFORE: %f\n", maxvalue); + } + */ + // solve Poisson equation solvePressurePre(_pressure, _divergence, _obstacles); + { + float maxvalue = 0; + for(unsigned int i = 0; i < _xRes * _yRes * _zRes; i++) + { + if(_pressure[i] > maxvalue) + maxvalue = _pressure[i]; + + /* HACK: Animated collision object sometimes result in a non converging solvePressurePre() */ + if(_pressure[i] > _dx * _dt) + _pressure[i] = _dx * _dt; + else if(_pressure[i] < -_dx * _dt) + _pressure[i] = -_dx * _dt; + + // if(_obstacle[i] && _pressure[i] != 0.0) + // printf("BAD PRESSURE i\n"); + + // if(_pressure[i]>1) + // printf("index: %d\n", i); + } + // printf("Max pressure: %f, dx: %f\n", maxvalue, _dx); + } + setObstaclePressure(_pressure, 0, _zRes); // project out solution @@ -848,12 +966,74 @@ void FLUID_3D::project() } } + setObstacleVelocity(0, _zRes); + if (_pressure) delete[] _pressure; if (_divergence) delete[] _divergence; } +////////////////////////////////////////////////////////////////////// +// calculate the obstacle velocity at boundary +////////////////////////////////////////////////////////////////////// +void FLUID_3D::setObstacleVelocity(int zBegin, int zEnd) +{ + + // completely TODO <-- who wrote this and what is here TODO? DG + const size_t index_ = _slabSize + _xRes + 1; + + //int vIndex=_slabSize + _xRes + 1; + + int bb=0; + int bt=0; + + if (zBegin == 0) {bb = 1;} + if (zEnd == _zRes) {bt = 1;} + // tag remaining obstacle blocks + for (int z = zBegin + bb; z < zEnd - bt; z++) + { + size_t index = index_ +(z-1)*_slabSize; + + for (int y = 1; y < _yRes - 1; y++, index += 2) + { + for (int x = 1; x < _xRes - 1; x++, index++) + { + if (!_obstacles[index]) + { + // if(_obstacles[index+1]) xright = - _xVelocityOb[index]; + if((_obstacles[index - 1] & 8) && abs(_xVelocityOb[index - 1]) > FLT_EPSILON ) + { + // printf("velocity x!\n"); + _xVelocity[index] = _xVelocityOb[index - 1]; + _xVelocity[index - 1] = _xVelocityOb[index - 1]; + } + // if(_obstacles[index+_xRes]) yup = - _yVelocityOb[index]; + if((_obstacles[index - _xRes] & 8) && abs(_yVelocityOb[index - _xRes]) > FLT_EPSILON) + { + // printf("velocity y!\n"); + _yVelocity[index] = _yVelocityOb[index - _xRes]; + _yVelocity[index - _xRes] = _yVelocityOb[index - _xRes]; + } + // if(_obstacles[index+_slabSize]) ztop = - _zVelocityOb[index]; + if((_obstacles[index - _slabSize] & 8) && abs(_zVelocityOb[index - _slabSize]) > FLT_EPSILON) + { + // printf("velocity z!\n"); + _zVelocity[index] = _zVelocityOb[index - _slabSize]; + _zVelocity[index - _slabSize] = _zVelocityOb[index - _slabSize]; + } + } + else + { + _density[index] = 0; + } + //vIndex++; + } // x loop + //vIndex += 2; + } // y loop + //vIndex += 2 * _xRes; + } // z loop +} ////////////////////////////////////////////////////////////////////// // diffuse heat @@ -892,7 +1072,7 @@ void FLUID_3D::addObstacle(OBSTACLE* obstacle) void FLUID_3D::setObstaclePressure(float *_pressure, int zBegin, int zEnd) { - // compleately TODO + // completely TODO <-- who wrote this and what is here TODO? DG const size_t index_ = _slabSize + _xRes + 1; @@ -914,7 +1094,7 @@ void FLUID_3D::setObstaclePressure(float *_pressure, int zBegin, int zEnd) for (int x = 1; x < _xRes - 1; x++, index++) { // could do cascade of ifs, but they are a pain - if (_obstacles[index]) + if (_obstacles[index] /* && !(_obstacles[index] & 8) DG TODO TEST THIS CONDITION */) { const int top = _obstacles[index + _slabSize]; const int bottom= _obstacles[index - _slabSize]; @@ -928,9 +1108,11 @@ void FLUID_3D::setObstaclePressure(float *_pressure, int zBegin, int zEnd) // const bool fully = (up && down); //const bool fullx = (left && right); + /* _xVelocity[index] = _yVelocity[index] = _zVelocity[index] = 0.0f; + */ _pressure[index] = 0.0f; // average pressure neighbors @@ -1253,7 +1435,35 @@ void FLUID_3D::advectMacCormackEnd2(int zBegin, int zEnd) 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 (zBegin == 0) {bb = 1;} + if (zEnd == _zRes) {bt = 1;} + + for (int z = zBegin + bb; z < zEnd - bt; z++) + { + size_t index = index_ +(z-1)*_slabSize; + 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; + } + } + } + } + } +#endif /*int begin=zBegin * _slabSize; int end=begin + (zEnd - zBegin) * _slabSize; diff --git a/intern/smoke/intern/FLUID_3D.h b/intern/smoke/intern/FLUID_3D.h index c9e18926fb2..5704cba3ed4 100644 --- a/intern/smoke/intern/FLUID_3D.h +++ b/intern/smoke/intern/FLUID_3D.h @@ -39,9 +39,6 @@ // #include "WTURBULENCE.h" #include "VEC3.h" -// timestep default value for nice appearance -#define DT_DEFAULT 0.1f; - using namespace std; using namespace BasicVector; class WTURBULENCE; @@ -49,7 +46,7 @@ class WTURBULENCE; class FLUID_3D { public: - FLUID_3D(int *res, /* int amplify, */ float *p0); + FLUID_3D(int *res, /* int amplify, */ float *p0, float dtdef); FLUID_3D() {}; virtual ~FLUID_3D(); @@ -72,7 +69,7 @@ class FLUID_3D int yRes() const { return _yRes; }; int zRes() const { return _zRes; }; - public: + public: // dimensions int _xRes, _yRes, _zRes, _maxRes; Vec3Int _res; @@ -89,6 +86,8 @@ class FLUID_3D void artificialDampingSL(int zBegin, int zEnd); void artificialDampingExactSL(int pos); + void setBorderObstacles(); + // fields float* _density; float* _densityOld; @@ -97,13 +96,17 @@ class FLUID_3D float* _xVelocity; float* _yVelocity; float* _zVelocity; + float* _xVelocityOb; + float* _yVelocityOb; + float* _zVelocityOb; float* _xVelocityOld; float* _yVelocityOld; float* _zVelocityOld; float* _xForce; float* _yForce; float* _zForce; - unsigned char* _obstacles; + unsigned char* _obstacles; /* only used (usefull) for static obstacles like domain boundaries */ + unsigned char* _obstaclesAnim; // Required for proper threading: float* _xVelocityTemp; @@ -137,6 +140,8 @@ class FLUID_3D // have to recalibrate borders if nothing has changed void setBorderCollisions(); + void setObstacleVelocity(int zBegin, int zEnd); + // WTURBULENCE object, if active // WTURBULENCE* _wTurbulence; diff --git a/intern/smoke/intern/OBSTACLE.h b/intern/smoke/intern/OBSTACLE.h index 61d47b727f0..da8ec6be024 100644 --- a/intern/smoke/intern/OBSTACLE.h +++ b/intern/smoke/intern/OBSTACLE.h @@ -27,9 +27,11 @@ #define OBSTACLE_H enum OBSTACLE_FLAGS { - EMPTY = 0, + EMPTY = 0, + /* 1 is used to flag an object cell */ MARCHED = 2, - RETIRED = 4 + RETIRED = 4, + ANIMATED = 8, }; class OBSTACLE diff --git a/intern/smoke/intern/WTURBULENCE.cpp b/intern/smoke/intern/WTURBULENCE.cpp index cd18cf7b344..83bec466c9f 100644 --- a/intern/smoke/intern/WTURBULENCE.cpp +++ b/intern/smoke/intern/WTURBULENCE.cpp @@ -431,8 +431,11 @@ void WTURBULENCE::decomposeEnergy(float *_energy, float *_highFreqEnergy) // compute velocity from energies and march into obstacles // for wavelet decomposition ////////////////////////////////////////////////////////////////////// -void WTURBULENCE::computeEnergy(float *_energy, float* xvel, float* yvel, float* zvel, unsigned char *obstacles) +void WTURBULENCE::computeEnergy(float *_energy, float* xvel, float* yvel, float* zvel, unsigned char *origObstacles) { + unsigned char *obstacles = new unsigned char[_totalCellsSm]; + memcpy(obstacles, origObstacles, sizeof(unsigned char) * _totalCellsSm); + // compute everywhere for (int x = 0; x < _totalCellsSm; x++) _energy[x] = 0.5f * (xvel[x] * xvel[x] + yvel[x] * yvel[x] + zvel[x] * zvel[x]); @@ -506,7 +509,9 @@ void WTURBULENCE::computeEnergy(float *_energy, float* xvel, float* yvel, float* for (int y = 1; y < _yResSm - 1; y++, index += 2) for (int x = 1; x < _xResSm - 1; x++, index++) if (obstacles[index]) - obstacles[index] = 1; + obstacles[index] = 1; // DG TODO ? animated obstacle flag? + + free(obstacles); } ////////////////////////////////////////////////////////////////////////////////////////// diff --git a/intern/smoke/intern/smoke_API.cpp b/intern/smoke/intern/smoke_API.cpp index a2f3c21bbbf..78f7d35360a 100644 --- a/intern/smoke/intern/smoke_API.cpp +++ b/intern/smoke/intern/smoke_API.cpp @@ -19,6 +19,7 @@ * All rights reserved. * * Contributor(s): Daniel Genrich + * Blender Foundation * * ***** END GPL LICENSE BLOCK ***** */ @@ -36,10 +37,10 @@ #include <math.h> // y in smoke is z in blender -extern "C" FLUID_3D *smoke_init(int *res, float *p0) +extern "C" FLUID_3D *smoke_init(int *res, float *p0, float dtdef) { // smoke lib uses y as top-bottom/vertical axis where blender uses z - FLUID_3D *fluid = new FLUID_3D(res, p0); + FLUID_3D *fluid = new FLUID_3D(res, p0, dtdef); // printf("xres: %d, yres: %d, zres: %d\n", res[0], res[1], res[2]); @@ -78,41 +79,9 @@ 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, size_t framenr, float fps) +extern "C" void smoke_step(FLUID_3D *fluid, float dtSubdiv) { - /* stability values copied from wturbulence.cpp */ - const int maxSubSteps = 25; - const float maxVel = 0.5f; /* TODO: maybe 0.5 is still too high, please confirm! -dg */ - - float dt = DT_DEFAULT; - float maxVelMag = 0.0f; - int totalSubsteps; - int substep = 0; - float dtSubdiv; - - /* get max velocity and lower the dt value if it is too high */ - size_t size= fluid->_xRes * fluid->_yRes * fluid->_zRes; - - for(size_t i = 0; i < size; i++) - { - float vtemp = (fluid->_xVelocity[i]*fluid->_xVelocity[i]+fluid->_yVelocity[i]*fluid->_yVelocity[i]+fluid->_zVelocity[i]*fluid->_zVelocity[i]); - if(vtemp > maxVelMag) - maxVelMag = vtemp; - } - - /* adapt timestep for different framerates, dt = 0.1 is at 25fps */ - dt *= (25.0f / fps); - - maxVelMag = sqrt(maxVelMag) * dt * (*(fluid->_dtFactor)); - totalSubsteps = (int)((maxVelMag / maxVel) + 1.0f); /* always round up */ - totalSubsteps = (totalSubsteps < 1) ? 1 : totalSubsteps; - totalSubsteps = (totalSubsteps > maxSubSteps) ? maxSubSteps : totalSubsteps; - dtSubdiv = (float)dt / (float)totalSubsteps; - - // printf("totalSubsteps: %d, maxVelMag: %f, dt: %f\n", totalSubsteps, maxVelMag, dt); - - for(substep = 0; substep < totalSubsteps; substep++) - fluid->step(dtSubdiv); + fluid->step(dtSubdiv); } extern "C" void smoke_turbulence_step(WTURBULENCE *wt, FLUID_3D *fluid) @@ -307,6 +276,18 @@ extern "C" unsigned char *smoke_get_obstacle(FLUID_3D *fluid) return fluid->_obstacles; } +extern "C" void smoke_get_ob_velocity(struct FLUID_3D *fluid, float **x, float **y, float **z) +{ + *x = fluid->_xVelocityOb; + *y = fluid->_yVelocityOb; + *z = fluid->_zVelocityOb; +} + +extern "C" unsigned char *smoke_get_obstacle_anim(FLUID_3D *fluid) +{ + return fluid->_obstaclesAnim; +} + extern "C" void smoke_turbulence_set_noise(WTURBULENCE *wt, int type) { wt->setNoise(type); |