[loki74]

Okay, so this is going to sound very VERY crazy... but just bear with me. I have several questions, if you can answer ANY that would be great! Depending on how much you know about this type of stuff, you may or may not need to read all my explaining.

I love 3D animation... but some aspects of it can be quite expensive... namely physics-based 3D dynamic fluid simulation. As such, I have set out to write my own, HIGHLY SIMPLIFIED fluid solver which will convincingly (note: not accurately) depict fluid phenomena such as gasses, liquids, etc.

I only have a high school Algebra 2 mathematics education. I KNOW ITS CRAZY. Please do NOT TELL ME TO GIVE UP! I've come pretty far already thanks to one Jos Stam. And if I was willing to wait I would just save up for RealFlow, a very expensive, industry standard SPH based fluid solver.

Right now I am working on the gas part of the solver.

A fellow named Jos Stam has provided a paper in which he implements a superfast fluid solver in less than 100 lines of readable C code.

This solver, like many of the solvers in CFD and CGI, is based on a voxel-grid. That is, the state of the fluid is known at only certain points across the domain, in this case, the centers of each grid cell. The solver involves two grids--a vector field on which velocity data is stored, and a scalar field in which density data is stored.

Right now I am working on the density part of the gas solver.

The governing equation for the density field is this:

∂ρ/∂t = - (u • ∇)ρ + κ∇²ρ + S

the variables are as follows:
∂t is some arbitrary timestep by which we advance the simulation each frame.
u is the corresponding cell in the aforementioned vector field.
ρ is the density of the cell in question
κ is the rate of diffusion
S represents a source of fluid

The three terms on the left hand side are as follows:
- (u • ∇)ρ means that the density advects according to a vector u
κ∇²ρ means that the density diffuses over time
S means that the density increases due to sources of fluid (ie the top of a smokestack)

Right now, this is all I care about--κ∇²ρ The DIFFUSION STEP!

NOW, the simplest approach would be that the current cell exchanges equally with its neighbouring cells. This is easy to code, except, if the cell diffuses beyond its neighbours, this is not accounted for, and the simulation "blows up." This happens when the grid is too coarse (cells too big), the diffusion rate kappa is too big, or the time step is too big.

Mr. Stam implemented an interative solver that uses Gauss-Seidel relaxation to solve a linear system for the diffusion. This is nice, because it is fast, easy to code, and will never blow up, regarless of grid cell size, timestep, or rate of diffusion. It's the perfect solution... but here's the catch--he patented the method!

Now I don't want to wait around for my solver to perform many many substeps to achieve proper diffusion... I realized the effects of diffusion were similar to that of the Gaussian Blur filter in Photoshop.

So I look up how to do this on Wikipedia, and it provides me with the following equation:

G(u,v) = (1/2πσ²)e^-(u²+v²)/(2σ²)

This will provide you with a convulsion kernel with Gaussian distribution, thus allowing you to convolve the image, resulting in a Gaussian Blur.

However, it is for only 2 dimensions--u and v. I need it in 3.

HERE IS MY FIRST QUESTION:
Would the 3D form of the above function be,
G(u,v,w) = (1/2πσ²)e^-(u²+v²+w²)/(2σ²)????
or would I need to do something more to make it work in 3D?

Moving on... the example Jos provides in his paper is 2D, and here is a sample of his diffusion code:
[FONT="Courier New"]
void diffuse_bad ( int N, int b, float * x, float * x0, float diff, float dt ) {

int i, j;
float a=dt*diff*N*N;

for ( i=1 ; i<=N ; i++ ) {

for ( j=1 ; j<=N ; j++ ) {

x[IX(i,j)] = x0[IX(i,j)] + a*(x0[IX(i-1,j)]+x0[IX(i+1,j)]+ x0[IX(i,j-1)]+x0[IX(i,j+1)]-4*x0[IX(i,j)]);

}

}
set_bnd ( N, b, x );

} [/FONT]

Now, he uses IX() to retrieve his variables two-dimensionally from a 1D array, so that's what that's about... what I'm really curious about is this line:

[FONT="Courier New"]float a=dt*diff*N*N;[/FONT]

a is what everything gets multiplied by. and this makes sense... dt is the timestep, diff is the rate of diffusion... but I notice the N*N. In his code, N is the size of one dimension of the domain. So the domain is N². The only thing I can see in the equation that this represents is the ∇² thing.

HERE IS MY SECOND QUESTION:
does ∇² = N², and if so, would it in 3D?

HERE IS MY THIRD QUESTION:
What would ∇² be in the case of a nonqubic 3D domain? (Say I wanted the domain to be 200x200x75)?

ANY help on this matter would be MUCH appreciated!

