random dense pack with clumps

Hello,

Yes, there is [1][2].

Cheers,

Robert

[1]https://yade-dem.org/doc/user.html#clumping-particles-together
[2]https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.BodyContainer.replaceByClumps

Thank you so much, it actually help me a lot geometrically but when running the model of clumps using the function " replaceclumps " for my model of spheres with jcpfm material type , the program geometrically makes the spheres can be replaced by clumps but once i run it, this error appears in the terminal.

In [1]: FATAL /build/yade-fDuCoe/yade-2018.02b/core/ThreadRunner.cpp:30 run: Exception occured:
Body #573: Body::material type JCFpmMat doesn't correspond to Body::state type State (should be JCFpmState instead).

Please, needing a hand with this, i dont understand what are the characteristics of "state type" and how to solve that error
i tried changing " JCFpmMat" to "JCFpmState" in all the code but it continues not working.

Just in case ,you need it, im sharing the code with you. Thank you in advance

from yade import pack, qt, plot
from numpy import arange
import numpy as np
import random

#DATOS DE ENTRADA
r=0.07
H=1 #altura del muro
B=.45 #base mayor
b=.25 #base menor
L=4 #longitud del muro
#EN SERVICIO
EmpujeTop=2400#Pa
EmpujeBottom=9100#Pa
damp=0.7
aceloutofplane=0
ii=0

# DEFINICION DEL MATERIAL ******************************************************
mat1 = JCFpmMat()
mat1.cohesion =0 #Pa32.37e3
mat1.density = 2821 #kg/m3
mat1.frictionAngle = radians(58.74) #rad
mat1.poisson = 0.18
mat1.young = 37.04e9 #Pa

# DEFINICION DE LA JUNTA ******************************************************
mat1.jointCohesion = 262000 #Pa
mat1.jointFrictionAngle = radians(40)
mat1.jointNormalStiffness= 1.75e8 #Pa/m
mat1.jointShearStiffness=0.7e8 #Pa/m
#mat1.jointDilationAngle = 0
#mat1.jointTensileStrength = 0
O.materials.append(mat1)

# GEOMETRIA ********************************************************************
a=(B-b)*.5
Vol1=pack.inParallelepiped(o=(0,0,0), a=(L,0,0), b=(0,B,0), c=(0,0,H))
Vol2=pack.inParallelepiped(o=(L,B+a,0), a=(0,B+a,0), b=(0,a+b,H), c=(L,a+b,H))
Vol3=pack.inParallelepiped(o=(L,-a,0), a=(0,-a,0), b=(0,a,H), c=(L,a,H))
Vol4=Vol1-Vol2-Vol3-Vol3
Esferas=pack.randomDensePack(Vol4,radius=r,rRelFuzz=0)
O.bodies.append(Esferas)

relRadList1 = [.04,.03]
relPosList1 = [[0.035,0,0],[.083,0,0]]
#peanut:
relRadList2 = [.05,.05]
relPosList2 = [[0,0,0],[.073,0,0]]

relRadList4 = [0.025,0.025,0.025,0.025]
relPosList4 = [[0,0,0],[.043,0,0],[0.02,0.045,0],[0.02,0,045]]
#stick:
relRadList3 = [0.05,0.025,0.05,0.025]
relPosList3 = [[.04,0,0],[.123,0,0.022],[0.14,0,0],[0.123,0,-0.022]]
templates= []
templates.append(clumpTemplate(relRadii=relRadList1,relPositions=relPosList1))
templates.append(clumpTemplate(relRadii=relRadList2,relPositions=relPosList2))
templates.append(clumpTemplate(relRadii=relRadList3,relPositions=relPosList3))
templates.append(clumpTemplate(relRadii=relRadList4,relPositions=relPosList4))
O.bodies.replaceByClumps(templates,[.5,.3,.1,.1])

# APLICACION DE FUERZAS ********************************************************
listSphereA=[]
listSphereB=[]
listSphereC=[]
listSphereD=[]
for x1 in O.bodies:
 if isinstance(x1.shape,Sphere) and x1.state.pos[1]<3*r and x1.state.pos[2]<=H and x1.state.pos[2]>3*H/4.000:
  listSphereA.append(x1.id)
for x2 in O.bodies:
 if isinstance(x2.shape,Sphere) and x2.state.pos[1]<3*r and x2.state.pos[2]<=3*H/4.000 and x2.state.pos[2]>H/2.000:
  listSphereB.append(x2.id)
for x3 in O.bodies:
 if isinstance(x3.shape,Sphere) and x3.state.pos[1]<3.5*r and x3.state.pos[2]<=H/2.000 and x3.state.pos[2]>H/4.000:
  listSphereC.append(x3.id)
for x4 in O.bodies:
 if isinstance(x4.shape,Sphere) and x4.state.pos[1]<3*r and x4.state.pos[2]<=H/4.000 and x4.state.pos[2]>0.000:
  listSphereD.append(x4.id)
ha=1*H/8.000
hb=3*H/8.000
hc=5*H/8.000
hd=7*H/8.000
m=EmpujeBottom-EmpujeTop
Ea=m*ha+EmpujeTop
Eb=m*hb+EmpujeTop
Ec=m*hc+EmpujeTop
Ed=m*hd+EmpujeTop

Fa=Ea*(H/4.000)*L/len(listSphereA)
Fb=Eb*(H/4.000)*L/len(listSphereB)
Fc=Ec*(H/4.000)*L/len(listSphereC)
Fd=Ed*(H/4.000)*L/len(listSphereD)

def aplicarFuerzaA():
 for i in listSphereA:
  O.forces.setPermF(i,(0,Fa,0))
