On Fri, Mar 13, 2009 at 11:16 AM, Robert Osfield
<[email protected]> wrote:
> 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.
Hi Robert,
yeah, looks good. I have tried to come up with something similar, but
I didn't like my solution.
I wonder if the double calls to length2() get optimized by the
compiler, maybe it would be more efficient to store them in a
temporary variable. Also depends on how often the condition is true,
but shouldn't hurt.
> 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.
This should definitely avoid the division by zero, as you say, but
since the problem only occurs randomly for us I can not directly
verify right away.
--
Thanks,
Csaba
#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;
double relative_wind_length2 = relative_wind.length2();
double wind_accel_length2 = wind_accel.length2();
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-users mailing list
[email protected]
http://lists.openscenegraph.org/listinfo.cgi/osg-users-openscenegraph.org
