smoke: decoupling high-res, introducing reset for heat+gravity, fixing another obstacle problem

This commit is contained in:
Daniel Genrich 2009-08-06 16:50:02 +00:00
parent 65cf500fdc
commit 2fd0263602
9 changed files with 194 additions and 107 deletions

@ -59,7 +59,10 @@ FLUID_3D::FLUID_3D(int *res, int amplify, float *p0, float dt) :
_maxRes = MAX3(_xRes, _yRes, _zRes);
// initialize wavelet turbulence
_wTurbulence = new WTURBULENCE(_res[0],_res[1],_res[2], amplify);
if(amplify)
_wTurbulence = new WTURBULENCE(_res[0],_res[1],_res[2], amplify);
else
_wTurbulence = NULL;
// scale the constants according to the refinement of the grid
_dx = 1.0f / (float)_maxRes;
@ -215,11 +218,11 @@ void FLUID_3D::step()
wipeBoundaries();
// run the solvers
addVorticity();
addBuoyancy(_heat, _density);
addVorticity();
addBuoyancy(_heat, _density);
addForce();
project();
diffuseHeat();
diffuseHeat();
// advect everything
advectMacCormack();
@ -249,6 +252,19 @@ void FLUID_3D::step()
*/
_totalTime += _dt;
_totalSteps++;
// clear obstacles
for (int i = 0; i < _totalCells; i++)
{
if(_obstacles[i])
{
_xVelocity[i] =
_yVelocity[i] =
_zVelocity[i] = 0.0f;
_pressure[i] = 0.0f;
}
_obstacles[i] = 0;
}
}
//////////////////////////////////////////////////////////////////////
@ -354,6 +370,7 @@ void FLUID_3D::addForce()
void FLUID_3D::project()
{
int index, x, y, z;
setObstacleBoundaries();
// copy out the boundaries
@ -397,6 +414,8 @@ void FLUID_3D::project()
// solve Poisson equation
solvePressure(_pressure, _divergence, _obstacles);
setObstaclePressure();
// project out solution
float invDx = 1.0f / _dx;
index = _slabSize + _xRes + 1;
@ -404,9 +423,12 @@ void FLUID_3D::project()
for (y = 1; y < _yRes - 1; y++, index += 2)
for (x = 1; x < _xRes - 1; x++, index++)
{
_xVelocity[index] -= 0.5f * (_pressure[index + 1] - _pressure[index - 1]) * invDx;
_yVelocity[index] -= 0.5f * (_pressure[index + _xRes] - _pressure[index - _xRes]) * invDx;
_zVelocity[index] -= 0.5f * (_pressure[index + _slabSize] - _pressure[index - _slabSize]) * invDx;
if(!_obstacles[index])
{
_xVelocity[index] -= 0.5f * (_pressure[index + 1] - _pressure[index - 1]) * invDx;
_yVelocity[index] -= 0.5f * (_pressure[index + _xRes] - _pressure[index - _xRes]) * invDx;
_zVelocity[index] -= 0.5f * (_pressure[index + _slabSize] - _pressure[index - _slabSize]) * invDx;
}
}
}
@ -443,6 +465,80 @@ void FLUID_3D::addObstacle(OBSTACLE* obstacle)
//////////////////////////////////////////////////////////////////////
// calculate the obstacle directional types
//////////////////////////////////////////////////////////////////////
void FLUID_3D::setObstaclePressure()
{
// tag remaining obstacle blocks
for (int z = 1, index = _slabSize + _xRes + 1;
z < _zRes - 1; z++, index += 2 * _xRes)
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])
{
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., vp = 0.;
if (left && !right) {
_pressure[index] += _pressure[index + 1];
pcnt += 1.;
}
if (!left && right) {
_pressure[index] += _pressure[index - 1];
pcnt += 1.;
}
if (up && !down) {
_pressure[index] += _pressure[index - _xRes];
pcnt += 1.;
}
if (!up && down) {
_pressure[index] += _pressure[index + _xRes];
pcnt += 1.;
}
if (top && !bottom) {
_pressure[index] += _pressure[index - _slabSize];
pcnt += 1.;
// _zVelocity[index] += - _zVelocity[index - _slabSize];
// vp += 1.0;
}
if (!top && bottom) {
_pressure[index] += _pressure[index + _slabSize];
pcnt += 1.;
// _zVelocity[index] += - _zVelocity[index + _slabSize];
// vp += 1.0;
}
if(pcnt > 0.000001f)
_pressure[index] /= pcnt;
// test - dg
// if(vp > 0.000001f)
// _zVelocity[index] /= vp;
// 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
}
}
}
void FLUID_3D::setObstacleBoundaries()
{
// cull degenerate obstacles , move to addObstacle?
@ -450,6 +546,7 @@ void FLUID_3D::setObstacleBoundaries()
z < _zRes - 1; z++, index += 2 * _xRes)
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];
@ -470,67 +567,14 @@ void FLUID_3D::setObstacleBoundaries()
if (counter < 3)
_obstacles[index] = EMPTY;
}
// tag remaining obstacle blocks
for (int z = 1, index = _slabSize + _xRes + 1;
z < _zRes - 1; z++, index += 2 * _xRes)
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] != 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];
// 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.;
if (_obstacles[index])
{
_xVelocity[index] =
_yVelocity[index] =
_zVelocity[index] = 0.0f;
_pressure[index] = 0.0f;
}
if (!left && right) {
_pressure[index] += _pressure[index - 1];
pcnt += 1.;
}
if (up && !down) {
_pressure[index] += _pressure[index - _xRes];
pcnt += 1.;
}
if (!up && down) {
_pressure[index] += _pressure[index + _xRes];
pcnt += 1.;
}
if (top && !bottom) {
_pressure[index] += _pressure[index - _xRes];
pcnt += 1.;
}
if (!top && bottom) {
_pressure[index] += _pressure[index + _xRes];
pcnt += 1.;
}
_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
}
}
}
//////////////////////////////////////////////////////////////////////