def aplicarFuerzaB():
 for i in listSphereB:
  O.forces.setPermF(i,(0,Fb,0))
def aplicarFuerzaC():
 for i in listSphereC:
  O.forces.setPermF(i,(0,Fc,0))
def aplicarFuerzaD():
 for i in listSphereD:
  O.forces.setPermF(i,(0,Fd,0.98))

# CONDICIONES DE BORDE *********************************************************

Borde1=utils.wall(0.43,axis=0, sense=0)
Borde2=utils.wall(L-0.43,axis=0, sense=0)
Borde3=utils.wall(0, axis=1, sense=0)
Borde4=utils.wall(0.05, axis=2, sense=0)

O.bodies.append(Borde1)
O.bodies.append(Borde2)
O.bodies.append(Borde3)
O.bodies.append(Borde4)

listSphereBase=[]
#listSphereBorde1=[]
#listSphereBorde2=[]
#restringir movimiento en la base
#for y1 in O.bodies:
# if isinstance(y1.shape,Sphere) and y1.state.pos[2]<=2.5*r:
# listSphereBase.append(y1.id)
#for y2 in O.bodies:
# if isinstance(y2.shape,Sphere) and y2.state.pos[0]<=2.5*r:
# listSphereBorde1.append(y2.id)
#for y3 in O.bodies:
# if isinstance(y3.shape,Sphere) and y3.state.pos[0]>L-2.5*r:
# listSphereBorde2.append(y3.id)
#for i in listSphereBase:
# O.bodies[i].state.blockedDOFs='xyz'
#for i in listSphereBorde1:
# O.bodies[i].state.blockedDOFs='xyzXYZ'
#for i in listSphereBorde2:
# O.bodies[i].state.blockedDOFs='xyzXYZ'

# DEFINICION DE FUNCIONES ******************************************************
listSphereTodo=[]
for m in O.bodies:
 if isinstance(m.shape,Sphere) and m.state.pos[2]<2*H:
  listSphereBase.append(m.id)

for i in listSphereTodo:
 O.bodies[i].state.blockedDOFs='XYZ'

def checkUnbalanced():
 print ('iter %d, unbalanced forces = %f'%(O.iter, utils.unbalancedForce()))
 print 'Fzas en c/ esfera (N)', Fa, Fb, Fc, Fd
# if unbalancedForce()<.00001:
 print ( 'Desplazamiento superior',O.bodies[w].state.pos[1],O.bodies[w].state.pos[2])

listTopBodies=[]
for z1 in O.bodies:
 if isinstance(z1.shape,Sphere) and z1.state.pos[2]<=H and z1.state.pos[2]>H-3.000*r and z1.state.pos[0]>(L/2.000)-1.1*r and z1.state.pos[0]<(L/2.000)+1.1*r and z1.state.pos[1]>B/2.000:
  listTopBodies.append(z1.id)
w=max(listTopBodies)

def addPlotData():
 plot.addData(Ypos=O.bodies[w].state.pos[1],Zpos=O.bodies[w].state.pos[2],iteracion=O.iter,unbal=unbalancedForce())
plot.plots={'Ypos':('Zpos')}

# ENGINE ***********************************************************************
O.engines=[
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Wall_Aabb()]),
  InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
   [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()],
   [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM()]
  ),
  NewtonIntegrator(gravity=(0,aceloutofplane,-9.81),damping=damp),
  PyRunner(firstIterRun=10000,command='aplicarFuerzaA()'),
  PyRunner(firstIterRun=10000,command='aplicarFuerzaB()'),
  PyRunner(firstIterRun=10000,command='aplicarFuerzaC()'),
  PyRunner(firstIterRun=10000,command='aplicarFuerzaD()'),
  PyRunner(iterPeriod=10000,command='checkUnbalanced()'),
  PyRunner(iterPeriod=1,command='addPlotData()'),
  VTKRecorder(iterPeriod=100,recorders=['spheres','jcfpm','facets','colors'],fileName='/tmp/p1-')
]

# DETALLES FINALES *************************************************************
O.dt=0.5*PWaveTimeStep() # establece el delta de tiempo como una fraccion del tiempo critico
plot.plot(subPlots=True) # llama a los sub plots
plot.live=True # ploteo en tiempo real
qt.Controller() # abre la ventana del controlador
V=qt.View() # abre la ventada de la vista
R=yade.qt.Renderer() # llama a la renderizacion
R.bgColor=(1.,1.,1.) # definel color blan
V.screenSize = (900,900) # tamano de pantalla

Thank you so much, it actually help me a lot geometrically but when running the model of clumps using the function " replaceclumps " for my model of spheres with jcpfm material type , the program geometrically makes the spheres can be replaced by clumps but once i run it, this error appears in the terminal.

In [1]: FATAL /build/yade-fDuCoe/yade-2018.02b/core/ThreadRunner.cpp:30 run: Exception occured:
Body #573: Body::material type JCFpmMat doesn't correspond to Body::state type State (should be JCFpmState instead).

Please, needing a hand with this, i dont understand what are the characteristics of "state type" and how to solve that error
i tried changing " JCFpmMat" to "JCFpmState" in all the code but it continues not working.

Just in case ,you need it, im sharing the code with you. Thank you in advance

from yade import pack, qt, plot
from numpy import arange
import numpy as np
import random

#DATOS DE ENTRADA
r=0.07
H=1 #altura del muro
B=.45 #base mayor
b=.25 #base menor
L=4 #longitud del muro
#EN SERVICIO
EmpujeTop=2400#Pa
EmpujeBottom=9100#Pa
damp=0.7
aceloutofplane=0
ii=0

