Shear force histogram

Well, could you please post a minimum working example? Doing so we can simply run it and have a look.
Thanks

Thanks Alexander,
Actually I'm using triaxial example but with some differences: 1- My loading is strain controlled, 2- I'm applying it with the same rate at all three directions ( Confining pressure) to a maximum stress value of 4 Mpa. Unfortunately I don't have access to my script at the moment but the above script should work with any given script because I put it in a function and calling it by pyrunner to produce shear force histogram plots. I just want to know whether the formulation is correct or not especially when it comes to the calculation of contact vector in tangential direction. thanks in advance for your help. :)

here is my full script:

# -*- coding: utf-8 -*-
#*************************************************************************
# Copyright (C) 2010 by Bruno Chareyre *
# bruno.chareyre_at_grenoble-inp.fr *
# *
# This program is free software; it is licensed under the terms of the *
# GNU General Public License v2 or later. See file LICENSE for details. *
#*************************************************************************

## This script details the simulation of a triaxial test on sphere packings using Yade
## See the associated pdf file for detailed exercises
## the algorithms presented here have been used in published papers, namely:
## * Chareyre et al. 2002 (http://www.geosyntheticssociety.org/Resources/Archive/GI/src/V9I2/GI-V9-N2-Paper1.pdf)
## * Chareyre and Villard 2005 (https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
## * Scholtès et al. 2009 (http://dx.doi.org/10.1016/j.ijengsci.2008.07.002)
## * Tong et al.2012 (http://dx.doi.org/10.2516/ogst/2012032)
##
## Most of the ideas were actually developped during my PhD.
## If you want to know more on micro-macro relations evaluated by triaxial simulations
## AND if you can read some french, it is here: http://tel.archives-ouvertes.fr/docs/00/48/68/07/PDF/Thesis.pdf

from yade import pack,plot,qt
import matplotlib; matplotlib.rc('axes',grid=True)
import pylab
############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################
key="K"
num_spheres=8000
psdSizes,psdCumm=[0.26*0.866,0.388*0.866,0.536*0.866,0.706*0.866,0.976*0.866,1.333*0.866,1.757*0.866,2.458*0.866,2.771*0.866,3.124*0.866,3.608*0.866,4.019*0.866,4.424*0.866,4.811*0.866,5.556*0.866,6.35*0.866],[0,1.865,3.657,6.12,9.25,14.4,18.88,29.18,35.45,42.4,52.2,61.857,71.716,80.898,90.522,100]
#targetPorosity = 0.387 #the porosity we want for the packing
compFricDegree = 26.5 # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
finalFricDegree = 26.5 # contact friction during the deviatoric loading
targetPorosity=0.3 ### INJA BARAYE MAX ACN va MIN e
rate=0.0002 # loading rate (strain rate)
damp=0.25 # damping coefficient!!!!!!!!!!
stabilityThreshold=0.01 # we test unbalancedForce against this value in different loops (see below)
young=520e6 # contact stiffness!!!CHANGED!!!
mn,mx=Vector3(0,0,0),Vector3(44.25,44.25,44.25) # corners of the initial packing

## create materials for spheres and plates
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(compFricDegree),density=2000,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))

# create walls around the packing
walls=aabbWalls([mn,mx],material='walls',oversizeFactor=1)
wallIds=O.bodies.append(walls)

## use a SpherePack object to generate a random loose particles packing
sp=pack.SpherePack()
sp.particleSD2(radii=psdSizes,passing=psdCumm,numSph=8000,cloudPorosity=0.57)
O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])
#walls=aabbWalls(material='walls',oversizeFactor=1)
#wallIds=O.bodies.append(walls)
#or alternatively (higher level function doing exactly the same):
#sp.toSimulation(material='spheres')

############################
### DEFINING ENGINES ###
############################

