Plate for three point bending load application

Dear all,

I'm trying to simulate three-point bending test in Yade.
I have prepared a sample with random dense pack.

In order to apply the load as a piston, I need to have a plate with 3mm in width and 50mm length parallel to xz plane in the middle of the beam.

and another two plates as rollers with the same dimensions at the bottom left and right locations?

How can I add a plate with a specific dimension?

Previously I have used facetcylinder to apply load as well as to use as rollers.

Now I need to use a single plate instead of this facet cylinder. could you please tell me how to do that?

this is my existing code,

# -*- coding: utf-8 -*-
# Code for three point bending test with crack mouth opening displacement and deflection
#Developed by S.M.C.U. Senanayake, Department of Civil Engineering, Monash University, Australia

from yade import pack

############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################

from yade import pack, plot,qt,ymport
import os
import csv

# Damping:
damp=0.7
stabilityThreshold=0.001 # we test unbalancedForce against this value in different loops (see below)
# material parameters
#young = 6.75e8
#poisson = .2
#frictionAngle = 0.6
#sigmaT = 1.5e6
#epsCrackOnset = 16.4e-3
#relDuctility = 1.05
#density = 2600

young = 3e9
poisson = .3
frictionAngle = 0.5
sigmaT = .3e6
epsCrackOnset = .9e-3
relDuctility = 30
density = 2600
# prestress
# axial strain rate
rate=0.01 # loading rate (strain rate)
 # stress path
'''
# 3 point bending test sample preparation parameters
L= length of the sample
W= width of the sample
H= height of the sample
C=clearance from two ends(canteliver length
A= width of the notch
B= height of the notch
R=radius of the rollers/piston(radius of facet clyinder) all mm s (10^-3) scale
Length along x axis, height along y axis and width along z axis
'''
#mn,mx=Vector3(0,0,0),Vector3(1,1,1) # corners of the initial packing

L=500e-3
W=50e-3
H=100e-3
C=50e-3
A=3e-3
B=50e-3
R=7e-3
psdmean=3e-3
rmin=1.5*psdmean
prepared=0
v=.1 # applied velocity on piston

concMat=O.materials.append(CpmMat(young=young,poisson=poisson,density=density,damLaw=1,frictionAngle=frictionAngle,
sigmaT=sigmaT,epsCrackOnset=epsCrackOnset,neverDamage=False,label='concrete',relDuctility=relDuctility,))

frictMat = O.materials.append(CpmMat(young=100*young,poisson=poisson,frictionAngle=0,sigmaT=0,epsCrackOnset=epsCrackOnset,relDuctility=1.0001,neverDamage=True,density=3000,label='walls'))

if prepared<0.5:
 pred=pack.inAlignedBox((0,0,0),(L,H,W))-pack.inAlignedBox((.5*(L-A),0,0),(.5*(L+A),B,W))
 #SpherePack=pack.randomDensePack(pred,spheresInCell=11000,radius=psdmean,rRelFuzz=0.2,color=(0,0,1),material=concMat,returnSpherePack=True)
 SpherePack=pack.randomDensePack(pred,spheresInCell=500,radius=psdmean,rRelFuzz=0.2,color=(0,0,1),material=concMat,returnSpherePack=True)
 sand=SpherePack.toSimulation() # assign all particles to sand object
else:
 O.bodies.append(ymport.textExt('tpb_sample.txt',format='x_y_z_r',material=concMat))
 print ('I have previously prepared packing!!')
for b in O.bodies:
 if isinstance(b.shape,Sphere):
  b.shape.color = (0.2,0.6,0.6)# give a colour
# In order to extract the ids of each particle in the left and right hand corner of the notch and make a list out of them
leftids=[]
for i in O.bodies:
 if i.state.pos[0]>(0.5*(L-A)-2*rmin) and i.state.pos[0]<(0.5*(L-A)):
  if i.state.pos[1]>0 and i.state.pos[1]<(2*rmin):
   i.shape.color=(1,0,0)
   leftids.append(i.id)
   #left=[O.bodies[i] in sand]

rightids=[]
c=0
for i in O.bodies:
 if i.state.pos[0]>(0.5*(L+A)) and i.state.pos[0]<(0.5*(L+A)+2*rmin):
  if i.state.pos[1]>0 and i.state.pos[1]<(2*rmin):
   i.shape.color=(1,0,0)
   c=c+1
   print i
   print i.id
   rightids.append(i.id)
   #right=[O.bodies[i] in sand]
print c
print rightids
#print right
# assign to a lh object all particles in the left hand side interested region of the notch
lh = [O.bodies[s] for s in sand if (O.bodies[s].state.pos[0]>(0.5*(L-A)-2*rmin) and O.bodies[s].state.pos[0]<(0.5*(L-A)) and O.bodies[s].state.pos[1]>0 and O.bodies[s].state.pos[1]<(2*rmin))]
rh=[O.bodies[s] for s in sand if (O.bodies[s].state.pos[0]>(0.5*(L+A)) and O.bodies[s].state.pos[0]<(0.5*(L+A)+2*rmin) and O.bodies[s].state.pos[1]>0 and O.bodies[s].state.pos[1]<(2*rmin))]
for i in rh:
 print i.id
