yade-dev team mailing list archive

[Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>]

> You can always derive from the geometry class and add your own data
> members. But you still have access to inherited data members and
> methods, should you need those. maxUn can (should) be stored in
> InteractionPhysics, the same for cell volume, no?
>   
(A) - maxUn (and other extrema of the variables contained in IntGeom, 
including initial orientations) in IntPhysics is ok with me. It will 
prevent the addition of excess parameters in IntGeometry.


> I have SpheresContactGeometry::splitToEpsTMax(Real epsTMax),
> which adjusts the original contact points such that |epsT|<epsTMax and
> returns increment of eps_plastic. I believe that most constitutive laws
> use similar notions of plasticity, so to me it made sense to add this to
> SpheresContactGeometry.
>   
You can't define this in ContGeom. This is in contradiction with (A), 
and also in contradiction with the fact the plastic part of a strain is 
defined for a given constitutive law.


> > An exemple of this : in creep equations (creep, like plasticity, is
> > changing the "neutral" position when it occures), the initial
> > orientations (from wich total bending/twisting are defined) are
> > updated to reflect creep. The way it is updated is valid only for the
> > specific creep law considered here, another constitutive law would
> > update the initial configuration differently.
> >
> > The conclusion is you can't avoid the manipulation/storage of
> > geometrical data in some constitutive laws/physicalParameters.
> >     
> Constitutive law can change geometry (as with the plasticity), in some
> special cases. I know that I say something more than just "compute
> stresses from strains". I would say that what is useful for multiple
> constitutive laws could be defined as SpheresContactGeometry method
> (plasticity), otherwise just that single constitutive law can manipulate
> that.
>   
Again, usefull for me is : normal, Un, dUt, rad1, rad2. I don't need 
(want) anything more in SpheresContact description.
It makes code harder to read and it wastes memory.

> If you say incremental shear, you still need access to 2 linear
> velocities and 2 angular velocities, not just penetrationDepth and radii.
>
>   
Well, dUt could be added in IntGeometry. But YOU don't need it...

rad1/rad2 should also be removed, this duplication caused troubles at 
least once when I was increasing the size of particles.
I was multiplying each body->radius by a constant, but I forgot to 
multiply each interaction->rad1/2 at the same time.
Spheres where increasing in size but contacts were computed based on the 
initial sizes. Consequence : no force were generated since the only 
"motion " was the radius evolution and it was not seen at the contact level.
This is a major drawback of the "let us copy what we need from bodies to 
interaction geometry" approach.

If the objective is only "good design", I'm not sure it is worth it.

Duplicating data is bad in itself. By trying to avoid accessing bodies 
in constitutive laws, we are generating duplications and extra 
members/local variables. The variable is somewhere, let us find it where 
and when we need it. Forecasting what people do/will/could need in 
constitutive laws is almost impossible.


> I personally don't trust the incremental code, although I admit that it
> can be my personal weakness that I wasn't able to derive it and verify
> myself. Besides that, it supposes that the rotation is such small that
> displacements can be linearized (if I understand what it does) and I
> think that could make for some problems in highly dynamic simulations.
>
>   
This incremental formulation needs a serious discussion.
I'm not really convinced by the global scale comparison with FEM. FEM 
nodes have rotational DOF's and radius?
Looking at your code, I suspect you wrote a definition of Ut which is 
not the classical one (so that the increment of your Ut is not my dUt). 
It is not "more precise", it is just "something else".
You are perhaps also mixing the "incremental formulation" (no assumption 
on small rotations) with the approximated method to rotate current Ft in 
Yade (which indeed needs this assumption).

For rotating Ft. In short : in a rotation, the exact displacement of P  
is dP=(M-Id).P. the operator "(M-Id)." can be replaced by "R x ()" when 
the rotation is small.  With R aligned on the rotation axis and scaled 
by the rotation angle, and "x" the vectorial cross product.

I've been reading "slipToDisplacementTMax" code, but I have now more 
questions than answers.
I still can't understand what it is supposed to do. I guess this is 
because we have very different situations in mind.

- How can you update the "initial" position based on a real instead of a 
3D vector?
- Not related but still interesting : 
if(displacementTMax<=Mathr::ZERO_TOLERANCE), why not if(dMax<=0)?

Even the comment seems obscure to me :
/*! Perform slip of the projected contact points so that their distance 
becomes equal (or remains smaller) than the given one.
* The slipped distance is returned.
*/

Why do you project material points?!

Quick check :

A sphere is in contact with a plane. The plane and the center of the 
sphere are fixed, the sphere can only rotate with an axis of rotation 
parallel to the plane.
IF the rotation of the sphere is pi/4, what is the tangential displacement?

