How can polyhedron use CpmMat 、Ip2_CpmMat_CpmMat_CpmPhys()and Law2_ScGeom_CpmPhys_Cpm() ?

Asked by weijie on 2021-04-26

Dear all,

I want to use polyhedron to simulate coarse aggregate of concrete. After I create interactions for non-overlapping polyhedron -polyhedron and sphere -polyhedron[1], I find that polyhedron can't use CpmMat. The script and error are as follows:

script:
##################
concrete = O.materials.append(CpmMat(
 young = young,
 poisson = poisson,
 frictionAngle = frictionAngle,
 epsCrackOnset = epsCrackOnset,
 relDuctility = relDuctility,
 sigmaT = sigmaT,
))
poly1= polyhedra_utils.polyhedra(material=concrete,size=(1,1,1),seed=5)
#########################
error:
#########################
poly1= polyhedra_utils.polyhedra(material=concrete,size=(1,1,1),seed=5)
  File "/home/wj/myYade/install/lib/x86_64-linux-gnu/yade-Unknown/py/yade/polyhedra_utils.py", line 44, in polyhedra
    b.mat = material
Boost.Python.ArgumentError: Python argument types in
    None.None(Body, int)
did not match C++ signature:
    None(yade::Body {lvalue}, boost::shared_ptr<yade::Material>)
#########################

1)How can I modify the source code so that polyhedron can use CpmMat correctly?
2)How can I modify the source code so that polyhedron can use Ip2_CpmMat_CpmMat_CpmPhys()and Law2_ScGeom_CpmPhys_Cpm()?

For law2, I found a place that may need to be modified[2]. The calculation of reflength, crosssection and refpd for polyhedron - polyhedron as well as polyhedron - sphere is added.I don't know whether this modification is correct or whether there are any other places that need to be modified. The added code is as follows:
###################################
  const int polyIndex=Polyhedra::getClassIndexStatic();
  if(b1index == polyIndex && b2index == polyIndex) {
   const Vector3r& pos1 = Body::byId(I->id1, scene)->state->pos;
   const Vector3r& pos2 = Body::byId(I->id2, scene)->state->pos;
   Real minRad = (geom->refR1 <= 0 ? geom->refR2 : (geom->refR2 <= 0 ? geom->refR1 : math::min(geom->refR1, geom->refR2)));
   Vector3r shift2 = scene->cell->hSize * I->cellDist.cast<Real>();
   phys->refLength = (pos2 - pos1 + shift2).norm();
   phys->crossSection = Mathr::PI * pow(minRad, 2);
   phys->refPD = geom->refR1 + geom->refR2 - phys->refLength;
  }
  if((b1index == sphereIndex && b2index == polyIndex)||(b1index == polyIndex && b2index == sphereIndex)) {
   const Vector3r& pos1 = Body::byId(I->id1, scene)->state->pos;
   const Vector3r& pos2 = Body::byId(I->id2, scene)->state->pos;
   Real minRad = (geom->refR1 <= 0 ? geom->refR2 : (geom->refR2 <= 0 ? geom->refR1 : math::min(geom->refR1, geom->refR2)));
   Vector3r shift2 = scene->cell->hSize * I->cellDist.cast<Real>();
   phys->refLength = (pos2 - pos1 + shift2).norm();
   phys->crossSection = Mathr::PI * pow(minRad, 2);
   phys->refPD = geom->refR1 + geom->refR2 - phys->refLength;
  }
######################################

Thanks in advance.
Jie

[1]https://answers.launchpad.net/yade/+question/696515
[2]https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/ConcretePM.cpp#L319

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
Last query:
Last reply:
Jan Stránský (honzik) said : #1

Hello,

> script: ...

Please try to make the codes MWE [3], W = working. Here:
NameError: name 'young' is not defined
NameError: name 'polyhedra_utils' is not defined

> 1)How can I modify the source code so that polyhedron can use CpmMat correctly?

