Creating a FCC structued packing and applying confining stress

Hey everyone!

I've been trying to create a cubic rock sample of FCC structure, and apply a 0.5 MPa confining stress to prepare the sample for other simulations. I tried following the periodic triaxial testing script but something is not correct in my setup as I'm getting an error message "body larger than half of the cell size encountered". I set up the FCC structure in the rock sample by specifying the centroid position and radius of every sphere particle, however I'm not sure if this is the correct way to create the packing. I have attached the my code below, if anyone could show me some pointers I would really appreciate it.

from yade import pack, qt, plot

r = 0.45
L = 2*sqrt(2)*r
sigmaIso = -5e5

# rock material
rockID = O.materials.append(CpmMat(young=15e6,poisson=.4,frictionAngle=radians(atan(0.5)),density=2650,label='spheres'))
#wall material
O.materials.append(FrictMat(young=15e6,poisson=.4,frictionAngle=0,density=0,label='frictionless'))
O.periodic = True

sphereArr = []

for m in range(0,46): #y axis
 for j in range(0,30): #z axis
  for i in range(0,15): #x axis
   if j%2==0:
    sphereArr.append((L*i,m*L/2,j*L/2))
   else:
    sphereArr.append((L*(i+(.5)),m*L/2,j*L/2))

O.bodies.append([sphere(center, r,material='spheres',color = (0.42,0.06,0.51) )for center in sphereArr])

## create walls around the packing
walls=aabbWalls([(0,0,0),(20,30,20)],thickness=1e-10,material='frictionless')
wallIds=O.bodies.append(walls)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 PeriTriaxController(label='triax',
 goal=(sigmaIso,sigmaIso,sigmaIso),stressMask=7,
 dynCell=True, maxStrainRate=(10,10,10),
 # wait until the unbalanced force goes below this value
 maxUnbalanced=.1,relStressTol=1e-3,
 # call this function when goal is reached and the packing is stable
 doneHook='compactionFinished()'
 ),
 NewtonIntegrator(damping=.2),
 PyRunner(command='GlobalData()',iterPeriod=500),
]

O.dt=.5*PWaveTimeStep()

def GlobalData():
 print 'porosity '+utils.porosity()
 print 'coord #' +utils.avgNumInteractions()

def compactionFinished():
 # set the current cell configuration to be the reference one
 O.cell.trsf=Matrix3.Identity
 # change control type: keep constant confinement in x,y, 20% compression in z
 triax.goal=(sigmaIso,sigmaIso,-.2)
 triax.stressMask=3
 # allow faster deformation along x,y to better maintain stresses
 triax.maxStrainRate=(1.,1.,.1)
 # next time, call triaxFinished instead of compactionFinished
 triax.doneHook='triaxFinished()'
 # do not wait for stabilization before calling triaxFinished
 triax.maxUnbalanced=10

def triaxFinished():
 print 'Finished'
 O.pause()