# DEFINICION DEL MATERIAL ******************************************************
mat1 = JCFpmMat()
mat1.cohesion =0 #Pa32.37e3
mat1.density = 2821 #kg/m3
mat1.frictionAngle = radians(58.74) #rad
mat1.poisson = 0.18
mat1.young = 37.04e9 #Pa

# DEFINICION DE LA JUNTA ******************************************************
mat1.jointCohesion = 262000 #Pa
mat1.jointFrictionAngle = radians(40)
mat1.jointNormalStiffness= 1.75e8 #Pa/m
mat1.jointShearStiffness=0.7e8 #Pa/m
#mat1.jointDilationAngle = 0
#mat1.jointTensileStrength = 0
O.materials.append(mat1)

# GEOMETRIA ********************************************************************
a=(B-b)*.5
Vol1=pack.inParallelepiped(o=(0,0,0), a=(L,0,0), b=(0,B,0), c=(0,0,H))
Vol2=pack.inParallelepiped(o=(L,B+a,0), a=(0,B+a,0), b=(0,a+b,H), c=(L,a+b,H))
Vol3=pack.inParallelepiped(o=(L,-a,0), a=(0,-a,0), b=(0,a,H), c=(L,a,H))
Vol4=Vol1-Vol2-Vol3-Vol3
Esferas=pack.randomDensePack(Vol4,radius=r,rRelFuzz=0)
O.bodies.append(Esferas)

relRadList1 = [.04,.03]
relPosList1 = [[0.035,0,0],[.083,0,0]]
#peanut:
relRadList2 = [.05,.05]
relPosList2 = [[0,0,0],[.073,0,0]]

relRadList4 = [0.025,0.025,0.025,0.025]
relPosList4 = [[0,0,0],[.043,0,0],[0.02,0.045,0],[0.02,0,045]]
#stick:
relRadList3 = [0.05,0.025,0.05,0.025]
relPosList3 = [[.04,0,0],[.123,0,0.022],[0.14,0,0],[0.123,0,-0.022]]
templates= []
templates.append(clumpTemplate(relRadii=relRadList1,relPositions=relPosList1))
templates.append(clumpTemplate(relRadii=relRadList2,relPositions=relPosList2))
templates.append(clumpTemplate(relRadii=relRadList3,relPositions=relPosList3))
templates.append(clumpTemplate(relRadii=relRadList4,relPositions=relPosList4))
O.bodies.replaceByClumps(templates,[.5,.3,.1,.1])

# APLICACION DE FUERZAS ********************************************************
listSphereA=[]
listSphereB=[]
listSphereC=[]
listSphereD=[]
for x1 in O.bodies:
 if isinstance(x1.shape,Sphere) and x1.state.pos[1]<3*r and x1.state.pos[2]<=H and x1.state.pos[2]>3*H/4.000:
  listSphereA.append(x1.id)
for x2 in O.bodies:
 if isinstance(x2.shape,Sphere) and x2.state.pos[1]<3*r and x2.state.pos[2]<=3*H/4.000 and x2.state.pos[2]>H/2.000:
  listSphereB.append(x2.id)
for x3 in O.bodies:
 if isinstance(x3.shape,Sphere) and x3.state.pos[1]<3.5*r and x3.state.pos[2]<=H/2.000 and x3.state.pos[2]>H/4.000:
  listSphereC.append(x3.id)
for x4 in O.bodies:
 if isinstance(x4.shape,Sphere) and x4.state.pos[1]<3*r and x4.state.pos[2]<=H/4.000 and x4.state.pos[2]>0.000:
  listSphereD.append(x4.id)
ha=1*H/8.000
hb=3*H/8.000
hc=5*H/8.000
hd=7*H/8.000
m=EmpujeBottom-EmpujeTop
Ea=m*ha+EmpujeTop
Eb=m*hb+EmpujeTop
Ec=m*hc+EmpujeTop
Ed=m*hd+EmpujeTop

Fa=Ea*(H/4.000)*L/len(listSphereA)
Fb=Eb*(H/4.000)*L/len(listSphereB)
Fc=Ec*(H/4.000)*L/len(listSphereC)
Fd=Ed*(H/4.000)*L/len(listSphereD)

def aplicarFuerzaA():
 for i in listSphereA:
  O.forces.setPermF(i,(0,Fa,0))
def aplicarFuerzaB():
 for i in listSphereB:
  O.forces.setPermF(i,(0,Fb,0))
def aplicarFuerzaC():
 for i in listSphereC:
  O.forces.setPermF(i,(0,Fc,0))
def aplicarFuerzaD():
 for i in listSphereD:
  O.forces.setPermF(i,(0,Fd,0.98))

# CONDICIONES DE BORDE *********************************************************

Borde1=utils.wall(0.43,axis=0, sense=0)
Borde2=utils.wall(L-0.43,axis=0, sense=0)
Borde3=utils.wall(0, axis=1, sense=0)
Borde4=utils.wall(0.05, axis=2, sense=0)

O.bodies.append(Borde1)
O.bodies.append(Borde2)
O.bodies.append(Borde3)
O.bodies.append(Borde4)

listSphereBase=[]
#listSphereBorde1=[]
#listSphereBorde2=[]
#restringir movimiento en la base
#for y1 in O.bodies:
# if isinstance(y1.shape,Sphere) and y1.state.pos[2]<=2.5*r:
# listSphereBase.append(y1.id)
#for y2 in O.bodies:
# if isinstance(y2.shape,Sphere) and y2.state.pos[0]<=2.5*r:
# listSphereBorde1.append(y2.id)
#for y3 in O.bodies:
# if isinstance(y3.shape,Sphere) and y3.state.pos[0]>L-2.5*r:
# listSphereBorde2.append(y3.id)
#for i in listSphereBase:
# O.bodies[i].state.blockedDOFs='xyz'
#for i in listSphereBorde1:
# O.bodies[i].state.blockedDOFs='xyzXYZ'
#for i in listSphereBorde2:
# O.bodies[i].state.blockedDOFs='xyzXYZ'