see above, just modify your script

> 2)How can I modify the source code so that polyhedron can use Ip2_CpmMat_CpmMat_CpmPhys()and Law2_ScGeom_CpmPhys_Cpm()?
> For law2, I found a place that may need to be modified[2]. ...

the modification seems ok, IMO it is the easiest approach.
It should work for "bouncing" = interactions without significant sliding/rotation.
It would behave strangely for "long sliding" or "largely rotating" interactions. The equivalent radii, area and refLength are set at the interaction creation time and kept constant for the interaction lifetime. If the situation changes significantly (sliding, rotating,...), this easy "hack" should be replaced by a more robust approach.

Just curious, are the two if branches (poly-poly and poly-sphere) different from each other and/or different from the sphere-sphere contact?

Cheers
Jan

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

weijie (amandajoe) said : #2

Hi Jan, and thank you again.

>Please try to make the codes MWE [3], W = working. Here:

Here are my codes:
#####
from __future__ import print_function
from yade import polyhedra_utils,qt

concrete = O.materials.append(CpmMat(
 young = 24e9,
 poisson = .2,
 frictionAngle = atan(0.8),
 epsCrackOnset = 1e-4,
 relDuctility =30,
 sigmaT = 3.5e6,
))
poly1= polyhedra_utils.polyhedra(material=concrete,size=(1,1,1),seed=5)
O.bodies.append(poly1)
poly1.state.pos = (1.5,0,0)
O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Polyhedra_Aabb()]),
    InteractionLoop(
        [Ig2_Polyhedra_Polyhedra_ScGeom()],
        [Ip2_CpmMat_CpmMat_CpmPhys()],
        [Law2_ScGeom_CpmPhys_Cpm()],
    ),
    NewtonIntegrator(),
]
#############

> If the situation changes significantly (sliding, rotating,...), this easy "hack" should be replaced by a more robust approach.

I found in law2_ ScGeom_ CpmPhys_ Cpm,for sphere's equivalent radius, area and reflength does not change when situation changes significantly (sliding, rotating,...). Do you have any suggestions or references for how to change these equivalent radii, area and refLength when situation changes significantly?

>Just curious, are the two if branches (poly-poly and poly-sphere) different from each other and/or different from the sphere-sphere contact?

Now it's the same. After that, I'll think if there is a better way.

Best regards
Jie

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

Hello,

>> 1)How can I modify the source code so that polyhedron can use CpmMat correctly?
> see above, just modify your script

sorry, I have included the corrected script, but then probably deleted it..
The problem in your code is that polyhedra_utils.polyhedra expects material to be actual Material, but your "concrete" is int (returned by O.materials.append).
The error says that. A bit cryptically thought..

Change it to
###
concrete = CpmMat(...) # concrete is Material
concreteId = O.materials.append(concrete)# concreteId is int
poly1= polyhedra_utils.polyhedra(material=concrete,...) # material=Material, ok
###
and it should work.
So no need to modify source code for this problem.

> Do you have any suggestions or references for how to change these equivalent radii, area and refLength when situation changes significantly?

No..
I **personally** would start with this simplified version, having in mind the simplification, testing from time to time if it is sufficient or not and deal with the more general solution only if necessary

Cheers
Jan

weijie (amandajoe) said : #4

Hi Jan, and thank you again.

I modified the script as you say,it works.But when I add walls to compress polyhedra, there's another error. The script and errors are as follows:

