triaxial test - flexible membrane

Hi everyone,
I am trying to modify an example (concrete triaxial) to do a triaxial test with a flexible membrane, but I want to apply the load using facets. I dont undertand why the facets at the top and bottom of the sample do not apply forces at the particles. The particles crossing the facets as it do not exits. I don't understand why the facets at the top and bottom of the sample do not apply forces at the particles. The particles acrossing the facets. Anyone can help me? Thanks.

################################################################################
#
# Triaxial test. Axial strain rate is prescribed and transverse prestress.
# Test is possible on prism or cylinder
# An independent c++ engine may be created from this script in the future.
#
################################################################################
from yade import pack, plot
import os

O.periodic=False
# default parameters or from table
readParamsFromTable(noTableOk=True,
 # type of test ['cyl','cube']
 testType = 'cyl',

 # material parameters
 young = 2.e9,
 youngWall=2e15,
 poisson = .4,
 frictionAngle = radians(5.),
 sigmaT = 1.5e6,
 epsCrackOnset = 1e-4,
 relDuctility = 30,

 # prestress
 preStress = -20.e6,
 # axial strain rate
 strainRate = -100,
 #velocity plates
 velPlates = 0.01,

 # assamlby parameters
 rParticle = 0.001,#0.0005, #
 width = 3.81e-2,
 height = 7.62e-2,
 bcCoeff = 5,

 # facets division
 nw = 10,#24,
 nh = 5,#15,

 # output specifications
 fileName = 'test',
 exportDir = '/tmp',
 runGnuplot = False,
 runInGui = True,
)
from yade.params.table import *
assert testType in ['cyl','cube']

# materials
sand = O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=frictionAngle,density=2600))
frictMat = O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=frictionAngle))
WallMat = O.materials.append(FrictMat(young=youngWall,poisson=0.2,frictionAngle=0.))

# spheres
pred = pack.inCylinder((0,0,0),(0,0,height),.5*width) if testType=='cyl' else pack.inAlignedBox((-.5*width,-.5*width,0),(.5*width,.5*width,height)) if testType=='cube' else None
sp=SpherePack()
sp = pack.randomDensePack(pred,spheresInCell=2000,radius=rParticle,rRelFuzz=0.8,returnSpherePack=True)
spheres=sp.toSimulation(color=(0,1,1),material=sand)

# bottom and top of specimen. Will have prescribed velocity
bot = [O.bodies[s] for s in spheres if O.bodies[s].state.pos[2]<rParticle*bcCoeff]
top = [O.bodies[s] for s in spheres if O.bodies[s].state.pos[2]>height-rParticle*bcCoeff]
vel = -strainRate*(height-rParticle*2*bcCoeff)

# facets
facets = []
if testType == 'cyl':
 rCyl2 = .5*width / cos(pi/float(nw))
 for r in xrange(nw):
  for h in xrange(nh):
   v1 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), rCyl2*sin(2*pi*(r+0)/float(nw)), height*(h+0)/float(nh) )
   v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), height*(h+0)/float(nh) )
   v3 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), height*(h+1)/float(nh) )
   v4 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), rCyl2*sin(2*pi*(r+0)/float(nw)), height*(h+1)/float(nh) )
   f1 = facet((v1,v2,v3),color=(0,0,1),material=frictMat)
   f2 = facet((v1,v3,v4),color=(0,0,1),material=frictMat)
   facets.extend((f1,f2))
elif testType == 'cube':
 nw2 = nw/4
 for r in xrange(nw2):
  for h in xrange(nh):
   v11 = Vector3( -.5*width + (r+0)*width/nw2, -.5*width, height*(h+0)/float(nh) )
   v12 = Vector3( -.5*width + (r+1)*width/nw2, -.5*width, height*(h+0)/float(nh) )
   v13 = Vector3( -.5*width + (r+1)*width/nw2, -.5*width, height*(h+1)/float(nh) )
   v14 = Vector3( -.5*width + (r+0)*width/nw2, -.5*width, height*(h+1)/float(nh) )
   f11 = facet((v11,v12,v13),color=(0,0,1),material=frictMat)
   f12 = facet((v11,v13,v14),color=(0,0,1),material=frictMat)
   v21 = Vector3( +.5*width, -.5*width + (r+0)*width/nw2, height*(h+0)/float(nh) )
   v22 = Vector3( +.5*width, -.5*width + (r+1)*width/nw2, height*(h+0)/float(nh) )
   v23 = Vector3( +.5*width, -.5*width + (r+1)*width/nw2, height*(h+1)/float(nh) )
   v24 = Vector3( +.5*width, -.5*width + (r+0)*width/nw2, height*(h+1)/float(nh) )
   f21 = facet((v21,v22,v23),color=(0,0,1),material=frictMat)
   f22 = facet((v21,v23,v24),color=(0,0,1),material=frictMat)
   v31 = Vector3( +.5*width - (r+0)*width/nw2, +.5*width, height*(h+0)/float(nh) )
   v32 = Vector3( +.5*width - (r+1)*width/nw2, +.5*width, height*(h+0)/float(nh) )
   v33 = Vector3( +.5*width - (r+1)*width/nw2, +.5*width, height*(h+1)/float(nh) )
   v34 = Vector3( +.5*width - (r+0)*width/nw2, +.5*width, height*(h+1)/float(nh) )
   f31 = facet((v31,v32,v33),color=(0,0,1),material=frictMat)
   f32 = facet((v31,v33,v34),color=(0,0,1),material=frictMat)
   v41 = Vector3( -.5*width, +.5*width - (r+0)*width/nw2, height*(h+0)/float(nh) )
   v42 = Vector3( -.5*width, +.5*width - (r+1)*width/nw2, height*(h+0)/float(nh) )
   v43 = Vector3( -.5*width, +.5*width - (r+1)*width/nw2, height*(h+1)/float(nh) )
   v44 = Vector3( -.5*width, +.5*width - (r+0)*width/nw2, height*(h+1)/float(nh) )
   f41 = facet((v41,v42,v43),color=(0,0,1),material=frictMat)
   f42 = facet((v41,v43,v44),color=(0,0,1),material=frictMat)
   facets.extend((f11,f12,f21,f22,f31,f32,f41,f42))