# DEFINICION DE FUNCIONES ******************************************************
listSphereTodo=[]
for m in O.bodies:
 if isinstance(m.shape,Sphere) and m.state.pos[2]<2*H:
  listSphereBase.append(m.id)

for i in listSphereTodo:
 O.bodies[i].state.blockedDOFs='XYZ'

def checkUnbalanced():
 print ('iter %d, unbalanced forces = %f'%(O.iter, utils.unbalancedForce()))
 print 'Fzas en c/ esfera (N)', Fa, Fb, Fc, Fd
# if unbalancedForce()<.00001:
 print ( 'Desplazamiento superior',O.bodies[w].state.pos[1],O.bodies[w].state.pos[2])

listTopBodies=[]
for z1 in O.bodies:
 if isinstance(z1.shape,Sphere) and z1.state.pos[2]<=H and z1.state.pos[2]>H-3.000*r and z1.state.pos[0]>(L/2.000)-1.1*r and z1.state.pos[0]<(L/2.000)+1.1*r and z1.state.pos[1]>B/2.000:
  listTopBodies.append(z1.id)
w=max(listTopBodies)

def addPlotData():
 plot.addData(Ypos=O.bodies[w].state.pos[1],Zpos=O.bodies[w].state.pos[2],iteracion=O.iter,unbal=unbalancedForce())
plot.plots={'Ypos':('Zpos')}

# ENGINE ***********************************************************************
O.engines=[
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Wall_Aabb()]),
  InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
   [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()],
   [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM()]
  ),
  NewtonIntegrator(gravity=(0,aceloutofplane,-9.81),damping=damp),
  PyRunner(firstIterRun=10000,command='aplicarFuerzaA()'),
  PyRunner(firstIterRun=10000,command='aplicarFuerzaB()'),
  PyRunner(firstIterRun=10000,command='aplicarFuerzaC()'),
  PyRunner(firstIterRun=10000,command='aplicarFuerzaD()'),
  PyRunner(iterPeriod=10000,command='checkUnbalanced()'),
  PyRunner(iterPeriod=1,command='addPlotData()'),
  VTKRecorder(iterPeriod=100,recorders=['spheres','jcfpm','facets','colors'],fileName='/tmp/p1-')
]

# DETALLES FINALES *************************************************************
O.dt=0.5*PWaveTimeStep() # establece el delta de tiempo como una fraccion del tiempo critico
plot.plot(subPlots=True) # llama a los sub plots
plot.live=True # ploteo en tiempo real
qt.Controller() # abre la ventana del controlador
V=qt.View() # abre la ventada de la vista
R=yade.qt.Renderer() # llama a la renderizacion
R.bgColor=(1.,1.,1.) # definel color blan
V.screenSize = (900,900) # tamano de pantalla

Hello,

Will you first please reduce the size of your script to a “minimal” working example? Thank you,

Robert

from yade import pack, qt, plot
from numpy import arange
import numpy as np
import random

#INPUT
r=0.07
H=1 #height
B=.45 #base
b=.25 #base
L=4 #lenght
#EN SERVICIO
EmpujeTop=2400#Pa
EmpujeBottom=9100#Pa
damp=0.7
aceloutofplane=0
ii=0

# DEFINICION DEL MATERIAL ******************************************************
mat1 = JCFpmMat()
mat1.cohesion =0 #Pa32.37e3
mat1.density = 2821 #kg/m3
mat1.frictionAngle = radians(58.74) #rad
mat1.poisson = 0.18
mat1.young = 37.04e9 #Pa

# JOINT ******************************************************
mat1.jointCohesion = 262000 #Pa
mat1.jointFrictionAngle = radians(40)
mat1.jointNormalStiffness= 1.75e8 #Pa/m
mat1.jointShearStiffness=0.7e8 #Pa/m
O.materials.append(mat1)

# GEOMETRY VOLUME ********************************************************************
a=(B-b)*.5
Vol1=pack.inParallelepiped(o=(0,0,0), a=(L,0,0), b=(0,B,0), c=(0,0,H))
Vol2=pack.inParallelepiped(o=(L,B+a,0), a=(0,B+a,0), b=(0,a+b,H), c=(L,a+b,H))
Vol3=pack.inParallelepiped(o=(L,-a,0), a=(0,-a,0), b=(0,a,H), c=(L,a,H))
Vol4=Vol1-Vol2-Vol3-Vol3
Esferas=pack.randomDensePack(Vol4,radius=r,rRelFuzz=0)
O.bodies.append(Esferas)

#RREPLACE CLUMPS********************************************************
relRadList1 = [.04,.03]
relPosList1 = [[0.035,0,0],[.083,0,0]]

templates= []
templates.append(clumpTemplate(relRadii=relRadList1,relPositions=relPosList1))
O.bodies.replaceByClumps(templates,[1])

# APLICACION DE FUERZAS ********************************************************
listSphereA=[]
for x1 in O.bodies:
 if isinstance(x1.shape,Sphere) and x1.state.pos[1]<3*r and x1.state.pos[2]<=H and x1.state.pos[2]>3*H/4.000:
  listSphereA.append(x1.id)

ha=1*H/8.000
m=EmpujeBottom-EmpujeTop
Ea=m*ha+EmpujeTop

