yade-dev team mailing list archive

[luc scholtes <lscholtes63@xxxxxxxxx>]

Hi Anton,

Thank you for the link.

After few tests, it seems that CohesivefrictionalPM does not work in its
actual version (is it related to Eigen?...). More particularly, it seems
that the problem comes from the part related to the moment rotation law (see
the concerned code below). I know that it is a pity to have copy/paste this
code from Bruno's work on CohesivefrictionalContactLaw, but I cannot see
where is the problem (I tried to follow all the wm3->Eigen conversions, and
it seems ok).

Does anibody have any suggestion? I am stucked...

  Luc


    /* Moment Rotation Law */
    // NOTE this part could probably be computed in ScGeom to avoid
copy/paste multiplication !!!
    Quaternionr delta( b1->state->ori *
phys->initialOrientation1.conjugate() *phys->initialOrientation2 *
b2->state->ori.conjugate()); //relative orientation
    AngleAxisr aa(angleAxisFromQuat(delta)); // axis of rotation - this is
the Moment direction UNIT vector; angle represents the power of resistant
ELASTIC moment
    if(aa.angle() > Mathr::PI) aa.angle() -= Mathr::TWO_PI; // angle is
between 0 and 2*pi, but should be between -pi and pi

    phys->cumulativeRotation = aa.angle();

    //Find angle*axis. That's all.  But first find angle about contact
normal. Result is scalar. Axis is contact normal.
    Real angle_twist(aa.angle() * aa.axis().dot(geom->normal) ); //rotation
about normal
    Vector3r axis_twist(angle_twist * geom->normal);
    Vector3r moment_twist(axis_twist * phys->kr);

    Vector3r axis_bending(aa.angle()*aa.axis() - axis_twist); //total
rotation minus rotation about normal
    Vector3r moment_bending(axis_bending * phys->kr);
    Vector3r moment = moment_twist + moment_bending;

    Real MomentMax = phys->maxBend*std::fabs(phys->normalForce.norm());
    Real scalarMoment = moment.norm();

    /*Plastic moment */
    if(scalarMoment > MomentMax)
    {
      Real ratio = 0;
      ratio *= MomentMax/scalarMoment; // to fix the moment to its yielding
value
      moment *= ratio;
      moment_twist *=  ratio;
      moment_bending *= ratio;
    }

    phys->moment_twist = moment_twist;
    phys->moment_bending = moment_bending;

    rootBody->forces.addTorque(id1,-moment);
    rootBody->forces.addTorque(id2, moment);

or maybe from the


2010/5/26 Anton Gladky <gladky.anton@xxxxxxxxx>

> Hi, Luc!
>
> Have a look here:
> https://yade-dem.org/wiki/Wm3%E2%86%92Eigen
>
> Eigen is a actively developed mathematical library. Wm3 is, as I know, does
> not develops now. Eigen gives a great functionality through API and
> well-supported.
> ______________________________
>
>