rigid boundaries

Hi everyone:

I am new to Yade and python and I am trying to simulate a pile as a rigid boundaries. One friend helps me to create a script to generate a sample of my soil. I'm studying geotechinical problems. But i don't know how to create a new rigid boundary to simulate a pile. Can someone tell me how to do this? This is my script:

# -*- coding: utf-8 -*-

from yade import pack,utils,plot,export

import numpy as np
import operator,fractions

# TODAS AS UNIDADES ESTÃO NO S.I.

###############################
### 1.PARÂMETROS DE ENTRADA ###
###############################

# 1.1 IMPORT PARAMETERS FROM TABLE
nRead=readParamsFromTable(
 DensidadeP = 2.7,
 PoissonP = 0.2,
 AnguloAtritoP = 0,
 CoesaoP = 0,
 PorosidadeS = 0.20,
 MatType = 'frict',
 TestType = '2D',
 PSD = [[0.0032,0.006],[0.6,1.]],
 verlet = 0.35,
 YoungP = 7.e4,
 unknownOk=True
)
from yade.params import table

# 1.2 PARAMETROS DO MEIO (MACROSCÓPICOS)
DensidadeS = table.DensidadeP #tf/m³
YoungP = table.YoungP #kN/m²
PoissonP = table.PoissonP #Adimensional
AnguloAtritoP = radians(table.AnguloAtritoP)
CoesaoP = table.CoesaoP
PorosidadeS = table.PorosidadeS

# 1.3 DIMENSÕES DOS ELEMENTOS

# 1.3.1 CURVA GRANULOMETRICA DO SOLO - PSD=[[diâmetros],[Quantidade passante Acumulada]]

PSD = table.PSD

# 1.3.2 DADOS MÉDIOS DOS ELEMENTOS
rMeanElemento = ((PSD[0][0]+PSD[0][-1]))/2
rRelFuzzElemento = PSD[0][-1]/(2*rMeanElemento)-1
rMaxElemento = PSD[0][-1]/2
rMinElemento = PSD[0][0]/2

# 1.4 TIPO DE MATERIAL
MatType = table.MatType # 'frict', 'cohfrict', 'concrete'

# 1.5 TIPO DE AMOSTRA
TestType = table.TestType

##################
### 2.MATERIAL ###
##################

# 2.1 Material Meio
if MatType == 'frict':
 Material=FrictMat(density=20400, frictionAngle=0)
elif MatType == 'cohfrict':
 Material=CohFrictMat(density=2600)
elif MatType == 'concrete':
# Material=CpmMat(density=2600)
 Material=FrictMat(density=2600, frictionAngle=0)
O.materials.append(Material)

# 2.2 Material Contorno
MatContorno = FrictMat(frictionAngle=0)
O.materials.append(MatContorno)

###################
### 3.GEOMETRIA ###
###################

# 3.1 Triaxial

# 3.1.1 Dimensões
width = 1
height = 0.5

if TestType == '2D':
 dx = 1.01*2*rMaxElemento
 dy = width
 dz = height
 CorteTransversal = dy
elif TestType == '3D':
 dx = width
 dy = width
 dz = height
 CorteTransversal = dy*dx

# 3.1.2 Faces
mn,mx=(-.5*dx,-.5*dy,0),(.5*dx,.5*dy,dz)
walls=aabbWalls([mn,mx],thickness=0,material=MatContorno)
wallIds=O.bodies.append(walls)

base=O.bodies[4]
topo=O.bodies[5]

##################
### 4.CÁLCULOS ###
##################

# 4.1. InteractionLoop
enlargeFactor=1.5

LoopFrict=InteractionLoop(
 [
  Ig2_Wall_Sphere_ScGeom(),
  Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=enlargeFactor),
  Ig2_Box_Sphere_ScGeom(),
  Ig2_Facet_Sphere_ScGeom(),
 ],
 [
  Ip2_FrictMat_FrictMat_FrictPhys(),
 ],
 [
  Law2_ScGeom_FrictPhys_CundallStrack(),
 ],
)

LoopCohFrict=InteractionLoop(
 [
  Ig2_Wall_Sphere_ScGeom(),
  Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=enlargeFactor),
  Ig2_Box_Sphere_ScGeom6D(),
  Ig2_Facet_Sphere_ScGeom6D(),
 ],
 [
  Ip2_FrictMat_FrictMat_FrictPhys(),
  Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=True,label='ipf')
 ],
 [
  Law2_ScGeom_FrictPhys_CundallStrack(),
  Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True)
 ],
)
LoopConcrete=InteractionLoop(
 [
  Ig2_Wall_Sphere_ScGeom(),
  Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=enlargeFactor),
  Ig2_Box_Sphere_ScGeom(),
  Ig2_Facet_Sphere_ScGeom(),
 ],
 [
  Ip2_FrictMat_FrictMat_FrictPhys(),
  Ip2_CpmMat_CpmMat_CpmPhys(cohesiveThresholdIter=10),
  Ip2_FrictMat_CpmMat_FrictPhys(),
 ],
 [
  Law2_ScGeom_FrictPhys_CundallStrack(),
  Law2_ScGeom_CpmPhys_Cpm(),
 ],
)

if MatType == 'frict':
 IntLoop = LoopFrict
elif MatType == 'cohfrict':
 IntLoop = LoopCohFrict
elif MatType == 'concrete':
 IntLoop = LoopConcrete

# 4.2 Factory Parameters

