"Law2_ScGeom_MindlinPhys_Mindlin ()" engine problem

Asked by Pengfei Tan on 2019-12-20

Hello,

I would like to consider the van Der van der Waals force between spherical particles using the DMT model in the engines "Ip2_FrictMat_FrictMat_MindlinPhys()" and "Law2_ScGeom_MindlinPhys_Mindlin ()". The code works well when only spherical particles exist. When I add the cylindrical particles into the code, the contact between spherical and cylindrical particles causes the following issues:

/build/yade-fDuCoe/yade-2018.02b/lib/multimethods/DynLibDispatcher.hpp:344: ambiguous 2d dispatch (arg1=GridCoGridCoGeom, arg2=MindlinPhys, distance=1), dispatch matrix:
AMBIGUOUS: 1+13 -> Law2_ScGeom_MindlinPhys_Mindlin
AMBIGUOUS: 2+5 -> Law2_ScGeom6D_CohFrictPhys_CohesionMoment
AMBIGUOUS: 6+3 -> Law2_GridCoGridCoGeom_FrictPhys_CundallStrack
AMBIGUOUS: 6+13 -> Law2_GridCoGridCoGeom_FrictPhys_CundallStrack
AMBIGUOUS: 7+5 -> Law2_ScGeom6D_CohFrictPhys_CohesionMoment
AMBIGUOUS: 11+3 -> Law2_ScGridCoGeom_FrictPhys_CundallStrack
AMBIGUOUS: 11+13 -> Law2_ScGridCoGeom_FrictPhys_CundallStrack
terminate called after throwing an instance of 'std::runtime_error'
  what(): Ambiguous dispatch.
Aborted (core dumped)

My code is as follows:

#################################################

from yade.gridpfacet import *
from yade import pack, plot,qt,export

fac=1000
radius=4e-6*fac
l=1e-3*fac
th=5e-6*fac
tsart=10
rad=30e-6*fac

#### Material properties ####
###cylinders###
Df=2590 #Density
Vf=0.21 #Poisson ratio
Ef=7.2e10 #Young's modulus

###powders###
Dp=1020 #Density
Vp=0.393 #Poisson ratio
Ep=1.935e9 #Young's modulus
#Friction coefficient 0.35

###walls###
Dw=8030 #Density
Vw=0.265 #Poisson ratio
Ew=1.90e10 #Young's modulus

O.engines=[
 ForceResetter(),
 InsertionSortCollider([
  Bo1_GridConnection_Aabb(),
  Bo1_PFacet_Aabb(),
  Bo1_Sphere_Aabb(),
 ]),
 InteractionLoop([
  Ig2_GridNode_GridNode_GridNodeGeom6D(),
  Ig2_Sphere_Sphere_ScGeom(),
  #Ig2_Sphere_Sphere_ScGeom6D(),
  Ig2_Sphere_GridConnection_ScGridCoGeom(),
  Ig2_GridConnection_GridConnection_GridCoGridCoGeom(), # cylinder-cylinder interaction
  Ig2_Sphere_PFacet_ScGridCoGeom(), # needed for GridNode-pFacet interaction (why is this not included in Ig2_GridConnection_PFacet_ScGeom???)
  Ig2_GridConnection_PFacet_ScGeom(), # Cylinder-pFcet interaction
 ],
 [
  Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=True),
  #Ip2_FrictMat_FrictMat_FrictPhys(),
  Ip2_FrictMat_FrictMat_MindlinPhys(gamma=40e-3),
 ],
 [
  #Law2_ScGeom_FrictPhys_CundallStrack(), # contact law for sphere-sphere
  Law2_ScGeom_MindlinPhys_Mindlin(includeAdhesion=True),
  Law2_ScGeom6D_CohFrictPhys_CohesionMoment(), # contact law for "internal" cylider forces
  Law2_ScGridCoGeom_FrictPhys_CundallStrack(), # contact law for Cylinder-pFacet interaction
  Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),# contact law for cylinder-cylinder interaction
 ]
 ),
 GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.01,label='ts'),
 NewtonIntegrator(gravity=(0.,0,-9.8),damping=0.2,label='newton'),
]
O.dt=0.1*PWaveTimeStep()

#cylindrical roller: the IDs of cylinder must be [0, 1, 2] which are used for translation and rotation engines.
O.materials.append( CohFrictMat( young=Ef,poisson=Vf,density=Df,frictionAngle=radians(18),normalCohesion=1e20,shearCohesion=1e20,momentRotationLaw=True,label='cMat' ) )
O.materials.append( FrictMat( young=Ef,poisson=Vf,density=Df,frictionAngle=radians(18),label='fMat' ) ) # material for general interactions