@ -132,6 +132,7 @@ class FLUID_3D
// handle obstacle boundaries
void setObstacleBoundaries();
void setObstaclePressure();
public:
// advection, accessed e.g. by WTURBULENCE class

@ -80,7 +80,7 @@ void FLUID_3D::addSmokeTestCase(float* field, Vec3Int res, float value)
void FLUID_3D::setNeumannX(float* field, Vec3Int res)
{
const int slabSize = res[0] * res[1];
int index;
size_t index;
for (int z = 0; z < res[2]; z++)
for (int y = 0; y < res[1]; y++)
{
@ -100,7 +100,7 @@ void FLUID_3D::setNeumannX(float* field, Vec3Int res)
void FLUID_3D::setNeumannY(float* field, Vec3Int res)
{
const int slabSize = res[0] * res[1];
int index;
size_t index;
for (int z = 0; z < res[2]; z++)
for (int x = 0; x < res[0]; x++)
{
@ -121,7 +121,7 @@ void FLUID_3D::setNeumannZ(float* field, Vec3Int res)
{
const int slabSize = res[0] * res[1];
const int totalCells = res[0] * res[1] * res[2];
int index;
size_t index;
for (int y = 0; y < res[1]; y++)
for (int x = 0; x < res[0]; x++)
{
@ -141,10 +141,11 @@ void FLUID_3D::setNeumannZ(float* field, Vec3Int res)
// top slab
index = x + y * res[0];
index += totalCells - slabSize;
if(field[index]<0.) field[index] = 0.;
if(field[index]<0.) field[index] = 0.0f;
index -= slabSize;
if(field[index]<0.) field[index] = 0.;
if(field[index]<0.) field[index] = 0.0f;
}
}
//////////////////////////////////////////////////////////////////////

@ -88,15 +88,18 @@ extern "C" float *smoke_get_velocity_z(FLUID_3D *fluid)
extern "C" float *smoke_get_bigdensity(FLUID_3D *fluid)
{
return fluid->_wTurbulence->getDensityBig();
return fluid->_wTurbulence ? fluid->_wTurbulence->getDensityBig() : NULL;
}
extern "C" void smoke_get_bigres(FLUID_3D *fluid, int *res)
{
Vec3Int r = fluid->_wTurbulence->getResBig();
res[0] = r[0];
res[1] = r[1];
res[2] = r[2];
if(fluid->_wTurbulence)
{
Vec3Int r = fluid->_wTurbulence->getResBig();
res[0] = r[0];
res[1] = r[1];
res[2] = r[2];
}
}
extern "C" unsigned char *smoke_get_obstacle(FLUID_3D *fluid)

@ -59,7 +59,10 @@ class PHYSICS_PT_smoke(PhysicButtonsPanel):
sub = split.column()
sub.itemL(text="Display:")
sub.itemR(md.domain_settings, "visibility")
sub.itemR(md.domain_settings, "color")
sub.itemR(md.domain_settings, "color", slider=True)
mysub = sub.column()
mysub.active = md.domain_settings.highres
mysub.itemR(md.domain_settings, "viewhighres")
layout.itemL(text="Noise Type:")
layout.itemR(md.domain_settings, "noise_type", expand=True)
@ -79,21 +82,20 @@ class PHYSICS_PT_smoke(PhysicButtonsPanel):
elif md.smoke_type == 'TYPE_FLOW':
layout.itemR(md.flow_settings, "outflow")
split = layout.split()
col = split.column()
col.itemR(md.flow_settings, "outflow")
col.itemL(text="Particle System:")
col.item_pointerR(md.flow_settings, "psys", ob, "particle_systems", text="")
if md.flow_settings.outflow:
col = split.column()
else:
col = split.column()
col.itemL(text="Behavior:")
col.itemR(md.flow_settings, "temperature")
col.itemR(md.flow_settings, "density")
sub = split.column()
sub.itemL(text="particle System:")
sub.item_pointerR(md.flow_settings, "psys", ob, "particle_systems", text="")
sub.itemL(text="Behavior:")
sub.itemR(md.flow_settings, "temperature")
sub.itemR(md.flow_settings, "density")
elif md.smoke_type == 'TYPE_COLL':
layout.itemS()

@ -283,7 +283,7 @@ int smokeModifier_init (SmokeModifierData *smd, Object *ob, Scene *scene, Derive
// printf("res[0]: %d, res[1]: %d, res[2]: %d\n", smd->domain->res[0], smd->domain->res[1], smd->domain->res[2]);
// dt max is 0.1
smd->domain->fluid = smoke_init(smd->domain->res, smd->domain->amplify, smd->domain->p0, smd->domain->p1, 2.5 / FPS);
smd->domain->fluid = smoke_init(smd->domain->res, (smd->domain->flags & MOD_SMOKE_HIGHRES) ? (smd->domain->amplify + 1) : 0, smd->domain->p0, smd->domain->p1, 2.5 / FPS);
smd->time = scene->r.cfra;
smd->domain->firstframe = smd->time;
@ -324,6 +324,10 @@ int smokeModifier_init (SmokeModifierData *smd, Object *ob, Scene *scene, Derive
size_t max_points = 0;
size_t quads = 0, facecounter = 0;
// copy obmat
Mat4CpyMat4(scs->mat, ob->obmat);
Mat4CpyMat4(scs->mat_old, ob->obmat);
// count quads
for(i = 0; i < dm->getNumFaces(dm); i++)
{
@ -451,7 +455,7 @@ int smokeModifier_init (SmokeModifierData *smd, Object *ob, Scene *scene, Derive
if(!smd->coll->bvhtree)
{
smd->coll->bvhtree = bvhtree_build_from_smoke ( ob->obmat, dm->getFaceArray(dm), dm->getNumFaces(dm), dm->getVertArray(dm), dm->getNumVerts(dm), 0.0 );
smd->coll->bvhtree = NULL; // bvhtree_build_from_smoke ( ob->obmat, dm->getFaceArray(dm), dm->getNumFaces(dm), dm->getVertArray(dm), dm->getNumVerts(dm), 0.0 );
}
}
@ -504,12 +508,12 @@ void calcTriangleDivs(Object *ob, MVert *verts, int numverts, MFace *faces, int
if(INPR(side1, side1) > fsTri*fsTri)
{
float tmp = Normalize(side1);
divs1 = (int)(tmp/fsTri);
divs1 = (int)ceil(tmp/fsTri);
}
if(INPR(side2, side2) > fsTri*fsTri)
{
float tmp = Normalize(side2);
divs2 = (int)(tmp/fsTri);
divs2 = (int)ceil(tmp/fsTri);
/*
// debug
@ -529,9 +533,11 @@ void calcTriangleDivs(Object *ob, MVert *verts, int numverts, MFace *faces, int
facecounter++;
VECCOPY(p1, verts[faces[i].v3].co);
VECCOPY(p0, verts[faces[i].v3].co);
Mat4MulVecfl (ob->obmat, p0);
VECCOPY(p1, verts[faces[i].v4].co);
Mat4MulVecfl (ob->obmat, p1);
VECCOPY(p2, verts[faces[i].v4].co);
VECCOPY(p2, verts[faces[i].v1].co);
Mat4MulVecfl (ob->obmat, p2);
VECSUB(side1, p1, p0);
@ -541,12 +547,12 @@ void calcTriangleDivs(Object *ob, MVert *verts, int numverts, MFace *faces, int
if(INPR(side1, side1) > fsTri*fsTri)
{
float tmp = Normalize(side1);
divs1 = (int)(tmp/fsTri);
divs1 = (int)ceil(tmp/fsTri);
}
if(INPR(side2, side2) > fsTri*fsTri)
{
float tmp = Normalize(side2);
divs2 = (int)(tmp/fsTri);
divs2 = (int)ceil(tmp/fsTri);
}
(*tridivs)[3 * facecounter + 0] = divs1;
@ -612,6 +618,10 @@ void smokeModifier_freeCollision(SmokeModifierData *smd)
smd->coll->bvhtree = NULL;
}
if(smd->coll->dm)
smd->coll->dm->release(smd->coll->dm);
smd->coll->dm = NULL;
MEM_freeN(smd->coll);
smd->coll = NULL;
}
@ -670,6 +680,11 @@ void smokeModifier_reset(struct SmokeModifierData *smd)
BLI_bvhtree_free(smd->coll->bvhtree);
smd->coll->bvhtree = NULL;
}
if(smd->coll->dm)
smd->coll->dm->release(smd->coll->dm);
smd->coll->dm = NULL;
}
}
}
@ -703,7 +718,7 @@ void smokeModifier_createType(struct SmokeModifierData *smd)
smd->domain->fluid_group = NULL;
smd->domain->coll_group = NULL;
smd->domain->maxres = 32;
smd->domain->amplify = 2;
smd->domain->amplify = 1;
smd->domain->omega = 1.0;
smd->domain->alpha = -0.001;
smd->domain->beta = 0.1;
@ -747,6 +762,7 @@ void smokeModifier_createType(struct SmokeModifierData *smd)
smd->coll->points = NULL;
smd->coll->numpoints = 0;
smd->coll->bvhtree = NULL;
smd->coll->dm = NULL;
}
}
}
@ -783,6 +799,15 @@ void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedM
bvhtree_update_from_smoke ( ob->obmat, smd->coll->bvhtree, dm->getFaceArray(dm), dm->getNumFaces(dm), dm->getVertArray(dm), dm->getNumVerts(dm));
else
printf("smoke coll with no bvh\n");
if(smd->coll->dm)
smd->coll->dm->release(smd->coll->dm);
smd->coll->dm = CDDM_copy(dm);
// rigid movement support
Mat4CpyMat4(smd->coll->mat_old, smd->coll->mat);
Mat4CpyMat4(smd->coll->mat, ob->obmat);
}
else if(scene->r.cfra < smd->time)
{
@ -883,9 +908,7 @@ void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedM
float *velocity_x = smoke_get_velocity_x(sds->fluid);
float *velocity_y = smoke_get_velocity_y(sds->fluid);
float *velocity_z = smoke_get_velocity_z(sds->fluid);
int bigres[3];
smoke_get_bigres(smd->domain->fluid, bigres);
int bigres[3];
// mostly copied from particle code
for(p=0, pa=psys->particles; p<psys->totpart; p++, pa++)
@ -926,15 +949,20 @@ void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedM
{
heat[index] = sfs->temp;
density[index] = sfs->density;
/*
velocity_x[index] = pa->state.vel[0];
velocity_y[index] = pa->state.vel[1];
velocity_z[index] = pa->state.vel[2];
*/
// we need different handling for the high-res feature
if(bigdensity)
{
// init all surrounding cells according to amplification, too
int i, j, k;
smoke_get_bigres(smd->domain->fluid, bigres);
for(i = 0; i < smd->domain->amplify; i++)
for(j = 0; j < smd->domain->amplify; j++)
for(k = 0; k < smd->domain->amplify; k++)
@ -957,6 +985,9 @@ void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedM
{
// init all surrounding cells according to amplification, too
int i, j, k;
smoke_get_bigres(smd->domain->fluid, bigres);
for(i = 0; i < smd->domain->amplify; i++)
for(j = 0; j < smd->domain->amplify; j++)
for(k = 0; k < smd->domain->amplify; k++)
@ -1039,16 +1070,14 @@ void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedM
{
// we got nice collision object
SmokeCollSettings *scs = smd2->coll;
int cell[3];
size_t index = 0;
size_t i, j;
unsigned char *obstacles = smoke_get_obstacle(smd->domain->fluid);
// int BLI_bvhtree_find_nearest(BVHTree *tree, const float *co, BVHTreeNearest *nearest, BVHTree_NearestPointCallback callback, void *userdata);
for(i = 0; i < scs->numpoints; i++)
{
int badcell = 0;
size_t index = 0;
int cell[3];
// 1. get corresponding cell
get_cell(smd, &scs->points[3 * i], cell, 0);

@ -5349,7 +5349,7 @@ void draw_object(Scene *scene, ARegion *ar, View3D *v3d, Base *base, int flag)
float *buffer = NULL;
int res[3];
float bigfactor = 1.0;
int big = smd->domain->flags & MOD_SMOKE_HIGHRES;
int big = (smd->domain->flags & MOD_SMOKE_HIGHRES) && (smd->domain->viewsettings & MOD_SMOKE_VIEW_USEBIG);
int new = 0;
// GUI sent redraw event

@ -44,6 +44,7 @@
#define MOD_SMOKE_VIEW_CHANGETOBIG (1<<5)
#define MOD_SMOKE_VIEW_REDRAWNICE (1<<6)
#define MOD_SMOKE_VIEW_REDRAWALL (1<<7)
#define MOD_SMOKE_VIEW_USEBIG (1<<8)
typedef struct SmokeDomainSettings {
struct SmokeModifierData *smd; /* for fast RNA access */
@ -103,6 +104,7 @@ typedef struct SmokeFlowSettings {
typedef struct SmokeCollSettings {
struct SmokeModifierData *smd; /* for fast RNA access */
struct BVHTree *bvhtree; /* bounding volume hierarchy for this cloth object */
struct DerivedMesh *dm;
float *points;
float *points_old;
float *vel;

@ -150,8 +150,13 @@ static void rna_def_smoke_domain_settings(BlenderRNA *brna)
prop= RNA_def_property(srna, "highres", PROP_BOOLEAN, PROP_NONE);
RNA_def_property_boolean_sdna(prop, NULL, "flags", MOD_SMOKE_HIGHRES);
RNA_def_property_ui_text(prop, "High res", "Show high resolution (using amplification).");
RNA_def_property_update(prop, NC_OBJECT|ND_DRAW, NULL);
RNA_def_property_ui_text(prop, "High res", "Enable high resolution (using amplification).");
RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset");
prop= RNA_def_property(srna, "viewhighres", PROP_BOOLEAN, PROP_NONE);
RNA_def_property_boolean_sdna(prop, NULL, "viewsettings", MOD_SMOKE_VIEW_USEBIG);
RNA_def_property_ui_text(prop, "Show high res", "Show high resolution (using amplification).");
RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, NULL);
prop= RNA_def_property(srna, "noise_type", PROP_ENUM, PROP_NONE);
RNA_def_property_enum_sdna(prop, NULL, "noise");
@ -171,14 +176,14 @@ static void rna_def_smoke_domain_settings(BlenderRNA *brna)
RNA_def_property_range(prop, -5.0, 5.0);
RNA_def_property_ui_range(prop, -5.0, 5.0, 0.02, 5);
RNA_def_property_ui_text(prop, "Gravity", "Higher value results in sinking smoke");
RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, NULL);
RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset");
prop= RNA_def_property(srna, "beta", PROP_FLOAT, PROP_NONE);
RNA_def_property_float_sdna(prop, NULL, "beta");
RNA_def_property_range(prop, -5.0, 5.0);
RNA_def_property_ui_range(prop, -5.0, 5.0, 0.02, 5);
RNA_def_property_ui_text(prop, "Heat", "Higher value results in faster rising smoke.");
RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, NULL);
RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset");
prop= RNA_def_property(srna, "coll_group", PROP_POINTER, PROP_NONE);
RNA_def_property_pointer_sdna(prop, NULL, "coll_group");