Hi Csaba,
Thanks for the fix. This helps me understand the nature and location
of the problem. Reviewing the changes and the original code I'm
thinking it might be simpler to just refactor the code block so the
division by wind_accel.length2() could be delayed, and have the if ()
statement refactored by moving the divisor to the other side of the <
expression :
double compenstated_dt = dt;
if (relative_wind.length2() < dt*dt*wind_accel.length2())
{
double critical_dt2 =
relative_wind.length2()/wind_accel.length2();
osg::notify(osg::NOTICE)<<"** Could be critical:
dt="<<dt<<" critical_dt="<<sqrtf(critical_dt2)<<std::endl;
compenstated_dt = sqrtf(critical_dt2)*0.8f;
}
Does this change make sense to you? If wind_accel.length2() is zero
then the division will never take place, chances of an overflow should
be avoided as well. The changed file is attached, could you try this
out on your models/app. I've done the standard example tests here and
they all work, but then they never highlighted a error so I can't
confirm a fix.
Cheers,
Robert.
On Thu, Mar 12, 2009 at 6:08 PM, Csaba Halász <[email protected]> wrote:
> Hi Robert and all,
>
> Attached modification works around a division by zero. NOTE: it does
> *not* fix potential overflow still lurking there if the wind_accel is
> a small vector (which it probably can be, if I have run into it being
> zero...)
>
> --
> Csaba
>
> _______________________________________________
> osg-submissions mailing list
> [email protected]
> http://lists.openscenegraph.org/listinfo.cgi/osg-submissions-openscenegraph.org
>
>
#include <osgParticle/FluidProgram>
osgParticle::FluidProgram::FluidProgram():
Program()
{
setFluidToAir();
}
osgParticle::FluidProgram::FluidProgram(const FluidProgram& copy, const osg::CopyOp& copyop):
Program(copy, copyop),
_acceleration(copy._acceleration),
_viscosity(copy._viscosity),
_density(copy._density),
_wind(copy._wind),
_viscosityCoefficient(copy._viscosityCoefficient),
_densityCoefficient(copy._densityCoefficient)
{
}
void osgParticle::FluidProgram::execute(double dt)
{
const float four_over_three = 4.0f/3.0f;
ParticleSystem* ps = getParticleSystem();
int n = ps->numParticles();
for (int i=0; i<n; ++i)
{
Particle* particle = ps->getParticle(i);
if (particle->isAlive())
{
float radius = particle->getRadius();
float Area = osg::PI*radius*radius;
float Volume = Area*radius*four_over_three;
// compute force due to gravity + boyancy of displacing the fluid that the particle is emersed in.
osg::Vec3 accel_gravity = _acceleration * ((particle->getMass() - _density*Volume) * particle->getMassInv());
// compute force due to friction
osg::Vec3 velBefore = particle->getVelocity();
osg::Vec3 relative_wind = particle->getVelocity()-_wind;
osg::Vec3 wind_force = - relative_wind * Area * (_viscosityCoefficient + _densityCoefficient*relative_wind.length());
osg::Vec3 wind_accel = wind_force * particle->getMassInv();
double compenstated_dt = dt;
if (relative_wind.length2() < dt*dt*wind_accel.length2())
{
double critical_dt2 = relative_wind.length2()/wind_accel.length2();
osg::notify(osg::NOTICE)<<"** Could be critical: dt="<<dt<<" critical_dt="<<sqrtf(critical_dt2)<<std::endl;
compenstated_dt = sqrtf(critical_dt2)*0.8f;
}
particle->addVelocity(accel_gravity*dt + wind_accel*compenstated_dt);
}
}
}
_______________________________________________
osg-submissions mailing list
[email protected]
http://lists.openscenegraph.org/listinfo.cgi/osg-submissions-openscenegraph.org