O.materials.append( FrictMat( young=Ep,poisson=Vp,density=Dp,frictionAngle=radians(18),label='sphereMat' ) ) # material for general interactions
O.materials.append( CohFrictMat( young=Ep,poisson=Vp,density=Dp,frictionAngle=radians(30),normalCohesion=1e10,shearCohesion=1e10,momentRotationLaw=True,label='sphereMatc' ) )

O.materials.append( CohFrictMat( young=Ew,poisson=Vw,density=Dw,frictionAngle=radians(30),normalCohesion=1e20,shearCohesion=1e20,momentRotationLaw=True,label='cMatw' ) )
O.materials.append( FrictMat( young=Ew,poisson=Vw,density=Dw,frictionAngle=radians(30),label='fMatw' ) ) # material for general interactions

nodesIds=[]
cylIds=[]

cylinder((-l/2,0,1*l),(l/2,0,1*l),radius=radius,nodesIds=nodesIds,cylIds=cylIds,color=[1,0,0],fixed=False,intMaterial='cMat',extMaterial='fMat')

s1=O.bodies.append(sphere([0,0,2*rad],radius=rad,fixed=True,material='sphereMat'))
s2=O.bodies.append(sphere([0,0,0.5*l],radius=rad,fixed=False,material='sphereMat'))

#Bottom surface
color=[255./255.,102./255.,0./255.]
n0=O.bodies.append( gridNode([-l,-l,0],th,wire=False,fixed=True,material='cMatw',color=color) )
n1=O.bodies.append( gridNode([l,-l,0],th,wire=False,fixed=True,material='cMatw',color=color) )
n2=O.bodies.append( gridNode([l,l,0],th,wire=False,fixed=True,material='cMatw',color=color) )
n3=O.bodies.append( gridNode([-l,l,0],th,wire=False,fixed=True,material='cMatw',color=color) )
O.bodies.append( gridConnection(n0,n1,th,color=color,material='fMatw') )
O.bodies.append( gridConnection(n1,n2,th,color=color,material='fMatw') )
O.bodies.append( gridConnection(n2,n0,th,color=color,material='fMatw') )
O.bodies.append( gridConnection(n2,n3,th,color=color,material='fMatw') )
O.bodies.append( gridConnection(n3,n0,th,color=color,material='fMatw') )
O.bodies.append( pfacet(n0,n1,n2,wire=False,material='fMatw',color=color) )
O.bodies.append( pfacet(n0,n2,n3,wire=False,material='fMatw',color=color) )

#### For viewing ####
qt.View()
Gl1_Sphere.stripes=True

#### Allows to reload the simulation ####
O.saveTmp()

##########################################################

Thank you so much for your help!

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2020-01-07
Last query:
2020-01-07
Last reply:
2020-01-07
Jérôme Duriez (jduriez) said : #2

See https://yade-dem.org/doc/user.html#functors-choice.

Is your Christmas list with plenty of Ig2* and Law2* functors thoughtfully designed for your needs ?

Pengfei Tan (tpf516) said : #3

Thank you for providing me the information about how to choose engine functors, Jérôme, but I still have the problem.

I want to consider the adhesive force between spherical particles. So, I used the "Law2_ScGeom_MindlinPhys_Mindlin ()" functor. My model also includes the cylindrical particles, which needs the following Law2 functors:
"Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),"
"Law2_ScGridCoGeom_FrictPhys_CundallStrack(),"
"Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),".

When I used the above four functors, I always encounter the issues mentioned in the original question. Could you help me check if the "Law2_ScGeom_MindlinPhys_Mindlin ()" functor (for considering the adhesive force between spherical particles) can be used together with cylinders?

Merry Christmas in advance!

Pengfei Tan (tpf516) said : #4

I still need help with this problem. Anyone can help me? Thank you!

Jérôme Duriez (jduriez) said : #5

Where did you see (Hertz-)Mindlin contact model included adhesive forces ?

Jan Stránský (honzik) said : #6

Hello,

> /build/yade-fDuCoe/yade-2018.02b/lib/multimethods/DynLibDispatcher.hpp:344: ambiguous 2d dispatch
(arg1=GridCoGridCoGeom, arg2=MindlinPhys, distance=1), dispatch matrix:

the error says, that you have an interaction with GridCoGridCoGeom and MindlinPhys, but no Law2 for this combination.

cheers
Jan

Jan Stránský (honzik) said : #7

Furthermore, I could not reproduced the problem with the code you have provides..
What version of Yade do you use?

Pengfei Tan (tpf516) said : #8

Thanks for your reply, Jérôme.

I saw that the MindlinPhys can consider the adhesive force from the Yade manual (Here, the adhesive force means van Der Waals force). "Law2_ScGeom_MindlinPhys_Mindlin" has the attribute "includeAdhesion" which can include the adhesion force by using the DMT model. And, "Ip2_FrictMat_FrictMat_MindlinPhys" has the attribute "gamma" to set the surface energy parameter which is needed for evaluating the adhesion force in the DMT model.

