[ODE] Hydrodynamic Forces Demo

[Ander Taylor]

I have posted this before.

Here you go!

void BoxGeom::ApplyHydrodynamicForces(dReal viscosity)
    {
        const dReal *lvel = dBodyGetLinearVel(this->_body->Id);
        const dReal *avel = dBodyGetAngularVel(this->_body->Id);
        const dReal *R = dBodyGetRotation(this->_body->Id);
        dReal ss[3];
        dGeomBoxGetLengths(this->_id, ss);

        dReal AreaX = ss[1] * ss[2];
        dReal AreaY = ss[0] * ss[2];
        dReal AreaZ = ss[0] * ss[1];

        dReal nx = (R[0] * lvel[0] + R[4] * lvel[1] + R[8] * lvel[2]) * 
AreaX;
        dReal ny = (R[1] * lvel[0] + R[5] * lvel[1] + R[9] * lvel[2]) * 
AreaY;
        dReal nz = (R[2] * lvel[0] + R[6] * lvel[1] + R[10] * lvel[2]) * 
AreaZ;

        dReal temp = -nx * viscosity;
        dBodyAddForce(this->_body->Id, temp * R[0], temp * R[4], temp * 
R[8]);

        temp = -ny * viscosity;
        dBodyAddForce(this->_body->Id, temp * R[1], temp * R[5], temp * 
R[9]);

        temp =-nz * viscosity;
        dBodyAddForce(this->_body->Id, temp * R[2], temp * R[6], temp * 
R[10]);

        nx = (R[0] * avel[0] + R[4] * avel[1] + R[8] * avel[2]) * AreaZ;
        ny = (R[1] * avel[0] + R[5] * avel[1] + R[9] * avel[2]) * AreaX;
        nz = (R[2] * avel[0] + R[6] * avel[1] + R[10] * avel[2]) * AreaY;

        temp = -nx * viscosity * 500;
        dBodyAddTorque(this->_body->Id, temp * R[0], temp * R[4], temp * 
R[8]);

        temp = -ny * viscosity * 500;
        dBodyAddTorque(this->_body->Id, temp * R[1], temp * R[5], temp * 
R[9]);

        temp = -nz * viscosity * 500;
        dBodyAddTorque(this->_body->Id, temp * R[2], temp * R[6], temp * 
R[10]);

        ////this might actually work now with the fabs()'s (may need to 
adjust scale):
        //dReal temp = -viscosity*fabs(nx+ny+nz);
        //dBodyAddForce(this->_body->Id, temp*lvel[0], temp*lvel[1], 
temp*lvel[2]);

    }

    void BoxGeom::ApplyHydrodynamicForces(dReal linear_viscosity, dReal 
angular_viscosity, dReal density, Vector3 flow)
    {
        const dReal *lvel = dBodyGetLinearVel(this->_body->Id);
        const dReal *avel = dBodyGetAngularVel(this->_body->Id);
        const dReal *R = dBodyGetRotation(this->_body->Id);

        dVector3 compFlow;
        compFlow[0] = lvel[0] - flow.x;
        compFlow[1] = lvel[1] - flow.y;
        compFlow[2] = lvel[2] - flow.z;

        dReal ss[3];
        dGeomBoxGetLengths(this->_id, ss);

        dReal AreaX = ss[1] * ss[2];
        dReal AreaY = ss[0] * ss[2];
        dReal AreaZ = ss[0] * ss[1];

        dReal nx = (R[0] * compFlow[0] + R[4] * compFlow[1] + R[8] * 
compFlow[2]) * AreaX;
        dReal ny = (R[1] * compFlow[0] + R[5] * compFlow[1] + R[9] * 
compFlow[2]) * AreaY;
        dReal nz = (R[2] * compFlow[0] + R[6] * compFlow[1] + R[10] * 
compFlow[2]) * AreaZ;

        dReal temp = -nx * linear_viscosity;
        dBodyAddForce(this->_body->Id, temp * R[0], temp * R[4], temp * 
R[8]);

        temp = -ny * linear_viscosity;
        dBodyAddForce(this->_body->Id, temp * R[1], temp * R[5], temp * 
R[9]);

        temp = -nz * linear_viscosity;
        dBodyAddForce(this->_body->Id, temp * R[2], temp * R[6], temp * 
R[10]);

        nx = (R[0] * avel[0] + R[4] * avel[1] + R[8] * avel[2]) * (AreaY + 
AreaZ);
        ny = (R[1] * avel[0] + R[5] * avel[1] + R[9] * avel[2]) * (AreaX + 
AreaZ);
        nz = (R[2] * avel[0] + R[6] * avel[1] + R[10] * avel[2]) * (AreaX + 
AreaY);

        temp = -nx * angular_viscosity;
        dBodyAddTorque(this->_body->Id, temp * R[0], temp * R[4], temp * 
R[8]);

        temp = -ny * angular_viscosity;
        dBodyAddTorque(this->_body->Id, temp * R[1], temp * R[5], temp * 
R[9]);

        temp = -nz * angular_viscosity;
        dBodyAddTorque(this->_body->Id, temp * R[2], temp * R[6], temp * 
R[10]);

        dReal gravity[3];
        dWorldGetGravity(this->_body->get_WorldId(), gravity);

        temp = -density * ss[0] * ss[1] * ss[2];
        dBodyAddForce(this->_body->Id, temp * gravity[0], temp * gravity[1], 
temp * gravity[2]);

        //which, unless you have a non-axis-aligned gravity, is probably 
just:
        //dBodyAddForce(body, 0, temp*gravity[1], 0);

    }

As I have said before, this is from a Managed C++ Wrapper so you will have 
to interpret it.

Cheers,

Ander

----Original Message Follows----
Hi Ander,

Sorry to ask, but it is also possible to have the source code? I would like 
to use it in my application, if it would be OK with you, obviously.

Best regards,

Artur Matos.

 >
 > -----Original Message-----
 > From: Ander Taylor [mailto:ander_taylor at hotmail.com]
 > Sent: miércoles, 12 de noviembre de 2003 13:28
 > To: amatos at create.human.nagoya-u.ac.jp; ode at q12.org
 > Subject: Re: [ODE] Hydrodynamic Forces Demo
 >
 > Hi Artur,
 >
 > You can get it from here.
 >
 > http://home.swiftdsl.com.au/~ander
 >
 > Cheers,
 >
 > Ander
 > _______________________________________________
 > ODE mailing list
 > ODE at q12.org
 > http://q12.org/mailman/listinfo/ode
 >
 > _________________________________________________________________
 > Hot chart ringtones and polyphonics. Go to
 > http://ninemsn.com.au/mobilemania/default.asp
 >
 > _______________________________________________
 > ODE mailing list
 > ODE at q12.org
 > http://q12.org/mailman/listinfo/ode
 >
 > _________________________________________________________________
 > Hot chart ringtones and polyphonics. Go to
 > http://ninemsn.com.au/mobilemania/default.asp
 >
 > _______________________________________________
 > ODE mailing list
 > ODE at q12.org
 > http://q12.org/mailman/listinfo/ode
 >



_______________________________________________
ODE mailing list
ODE at q12.org
http://q12.org/mailman/listinfo/ode

_________________________________________________________________
Hot chart ringtones and polyphonics. Go to  
http://ninemsn.com.au/mobilemania/default.asp