Polyhedra contact instability (potential Bug?)

Hi,
Thanks for your report (I do not confirm, didn't check details).

> I wish, I could be able to verify this by myself, but it is not my league

You don't seem far from it. If you have been able to spot the potential issue you should be able to edit/recompile - this is mainly about changing some algebraic expressions (ok, a bit more than that ;) ). That would be a perfect training.

So if no one pick up the task soon enough I suggest you give it a try. We can help.
I opened an issue on gitlab to track progress. [1]

Bruno

[1] https://gitlab.com/yade-dev/trunk/-/issues/192

Thank you for your encouragement! I have to admit there is a long way for me before I would be able to 'give it a try'. The last time, I compiled something, I was in high school (half of my lifetime ago). And I only started using Ubuntu, because of Yade. I started learning C++ from basics recently, and I am on 'while loop' right now ;)
I do not intend to learn everything before I get my hands on the code, but it is still a matter of weeks or even months. If the issue that I raised will still be around at that time, I will try to tackle this.

Hello,

I can confirm this problem. The contact normal together with contact force magnitude might be the source of the problem (I was not solving the problem, just detecting it).

Apart from Bruno's option, have a try potential blocks [1]. It models the same shape, but different contact logic. The results seems more "natural" for these cases. Have look at the testing code below.

Disclainmer: I have ever used potential blocks for this testing script only :-) it was some time ago for my colleague playing with polyhedrons, running into the same problems you describe.

cheers
Jan

[1] http://yade-dem.org/doc/potentialparticles.html

###
O.engines=[
    ForceResetter(),
    InsertionSortCollider([PotentialBlock2AABB()]),
    InteractionLoop(
        [Ig2_PB_PB_ScGeom(unitWidth2D=1.0, calContactArea=True)],
        [Ip2_FrictMat_FrictMat_KnKsPBPhys(kn_i=1e8, ks_i=1e7, Knormal=1e8, Kshear=1e7, viscousDamping=0.2)],
        [Law2_SCG_KnKsPBPhys_KnKsPBLaw(label='law', neverErase=False, allowViscousAttraction=False)]
    ),
    NewtonIntegrator(damping=0.2, exactAsphericalRot=True),
    PotentialBlockVTKRecorder(fileName='/tmp/pb', iterPeriod=50, sampleX=60, sampleY=60, sampleZ=60, label='vtkRecorder')
]

def cube(pos,size):
    ret = Body()
    r = .5*size
    a = [1,-1,0,0,0,0]
    b = [0,0,1,-1,0,0]
    c = [0,0,0,0,1,-1]
    d = [.5*size-r for _ in range(6)]
    minmaxAabb = 1.05*Vector3( d[0], d[2], d[4] )
    ret.shape=PotentialBlock(r=r, R=0.0, a=a, b=b, c=c, d=d, id=len(O.bodies))
    utils._commonBodySetup(ret,ret.shape.volume,ret.shape.inertia,material=-1,pos=pos)
    ret.state.ori = ret.shape.orientation
    return ret

b1 = cube((0,0,0),1)
b2 = cube((2,0,0),1)
O.bodies.append((b1,b2))
b1.state.vel = (1,0,0)
b1.state.angVel = (0,0,2)
O.dt = 1e-3
O.run(3000,True)
###

Hi Jan,

thank you for this example. This can be a good starting point for playing with potential blocks. Right now, I am not able to run it because my version of Yade doesn't have 'misc.py' and 'PotentialBlock2AABB()'. If I knew it is the only problem, I would dig into the Potential Blocks, but I want to split those polyhedra later.

Referring to the original problem of normal force direction, I think that I have a conceptual solution. I do not have any proof, it's kind of intuitive.
In general, all the systems tend to have minimum energy. -> So the contact force should act in the direction that causes a certain force change (volume change) with minimum displacement. -> Movement should be perpendicular to 'the widest cross-section'. -> The projection of hull intersection on the sought plane (normal to the force) should give a polygon of the biggest area.

Implementing it in Yade is out of my reach, but maybe I will test it in Python (for geometry from my example).

Best wishes,
Karol

> Right now, I am not able to run it

I have tried Ubuntu 18.04 and yadedaily and it worked.

Concerning the biggest area polygon, I would start with some simple "testing" 2D situations.

Consider a "ground", large enough polyhedron with xy surface.

Consider a prism coming to ground. In the beginning, The intersection is "flat square", the largest section is perpendicular to z axis, as expected.
    | |
__|______|__
    |______|

Consider a "needle", a very sharp pyramid or a "long" prism with z axis. Then the intersection is a subshape with z axis and the contact normal is expected to be also z axis.
But it is the direction with **least** section area..
The largest area is perpendicular to z axis..
__\__/___
      \/