I'd say Ut here is _exactly_  (radius x pi/4). Do you agree (I 
anticipate you don't, with this "projection" feature...)?
In Yade, the incremental formulation gives exactly the same result. I 
mean, if you take the algorithm and write it in a mathematical form, you 
can proove that if the rotation pi/4 is discretised in N small steps, 
the sum of dUt's is exactly radius x pi/4.


> If you think it is realiable and useful, I would be all for extracting
> it to a stand-alone function that would get all the parameters it needs
> and return what it computes.
It is not just "reliable" : it is what most people knowing DEM will look 
for in Yade, it is the very root of the DEM. Really.
So, yes, it will need to be adapted to the new design.
And if there is an error in the code, I need to spot it.
>  Or shear strain could be added to
> SpheresContactGeometry (or a new class, if you like), but computed by
> InteractingSphere2InteractingSphere4SpheresContactGeometry. That would
> make more sense to me.
>   
Sorry to insist, but, what I'd like ideally is SpheresContactGeometry 
with Un, dUt, normal (dUt is actually not there, perhaps it should, 
rad1/rad2 are there - they should not).
Anything else (rotation, etc.) should be in other classes inheriting 
from SpheresContactGeometry, one for each contact law.

I'd like to keep the basic and common Yade classes very simple, with 
nothing more than the classical dem. Same for IntPhysics : kn, ks, 
friction. Nothing more (this useless "intitials" kn/ks caused troubles 
in the past, a bit like the duplicated radii).

> Oh, I didn't add those (except facetContactFace, which I'd like to keep
> until SphresFacetGeometry is done). "svn blame" gives me "richefeu". But
> I can remove that, no problem at all ;-)
>   
Hehe. You started the "abusive storage" and others filled it with their 
stuff.
You are only "partly" to blame then. ;-)

> I put those extra things under condition 'hasShear', they are not
> assigned and should "just" take memory if you don't use them (29 bytes,
> if I calculate correctly: 5 quaternions, 2 Vector3, 3 Reals).
Ok, all this instead of 2 vectors and 1 real makes, what? 4 times more?
It reminds me the question of Luc about the size of xml's... 70% of an 
xml is for interactions...

+ It is harder to read.

>  If that is
> not OK and if other peoples will not use it, I can move it to a separate
> class.
Yes please, move it in a new inheriting class.

Currently, the rotation variables for CohesiveFrictional are in 
CFInteractionPhysics.hpp, because it was simpler for Janek (see below). 
They are very similar to yours. If there is a new geometry class to put 
them, it is all good.

Currently in the CFI.hpp :


/Quaternionr    initialOrientation1,initialOrientation2,
               orientationToContact1,orientationToContact2,
               currentContactOrientation,initialContactOrientation,
               twistCreep;
Vector3r    initialPosition1,initialPosition2; /

I think initial's should not move. They would store in IPhys. the 
extrema of some IGeom. variables (still, I wonder if kinematic variables 
of contact laws will always be extrema of standard geometric variables. 
Wait and see...)

> It is already stored there (iff hasShear): pos1, pos2, ori1, ori2. But I
> have yet to see non-spherical particles in yade (if someone is serious
> about that, he'd better fix global/local mismatch of inertia in Newton's
> law first)
>   
"Yet to see"?! Open your eyes then! :-P
Janek has been developping non-spherical grains for a while. I did not 
invent this position example, I got it from him.
Your clumps are a great exemple too, even with some mismatches left.

> > So, fundamentally, with you first statement in mind : "Constitutive
> > law should operate just on interaction geometry and interaction
> > physics. If it accesses rootBody->bodies, it is a design problem", all
> > information from bodies must be duplicated in interactionGeometry,
> > which is a serious design problem too...
> >     
> I know it is a hassle, it could be that SpheresContactGeometry would
> contain shared_ptr<Se3> shared with Body::physicalParameters, but
> current serialization wouldn't handle that correctly (it would create 2
> copies of the same Se3 instead of the shared one). I still have the
> parallel stuff in mind and it is advantageous if the computation is as
> local as possible (i.e. no access to rootBody->bodies and such).
>
>   
Yes, I was suspecting potential conflicts with parallelisation... so ok.

Appart from that, yes, it is obvious that "interactions" should link to 
the particles in interaction.
In PFC, you can get the pointer to grains from the pointer to an 
interaction, and you can get the pointer to contacts from the pointer to 
a particle. This is very usefull to write a number of algorithms, and I 
always wondered why we don't have this double-linked lists in Yade. 
Parallelisation is one reason now.

Bruno






--

_______________
Chareyre Bruno
Maitre de conference

Grenoble INP
Laboratoire 3SR - bureau E145
BP 53 - 38041, Grenoble cedex 9 - France
Tél : 33 4 56 52 86 21
Fax : 33 4 76 82 70 43
________________