Asked by kelaogui on 2013-06-19


I would like to make some modification to updateProperties function in clump.cpp line 99. First I wonder where the function is called. It seems not in NewtonIntegrator.cpp.

I guess if this function works, the properties of the clump, including inertia and mass, will be computated in each step based on clump members' (sphere) mass and position. But I found if I assigned the inertia and mass directly to the clump using O.bodies[clumpid].state.inertia=xxx, it keeps constant in later steps, won't be updated via sphere mass and position. Does that mean we can just give a right mass and inertia to the clump and it will be used in the movement computation always?

Question information

English Edit question
Yade Edit question
No assignee Edit question
Solved by:
Last query:
Last reply:
Christian Jakob (jakob-ifgt) said : #1

I do not really understand, why you want to compute inertia "in each step" ...
Inertia of a body is constant since mass (or shape) does not change. So it is calculated once after creating or modifiing a clump*. There is no need to calculate it in each step.

Inertia tensor is actually just approximated in Yade [1]. It is simply the sum of all inertia tensors of clump members. If you need a more exactly inertia tensor, you can try to improve the section between line 121 [2] and 124. For this you can have a look how you can get inertia tensor by integration [3]. I just found a german entry in wiki (no english, sorry) about inertia tensor.

Feel free if you want to improve this section.




*Clump::updateProperties is called by following methods:
- clump()
- appendClumped()
- replaceByClumps()
- addToClump()
- releaseFromClump()
- adaptClumpMasses()

kelaogui (kelaogui91) said : #2

Oh, yes. I think I made a mistake in understanding the whole thing. Now I get your idea.
Thanks very much, Christian.

Christian, what about the Monte-Carlo method?
Isn't it a better approximation that sum of members inertia?

kelaogui (kelaogui91) said : #4

I think Monte Carlo is a good way. But for 3d clump and large sample, the computation consuming will be a problem and also, it's still an approximation,right?

Christian Jakob (jakob-ifgt) said : #5

@bruno: yes, it should basically look like getClumpVolume()

@kelaogui: Of course it is time consuming, but it is just calculated once. Also the user can choose the accurateness by number of points for integration/summation. The higher numTries, the better the result but also the longer the calculation ...

Jan Stránský (honzik) said : #6

It depends on the fact, if the bodies of one clump are overlapping or not.
If they overlap, the Monte Carlo method is the easiest way from
implementation point of view. On the other hand, If they do not overlap,
than the "sum of members inertia" is the fastest as well as exact approach.

2013/6/19 Christian Jakob <email address hidden>

> Question #230986 on Yade changed:
> Christian Jakob posted a new comment:
> @bruno: yes, it should basically look like getClumpVolume()
> @kelaogui: Of course it is time consuming, but it is just calculated
> once. Also the user can choose the accurateness by number of points for
> integration/summation. The higher numTries, the better the result but
> also the longer the calculation ...
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
> _______________________________________________
> Mailing list:
> Post to : <email address hidden>
> Unsubscribe :
> More help :

> yes, it should basically look like getClumpVolume()

My bad, sorry. I thought it was already implemented for inertia to.
Probably only one more line in the MC algorithm, something like

V = V + 1
+ I = I + 1*dist*dist*density

where "*" is diadic product?

Christian Jakob (jakob-ifgt) said : #8

already working on that ;)

i think dyadic product is not correct, please have a look at equation 2:

Still one line: :)

I = I + (dist.dist x identity - dist*dist)*density

where "." is dot, "x" is scalar*matrix, and "*" is dyadic!

Václav Šmilauer (eudoxos) said : #10

Hi, coming accross this discussion... Regular grid sampling is not Monte-Carlo (random sampling). For the dyadic product, you can write vector*vector.transpose() for dyadic (outer) product in Eigen. Cheers!

Christian Jakob (jakob-ifgt) said : #11

thanks bruno and vaclav, i think i have got it (not tested yet).
but it also leads to some changes in wrapper and examples since adaptClumpVolume() is no more needed.

Václav Šmilauer (eudoxos) said : #13

BTW someone can derive a formula for critical timestep of a clump?

As it is currently, the time stepper computes the critical timestep for clumps. The formula is the same as for standalone bodies.
There is still a small thing wrong in the rotational stiffness, as it should not be simply the sum of RStiffness of each member but should be translated to the clump's center (just like rotational inertia changes when expressed wrt a different point).

Well done!

Václav Šmilauer (eudoxos) said : #18

For the formula for dt: I guess you can be way off. "They are only approximated assuming than DOF’s are uncoupled" - they are not uncoupled in clumps at all, on the contrary particles in the clump are linked with (quasi-)infinite stiffness. So I think the whole clump needs to be considered, not just individual particles.

The whole clump is indeed considered as a unique particle with a total mass and stiffness.
Uncoupled DOFs means that a displacement in the direction x generates a dominant force along x (i.e. the generalized stiffness matrix is assumed diagonal). Be it for a clump or a sphere, this approximation is the same.
It is possible to compute without this assumption: compute eigen values of a 6x6 matrix for each body instead of using the diagonal values.
It would maybe let a higher coefficient be used: 0.9dt or so instead of 0.8dt, if you think it is worth it.

Václav Šmilauer (eudoxos) said : #20

@Bruno: my bad, I did not read the code carefully, it is going through all clumped particles (I though they were handled separately, in which case x-force on one particle would generate dominantly x-force on another particle - that's what I meant).