Fa=Ea*(H/4.000)*L/len(listSphereA)

def aplicarFuerzaA():
 for i in listSphereA:
  O.forces.setPermF(i,(0,Fa,0))

# CONDICIONES DE BORDE *********************************************************

Borde1=utils.wall(0.43,axis=0, sense=0)
Borde2=utils.wall(L-0.43,axis=0, sense=0)
Borde3=utils.wall(0, axis=1, sense=0)
Borde4=utils.wall(0.05, axis=2, sense=0)

O.bodies.append(Borde1)
O.bodies.append(Borde2)
O.bodies.append(Borde3)
O.bodies.append(Borde4)

# BLOCKING ROTATION******************************************************
listSphereTodo=[]
for m in O.bodies:
 if isinstance(m.shape,Sphere) and m.state.pos[2]<2*H:
  listSphereBase.append(m.id)

for i in listSphereTodo:
 O.bodies[i].state.blockedDOFs='XYZ'

# ENGINE ***********************************************************************
O.engines=[
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Wall_Aabb()]),
  InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
   [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()],
   [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM()]
  ),
  NewtonIntegrator(gravity=(0,aceloutofplane,-9.81),damping=damp),
  PyRunner(firstIterRun=10000,command='aplicarFuerzaA()'),
   VTKRecorder(iterPeriod=100,recorders=['spheres','jcfpm','facets','colors'],fileName='/tmp/p1-')
]

# DETALLES FINALES *************************************************************
O.dt=0.5*PWaveTimeStep() # establece el delta de tiempo como una fraccion del tiempo critico
plot.plot(subPlots=True) # llama a los sub plots
plot.live=True # ploteo en tiempo real
qt.Controller() # abre la ventana del controlador
V=qt.View() # abre la ventada de la vista
R=yade.qt.Renderer() # llama a la renderizacion
R.bgColor=(1.,1.,1.) # definel color blan
V.screenSize = (900,900) # tamano de pantalla

There is a bug in replaceByClumps, see [1,2,3]
cheers
Jan

[1] https://gitlab.com/yade-dev/trunk/issues/46
[2] https://answers.launchpad.net/yade/+question/676746
[3] https://answers.launchpad.net/yade/+question/676382

Hello Jan,
Thank you for your answer, as i read in those links , i understand that its all about a bug "I think it is a bug for this contact model and it should be modified from the source." but the reason why i wrote this question its because i dont know how to delete that error. The model with spheres works correctly but once, i tried the replace clumps, that error appears. I dont have any idea how to " modify the source" as the comments in the links mention.

Thank you in advance, hope i can get a solution, i was trying to replace the word "JCFpmMat" to "JCFpmState" but the problem is that with " state",the code doesnt work . I truly dont know how to solve it. How can i solve the bug in replaceByClumps? Please

There are several solutions, e.g.:

1) you or somebody else will fix the bug in Yade source code (unlikely according to your reaction and the time the issue is there unfixed)
2) use another method than using replaceByClumps directly, e.g. rewriting it to python
3) reorder the script not to use JCFpmMat before you need it

I think 3) is the way to go. In the script, you actually do not need the JCFpmMat (which is the source of the error in replaceByClumps) at all, so you can define / set it to bodies just before O.step() / O.run() / clicking play button.
something like:
###
... # no JCFpmMat yet!
Esferas=pack.randomDensePack(Vol4,radius=r,rRelFuzz=0)
O.bodies.append(Esferas)
...
O.bodies.replaceByClumps(templates,[.5,.3,.1,.1])
mat1 = JCFpmMat()
...
O.materials.append(mat1)
...
... # the end of your script
for b in O.bodies:
   b.mat = mat1
###

cheers
Jan

thank you so much for your answer, i try to incorporate what you suggest . Even though, it continues appearing " In [1]: FATAL /build/yade-fDuCoe/yade-2018.02b/core/ThreadRunner.cpp:30 run: Exception occured:
Body #573: Body::material type JCFpmMat doesn't correspond to Body::state type State (should be JCFpmState instead)."

I would like to know if i have to change also the engine? or the problem its related to the engine with JCFPm type? or is there any other thing im doing wrong in the script? please, hope i can solve this error. Thank you in advance, im very thankful for the responses.

im copying the script with the modifications
from yade import pack, qt, plot
from numpy import arange
import numpy as np
import random

#DATOS DE ENTRADA
r=0.07
H=1 #altura del muro
B=.45 #base mayor
b=.25 #base menor
L=4 #longitud del muro
#EN SERVICIO
EmpujeTop=2400#Pa2400
EmpujeBottom=9100#Pa9100
damp=0.7
aceloutofplane=0
ii=0

# GEOMETRIA ********************************************************************
a=(B-b)*.5
Vol1=pack.inParallelepiped(o=(0,0,0), a=(L,0,0), b=(0,B,0), c=(0,0,H))
Vol2=pack.inParallelepiped(o=(L,B+a,0), a=(0,B+a,0), b=(0,a+b,H), c=(L,a+b,H))
Vol3=pack.inParallelepiped(o=(L,-a,0), a=(0,-a,0), b=(0,a,H), c=(L,a,H))
Vol4=Vol1-Vol2-Vol3-Vol3
Esferas=pack.randomDensePack(Vol4,radius=r,rRelFuzz=0)
O.bodies.append(Esferas)

relRadList1 = [.04,.03]
relPosList1 = [[0.035,0,0],[.083,0,0]]
#peanut:
relRadList2 = [.05,.05]
relPosList2 = [[0,0,0],[.073,0,0]]

