/*?Author:?Johannes?Schmid?2006?johnny@grob.org?*/
#include?“fluid.h“
#include?
#include?

Fluid::Fluid（）
{
	int?i;

	for?（i=0;?i<10;?i++）
		clear_buffer（buffers[i]）;

	i=0;
	d=buffers[i++];?d0=buffers[i++];
	u=buffers[i++];?u0=buffers[i++];
	v=buffers[i++];?v0=buffers[i++];
	w=buffers[i++];?w0=buffers[i++];

	clear_sources（）;
}

Fluid::~Fluid（）
{
}

#define?NO_BOUNDARY
void?Fluid::set_bnd（int?b?float*?x）
{
#ifndef?NO_BOUNDARY
	int?i?j;
	for?（i=1;?i<=N;?i++）
	{
		for?（j=1;?j<=N;?j++）?{
			x[_I（0ij）]????=?（b==1）???-x[_I（1ij）]?:?x[_I（1ij）];
			x[_I（N+1ij）]??=?（b==1）???-x[_I（Nij）]?:?x[_I（Nij）];
			x[_I（i0j）]????=?（b==2）???-x[_I（i1j）]?:?x[_I（i1j）];
			x[_I（iN+1j）]??=?（b==2）???-x[_I（iNj）]?:?x[_I（iNj）];
			x[_I（ij0）]????=?（b==3）???-x[_I（ij1）]?:?x[_I（ij1）];
			x[_I（ijN+1）]??=?（b==3）???-x[_I（ijN）]?:?x[_I（ijN）];
		}
	}

	x[_I（000）]???????=?（x[_I（100）]????+x[_I（010）]????+x[_I（001）]）????/3;
	x[_I（0N+10）]?????=?（x[_I（1N+10）]??+x[_I（0N0）]????+x[_I（0N+11）]）??/3;
	x[_I（N+100）]?????=?（x[_I（N00）]????+x[_I（N+110）]??+x[_I（N+101）]）??/3;
	x[_I（N+1N+10）]???=?（x[_I（NN+10）]??+x[_I（N+1N0）]??+x[_I（N+1N+11）]）/3;
	x[_I（00N+1）]?????=?（x[_I（10N+1）]??+x[_I（01N+1）]??+x[_I（00N）]）????/3;
	x[_I（0N+1N+1）]???=?（x[_I（1N+1N+1）]+x[_I（0NN+1）]??+x[_I（0N+1N）]）??/3;
	x[_I（N+10N+1）]???=?（x[_I（N0N+1）]??+x[_I（N+11N+1）]+x[_I（N+10N）]）??/3;
	x[_I（N+1N+1N+1）]?=?（x[_I（NN+1N+1）]+x[_I（N+1NN+1）]+x[_I（N+1N+1N）]）/3;
#endif
}

void?Fluid::add_source（float*?src?float?*dst?float?dt）
{
	int?i?size=（N+2）*（N+2）*（N+2）;

	for?（i=0;?i		dst[i]?+=?src[i]*dt;
}

void?Fluid::add_buoyancy（float?dt）
{
	int?i?size=（N+2）*（N+2）*（N+2）;

	for?（i=0;?i		v[i]?+=?-d[i]*buoyancy*dt;
}

inline?void?Fluid::diffuse（int?b?float*?x0?float*?x?float?diff?float?dt）
{
	int?i?j?k?l;
	float?a=dt*diff*N*N*N;
	for?（l=0;?l<20;?l++）?
	{
		for?（k=1;?k<=N;?k++）
		{
			for?（j=1;?j<=N;?j++）
			{
				for?（i=1;?i<=N;?i++）
				{
					x[_I（ijk）]?=?（x0[_I（ijk）]?+?a*（
						x[_I（i-1jk）]+x[_I（i+1jk）]+
						x[_I（ij-1k）]+x[_I（ij+1k）]+
						x[_I（ijk-1）]+x[_I（ijk+1）]））/（1+6*a）;
				}
			}
		}
		set_bnd（bx）;
	}
}

inline?void?Fluid::advect（int?b?float*?x0?float*?x?float*?uu?float*?vv?float*?ww?float?dt）
{
	int?i?j?k?i0?j0?k0?i1?j1?k1;
	float?sx0?sx1?sy0?sy1?sz0?sz1?v0?v1;
	float?xx?yy?zz?dt0;
	dt0?=?dt*N;
	for?（k=1;?k<=N;?k++）
	{
		for?（j=1;?j<=N;?j++）
		{
			for?（i=1;?i<=N;?i++）
			{
				xx?=?i-dt0*uu[_I（ijk）];
				yy?=?j-dt0*vv[_I（ijk）];
				zz?=?k-dt0*ww[_I（ijk）];
				if?（xx<0.5）?xx=0.5f;?if?（xx>N+0.5）?xx=N+0.5f;?i0=（int）xx;?i1=i0+1;
				if?（yy<0.5）?yy=0.5f;?if?（yy>N+0.5）?yy=N+0.5f;?j0=（int）yy;?j1=j0+1;
				if?（zz<0.5）?zz=0.5f;?if?（zz>N+0.5）?zz=N+0.5f;?k0=（int）zz;?k1=k0+1;
				sx1?=?xx-i0;?sx0?=?1-sx1;
				sy1?=?yy-j0;?sy0?=?1-sy1;
				sz1?=?zz-k0;?sz0?=?1-sz1;
				v0?=?sx0*（sy0*x0[_I（i0j0k0）]+sy1*x0[_I（i0j1k0）]）+sx1*（sy0*x0[_I（i1j0k0）]+sy1*x0[_I