Thanks in advance,
-loki

PS I know that if I only use this for personal use I can use the patented method--but I may want to distribute this to a few select friends, and I'm not sure if that would be legal, given that case. And in ANY case, I wouldn't want to risk it AT ALL seeing as how I don't have an attorney... Thank you all again!

[olePigeon]

42?

[The Godfather]

Can't you use Core Image to diffuse your space slices, once in each direction?

http://www.vis.uni-stuttgart.de/ger/.../vis99hopf.pdf
http://www.apple.com/macosx/features/coreimage/

It is totally cool that you are doing this on your own initiative, for the advancement of the art.

I don't claim to know the answer, but G(u,v,w) = (1/2πσ²)e^-(u²+v²+w²)/(2σ²) needs a different scaling constant. (1/2πσ²) works for 2D because if you integrate the kernel it over the UV plane, it yields 1.

I assume that any Gaussian-like 3D convolution kernel will diffuse your voxels. Why don't you try your formula even if it is most likely wrong? Then tweak with the scaling constant.

Do you have pics?

[rickey939]

Paging Ghoser777 to the Lounge....Ghoser777, please come to the Lounge immediately.

[tie]

Using a Gaussian blur will be more expensive computationally (a big deal in 3D), and will give inferior results. For example, it can't naturally handle fluid boundaries.

To answer your questions:
If you did you use Gaussian blur, the normalization would be 1/(2 pi sigma^2)^{3/2}. (The integral factors: e^{x^2+y^2+z^2} = e^{x^2} e^{y^2} e^{z^2}.)

The multiplication by N*N is because he is approximating the second differential (d^2X/dx^2+d^2X/dy^2) -- that is, he is dividing by the square of the cell width, dx=1/N. This won't change in 3D -- divide by the square of the cell dimension.

In other words,

Code:
x0[IX(i,j)] + a*(x0[IX(i-1,j)]+x0[IX(i+1,j)]+ x0[IX(i,j-1)]+x0[IX(i,j+1)]-4*x0[IX(i,j)])

should be, in 3D,

Code:
 x0[IX(i,j,k)] + a*(x0[IX(i-1,j,k)]+x0[IX(i+1,j,k)]+ x0[IX(i,j-1,k)]+x0[IX(i,j+1,k)]+x0[IX(i,j,k-1)]+x0[IX(i,j,k+1)]-6*x0[IX(i,j,k)])

Use Stam's code and save yourself the trouble. The patent isn't an issue for you.

[Ulrich Kinbote]

*Innumerate enters thread, and slowly backs out again*

[The Godfather]

Originally Posted by rickey939

Paging Ghoser777 to the Lounge....Ghoser777, please come to the Lounge immediately.

You should PM him.

[rickey939]

Originally Posted by The Godfather

You should PM him.

Nah, PM is o'skool. I caught him on iChat.

[tavilach]

Whoa. You lost me on fluid .

[loki74]

Originally Posted by tie

Using a Gaussian blur will be more expensive computationally (a big deal in 3D), and will give inferior results. For example, it can't naturally handle fluid boundaries.

To answer your questions:
If you did you use Gaussian blur, the normalization would be 1/(2 pi sigma^2)^{3/2}. (The integral factors: e^{x^2+y^2+z^2} = e^{x^2} e^{y^2} e^{z^2}.)

The multiplication by N*N is because he is approximating the second differential (d^2X/dx^2+d^2X/dy^2) -- that is, he is dividing by the square of the cell width, dx=1/N. This won't change in 3D -- divide by the square of the cell dimension.

In other words,

Code:
x0[IX(i,j)] + a*(x0[IX(i-1,j)]+x0[IX(i+1,j)]+ x0[IX(i,j-1)]+x0[IX(i,j+1)]-4*x0[IX(i,j)])

should be, in 3D,

Code:
 x0[IX(i,j,k)] + a*(x0[IX(i-1,j,k)]+x0[IX(i+1,j,k)]+ x0[IX(i,j-1,k)]+x0[IX(i,j+1,k)]+x0[IX(i,j,k-1)]+x0[IX(i,j,k+1)]-6*x0[IX(i,j,k)])

Use Stam's code and save yourself the trouble. The patent isn't an issue for you.

Ah... so basically for any number of dimensions n, we take the 2*pi to the n/2 power, and the sigma to the n power. (And of course change the power e is raised to accordingly). That makes sense.

Anyway, I have taken into consideration the boundary handling issue... it is highly experimental, I know. I did post the question on Wikipedia's help desk, and the people there say its pretty sound.

so in 3d, a would still equal dt*diff*N*N?

thanks!

[residentEvil]

where is the button you click that just does it?

my head hurts.