//--------------------------------------------------------------------------------------------------
// Melax, Volint.cpp, http://melax.com
void lmPoly::ComputeMass( lmMassData* md ) const
{
    lmVec3 diag;
    lmVec3 offDiag;
    lmVec3 centroid;
    lmIdentity( diag );
    lmIdentity( offDiag );
    lmIdentity( centroid );
    r32 volume = r32( 0.0 );

    for ( i32 i = 0; i < hull->faceCount; ++i )
    {
        const lmHalfEdge *start = hull->GetEdge( hull->GetFace( i )->edge );
        const lmHalfEdge *middle = hull->GetEdge( start->next );
        const lmHalfEdge *last = hull->GetEdge( middle->next );

        while ( last != start )
        {
            lmVec3 u = hull->GetVertex( start->vert );
            lmVec3 v = hull->GetVertex( middle->vert );
            lmVec3 w = hull->GetVertex( last->vert );

            u = lmMul( local, u );
            v = lmMul( local, v );
            w = lmMul( local, w );

            r32 currentVolume = lmScalarTripleProduct( u, v, w );
            volume += currentVolume;
            centroid += (u + v + w) * currentVolume;

            for ( i32 j = 0; j < 3; ++j )
            {
                i32 j1 = (j + 1) % 3;
                i32 j2 = (j + 2) % 3;

                diag[ j ] += (
                    u[ j ] * v[ j ] + v[ j ] * w[ j ] + w[ j ] * u[ j ] +
                    u[ j ] * u[ j ] + v[ j ] * v[ j ] + w[ j ] * w[ j ]) * currentVolume;

                offDiag[ j ] += (
                    u[ j1 ] * v[ j2 ]        + v[ j1 ] * w[ j2 ]        + w[ j1 ] * u[ j2 ]  +
                    u[ j1 ] * w[ j2 ]        + v[ j1 ] * u[ j2 ]        + w[ j1 ] * v[ j2 ]  +
                    u[ j1 ] * u[ j2 ] * 2.0f + v[ j1 ] * v[ j2 ] * 2.0f + w[ j1 ] * w[ j2 ] * 2.0f ) * currentVolume;
            }

            middle = last;
            last = hull->GetEdge( last->next );
        }
    }

    centroid /= (volume * 4.0f);
    volume *= (1.0f / 6.0f);
    diag /= volume * 60.0f;
    offDiag /= volume * 120.0f;
    md->mass = density * volume;

    if ( md->mass != 0 )
    {
        md->inertia.Set(
            diag.y + diag.z, -offDiag.z     , -offDiag.y,
            -offDiag.z     , diag.x + diag.z, -offDiag.x,
            -offDiag.y     , -offDiag.x     , diag.x + diag.y
            );
        md->inertia *= md->mass;
    }
    else
        md->inertia = lmDiagonal( 0 );

    md->center = centroid;
}