triax=TriaxialStressController(
 ## ThreeDTriaxialEngine will be used to control stress and strain. It controls particles size and plates positions.
 ## this control of boundary conditions was used for instance in http://dx.doi.org/10.1016/j.ijengsci.2008.07.002
 maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)!!!!!!
 finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)!!!!!!!!!!
 thickness = 0,
 ## switch stress/strain control using a bitmask. What is a bitmask, huh?!
 ## Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0.
 ## Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4
 ## to put it differently, the mask is the integer whose binary representation is xyz, i.e.
 ## "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x and y and z", etc.
 stressMask = 7,
 internalCompaction=True, # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_MindlinPhys()],
  [Law2_ScGeom_MindlinPhys_Mindlin()]
 ),
 ## We will use the global stiffness of each body to determine an optimal timestep (see https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.7),
 triax,
 TriaxialStateRecorder(iterPeriod=100,file='WallStresses.dat'),
 qt.SnapshotEngine(fileBase="snap",iterPeriod=0,label='snapshoter'),
 newton,

]

########################################
#### APPLYING CONFINING PRESSURE ###
########################################

#the value of (isotropic) confining stress defines the target stress to be applied in all three directions
triax.goal1=triax.goal2=triax.goal3=100000
while 1:
  O.run(1000, True)
  #the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
  unb=unbalancedForce()
  print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
  if unb<stabilityThreshold and abs(100000-triax.meanStress)/100000<0.01:
    break

O.save('confinedState'+'.yade.gz')
print "### Isotropic state saved ###"
print 'ACN=',utils.avgNumInteractions(),'Porosity=',utils.voxelPorosityTriaxial(triax),'Calculation Time(Sec)=',O.realtime

##############################
### DEVIATORIC LOADING ###
##############################
#We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
triax.internalCompaction=False

# Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction(radians(finalFricDegree))

triax.wall_bottom_activated=True
triax.wall_top_activated=True
triax.wall_left_activated=True
triax.wall_right_activated=True
triax.wall_back_activated=True
triax.wall_front_activated=True

#set stress control on x and z, we will impose strain rate on y
triax.stressMask = 0
#now goal2 is the target strain rate
triax.goal2=-rate
# we define three lateral stresses during the test, here the same 10kPa as for the initial confinement.
triax.goal1=-rate
triax.goal3=-rate

#we can change damping here. What is the effect in your opinion?
#newton.damping=0.1

#Save temporary state in live memory. This state will be reloaded from the interface with the "reload" button.
O.saveTmp()
#####################################################
### Example of how to record and plot data ###
#####################################################

from yade import plot
import pylab
import numpy as np
import os

#########################################################################
#########################################################################
#########################################################################

class StressChecker():
 dStress=nextStress=100000
 def Sintrhisto(self):
  stress=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
  axis=2
  ax1,ax2=(axis+0)%3,(axis+1)%3
  angles,forces=[],[]
  for z in O.interactions:
   if not z.isReal: continue
   if z.id1<6 or z.id2<6: continue
   norm=z.geom.normal
   if norm[ax1]==0:
    angle=0
    force=z.phys.shearForce.norm()
   else:
    angle=atan(norm[ax2]/norm[ax1])
    force=z.phys.shearForce.norm()
   angles.append(angle+pi/2)
   forces.append(force*10e-6)
  pylab.figure()
  values,bins=numpy.histogram(angles,weights=forces,bins=30)
  subp=pylab.subplot(111,polar=True)
  pylab.bar(left=bins[:-1],height=values,width=np.pi/(1.05*30),alpha=.7,label=['xy'])
  pylab.xlabel('Shear Force Histogram-XY Plane')
  pylab.plot()
  pylab.savefig("interaction histogram-shear")
  FName = "interaction histogram-shear.png"
  NFName = "interaction histogram-shear-%.2f kPa.png"%float(stress/1000)
  os.rename(FName,NFName)

 def output(self):
  stress=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
  if abs(stress) > self.nextStress:
   self.nextStress += self.dStress
   self.Sintrhisto()

checker=StressChecker()
# include a periodic engine calling that function in the simulation loop
O.engines=O.engines[0:6]+[PyRunner(iterPeriod=20,command='checker.output()')]+O.engines[6:8]
#plot.plot()
#O.run(100,True)
#O.run(1000000)
#### PLAY THE SIMULATION HERE WITH "PLAY" BUTTON OR WITH THE COMMAND O.run(N) #####

I figured out that if you find the average of forces in each direction, you can draw polar histogram of forces easily.