relRadList4 = [0.025,0.025,0.025,0.025]
relPosList4 = [[0,0,0],[.043,0,0],[0.02,0.045,0],[0.02,0,045]]
#stick:
relRadList3 = [0.05,0.025,0.05,0.025]
relPosList3 = [[.04,0,0],[.123,0,0.022],[0.14,0,0],[0.123,0,-0.022]]
templates= []
templates.append(clumpTemplate(relRadii=relRadList1,relPositions=relPosList1))
templates.append(clumpTemplate(relRadii=relRadList2,relPositions=relPosList2))
templates.append(clumpTemplate(relRadii=relRadList3,relPositions=relPosList3))
templates.append(clumpTemplate(relRadii=relRadList4,relPositions=relPosList4))
O.bodies.replaceByClumps(templates,[.5,.3,.1,.1])

# DEFINICION DEL MATERIAL ******************************************************
mat1 = JCFpmMat()
mat1.cohesion =0 #Pa32.37e3
mat1.density = 2821 #kg/m3
mat1.frictionAngle = radians(58.74) #rad
mat1.poisson = 0.18
mat1.young = 37.04e9 #Pa

# DEFINICION DE LA JUNTA ******************************************************
mat1.jointCohesion = 262000 #Pa
mat1.jointFrictionAngle = radians(40)
mat1.jointNormalStiffness= 1.75e8 #Pa/m
mat1.jointShearStiffness=0.7e8 #Pa/m
#mat1.jointDilationAngle = 0
#mat1.jointTensileStrength = 0
O.materials.append(mat1)

# APLICACION DE FUERZAS ********************************************************
listSphereA=[]
listSphereB=[]
listSphereC=[]
listSphereD=[]
for x1 in O.bodies:
 if isinstance(x1.shape,Sphere) and x1.state.pos[1]<3*r and x1.state.pos[2]<=H and x1.state.pos[2]>3*H/4.000:
  listSphereA.append(x1.id)
for x2 in O.bodies:
 if isinstance(x2.shape,Sphere) and x2.state.pos[1]<3*r and x2.state.pos[2]<=3*H/4.000 and x2.state.pos[2]>H/2.000:
  listSphereB.append(x2.id)
for x3 in O.bodies:
 if isinstance(x3.shape,Sphere) and x3.state.pos[1]<3.5*r and x3.state.pos[2]<=H/2.000 and x3.state.pos[2]>H/4.000:
  listSphereC.append(x3.id)
for x4 in O.bodies:
 if isinstance(x4.shape,Sphere) and x4.state.pos[1]<3*r and x4.state.pos[2]<=H/4.000 and x4.state.pos[2]>0.000:
  listSphereD.append(x4.id)
ha=1*H/8.000
hb=3*H/8.000
hc=5*H/8.000
hd=7*H/8.000
m=EmpujeBottom-EmpujeTop
Ea=m*ha+EmpujeTop
Eb=m*hb+EmpujeTop
Ec=m*hc+EmpujeTop
Ed=m*hd+EmpujeTop

Fa=Ea*(H/4.000)*L/len(listSphereA)
Fb=Eb*(H/4.000)*L/len(listSphereB)
Fc=Ec*(H/4.000)*L/len(listSphereC)
Fd=Ed*(H/4.000)*L/len(listSphereD)

def aplicarFuerzaA():
 for i in listSphereA:
  O.forces.setPermF(i,(0,Fa,0))
def aplicarFuerzaB():
 for i in listSphereB:
  O.forces.setPermF(i,(0,Fb,0))
def aplicarFuerzaC():
 for i in listSphereC:
  O.forces.setPermF(i,(0,Fc,0))
def aplicarFuerzaD():
 for i in listSphereD:
  O.forces.setPermF(i,(0,Fd,0.98))

# CONDICIONES DE BORDE *********************************************************

Borde1=utils.wall(0.43,axis=0, sense=0)
Borde2=utils.wall(L-0.43,axis=0, sense=0)
Borde3=utils.wall(0, axis=1, sense=0)
Borde4=utils.wall(0.05, axis=2, sense=0)

O.bodies.append(Borde1)
O.bodies.append(Borde2)
O.bodies.append(Borde3)
O.bodies.append(Borde4)

listSphereBase=[]
#listSphereBorde1=[]
#listSphereBorde2=[]
#restringir movimiento en la base
#for y1 in O.bodies:
# if isinstance(y1.shape,Sphere) and y1.state.pos[2]<=2.5*r:
# listSphereBase.append(y1.id)
#for y2 in O.bodies:
# if isinstance(y2.shape,Sphere) and y2.state.pos[0]<=2.5*r:
# listSphereBorde1.append(y2.id)
#for y3 in O.bodies:
# if isinstance(y3.shape,Sphere) and y3.state.pos[0]>L-2.5*r:
# listSphereBorde2.append(y3.id)
#for i in listSphereBase:
# O.bodies[i].state.blockedDOFs='xyz'
#for i in listSphereBorde1:
# O.bodies[i].state.blockedDOFs='xyzXYZ'
#for i in listSphereBorde2:
# O.bodies[i].state.blockedDOFs='xyzXYZ'

# DEFINICION DE FUNCIONES ******************************************************
listSphereTodo=[]
for m in O.bodies:
 if isinstance(m.shape,Sphere) and m.state.pos[2]<2*H:
  listSphereBase.append(m.id)

for i in listSphereTodo:
 O.bodies[i].state.blockedDOFs='XYZ'

def checkUnbalanced():
 print ('iter %d, unbalanced forces = %f'%(O.iter, utils.unbalancedForce()))
 print 'Fzas en c/ esfera (N)', Fa, Fb, Fc, Fd
