Proper Law2 for Mindlin Physics w/ Poly, Sphere and Facets

Hello,

I am not a polyhedra expert so I may be of little use. I will try to help but first I should give you tips to improve your question.

Please provide an MWE [1].

MWE means you essentially say "based on the polyhedra example script, when I change *BLANK*, the script throws this error." I.e. at what level of complexity does your script work? E.g. does it work when you omit one of geometries? Such as just omitting the walls or spheres? Does it even work with just spheres?

MWE also means someone like me can copy and paste your code into a single script and run, to come up with the same error. It means you probably need to spend some time making a dedicated script that recreates your problem in a more simplistic way. This helps you and us solve the problem simultaneously in the most efficient way possible.

Thanks for following our forum guidelines.

Cheers,

Robert

[1] https://www.yade-dem.org/wiki/Howtoask

Hello Robert,

Thanks for responding so quickly! I am new to open source DEM. Hopefully my company can start submitting valuable contributions soon.

Here is an easy to run example of my problem. First, a working script with only spheres. It is a slightly modified version of the gravity deposition in box tutorial:

from yade import pack,plot,ymport

O.bodies.append(geom.facetBox((.5,.5,.5),(.5,.5,.5),wallMask=31))

sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1,1,1),rMean=.05,rRelFuzz=.5)

id_matSpheres = O.materials.append(ViscElMat(
 young=1.0e8,
 frictionAngle=radians(25),
 density=2160,
 poisson=0.3,
))
sp.toSimulation(material=id_matSpheres)

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()],verletDist=(2.0e-3),label='collider',ompThreads=1),
    InteractionLoop(
        ## handle sphere+sphere and facet+sphere collisions
        [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_MindlinPhys(en=0.3, gamma=0.4, krot=0.3)],
        [Law2_ScGeom_MindlinPhys_Mindlin(includeMoment=True, includeAdhesion=True)],
    ),
    NewtonIntegrator(damping=0,gravity=[0,0,-9.81]),

    VTKRecorder(fileName='Sph-',recorders=['spheres', 'facets', 'velocity', 'mass'],iterPeriod=10000),
    ]

O.dt=.5*PWaveTimeStep()

O.stopAtIter=30001
O.run()

I am running all my scripts in the cloud and using Paraview to view the results so all my scripts so I removed all the code around plots. Now, here is a script that replaces some of the spheres in the previous script with polyhedra. The script exports the new packing in two vtk files which you can load in paraview and see that the code successfully generated a mix of spheres and polyhedra. After the export, the reports a fatal error when it tries to find an acceptable Law2 for this interaction.

from yade import pack,plot,ymport,export,polyhedra_utils

#define materials
id_matSpheres = O.materials.append(ViscElMat(
 young=1.0e8,
 frictionAngle=radians(25),
 density=2160,
 poisson=0.3,
))

matPoly = PolyhedraMat(
 young=1.0e8,
 frictionAngle=radians(25),
 density=2160,
 poisson=0.3,
)
id_matPoly = O.materials.append(matPoly)

# Add box to hold our particles
O.bodies.append(geom.facetBox((.5,.5,.5),(.5,.5,.5),wallMask=31))

# Create cloud of spheres
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1,1,1),rMean=.05,rRelFuzz=.5)

# Replace large spheres by polys
for pos,radius in sp:
    if radius > 0.04:
        t = polyhedra_utils.polyhedra(matPoly,(radius,radius,radius))
        t.state.pos = pos
        O.bodies.append(t)
    else:
        O.bodies.append(utils.sphere(pos,radius))

#export the initial packing
vtkExporter = export.VTKExporter('Mixed.vtk')
vtkExporter.exportSpheres()
vtkExporter.exportPolyhedra()

O.engines=[
ForceResetter(),
 InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()],verletDist=(2.0e-3),label='collider',ompThreads=1),
 InteractionLoop(
  ## handle sphere+sphere and facet+sphere collisions
  [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom(), Ig2_Sphere_Polyhedra_ScGeom(), Ig2_Facet_Polyhedra_PolyhedraGeom(), Ig2_Polyhedra_Polyhedra_PolyhedraGeom()],
  [Ip2_FrictMat_FrictMat_MindlinPhys(en=0.3, gamma=0.4, krot=0.3)],
  [Law2_ScGeom_MindlinPhys_Mindlin(includeMoment=True, includeAdhesion=True), Law2_ScGeom_FrictPhys_CundallStrack()],
 ),
    NewtonIntegrator(damping=0,gravity=[0,0,-9.81]),

    VTKRecorder(fileName='Sph-',recorders=['spheres', 'facets', 'velocity', 'mass'],iterPeriod=10000),
    ]

O.dt=.5*PWaveTimeStep()

O.stopAtIter=30001
O.run()

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

I believe I solved my problem. Here is the proper loop:

O.engines=[
    ForceResetter(),
    InsertionSortCollider([
        Bo1_Polyhedra_Aabb(),
        Bo1_Sphere_Aabb(),
        Bo1_Facet_Aabb()
    ],verletDist=(2.0e-3),label='collider',ompThreads=1),

    InteractionLoop(
        [
            Ig2_Sphere_Polyhedra_ScGeom(),
            Ig2_Sphere_Sphere_ScGeom(),
            Ig2_Polyhedra_Polyhedra_PolyhedraGeom(),
            Ig2_Facet_Polyhedra_PolyhedraGeom(),
            Ig2_Facet_Sphere_ScGeom(),
        ],
        [
            Ip2_FrictMat_FrictMat_MindlinPhys(en=0.3, gamma=0.4, krot=0.3),
            Ip2_FrictMat_PolyhedraMat_FrictPhys(),
            Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys(),
        ],
        [
            Law2_ScGeom_MindlinPhys_Mindlin(includeMoment=True, includeAdhesion=True),
            Law2_ScGeom_FrictPhys_CundallStrack(),
            Law2_PolyhedraGeom_PolyhedraPhys_Volumetric(),
        ],
    ),
    NewtonIntegrator(damping=0,gravity=[0,0,-9.81]),
]