How to plot the total Force on a Clump

Hi everyone,

I am new to Yade and python and I am trying to simulate a pile and a shallow as a clump. One friend helps me to create a script to generate a sample of my soil and my foundation. I'm studying geotechinical problems. But i don't know how to obtain the total Force acting in the clump, because of a velocity imposed. Can someone tell me how to do this? This is my script:

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

from yade import pack,utils,plot,ymport
import math

################
### 1.INPUTS ###
################

# 1.1 IMPORT PARAMETERS FROM TABLE
nRead=readParamsFromTable(
 DensidadeP = 2.70,
 DensidadeS = 3.50,
 YoungP = 1200e6,
 PoissonP = 0.20,
 AnguloAtritoP = 39.5, #variar angulo de atrito para calibrar com a curva do professor Nakai
 Porosidade = 0.20,
 MatType = 'frict',
 TestDim = '2D',
 rate = 0.1,
 PSD = [[0.0032,0.006],[0.6,1.]],
 TensaoIsotropica = 1.9e4,
 verlet = 0.035,
 unknownOk=True
)
from yade.params import table

# 1.2 PARTICLE/MEDIUM PARAMETERS
mediumPorosity = table.Porosidade
particleDensity = table.DensidadeP #kg/m³
particleYoung = table.YoungP #Pa
particlePoisson = table.PoissonP #Adimensional
particleFricAng = radians(table.AnguloAtritoP)

# 1.3 ELEMENTS DIMENSIONS

 # 1.3.1 MEDIUM PARTICLE SIZE DISTRIBUTION
PSD = table.PSD

  # 1.3.2 SPHERES RADIUS
meanElementRadus = ((PSD[0][0]+PSD[0][-1])/2)/2
maxElementRadus=PSD[0][-1]/2
minElementRadus=PSD[0][0]/2

# 1.4 MATERIAL TYPE
materialType = 'frict' # 'frict', 'cohfrict', 'concrete'

# 1.5 DIMENSIONS (2D OR 3D)
testDimensions = table.TestDim

# 1.6 SHALLOW AND PILE (2D)
AlturaBloco=0.03 #Altura do bloco de fundação
BordasBloco=0.055 #Exedente dos bordos do bloco com relação à(s) estaca(s)
ProfEstaca=0.24 #Profundidade da estaca
CompPonta=0 #Comprimento da ponta da estaca (estacas com ponta ou tubulões)
DiametroEstaca=0.01 #Diâmetro da estaca
Afastamentos=5*DiametroEstaca #Região do solo nao influenciada pela estaca
XCentroE=0. #Coordenada 'x' do centro da estaca
YCentroE1=0.
YCentroE2=0.04
YCentroE3=-0.04 #Coordenada 'y' do centro da estaca
ZCentroE= (0.5-ProfEstaca/2) #Coordenada 'z' do centro da estaca
CentroEstaca1=(XCentroE,YCentroE1,ZCentroE) #Vetor de coordenadas do centro da estaca
TopoEstaca1=(XCentroE,YCentroE1,ZCentroE-ProfEstaca/2) #Vetor de coordenadas do topo da estaca
BaseEstaca1=(XCentroE,YCentroE1,ZCentroE+ProfEstaca/2) #Vetor de coordenadas da ponta da estaca
CentroEstaca2=(XCentroE,YCentroE2,ZCentroE) #Vetor de coordenadas do centro da estaca
TopoEstaca2=(XCentroE,YCentroE2,ZCentroE-ProfEstaca/2) #Vetor de coordenadas do topo da estaca
BaseEstaca2=(XCentroE,YCentroE2,ZCentroE+ProfEstaca/2) #Vetor de coordenadas da ponta da estaca
TopoEstaca3=(XCentroE,YCentroE3,ZCentroE-ProfEstaca/2) #Vetor de coordenadas do topo da estaca
BaseEstaca3=(XCentroE,YCentroE3,ZCentroE+ProfEstaca/2)
XCentroB=XCentroE #Coordenada 'x' do centro do bloco
YCentroB=YCentroE1 #Coordenada 'y' do centro do bloco
ZCentroB=0.5+AlturaBloco/2 #Coordenada 'z' do centro do bloco
CentroBloco=(XCentroB,YCentroB,ZCentroB) #Vetor de coordenadas do centro do bloco
Vertice1Bloco=(-1.01*1*0.003,YCentroB-DiametroEstaca-BordasBloco,ZCentroB-AlturaBloco/2)
Vertice2Bloco=(1.01*1*0.003,YCentroB+DiametroEstaca+BordasBloco,ZCentroB+AlturaBloco/2)

# 1.7 IMPORT/EXPORT NAMES
ImportName = 'Packing' + testDimensions + str(materialType)+ '-' + str(table.Porosidade)
ExportName = 'Triaxial' + testDimensions + str(materialType) + '-Phi' + str(table.AnguloAtritoP) + '-Y' + str(particleYoung/1.e6)+'MPa'+ '-Poisson' + str(particlePoisson) + '-' + str(mediumPorosity) + ' - Iso' + str(table.TensaoIsotropica)

# 1.8. USER DEFINED FUNCTIONS

 # 1.8.1. POROSITY

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

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

 # 2.1 MEDIUM MATERIAL

if materialType == 'frict':
 spheresMaterial=FrictMat(density=particleDensity, frictionAngle=particleFricAng, young=particleYoung, poisson=particlePoisson)

O.materials.append(spheresMaterial)

# 2.2 BOUNDARY MATERIAL
boundaryMaterial=FrictMat(density=particleDensity, frictionAngle=particleFricAng, young=10.*particleYoung, poisson=2*particlePoisson)
O.materials.append(boundaryMaterial)

