# Clumps pack sphere in a irregular volume

Asked by Grace Mejico on 2019-10-25

Good evening,
I'm from Peru and I'm part of an investigation that tries to represent the pirca performance at seismic action. The pirca is a dry stone retaining very used here.
For making the model, I was trying to use spherepack with clumps to represent the blocks of rock form.

First, If i use a volume like a regular prism using vector,it works perfectly but what i want is a trapezoidal form (cross section) of the dry stone retaining wall. I tried many times but i couldnt find a way to make it work,thats why im sharing the code i was working on and hopefully i can have a hand with this problem. Is there a way to put the clumps in a volume with a trapezoidal form as cross section?

Second,I also want to recognize the clump that has the greatest displacement thats why i wrote that code . I would like to know if there is another way to do it.

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

plot.plots={'Ypos':('Zpos'),'iteracion':('unbal')}

--------------------------------------------------------------------------------------------------------------------------------------------------------------------
Here i share the complete code for the model:

# MODELO NUMERICO CON CLUMPS *************************
from yade import pack, qt, plot
from numpy import arange
import numpy as np
import random
r=0.07
H=1 #altura del muro
B=.45 #base mayor
b=.25 #base menor
L=2 #longitud del muro
#EN SERVICIO
EmpujeTop=2400 #Pa
EmpujeBottom=9100 #Pa
damp=0.7
aceloutofplane=0
ii=0
# Definicion del material mat1 para la pirca **********
mat1 = JCFpmMat()
mat1.cohesion = 0 #32.37e3 #Pa
mat1.density = 2233.2 #kg/m3
mat1.poisson = 0.18
mat1.young = 37.04e9 #Pa

# DEFINICION DE LA JUNTA ******************************************************
mat1.jointCohesion = 232500
mat1.jointNormalStiffness= 1.75e8
mat1.jointShearStiffness=0.875e8
#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

sp=pack.SpherePack()
biesfera=pack.SpherePack([((0,0,0),.033),((.033,0,0),.033)]) # clump de dos esferas
triesfera=pack.SpherePack([((0,0,0),.04),((.02,.02,0),.04),((0,.02,.02),.04)]) # clump de tres esferas
tetraesfera=pack.SpherePack([((0,0,0),.033),((.033,0,0),.033),((0,.033,0),.033),((0,0,.033),.033)]) # clump decuatro esferas
sp.makeClumpCloud(Vol4,[biesfera,triesfera, tetraesfera]) # Empaquetamiento suelto y aleatoriode clumps
sp.toSimulation() # llama las esferas a la simulacion
#Imprime la cantidad total de cuerpos
Q=len(O.bodies)
print 'numero de cuerpos', Q
# Definicion del piso
piso=utils.wall(position=Vector3(0,0,0), axis=2, sense=0, material=mat1, color=(0.5,0.5,0.5))
O.bodies.append(piso)

# Definicion de las fuerzas en aplicacion
listSphereA=[]
listSphereB=[]
listSphereC=[]
listSphereD=[]
for x1 in O.bodies:
if isinstance(x1.shape,Sphere) and x1.state.pos<3*r and x1.state.pos<=H and x1.state.pos>3*H/4.000:
listSphereA.append(x1.id)
for x2 in O.bodies:
if isinstance(x2.shape,Sphere) and x2.state.pos<3*r and x2.state.pos<=3*H/4.000+0.1 and x2.state.pos>H/2.000+0.1:
listSphereB.append(x2.id)
for x3 in O.bodies:
if isinstance(x3.shape,Sphere) and x3.state.pos<3*r and x3.state.pos<=H/2.000+0.1 and x3.state.pos>H/4.000+0.1:
listSphereC.append(x3.id)
for x4 in O.bodies:
if isinstance(x4.shape,Sphere) and x4.state.pos<3*r and x4.state.pos<=H/4.000+0.1 and x4.state.pos>0.000+0.1:
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)

print 'Fzas en c/ esfera (N)', Fa, Fb, Fc, Fd

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))
for i in listSphereD:
O.forces.setPermF(i,(0,Fd,0))

# CONDICIONES DE BORDE *********************************************************
Base=utils.wall(0.1,axis=2,sense=0)
Borde1=utils.wall(0.22, axis=0, sense=0)
Borde2=utils.wall(L-.22, axis=0, sense=0)
Borde3=utils.wall(-.0, axis=1, sense=0)
O.bodies.append(Base)
O.bodies.append(Borde1)
O.bodies.append(Borde2)
O.bodies.append(Borde3)

listSphereBase=[]

for y1 in O.bodies:
if isinstance(y1.shape,Sphere) and y1.state.pos<=3*r+0.1:
listSphereBase.append(y1.id)
#for y2 in O.bodies:
# if isinstance(y2.shape,Sphere) and y2.state.pos<=2.5*r:
# listSphereBorde1.append(y2.id)
#for y3 in O.bodies:
# if isinstance(y3.shape,Sphere) and y3.state.pos>L-2.5*r:
# listSphereBorde2.append(y3.id)
for i in listSphereBase:
O.bodies[i].state.blockedDOFs='xyzXYZ'
#for i in listSphereBorde1:
# O.bodies[i].state.blockedDOFs='xyzXYZ'
#for i in listSphereBorde2:
# O.bodies[i].state.blockedDOFs='xyzXYZ'

def checkUnbalanced():
print ('iter %d, unbalanced forces = %f'%(O.iter, utils.unbalancedForce()))
# if unbalancedForce()<.00001:

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

plot.plots={'Ypos':('Zpos'),'iteracion':('unbal')}

# MOTOR DEL ALGORITMO DEM *******************************
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(iterPeriod=10000,command='checkUnbalanced()'),
]

# 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,500) # tamano de pantalla

## Question information

Language:
English Edit question
Status:
Expired
For:
Assignee:
No assignee Edit question
Last query:
2019-10-25
2019-11-09
 Jérôme Duriez (jduriez) said on 2019-10-25: #1

Hi,

> Is there a way to put the clumps in a volume with a trapezoidal form as cross section?

Did you look at https://yade-dem.org/doc/user.html#volume-representation and the use of predicate ? Maybe it would make sense to try first with regular Spheres, before trying with clumps.

 Grace Mejico (mejicograce) said on 2019-10-25: #2

Yes, ive already done it. And it works, thats why now im trying with clumps.