script:
######################
from __future__ import print_function
from yade import polyhedra_utils,qt
concrete = CpmMat(
 young = 24e9,
 poisson = .2,
 frictionAngle = atan(0.8),
 epsCrackOnset = 1e-4,
 relDuctility =30,
 sigmaT = 3.5e6,
)
concreteId = O.materials.append(concrete)
poly1= polyhedra_utils.polyhedra(material=concrete,size=(1,1,1),seed=5)
O.bodies.append(poly1)
poly1.state.pos = (0,0,0)
w1 = utils.wall((-0.75,0,0), axis=0, sense=0)
w2 = utils.wall((0.75,0,0), axis=0, sense=0)
v=0.1
w1.state.vel = (+v,0,0)
w2.state.vel = (-v,0,0)
O.bodies.append((w1,w2))
O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Wall_Aabb()]),
    InteractionLoop(
        [Ig2_Polyhedra_Polyhedra_ScGeom(),Ig2_Wall_Polyhedra_PolyhedraGeom()],
        [Ip2_CpmMat_CpmMat_CpmPhys()],
        [Law2_ScGeom_CpmPhys_Cpm()],
    ),
    NewtonIntegrator(),
]

O.dt = 1e-6
O.run()
###############

error:

<FATAL ERROR> ThreadRunner:35 void yade::ThreadRunner::run(): Exception occured:
Body #0: Body::material type CpmMat doesn't correspond to Body::state type State (should be CpmState instead).

How to solve it?

Best regards
Jie

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

One approach is this:
- add "basic" Frict... stuff to InteractionLoop
- create the packing, spheres/polyhedrons with CpmMat, walls with FrictMat, set interactionDetectionFactor, aabbEnlargeFactor etc.
- run one iteration to create cohesive bonds
- reset interactionDetectionFactor, aabbEnlargeFactor..
- set material to all particles to FrictMat

This way, cohesive bonds would remain cohesive Cpm interactions (modifying particle material ha no effect on existing interactions) and new contacts would be of FrcitMat type.
This IMO makes physical sense, Cpm gives no extra features for non-cohesive interactions plus you can independently set friction angle, which has different meaning for cohesive and non-cohesive interactions in Cpm.

Cheers
Jan

weijie (amandajoe) said : #6

Hi Jan, and thank you again.

I tried to make a simple version first. Based on the above script(#4), firstly, the polyhedron material was set to CPM, the material of the wall was set to FrictMat, and then I ran a step to set the polyhedron material to frictmat, and an error occurred. The script and error are as follows:

script:
###########
from __future__ import print_function
from yade import polyhedra_utils,qt
concrete = CpmMat(
 young = 24e9,
 poisson = .2,
 frictionAngle = atan(0.8),
 epsCrackOnset = 1e-4,
 relDuctility =30,
 sigmaT = 3.5e6,
)
concreteId = O.materials.append(concrete)
poly1= polyhedra_utils.polyhedra(material=concrete,size=(1,1,1),seed=5)
O.bodies.append(poly1)
poly1.state.pos = (0,0,0)
m=FrictMat(density=1000,young=1e5,poisson=0.5,frictionAngle=radians(20))
O.materials.append(m)
w1 = utils.wall((-0.75,0,0), axis=0, sense=0,material=m)
w2 = utils.wall((0.75,0,0), axis=0, sense=0,material=m)
v=0.1
w1.state.vel = (+v,0,0)
w2.state.vel = (-v,0,0)
O.bodies.append((w1,w2))

O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Wall_Aabb()]),
    InteractionLoop(
        [Ig2_Polyhedra_Polyhedra_ScGeom(),Ig2_Wall_Polyhedra_PolyhedraGeom()],
        [Ip2_CpmMat_CpmMat_CpmPhys(),Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_CpmPhys_Cpm(),Law2_ScGeom_FrictPhys_CundallStrack()],
    ),
    NewtonIntegrator(),
]

O.dt = 0. #run one iteration
poly1.material=m #change material
O.dt = 1e-6
##########################
error:

<FATAL ERROR> InteractionLoop:182 virtual void yade::InteractionLoop::action(): None of given Law2 functors can handle interaction #1+0, types geom:PolyhedraGeom=10 and phys:FrictPhys=3 (LawDispatcher::getFunctor2D returned empty functor)
QObject::~QObject: Timers cannot be stopped from another thread

