yade-dev team mailing list archive

[Gaël Lorieul <gael.lorieul@xxxxxxxxxxxxxx>]

Hello all,

I think I fixed issue 694548.

Contrary to my initial belief, it was not related to what I had
observed for gridConnection objects (this will be the subject of a
different thread). Instead, it was a bug regarding the calculation of
weights. I'll show that in the case of the MWE posted by Jeffrey
Knowles (issue 694548's original poster).

The motion of PFacet is dictated by that of the 3 gridNode objects it
connects. Each gridNode recieves a different force value calculated
from the force of the interaction. The function
`Law2_ScGridCoGeom_FrictPhys_CundallStrack::go`
(`pkg/common/Grid.cpp:537`) does that using weights:
```
scene->forces.addForce(geom->id3, geom->weight[0] * -force);
scene->forces.addForce(geom->id4, geom->weight[1] * -force);
scene->forces.addForce(geom->id5, geom->weight[2] * -force);
```

Those weights are calculated in `pkg/common/PFacet.cpp`
(`Ig2_Sphere_PFacet_ScGridCoGeom::projection()`) as

         ‖v₁‖²(v₀·v₂) - (v₁·v₂)(v₀·v₁)
    p₁ = —————————————————————————————        (1a)
             ‖v₀‖²‖v₁‖² - (v₀·v₁)²
         ‖v₀‖²(v₁·v₂) - (v₀·v₂)(v₀·v₁)
    p₂ = —————————————————————————————        (1b)
             ‖v₀‖²‖v₁‖² - (v₀·v₁)²
    p₃ = 1 - p₂ - p₃                          (1c)

But actually, it should be:
         ‖v₁‖²(v₀·v₂) - (v₁·v₂)(v₀·v₁)
    p₂ = —————————————————————————————        (2a)
             ‖v₀‖²‖v₁‖² - (v₀·v₁)²
         ‖v₀‖²(v₁·v₂) - (v₀·v₂)(v₀·v₁)
    p₃ = —————————————————————————————        (2b)
             ‖v₀‖²‖v₁‖² - (v₀·v₁)²
    p₁ = 1 - p₂ - p₃                          (2c)

Hence the bug (MWE "passes" after making the permutation).

In the end, it is only a permutation of p₁,p₂ and p₃. However, let me
stress that the fix I propose is not just a "hack" (me randomly
interchanging p1,p2,p3 until I get the result I expected). Instead, I
did the calculations from scratch again, and the equations (2) are the
result of that. The scans of my calculations are available here [1]. If
need be, I can make a LaTeX transcript of those calculations.

However, it uneases me that the consequences of this bug can be seen
fairly easily, and most certainly the original author of the code would
have noticed it straight away. Especially since the rest of the
function seems well written. Hence perhaps things are more complicated
than that? Could I be breaking something else? To make sure of that, I
attempted to compile the first version of the code to feature this
calculation (commit dfec202), but unfortunately without success 🙁️
(and I don't really want to spend more than the 2.5 days I already
spent just to try to compile (old) stuff). Note that commit dfec202 is
the commit that introduced PFacet objects in Yade.

Final point: although the changes to be made are small, I would really
like to submit it personnaly to the git repository so that I learn how
it is done. I would need a bit of guidance for that though 😉️. I also
think I might need the confirmation that the fix I propose indeed is
relevant and safe to be implemented.

Gaël

[1] http://dl.free.fr/mbGj9Moc0