From 967d25ac1c269f28ab793f7c63edb9fd540dd5cf Mon Sep 17 00:00:00 2001 From: Daniel Genrich Date: Tue, 27 Jul 2010 14:53:20 +0000 Subject: =?UTF-8?q?Smoke=20Patch=20+=20additions:=20a)=20Applying=20patch?= =?UTF-8?q?=20#22765=20by=20Miika=20H=C3=A4m=C3=A4l=C3=A4inen=20(domain=20?= =?UTF-8?q?border=20collision=20settings,=20vorticity=20settings,=20time?= =?UTF-8?q?=20scale,=20non=20absolute=20density,=20smooth=20high=20res=20e?= =?UTF-8?q?mitter,=20initial=20velocity=20multiplier,=20high=20res=20stren?= =?UTF-8?q?gth=20available=20to=20be=20set=20to=200),=20b)=20Additions=20b?= =?UTF-8?q?y=20me:=20--Initial=20velocity=20is=20now=20per=20flow=20object?= =?UTF-8?q?,=20not=20per=20domain;=20--Using=20boundingbox=20as=20standard?= =?UTF-8?q?=20display=20mode=20for=20domains=20(was=20wire=20before);=20--?= =?UTF-8?q?When=20adding=20a=20flow=20object,=20an=20initial=20nice=20Smok?= =?UTF-8?q?eParticle=20system=20is=20added=20too=20with=20nice=20initial?= =?UTF-8?q?=20settings=20(life=3D1,=20no=5Frender,=20unborn,=20etc)=20fitt?= =?UTF-8?q?ing=20smoke=20simulation;=20--Adaptive=20timesteps=20introduced?= =?UTF-8?q?=20to=20the=20smoke=20sim=20(depending=20on=20the=20magnitude?= =?UTF-8?q?=20of=20the=20velocity)=20because=20it=20was=20quite=20unstable?= =?UTF-8?q?=20when=20used=20for=20fire=20simulations,=20still=20needs=20to?= =?UTF-8?q?=20be=20tested=20and=20will=20also=20slow=20down=20some=20simul?= =?UTF-8?q?ations.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- intern/smoke/extern/smoke_API.h | 4 +- intern/smoke/intern/FLUID_3D.cpp | 197 +++++++++++++++++++++++++++++++------- intern/smoke/intern/FLUID_3D.h | 23 ++++- intern/smoke/intern/smoke_API.cpp | 39 +++++++- 4 files changed, 217 insertions(+), 46 deletions(-) (limited to 'intern/smoke') diff --git a/intern/smoke/extern/smoke_API.h b/intern/smoke/extern/smoke_API.h index 3e296fd7c90..d1012eef76e 100644 --- a/intern/smoke/extern/smoke_API.h +++ b/intern/smoke/extern/smoke_API.h @@ -38,10 +38,10 @@ struct FLUID_3D; 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); +struct FLUID_3D *smoke_init(int *res, float *p0); void smoke_free(struct FLUID_3D *fluid); -void smoke_initBlenderRNA(struct FLUID_3D *fluid, float *alpha, float *beta); +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, size_t framenr); float *smoke_get_density(struct FLUID_3D *fluid); diff --git a/intern/smoke/intern/FLUID_3D.cpp b/intern/smoke/intern/FLUID_3D.cpp index 4ac960b5ef0..05fbb918d24 100644 --- a/intern/smoke/intern/FLUID_3D.cpp +++ b/intern/smoke/intern/FLUID_3D.cpp @@ -35,23 +35,15 @@ #include #endif // PARALLEL -// boundary conditions of the fluid domain -#define DOMAIN_BC_FRONT 0 // z -#define DOMAIN_BC_TOP 1 // y -#define DOMAIN_BC_LEFT 1 // x -#define DOMAIN_BC_BACK DOMAIN_BC_FRONT -#define DOMAIN_BC_BOTTOM DOMAIN_BC_TOP -#define DOMAIN_BC_RIGHT DOMAIN_BC_LEFT - ////////////////////////////////////////////////////////////////////// // Construction/Destruction ////////////////////////////////////////////////////////////////////// -FLUID_3D::FLUID_3D(int *res, float *p0, float dt) : - _xRes(res[0]), _yRes(res[1]), _zRes(res[2]), _res(0.0f), _dt(dt) +FLUID_3D::FLUID_3D(int *res, float *p0) : + _xRes(res[0]), _yRes(res[1]), _zRes(res[2]), _res(0.0f) { // set simulation consts - // _dt = dt; // 0.10 + _dt = DT_DEFAULT; // just in case. set in step from a RNA factor // start point of array _p0[0] = p0[0]; @@ -61,7 +53,6 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dt) : _iterations = 100; _tempAmb = 0; _heatDiffusion = 1e-3; - _vorticityEps = 2.0; _totalTime = 0.0f; _totalSteps = 0; _res = Vec3Int(_xRes,_yRes,_zRes); @@ -77,9 +68,9 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dt) : // scale the constants according to the refinement of the grid _dx = 1.0f / (float)_maxRes; - float scaling = 64.0f / _maxRes; - scaling = (scaling < 1.0f) ? 1.0f : scaling; - _vorticityEps /= scaling; + _constantScaling = 64.0f / _maxRes; + _constantScaling = (_constantScaling < 1.0f) ? 1.0f : _constantScaling; + _vorticityEps = 2.0f / _constantScaling; // Just in case set a default value // allocate arrays _totalCells = _xRes * _yRes * _zRes; @@ -126,30 +117,42 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dt) : _obstacles[x] = false; } + // boundary conditions of the fluid domain + // set default values -> vertically non-colliding + _domainBcFront = true; + _domainBcTop = false; + _domainBcLeft = true; + _domainBcBack = _domainBcFront; + _domainBcBottom = _domainBcTop; + _domainBcRight = _domainBcLeft; + + _colloPrev = 1; // default value + + // set side obstacles int index; for (int y = 0; y < _yRes; y++) for (int x = 0; x < _xRes; x++) { - // front slab + // bottom slab index = x + y * _xRes; - if(DOMAIN_BC_FRONT==1) _obstacles[index] = 1; + if(_domainBcBottom==1) _obstacles[index] = 1; - // back slab + // top slab index += _totalCells - _slabSize; - if(DOMAIN_BC_BACK==1) _obstacles[index] = 1; + if(_domainBcTop==1) _obstacles[index] = 1; } for (int z = 0; z < _zRes; z++) for (int x = 0; x < _xRes; x++) { - // bottom slab + // front slab index = x + z * _slabSize; - if(DOMAIN_BC_BOTTOM==1) _obstacles[index] = 1; + if(_domainBcFront==1) _obstacles[index] = 1; - // top slab + // back slab index += _slabSize - _xRes; - if(DOMAIN_BC_TOP==1) _obstacles[index] = 1; + if(_domainBcBack==1) _obstacles[index] = 1; } for (int z = 0; z < _zRes; z++) @@ -157,12 +160,13 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dt) : { // left slab index = y * _xRes + z * _slabSize; - if(DOMAIN_BC_LEFT==1) _obstacles[index] = 1; + if(_domainBcLeft==1) _obstacles[index] = 1; // right slab index += _xRes - 1; - if(DOMAIN_BC_RIGHT==1) _obstacles[index] = 1; + if(_domainBcRight==1) _obstacles[index] = 1; } + } FLUID_3D::~FLUID_3D() @@ -193,17 +197,32 @@ FLUID_3D::~FLUID_3D() } // init direct access functions from blender -void FLUID_3D::initBlenderRNA(float *alpha, float *beta) +void FLUID_3D::initBlenderRNA(float *alpha, float *beta, float *dt_factor, float *vorticity, int *borderCollision) { _alpha = alpha; _beta = beta; + _dtFactor = dt_factor; + _vorticityRNA = vorticity; + _borderColli = borderCollision; } ////////////////////////////////////////////////////////////////////// // step simulation once ////////////////////////////////////////////////////////////////////// -void FLUID_3D::step() +void FLUID_3D::step(float dt) { + // If border rules have been changed + if (_colloPrev != *_borderColli) { + setBorderCollisions(); + } + + + // set delta time by dt_factor + _dt = (*_dtFactor) * dt; + // set vorticity from RNA value + _vorticityEps = (*_vorticityRNA)/_constantScaling; + + #if PARALLEL==1 int threadval = 1; threadval = omp_get_max_threads(); @@ -246,6 +265,13 @@ void FLUID_3D::step() #pragma omp single { #endif + /* + * addForce() changed Temp values to preserve thread safety + * (previous functions in per thread loop still needed + * original velocity data) + * + * So swap temp values to velocity + */ SWAP_POINTERS(_xVelocity, _xVelocityTemp); SWAP_POINTERS(_yVelocity, _yVelocityTemp); SWAP_POINTERS(_zVelocity, _zVelocityTemp); @@ -276,6 +302,10 @@ void FLUID_3D::step() #pragma omp single { #endif + /* + * For thread safety use "Old" to read + * "current" values but still allow changing values. + */ SWAP_POINTERS(_xVelocity, _xVelocityOld); SWAP_POINTERS(_yVelocity, _yVelocityOld); SWAP_POINTERS(_zVelocity, _zVelocityOld); @@ -334,6 +364,10 @@ void FLUID_3D::step() } #endif + /* + * swap final velocity back to Velocity array + * from temp xForce storage + */ SWAP_POINTERS(_xVelocity, _xForce); SWAP_POINTERS(_yVelocity, _yForce); SWAP_POINTERS(_zVelocity, _zForce); @@ -351,6 +385,88 @@ void FLUID_3D::step() } + +// Set border collision model from RNA setting + +void FLUID_3D::setBorderCollisions() { + + + _colloPrev = *_borderColli; // saving the current value + + // boundary conditions of the fluid domain + if (_colloPrev == 0) + { + // No collisions + _domainBcFront = false; + _domainBcTop = false; + _domainBcLeft = false; + } + else if (_colloPrev == 2) + { + // Collide with all sides + _domainBcFront = true; + _domainBcTop = true; + _domainBcLeft = true; + } + else + { + // Default values: Collide with "walls", but not top and bottom + _domainBcFront = true; + _domainBcTop = false; + _domainBcLeft = true; + } + + _domainBcBack = _domainBcFront; + _domainBcBottom = _domainBcTop; + _domainBcRight = _domainBcLeft; + + + + // set side obstacles + int index; + for (int y = 0; y < _yRes; y++) + for (int x = 0; x < _xRes; x++) + { + // front slab + index = x + y * _xRes; + if(_domainBcBottom==1) _obstacles[index] = 1; + else _obstacles[index] = 0; + + // back slab + index += _totalCells - _slabSize; + if(_domainBcTop==1) _obstacles[index] = 1; + else _obstacles[index] = 0; + } + + for (int z = 0; z < _zRes; z++) + for (int x = 0; x < _xRes; x++) + { + // bottom slab + index = x + z * _slabSize; + if(_domainBcFront==1) _obstacles[index] = 1; + else _obstacles[index] = 0; + + // top slab + index += _slabSize - _xRes; + if(_domainBcBack==1) _obstacles[index] = 1; + else _obstacles[index] = 0; + } + + for (int z = 0; z < _zRes; z++) + for (int y = 0; y < _yRes; y++) + { + // left slab + index = y * _xRes + z * _slabSize; + if(_domainBcLeft==1) _obstacles[index] = 1; + else _obstacles[index] = 0; + + // right slab + index += _xRes - 1; + if(_domainBcRight==1) _obstacles[index] = 1; + else _obstacles[index] = 0; + } +} + ////////////////////////////////////////////////////////////////////// // helper function to dampen co-located grid artifacts of given arrays in intervals // (only needed for velocity, strength (w) depends on testcase... @@ -428,6 +544,10 @@ void FLUID_3D::artificialDampingExactSL(int pos) { for (y = 1; y < _res[1]-1; y++) for (x = 1+(y+z)%2; x < _res[0]-1; x+=2) { index = x + y*_res[0] + posslab; + /* + * Uses xForce as temporary storage to allow other threads to read + * old values from xVelocityTemp + */ _xForce[index] = (1-w)*_xVelocityTemp[index] + 1./6. * w*( _xVelocityTemp[index+1] + _xVelocityTemp[index-1] + _xVelocityTemp[index+_res[0]] + _xVelocityTemp[index-_res[0]] + @@ -450,6 +570,11 @@ void FLUID_3D::artificialDampingExactSL(int pos) { for (y = 1; y < _res[1]-1; y++) for (x = 1+(y+z+1)%2; x < _res[0]-1; x+=2) { index = x + y*_res[0] + posslab; + + /* + * Uses xForce as temporary storage to allow other threads to read + * old values from xVelocityTemp + */ _xForce[index] = (1-w)*_xVelocityTemp[index] + 1./6. * w*( _xVelocityTemp[index+1] + _xVelocityTemp[index-1] + _xVelocityTemp[index+_res[0]] + _xVelocityTemp[index-_res[0]] + @@ -661,13 +786,13 @@ void FLUID_3D::project() setObstacleBoundaries(_pressure, 0, _zRes); // copy out the boundaries - if(DOMAIN_BC_LEFT == 0) setNeumannX(_xVelocity, _res, 0, _zRes); + if(_domainBcLeft == 0) setNeumannX(_xVelocity, _res, 0, _zRes); else setZeroX(_xVelocity, _res, 0, _zRes); - if(DOMAIN_BC_TOP == 0) setNeumannY(_yVelocity, _res, 0, _zRes); + if(_domainBcFront == 0) setNeumannY(_yVelocity, _res, 0, _zRes); else setZeroY(_yVelocity, _res, 0, _zRes); - if(DOMAIN_BC_FRONT == 0) setNeumannZ(_zVelocity, _res, 0, _zRes); + if(_domainBcTop == 0) setNeumannZ(_zVelocity, _res, 0, _zRes); else setZeroZ(_zVelocity, _res, 0, _zRes); // calculate divergence @@ -1060,13 +1185,13 @@ void FLUID_3D::advectMacCormackBegin(int zBegin, int zEnd) { Vec3Int res = Vec3Int(_xRes,_yRes,_zRes); - if(DOMAIN_BC_LEFT == 0) copyBorderX(_xVelocityOld, res, zBegin, zEnd); + if(_domainBcLeft == 0) copyBorderX(_xVelocityOld, res, zBegin, zEnd); else setZeroX(_xVelocityOld, res, zBegin, zEnd); - if(DOMAIN_BC_TOP == 0) copyBorderY(_yVelocityOld, res, zBegin, zEnd); + if(_domainBcFront == 0) copyBorderY(_yVelocityOld, res, zBegin, zEnd); else setZeroY(_yVelocityOld, res, zBegin, zEnd); - if(DOMAIN_BC_FRONT == 0) copyBorderZ(_zVelocityOld, res, zBegin, zEnd); + if(_domainBcTop == 0) copyBorderZ(_zVelocityOld, res, zBegin, zEnd); else setZeroZ(_zVelocityOld, res, zBegin, zEnd); } @@ -1114,13 +1239,13 @@ void FLUID_3D::advectMacCormackEnd2(int zBegin, int 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); - if(DOMAIN_BC_LEFT == 0) copyBorderX(_xVelocityTemp, res, zBegin, zEnd); + if(_domainBcLeft == 0) copyBorderX(_xVelocityTemp, res, zBegin, zEnd); else setZeroX(_xVelocityTemp, res, zBegin, zEnd); - if(DOMAIN_BC_TOP == 0) copyBorderY(_yVelocityTemp, res, zBegin, zEnd); + if(_domainBcFront == 0) copyBorderY(_yVelocityTemp, res, zBegin, zEnd); else setZeroY(_yVelocityTemp, res, zBegin, zEnd); - if(DOMAIN_BC_FRONT == 0) copyBorderZ(_zVelocityTemp, res, zBegin, zEnd); + if(_domainBcTop == 0) copyBorderZ(_zVelocityTemp, res, zBegin, zEnd); else setZeroZ(_zVelocityTemp, res, zBegin, zEnd); setZeroBorder(_density, res, zBegin, zEnd); diff --git a/intern/smoke/intern/FLUID_3D.h b/intern/smoke/intern/FLUID_3D.h index 6ca50b2c032..c244180ee57 100644 --- a/intern/smoke/intern/FLUID_3D.h +++ b/intern/smoke/intern/FLUID_3D.h @@ -36,6 +36,9 @@ // #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; @@ -43,11 +46,11 @@ class WTURBULENCE; class FLUID_3D { public: - FLUID_3D(int *res, /* int amplify, */ float *p0, float dt); + FLUID_3D(int *res, /* int amplify, */ float *p0); FLUID_3D() {}; virtual ~FLUID_3D(); - void initBlenderRNA(float *alpha, float *beta); + void initBlenderRNA(float *alpha, float *beta, float *dt_factor, float *vorticity, int *border_colli); // create & allocate vector noise advection void initVectorNoise(int amplify); @@ -55,7 +58,7 @@ class FLUID_3D void addSmokeColumn(); static void addSmokeTestCase(float* field, Vec3Int res); - void step(); + void step(float dt); void addObstacle(OBSTACLE* obstacle); const float* xVelocity() { return _xVelocity; }; @@ -111,11 +114,25 @@ class FLUID_3D // simulation constants float _dt; + float *_dtFactor; float _vorticityEps; float _heatDiffusion; + float *_vorticityRNA; // RNA-pointer. float *_alpha; // for the buoyancy density term <-- as pointer to get blender RNA in here float *_beta; // was _buoyancy <-- as pointer to get blender RNA in here float _tempAmb; /* ambient temperature */ + float _constantScaling; + + bool _domainBcFront; // z + bool _domainBcTop; // y + bool _domainBcLeft; // x + bool _domainBcBack; // DOMAIN_BC_FRONT + bool _domainBcBottom; // DOMAIN_BC_TOP + bool _domainBcRight; // DOMAIN_BC_LEFT + int *_borderColli; // border collision rules <-- as pointer to get blender RNA in here + int _colloPrev; // To track whether value has been changed (to not + // have to recalibrate borders if nothing has changed + void setBorderCollisions(); // WTURBULENCE object, if active // WTURBULENCE* _wTurbulence; diff --git a/intern/smoke/intern/smoke_API.cpp b/intern/smoke/intern/smoke_API.cpp index 427f4d53171..23c12c17014 100644 --- a/intern/smoke/intern/smoke_API.cpp +++ b/intern/smoke/intern/smoke_API.cpp @@ -33,10 +33,10 @@ #include // y in smoke is z in blender -extern "C" FLUID_3D *smoke_init(int *res, float *p0, float dt) +extern "C" FLUID_3D *smoke_init(int *res, float *p0) { // smoke lib uses y as top-bottom/vertical axis where blender uses z - FLUID_3D *fluid = new FLUID_3D(res, p0, dt); + FLUID_3D *fluid = new FLUID_3D(res, p0); // printf("xres: %d, yres: %d, zres: %d\n", res[0], res[1], res[2]); @@ -77,7 +77,36 @@ extern "C" size_t smoke_get_index2d(int x, int max_x, int y /*, int max_y, int z extern "C" void smoke_step(FLUID_3D *fluid, size_t framenr) { - fluid->step(); + /* 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 */ + + const 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; + } + + 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); } extern "C" void smoke_turbulence_step(WTURBULENCE *wt, FLUID_3D *fluid) @@ -85,9 +114,9 @@ extern "C" void smoke_turbulence_step(WTURBULENCE *wt, FLUID_3D *fluid) 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) +extern "C" void smoke_initBlenderRNA(FLUID_3D *fluid, float *alpha, float *beta, float *dt_factor, float *vorticity, int *border_colli) { - fluid->initBlenderRNA(alpha, beta); + fluid->initBlenderRNA(alpha, beta, dt_factor, vorticity, border_colli); } extern "C" void smoke_dissolve(FLUID_3D *fluid, int speed, int log) -- cgit v1.2.3