烟雾模拟源工程代码

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