Unexpected PFacet Collision Physics

Reported on launchpad by Jeffrey Knowles. Apparently an erroneous assignement of nodal forces.

https://answers.launchpad.net/yade/+question/694548

import sys
p = os.getcwd() # current working directory
sys.path.append(p)
from yade import pack, geom, PyRunner, qt
from yade.gridpfacet import *

############################################# engines ###########################################
O.engines=[

 ForceResetter(),
 InsertionSortCollider(

  [
  Bo1_PFacet_Aabb(),
  Bo1_GridConnection_Aabb()
  ]),

 InteractionLoop(

  [
  Ig2_GridConnection_PFacet_ScGeom(),
  Ig2_Sphere_PFacet_ScGridCoGeom(),
  Ig2_GridNode_GridNode_GridNodeGeom6D(),
  Ig2_Sphere_GridConnection_ScGridCoGeom(),
  Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
  Ig2_PFacet_PFacet_ScGeom()
  ],

  [
  Ip2_FrictMat_FrictMat_FrictPhys(),
  Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False),
  ],

  [
  Law2_ScGeom_FrictPhys_CundallStrack(),
  Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm = True),
  Law2_ScGridCoGeom_FrictPhys_CundallStrack(),
  Law2_ScGridCoGeom_CohFrictPhys_CundallStrack(),
  Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),
  ]),

  NewtonIntegrator(gravity=(0,0,0), damping=0.0)
  ]

O.dt = 1.0e-9 # times step

############################################# materials ###########################################
cmatName = 'cmat'
fmatName = 'fmat'

O.materials.append( CohFrictMat( young=1.0e14,poisson=0.3,density=1000.0, normalCohesion=1.0e12 ,shearCohesion=1.0e12, fragile = False, momentRotationLaw=True, isCohesive = True, etaRoll=1000.0, etaTwist = 1000.0, frictionAngle=0.0, alphaKr = 100.0, alphaKtw = 100.0, label=cmatName ) )

O.materials.append( FrictMat( young=1.0e14,poisson=0.3,density=1000.0,frictionAngle=0.0,label=fmatName ) )

############################################## pfacets #############################################

# pfacet 1
gnh1_0 = O.bodies.append( gridNode([0.0, 0.0, 0.3],radius = 0.1,fixed=False,material=cmatName,color=[0, 0, 0]) )
gnh1_2 = O.bodies.append( gridNode([-1.0, 1.0, 0.3],radius = 0.1,fixed=False,material=cmatName,color=[0, 0, 0]) )
gnh1_3 = O.bodies.append( gridNode([-1.0, -1.0, 0.3],radius = 0.1,fixed=False,material=cmatName,color=[0, 0, 0]) )

O.bodies.append( gridConnection(gnh1_2,gnh1_3,radius = 0.1,color=[0, 0, 0],material=fmatName, wire=False) )
O.bodies.append( gridConnection(gnh1_0,gnh1_2,radius = 0.1,color=[0, 0, 0],material=fmatName, wire=False) )
O.bodies.append( gridConnection(gnh1_0,gnh1_3,radius = 0.1,color=[0, 0, 0],material=fmatName, wire=False) )

O.bodies.append( pfacet(gnh1_0, gnh1_2, gnh1_3,material=fmatName,color=[1, 0, 0]) )

# assign initial velocity to pfacet 1
vi = Vector3(10.0, 0.0, 0.0)
for b in O.bodies:
 b.state.vel = vi

# pfacet 2
gnh2_0 = O.bodies.append( gridNode([0.21, 0.0, 0.0],radius = 0.1,fixed=False,material=cmatName,color=[0, 0, 0]) )
gnh2_1 = O.bodies.append( gridNode([0.21, 1.0, 1.0],radius = 0.1,fixed=False,material=cmatName,color=[0, 0, 0]) )
gnh2_2 = O.bodies.append( gridNode([0.21, -1.0, 1.0],radius = 0.1,fixed=False,material=cmatName,color=[0, 0, 0]) )

O.bodies.append( gridConnection(gnh2_0,gnh2_1,radius = 0.1,color=[0, 0, 0],material=fmatName, wire=False) )
O.bodies.append( gridConnection(gnh2_0,gnh2_2,radius = 0.1,color=[0, 0, 0],material=fmatName, wire=False) )
O.bodies.append( gridConnection(gnh2_1,gnh2_2,radius = 0.1,color=[0, 0, 0],material=fmatName, wire=False) )

O.bodies.append( pfacet(gnh2_0, gnh2_1, gnh2_2,material=fmatName,color=[1, 0, 0]) )

#################################### appearance and view ##########################################
qtr = qt.Renderer()
qt.Controller()
v = qt.View()

qtr.bgColor = [1,1,1]