Hey all,

I'm seeing some strange behavior from some SPH fluid code I'm working on. I modified it to have a cylinder as a container/collision shape and that works fine. However, I then change it so that it slopes inward as it gets taller, and the particles start to move around constantly like the fluid is boiling and they never settle down. Also, the fluid level doesn't get higher as the container slopes inward as I'd expect. See image:

[attachment=6419:fluid.gif]

Here is the code I'm using to calculate positions (it's a modified form of the code at http://www.rchoetzle...hics/fluids.htm)



accel = p->sph_force;
accel *= m_Param[SPH_PMASS];
.....
// p is the particle structure, with z being up
// radius is the particle radius
// norm is a 3-vector

float boundary = abs((-.4 * p->pos.z) + 9.8);

diff = 2.0f * radius - (boundary - (float)fastSqrt_2(p->pos.x*p->pos.x + p->pos.y*p->pos.y)) * ss;

if (((diff < 8 *radius) && (diff > EPSILON)))
{
norm.Set ( -p->pos.x, -p->pos.y, 0 );
norm.Normalize();

adj = stiff * diff - damp * norm.Dot ( p->vel_eval );
accel.x += adj * norm.x; accel.y += adj * norm.y; accel.z += adj * norm.z;
}

....

// Leapfrog Integration ----------------------------
vnext = accel;
vnext *= m_DT;
vnext += p->vel; // v(t+1/2) = v(t-1/2) + a(t) dt
p->vel_eval = p->vel;
p->vel_eval += vnext;
p->vel_eval *= 0.5; // v(t+1) = [v(t-1/2) + v(t+1/2)] * 0.5 used to compute forces later
p->vel = vnext;
vnext *= m_DT/ss;
p->pos += vnext; // p(t+1) = p(t) + v(t+1/2) dt



I'm pretty sure the way the normal is being added to the acceleration is causing the "boiling" of the particles, but I'm not sure. Thanks for any help

Here I am assuming the norm variable holds the normal vector of the container, which would indicate to me that it needs to have a non-zero z value if the container is a frustrum instead of a cylinder.

Here I am assuming the norm variable holds the normal vector of the container, which would indicate to me that it needs to have a non-zero z value if the container is a frustrum instead of a cylinder.

Thanks for the reply, you are correct the norm variable holds the normal vector of the container. Unfortunately when I change it to have a non-zero z value I still get the erratic particle movement, i.e :


norm.Set ( -p->pos.x, -p->pos.y, p->pos.z );
norm.Normalize();


Maybe I need something different for the z value?

do you understand what the normal to a surface means?

you've been calculating the normal from the particle position, which is a little strange but I guess it works if you know the centre of the cylinder is at (0,0).

[font="Arial"]The equation for a cone with it's apex at the origin and axis parallel to z with aperture 2? is:
f(u) = u.z - |u|.cos ? = 0

Since your cone won't have it's apex at the origin, you can simply offset u.z by -h for apex at height 'h'.

Letting d = sqrt(x[sup]2[/sup]+y[sup]2[/sup]+(z-h)[sup]2[/sup])
the partial derivatives of f(u+h.z) are:
pf/px = -(x/d).cos ?, pf/py = -(y/d).cos ?, pf/pz = 1 - ((z-h)/d).cos ?

which are the components of your normal (to be normalised)[/font]


you can speed up the normalisation as we can determine the magnitude of the vector (pf/px, pf/py, pf/pz) to be sin [font=Arial][size=2]? as long as you know the point is on the surface of the cone[/font]

do you understand what the normal to a surface means?

I do, but this is my first attempt at coding collision detection, probably why I got confused on the implementation

Thanks for the help, it was the z component that was causing problems, I finally ended up using this based on your code:

norm.Set ( -p->pos.x/diff, -p->pos.y/diff, 1.0-(p->pos.z/diff) );

along with the initial boundary slope code:
float boundary = abs((-.4 * p->pos.z) + 9.8);

Thanks again!