I also found the formulation which calculates the adhesion force in the "HertzMindlin.cpp" as follows:
Real Adhesion = 4.*Mathr::PI*R*gamma; // calculate adhesion force as predicted by DMT theory.

I read Chiara Modenese's thesis "Numerical study of the mechanical properties of lunar soil by the discrete element method". He used the DMT model to calculate the van Der Waals force.

Pengfei Tan (tpf516) said : #9

Thanks, Jan.

I used the Yade 2018.02b + Ubuntu 18.04.3 LTS. The error would occur when the particles contact each other.

I understand the meaning of this error. In my code, I used " Ip2_FrictMat_FrictMat_MindlinPhys()" functor. So, as two particles whose materials are defined by the "FrictMat" function collide each other, Yade would assign the MindlinPhys theory to deal with the contact. My model has the interaction between the cylinder (gridconnection) and PFacet (gridconnection) which have friction properties. No law2 can handle this contact which includes "GridCoGridCoGeom“” geomery and “”MindlinPhys" thoery.

I would like to know if I can realize the following requirements:

For the contact between spherical particles, “ Law2_ScGeom_MindlinPhys_Mindlin()" would be used. For contacts including the gridconnection, law2 functors such as "Law2_ScGridCoGeom_FrictPhys_CundallStrack()" and "Law2_GridCoGridCoGeom_FrictPhys_CundallStrack()" would be chosen by Yade.

Jan Stránský (honzik) said : #10

> I used the Yade 2018.02b + Ubuntu 18.04.3 LTS.

thanks

> The error would occur when the particles contact each other.

I have run 3M timesteps without any problem.
Please update the script to be MWE also w.r.t time (making the problematic interaction in the first or a few time steps).

> I would like to know if I can realize the following requirements:

probably not "as it is".
If you had fewer Law2s, it could be feasible.
In the case Yade could not find the Law2 exactly matching given IGeom and IPhys, it tries to use one based on igeom and iphys base classes (at least I think so, but I am not sure :-).
In your case, it finds a lot of matches (since ScGeom is base of GridCoGridCoGeom and ScGeom6D and FrictPhys is base of CohFrictPhys and MindlinPhys) and do not know what to do, therefore ending with the ambiguity error.

A solution is to add another Law2, exactly matching GridCoGridCoGeom and MindlinPhys.
Internally, it can just call any other Law2.
This way, the dispatch would not be ambiguous.
However, this solution needs modification of source code..

cheers
Jan

Pengfei Tan (tpf516) said : #11

Thank you for your patient reply, Jan.

>Please update the script to be MWE also w.r.t time (making the problematic interaction in the first or a few time steps).

 Sorry for not providing a suitable code in which the error can occur fast.

Please try the following code:

############################################################

from yade.gridpfacet import *
from yade import pack, plot,qt,export

radius=4e-2
l=10
th=5e-2
rad=30e-2

#### Material properties ####
###cylinders###
Df=2590 #Density
Vf=0.21 #Poisson ratio
Ef=1e8 #Young's modulus

###powders###
Dp=1020
Vp=0.393
Ep=1e8

###walls###
Dw=8030
Vw=0.265
Ew=1e8

O.engines=[
 ForceResetter(),
 InsertionSortCollider([
  Bo1_GridConnection_Aabb(),
  Bo1_PFacet_Aabb(),
  Bo1_Sphere_Aabb(),
 ]),
 InteractionLoop([
  Ig2_GridNode_GridNode_GridNodeGeom6D(),
  Ig2_Sphere_Sphere_ScGeom(),
  #Ig2_Sphere_Sphere_ScGeom6D(),
  Ig2_Sphere_GridConnection_ScGridCoGeom(),
  Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
  Ig2_Sphere_PFacet_ScGridCoGeom(),
  Ig2_GridConnection_PFacet_ScGeom(),
 ],
 [
  Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=True),
  Ip2_FrictMat_FrictMat_MindlinPhys(gamma=40e-3),
  #Ip2_FrictMat_FrictMat_FrictPhys(),
 ],
 [
  #Law2_ScGeom_FrictPhys_CundallStrack(),
  Law2_ScGeom_MindlinPhys_Mindlin(includeAdhesion=True),
  Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),
  Law2_ScGridCoGeom_FrictPhys_CundallStrack(),
  Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),
 ]
 ),
 GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.01,label='ts'),
 NewtonIntegrator(gravity=(0.,0,-9.8),damping=0.5,label='newton'),
]
O.dt=0.1*PWaveTimeStep()

