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 _______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : yade-dev@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp