yade-dev team mailing list archive

[Václav S(milauer <eudoxos@xxxxxxxx>]

====
Again, without attachments. Find them here:
http://beta.arcig.cz/~eudoxos/yade/shear-displacement.pdf and
http://beta.arcig.cz/~eudoxos/yade/shear-slip.pdf
===
> (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.
slipToEpsTMax(..) doesn't add any data to SpheresContactGeometry. It
merely adjusts initial contact points (which is really their relative
orientations WRT each sphere).

I hope it will be clear from the attached figure. c1,c2 are sphere
centers, contact plane is perpendicular to line connecting those,
passing through the contact point. p_01 and p_02 are initial contact
points on sphere1 and sphere2; because spheres rotate (and move), these
points go apart. To compute shear displacement, you have to get the
length and direction of their trajectories (-.- on the figure), which is
done by projecting (dotted lines; unroll* method) these arcs to the
contact plane. That gives you projected initial contact points p'_01 and
p'_02 and their difference is the shear displacement u_T (a Vector3, but
always in the plane; oriented in global coords, origin shifted to
contact point).

In the case I get u_T outside the plasticity surface, I can say what is
the maximum |u_T| in that direction. In that case, p'_01 and p'_02 will
be shifted one towards another, such that the |u_T| is as desired, while
keeping its direction (radial stress return, non-associated law, since
it doesn't affect normal strain). p_01 and p_02 are adjusted such that
this change is permament. This way, plasticity adds no extra data to the
geometry, it merely changes initial contact points. Also, I don't need
to remember the history, as I am always within/on the plastic surface.
>> 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...
True, but I don't care about 3 extra numbers per interaction. As said, I
can separate my code completely, it is just pain in the neck to have too
much classes, some of them with just 1 user; why not, Yade is a
platform, after all...
> 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.
Duplicating data is not bad in itself; always the situation matters. For
SpheresContactGeometry, for example, it is IS2IS4SCG that should update
all that, no user intervention at all. If you fiddle with that, however,
you should know what you're doing. Some data will always be duplicated;
I think that locality makes it possible for the compiler to better
optimize the code. And debugging is also easier, of course.
> This incremental formulation needs a serious discussion.
Agreed. Some can derive those formulas and post a PDF with that and a
few images?
> I'm not really convinced by the global scale comparison with FEM. FEM
> nodes have rotational DOF's and radius?
Why not? In this model, they have 6 DOFs. Nodes (at sphere centers)
don't have radii, but the elements ("interactions") are generated by
yade during the first iteration and then exported with all their
parameters to oofem.
> 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".
I hope the figure makes it clear. To me it seems that this formulation
is that self-explanatory that I don't see a simpler way (as far as the
idea goes) to compute it. We can try comparing these two. They should
really be the same. What makes thing complicate is that interaction
don't have local coordinate system, therefore everything has to be
transformed at every iteration.
> 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).
True, sorry about that. I suspect the approximated method that it will
accumulate errors over time. Besides that, it is not necessarily true
that angularVelocity*dt=rotation --
> 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?
Attached shear-slip.pdf. There is the old u_T (from the first figure)
and I request its maximum length (so that it is on the plasticity
surface); projected points (original contact points in the contact
plane) are moved in that way, then original contact points on spheres
are moved as well.
> - Not related but still interesting :
> if(displacementTMax<=Mathr::ZERO_TOLERANCE), why not if(dMax<=0)?
My fault. I wasn't sure, I corrected that now.
> 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...)?
Of course I agree. I project radius*pi/4 to the contact plane (well, not
really radius, but the intial distance from sphere center to contact
point) in the right direction. I get radius*pi/4 exactly.
>> 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.
So someone go ahead and factor new function out of that (I don't use
that code and the idea of debugging 20 lines of ununderstandable
formulas makes me sick already); it needs just linear and angular
velocities, positions, contact point, normal, dt, isDynamic1, isDynamic2
for parameters, as far as I can see.

> 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.
Not one for each contact law (many need same geometry info -- then you
would need also
InteractingThing2InteractingThing4YetAnotherContactGeometry for every
constitutive law!!), but otherwise I agree.
> Hehe. You started the "abusive storage" and others filled it with
> their stuff.
> You are only "partly" to blame then. ;-)
No, I just _marked_ that thing "abusive storage". It was there before I
touched that code IIRC. I don't know if Vincent really uses that...
Vincent, ping?

>> 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.
That's very true, I will work on that.
> 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.
I keep that in mind. Janek, what do you need for your law? Normal displ,
shear displ, bending and torsion is enough?
> 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...)
I don't get the point here, sorry.
> "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.
Clumps don't interact, only their constitutents do. For snow, I never
looked.
> 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.
We could have list of interaction within body, but what are you going to
store? Shared_ptr's wouldn't work because of broken serialization. You
could store id of the other body, but it would need lookup in the
interaction container, but it can be quite fast. The list would be
maintained just by the collider, since only collider may add or delete
interactions. That needs to be fixed in some of the constitutive law, too.

Vaclav