___|_|___
      | |
      |_|

So it seems the biggest area polygon is not universally good..
BUT, it just a quick brainstorming, maybe there is some mistake in my assumptions or there exist some easy fix :-)

cheers
Jan

Hi to everyone,

I think you might find two recent papers of Prof. Feng useful [1,2], as he proposes that volumetric contact law is the way to go, having energy conservation in mind, while a calculation for the contact normals is also proposed.

In the Potential Blocks we have stable/reasonable contact normal values for face-to-face contacts, but personally, I consider a volumetric law (and volume-based normals) much more accurate from an energy point of view, for complex contact geometries (like face-to-vertex), if the contact normals are somehow made more stable.

Hope this helps,
Vasileios

[1] https://www.sciencedirect.com/science/article/pii/S0045782520306782#fig6
[2] https://www.sciencedirect.com/science/article/pii/S0045782520306393

Vasileios,

thank you for providing those top-notch papers! They are very general and formalized, so I have only skim-read them. My intuitive idea seems to be similar to the solution presented in the paper of Prof. Feng [1]: "(...) The above conclusion has a sound physical explanation: if the current contact state is fully described by a contact tenergy function w, the true contact force and moment at the state will reduce the contact energy most effectively or at the greatest rate, i.e. along the negative gradient direction of w."

The main problem is the calculation of the contact energy function - w. As I said, his approach is more general and can be applied to 'arbitrarily shaped particles'. Still, I think that the 'maximum area polygon approach' presented above, can be a solution of a special case: convex-shaped polyhedra :)

I also agree that if only contacts were more stable, the polyhedra would be unbeatable for many applications.

Jan,

I do not know if the proposed solution can be easily interpreted in 2D. However, we can benefit from the simplicity of your 2D sketches and stay in 3D space, by assuming axisymmetry. In all of the proposed cases, the hull intersection of the contact will be a circle. This circle projected on any plane other than horizontal would be an ellipse with an area smaller than the area of the circle. So, it leaves only one choice of the normal contact direction :)

Best wishes,
Karol

[1] https://www.sciencedirect.com/science/article/pii/S0045782520306393

I just wanted to comment that, I don't deny the conceptual idea presented by Eliáš [1]. Perhaps, I talk too much about a potential new solution and draw attention back from the problem stated in the first post, that there is probably an implementation issue.

Best wishes,
Karol

[1] Eliáš, J. (2014). Simulation of railway ballast using crushable polyhedral particles. Powder Technology, 264, 458-465.

Dear all,

Vasileios directed me to this conversation, because I used Eliáš' model in the past and he remembered once I told him that I'm using and older Ubuntu version for this particular model since I experienced stability issues with newer versions. This version is Ubuntu 16.04 xenial with yadedaily (2018.02b-28948d72e6~xenial). Older Yade versions doesn't support every utility for the polyhedral model I need and the newer ones are unstable, as Karol described.

I've tried Karol's script and I can confirm that there's a sudden change in contact forces in case of the newer versions (even with Yade 2018.02b on Ubuntu 18.04). However, Ubuntu 16.04 with 2018.02b works completely fine and the contact forces are smooth.

I also had a conversation with Jan Eliáš and he pointed out that the bug might only be in Yade, but in CGAL as well, as the Ubuntu versions use different versions of the library.

I hope I could help and sorry for not letting you these information sooner,
Ákos

I've also added a comment with figures to the issue on gitlab, confirming the bug. [1]

Ákos

[1] https://gitlab.com/yade-dev/trunk/-/issues/192

Dear Ákos,

thank you for adding the description of the problem o the Gitlab and your help. Taking into account all your advice, including Bruno's encouragement for learning by practice and Jan's proposal for searching for an easy fix, I did a crazy attempt.

I've learned how to compile Yade and just tried different small options and modifications (pure alchemy, not really knowing what I am doing). One of my first guesses led to the potential "solution".
I have changed the Dimension_tag parameter in line 515 of the Polyhedra_support.cpp (dev trunk) [1] from '1' to '0'.

"linear_least_squares_fitting_3(segments.begin(), segments.end(), fit, CGAL::Dimension_tag<0>());"

Now I get results similar to yours and way more stable simulation in general. I think that may be the solution, but please don't ask me why it works :)

Best wishes,
Karol

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/Polyhedra_support.cpp

Hi Karol,

I've tried your solution and can confirm that it works, good find. However I'm not perfectly sure that this method works for every simulation on long term, but it is currently out of my knowledge as I'm also mainly on the user side. I've checked a few source code of older versions of Yade and this flag was always set to 1. Although it makes sense that if the forces change their direction suddenly, the problem might be related to the change of contact normals.

Ákos