yade-users team mailing list archive

[Jan Stránský <question692998@xxxxxxxxxxxxxxxxxxxxx>]

Question #692998 on Yade changed:
https://answers.launchpad.net/yade/+question/692998

    Status: Open => Answered

Jan Stránský proposed the following answer:
Thanks, now the question is finally nearly perfect* :-)
*nearly - there are "IndentationError: unexpected indent" in the code, but easy to fix

What you encounter is actually a bug in Yade [6], thanks for pointing it out.
I am not sure if the bug influences the computation with clumps or not..

Anyway, to compute true principle axes of inertia, you can "mirror" the C++ code [2]:
###
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]

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()
eigvecMin = eigvecs.col(0)
print("V",eigvecMin)
###

But, note that the bug might influence the clump computation itself (and
maybe not, I don't have time currently to investigate)

cheers
Jan

[2] https://gitlab.com/yade-dev/trunk/-/blob/master/core/Clump.cpp
[6] https://gitlab.com/yade-dev/trunk/-/issues/167

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.