Question #683526 on Yade changed:
https://answers.launchpad.net/yade/+question/683526
Status: Needs information => Open
Akm gave more information on the question:
from yade import pack, plot, qt
O.reset()
#Creating the upper box
U_box=geom.facetBox((.2,.1,.8),(.2,.1,.2),wallMask=51)
O.bodies.append(U_box)
#Creating the lower box
L_box=geom.facetBox((.2,-.1,.8),(.2,.1,.2),wallMask=55)
O.bodies.append(L_box)
'''
#Creating the upper box
U_left=utils.box((0,.1,.8),(0,.1,.2),dynamic=False,fixed=True,wire=True,color=None)
O.bodies.append(U_left)
U_right=utils.box((.4,.1,.8),(0,.1,.2),dynamic=False,fixed=True,wire=True,color=None)
O.bodies.append(U_right)
U_front=utils.box((.2,.1,1),(.2,.1,0),dynamic=False,fixed=True,wire=True,color=None)
O.bodies.append(U_front)
U_back=utils.box((.2,.1,.6),(.2,.1,0),dynamic=False,fixed=True,wire=True,color=None)
O.bodies.append(U_back)
#Creating the lower box
L_left=utils.box((0,-.1,.8),(0,.1,.2),dynamic=False,fixed=True,wire=True,color=None)
O.bodies.append(L_left)
L_right=utils.box((.4,-.1,.8),(0,.1,.2),dynamic=False,fixed=True,wire=True,color=None)
O.bodies.append(L_right)
L_front=utils.box((.2,-.1,1),(.2,.1,0),dynamic=False,fixed=True,wire=True,color=None)
O.bodies.append(L_front)
L_back=utils.box((.2,-.1,.6),(.2,.1,0),dynamic=False,fixed=True,wire=True,color=None)
O.bodies.append(L_back)
L_bottom
=utils.box((.2,-0.05,.8),(.2,0,.2),dynamic=False,fixed=True,wire=True,color=None)
O.bodies.append(L_bottom)
#Creating the extension for sliding plane so particles don't fall out
#Give extra .01m in Y dir. so that the two facet walls...
#...at the shear interface do not rub against each other
U_box_extension1=geom.facetBox((-.2,.5,.8),(.2,.5,.2),wallMask=4)
O.bodies.append(U_box_extension1)
U_box_extension2=geom.facetBox((.6,.5,.8),(.2,.5,.2),wallMask=4)
O.bodies.append(U_box_extension2)
#Creating the extension for sliding plane so particles don't fall out
L_box_extension1=geom.facetBox((.6,-.5001,.8),(.2,.5,.2),wallMask=8)
O.bodies.append(L_box_extension1)
L_box_extension2=geom.facetBox((-.2,-.5001,.8),(.2,.5,.2),wallMask=8)
O.bodies.append(L_box_extension2)
#Material assignment- normalcohesion = Tensile strength, shearCohesion = Shear
strength
sample_material=FrictMat(
young=50e6 ,poisson=0.25 ,density=1800
,frictionAngle=radians(45),label='sample_mat')
O.materials.append(sample_material)
wall_material= FrictMat(
young=4e9 ,poisson=0.25 ,density=5000 ,frictionAngle=radians(50)
,label='wall_mat')
O.materials.append(wall_material)
#Speed of the translation engine
Vel=1
#Sphere Pack Creation
sp=pack.SpherePack()
sp.makeCloud((0,-0.05,0.6),(0.4,0.2,1.0),rMean=0.008,rRelFuzz=0.3,porosity=0.5,
num = 6000)
sp.toSimulation(material=sample_material) #already one cloud of particles here
#Collecting the ids for translation engine run
id_list=[]
#for i in L_box: id_list.append(i.id)#for facet box
id_list.append(L_left.id)
id_list.append(L_right.id)
id_list.append(L_front.id)
id_list.append(L_back.id)
id_list.append(L_bottom.id)
for i in L_box_extension1: id_list.append(i.id)
for i in L_box_extension2: id_list.append(i.id)
#Simulation
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom6D()],
[Ip2_FrictMat_FrictMat_FrictPhys(),
Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=False,label="cohesiveIp")
],
[Law2_ScGeom_FrictPhys_CundallStrack(),
Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
),
NewtonIntegrator(gravity=(0,-1,0),damping=0.3),
GlobalStiffnessTimeStepper(active=True,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.8),
#using global timestep for faster simulation
# (command='checkUnbalanced()',realPeriod=0.5,label='checker'),
PyRunner(command='addplate()',iterPeriod=10,label='runner') #real period is
actually not good
]
maxLoad=1e5
minLoad=1e1
inst_position_lower=min([b.state.pos[1]-(b.shape.radius) for b in
O.bodies if isinstance(b.shape,Sphere)])
global plate_b
plate_b=wall(min([b.state.pos[1]-(b.shape.radius) for b in O.bodies if
isinstance(b.shape,Sphere)]),axis=1,sense=+1)
#O.bodies.append(plate_b)
print('Plate in -Y added')
#Plate is being added to the simulation to consolidate the sample
def addplate():
#if unbalancedForce()<0.5 and O.iter>15000:
global plate
inst_position1=max([b.state.pos[1]+(b.shape.radius) for b in O.bodies
if isinstance(b.shape,Sphere)])
plate=utils.box((.2,inst_position1,.8),(.199,0,.199),dynamic=False,fixed=False,wire=False,color=(0,0,1))#,material=wall_material)
O.bodies.append(plate)
print('Plate added')
plate.state.vel=(0,-2,0)
runner.command='unloadPlate()'
#Plate rebounding after consolidation upto a maximum load
def unloadPlate():