# assign global variable to total numbe rof particles in each block
global left_part
left_part=len(leftids)
global right_part
right_part=len(rightids)

piston=O.bodies.append(geom.facetCylinder(center=(.5*L,H+R,.5*W),radius=R,height=W,orientation=Quaternion((0, 0, 1), pi/2),segmentsNumber=20,wire=False,material=frictMat)) #piston radius is 15mm

left_roller=O.bodies.append(geom.facetCylinder(center=(C,-R,.5*W),radius=R,height=W,orientation=Quaternion((0, 0, 1), pi/2),segmentsNumber=20,wire=False,material=frictMat))

right_roller=O.bodies.append(geom.facetCylinder(center=(L-C,-R,.5*W),radius=R,height=W,orientation=Quaternion((0, 0, 1), pi/2),segmentsNumber=20,wire=False,material=frictMat))
newton=NewtonIntegrator(damping=damp)
enlargeFactor=1.5
for b in O.bodies:
 if isinstance(b.shape,Sphere):
  b.shape.color = (0.2,0.6,0.6)
print len(O.bodies)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([
  Bo1_Sphere_Aabb(aabbEnlargeFactor=enlargeFactor,label='bo1s'),
  Bo1_Facet_Aabb()
 ]),
 InteractionLoop(
  [
   Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=enlargeFactor,label='ig2ss'),
   Ig2_Facet_Sphere_ScGeom()
  ],
  [
   Ip2_CpmMat_CpmMat_CpmPhys(cohesiveThresholdIter=O.iter+5),
   Ip2_FrictMat_CpmMat_FrictPhys(),
   Ip2_FrictMat_FrictMat_FrictPhys(),
  ],
  [
   Law2_ScGeom_CpmPhys_Cpm(),
   Law2_ScGeom_FrictPhys_CundallStrack(),
  ],
 ),
 ## We will use the global stiffness of each body to determine an optimal timestep (see https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
 #TranslationEngine(ids=[piston.id],translationAxis=[0,0,1],velocity=1),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.3),
 #VTKRecorder(fileName='/home/ud/vtkdata/3pb_test', recorders=['all'], iterPeriod=100),
 PyRunner(iterPeriod=1,command='apply()',label='apply'),
 PyRunner(iterPeriod=1,command='addPlotData()',label='graph'),
 #PyRunner(iterPeriod=50,command='pr()',label='ph'),
 newton,
 #CpmStateUpdater(iterPeriod=100),
]

#Display spheres with 2 colors for seeing rotations better
Gl1_Sphere.stripes=0
qtr=yade.qt.Renderer()
qtr.bgColor=(1,1,1)
yade.qt.Controller(), yade.qt.View()
O.step()
bo1s.aabbEnlargeFactor=ig2ss.interactionDetectionFactor=1.
cylinderIds=[]
for i in range(0,len(piston)):
 cylinderIds.append(piston[i])
 #print piston[i]
#print cylinderIds

def apply():
 for i in range(0,len(piston)):
  O.bodies[piston[i]].state.vel[1]=-v
# a function to find crack mouth opening displacement
def cmod():
 sumlh=0
 sumrh=0
 for i in lh:
  dp_left=i.state.pos[0]-i.state.refPos[0]# left displacment
  sumlh+=dp_left # addig all the displacement in x direction components

 for j in rh:
  dp_right=j.state.pos[0]-j.state.refPos[0]
  sumrh+=dp_right

 d1=sumlh/left_part # average displacment in left hand side of the notch
 d2=sumrh/right_part
 cmodi=abs(d1)+abs(d2) # cmod by adding right ahand side and left hand side
 return (cmodi)

def addPlotData():
 f_2=utils.sumForces(ids=cylinderIds,direction =(0,1,0))# measure total force induced on the cylinder
 delta_z=v*O.time
 crackopen=cmod()
 #print f_2
 plot.addData(t = O.time,f_2 = f_2,delta_z=delta_z,crackopen=crackopen)

#O.engines=O.engines+[PyRunner(command='pr()',iterPeriod=50)]
def pr():
 presentcohesive_count = 0
 for i in O.interactions:
  if hasattr(i.phys, 'isCohesive'):
   if i.phys.isCohesive == 0:
    presentcohesive_count+=1
    a=i.id1
    b=i.id2
    O.bodies[a].shape.color = (1,0,0)
    O.bodies[b].shape.color = (1,0,0)
    #print (i.id1, "-", i.id2)
 noncohesive_count=presentcohesive_count

#O.engines=O.engines+[TranslationEngine(ids=[piston.id],translationAxis=[0,0,1],velocity=1)]
#plot.plots={'delta_z':('f_2',), 'crackopen':('f_2 ')}
plot.plots={'crackopen':('f_2',)}
plot.plot()

#print (len(O.bodies))
from yade import qt
qt.Controller()
qt.View()