principal axis of clump

Hello,

> the principal axis of clump

What is "the principal axis"? principle axes of inertia?
What is "clump"? clump of spheres? intersecting or not? clump of non-spherical particles?

the implementation is in core/Clump.*pp files [1,2].
The geometric properties (center of mass, inertia and orientation) is computed in Clump::updateProperties.

There are two main approaches:
- using approximate "voxelization" for the case of intersecting bodies
- using standard analytical formulas for non-intersecting case

cheers
Jan

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/core/Clump.hpp
[2] https://gitlab.com/yade-dev/trunk/-/blob/master/core/Clump.cpp

Hi,

By the way, would it be possible to elaborate about the distinction between Clump::updateProperties and Clump::updatePropertiesNonSpherical ?

@Jerome: of course :-) actually there are two Clump::updatePropertiesNonSpherical methods. I just did not want to go much in detail (mainly because I do not know the details myself :-D before the question is focused enough.
Another option is to start a separate question on this updatePropertiesNonSpherical topic (or a discussion on yade-dev mailing list or gitlab issue?).
cheers
Jan

Hi,
   sorry for reply late..I mean the clump of spheres and principle axes of inertia..did yade have any functions to calculate it directly??

Thanks!

> I wonder how yade calculate the principal axis of clump? is there any functions or any sourse codes?

there are functions in the source code [1,2]

> did yade have any functions to calculate it directly??

Yade does have functions to calculate it directly, see [1,2]

cheers
Jan

PS: to get a better answer, please read [3] and provide a good question

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/core/Clump.hpp
[2] https://gitlab.com/yade-dev/trunk/-/blob/master/core/Clump.cpp
[3] https://www.yade-dem.org/wiki/Howtoask

hey, ??
  sorry for reply late..I mean the orientation long axis of clump of spheres..did the function in yade is i.state.ori??
Thanks!

Hello,

> hey, ??

what is this salutation?

> I mean the orientation long axis

What is "orientation long axis"?

It is still the same. Please read [3].
Then read [3] again.

Then, to get a better answer, please read [3] once more and provide a good question (as requested in #5). I.e.:
- explain, how/why the previous answers did not help.
- provide requested information. e.g. is the particles in the clump intersecting, or not - asked in #1.
- be as clear as possible. What is "orientation long axis"?? orientation OF long axis? Then still "long axis" is nor clear. Do you have some reference / definition?
- provide a MWE [3] explaining what you want to achieve - illustrating "long axis", illustrating a typical investigated clump (spheres with overlaps or not, ...), ideally also expected results, ...

You can compute all the clump-related values:
- by hand (after all it is "just" not too difficult math)
- by your own function(s) ("automated" by hand)
- by Yade predefined function(s)

We can help you with all three approaches, but first we need to know with what you want to help and where is the problem (which is not the case of this thread).
If you want us to spend time on your problem, you spend time first, at least to formulate reasonable, understandable and clear question..

cheers
Jan

[3] https://www.yade-dem.org/wiki/Howtoask

Hello,
    sorry for didn't ask the questions clearly.. I am conducting DEM of flat sands with clumps consists of several spheres, some spheres are intersecting. I have checked yade and found that the function "inertia" could calculate the moment of inertia in the space of principal axes, but didn't give the orientation of principle axes .
    And my question is : how can we know the orientation of the longest axis of clump?? because sands are always flat..did the orientation of longest axis have some relation with the above principal axes of moment of inertia?

Hello,

> how can we know the orientation of the longest axis of clump??

how can we know what you mean by "the longest axis of clump"??

> did the orientation of longest axis have some relation with the above principal axes of moment of inertia?

It depends on what "longest axis" is.

cheers
Jan

Hello,
          the longest axis of clump is defined as the longest line segment that links two points on the edge of the particle. Its inclination angle (0–180°) in the plane is defined as the long-axis direction..

Thanks

You have clump of spheres.
Does such clump has any edge (I would say no)? If yes, how you define such edge?
If it does not have edge, how a line segment can link point on the edge?
Did you mean surface instead of edge?

Anyway, such geometric characteristics are not directly related to inertia.
You will have to use another approach. I am not sure if Yade has predefined tool for this (I would say it does not)..

cheers
Jan

hey,
    sorry, I mean on the surface, and how yade calculate the orientation of principle axes( moment of inertia ) of clump and what dose it mean??

thanks

> how yade calculate the orientation of principle axes( moment of inertia )

as eigenvectors of inertia tensor, see [1,2] for the implementation.

> and what dose it mean??

what is "it"?

cheers
Jan

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/core/Clump.hpp
[2] https://gitlab.com/yade-dev/trunk/-/blob/master/core/Clump.cpp

hey
    What does the eigenvector corresponding to the maximum moment of inertia of clump mean? and Is this eigenvector independent of the global coordinate system or the local coordinate system?

Thanks

> What does the eigenvector corresponding to the maximum moment of inertia of clump mean?

It is the orientation (unit direction vector) of principle axis corresponding to the maximum moment of inertia

Note that, the principle axis corresponding to minimum (not maximum) moment of inertia would normally be closest to the "longest axis".

> and Is this eigenvector independent of the global coordinate system or the local coordinate system?

It is a vector. As such, it is independent of any coordinate system.

But once again note that principle axes of inertia has no direct relation to the "longest axis".

cheers
Jan

sure,but seems function "state.inertia" could only calculate the principle moment of inertia, and is there any functions which could calculate the inertia tensor??

thanks

> is there any functions which could calculate the inertia tensor??

not directly, but e.g a matrix representation of inertia tensor could be obtained by
inertiaTensor = state.inertia.asDiagonal()
or
r = state.ori.toRotationMatrix()
inertiaTensor = r*state.inertia.asDiagonal()*r.transpose()
or
r = state.refOri.toRotationMatrix()
inertiaTensor = r*state.inertia.asDiagonal()*r.transpose()

Which one to use depends on circumstances.

cheers
Jan

I had a look in documentation. All we discussed (expect long axis) is there.

How inertia is treated for non-spherical particles (including clumps) is explained at [4].
How properties of clumps is evaluated is explained at [5] (and of course in [1,2]).
Please read it first.

Also try to be consistent and clear. Almost every new post, you are jumping to different topics:

- principle axes of inertia
- long(est) axis
- inertia tensor

- source code
- functions for direct calculations
- what something does mean
- how yade does something
- ...

Every second post is unclear and needs clarification..

It is then really difficult to give some reasonable answer..

If you want to continue, please:
- read [3] again and start to follow the suggestions ("Try to describe your problem in a concise and accurate way", "Do not merge different questions in one single post", "Train your empathy on this example")
- provide a good and clear question (requested in #5)
- provide an example (requested in #7), sample input and expected output. For this case, an example would solve all the confusion, something like:
"""
I have this simple artificial clump. I want to compute ......, the result should be ......
###
s1 = sphere((0,0,0),1)
s2 = sphere((1,1,1),2)
O.bodies.append((s1,s2))
clump = O.bodies[2]
# the result should be the inertia matrix Matrix3(......), how to do it in Yade?
"""

cheers
Jan

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/core/Clump.hpp
[2] https://gitlab.com/yade-dev/trunk/-/blob/master/core/Clump.cpp
[3] https://www.yade-dem.org/wiki/Howtoask
[4] https://yade-dem.org/doc/formulation.html#orientation-aspherical
[5] https://yade-dem.org/doc/formulation.html#clumps-rigid-aggregates

hello:
       very thanks for your advice, I choose this simple artificial clump consists of three particles. I want to compute principle axes of minimum moment of inertia, I think the result should be (0.57735,0.57735,0.57735), but I always get (-0.80473785, -0.50587936,0.31061722)..could you help check it?

Here are the codes:

from __future__ import print_function
import numpy as np
from builtins import range
from yade import pack,export,qt
from numpy import linalg as LA
import math

#define material for all bodies:
id_Mat=O.materials.append(FrictMat(young=1e6,poisson=0.3,density=1000,frictionAngle=1))
Mat=O.materials[id_Mat]

#define engines:
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 NewtonIntegrator(damping=0.7,gravity=[0,0,-10])
]

#### show how to use appendClumped():

#create 3 clumps:
clump1=O.bodies.appendClumped([\
sphere([0,0,0], material=Mat, radius=1),\
sphere([1,1,1], material=Mat, radius=1),\
sphere([2,2,2], material=Mat, radius=1)\
])

#definition for getting informations from all clumps:
def getClumpInfo():
 for b in O.bodies:
  if b.isClump:
   print('Clump ',b.id,' has following members:')
   keys = list(b.shape.members.keys())
   for ii in range(0,len(keys)):
    print('- Body ',keys[ii])
   r = b.state.ori.toRotationMatrix()
   inertiaTensor = r*b.state.inertia.asDiagonal()*r.transpose()
     w,v = LA.eig(inertiaTensor)
     emax=max(abs(w[0]),abs(w[1]),abs(w[2]))
   print(inertiaTensor)
   print(w)
     print(v)
     emin=min(abs(w[0]),abs(w[1]),abs(w[2]))
     o=np.zeros(3)
   for j in range(0,3):
     if abs(w[j])==emin:
       o[0]=v.item((0,j))
       print(v.item((0,j)))
       o[1]=v.item((1,j))
       o[2]=v.item((2,j))
   print(o)
                        print('inertia:',b.state.inertia)
   print('mass:',b.state.mass,'\n')
   print('ori:',b.state.ori,'\n')

getClumpInfo()
renderer = qt.Renderer()
qt.View()

Thanks

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

Considering the "longest axis", I **think** (am not sure, you can test it), that if you check all sphere pairs sph1,sph2 and take the pair with maximum
length = (sph1.state.pos-sph2.state.pos).norm() + sph1.shape.radius + sph2.shape.radius
you then have the maximum clump size (quickly I cannot imagine a longer surface connection) and the axis is given by sph1.state.pos-sph2.state.pos

sph1.state.pos-sph2.state.pos is vector connecting sphere centers
(sph1.state.pos-sph2.state.pos).norm() is distance of centers

cheers
Jan