test/intern/smoke/intern/FLUID_3D_STATIC.cpp

/** \file
 * \ingroup smoke
 */
//////////////////////////////////////////////////////////////////////
// This file is part of Wavelet Turbulence.
//
// Wavelet Turbulence is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// Wavelet Turbulence is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Wavelet Turbulence.  If not, see <http://www.gnu.org/licenses/>.
//
// Copyright 2008 Theodore Kim and Nils Thuerey
//
// FLUID_3D.cpp: implementation of the static functions 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 <zlib.h>
#include "FLUID_3D.h"
#include "IMAGE.h"
#include "WTURBULENCE.h"
#include "INTERPOLATE.h"

//////////////////////////////////////////////////////////////////////
// add a test cube of density to the center
//////////////////////////////////////////////////////////////////////
/*
void FLUID_3D::addSmokeColumn() {
	addSmokeTestCase(_density, _res, 1.0);
	// addSmokeTestCase(_zVelocity, _res, 1.0);
	addSmokeTestCase(_heat, _res, 1.0);
	if (_wTurbulence) {
		addSmokeTestCase(_wTurbulence->getDensityBig(), _wTurbulence->getResBig(), 1.0);
	}
}
*/

//////////////////////////////////////////////////////////////////////
// generic static version, so that it can be applied to the
// WTURBULENCE grid as well
//////////////////////////////////////////////////////////////////////

void FLUID_3D::addSmokeTestCase(float* field, Vec3Int res)
{
	const int slabSize = res[0]*res[1]; int maxRes = (int)MAX3V(res);
	float dx = 1.0f / (float)maxRes;

	float xTotal = dx * res[0];
	float yTotal = dx * res[1];

	float heighMin = 0.05;
	float heighMax = 0.10;

	for (int y = 0; y < res[2]; y++)
		for (int z = (int)(heighMin*res[2]); z <= (int)(heighMax * res[2]); z++)
			for (int x = 0; x < res[0]; x++) {
				float xLength = x * dx - xTotal * 0.4f;
				float yLength = y * dx - yTotal * 0.5f;
				float radius = sqrtf(xLength * xLength + yLength * yLength);

				if (radius < 0.075f * xTotal) {
					int index = x + y * res[0] + z * slabSize;
					field[index] = 1.0f;
				}
			}
}


//////////////////////////////////////////////////////////////////////
// set x direction to Neumann boundary conditions
//////////////////////////////////////////////////////////////////////
void FLUID_3D::setNeumannX(float* field, Vec3Int res, int zBegin, int zEnd)
{
	const int slabSize = res[0] * res[1];
	int index;
	for (int z = zBegin; z < zEnd; z++)
		for (int y = 0; y < res[1]; y++)
		{
			// left slab
			index = y * res[0] + z * slabSize;
			field[index] = field[index + 2];

			// right slab
			index = y * res[0] + z * slabSize + res[0] - 1;
			field[index] = field[index - 2];
		}
 }

//////////////////////////////////////////////////////////////////////
// set y direction to Neumann boundary conditions
//////////////////////////////////////////////////////////////////////
void FLUID_3D::setNeumannY(float* field, Vec3Int res, int zBegin, int zEnd)
{
	const int slabSize = res[0] * res[1];
	int index;
	for (int z = zBegin; z < zEnd; z++)
		for (int x = 0; x < res[0]; x++)
		{
			// front slab
			index = x + z * slabSize;
			field[index] = field[index + 2 * res[0]];

			// back slab
			index = x + z * slabSize + slabSize - res[0];
			field[index] = field[index - 2 * res[0]];
		}
}