Best regards
Jie

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

Hello,

sorry for late answer

> and then I ran a step to set the polyhedron material to frictmat
>
> O.dt = 0. #run one iteration
> poly1.material=m #change material
> O.dt = 1e-6

If you insist on automatic interaction creation, you need O.step (or similar) before changing poly1.material.

Concerning the error, you probably need a new Ig2_Wall_Polyhedra_ScGeom..
You can hack it by creating the walls made of polyhedron.

Cheers
Jan

weijie (amandajoe) said : #8

Hi Jan, and thank you again.

>If you insist on automatic interaction creation, you need O.step (or similar) before changing poly1.material.

After I added O. step (), I found the error of material type doesn't respond to body.The script and error see[*].
If I want to create the initial cohesive contact manually and then reset material ,the method is like below?
##
createInteraction(poly1.id,poly2.id)
O.run(1,True)
poly1.material=FrictMat
ss2sc.interactionDetectionFactor=1.

>You can hack it by creating the walls made of polyhedron.
Walls created by utils.wall will still move at the same speed after contacting with other bodies.How to make the wall created by polyhedron move in the same speed after colliding?

Best regards
Jie

[*]script:
###########
from __future__ import print_function
from yade import polyhedra_utils,qt
concrete = CpmMat(
 young = 24e9,
 poisson = .2,
 frictionAngle = atan(0.8),
 epsCrackOnset = 1e-4,
 relDuctility =30,
 sigmaT = 3.5e6,
)
concreteId = O.materials.append(concrete)
poly1= polyhedra_utils.polyhedra(material=concrete,size=(1,1,1),seed=5)
O.bodies.append(poly1)
poly1.state.pos = (0,0,0)
m=FrictMat(density=1000,young=1e5,poisson=0.5,frictionAngle=radians(20))
O.materials.append(m)

######### Create wall
wallThickness = 0.1
length = 1.5
height = 1.5
vC =((-0.5*length,0.5*wallThickness,0.5*height),(0.5*length,0.5*wallThickness,0.5*height),(0.5*length,0.5*wallThickness,-0.5*height),(-0.5*length,0.5*wallThickness,-0.5*height),(-0.5*length,-0.5*wallThickness,0.5*height),(0.5*length,-0.5*wallThickness,0.5*height),(0.5*length,-0.5*wallThickness,-0.5*height),(-0.5*length,-0.5*wallThickness,-0.5*height))
vD=((-0.5*length,0.5*wallThickness,0.5*height),(0.5*length,0.5*wallThickness,0.5*height),(0.5*length,0.5*wallThickness,-0.5*height),(-0.5*length,0.5*wallThickness,-0.5*height),(-0.5*length,-0.5*wallThickness,0.5*height),(0.5*length,-0.5*wallThickness,0.5*height),(0.5*length,-0.5*wallThickness,-0.5*height),(-0.5*length,-0.5*wallThickness,-0.5*height))
bC = polyhedra_utils.polyhedra(material=m,v=vC,color=[0,0.5,1])
bC.state.pos = [0,0.45*length,0]
bD = polyhedra_utils.polyhedra(material=m,v=vD,color=[0,0.5,1])
bD.state.pos = [0,-0.45*length,0]
O.bodies.append((bC,bD))
v=0.1
bC.state.vel = (0,-v,0)
bD.state.vel = (0,+v,0)

O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Polyhedra_Aabb()]),
    InteractionLoop(
        [Ig2_Polyhedra_Polyhedra_ScGeom()],
        [Ip2_CpmMat_CpmMat_CpmPhys(),Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_CpmPhys_Cpm(),Law2_ScGeom_FrictPhys_CundallStrack()],
    ),
    NewtonIntegrator(),
]

O.dt = 0. #run one iteration
O.step()
poly1.material=m #change material
O.dt = 1e-6
##########################
error:
    O.step()