# if unbalancedForce()<.00001:
 print ( 'Desplazamiento superior',O.bodies[w].state.pos[1],O.bodies[w].state.pos[2])

listTopBodies=[]
for z1 in O.bodies:
 if isinstance(z1.shape,Sphere) and z1.state.pos[2]<=H and z1.state.pos[2]>H-3.000*r and z1.state.pos[0]>(L/2.000)-1.1*r and z1.state.pos[0]<(L/2.000)+1.1*r and z1.state.pos[1]>B/2.000:
  listTopBodies.append(z1.id)
w=max(listTopBodies)

def addPlotData():
 plot.addData(Ypos=O.bodies[w].state.pos[1],Zpos=O.bodies[w].state.pos[2],iteracion=O.iter,unbal=unbalancedForce())
plot.plots={'Ypos':('Zpos')}

# ENGINE ***********************************************************************
O.engines=[
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Wall_Aabb()]),
  InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
   [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()],
   [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM()]
  ),
  NewtonIntegrator(gravity=(0,aceloutofplane,-9.81),damping=damp),
  PyRunner(firstIterRun=10000,command='aplicarFuerzaA()'),
  PyRunner(firstIterRun=10000,command='aplicarFuerzaB()'),
  PyRunner(firstIterRun=10000,command='aplicarFuerzaC()'),
  PyRunner(firstIterRun=10000,command='aplicarFuerzaD()'),
  PyRunner(iterPeriod=10000,command='checkUnbalanced()'),
  PyRunner(iterPeriod=1,command='addPlotData()'),
  VTKRecorder(iterPeriod=100,recorders=['spheres','jcfpm','facets','colors'],fileName='/tmp/p1-')
]

for b in O.bodies:
 b.mat=mat1
# DETALLES FINALES *************************************************************

O.dt=0.5*PWaveTimeStep() # establece el delta de tiempo como una fraccion del tiempo critico
plot.plot(subPlots=True) # llama a los sub plots
plot.live=True # ploteo en tiempo real
qt.Controller() # abre la ventana del controlador
V=qt.View() # abre la ventada de la vista
R=yade.qt.Renderer() # llama a la renderizacion
R.bgColor=(1.,1.,1.) # definel color blan
V.screenSize = (900,900) # tamano de pantalla

Sorry, I missed that you also have to change state, not only material:
###
for b in O.bodies:
   b.mat=mat1
   b.state=JCFpmState()
###

cheers
Jan

Jan, thank you for the answers . Ive already change that part, but when running it , what appears is only one sphere and not with the packing in the volume i made, not like before when geometrically it was working well but appeared the "fatal error" . So sorry for the replies. Dont know how to make it work.

Im copying the script in case i can have a hand with it.

from yade import pack, qt, plot
from numpy import arange
import numpy as np
import random

#DATOS DE ENTRADA
r=0.07
H=1 #altura del muro
B=.45 #base mayor
b=.25 #base menor
L=4 #longitud del muro
#EN SERVICIO
EmpujeTop=2400#Pa2400
EmpujeBottom=9100#Pa9100
damp=0.7
aceloutofplane=0
ii=0

# GEOMETRIA ********************************************************************
a=(B-b)*.5
Vol1=pack.inParallelepiped(o=(0,0,0), a=(L,0,0), b=(0,B,0), c=(0,0,H))
Vol2=pack.inParallelepiped(o=(L,B+a,0), a=(0,B+a,0), b=(0,a+b,H), c=(L,a+b,H))
Vol3=pack.inParallelepiped(o=(L,-a,0), a=(0,-a,0), b=(0,a,H), c=(L,a,H))
Vol4=Vol1-Vol2-Vol3-Vol3
Esferas=pack.randomDensePack(Vol4,radius=r,rRelFuzz=0)
O.bodies.append(Esferas)

relRadList1 = [.04,.03]
relPosList1 = [[0.035,0,0],[.083,0,0]]
#peanut:
relRadList2 = [.05,.05]
relPosList2 = [[0,0,0],[.073,0,0]]

relRadList4 = [0.025,0.025,0.025,0.025]
relPosList4 = [[0,0,0],[.043,0,0],[0.02,0.045,0],[0.02,0,045]]
#stick:
relRadList3 = [0.05,0.025,0.05,0.025]
relPosList3 = [[.04,0,0],[.123,0,0.022],[0.14,0,0],[0.123,0,-0.022]]
templates= []
templates.append(clumpTemplate(relRadii=relRadList1,relPositions=relPosList1))
templates.append(clumpTemplate(relRadii=relRadList2,relPositions=relPosList2))
templates.append(clumpTemplate(relRadii=relRadList3,relPositions=relPosList3))
templates.append(clumpTemplate(relRadii=relRadList4,relPositions=relPosList4))
O.bodies.replaceByClumps(templates,[.5,.3,.1,.1])

# DEFINICION DEL MATERIAL ******************************************************
mat1 = JCFpmMat()
mat1.cohesion =0 #Pa32.37e3
mat1.density = 2821 #kg/m3
mat1.frictionAngle = radians(58.74) #rad
mat1.poisson = 0.18
mat1.young = 37.04e9 #Pa

# DEFINICION DE LA JUNTA ******************************************************
mat1.jointCohesion = 262000 #Pa
mat1.jointFrictionAngle = radians(40)
mat1.jointNormalStiffness= 1.75e8 #Pa/m
mat1.jointShearStiffness=0.7e8 #Pa/m
O.materials.append(mat1)