//////////////////////////////////////////////////////////////////////
// set z direction to Neumann boundary conditions
//////////////////////////////////////////////////////////////////////
void FLUID_3D::setNeumannZ(float* field, Vec3Int res, int zBegin, int zEnd)
{
	const int slabSize = res[0] * res[1];
	const int totalCells = res[0] * res[1] * res[2];
	const int cellsslab = totalCells - slabSize;
	int index;

	if (zBegin == 0) {
		for (int y = 0; y < res[1]; y++)
			for (int x = 0; x < res[0]; x++)
			{
				// front slab
				index = x + y * res[0];
				field[index] = field[index + 2 * slabSize];
			}
	}

	if (zEnd == res[2]) {
		for (int y = 0; y < res[1]; y++)
			for (int x = 0; x < res[0]; x++)
			{
				// back slab
				index = x + y * res[0] + cellsslab;
				field[index] = field[index - 2 * slabSize];
			}
	}
		
}

//////////////////////////////////////////////////////////////////////
// set x direction to zero
//////////////////////////////////////////////////////////////////////
void FLUID_3D::setZeroX(float* field, Vec3Int res, int zBegin, int zEnd)
{
	const int slabSize = res[0] * res[1];
	int index;
	for (int z = zBegin; z < zEnd; z++)
		for (int y = 0; y < res[1]; y++)
		{
			// left slab
			index = y * res[0] + z * slabSize;
			field[index] = 0.0f;

			// right slab
			index += res[0] - 1;
			field[index] = 0.0f;
		}
}

//////////////////////////////////////////////////////////////////////
// set y direction to zero
//////////////////////////////////////////////////////////////////////
void FLUID_3D::setZeroY(float* field, Vec3Int res, int zBegin, int zEnd)
{
	const int slabSize = res[0] * res[1];
	int index;
	for (int z = zBegin; z < zEnd; z++)
		for (int x = 0; x < res[0]; x++)
		{
			// bottom slab
			index = x + z * slabSize;
			field[index] = 0.0f;

			// top slab
			index += slabSize - res[0];
			field[index] = 0.0f;
		}
}

//////////////////////////////////////////////////////////////////////
// set z direction to zero
//////////////////////////////////////////////////////////////////////
void FLUID_3D::setZeroZ(float* field, Vec3Int res, int zBegin, int zEnd)
{
	const int slabSize = res[0] * res[1];
	const int totalCells = res[0] * res[1] * res[2];

	int index = 0;
	if (zBegin == 0)
	{
		for (int y = 0; y < res[1]; y++)
			for (int x = 0; x < res[0]; x++, index++)
			{
				// front slab
				field[index] = 0.0f;
		}
	}

	if (zEnd == res[2])
	{
		index=0;
		int indexx=0;
		const int cellsslab = totalCells - slabSize;

		for (int y = 0; y < res[1]; y++)
			for (int x = 0; x < res[0]; x++, index++)
			{

				// back slab
				indexx = index + cellsslab;
				field[indexx] = 0.0f;
			}
	}
 }
//////////////////////////////////////////////////////////////////////
// copy grid boundary
//////////////////////////////////////////////////////////////////////
void FLUID_3D::copyBorderX(float* field, Vec3Int res, int zBegin, int zEnd)
{
	const int slabSize = res[0] * res[1];
	int index;
	for (int z = zBegin; z < zEnd; z++)
		for (int y = 0; y < res[1]; y++)
		{
			// left slab
			index = y * res[0] + z * slabSize;
			field[index] = field[index + 1];

			// right slab
			index += res[0] - 1;
			field[index] = field[index - 1];
		}
}
void FLUID_3D::copyBorderY(float* field, Vec3Int res, int zBegin, int zEnd)
{
	const int slabSize = res[0] * res[1];
	//const int totalCells = res[0] * res[1] * res[2];
	int index;
	for (int z = zBegin; z < zEnd; z++)
		for (int x = 0; x < res[0]; x++)
		{
			// bottom slab
			index = x + z * slabSize;
			field[index] = field[index + res[0]]; 
			// top slab
			index += slabSize - res[0];
			field[index] = field[index - res[0]];
		}
}
void FLUID_3D::copyBorderZ(float* field, Vec3Int res, int zBegin, int zEnd)
{
	const int slabSize = res[0] * res[1];
	const int totalCells = res[0] * res[1] * res[2];
	int index=0;

	if (zBegin == 0)
	{
		for (int y = 0; y < res[1]; y++)
			for (int x = 0; x < res[0]; x++, index++)
			{
				field[index] = field[index + slabSize]; 
			}
	}

	if (zEnd == res[2])
	{

		index=0;
		int indexx=0;
		const int cellsslab = totalCells - slabSize;

		for (int y = 0; y < res[1]; y++)
			for (int x = 0; x < res[0]; x++, index++)
			{
				// back slab
				indexx = index + cellsslab;
				field[indexx] = field[indexx - slabSize];
			}
	}
}

