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

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

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

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

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

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

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

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).

> 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

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

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

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