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