/////////////////////////////////////////////////////////////////////
// advect field with the semi lagrangian method
//////////////////////////////////////////////////////////////////////
void FLUID_3D::advectFieldSemiLagrange(const float dt, const float* velx, const float* vely,  const float* velz,
		float* oldField, float* newField, Vec3Int res, int zBegin, int zEnd)
{
	const int xres = res[0];
	const int yres = res[1];
	const int zres = res[2];
	const int slabSize = res[0] * res[1];


	for (int z = zBegin; z < zEnd; z++)
		for (int y = 0; y < yres; y++)
			for (int x = 0; x < xres; x++)
			{
				const int index = x + y * xres + z * xres*yres;
				
        // backtrace
				float xTrace = x - dt * velx[index];
				float yTrace = y - dt * vely[index];
				float zTrace = z - dt * velz[index];

				// clamp backtrace to grid boundaries
				if (xTrace < 0.5f) xTrace = 0.5f;
				if (xTrace > xres - 1.5f) xTrace = xres - 1.5f;
				if (yTrace < 0.5f) yTrace = 0.5f;
				if (yTrace > yres - 1.5f) yTrace = yres - 1.5f;
				if (zTrace < 0.5f) zTrace = 0.5f;
				if (zTrace > zres - 1.5f) zTrace = zres - 1.5f;

				// locate neighbors to interpolate
				const int x0 = (int)xTrace;
				const int x1 = x0 + 1;
				const int y0 = (int)yTrace;
				const int y1 = y0 + 1;
				const int z0 = (int)zTrace;
				const int z1 = z0 + 1;

				// get interpolation weights
				const float s1 = xTrace - x0;
				const float s0 = 1.0f - s1;
				const float t1 = yTrace - y0;
				const float t0 = 1.0f - t1;
				const float u1 = zTrace - z0;
				const float u0 = 1.0f - u1;

				const int i000 = x0 + y0 * xres + z0 * slabSize;
				const int i010 = x0 + y1 * xres + z0 * slabSize;
				const int i100 = x1 + y0 * xres + z0 * slabSize;
				const int i110 = x1 + y1 * xres + z0 * slabSize;
				const int i001 = x0 + y0 * xres + z1 * slabSize;
				const int i011 = x0 + y1 * xres + z1 * slabSize;
				const int i101 = x1 + y0 * xres + z1 * slabSize;
				const int i111 = x1 + y1 * xres + z1 * slabSize;

				// interpolate
				// (indices could be computed once)
				newField[index] = u0 * (s0 * (t0 * oldField[i000] +
							t1 * oldField[i010]) +
						s1 * (t0 * oldField[i100] +
							t1 * oldField[i110])) +
					u1 * (s0 * (t0 * oldField[i001] +
								t1 * oldField[i011]) +
							s1 * (t0 * oldField[i101] +
								t1 * oldField[i111]));
			}
}


/////////////////////////////////////////////////////////////////////
// advect field with the maccormack method
//
// comments are the pseudocode from selle's paper
//////////////////////////////////////////////////////////////////////
void FLUID_3D::advectFieldMacCormack1(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity, 
				float* oldField, float* tempResult, Vec3Int res, int zBegin, int zEnd)
{
	/*const int sx= res[0];
	const int sy= res[1];
	const int sz= res[2];

	for (int x = 0; x < sx * sy * sz; x++)
		phiHatN[x] = phiHatN1[x] = oldField[x];*/	// not needed as all the values are written first

	float*& phiN    = oldField;
	float*& phiN1   = tempResult;



	// phiHatN1 = A(phiN)
	advectFieldSemiLagrange(  dt, xVelocity, yVelocity, zVelocity, phiN, phiN1, res, zBegin, zEnd);		// uses wide data from old field and velocities (both are whole)
}



