Fix for smoke boundary conditions: "Open borders" behaved differently for positive and negative axises.

+ Some code cleanup.
This commit is contained in:
Miika Hamalainen 2012-05-26 21:36:19 +00:00
parent ebdeed07e5
commit eda0f3b186
2 changed files with 56 additions and 116 deletions

@ -153,11 +153,11 @@ void FLUID_3D::setBorderObstacles()
{ {
// bottom slab // bottom slab
index = x + y * _xRes; index = x + y * _xRes;
if(_domainBcBottom==1) _obstacles[index] = 1; if(_domainBcBottom) _obstacles[index] = 1;
// top slab // top slab
index += _totalCells - _slabSize; index += _totalCells - _slabSize;
if(_domainBcTop==1) _obstacles[index] = 1; if(_domainBcTop) _obstacles[index] = 1;
} }
for (int z = 0; z < _zRes; z++) for (int z = 0; z < _zRes; z++)
@ -165,11 +165,11 @@ void FLUID_3D::setBorderObstacles()
{ {
// front slab // front slab
index = x + z * _slabSize; index = x + z * _slabSize;
if(_domainBcFront==1) _obstacles[index] = 1; if(_domainBcFront) _obstacles[index] = 1;
// back slab // back slab
index += _slabSize - _xRes; index += _slabSize - _xRes;
if(_domainBcBack==1) _obstacles[index] = 1; if(_domainBcBack) _obstacles[index] = 1;
} }
for (int z = 0; z < _zRes; z++) for (int z = 0; z < _zRes; z++)
@ -177,11 +177,11 @@ void FLUID_3D::setBorderObstacles()
{ {
// left slab // left slab
index = y * _xRes + z * _slabSize; index = y * _xRes + z * _slabSize;
if(_domainBcLeft==1) _obstacles[index] = 1; if(_domainBcLeft) _obstacles[index] = 1;
// right slab // right slab
index += _xRes - 1; index += _xRes - 1;
if(_domainBcRight==1) _obstacles[index] = 1; if(_domainBcRight) _obstacles[index] = 1;
} }
} }
@ -449,48 +449,7 @@ void FLUID_3D::setBorderCollisions() {
// set side obstacles // set side obstacles
int index; setBorderObstacles();
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;
}
} }
////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////
@ -813,13 +772,13 @@ void FLUID_3D::project()
setObstacleBoundaries(_pressure, 0, _zRes); setObstacleBoundaries(_pressure, 0, _zRes);
// copy out the boundaries // copy out the boundaries
if(_domainBcLeft == 0) setNeumannX(_xVelocity, _res, 0, _zRes); if(!_domainBcLeft) setNeumannX(_xVelocity, _res, 0, _zRes);
else setZeroX(_xVelocity, _res, 0, _zRes); else setZeroX(_xVelocity, _res, 0, _zRes);
if(_domainBcFront == 0) setNeumannY(_yVelocity, _res, 0, _zRes); if(!_domainBcFront) setNeumannY(_yVelocity, _res, 0, _zRes);
else setZeroY(_yVelocity, _res, 0, _zRes); else setZeroY(_yVelocity, _res, 0, _zRes);
if(_domainBcTop == 0) setNeumannZ(_zVelocity, _res, 0, _zRes); if(!_domainBcTop) setNeumannZ(_zVelocity, _res, 0, _zRes);
else setZeroZ(_zVelocity, _res, 0, _zRes); else setZeroZ(_zVelocity, _res, 0, _zRes);
/* /*
@ -1369,13 +1328,13 @@ void FLUID_3D::advectMacCormackBegin(int zBegin, int zEnd)
{ {
Vec3Int res = Vec3Int(_xRes,_yRes,_zRes); Vec3Int res = Vec3Int(_xRes,_yRes,_zRes);
if(_domainBcLeft == 0) copyBorderX(_xVelocityOld, res, zBegin, zEnd); if(!_domainBcLeft) copyBorderX(_xVelocityOld, res, zBegin, zEnd);
else setZeroX(_xVelocityOld, res, zBegin, zEnd); else setZeroX(_xVelocityOld, res, zBegin, zEnd);
if(_domainBcFront == 0) copyBorderY(_yVelocityOld, res, zBegin, zEnd); if(!_domainBcFront) copyBorderY(_yVelocityOld, res, zBegin, zEnd);
else setZeroY(_yVelocityOld, res, zBegin, zEnd); else setZeroY(_yVelocityOld, res, zBegin, zEnd);
if(_domainBcTop == 0) copyBorderZ(_zVelocityOld, res, zBegin, zEnd); if(!_domainBcTop) copyBorderZ(_zVelocityOld, res, zBegin, zEnd);
else setZeroZ(_zVelocityOld, res, zBegin, zEnd); else setZeroZ(_zVelocityOld, res, zBegin, zEnd);
} }
@ -1423,13 +1382,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, _yVelocityOld, _yVelocityTemp, _yVelocity, t1, res, _obstacles, zBegin, zEnd);
advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _zVelocityOld, _zVelocityTemp, _zVelocity, t1, res, _obstacles, zBegin, zEnd); advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _zVelocityOld, _zVelocityTemp, _zVelocity, t1, res, _obstacles, zBegin, zEnd);
if(_domainBcLeft == 0) copyBorderX(_xVelocityTemp, res, zBegin, zEnd); if(!_domainBcLeft) copyBorderX(_xVelocityTemp, res, zBegin, zEnd);
else setZeroX(_xVelocityTemp, res, zBegin, zEnd); else setZeroX(_xVelocityTemp, res, zBegin, zEnd);
if(_domainBcFront == 0) copyBorderY(_yVelocityTemp, res, zBegin, zEnd); if(!_domainBcFront) copyBorderY(_yVelocityTemp, res, zBegin, zEnd);
else setZeroY(_yVelocityTemp, res, zBegin, zEnd); else setZeroY(_yVelocityTemp, res, zBegin, zEnd);
if(_domainBcTop == 0) copyBorderZ(_zVelocityTemp, res, zBegin, zEnd); if(!_domainBcTop) copyBorderZ(_zVelocityTemp, res, zBegin, zEnd);
else setZeroZ(_zVelocityTemp, res, zBegin, zEnd); else setZeroZ(_zVelocityTemp, res, zBegin, zEnd);
setZeroBorder(_density, res, zBegin, zEnd); setZeroBorder(_density, res, zBegin, zEnd);

@ -92,19 +92,15 @@ void FLUID_3D::setNeumannX(float* field, Vec3Int res, int zBegin, int zEnd)
// left slab // left slab
index = y * res[0] + z * slabSize; index = y * res[0] + z * slabSize;
field[index] = field[index + 2]; field[index] = field[index + 2];
/* only allow outwards flux */
if(field[index]>0.) field[index] = 0.;
index += 1;
if(field[index]>0.) field[index] = 0.;
// right slab // right slab
index += res[0] - 1; index = y * res[0] + z * slabSize + res[0] - 1;
field[index] = field[index - 2]; field[index] = field[index - 2];
} /* only allow outwards flux */
// fix, force top slab to only allow outwards flux
for (int y = 0; y < res[1]; y++)
for (int z = zBegin; z < zEnd; z++)
{
// top slab
index = y * res[0] + z * slabSize;
index += res[0] - 1;
if(field[index]<0.) field[index] = 0.; if(field[index]<0.) field[index] = 0.;
index -= 1; index -= 1;
if(field[index]<0.) field[index] = 0.; if(field[index]<0.) field[index] = 0.;
@ -121,27 +117,22 @@ void FLUID_3D::setNeumannY(float* field, Vec3Int res, int zBegin, int zEnd)
for (int z = zBegin; z < zEnd; z++) for (int z = zBegin; z < zEnd; z++)
for (int x = 0; x < res[0]; x++) for (int x = 0; x < res[0]; x++)
{ {
// bottom slab // front slab
index = x + z * slabSize; index = x + z * slabSize;
field[index] = field[index + 2 * res[0]]; field[index] = field[index + 2 * res[0]];
/* only allow outwards flux */
if(field[index]>0.) field[index] = 0.;
index += res[0];
if(field[index]>0.) field[index] = 0.;
// top slab // back slab
index += slabSize - res[0]; index = x + z * slabSize + slabSize - res[0];
field[index] = field[index - 2 * res[0]]; field[index] = field[index - 2 * res[0]];
} /* only allow outwards flux */
// fix, force top slab to only allow outwards flux
for (int z = zBegin; z < zEnd; z++)
for (int x = 0; x < res[0]; x++)
{
// top slab
index = x + z * slabSize;
index += slabSize - res[0];
if(field[index]<0.) field[index] = 0.; if(field[index]<0.) field[index] = 0.;
index -= res[0]; index -= res[0];
if(field[index]<0.) field[index] = 0.; if(field[index]<0.) field[index] = 0.;
} }
} }
////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////
@ -154,43 +145,33 @@ void FLUID_3D::setNeumannZ(float* field, Vec3Int res, int zBegin, int zEnd)
const int cellsslab = totalCells - slabSize; const int cellsslab = totalCells - slabSize;
int index; int index;
index = 0; if (zBegin == 0) {
if (zBegin == 0) for (int y = 0; y < res[1]; y++)
for (int y = 0; y < res[1]; y++) for (int x = 0; x < res[0]; x++)
for (int x = 0; x < res[0]; x++, index++) {
{ // front slab
// front slab index = x + y * res[0];
field[index] = field[index + 2 * slabSize]; field[index] = field[index + 2 * slabSize];
} /* only allow outwards flux */
if(field[index]>0.) field[index] = 0.;
index += slabSize;
if(field[index]>0.) field[index] = 0.;
}
}
if (zEnd == res[2]) if (zEnd == res[2]) {
{ for (int y = 0; y < res[1]; y++)
index = 0; for (int x = 0; x < res[0]; x++)
int indexx = 0; {
// back slab
for (int y = 0; y < res[1]; y++) index = x + y * res[0] + cellsslab;
for (int x = 0; x < res[0]; x++, index++) field[index] = field[index - 2 * slabSize];
{ /* only allow outwards flux */
if(field[index]<0.) field[index] = 0.;
// back slab index -= slabSize;
indexx = index + cellsslab; if(field[index]<0.) field[index] = 0.;
field[indexx] = field[indexx - 2 * slabSize]; }
} }
// fix, force top slab to only allow outwards flux
for (int y = 0; y < res[1]; y++)
for (int x = 0; x < res[0]; x++)
{
// top slab
index = x + y * res[0];
index += cellsslab;
if(field[index]<0.) field[index] = 0.;
index -= slabSize;
if(field[index]<0.) field[index] = 0.;
}
} // zEnd == res[2]
} }