Skip to content

Wrong state.ori of clumps

Consider a clump below. Its computed state.ori is wrong. (using Ubuntu 18.04 and Yade 2020-06-23.git-f03a37c0, but looking to current source code, it should be the same with the most recent version)

Does wrong state.ori influence clump computation itself or not?

O.bodies.appendClumped([
    sphere([0,0,0], radius=1),
    sphere([1,1,1], radius=1),
    sphere([2,2,2], radius=1),
    sphere([0,2,1], radius=0.5),
    sphere([2,0,1], radius=0.5),
])
clump = O.bodies[-1]

moment of inertia computed "by hand" (based on Cump::updateProperties):

def inertiaTensorTranslate(I,m,off):
    return I + m * (off.dot(off) * Matrix3.Identity - off.outer(off))
M = 0
Sg = Vector3.Zero
Ig = Matrix3.Zero
members = [O.bodies[i] for i in clump.shape.members.keys()]
for mm in members:
    dens = mm.material.density
    subState = mm.state
    sphere = mm.shape
    vol = 4./3.*pi*pow(sphere.radius,3)
    m = dens*vol
    M += m
    Sg += m*subState.pos
    Ig += inertiaTensorTranslate((Vector3.Ones*(2 / 5.) * m * pow(sphere.radius, 2)).asDiagonal(), m, -1. * subState.pos);
pos = Sg / M
Ic_orientG = inertiaTensorTranslate(Ig,-M,pos)
Ic_orientG[1, 0] = Ic_orientG[0, 1]
Ic_orientG[2, 0] = Ic_orientG[0, 2]
Ic_orientG[2, 1] = Ic_orientG[1, 2]
eigvecs,eigvals = Ic_orientG.spectralDecomposition()
#### START MISSING IN C++ SOURCE CODE
if 1: # switch to get good/wrong results
    if eigvecs.determinant() < 0:
        c1,c2,c3 = [eigvecs.col(i) for i in range(3)]
        eigvecs = Matrix3(c1,c2,-c3).transpose() # make determinant +1 instead of -1, some better approach?
#### END MISSING
ori = Quaternion(eigvecs)
ori.normalize()
inertia = eigvals
print(Ic_orientG)
r = ori.toRotationMatrix()
print(r*inertia.asDiagonal()*r.transpose())

If the determinant of eigenvector matrix is +1, then it is OK (Ic_orientG is almost equal to r*inertia.asDiagonal()*r.transpose() "composition").

However, if the if eigvecs.determinant() < 0 check is not done, the determinant of eigenvector matrix can be -1. Then, Quaternion(orthogonalMatrixButWithDeterminantMinusOne) gives nonsense orientation. Then, naturally, Ic_orientG is different from r*inertia.asDiagonal()*r.transpose() "composition".

Based on https://answers.launchpad.net/yade/+question/692998