MassFlowRate = 5000/(CorteTransversal*dz)
VMAX = 5
VMIN = 5
FCenter = (0.,0.,height/2)
FExtents2D = (0,width/2.,height/2)
FExtents3D = (width/2.,width/2.,height/2)

factoy2D = BoxFactory(
  center = FCenter, extents = FExtents2D,
  PSDsizes = PSD[0], PSDcum = PSD[1], PSDcalculateMass = True,
  vMax = VMAX, vMin = VMIN, vAngle = 0., normal = (0.,0.,-1.),
  maxParticles = -1,
  massFlowRate = MassFlowRate,
# exactDiam = False,
  label = 'factory',
  blockedDOFs = 'xYZ',
  materialId = 0
)

factoy3D = BoxFactory(
  center = FCenter, extents = FExtents3D,
  PSDsizes = PSD[0], PSDcum = PSD[1], PSDcalculateMass = True,
  vMax = 5*VMAX, vMin = 5*VMIN, vAngle = 0., normal = (0.,0.,-1.),
  maxParticles = -1,
  massFlowRate = MassFlowRate,
# exactDiam = False,
  label = 'factory',
  materialId = 0
)

if TestType == '2D':
 factory = factoy2D
elif TestType == '3D':
 factory = factoy3D

#4.3 Engines
O.engines=[
 ForceResetter(),
 InsertionSortCollider([
  Bo1_Sphere_Aabb(aabbEnlargeFactor=enlargeFactor),
  Bo1_Box_Aabb(),
  Bo1_Wall_Aabb(),
  Bo1_Facet_Aabb()
 ], verletDist=0.07*2*rMeanElemento),
 IntLoop,
 factory,
 DomainLimiter(lo=(-dx/2,-dy/2.,0.),hi=(dx/2,dy/2.,dz+5*rMeanElemento),iterPeriod=100),
 GlobalStiffnessTimeStepper(
  active=True,
  defaultDt=SpherePWaveTimeStep(radius=rMeanElemento, density=O.materials[0].density, young=O.materials[0].young),
  timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8
 ),
 NewtonIntegrator(damping=0.3, gravity=[0.,0.,-9.81], label='newton'),
 PyRunner(command='PreencherTriaxial()',iterPeriod=10,label='PrepararAmostra'),
# PyRunner(command='addPlotData()',iterPeriod=10),
]

##################
### 5.PyRunner ###
##################

def Porosity(Region):
 if TestType == '3D':
  return porosity(Region)
 elif TestType == '2D':
  soma=0.
  for b in O.bodies:
   if isinstance(b.shape,Sphere):
    soma+=pi*(b.shape.radius**2)
  return (Region-soma)/Region

if TestType == '2D':
 acelerador = 1
elif TestType == '3D':
 acelerador = 5

Compressao=False
Descompressao=False
AlturaPreenchida=0

def PreencherTriaxial():
 global Compressao
 global Descompressao
 global topo
 global AlturaPreenchida
 if Compressao:
  if Porosity(CorteTransversal*topo.state.pos[2])<PorosidadeS:
   topo.state.vel = (0.,0.,1.*acelerador)
   Descompressao = True
   Compressao = False
   newton.damping = 0.7
  else:
   topo.state.vel = (0,0,-10*acelerador)
 elif Descompressao:
  if abs(O.forces.f(topo.id)[2])==0:
   if topo.state.vel==(0,0,1*acelerador):
    print Porosity(CorteTransversal*topo.state.pos[2])
    AlturaPreenchida=topo.state.pos[2]
    print AlturaPreenchida/height
   if topo.state.pos==(0,0,dz+rMaxElemento):
    Material.density=DensidadeS*100/(1-Porosity(CorteTransversal*AlturaPreenchida))
    if AlturaPreenchida>height-rMeanElemento:
     topo.state.vel=(0,0,0)
     calm()
     Descompressao=False
     Material.frictionAngle=radians(35)
     newton.damping=0.4
     PrepararAmostra.command='EstabilizarTriaxial()'
     PrepararAmostra.iterPeriod=1
    else:
     topo.state.vel = (0,0,0)
     Descompressao = False
     factory.massFlowRate = 220/(CorteTransversal*(height-AlturaPreenchida))
# factory.center[2] = (AlturaPreenchida+height)/2
# factory.extents[2] = (height-AlturaPreenchida)/2
     newton.damping = 0.1
# if AlturaPreenchida > 0.95*height:
# newton.gravity = [0.,0.,-9.81]
   else:
# topo.state.vel=(0,0,30*acelerador)
    topo.state.vel=(0,0,0)
    topo.state.pos=(0,0,dz+rMaxElemento)
 elif factory.massFlowRate==0:
  Compressao=True

def EstabilizarTriaxial():
 if topo.state.pos[2]>0.98*height:
  topo.state.vel=(0,0,-25)
 else:
  topo.state.vel=(0,0,0)
  if utils.unbalancedForce()<0.005: #Checar isso para concreto
   PrepararAmostra.command='FinalizarPreparo()'

def FinalizarPreparo():
 print "Amostra Preparada"
 print "Densidade"
 print Material.density
 print "Porosidade"
 print Porosity(CorteTransversal*AlturaPreenchida)
 print "Tempo"
 print O.realtime
 export.textExt('Packing' + TestType + MatType + '-' + str(PorosidadeS),'x_y_z_r_matId')
 O.pause()

O.run()
utils.waitIfBatch()

Thank you so much and wish have a good day!

Regards,
Jessica.