void FLUID_3D::advectFieldMacCormack2(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity, 
				float* oldField, float* newField, float* tempResult, float* temp1, Vec3Int res, const unsigned char* obstacles, int zBegin, int zEnd)
{
	float* phiHatN  = tempResult;
	float* t1  = temp1;
	const int sx= res[0];
	const int sy= res[1];

	float*& phiN    = oldField;
	float*& phiN1   = newField;



	// phiHatN = A^R(phiHatN1)
	advectFieldSemiLagrange( -1.0f*dt, xVelocity, yVelocity, zVelocity, phiHatN, t1, res, zBegin, zEnd);		// uses wide data from old field and velocities (both are whole)

	// phiN1 = phiHatN1 + (phiN - phiHatN) / 2
	const int border = 0; 
	for (int z = zBegin+border; z < zEnd-border; z++)
		for (int y = border; y < sy-border; y++)
			for (int x = border; x < sx-border; x++) {
				int index = x + y * sx + z * sx*sy;
				phiN1[index] = phiHatN[index] + (phiN[index] - t1[index]) * 0.50f;
				//phiN1[index] = phiHatN1[index]; // debug, correction off
			}
	copyBorderX(phiN1, res, zBegin, zEnd);
	copyBorderY(phiN1, res, zBegin, zEnd);
	copyBorderZ(phiN1, res, zBegin, zEnd);

	// clamp any newly created extrema
	clampExtrema(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res, zBegin, zEnd);		// uses wide data from old field and velocities (both are whole)

	// if the error estimate was bad, revert to first order
	clampOutsideRays(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res, obstacles, phiHatN, zBegin, zEnd);	// phiHatN is only used at cells within thread range, so its ok

} 


//////////////////////////////////////////////////////////////////////
// Clamp the extrema generated by the BFECC error correction
//////////////////////////////////////////////////////////////////////
void FLUID_3D::clampExtrema(const float dt, const float* velx, const float* vely,  const float* velz,
		float* oldField, float* newField, Vec3Int res, int zBegin, int zEnd)
{
	const int xres= res[0];
	const int yres= res[1];
	const int zres= res[2];
	const int slabSize = res[0] * res[1];

	int bb=0;
	int bt=0;

	if (zBegin == 0) {bb = 1;}
	if (zEnd == res[2]) {bt = 1;}


	for (int z = zBegin+bb; z < zEnd-bt; z++)
		for (int y = 1; y < yres-1; y++)
			for (int x = 1; x < xres-1; x++)
			{
				const int index = x + y * xres+ z * xres*yres;
				// backtrace
				float xTrace = x - dt * velx[index];
				float yTrace = y - dt * vely[index];
				float zTrace = z - dt * velz[index];

				// clamp backtrace to grid boundaries
				if (xTrace < 0.5f) xTrace = 0.5f;
				if (xTrace > xres - 1.5f) xTrace = xres - 1.5f;
				if (yTrace < 0.5f) yTrace = 0.5f;
				if (yTrace > yres - 1.5f) yTrace = yres - 1.5f;
				if (zTrace < 0.5f) zTrace = 0.5f;
				if (zTrace > zres - 1.5f) zTrace = zres - 1.5f;

				// locate neighbors to interpolate
				const int x0 = (int)xTrace;
				const int x1 = x0 + 1;
				const int y0 = (int)yTrace;
				const int y1 = y0 + 1;
				const int z0 = (int)zTrace;
				const int z1 = z0 + 1;

				const int i000 = x0 + y0 * xres + z0 * slabSize;
				const int i010 = x0 + y1 * xres + z0 * slabSize;
				const int i100 = x1 + y0 * xres + z0 * slabSize;
				const int i110 = x1 + y1 * xres + z0 * slabSize;
				const int i001 = x0 + y0 * xres + z1 * slabSize;
				const int i011 = x0 + y1 * xres + z1 * slabSize;
				const int i101 = x1 + y0 * xres + z1 * slabSize;
				const int i111 = x1 + y1 * xres + z1 * slabSize;

				float minField = oldField[i000];
				float maxField = oldField[i000];

				minField = (oldField[i010] < minField) ? oldField[i010] : minField;
				maxField = (oldField[i010] > maxField) ? oldField[i010] : maxField;

				minField = (oldField[i100] < minField) ? oldField[i100] : minField;
				maxField = (oldField[i100] > maxField) ? oldField[i100] : maxField;

				minField = (oldField[i110] < minField) ? oldField[i110] : minField;
				maxField = (oldField[i110] > maxField) ? oldField[i110] : maxField;

				minField = (oldField[i001] < minField) ? oldField[i001] : minField;
				maxField = (oldField[i001] > maxField) ? oldField[i001] : maxField;

				minField = (oldField[i011] < minField) ? oldField[i011] : minField;
				maxField = (oldField[i011] > maxField) ? oldField[i011] : maxField;

				minField = (oldField[i101] < minField) ? oldField[i101] : minField;
				maxField = (oldField[i101] > maxField) ? oldField[i101] : maxField;

				minField = (oldField[i111] < minField) ? oldField[i111] : minField;
				maxField = (oldField[i111] > maxField) ? oldField[i111] : maxField;

				newField[index] = (newField[index] > maxField) ? maxField : newField[index];
				newField[index] = (newField[index] < minField) ? minField : newField[index];
			}
}

