FLUID_3D.cpp

Go to the documentation of this file.

// 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>

#if PARALLEL==1
#include <omp.h>
#endif // PARALLEL 

// Construction/Destruction

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_DEFAULT;   // just in case. set in step from a RNA factor
    
    // start point of array
    _p0[0] = p0[0];
    _p0[1] = p0[1];
    _p0[2] = p0[2];

    _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
    _dx = 1.0f / (float)_maxRes;
    _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];
    _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];
    _heat         = new float[_totalCells];
    _heatOld      = 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];
    _heatTemp      = 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;
        _heat[x]         = 0.0f;
        _heatOld[x]      = 0.0f;
        _xVelocity[x]    = 0.0f;
        _yVelocity[x]    = 0.0f;
        _zVelocity[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;
    }

    // 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++)
    {
        // bottom slab
        index = x + y * _xRes;
        if(_domainBcBottom==1) _obstacles[index] = 1;

        // top slab
        index += _totalCells - _slabSize;
        if(_domainBcTop==1) _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==1) _obstacles[index] = 1;

        // back slab
        index += _slabSize - _xRes;
        if(_domainBcBack==1) _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==1) _obstacles[index] = 1;

        // right slab
        index += _xRes - 1;
        if(_domainBcRight==1) _obstacles[index] = 1;
    }

}

FLUID_3D::~FLUID_3D()
{
    if (_xVelocity) delete[] _xVelocity;
    if (_yVelocity) delete[] _yVelocity;
    if (_zVelocity) delete[] _zVelocity;
    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 (_wTurbulence) delete _wTurbulence;

    if (_xVelocityTemp) delete[] _xVelocityTemp;
    if (_yVelocityTemp) delete[] _yVelocityTemp;
    if (_zVelocityTemp) delete[] _zVelocityTemp;
    if (_densityTemp) delete[] _densityTemp;
    if (_heatTemp) delete[] _heatTemp;

    // 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)
{
    _alpha = alpha;
    _beta = beta;
    _dtFactor = dt_factor;
    _vorticityRNA = vorticity;
    _borderColli = borderCollision;
}

// step simulation once
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();

    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


#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

        wipeBoundariesSL(zBegin, zEnd);
        addVorticity(zBegin, zEnd);
        addBuoyancy(_heat, _density, 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
        {
#endif
        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);

    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
    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...


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./6. * 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./6. * 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./6. * 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./6. * 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./6. * 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./6. * 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./6. * 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./6. * 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./6. * 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./6. * 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./6. * 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./6. * 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);
}

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;

            // right slab
            index += _xRes - 1;
            _xVelocity[index] = 0.0f;
            _yVelocity[index] = 0.0f;
            _zVelocity[index] = 0.0f;
            _density[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;

            // top slab
            index += slabSize - _xRes;
            _xVelocity[index] = 0.0f;
            _yVelocity[index] = 0.0f;
            _zVelocity[index] = 0.0f;
            _density[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 (zEnd == _zRes)
    {
        index=0;
        int indexx=0;
        const int cellsslab = totalCells - slabSize;

        for (y = 0; y < _yRes; y++)
            for (x = 0; x < _xRes; x++, index++)
            {

                // back slab
                indexx = index + cellsslab;
                _xVelocity[indexx] = 0.0f;
                _yVelocity[indexx] = 0.0f;
                _zVelocity[indexx] = 0.0f;
                _density[indexx] = 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);
    
    setObstacleBoundaries(_pressure, 0, _zRes);
    
    // copy out the boundaries
    if(_domainBcLeft == 0)  setNeumannX(_xVelocity, _res, 0, _zRes);
    else setZeroX(_xVelocity, _res, 0, _zRes); 

    if(_domainBcFront == 0)   setNeumannY(_yVelocity, _res, 0, _zRes);
    else setZeroY(_yVelocity, _res, 0, _zRes); 

    if(_domainBcTop == 0) 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++)
            {
                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];
                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];

                _divergence[index] = -_dx * 0.5f * (
                        xright - xleft +
                        yup - ydown +
                        ztop - zbottom );

                // DG: commenting this helps CG to get a better start, 10-20% speed improvement
                // _pressure[index] = 0.0f;
            }
    copyBorderAll(_pressure, 0, _zRes);

    // solve Poisson equation
    solvePressurePre(_pressure, _divergence, _obstacles);

    setObstaclePressure(_pressure, 0, _zRes);

    // project out solution
    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++)
            {
                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;
                }
            }

    if (_pressure) delete[] _pressure;
    if (_divergence) delete[] _divergence;
}




// 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)
{

    // compleately TODO

    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])
            {
                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 (!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.;
                }
                if (!top && bottom) {
                    _pressure[index] += _pressure[index + _slabSize];
                    pcnt += 1.;
                }
                
                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
}

