determine the force applied on the particle

Hello everyone,

I want to calculate the potential energy of the beam and the work performed on the beam in 3-points bending test.
To calculate the energy potential in each particle i use this formula:
'Ub(particle) = m(mass of particle)*g*d(displacement of particle)
and the work performed in each particle on the boundary of the model, i use this formula:
'W(particle)=F(force exerced on the particle)*d(displacement of particle)

My question is how can i determine the force exerced on each particle and the displacement of each particle?

Thanks so much
Jabrane

from yade import ymport, utils, pack, export
from yade import plot
from pylab import *
import math

#---------------- SIMULATIONS DEFINED HERE (assembly, material, boundary conditions)

#### packing (previously constructed)
PACKING='X240Y80Z30_400'
OUT=PACKING+'_flexion_3_points'

#### Simulation Control
DAMP=0.4 # numerical damping
saveData=1000 # data record interval
iterMax=10000 # maximum number of iteration (to be adjusted)
saveVTK=1000 # Vtk files record interval

#### Material microproperties -> Lac du Bonnet granite (cf. A DEM model for soft and hard rock, Scholtes & Donze, JMPS 2013)
intR=1.6684# allows near neighbour interaction and coordination number K=13 (determined with coordinationNumber.py -> to be adjusted for every packing)
DENS=4000 # could be adapted to match material density: dens_DEM=dens_rock*(V_rock/V_particles)=dens_rock*1/(1-poro_DEM) -> poro?
YOUNG=65e9
FRICT=10
ALPHA=0.4
TENS=8e6
COH=160e6
#### material definition
sphereMat = JCFpmMat(type=1,density=DENS,young=YOUNG,poisson = ALPHA,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH)
wallMat = JCFpmMat(type=0,density=DENS,young=YOUNG,frictionAngle=radians(0))

for mat in (sphereMat,wallMat):
   O.materials.append(mat) # then wallMat will be used if material is not specified

#### preprocessing to get dimensions
O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0)))

dim=utils.aabbExtrema()
xinf=dim[0][0]
xsup=dim[1][0]
X=xsup-xinf
yinf=dim[0][1]
ysup=dim[1][1]
Y=ysup-yinf
zinf=dim[0][2]
zsup=dim[1][2]
Z=zsup-zinf

r=X/15.
#### preprocessing to get dimensions
R=0
Rmax=0
numSpheres=0.
for o in O.bodies:
 if isinstance(o.shape,Sphere):
   numSpheres+=1
   R+=o.shape.radius
   if o.shape.radius>Rmax:
     Rmax=o.shape.radius
Rmean=R/numSpheres

O.reset() # all previous lines were for getting dimensions of the packing to create walls at the right positions (below) because walls have to be genrated after spheres for FlowEngine

O.bodies.append(geom.facetCylinder(center=(xinf+X/12.,yinf-0.9*r,Z/2.),radius=r,height=Z,orientation=Quaternion((1, 0, 0), 0),segmentsNumber=20,wire=False,material=wallMat)) # bas gauche
O.bodies.append(geom.facetCylinder(center=(xinf+11*X/12.,yinf-0.9*r,Z/2.),radius=r,height=Z,orientation=Quaternion((1, 0, 0), 0),segmentsNumber=20,wire=False,material=wallMat)) # bas droite
piston=O.bodies.append(geom.facetCylinder(center=(0.5*X,Y+0.9*r,Z/2.),radius=r,height=Z,dynamic=False,orientation=Quaternion((1, 0, 0), 0),segmentsNumber=20,wire=False,material=wallMat)) # haut

p0=O.bodies[piston[0]].state.pos[1]

### packing
beam=O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat))

# Definition of the facets for joint's geometry
alpha = 90 #angle of the crack with horizontal (pi/4.)
c = 0.375*Y # length of the crack

# 4 points pour l'entaille
ptA = gts.Vertex( X/2., c, 8./7.*Z)
ptB = gts.Vertex( X/2., 0, 8./7.*Z)
ptApr = gts.Vertex(X/2., c, -1./7.*Z)
ptBpr = gts.Vertex(X/2., 0, -1./7.*Z)

e1 = gts.Edge(ptA,ptB)
e2 = gts.Edge(ptA,ptApr)
e3 = gts.Edge(ptApr,ptB)
f1 = gts.Face(e1,e2,e3)