O.bodies.append(facets)
# walls top and bottom
#---------------------
walls = []
v1t=Vector3(0.6*width,-0.6*width,height)
v2t=Vector3(0.6*width,0.6*width,height)
v3t=Vector3(-0.6*width,0.6*width,height)
v4t=Vector3(-0.6*width,-0.6*width,height)
v1b=Vector3(0.6*width,-0.6*width,0.)
v2b=Vector3(0.6*width,0.6*width,0.)
v3b=Vector3(-0.6*width,0.6*width,0.)
v4b=Vector3(-0.6*width,-0.6*width,0.)
#-------------------
#-------------------
ftop1=facet((v2t,v4t,v1t),color=(1,0,0),material=frictMat,wire=False)
ftop2=facet((v2t,v3t,v4t),color=(1,0,0),material=frictMat,wire=False)
fbot1=facet((v1b,v2b,v4b),color=(0,1,0),material=frictMat,wire=False)
fbot2=facet((v2b,v3b,v4b),color=(0,1,0),material=frictMat,wire=False)
walls.extend((ftop1,ftop2,fbot1,fbot2))
O.bodies.append(walls)
#---------------------
mass = O.bodies[0].state.mass
for f in facets:
 f.state.mass = mass
 f.state.blockedDOFs = 'XYZz'
for w in walls:
 w.state.mass = mass
 w.state.blockedDOFs='XYZxy'
 print 'mass:',mass

# plots
plot.plots = { 'e':('s',), }
def plotAddData():
 f1 = sum(O.forces.f(b.id)[2] for b in top)
 f2 = sum(O.forces.f(b.id)[2] for b in bot)
 f = .5*(f2-f1)
 s = f/(pi*.25*width*width) if testType=='cyl' else f/(width*width) if testType=='cube' else None
 e = (top[0].state.displ()[2] - bot[0].state.displ()[2]) / (height-rParticle*2*bcCoeff)
 plot.addData(
  i = O.iter,
  s = -s,
  e = -e,
 )

# apply prestress to facets
def addForces():
 for f in facets:
  n = f.shape.normal
  a = f.shape.area
  O.forces.addF(f.id,preStress*a*n)

#apply velocity on walls (Plates)
def velWall():
 for w in walls:
  wID = w.id
  print 'gwID',wID
  if O.bodies[wID].state.pos[2] > 0.5*height:
   O.bodies[wID].state.blockedDOFs='z'
   O.bodies[wID].state.vel=(0.,0.,-velPlates)
  if O.bodies[wID].state.pos[2] < 0.5*height:
   O.bodies[wID].state.blockedDOFs='z'
   O.bodies[wID].state.vel=(0.,0.,velPlates)

# stop condition and exit of the simulation
def stopIfDamaged(maxEps=5e-2):
 extremum = max(abs(s) for s in plot.data['s'])
 s = abs(plot.data['s'][-1])
 e = abs(plot.data['e'][-1])
 if O.iter < 1000 or s > .5*extremum and e < maxEps:
  return
 f = os.path.join(exportDir,fileName)
 print 'gnuplot',plot.saveGnuplot(f,term='png')
 if runGnuplot:
  import subprocess
  os.chdir(exportDir)
  subprocess.Popen(['gnuplot',f+'.gnuplot']).wait()
 print 'Simulation finished'
 O.pause()
 #sys.exit(0) # results in some threading exception

O.dt=.1*utils.PWaveTimeStep()
enlargeFactor=1.5
O.engines=[
 ForceResetter(),
 InsertionSortCollider([
  Bo1_Sphere_Aabb(aabbEnlargeFactor=enlargeFactor,label='bo1s'),
  Bo1_Facet_Aabb()
 ]),
 InteractionLoop(
  [
   Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=enlargeFactor,label='ss2d3dg'),
   Ig2_Facet_Sphere_ScGeom(),
  ],
  [
   Ip2_FrictMat_FrictMat_FrictPhys(),
  ],
  [
   Law2_ScGeom_FrictPhys_CundallStrack(),
  ],
 ),
 PyRunner(iterPeriod=1,command="addForces()"),
 NewtonIntegrator(damping=.3),
 PyRunner(command='plotAddData()',iterPeriod=10),
 #PyRunner(iterPeriod=50,command='stopIfDamaged()'),
 #velWall(),
]

# run one step
velWall()
O.step()

# reset interaction detection enlargement
bo1s.aabbEnlargeFactor=ss2d3dg.interactionDetectionFactor=1.0

# initialize auto-updated plot
if runInGui:
 plot.plot()
 try:
  from yade import qt
  renderer=qt.Renderer()
  # uncomment following line to exagerate displacement
  #renderer.dispScale=(100,100,100)
 except:
  pass

# run
#O.run()