Alexander Gavrilov
f593333099
The problem happens because smoke collides only with the surface of the collider and uses incompressible fluid solver. This means that scaling the collider tries to compress or decompress fluid within the volume of the collider, which can't be handled by the simulation. Fast rotation likely also causes transient scaling due to emulation of arcs by chords. This can be fixed by finding compartments completely isolated by obstacles from the rest of the domain, and forcing total divergence within each one to be zero so that equations are solvable. Physical validity is somewhat dubious, but without this the solver simply breaks down. From the physics point of view, the effect of the correction should be similar to opening a hole from every cell to another dimension that lets an equal amount of air to pass through to balance the change in volume. Reviewers: miikah, lukastoenne Reviewed By: lukastoenne Subscribers: dafassi, scorpion81, #physics Maniphest Tasks: T43220, T47551 Differential Revision: https://developer.blender.org/D2112
1793 lines
52 KiB
C++
1793 lines
52 KiB
C++
/** \file smoke/intern/FLUID_3D.cpp
|
|
* \ingroup smoke
|
|
*/
|
|
//////////////////////////////////////////////////////////////////////
|
|
// This file is part of Wavelet Turbulence.
|
|
//
|
|
// Wavelet Turbulence is free software: you can redistribute it and/or modify
|
|
// it under the terms of the GNU General Public License as published by
|
|
// the Free Software Foundation, either version 3 of the License, or
|
|
// (at your option) any later version.
|
|
//
|
|
// Wavelet Turbulence is distributed in the hope that it will be useful,
|
|
// but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
// GNU General Public License for more details.
|
|
//
|
|
// You should have received a copy of the GNU General Public License
|
|
// along with Wavelet Turbulence. If not, see <http://www.gnu.org/licenses/>.
|
|
//
|
|
// Copyright 2008 Theodore Kim and Nils Thuerey
|
|
//
|
|
// FLUID_3D.cpp: implementation of the FLUID_3D class.
|
|
//
|
|
//////////////////////////////////////////////////////////////////////
|
|
// Heavy parallel optimization done. Many of the old functions now
|
|
// take begin and end parameters and process only specified part of the data.
|
|
// Some functions were divided into multiple ones.
|
|
// - MiikaH
|
|
//////////////////////////////////////////////////////////////////////
|
|
|
|
#include "FLUID_3D.h"
|
|
#include "IMAGE.h"
|
|
#include <INTERPOLATE.h>
|
|
#include "SPHERE.h"
|
|
#include <zlib.h>
|
|
|
|
#include "float.h"
|
|
|
|
#if PARALLEL==1
|
|
#include <omp.h>
|
|
#endif // PARALLEL
|
|
|
|
//////////////////////////////////////////////////////////////////////
|
|
// Construction/Destruction
|
|
//////////////////////////////////////////////////////////////////////
|
|
|
|
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
|
|
|
|
_iterations = 100;
|
|
_tempAmb = 0;
|
|
_heatDiffusion = 1e-3;
|
|
_totalTime = 0.0f;
|
|
_totalSteps = 0;
|
|
_res = Vec3Int(_xRes,_yRes,_zRes);
|
|
_maxRes = MAX3(_xRes, _yRes, _zRes);
|
|
|
|
// initialize wavelet turbulence
|
|
/*
|
|
if(amplify)
|
|
_wTurbulence = new WTURBULENCE(_res[0],_res[1],_res[2], amplify, noisetype);
|
|
else
|
|
_wTurbulence = NULL;
|
|
*/
|
|
|
|
// scale the constants according to the refinement of the grid
|
|
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
|
|
|
|
// allocate arrays
|
|
_totalCells = _xRes * _yRes * _zRes;
|
|
_slabSize = _xRes * _yRes;
|
|
_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];
|
|
_xForce = new float[_totalCells];
|
|
_yForce = new float[_totalCells];
|
|
_zForce = new float[_totalCells];
|
|
_density = new float[_totalCells];
|
|
_densityOld = new float[_totalCells];
|
|
_obstacles = new unsigned char[_totalCells]; // set 0 at end of step
|
|
|
|
// For threaded version:
|
|
_xVelocityTemp = new float[_totalCells];
|
|
_yVelocityTemp = new float[_totalCells];
|
|
_zVelocityTemp = new float[_totalCells];
|
|
_densityTemp = new float[_totalCells];
|
|
|
|
// DG TODO: check if alloc went fine
|
|
|
|
for (int x = 0; x < _totalCells; x++)
|
|
{
|
|
_density[x] = 0.0f;
|
|
_densityOld[x] = 0.0f;
|
|
_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;
|
|
_xForce[x] = 0.0f;
|
|
_yForce[x] = 0.0f;
|
|
_zForce[x] = 0.0f;
|
|
_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;
|
|
_domainBcTop = false;
|
|
_domainBcLeft = true;
|
|
_domainBcBack = _domainBcFront;
|
|
_domainBcBottom = _domainBcTop;
|
|
_domainBcRight = _domainBcLeft;
|
|
|
|
_colloPrev = 1; // default value
|
|
}
|
|
|
|
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()
|
|
{
|
|
|
|
// set side obstacles
|
|
unsigned int index;
|
|
for (int y = 0; y < _yRes; y++)
|
|
for (int x = 0; x < _xRes; x++)
|
|
{
|
|
// bottom slab
|
|
index = x + y * _xRes;
|
|
if(_domainBcBottom) _obstacles[index] = 1;
|
|
|
|
// top slab
|
|
index += _totalCells - _slabSize;
|
|
if(_domainBcTop) _obstacles[index] = 1;
|
|
}
|
|
|
|
for (int z = 0; z < _zRes; z++)
|
|
for (int x = 0; x < _xRes; x++)
|
|
{
|
|
// front slab
|
|
index = x + z * _slabSize;
|
|
if(_domainBcFront) _obstacles[index] = 1;
|
|
|
|
// back slab
|
|
index += _slabSize - _xRes;
|
|
if(_domainBcBack) _obstacles[index] = 1;
|
|
}
|
|
|
|
for (int z = 0; z < _zRes; z++)
|
|
for (int y = 0; y < _yRes; y++)
|
|
{
|
|
// left slab
|
|
index = y * _xRes + z * _slabSize;
|
|
if(_domainBcLeft) _obstacles[index] = 1;
|
|
|
|
// right slab
|
|
index += _xRes - 1;
|
|
if(_domainBcRight) _obstacles[index] = 1;
|
|
}
|
|
}
|
|
|
|
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;
|
|
if (_xForce) delete[] _xForce;
|
|
if (_yForce) delete[] _yForce;
|
|
if (_zForce) delete[] _zForce;
|
|
if (_density) delete[] _density;
|
|
if (_densityOld) delete[] _densityOld;
|
|
if (_heat) delete[] _heat;
|
|
if (_heatOld) delete[] _heatOld;
|
|
if (_obstacles) delete[] _obstacles;
|
|
|
|
if (_xVelocityTemp) delete[] _xVelocityTemp;
|
|
if (_yVelocityTemp) delete[] _yVelocityTemp;
|
|
if (_zVelocityTemp) delete[] _zVelocityTemp;
|
|
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, 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, float gravity[3])
|
|
{
|
|
#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
|
|
_dt = (*_dtFactor) * dt;
|
|
// set vorticity from RNA value
|
|
_vorticityEps = (*_vorticityRNA)/_constantScaling;
|
|
|
|
#if PARALLEL==1
|
|
int threadval = 1;
|
|
threadval = omp_get_max_threads();
|
|
|
|
int stepParts = 1;
|
|
float partSize = _zRes;
|
|
|
|
stepParts = threadval*2; // Dividing parallelized sections into numOfThreads * 2 sections
|
|
partSize = (float)_zRes/stepParts; // Size of one part;
|
|
|
|
if (partSize < 4) {stepParts = threadval; // If the slice gets too low (might actually slow things down, change it to larger
|
|
partSize = (float)_zRes/stepParts;}
|
|
if (partSize < 4) {stepParts = (int)(ceil((float)_zRes/4.0f)); // If it's still too low (only possible on future systems with +24 cores), change it to 4
|
|
partSize = (float)_zRes/stepParts;}
|
|
#else
|
|
int zBegin=0;
|
|
int zEnd=_zRes;
|
|
#endif
|
|
|
|
wipeBoundariesSL(0, _zRes);
|
|
|
|
#if PARALLEL==1
|
|
#pragma omp parallel
|
|
{
|
|
#pragma omp for schedule(static,1)
|
|
for (int i=0; i<stepParts; i++)
|
|
{
|
|
int zBegin = (int)((float)i*partSize + 0.5f);
|
|
int zEnd = (int)((float)(i+1)*partSize + 0.5f);
|
|
#endif
|
|
|
|
addVorticity(zBegin, zEnd);
|
|
addBuoyancy(_heat, _density, gravity, zBegin, zEnd);
|
|
addForce(zBegin, zEnd);
|
|
|
|
#if PARALLEL==1
|
|
} // end of parallel
|
|
#pragma omp barrier
|
|
|
|
#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);
|
|
#if PARALLEL==1
|
|
} // end of single
|
|
|
|
#pragma omp barrier
|
|
|
|
#pragma omp for
|
|
for (int i=0; i<2; i++)
|
|
{
|
|
if (i==0)
|
|
{
|
|
#endif
|
|
project();
|
|
#if PARALLEL==1
|
|
}
|
|
else if (i==1)
|
|
{
|
|
#endif
|
|
if (_heat) {
|
|
diffuseHeat();
|
|
}
|
|
#if PARALLEL==1
|
|
}
|
|
}
|
|
|
|
#pragma omp barrier
|
|
|
|
#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);
|
|
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
|
|
} // end of single
|
|
|
|
#pragma omp barrier
|
|
|
|
|
|
#pragma omp for schedule(static,1)
|
|
for (int i=0; i<stepParts; i++)
|
|
{
|
|
|
|
int zBegin = (int)((float)i*partSize + 0.5f);
|
|
int zEnd = (int)((float)(i+1)*partSize + 0.5f);
|
|
#endif
|
|
|
|
advectMacCormackEnd1(zBegin, zEnd);
|
|
|
|
#if PARALLEL==1
|
|
} // end of parallel
|
|
|
|
#pragma omp barrier
|
|
|
|
#pragma omp for schedule(static,1)
|
|
for (int i=0; i<stepParts; i++)
|
|
{
|
|
|
|
int zBegin = (int)((float)i*partSize + 0.5f);
|
|
int zEnd = (int)((float)(i+1)*partSize + 0.5f);
|
|
#endif
|
|
|
|
advectMacCormackEnd2(zBegin, zEnd);
|
|
|
|
artificialDampingSL(zBegin, zEnd);
|
|
|
|
// Using forces as temp arrays
|
|
|
|
#if PARALLEL==1
|
|
}
|
|
}
|
|
|
|
|
|
|
|
for (int i=1; i<stepParts; i++)
|
|
{
|
|
int zPos=(int)((float)i*partSize + 0.5f);
|
|
|
|
artificialDampingExactSL(zPos);
|
|
|
|
}
|
|
#endif
|
|
|
|
/*
|
|
* swap final velocity back to Velocity array
|
|
* from temp xForce storage
|
|
*/
|
|
SWAP_POINTERS(_xVelocity, _xForce);
|
|
SWAP_POINTERS(_yVelocity, _yForce);
|
|
SWAP_POINTERS(_zVelocity, _zForce);
|
|
|
|
_totalTime += _dt;
|
|
_totalSteps++;
|
|
|
|
for (int i = 0; i < _totalCells; i++)
|
|
{
|
|
_xForce[i] = _yForce[i] = _zForce[i] = 0.0f;
|
|
}
|
|
|
|
}
|
|
|
|
|
|
// 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
|
|
setBorderObstacles();
|
|
}
|
|
|
|
//////////////////////////////////////////////////////////////////////
|
|
// helper function to dampen co-located grid artifacts of given arrays in intervals
|
|
// (only needed for velocity, strength (w) depends on testcase...
|
|
//////////////////////////////////////////////////////////////////////
|
|
|
|
|
|
void FLUID_3D::artificialDampingSL(int zBegin, int zEnd) {
|
|
const float w = 0.9;
|
|
|
|
memmove(_xForce+(_slabSize*zBegin), _xVelocityTemp+(_slabSize*zBegin), sizeof(float)*_slabSize*(zEnd-zBegin));
|
|
memmove(_yForce+(_slabSize*zBegin), _yVelocityTemp+(_slabSize*zBegin), sizeof(float)*_slabSize*(zEnd-zBegin));
|
|
memmove(_zForce+(_slabSize*zBegin), _zVelocityTemp+(_slabSize*zBegin), sizeof(float)*_slabSize*(zEnd-zBegin));
|
|
|
|
|
|
if(_totalSteps % 4 == 1) {
|
|
for (int z = zBegin+1; z < zEnd-1; z++)
|
|
for (int y = 1; y < _res[1]-1; y++)
|
|
for (int x = 1+(y+z)%2; x < _res[0]-1; x+=2) {
|
|
const int index = x + y*_res[0] + z * _slabSize;
|
|
_xForce[index] = (1-w)*_xVelocityTemp[index] + 1.0f/6.0f * w*(
|
|
_xVelocityTemp[index+1] + _xVelocityTemp[index-1] +
|
|
_xVelocityTemp[index+_res[0]] + _xVelocityTemp[index-_res[0]] +
|
|
_xVelocityTemp[index+_slabSize] + _xVelocityTemp[index-_slabSize] );
|
|
|
|
_yForce[index] = (1-w)*_yVelocityTemp[index] + 1.0f/6.0f * w*(
|
|
_yVelocityTemp[index+1] + _yVelocityTemp[index-1] +
|
|
_yVelocityTemp[index+_res[0]] + _yVelocityTemp[index-_res[0]] +
|
|
_yVelocityTemp[index+_slabSize] + _yVelocityTemp[index-_slabSize] );
|
|
|
|
_zForce[index] = (1-w)*_zVelocityTemp[index] + 1.0f/6.0f * w*(
|
|
_zVelocityTemp[index+1] + _zVelocityTemp[index-1] +
|
|
_zVelocityTemp[index+_res[0]] + _zVelocityTemp[index-_res[0]] +
|
|
_zVelocityTemp[index+_slabSize] + _zVelocityTemp[index-_slabSize] );
|
|
}
|
|
}
|
|
|
|
if(_totalSteps % 4 == 3) {
|
|
for (int z = zBegin+1; z < zEnd-1; z++)
|
|
for (int y = 1; y < _res[1]-1; y++)
|
|
for (int x = 1+(y+z+1)%2; x < _res[0]-1; x+=2) {
|
|
const int index = x + y*_res[0] + z * _slabSize;
|
|
_xForce[index] = (1-w)*_xVelocityTemp[index] + 1.0f/6.0f * w*(
|
|
_xVelocityTemp[index+1] + _xVelocityTemp[index-1] +
|
|
_xVelocityTemp[index+_res[0]] + _xVelocityTemp[index-_res[0]] +
|
|
_xVelocityTemp[index+_slabSize] + _xVelocityTemp[index-_slabSize] );
|
|
|
|
_yForce[index] = (1-w)*_yVelocityTemp[index] + 1.0f/6.0f * w*(
|
|
_yVelocityTemp[index+1] + _yVelocityTemp[index-1] +
|
|
_yVelocityTemp[index+_res[0]] + _yVelocityTemp[index-_res[0]] +
|
|
_yVelocityTemp[index+_slabSize] + _yVelocityTemp[index-_slabSize] );
|
|
|
|
_zForce[index] = (1-w)*_zVelocityTemp[index] + 1.0f/6.0f * w*(
|
|
_zVelocityTemp[index+1] + _zVelocityTemp[index-1] +
|
|
_zVelocityTemp[index+_res[0]] + _zVelocityTemp[index-_res[0]] +
|
|
_zVelocityTemp[index+_slabSize] + _zVelocityTemp[index-_slabSize] );
|
|
}
|
|
|
|
}
|
|
}
|
|
|
|
|
|
|
|
void FLUID_3D::artificialDampingExactSL(int pos) {
|
|
const float w = 0.9;
|
|
int index, x,y,z;
|
|
|
|
|
|
size_t posslab;
|
|
|
|
for (z=pos-1; z<=pos; z++)
|
|
{
|
|
posslab=z * _slabSize;
|
|
|
|
if(_totalSteps % 4 == 1) {
|
|
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.0f/6.0f * w*(
|
|
_xVelocityTemp[index+1] + _xVelocityTemp[index-1] +
|
|
_xVelocityTemp[index+_res[0]] + _xVelocityTemp[index-_res[0]] +
|
|
_xVelocityTemp[index+_slabSize] + _xVelocityTemp[index-_slabSize] );
|
|
|
|
_yForce[index] = (1-w)*_yVelocityTemp[index] + 1.0f/6.0f * w*(
|
|
_yVelocityTemp[index+1] + _yVelocityTemp[index-1] +
|
|
_yVelocityTemp[index+_res[0]] + _yVelocityTemp[index-_res[0]] +
|
|
_yVelocityTemp[index+_slabSize] + _yVelocityTemp[index-_slabSize] );
|
|
|
|
_zForce[index] = (1-w)*_zVelocityTemp[index] + 1.0f/6.0f * w*(
|
|
_zVelocityTemp[index+1] + _zVelocityTemp[index-1] +
|
|
_zVelocityTemp[index+_res[0]] + _zVelocityTemp[index-_res[0]] +
|
|
_zVelocityTemp[index+_slabSize] + _zVelocityTemp[index-_slabSize] );
|
|
|
|
}
|
|
}
|
|
|
|
if(_totalSteps % 4 == 3) {
|
|
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.0f/6.0f * w*(
|
|
_xVelocityTemp[index+1] + _xVelocityTemp[index-1] +
|
|
_xVelocityTemp[index+_res[0]] + _xVelocityTemp[index-_res[0]] +
|
|
_xVelocityTemp[index+_slabSize] + _xVelocityTemp[index-_slabSize] );
|
|
|
|
_yForce[index] = (1-w)*_yVelocityTemp[index] + 1.0f/6.0f * w*(
|
|
_yVelocityTemp[index+1] + _yVelocityTemp[index-1] +
|
|
_yVelocityTemp[index+_res[0]] + _yVelocityTemp[index-_res[0]] +
|
|
_yVelocityTemp[index+_slabSize] + _yVelocityTemp[index-_slabSize] );
|
|
|
|
_zForce[index] = (1-w)*_zVelocityTemp[index] + 1.0f/6.0f * w*(
|
|
_zVelocityTemp[index+1] + _zVelocityTemp[index-1] +
|
|
_zVelocityTemp[index+_res[0]] + _zVelocityTemp[index-_res[0]] +
|
|
_zVelocityTemp[index+_slabSize] + _zVelocityTemp[index-_slabSize] );
|
|
|
|
}
|
|
|
|
}
|
|
}
|
|
}
|
|
|
|
//////////////////////////////////////////////////////////////////////
|
|
// copy out the boundary in all directions
|
|
//////////////////////////////////////////////////////////////////////
|
|
void FLUID_3D::copyBorderAll(float* field, int zBegin, int zEnd)
|
|
{
|
|
int index, x, y, z;
|
|
int zSize = zEnd-zBegin;
|
|
int _blockTotalCells=_slabSize * zSize;
|
|
|
|
if (zBegin==0)
|
|
for (int y = 0; y < _yRes; y++)
|
|
for (int x = 0; x < _xRes; x++)
|
|
{
|
|
// front slab
|
|
index = x + y * _xRes;
|
|
field[index] = field[index + _slabSize];
|
|
}
|
|
|
|
if (zEnd==_zRes)
|
|
for (y = 0; y < _yRes; y++)
|
|
for (x = 0; x < _xRes; x++)
|
|
{
|
|
|
|
// back slab
|
|
index = x + y * _xRes + _blockTotalCells - _slabSize;
|
|
field[index] = field[index - _slabSize];
|
|
}
|
|
|
|
for (z = 0; z < zSize; z++)
|
|
for (x = 0; x < _xRes; x++)
|
|
{
|
|
// bottom slab
|
|
index = x + z * _slabSize;
|
|
field[index] = field[index + _xRes];
|
|
|
|
// top slab
|
|
index += _slabSize - _xRes;
|
|
field[index] = field[index - _xRes];
|
|
}
|
|
|
|
for (z = 0; z < zSize; z++)
|
|
for (y = 0; y < _yRes; y++)
|
|
{
|
|
// left slab
|
|
index = y * _xRes + z * _slabSize;
|
|
field[index] = field[index + 1];
|
|
|
|
// right slab
|
|
index += _xRes - 1;
|
|
field[index] = field[index - 1];
|
|
}
|
|
}
|
|
|
|
//////////////////////////////////////////////////////////////////////
|
|
// wipe boundaries of velocity and density
|
|
//////////////////////////////////////////////////////////////////////
|
|
void FLUID_3D::wipeBoundaries(int zBegin, int zEnd)
|
|
{
|
|
setZeroBorder(_xVelocity, _res, zBegin, 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)
|
|
{
|
|
|
|
/////////////////////////////////////
|
|
// setZeroBorder to all:
|
|
/////////////////////////////////////
|
|
|
|
/////////////////////////////////////
|
|
// setZeroX
|
|
/////////////////////////////////////
|
|
|
|
const int slabSize = _xRes * _yRes;
|
|
int index, x,y,z;
|
|
|
|
for (z = zBegin; z < zEnd; z++)
|
|
for (y = 0; y < _yRes; y++)
|
|
{
|
|
// left slab
|
|
index = y * _xRes + z * slabSize;
|
|
_xVelocity[index] = 0.0f;
|
|
_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;
|
|
_xVelocity[index] = 0.0f;
|
|
_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;
|
|
}
|
|
}
|
|
|
|
/////////////////////////////////////
|
|
// setZeroY
|
|
/////////////////////////////////////
|
|
|
|
for (z = zBegin; z < zEnd; z++)
|
|
for (x = 0; x < _xRes; x++)
|
|
{
|
|
// bottom slab
|
|
index = x + z * slabSize;
|
|
_xVelocity[index] = 0.0f;
|
|
_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;
|
|
_xVelocity[index] = 0.0f;
|
|
_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;
|
|
}
|
|
|
|
}
|
|
|
|
/////////////////////////////////////
|
|
// setZeroZ
|
|
/////////////////////////////////////
|
|
|
|
|
|
const int totalCells = _xRes * _yRes * _zRes;
|
|
|
|
index = 0;
|
|
if (zBegin == 0)
|
|
for (y = 0; y < _yRes; y++)
|
|
for (x = 0; x < _xRes; x++, index++)
|
|
{
|
|
// front slab
|
|
_xVelocity[index] = 0.0f;
|
|
_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)
|
|
{
|
|
index=0;
|
|
int index_top=0;
|
|
const int cellsslab = totalCells - slabSize;
|
|
|
|
for (y = 0; y < _yRes; y++)
|
|
for (x = 0; x < _xRes; x++, index++)
|
|
{
|
|
|
|
// back slab
|
|
index_top = index + cellsslab;
|
|
_xVelocity[index_top] = 0.0f;
|
|
_yVelocity[index_top] = 0.0f;
|
|
_zVelocity[index_top] = 0.0f;
|
|
_density[index_top] = 0.0f;
|
|
if (_fuel) {
|
|
_fuel[index_top] = 0.0f;
|
|
_react[index_top] = 0.0f;
|
|
}
|
|
if (_color_r) {
|
|
_color_r[index_top] = 0.0f;
|
|
_color_g[index_top] = 0.0f;
|
|
_color_b[index_top] = 0.0f;
|
|
}
|
|
}
|
|
}
|
|
|
|
}
|
|
//////////////////////////////////////////////////////////////////////
|
|
// add forces to velocity field
|
|
//////////////////////////////////////////////////////////////////////
|
|
void FLUID_3D::addForce(int zBegin, int zEnd)
|
|
{
|
|
int begin=zBegin * _slabSize;
|
|
int end=begin + (zEnd - zBegin) * _slabSize;
|
|
|
|
for (int i = begin; i < end; i++)
|
|
{
|
|
_xVelocityTemp[i] = _xVelocity[i] + _dt * _xForce[i];
|
|
_yVelocityTemp[i] = _yVelocity[i] + _dt * _yForce[i];
|
|
_zVelocityTemp[i] = _zVelocity[i] + _dt * _zForce[i];
|
|
}
|
|
}
|
|
//////////////////////////////////////////////////////////////////////
|
|
// project into divergence free field
|
|
//////////////////////////////////////////////////////////////////////
|
|
void FLUID_3D::project()
|
|
{
|
|
int x, y, z;
|
|
size_t index;
|
|
|
|
float *_pressure = new float[_totalCells];
|
|
float *_divergence = new float[_totalCells];
|
|
|
|
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
|
|
if(!_domainBcLeft) setNeumannX(_xVelocity, _res, 0, _zRes);
|
|
else setZeroX(_xVelocity, _res, 0, _zRes);
|
|
|
|
if(!_domainBcFront) setNeumannY(_yVelocity, _res, 0, _zRes);
|
|
else setZeroY(_yVelocity, _res, 0, _zRes);
|
|
|
|
if(!_domainBcTop) setNeumannZ(_zVelocity, _res, 0, _zRes);
|
|
else setZeroZ(_zVelocity, _res, 0, _zRes);
|
|
|
|
// 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];
|
|
float ydown = _yVelocity[index - _xRes];
|
|
float ztop = _zVelocity[index + _slabSize];
|
|
float zbottom = _zVelocity[index - _slabSize];
|
|
|
|
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 );
|
|
|
|
// Pressure is zero anyway since now a local array is used
|
|
_pressure[index] = 0.0f;
|
|
}
|
|
|
|
copyBorderAll(_pressure, 0, _zRes);
|
|
|
|
// fix fluid compression caused in isolated components by obstacle movement
|
|
fixObstacleCompression(_divergence);
|
|
|
|
// solve Poisson equation
|
|
solvePressurePre(_pressure, _divergence, _obstacles);
|
|
|
|
setObstaclePressure(_pressure, 0, _zRes);
|
|
|
|
// project out solution
|
|
// New idea for code from NVIDIA graphic gems 3 - DG
|
|
float invDx = 1.0f / _dx;
|
|
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++)
|
|
{
|
|
float vMask[3] = {1.0f, 1.0f, 1.0f}, vObst[3] = {0, 0, 0};
|
|
// float vR = 0.0f, vL = 0.0f, vT = 0.0f, vB = 0.0f, vD = 0.0f, vU = 0.0f; // UNUSED
|
|
|
|
float pC = _pressure[index]; // center
|
|
float pR = _pressure[index + 1]; // right
|
|
float pL = _pressure[index - 1]; // left
|
|
float pU = _pressure[index + _xRes]; // Up
|
|
float pD = _pressure[index - _xRes]; // Down
|
|
float pT = _pressure[index + _slabSize]; // top
|
|
float pB = _pressure[index - _slabSize]; // bottom
|
|
|
|
if(!_obstacles[index])
|
|
{
|
|
// DG TODO: What if obstacle is left + right and one of them is moving?
|
|
if(_obstacles[index+1]) { pR = pC; vObst[0] = _xVelocityOb[index + 1]; vMask[0] = 0; }
|
|
if(_obstacles[index-1]) { pL = pC; vObst[0] = _xVelocityOb[index - 1]; vMask[0] = 0; }
|
|
if(_obstacles[index+_xRes]) { pU = pC; vObst[1] = _yVelocityOb[index + _xRes]; vMask[1] = 0; }
|
|
if(_obstacles[index-_xRes]) { pD = pC; vObst[1] = _yVelocityOb[index - _xRes]; vMask[1] = 0; }
|
|
if(_obstacles[index+_slabSize]) { pT = pC; vObst[2] = _zVelocityOb[index + _slabSize]; vMask[2] = 0; }
|
|
if(_obstacles[index-_slabSize]) { pB = pC; vObst[2] = _zVelocityOb[index - _slabSize]; vMask[2] = 0; }
|
|
|
|
_xVelocity[index] -= 0.5f * (pR - pL) * invDx;
|
|
_yVelocity[index] -= 0.5f * (pU - pD) * invDx;
|
|
_zVelocity[index] -= 0.5f * (pT - pB) * invDx;
|
|
|
|
_xVelocity[index] = (vMask[0] * _xVelocity[index]) + vObst[0];
|
|
_yVelocity[index] = (vMask[1] * _yVelocity[index]) + vObst[1];
|
|
_zVelocity[index] = (vMask[2] * _zVelocity[index]) + vObst[2];
|
|
}
|
|
else
|
|
{
|
|
_xVelocity[index] = _xVelocityOb[index];
|
|
_yVelocity[index] = _yVelocityOb[index];
|
|
_zVelocity[index] = _zVelocityOb[index];
|
|
}
|
|
}
|
|
|
|
// DG: was enabled in original code but now we do this later
|
|
// 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
|
|
//////////////////////////////////////////////////////////////////////
|
|
void FLUID_3D::diffuseHeat()
|
|
{
|
|
SWAP_POINTERS(_heat, _heatOld);
|
|
|
|
copyBorderAll(_heatOld, 0, _zRes);
|
|
solveHeat(_heat, _heatOld, _obstacles);
|
|
|
|
// zero out inside obstacles
|
|
for (int x = 0; x < _totalCells; x++)
|
|
if (_obstacles[x])
|
|
_heat[x] = 0.0f;
|
|
}
|
|
|
|
//////////////////////////////////////////////////////////////////////
|
|
// stamp an obstacle in the _obstacles field
|
|
//////////////////////////////////////////////////////////////////////
|
|
void FLUID_3D::addObstacle(OBSTACLE* obstacle)
|
|
{
|
|
int index = 0;
|
|
for (int z = 0; z < _zRes; z++)
|
|
for (int y = 0; y < _yRes; y++)
|
|
for (int x = 0; x < _xRes; x++, index++)
|
|
if (obstacle->inside(x * _dx, y * _dx, z * _dx)) {
|
|
_obstacles[index] = true;
|
|
}
|
|
}
|
|
|
|
//////////////////////////////////////////////////////////////////////
|
|
// calculate the obstacle directional types
|
|
//////////////////////////////////////////////////////////////////////
|
|
void FLUID_3D::setObstaclePressure(float *_pressure, 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++)
|
|
{
|
|
// could do cascade of ifs, but they are a pain
|
|
if (_obstacles[index] /* && !(_obstacles[index] & 8) DG TODO TEST THIS CONDITION */)
|
|
{
|
|
const int top = _obstacles[index + _slabSize];
|
|
const int bottom= _obstacles[index - _slabSize];
|
|
const int up = _obstacles[index + _xRes];
|
|
const int down = _obstacles[index - _xRes];
|
|
const int left = _obstacles[index - 1];
|
|
const int right = _obstacles[index + 1];
|
|
|
|
// unused
|
|
// const bool fullz = (top && bottom);
|
|
// 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
|
|
float pcnt = 0.;
|
|
if (left && !right) {
|
|
_pressure[index] += _pressure[index + 1];
|
|
pcnt += 1.0f;
|
|
}
|
|
if (!left && right) {
|
|
_pressure[index] += _pressure[index - 1];
|
|
pcnt += 1.0f;
|
|
}
|
|
if (up && !down) {
|
|
_pressure[index] += _pressure[index - _xRes];
|
|
pcnt += 1.0f;
|
|
}
|
|
if (!up && down) {
|
|
_pressure[index] += _pressure[index + _xRes];
|
|
pcnt += 1.0f;
|
|
}
|
|
if (top && !bottom) {
|
|
_pressure[index] += _pressure[index - _slabSize];
|
|
pcnt += 1.0f;
|
|
}
|
|
if (!top && bottom) {
|
|
_pressure[index] += _pressure[index + _slabSize];
|
|
pcnt += 1.0f;
|
|
}
|
|
|
|
if(pcnt > 0.000001f)
|
|
_pressure[index] /= pcnt;
|
|
|
|
// TODO? set correct velocity bc's
|
|
// velocities are only set to zero right now
|
|
// this means it's not a full no-slip boundary condition
|
|
// but a "half-slip" - still looks ok right now
|
|
}
|
|
//vIndex++;
|
|
} // x loop
|
|
//vIndex += 2;
|
|
} // y loop
|
|
//vIndex += 2 * _xRes;
|
|
} // z loop
|
|
}
|
|
|
|
void FLUID_3D::setObstacleBoundaries(float *_pressure, int zBegin, int zEnd)
|
|
{
|
|
// cull degenerate obstacles , move to addObstacle?
|
|
|
|
// r = b - Ax
|
|
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++)
|
|
{
|
|
if (_obstacles[index] != EMPTY)
|
|
{
|
|
const int top = _obstacles[index + _slabSize];
|
|
const int bottom= _obstacles[index - _slabSize];
|
|
const int up = _obstacles[index + _xRes];
|
|
const int down = _obstacles[index - _xRes];
|
|
const int left = _obstacles[index - 1];
|
|
const int right = _obstacles[index + 1];
|
|
|
|
int counter = 0;
|
|
if (up) counter++;
|
|
if (down) counter++;
|
|
if (left) counter++;
|
|
if (right) counter++;
|
|
if (top) counter++;
|
|
if (bottom) counter++;
|
|
|
|
if (counter < 3)
|
|
_obstacles[index] = EMPTY;
|
|
}
|
|
if (_obstacles[index])
|
|
{
|
|
_xVelocity[index] =
|
|
_yVelocity[index] =
|
|
_zVelocity[index] = 0.0f;
|
|
_pressure[index] = 0.0f;
|
|
}
|
|
//vIndex++;
|
|
} // x-loop
|
|
//vIndex += 2;
|
|
} // y-loop
|
|
//vIndex += 2* _xRes;
|
|
} // z-loop
|
|
}
|
|
|
|
void FLUID_3D::floodFillComponent(int *buffer, size_t *queue, size_t limit, size_t pos, int from, int to)
|
|
{
|
|
/* Flood 'from' cells with 'to' in the grid. Rely on (from != 0 && from != to && edges == 0) to stop. */
|
|
int offsets[] = { -1, +1, -_xRes, +_xRes, -_slabSize, +_slabSize };
|
|
size_t qend = 0;
|
|
|
|
buffer[pos] = to;
|
|
queue[qend++] = pos;
|
|
|
|
for (size_t qidx = 0; qidx < qend; qidx++)
|
|
{
|
|
pos = queue[qidx];
|
|
|
|
for (int i = 0; i < 6; i++)
|
|
{
|
|
size_t next = pos + offsets[i];
|
|
|
|
if (next < limit && buffer[next] == from)
|
|
{
|
|
buffer[next] = to;
|
|
queue[qend++] = next;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
void FLUID_3D::mergeComponents(int *buffer, size_t *queue, size_t cur, size_t other)
|
|
{
|
|
/* Replace higher value with lower. */
|
|
if (buffer[other] < buffer[cur])
|
|
{
|
|
floodFillComponent(buffer, queue, cur, cur, buffer[cur], buffer[other]);
|
|
}
|
|
else if (buffer[cur] < buffer[other])
|
|
{
|
|
floodFillComponent(buffer, queue, cur, other, buffer[other], buffer[cur]);
|
|
}
|
|
}
|
|
|
|
void FLUID_3D::fixObstacleCompression(float *divergence)
|
|
{
|
|
int x, y, z;
|
|
size_t index;
|
|
|
|
/* Find compartments completely separated by obstacles.
|
|
* Edge of the domain is automatically component 0. */
|
|
int *component = new int[_totalCells];
|
|
size_t *queue = new size_t[_totalCells];
|
|
|
|
memset(component, 0, sizeof(int) * _totalCells);
|
|
|
|
int next_id = 1;
|
|
|
|
for (z = 1, index = _slabSize + _xRes + 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])
|
|
{
|
|
/* Check for connection to the domain edge at iteration end. */
|
|
if ((x == _xRes-2 && !_obstacles[index + 1]) ||
|
|
(y == _yRes-2 && !_obstacles[index + _xRes]) ||
|
|
(z == _zRes-2 && !_obstacles[index + _slabSize]))
|
|
{
|
|
component[index] = 0;
|
|
}
|
|
else {
|
|
component[index] = next_id;
|
|
}
|
|
|
|
if (!_obstacles[index - 1])
|
|
mergeComponents(component, queue, index, index - 1);
|
|
if (!_obstacles[index - _xRes])
|
|
mergeComponents(component, queue, index, index - _xRes);
|
|
if (!_obstacles[index - _slabSize])
|
|
mergeComponents(component, queue, index, index - _slabSize);
|
|
|
|
if (component[index] == next_id)
|
|
next_id++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
delete[] queue;
|
|
|
|
/* Compute average divergence within each component. */
|
|
float *total_divergence = new float[next_id];
|
|
int *component_size = new int[next_id];
|
|
|
|
memset(total_divergence, 0, sizeof(float) * next_id);
|
|
memset(component_size, 0, sizeof(int) * next_id);
|
|
|
|
for (z = 1, index = _slabSize + _xRes + 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])
|
|
{
|
|
int ci = component[index];
|
|
|
|
component_size[ci]++;
|
|
total_divergence[ci] += divergence[index];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/* Adjust divergence to make the average zero in each component except the edge. */
|
|
total_divergence[0] = 0.0f;
|
|
|
|
for (z = 1, index = _slabSize + _xRes + 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])
|
|
{
|
|
int ci = component[index];
|
|
|
|
divergence[index] -= total_divergence[ci] / component_size[ci];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
delete[] component;
|
|
delete[] component_size;
|
|
delete[] total_divergence;
|
|
}
|
|
|
|
//////////////////////////////////////////////////////////////////////
|
|
// add buoyancy forces
|
|
//////////////////////////////////////////////////////////////////////
|
|
void FLUID_3D::addBuoyancy(float *heat, float *density, float gravity[3], int zBegin, int zEnd)
|
|
{
|
|
int index = zBegin*_slabSize;
|
|
|
|
for (int z = zBegin; z < zEnd; z++)
|
|
for (int y = 0; y < _yRes; y++)
|
|
for (int x = 0; x < _xRes; x++, index++)
|
|
{
|
|
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;
|
|
}
|
|
}
|
|
|
|
|
|
//////////////////////////////////////////////////////////////////////
|
|
// add vorticity to the force field
|
|
//////////////////////////////////////////////////////////////////////
|
|
#define VORT_VEL(i, j) \
|
|
((_obstacles[obpos[(i)]] & 8) ? ((abs(objvelocity[(j)][obpos[(i)]]) > FLT_EPSILON) ? objvelocity[(j)][obpos[(i)]] : velocity[(j)][index]) : velocity[(j)][obpos[(i)]])
|
|
|
|
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+flame_vorticity<=0.0f) return;
|
|
|
|
int _blockSize=zEnd-zBegin;
|
|
int _blockTotalCells = _slabSize * (_blockSize+2);
|
|
|
|
float *_xVorticity, *_yVorticity, *_zVorticity, *_vorticity;
|
|
|
|
int bb=0;
|
|
int bt=0;
|
|
int bb1=-1;
|
|
int bt1=-1;
|
|
|
|
if (zBegin == 0) {bb1 = 1; bb = 1; _blockTotalCells-=_blockSize;}
|
|
if (zEnd == _zRes) {bt1 = 1;bt = 1; _blockTotalCells-=_blockSize;}
|
|
|
|
_xVorticity = new float[_blockTotalCells];
|
|
_yVorticity = new float[_blockTotalCells];
|
|
_zVorticity = new float[_blockTotalCells];
|
|
_vorticity = new float[_blockTotalCells];
|
|
|
|
memset(_xVorticity, 0, sizeof(float)*_blockTotalCells);
|
|
memset(_yVorticity, 0, sizeof(float)*_blockTotalCells);
|
|
memset(_zVorticity, 0, sizeof(float)*_blockTotalCells);
|
|
memset(_vorticity, 0, sizeof(float)*_blockTotalCells);
|
|
|
|
//const size_t indexsetupV=_slabSize;
|
|
const size_t index_ = _slabSize + _xRes + 1;
|
|
|
|
// calculate vorticity
|
|
float gridSize = 0.5f / _dx;
|
|
//index = _slabSize + _xRes + 1;
|
|
|
|
float *velocity[3];
|
|
float *objvelocity[3];
|
|
|
|
velocity[0] = _xVelocity;
|
|
velocity[1] = _yVelocity;
|
|
velocity[2] = _zVelocity;
|
|
|
|
objvelocity[0] = _xVelocityOb;
|
|
objvelocity[1] = _yVelocityOb;
|
|
objvelocity[2] = _zVelocityOb;
|
|
|
|
size_t vIndex=_xRes + 1;
|
|
for (int z = zBegin + bb1; z < (zEnd - bt1); z++)
|
|
{
|
|
size_t index = index_ +(z-1)*_slabSize;
|
|
vIndex = index-(zBegin-1+bb)*_slabSize;
|
|
|
|
for (int y = 1; y < _yRes - 1; y++, index += 2)
|
|
{
|
|
for (int x = 1; x < _xRes - 1; x++, index++)
|
|
{
|
|
if (!_obstacles[index])
|
|
{
|
|
int obpos[6];
|
|
|
|
obpos[0] = (_obstacles[index + _xRes] == 1) ? index : index + _xRes; // up
|
|
obpos[1] = (_obstacles[index - _xRes] == 1) ? index : index - _xRes; // down
|
|
float dy = (obpos[0] == index || obpos[1] == index) ? 1.0f / _dx : gridSize;
|
|
|
|
obpos[2] = (_obstacles[index + _slabSize] == 1) ? index : index + _slabSize; // out
|
|
obpos[3] = (_obstacles[index - _slabSize] == 1) ? index : index - _slabSize; // in
|
|
float dz = (obpos[2] == index || obpos[3] == index) ? 1.0f / _dx : gridSize;
|
|
|
|
obpos[4] = (_obstacles[index + 1] == 1) ? index : index + 1; // right
|
|
obpos[5] = (_obstacles[index - 1] == 1) ? index : index - 1; // left
|
|
float dx = (obpos[4] == index || obpos[5] == index) ? 1.0f / _dx : gridSize;
|
|
|
|
float xV[2], yV[2], zV[2];
|
|
|
|
zV[1] = VORT_VEL(0, 2);
|
|
zV[0] = VORT_VEL(1, 2);
|
|
yV[1] = VORT_VEL(2, 1);
|
|
yV[0] = VORT_VEL(3, 1);
|
|
_xVorticity[vIndex] = (zV[1] - zV[0]) * dy + (-yV[1] + yV[0]) * dz;
|
|
|
|
xV[1] = VORT_VEL(2, 0);
|
|
xV[0] = VORT_VEL(3, 0);
|
|
zV[1] = VORT_VEL(4, 2);
|
|
zV[0] = VORT_VEL(5, 2);
|
|
_yVorticity[vIndex] = (xV[1] - xV[0]) * dz + (-zV[1] + zV[0]) * dx;
|
|
|
|
yV[1] = VORT_VEL(4, 1);
|
|
yV[0] = VORT_VEL(5, 1);
|
|
xV[1] = VORT_VEL(0, 0);
|
|
xV[0] = VORT_VEL(1, 0);
|
|
_zVorticity[vIndex] = (yV[1] - yV[0]) * dx + (-xV[1] + xV[0])* dy;
|
|
|
|
_vorticity[vIndex] = sqrtf(_xVorticity[vIndex] * _xVorticity[vIndex] +
|
|
_yVorticity[vIndex] * _yVorticity[vIndex] +
|
|
_zVorticity[vIndex] * _zVorticity[vIndex]);
|
|
|
|
}
|
|
vIndex++;
|
|
}
|
|
vIndex+=2;
|
|
}
|
|
//vIndex+=2*_xRes;
|
|
}
|
|
|
|
// calculate normalized vorticity vectors
|
|
float eps = _vorticityEps;
|
|
|
|
//index = _slabSize + _xRes + 1;
|
|
vIndex=_slabSize + _xRes + 1;
|
|
|
|
for (int z = zBegin + bb; z < (zEnd - bt); z++)
|
|
{
|
|
size_t index = index_ +(z-1)*_slabSize;
|
|
vIndex = index-(zBegin-1+bb)*_slabSize;
|
|
|
|
for (int y = 1; y < _yRes - 1; y++, index += 2)
|
|
{
|
|
for (int x = 1; x < _xRes - 1; x++, index++)
|
|
{
|
|
//
|
|
|
|
if (!_obstacles[index])
|
|
{
|
|
float N[3];
|
|
|
|
int up = (_obstacles[index + _xRes] == 1) ? vIndex : vIndex + _xRes;
|
|
int down = (_obstacles[index - _xRes] == 1) ? vIndex : vIndex - _xRes;
|
|
float dy = (up == vIndex || down == vIndex) ? 1.0f / _dx : gridSize;
|
|
|
|
int out = (_obstacles[index + _slabSize] == 1) ? vIndex : vIndex + _slabSize;
|
|
int in = (_obstacles[index - _slabSize] == 1) ? vIndex : vIndex - _slabSize;
|
|
float dz = (out == vIndex || in == vIndex) ? 1.0f / _dx : gridSize;
|
|
|
|
int right = (_obstacles[index + 1] == 1) ? vIndex : vIndex + 1;
|
|
int left = (_obstacles[index - 1] == 1) ? vIndex : vIndex - 1;
|
|
float dx = (right == vIndex || left == vIndex) ? 1.0f / _dx : gridSize;
|
|
|
|
N[0] = (_vorticity[right] - _vorticity[left]) * dx;
|
|
N[1] = (_vorticity[up] - _vorticity[down]) * dy;
|
|
N[2] = (_vorticity[out] - _vorticity[in]) * dz;
|
|
|
|
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 + 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++;
|
|
} // x loop
|
|
vIndex+=2;
|
|
} // y loop
|
|
//vIndex+=2*_xRes;
|
|
} // z loop
|
|
|
|
if (_xVorticity) delete[] _xVorticity;
|
|
if (_yVorticity) delete[] _yVorticity;
|
|
if (_zVorticity) delete[] _zVorticity;
|
|
if (_vorticity) delete[] _vorticity;
|
|
}
|
|
|
|
|
|
void FLUID_3D::advectMacCormackBegin(int zBegin, int zEnd)
|
|
{
|
|
Vec3Int res = Vec3Int(_xRes,_yRes,_zRes);
|
|
|
|
setZeroX(_xVelocityOld, res, zBegin, zEnd);
|
|
setZeroY(_yVelocityOld, res, zBegin, zEnd);
|
|
setZeroZ(_zVelocityOld, res, zBegin, zEnd);
|
|
}
|
|
|
|
//////////////////////////////////////////////////////////////////////
|
|
// Advect using the MacCormack method from the Selle paper
|
|
//////////////////////////////////////////////////////////////////////
|
|
void FLUID_3D::advectMacCormackEnd1(int zBegin, int zEnd)
|
|
{
|
|
Vec3Int res = Vec3Int(_xRes,_yRes,_zRes);
|
|
|
|
const float dt0 = _dt / _dx;
|
|
|
|
int begin=zBegin * _slabSize;
|
|
int end=begin + (zEnd - zBegin) * _slabSize;
|
|
for (int x = begin; x < end; x++)
|
|
_xForce[x] = 0.0;
|
|
|
|
// advectFieldMacCormack1(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res)
|
|
|
|
advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _densityOld, _densityTemp, 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);
|
|
|
|
// Have to wait untill all the threads are done -> so continuing in step 3
|
|
}
|
|
|
|
//////////////////////////////////////////////////////////////////////
|
|
// Advect using the MacCormack method from the Selle paper
|
|
//////////////////////////////////////////////////////////////////////
|
|
void FLUID_3D::advectMacCormackEnd2(int zBegin, int zEnd)
|
|
{
|
|
const float dt0 = _dt / _dx;
|
|
Vec3Int res = Vec3Int(_xRes,_yRes,_zRes);
|
|
|
|
// 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);
|
|
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);
|
|
|
|
if(!_domainBcFront) copyBorderY(_yVelocityTemp, res, zBegin, zEnd);
|
|
else setZeroY(_yVelocityTemp, res, zBegin, zEnd);
|
|
|
|
if(!_domainBcTop) copyBorderZ(_zVelocityTemp, res, zBegin, zEnd);
|
|
else setZeroZ(_zVelocityTemp, res, zBegin, zEnd);
|
|
|
|
/* clear data boundaries */
|
|
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::processBurn(float *fuel, float *smoke, float *react, 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++)
|
|
{
|
|
float orig_fuel = fuel[index];
|
|
float orig_smoke = smoke[index];
|
|
float smoke_emit = 0.0f;
|
|
float flame = 0.0f;
|
|
|
|
/* process fuel */
|
|
fuel[index] -= burning_rate * dt;
|
|
if (fuel[index] < 0.0f) fuel[index] = 0.0f;
|
|
/* process reaction coordinate */
|
|
if (orig_fuel > FLT_EPSILON) {
|
|
react[index] *= fuel[index]/orig_fuel;
|
|
flame = pow(react[index], 0.5f);
|
|
}
|
|
else {
|
|
react[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);
|
|
|
|
/* set fluid temperature from the flame temperature profile */
|
|
if (heat && flame)
|
|
heat[index] = (1.0f - flame)*ignition_point + flame*temp_max;
|
|
|
|
/* mix new color */
|
|
if (r && smoke_emit > FLT_EPSILON) {
|
|
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;
|
|
}
|
|
}
|
|
}
|
|
|
|
void FLUID_3D::updateFlame(float *react, float *flame, int total_cells)
|
|
{
|
|
for (int index = 0; index < total_cells; index++)
|
|
{
|
|
/* model flame temperature curve from the reaction coordinate (fuel)
|
|
* TODO: Would probably be best to get rid of whole "flame" data field.
|
|
* Currently it's just sqrt mirror of reaction coordinate, and therefore
|
|
* basically just waste of memory and disk space...
|
|
*/
|
|
if (react[index]>0.0f) {
|
|
/* do a smooth falloff for rest of the values */
|
|
flame[index] = pow(react[index], 0.5f);
|
|
}
|
|
else
|
|
flame[index] = 0.0f;
|
|
}
|
|
}
|