# APLICACION DE FUERZAS ********************************************************
listSphereA=[]
listSphereB=[]
listSphereC=[]
listSphereD=[]
for x1 in O.bodies:
 if isinstance(x1.shape,Sphere) and x1.state.pos[1]<3*r and x1.state.pos[2]<=H and x1.state.pos[2]>3*H/4.000:
  listSphereA.append(x1.id)
for x2 in O.bodies:
 if isinstance(x2.shape,Sphere) and x2.state.pos[1]<3*r and x2.state.pos[2]<=3*H/4.000 and x2.state.pos[2]>H/2.000:
  listSphereB.append(x2.id)
for x3 in O.bodies:
 if isinstance(x3.shape,Sphere) and x3.state.pos[1]<3.5*r and x3.state.pos[2]<=H/2.000 and x3.state.pos[2]>H/4.000:
  listSphereC.append(x3.id)
for x4 in O.bodies:
 if isinstance(x4.shape,Sphere) and x4.state.pos[1]<3*r and x4.state.pos[2]<=H/4.000 and x4.state.pos[2]>0.000:
  listSphereD.append(x4.id)
ha=1*H/8.000
hb=3*H/8.000
hc=5*H/8.000
hd=7*H/8.000
m=EmpujeBottom-EmpujeTop
Ea=m*ha+EmpujeTop
Eb=m*hb+EmpujeTop
Ec=m*hc+EmpujeTop
Ed=m*hd+EmpujeTop

Fa=Ea*(H/4.000)*L/len(listSphereA)
Fb=Eb*(H/4.000)*L/len(listSphereB)
Fc=Ec*(H/4.000)*L/len(listSphereC)
Fd=Ed*(H/4.000)*L/len(listSphereD)

def aplicarFuerzaA():
 for i in listSphereA:
  O.forces.setPermF(i,(0,Fa,0))
def aplicarFuerzaB():
 for i in listSphereB:
  O.forces.setPermF(i,(0,Fb,0))
def aplicarFuerzaC():
 for i in listSphereC:
  O.forces.setPermF(i,(0,Fc,0))
def aplicarFuerzaD():
 for i in listSphereD:
  O.forces.setPermF(i,(0,Fd,0.98))

# CONDICIONES DE BORDE *********************************************************

Borde1=utils.wall(0.43,axis=0, sense=0)
Borde2=utils.wall(L-0.43,axis=0, sense=0)
Borde3=utils.wall(0, axis=1, sense=0)
Borde4=utils.wall(0.05, axis=2, sense=0)

O.bodies.append(Borde1)
O.bodies.append(Borde2)
O.bodies.append(Borde3)
O.bodies.append(Borde4)

listSphereBase=[]
for i in listSphereBase:
  O.bodies[i].state.blockedDOFs='xyz'

# DEFINICION DE FUNCIONES ******************************************************
for m in O.bodies:
 if isinstance(m.shape,Sphere) and m.state.pos[2]<2*H:
  listSphereBase.append(m.id)
listSphereTodo=[]
for i in listSphereTodo:
 O.bodies[i].state.blockedDOFs='XYZ'

def checkUnbalanced():
 print ('iter %d, unbalanced forces = %f'%(O.iter, utils.unbalancedForce()))
 print 'Fzas en c/ esfera (N)', Fa, Fb, Fc, Fd
 print ( 'Desplazamiento superior',O.bodies[w].state.pos[1],O.bodies[w].state.pos[2])

listTopBodies=[]
for z1 in O.bodies:
 if isinstance(z1.shape,Sphere) and z1.state.pos[2]<=H and z1.state.pos[2]>H-3.000*r and z1.state.pos[0]>(L/2.000)-1.1*r and z1.state.pos[0]<(L/2.000)+1.1*r and z1.state.pos[1]>B/2.000:
  listTopBodies.append(z1.id)
w=max(listTopBodies)

def addPlotData():
 plot.addData(Ypos=O.bodies[w].state.pos[1],Zpos=O.bodies[w].state.pos[2],iteracion=O.iter,unbal=unbalancedForce())
plot.plots={'Ypos':('Zpos')}

# ENGINE ***********************************************************************
O.engines=[
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Wall_Aabb()]),
  InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
   [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()],
   [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM()]
  ),
  NewtonIntegrator(gravity=(0,aceloutofplane,-9.81),damping=damp),
  PyRunner(firstIterRun=10000,command='aplicarFuerzaA()'),
  PyRunner(firstIterRun=10000,command='aplicarFuerzaB()'),
  PyRunner(firstIterRun=10000,command='aplicarFuerzaC()'),
  PyRunner(firstIterRun=10000,command='aplicarFuerzaD()'),
  PyRunner(iterPeriod=10000,command='checkUnbalanced()'),
  PyRunner(iterPeriod=1,command='addPlotData()'),
  VTKRecorder(iterPeriod=100,recorders=['spheres','jcfpm','facets','colors'],fileName='/tmp/p1-')
]

for b in O.bodies:
 b.mat=mat1
 b.state=JCFpmState()

# DETALLES FINALES *************************************************************

O.dt=0.5*PWaveTimeStep() # establece el delta de tiempo como una fraccion del tiempo critico
plot.plot(subPlots=True) # llama a los sub plots
plot.live=True # ploteo en tiempo real
qt.Controller() # abre la ventana del controlador
V=qt.View() # abre la ventada de la vista
R=yade.qt.Renderer() # llama a la renderizacion
R.bgColor=(1.,1.,1.) # definel color blan
V.screenSize = (900,900) # tamano de pantalla

Sorry, I have almost no experience with clumps, you have to wait for somebody else to answer..
Or maybe you can try to create a new question (this one is already a bit lengthy..)
cheers
Jan

anyway, thanks Jan for the help and for your constant reply!

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