#2.3 SHALLOW AND PILE MATERIAL
AngAtritoA=atan(0)
DensidadeA=27000.
YoungA=70.e9
PoissonA=.33
MatConcreto=O.materials.append(FrictMat(density=DensidadeA,young=YoungA,poisson=PoissonA,frictionAngle=AngAtritoA,label="aluminio"))

 ################## ### 3.GEOMETRY ### ##################

# 3.1 PACKING DIMENSIONS

width = 1
height = 0.5
dx = 1.01*2*maxElementRadus
dy = width
dz = height
CorteTransversal = dy
baseArea = dy
lateralAreaX = dy
lateralAreaY = 1.

# 3.2 BOUNDARIES

corner1,corner2 = (-.5*dx,-.5*dy,0.),(.5*dx,.5*dy,dz)
walls = aabbWalls([corner1, corner2], thickness=0., oversizeFactor=2, material=boundaryMaterial)
wallIds = O.bodies.append(walls)

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

#3.3 SHALLOW AND PILE
GeomEst1=pack.inCylinder(BaseEstaca1,TopoEstaca1,DiametroEstaca/2)
GeomEst2=pack.inCylinder(BaseEstaca2,TopoEstaca2,DiametroEstaca/2)
GeomEst3=pack.inCylinder(BaseEstaca3,TopoEstaca3,DiametroEstaca/2)
GeomBloco=pack.inAlignedBox(Vertice1Bloco,Vertice2Bloco)

GeomBloco=pack.regularOrtho(GeomBloco,radius=0.002,gap=0.,color=(0,0,1),material="aluminio")
GeomEst1=pack.regularOrtho(GeomEst1,radius=0.002,gap=0.,color=(0,0,1),material="aluminio")
GeomEst2=pack.regularOrtho(GeomEst2,radius=0.002,gap=0.,color=(0,0,1),material="aluminio")
GeomEst3=pack.regularOrtho(GeomEst3,radius=0.002,gap=0.,color=(0,0,1),material="aluminio")

O.bodies.appendClumped(GeomBloco)
O.bodies.appendClumped(GeomEst1)
O.bodies.appendClumped(GeomEst2)
O.bodies.appendClumped(GeomEst3)

# 3.4 IMPORT SPHERE PACKING

O.bodies.append(ymport.textExt(fileName = ImportName, format='x_y_z_r_matId'))

# 3.5. ASSIGN MATERIAL PROPERTIES TO THE SPHERES

for b in O.bodies:
 if isinstance(b.shape,Sphere):
  b.material=spheresMaterial

# 3.6. DEFINE BLOCKED DOFs
for b in O.bodies:
 if isinstance(b.shape,Sphere):
  b.state.blockedDOFs='xYZ'
 elif b.isClump:
  b.state.blockedDOFs='xXYZ'

# 4.1. TEST PROCEDURE

def Analise():
 O.forces.f(topo.id)[2]

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

if testDimensions == '2D':
 acelerador = 2
elif testDimensions == '3D':
 acelerador = 5

topo.state.pos[2] = topo.state.pos[2]-0.24

def PreencherTriaxial():
 if Porosity(CorteTransversal*topo.state.pos[2])>0.13:
  for b in O.bodies:
   if b.isClump:
    b.state.vel = (0,0,-10*acelerador)
  newton.damping = 0.7

 else:
  print 'fim'

#def force():
# global force
# if O.iter>0.:
# i=6.
# forca=0.
# while i>5 and i<352:
# forca= forca +O.forces.f(i)[2]
# i=i+1
# return list(force)

def myAddData():
        b=O.bodies[351]
        plot.addData(z1=0.53-b.state.pos[2], f1=O.forces.f(6)[2], i=O.iter, t=O.time)

plot.plots={'t':('f1','z1'),'f1':'z1'}
plot.plot()

O.run()
utils.waitIfBatch()

  ############################### ### 5. RECORD AND PLOT DATA ### ###############################

#def Displacement():
# global Displacement
# for b in b.isClump:
# print b.state.pos[2]

 ############################## ### 6.SIMULATION PROCEDURE ### ##############################

enlargeFactor = 1.5

O.engines=[
 ForceResetter(),
 InsertionSortCollider([
  Bo1_Sphere_Aabb(),
  Bo1_Box_Aabb(),
  Bo1_Wall_Aabb(),
  Bo1_Facet_Aabb()
 ],verletDist=table.verlet*2.*meanElementRadus),
 InteractionLoop(
  [Ig2_Wall_Sphere_ScGeom(),
  Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=enlargeFactor),
  Ig2_Box_Sphere_ScGeom(),
  Ig2_Facet_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_MindlinPhys(en=0.8, es=0.8)],
  [Law2_ScGeom_MindlinPhys_Mindlin()],
),

 GlobalStiffnessTimeStepper(
  active=True,
  defaultDt=SpherePWaveTimeStep(radius=meanElementRadus, density=O.materials[0].density, young=O.materials[0].young),
  timeStepUpdateInterval=1000,timestepSafetyCoefficient=0.5
 ),
# PyRunner(command='force()',iterPeriod=10,label='PrepararAmostra'),
 PyRunner(command='PreencherTriaxial()',iterPeriod=10,label='PrepararAmostra'),
 PyRunner(command='myAddData()',iterPeriod=10,label='PrepararAmostra'),
 NewtonIntegrator(damping=0.7, gravity=[0.,0.,-9.81], label='newton'),
]

 ########################## ### 7. RUN SIMULATION ### ##########################

O.run()
utils.waitIfBatch()