ESyS-Particle

Possible identification of Yucang's new implementation bug

Asked by Diego Peinado Martín

Hello, I do not know if I should use the bug report tool or this is the adecuate channel.

Reading the code of RotFricInteraction.cpp, I've seen that a part of code that do not compile regarding a new implementation can have a bug.

The piece of code is (I will insert comments in the code)

/**
 * Yucang's new implementation which takes into
 * account rigid body rotaton of two particles.
 * Unit tests failing for this implementation.
 */
void CRotFrictionInteraction::calcForces()
{ //cout << "wyc in RotFric:: calcF " <<endl;
  Vec3 pos;
  Vec3 force;
  Vec3 dv,ds ;
  Vec3 d_Ffric ;
  // calculate distance
//**********************************************
  Vec3 D=m_p2->getPos()-m_p1->getPos();
//**** this line is different than in other modules. It defines D as a vector that goes from particle 1 to particle 2
//**** In other parts including the working implementation is defined as the opposite: particle 2 to particle 1
  double dist=D*D;
  double eq_dist=m_p1->getRad()+m_p2->getRad();
  // check if there is contact
  if(dist<(eq_dist*eq_dist)){ // contact -> calculate forces
    //--- elastic force ---
    dist=sqrt(dist);
    force=D*(m_k*(dist-eq_dist)/dist);
    m_normal_force=force;
//*********************************************
    pos=m_p2->getPos()+(m_p2->getRad()/eq_dist)*D;
//** Here is the problem ... as D goes from 1 to 2, the vector D/eq_dist is a unit vector from 1 to 2, so the vector
//** pos points to one point in the opposite part of the sphere than the contact.

Then as the rest of the code seems to work like the others, the forces are applied in a wrong point.

By the way, I've tested the interaction law with a relative velocity dependent damping and it works. Gives the 'correct' restitution coefficient. I say 'correct' because as indicated in the book of Poschel, in eliminating the attractive forces during the contact, the restitution coefficient is somewhat bigger than the theoretical.

I'm waiting impatiently the tutorial for implementing these new interaction laws because unfortulately I can not assist to the course in Aachen.

Regards,

Question information

Language:
English Edit question
Status:
Solved
For:
ESyS-Particle Edit question
Assignee:
No assignee Edit question
Solved by:
Diego Peinado Martín
Solved:
Last query:
Last reply:
Revision history for this message
Dion Weatherley (d-weatherley) said : 		#1

Hi Diego,

Thanks for the message. I assume that you noticed that this implementation had been excluded from compilation and replaced with another implementation of RotFriction via a #if 1 - #else statement. If you wish to test the implementation, you need to change #if 1 to #if 0.

I had a close look at the code and I'm not sure there is a bug with the direction of forces. In Yucang's implementation, the elastic forces are applied with these commands:
    m_p2->applyForce(-1.0*force,pos);
    m_p1->applyForce(force,pos);

whereas in the other implementation, the relevant commands are:
    m_p2->applyForce(force,pos);
    m_p1->applyForce(-1.0*force,pos);

Notice that the force applied to particle 2 (and particle 1) is in the opposite direction in Yucang's implementation. I think this corrects the sign-error you found.

Congratulations on getting the velocity-dependent damping interaction working! Steffen is still working on the tutorial to add new interactions but it will hopefully be ready for release by the end of next week. Sorry for the delay. Once you have created a new interaction group for this law, would you be willing to contribute that interaction law for inclusion in the next release of ESyS-Particle so other users can benefit from it also? I would greatly appreciate it if you are willing to do so.

Cheers,

Dion.

Revision history for this message
Diego Peinado Martín (diego-peinado) said : 		#2

Hi Dion,
Yes, I saw how the direction of forces is changed from the other implementation. The problem, I think in in the position of the application of forces.

 pos=m_p2->getPos()+(m_p2->getRad()/eq_dist)*D;

So, the correction I think is to change the sign affected by the unit vector. The pos should be:

 pos=m_p2->getPos() **** - ****(m_p2->getRad()/eq_dist)*D;

So the pos vector points to the contact point that is between the particles, and not to a point that is outside the centers of the particles.

Regarding the contribution, of course I will. After implementing the velocity-dependent interaction for the friction particles, and rotating elastic and rotating friction, I'm thinking on implement also a Hertz-Midlin type of interacction.

If I everything goes smooth, and I have time, then I would like to contribute in the solid geometry-particles boundary. I've seen a opensouce library called opencascade, so if possible to be able to import solid or surface geometry from CAD (not only meshes), and to use it for boundary conditions. Well these are plans at long time. First, I have to do the simplest step.

Thanks a lot, best regards,

Revision history for this message
Dion Weatherley (d-weatherley) said : 		#3

Hi Diego,

Thanks for the explanation. I misunderstood your first post. I can see your point about the 'pos' vector. Since this is used to compute frictional forces and moments, this might be the reason Yucang's new implemention caused the Unit Tests to fail. It would appear that this implementation has not been reviewed since it was commented out about 4 years ago. I will take a closer look at a later stage to see if the bug you found might be the reason for the problems.

Thanks for your willingness to contribute to ESyS-Particle. Your idea w.r.t to importing solid geometries and using them for boundary conditions is very interesting. It would be nice to have some ways to import geometries from CAD packages for BCs.

Cheers,

Dion.