// add buoyancy forces
void FLUID_3D::addBuoyancy(float *heat, float *density, 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++)
            {
                _zForce[index] += *_alpha * density[index] + (*_beta * (heat[index] - _tempAmb)); // DG: was _yForce, changed for Blender
            }
}

// add vorticity to the force field
void FLUID_3D::addVorticity(int zBegin, int zEnd)
{
    //int x,y,z,index;
    if(_vorticityEps<=0.) 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;


    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++)
            {

                int up    = _obstacles[index + _xRes] ? index : index + _xRes;
                int down  = _obstacles[index - _xRes] ? index : index - _xRes;
                float dy  = (up == index || down == index) ? 1.0f / _dx : gridSize;
                int out   = _obstacles[index + _slabSize] ? index : index + _slabSize;
                int in    = _obstacles[index - _slabSize] ? index : index - _slabSize;
                float dz  = (out == index || in == index) ? 1.0f / _dx : gridSize;
                int right = _obstacles[index + 1] ? index : index + 1;
                int left  = _obstacles[index - 1] ? index : index - 1;
                float dx  = (right == index || right == index) ? 1.0f / _dx : gridSize;

                _xVorticity[vIndex] = (_zVelocity[up] - _zVelocity[down]) * dy + (-_yVelocity[out] + _yVelocity[in]) * dz;
                _yVorticity[vIndex] = (_xVelocity[out] - _xVelocity[in]) * dz + (-_zVelocity[right] + _zVelocity[left]) * dx;
                _zVorticity[vIndex] = (_yVelocity[right] - _yVelocity[left]) * dx + (-_xVelocity[up] + _xVelocity[down])* 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] ? vIndex : vIndex + _xRes;
                    int down  = _obstacles[index - _xRes] ? vIndex : vIndex - _xRes;
                    float dy  = (up == vIndex || down == vIndex) ? 1.0f / _dx : gridSize;
                    int out   = _obstacles[index + _slabSize] ? vIndex : vIndex + _slabSize;
                    int in    = _obstacles[index - _slabSize] ? vIndex : vIndex - _slabSize;
                    float dz  = (out == vIndex || in == vIndex) ? 1.0f / _dx : gridSize;
                    int right = _obstacles[index + 1] ? vIndex : vIndex + 1;
                    int left  = _obstacles[index - 1] ? vIndex : vIndex - 1;
                    float dx  = (right == vIndex || right == 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 > 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;
                        _yForce[index] -= (N[0] * _zVorticity[vIndex] - N[2] * _xVorticity[vIndex]) * _dx * eps;
                        _zForce[index] += (N[0] * _yVorticity[vIndex] - N[1] * _xVorticity[vIndex]) * _dx * eps;
                    }
                    }   // 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);

    if(_domainBcLeft == 0) copyBorderX(_xVelocityOld, res, zBegin, zEnd);
    else setZeroX(_xVelocityOld, res, zBegin, zEnd);

    if(_domainBcFront == 0) copyBorderY(_yVelocityOld, res, zBegin, zEnd);
    else setZeroY(_yVelocityOld, res, zBegin, zEnd); 

    if(_domainBcTop == 0) copyBorderZ(_zVelocityOld, res, zBegin, zEnd);
    else 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);
    advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _heatOld, _heatTemp, 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 arrays
    float* t1 = _xForce;

    // advectFieldMacCormack2(dt, xVelocity, yVelocity, zVelocity, oldField, newField, tempfield, temp, res, obstacles)

    advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _densityOld, _density, _densityTemp, t1, res, _obstacles, zBegin, zEnd);
    advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _heatOld, _heat, _heatTemp, 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);

    if(_domainBcLeft == 0) copyBorderX(_xVelocityTemp, res, zBegin, zEnd);
    else setZeroX(_xVelocityTemp, res, zBegin, zEnd);               

    if(_domainBcFront == 0) copyBorderY(_yVelocityTemp, res, zBegin, zEnd);
    else setZeroY(_yVelocityTemp, res, zBegin, zEnd); 

    if(_domainBcTop == 0) copyBorderZ(_zVelocityTemp, res, zBegin, zEnd);
    else setZeroZ(_zVelocityTemp, res, zBegin, zEnd);

    setZeroBorder(_density, res, zBegin, zEnd);
    setZeroBorder(_heat, res, zBegin, zEnd);


    /*int begin=zBegin * _slabSize;
    int end=begin + (zEnd - zBegin) * _slabSize;
  for (int x = begin; x < end; x++)
    _xForce[x] = _yForce[x] = 0.0f;*/
}