O.materials.append( CohFrictMat( young=Ef,poisson=Vf,density=Df,frictionAngle=radians(18),normalCohesion=1e20,shearCohesion=1e20,momentRotationLaw=True,label='cMat' ) )
O.materials.append( FrictMat( young=Ef,poisson=Vf,density=Df,frictionAngle=radians(18),label='fMat' ) )

O.materials.append( FrictMat( young=Ep,poisson=Vp,density=Dp,frictionAngle=radians(18),label='sphereMat' ) )

O.materials.append( CohFrictMat( young=Ew,poisson=Vw,density=Dw,frictionAngle=radians(30),normalCohesion=1e20,shearCohesion=1e20,momentRotationLaw=True,label='cMatw' ) )
O.materials.append( FrictMat( young=Ew,poisson=Vw,density=Dw,frictionAngle=radians(30),label='fMatw' ) )

nodesIds=[]
cylIds=[]
cylinder((-l/2,0,1*l),(l/2,0,1*l),radius=radius,nodesIds=nodesIds,cylIds=cylIds,color=[1,0,0],fixed=False,intMaterial='cMat',extMaterial='fMat')

s1=O.bodies.append(sphere([0,0,2*rad],radius=rad,fixed=True,material='sphereMat'))
s2=O.bodies.append(sphere([0,0,0.5*l],radius=rad,fixed=False,material='sphereMat'))

#Bottom surface
color=[255./255.,102./255.,0./255.]
n0=O.bodies.append( gridNode([-l,-l,0],th,wire=False,fixed=True,material='cMatw',color=color) )
n1=O.bodies.append( gridNode([l,-l,0],th,wire=False,fixed=True,material='cMatw',color=color) )
n2=O.bodies.append( gridNode([l,l,0],th,wire=False,fixed=True,material='cMatw',color=color) )
n3=O.bodies.append( gridNode([-l,l,0],th,wire=False,fixed=True,material='cMatw',color=color) )
O.bodies.append( gridConnection(n0,n1,th,color=color,material='fMatw') )
O.bodies.append( gridConnection(n1,n2,th,color=color,material='fMatw') )
O.bodies.append( gridConnection(n2,n0,th,color=color,material='fMatw') )
O.bodies.append( gridConnection(n2,n3,th,color=color,material='fMatw') )
O.bodies.append( gridConnection(n3,n0,th,color=color,material='fMatw') )
O.bodies.append( pfacet(n0,n1,n2,wire=False,material='fMatw',color=color) )
O.bodies.append( pfacet(n0,n2,n3,wire=False,material='fMatw',color=color) )

#### For viewing ####
qt.View()
Gl1_Sphere.stripes=True

#### Allows to reload the simulation ####
O.saveTmp()

############################################################

Ps: In the new code, I just reduce Young's modulus and increase geometrical dimensions so that the time step can be increased.

>A solution is to add another Law2, exactly matching GridCoGridCoGeom and MindlinPhys.

In my case, there are PFacet, cylindrical and spherical particles. Actually, I only want to use the MindlinPhys (which includes the adhesion force) for the contact between spherical particles, which is realized by "Law2_ScGeom_MindlinPhys_Mindlin".

Would it be possible to modify the source code to let Yade choose the "Law2_ScGeom_MindlinPhys_Mindlin" for the contact between spherical particles and other Law2s (such as "Law2_ScGridCoGeom_FrictPhys_CundallStrack" and "Law2_GridCoGridCoGeom_FrictPhys_CundallStrack" ) for the contact associated with PFacet and cylindrical particles?

Best Jan Stránský (honzik) said : #12

Thanks for updating the code

> Would it be possible to modify the source code ...

of course yes :-)
the right questions are who would do it, how long would it take, how difficult it would be ...

Currently, the Law2 is chosen merely by IGeom and IPhys, but in principle, it could take / consider other parameters, like Shape, groupMask or whatever.

E.g. you could have some new Law2_ScGeom_FrictPhys_furtherDispacth, which would decide the specific behavior according to Shape (Law2::go has Interaction argument, which has access to bodies and their shapes).

For **this specific case**, the easiest way, as I pointed, is to implement something like
Law2_GridCoGridCoGeom_MindlinPhys_asCundallStrack
such that for grid connections it internally calls another existing Law2 (possibly just Law2_ScGeom_FrictPhys_CundallStrack) and for sphere-sphere you let the default Mindlin behavior..

cheers
Jan

Pengfei Tan (tpf516) said : #13

Sincerely appreciate your help, Jan.

I will follow your suggestions to modify the source code.

Best,

Pengfei

Jan Stránský (honzik) said : #14

let us know if you need some suggestion :-)
cheers
Jan

Pengfei Tan (tpf516) said : #15

Thanks Jan Stránský, that solved my question.