e4 = gts.Edge(ptB,ptBpr)
e5 = gts.Edge(ptBpr,ptApr)
f2 = gts.Face(e4,e5,e3)

s1 = gts.Surface()
s1.add(f1)
s1.add(f2)
facet = gtsSurface2Facets(s1,wire = False,material=wallMat)
O.bodies.append(facet)

## Determinate the particle in the boundary
limite=[]
limite=range(0, len(beam))

for i in range(0, len(beam)):
 if X/3. <= O.bodies[beam[i]].state.pos[0] <= 2*X/3. or Y/3. <= O.bodies[beam[i]].state.pos[1] <= 2*Y/3. or Z/3. <= O.bodies[beam[i]].state.pos[2] <= 2*Z/3. :
   limite[i] = 0
 else:
   limite[i] = beam[i]

for i in range(0, len(limite)):
 try:
      limite.remove(0)
 except ValueError:
      pass

## Determinate the initial position of the particle on the boundary

poslimx=[]
poslimy=[]
poslimz=[]
poslimx=range(0,len(limite))
poslimy=range(0,len(limite))
poslimz=range(0,len(limite))

for i in range(0,len(limite)):
 poslimx[i]=O.bodies[limite[i]].state.pos[0]

for i in range(0,len(limite)):
 poslimy[i]=O.bodies[limite[i]].state.pos[1]

for i in range(0,len(limite)):
 poslimz[i]=O.bodies[limite[i]].state.pos[2]

### set a color to the spheres
for o in O.bodies:
 if isinstance(o.shape,Sphere):
   o.shape.color=(0.7,0.5,0.3)

#---------------- ENGINES DEFINED HERE

#### simulation is defined here (DEM loop, interaction law, servo control, recording, etc...)
##### simulation piston's motion

for i in range(0,len(piston)):
 O.bodies[piston[i]].state.vel[1]=-1

##### simulation beam's grains motion
O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb'),Bo1_Facet_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom'),Ig2_Facet_Sphere_ScGeom()],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT,label='interactionLaw')]
 ),

        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.4, defaultDt=0.1*utils.PWaveTimeStep()),
        NewtonIntegrator(damping=DAMP,label="newton"),
        PyRunner(iterPeriod=int(saveData),initRun=True,command='recorder()',label='data'),
        VTKRecorder(iterPeriod=int(saveVTK),initRun=True,fileName=OUT+'-',recorders=['spheres','bstresses','jcfpm','cracks'],Key=OUT,label='vtk')
]

## Determinate the initial position of all particles of the model
posx=[]
posy=[]
posz=[]
posx=range(0,len(beam))
posy=range(0,len(beam))
posz=range(0,len(beam))

for i in range(0,len(beam)):
 posx[i]=O.bodies[beam[i]].state.pos[0]

for i in range(0,len(beam)):
 posy[i]=O.bodies[beam[i]].state.pos[1]

for i in range(0,len(beam)):
 posz[i]=O.bodies[beam[i]].state.pos[2]

## Calculate the work on the boundary of model
W=[]
W=range(0,len(limite))
## W est le travail effectue sur les billes du limite de modele
#for i in range(0, len(W)):
# W[i]=utils.sumFacetNormalForces(ids=limite[i],axis=0)*abs(O.bodies[limite[i]].state.pos[0]-poslimx[i])+utils.sumFacetNormalForces(ids=limite[i],axis=1)*abs(O.bodies[limite[i]].state.pos[1]-poslimy[i])+utils.sumFacetNormalForces(ids=limite[i],axis=2)*abs(O.bodies[limite[i]].state.pos[2]-poslimz[i])

## gravitational acceleration
g=[0,0,-9.8]

## calculate the potential energy of model
Ub=[]
Ub=range(0,len(beam))
## Ub is the potential energy of the model

for i in range(0, len(Ub)):
 Ub[i]=O.bodies[beam[i]].state.mass*abs(g[2])*abs(O.bodies[beam[i]].state.pos[2]-posz[i])

# if you want to plot during simulation
plot.plots={'i':('f')}
#plot.plot()

#---------------- SIMULATION STARTS HERE

#### manage interaction detection factor during the first timestep and then set default interaction range ((cf. A DEM model for soft and hard rock, Scholtes & Donze, JMPS 2013))
O.step();
### initializes the interaction detection factor
SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.

O.run(iterMax)