Skip to content

fix condition that would result in missing virtual interactions when PFacets have parallel edges

bchareyre requested to merge anotherPFacetFix into master

BUG REPORT:

The following script (inspired by this [1]) results in two facets crossing each other with no contact before this MR.

The solution is to remove some code which is a bit suspicious (e.g. checking is Vector3.dot(anotherVector3) is equal to 1. Which, I suspect, only works by chance if the vectors are axis-aligned.) and would lead to ignore most interactions to handle just one special case.

[1] https://answers.launchpad.net/yade/+question/694548

See also !628 (merged)

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=(10,0,0), damping=0.5)
  ]

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

O.materials.append( CohFrictMat( young=1.0e10,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.2, alphaKr = 100.0, alphaKtw = 100.0, label=cmatName ) )

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

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

# pfacet 1
gnh1_0 = O.bodies.append( gridNode([0.0, 0.0, 1.02],radius = 0.1,fixed=False,material=cmatName,color=[0, 0, 0]) )
gnh1_2 = O.bodies.append( gridNode([-0.3, 1.0, 0.5],radius = 0.1,fixed=False,material=cmatName,color=[0, 0, 0]) )
gnh1_3 = O.bodies.append( gridNode([-0.3, -1.0, 0.5],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=True,material=cmatName,color=[0, 0, 0]) )
gnh2_1 = O.bodies.append( gridNode([0.21, 1.0, 1.0],radius = 0.1,fixed=True,material=cmatName,color=[0, 0, 0]) )
gnh2_2 = O.bodies.append( gridNode([0.21, -1.0, 1.0],radius = 0.1,fixed=True,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]
O.dt=1e-5
O.saveTmp()

Merge request reports