//////////////////////////////////////////////////////////////////////
// Reverts any backtraces that go into boundaries back to first 
// order -- in this case the error correction term was totally
// incorrect
//////////////////////////////////////////////////////////////////////
void FLUID_3D::clampOutsideRays(const float dt, const float* velx, const float* vely,  const float* velz,
				float* oldField, float* newField, Vec3Int res, const unsigned char* obstacles, const float *oldAdvection, int zBegin, int zEnd)
{
	const int sx= res[0];
	const int sy= res[1];
	const int sz= res[2];
	const int slabSize = res[0] * res[1];

	int bb=0;
	int bt=0;

	if (zBegin == 0) {bb = 1;}
	if (zEnd == res[2]) {bt = 1;}

	for (int z = zBegin+bb; z < zEnd-bt; z++)
		for (int y = 1; y < sy-1; y++)
			for (int x = 1; x < sx-1; x++)
			{
				const int index = x + y * sx+ z * slabSize;
				// backtrace
				float xBackward = x + dt * velx[index];
				float yBackward = y + dt * vely[index];
				float zBackward = z + dt * velz[index];
				float xTrace    = x - dt * velx[index];
				float yTrace    = y - dt * vely[index];
				float zTrace    = z - dt * velz[index];

				// see if it goes outside the boundaries
				bool hasObstacle = 
					(zTrace < 1.0f)    || (zTrace > sz - 2.0f) ||
					(yTrace < 1.0f)    || (yTrace > sy - 2.0f) ||
					(xTrace < 1.0f)    || (xTrace > sx - 2.0f) ||
					(zBackward < 1.0f) || (zBackward > sz - 2.0f) ||
					(yBackward < 1.0f) || (yBackward > sy - 2.0f) ||
					(xBackward < 1.0f) || (xBackward > sx - 2.0f);
				// reuse old advection instead of doing another one...
				if(hasObstacle) { newField[index] = oldAdvection[index]; continue; }

				// clamp to prevent an out of bounds access when looking into
				// the _obstacles array
				zTrace = (zTrace < 0.5f) ? 0.5f : zTrace;
				zTrace = (zTrace > sz - 1.5f) ? sz - 1.5f : zTrace;
				yTrace = (yTrace < 0.5f) ? 0.5f : yTrace;
				yTrace = (yTrace > sy - 1.5f) ? sy - 1.5f : yTrace;
				xTrace = (xTrace < 0.5f) ? 0.5f : xTrace;
				xTrace = (xTrace > sx - 1.5f) ? sx - 1.5f : xTrace;

				// locate neighbors to interpolate,
				// do backward first since we will use the forward indices if a
				// reversion is actually necessary
				zBackward = (zBackward < 0.5f) ? 0.5f : zBackward;
				zBackward = (zBackward > sz - 1.5f) ? sz - 1.5f : zBackward;
				yBackward = (yBackward < 0.5f) ? 0.5f : yBackward;
				yBackward = (yBackward > sy - 1.5f) ? sy - 1.5f : yBackward;
				xBackward = (xBackward < 0.5f) ? 0.5f : xBackward;
				xBackward = (xBackward > sx - 1.5f) ? sx - 1.5f : xBackward;

				int x0 = (int)xBackward;
				int x1 = x0 + 1;
				int y0 = (int)yBackward;
				int y1 = y0 + 1;
				int z0 = (int)zBackward;
				int z1 = z0 + 1;
				if(obstacles && !hasObstacle) {
					hasObstacle = hasObstacle || 
						obstacles[x0 + y0 * sx + z0*slabSize] ||
						obstacles[x0 + y1 * sx + z0*slabSize] ||
						obstacles[x1 + y0 * sx + z0*slabSize] ||
						obstacles[x1 + y1 * sx + z0*slabSize] ||
						obstacles[x0 + y0 * sx + z1*slabSize] ||
						obstacles[x0 + y1 * sx + z1*slabSize] ||
						obstacles[x1 + y0 * sx + z1*slabSize] ||
						obstacles[x1 + y1 * sx + z1*slabSize] ;
				}
				// reuse old advection instead of doing another one...
				if(hasObstacle) { newField[index] = oldAdvection[index]; continue; }

				x0 = (int)xTrace;
				x1 = x0 + 1;
				y0 = (int)yTrace;
				y1 = y0 + 1;
				z0 = (int)zTrace;
				z1 = z0 + 1;
				if(obstacles && !hasObstacle) {
					hasObstacle = hasObstacle || 
						obstacles[x0 + y0 * sx + z0*slabSize] ||
						obstacles[x0 + y1 * sx + z0*slabSize] ||
						obstacles[x1 + y0 * sx + z0*slabSize] ||
						obstacles[x1 + y1 * sx + z0*slabSize] ||
						obstacles[x0 + y0 * sx + z1*slabSize] ||
						obstacles[x0 + y1 * sx + z1*slabSize] ||
						obstacles[x1 + y0 * sx + z1*slabSize] ||
						obstacles[x1 + y1 * sx + z1*slabSize] ;
				} // obstacle array
				// reuse old advection instead of doing another one...
				if(hasObstacle) { newField[index] = oldAdvection[index]; continue; }

				// see if either the forward or backward ray went into
				// a boundary
				if (hasObstacle) {
					// get interpolation weights
					float s1 = xTrace - x0;
					float s0 = 1.0f - s1;
					float t1 = yTrace - y0;
					float t0 = 1.0f - t1;
					float u1 = zTrace - z0;
					float u0 = 1.0f - u1;

					const int i000 = x0 + y0 * sx + z0 * slabSize;
					const int i010 = x0 + y1 * sx + z0 * slabSize;
					const int i100 = x1 + y0 * sx + z0 * slabSize;
					const int i110 = x1 + y1 * sx + z0 * slabSize;
					const int i001 = x0 + y0 * sx + z1 * slabSize;
					const int i011 = x0 + y1 * sx + z1 * slabSize;
					const int i101 = x1 + y0 * sx + z1 * slabSize;
					const int i111 = x1 + y1 * sx + z1 * slabSize;

					// interpolate, (indices could be computed once)
					newField[index] = u0 * (s0 * (
								t0 * oldField[i000] +
								t1 * oldField[i010]) +
							s1 * (t0 * oldField[i100] +
								t1 * oldField[i110])) +
						u1 * (s0 * (t0 * oldField[i001] +
									t1 * oldField[i011]) +
								s1 * (t0 * oldField[i101] +
									t1 * oldField[i111])); 
				}
			} // xyz
}