RuntimeError: Body #0: Body::material type CpmMat doesn't correspond to Body::state type State (should be CpmState instead).

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

> Walls created by utils.wall will still move at the same speed after contacting with other bodies.How to make the wall created by polyhedron move in the same speed after colliding?

poly.state.blockedDOFs = "xyzXYZ"

> RuntimeError: Body #0: Body::material type CpmMat doesn't correspond to Body::state type State (should be CpmState instead).

poly1.state = CpmState()

Cheers
Jan

weijie (amandajoe) said : #10

Hi Jan, and thank you again.

>poly1.state = CpmState()
After I add poly1. state = CpmState(), I find poly1 disappeared in the scene.I've tried for a long time and it's still like this.script see[1].

Best regards
Jie

################################
from __future__ import print_function
from yade import polyhedra_utils,qt
concrete = CpmMat(
 young = 24e9,
 poisson = .2,
 frictionAngle = atan(0.8),
 epsCrackOnset = 1e-4,
 relDuctility =30,
 sigmaT = 3.5e6,
)
concreteId = O.materials.append(concrete)
poly1= polyhedra_utils.polyhedra(material=concrete,size=(1,1,1),seed=5)
poly1.state = CpmState()
poly1.state.pos = (0,0,0)
O.bodies.append(poly1)
m=FrictMat(density=1000,young=1e5,poisson=0.5,frictionAngle=radians(20))
O.materials.append(m)

######### Create wall
wallThickness = 0.1
length = 1.5
height = 1.5
vC =((-0.5*length,0.5*wallThickness,0.5*height),(0.5*length,0.5*wallThickness,0.5*height),(0.5*length,0.5*wallThickness,-0.5*height),(-0.5*length,0.5*wallThickness,-0.5*height),(-0.5*length,-0.5*wallThickness,0.5*height),(0.5*length,-0.5*wallThickness,0.5*height),(0.5*length,-0.5*wallThickness,-0.5*height),(-0.5*length,-0.5*wallThickness,-0.5*height))
vD=((-0.5*length,0.5*wallThickness,0.5*height),(0.5*length,0.5*wallThickness,0.5*height),(0.5*length,0.5*wallThickness,-0.5*height),(-0.5*length,0.5*wallThickness,-0.5*height),(-0.5*length,-0.5*wallThickness,0.5*height),(0.5*length,-0.5*wallThickness,0.5*height),(0.5*length,-0.5*wallThickness,-0.5*height),(-0.5*length,-0.5*wallThickness,-0.5*height))
bC = polyhedra_utils.polyhedra(material=m,v=vC,color=[0,0.5,1])
bC.state.pos = [0,0.5*length,0]
bD = polyhedra_utils.polyhedra(material=m,v=vD,color=[0,0.5,1])
bD.state.pos = [0,-0.5*length,0]
O.bodies.append((bC,bD))
v=0.1
bC.state.vel = (0,-v,0)
bD.state.vel = (0,+v,0)
bC.state.blockedDOFs = "xyzXYZ"
bD.state.blockedDOFs = "xyzXYZ"
O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Polyhedra_Aabb()]),
    InteractionLoop(
        [Ig2_Polyhedra_Polyhedra_ScGeom()],
        [Ip2_CpmMat_CpmMat_CpmPhys(),Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_CpmPhys_Cpm(),Law2_ScGeom_FrictPhys_CundallStrack()],
    ),
    NewtonIntegrator(),
]

O.dt = 0. #run one iteration
O.step()
poly1.material=m #change material
O.dt = 1e-6

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

CpmState() is "fresh" e.g. with zero mass (after O.step, pos is NaN because of zero mass and therefore "disappeared"), you have to "clone" the state:
###
oldState = poly1.state
poly1.state = CpmState()
poly1.state.mass = oldState.mass
# same for inertia, ori, ...
###

cheers
Jan

weijie (amandajoe) said : #12

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