Direct shear test simulation

Hello all,

I am trying to set-up a direct shear test box simulation to retrieve the soil resistance parameters of a soil. I am using the PeriTriaxController, with uniaxial compression first and then shearing. I want to force the failure in the middle of the box, so I am blocking the displacements and rotations in the lower half and applying a velocity in one of the directions, but it seems it is not working. I think it is in part because of the periodic boundary conditions. Is there a way to apply a periodic boundary conditions to only one of the boundaries, let us say the top one? I have looked for something similar in the mailing list and the documentation and I have not found anything. If my approach is wrong please let me know if I should try something else. Here is my script. Thanks in advance.

Cheers,

Mario

#coding: utf-8
#Direct shear test simulation
#
from yade import pack,qt,plot

#set periodic boundaries
O.periodic=True
#
O.cell.setBox(.1,.1,.1)
#O.cell.Hsize=Matrix3(0.1,0,0, 0,0.1,0, 0,0,0.1)
sp=pack.SpherePack()
radius=5e-3
num=sp.makeCloud((0,0,0),(.1,.1,.1),radius,.2,500,periodic=True) # min,max,radius,rRelFuzz,spheresInCell,periodic
O.bodies.append([sphere(s[0],s[1]) for s in sp])

O.engines=[
#reset forces
        ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb()],verletDist=.05*radius),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
                ),
        #stressMask, 1=stress, 0=strain. For x,y,z in binary format, 110 means stress in x and y and strain in z.
        #In python the format is 0b110, can be converted to decimal format.
 PeriTriaxController(dynCell=True,mass=0.2,maxUnbalanced=0.01,relStressTol=0.02,goal=[0,-0.25e4,0],
                            stressMask=0b010,globUpdate=5,maxStrainRate=[1.,1.,1.],doneHook='shear()',label='dst'),
 NewtonIntegrator(damping=.2),
        # record data for plotting every 100 steps; addData function is defined below
        PyRunner(command='addPlotData()',iterPeriod=500),
        # saves data to a txt file, calls the function 'savenormaltxt', defined below
        PyRunner(command='savenormaltxt()',iterPeriod=500)
          ]

phase=0
def shear():
 global phase
 if phase==0:
  print 'Here we are: stress',dst.stress,'strain',dst.strain,'stiffness',dst.stiff
  print 'Now shearing.'
  #O.cell.velGrad=Matrix3(6,0,0, 0,0,0, 0,0,0)
  #dst.stressMask=0b010
                goalx=0.5
  dst.goal=[goalx,-0.25e4,0]
                dst.doneHook='endSim()'
                for b in O.bodies:
                    b.state.vel=(10,0,0)
                    # block X,Y,Z rotations and translations for the lower half of the box
                    if b.state.pos[1]<0.05:
                        b.dynamic=False
  phase+=1
 #elif phase==1:
  #print 'Here we are: stress',triax.stress,'strain',triax.strain,'stiffness',triax.stiff
  #phase+=1
  ##print 'Done, pausing now.'
  ##O.pause()

def endSim():
    print 'End of simulation'
    print 'stress',dst.stress
    print 'strain',dst.strain
    O.pause()

O.dt=PWaveTimeStep()
qt.View()

#uses the function addData to assign labels and plot the values of the stresses of the stress tensor.
def addPlotData():
    plot.addData(sx=dst.stress[0],sy=dst.stress[1],sz=dst.stress[2],ex=dst.strain[0],ey=dst.strain[1],ez=dst.strain[2],
                        i=O.iter,unbalanced=utils.unbalancedForce())
def savenormaltxt():
    plot.saveDataTxt('dst-normal.txt')

#defines the values of x and y for every plot
plot.plots={'ey':('sy')}
#shows the plot
plot.plot()

O.run()