YADE Tutorials

Tutorials

YADE uses python scripts as inputs for all simulations. This allows considerable flexibility as helper functions (for particle generation, analysis, plotting, etc.) can be included in each script.

Installation

YADE is supported on Debian and Ubuntu systems. On these supported systems, YADE can typically be installed as a package with:

sudo apt-get install yade

Full installation instructions for installing the latest daily version or building from source can be found at the YADE installation documentation.

Tutorial 1 - Bouncing Sphere

A free sphere bounces on top of a fixed sphere under gravity.

# Script to simulate a free sphere bouncing under gravity on a fixed sphere
# 2025 © Vasileios Angelidakis <v.angelidakis@qub.ac.uk>
# Developed based on https://yade-dem.org/doc/tutorial-examples.html with modifications to adopt the viscoelastic Hertz-Mindlin contact law

# Create results directory
import os
import errno
try:
	os.mkdir('./results/')
except OSError as exc:
	if exc.errno != errno.EEXIST:
		raise
	pass

# Create materials
O.materials.append(FrictMat(young=5e9, poisson=.25, frictionAngle=0.0, density=2500, label='particles'))

# Create particles
O.bodies.append([
#	box(center=(0,0,0), extents=(1,1,0.25), material='particles', fixed=True),
	sphere(center=(0, 0, 0), radius=.5, material='particles', fixed=True), # fixed sphere
	sphere(center=(0, 0, 2), radius=.5, material='particles') # free sphere
])

# Engines - Simulation loop
O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(hertzian=True)],  # collision geometry
                [Ip2_FrictMat_FrictMat_MindlinPhys(en=1.0,label='ip2')],  # collision "physics"
                [Law2_ScGeom_MindlinPhys_Mindlin()]  # contact law -- apply forces
        ),
        # Apply gravity force to particles. damping: numerical dissipation of energy.
        NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.0) # No numerical damping
]

# Timestep
O.dt = 0.0005 * RayleighWaveTimeStep() # set timestep to a fraction of the Rayleigh Wave timestep

# save the simulation, so that it can be reloaded later, for experimentation
O.saveTmp()

Gl1_Sphere.quality=5 # Make the visualisation of spheres nicer


# Live plotting of results
showPlots = True

if showPlots:
	from yade import plot

	# Function to plot data periodically
	def addPlotData():
		if O.interactions.countReal()>0: # If the system has "real" interactions, then plot force vs displacement
			i=O.interactions[0,1]
			a=i.geom.penetrationDepth	# Normal deformation
			Fn=i.phys.normalForce[2] 	# Normal force -> Z component
		else:
			a=0
			Fn=0
		plot.addData(deformation=a, force=Fn, position=O.bodies[1].state.pos[2], iteration=O.iter)
#		plot.saveDataTxt('results/bouncing-sphere.txt')

	# besides unbalanced force evolution, also plot the displacement-force diagram
	plot.plots = {'deformation': ('force'), 'iteration': ('position')}
	plot.plot()

	O.engines = O.engines + [PyRunner(command='addPlotData()', iterPeriod=200)]

Tutorial 2 - Angle of Repose (AoR)

A granular packing is generated within a cylinder made of facets. It is let settle under gravity, and then the top part of the cylinder is lifted away to make the material form a repose state.

# Script to perform an Angle of Repose simulation
# 2025 © Vasileios Angelidakis <v.angelidakis@qub.ac.uk>

from yade import qt, plot

# Create results directories
import os
import errno
try:
	os.mkdir('./results/')
	os.mkdir('./results/angle-of-repose-vtk/')
except OSError as exc:
	if exc.errno != errno.EEXIST:
		raise
	pass

# Define materials
O.materials.append(FrictMat(young=5e5, poisson=.25, frictionAngle=atan(0.3), density=2500, label='walls'))
O.materials.append(FrictMat(young=5e5, poisson=.25, frictionAngle=atan(0.4), density=2500, label='spheres'))

# Generate facets
radius=0.2
height=0.1

oriBody = Quaternion(Vector3(0, 0, 0), 0)
kwMeshes = {'color': [0, 1, 0], 'wire': True, 'dynamic': False, 'material': 'walls', 'wallMask': 62}
base_cylinder=O.bodies.append(geom.facetCylinder(center=(0,0,height/2), radius=radius, height=height, orientation=oriBody, segmentsNumber=20, **kwMeshes))

kwMeshes = {'color': [1, 0, 0], 'wire': True, 'dynamic': False, 'material': 'walls', 'wallMask': 4}
top_cylinder=O.bodies.append(geom.facetCylinder(center=(0,0,height+1.5*height), radius=radius, height=3*height, orientation=oriBody, segmentsNumber=20, **kwMeshes))


# Generate sphere packing in random positions within a cuboidal domain
sp=pack.SpherePack()
sp.makeCloud(minCorner=(-radius/sqrt(2),-radius/sqrt(2), 0.0),	# min corner of the domain
			 maxCorner=( radius/sqrt(2), radius/sqrt(2), 5*radius),	# max corner of the domain
			 rMean=0.01,					# average radius
			 rRelFuzz=0.3,					# relative deviation from average radius, e.g. +-30%
#			 num=100,					# number of spheres to be generated (optional)
			 periodic=False)				# whether to generate periodic packing or not 
sp.toSimulation(material='spheres')					# add spheres to the simulation


# Engines
O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(hertzian=True),Ig2_Wall_Sphere_ScGeom(hertzian=True)],  # collision "geometry"
                [Ip2_FrictMat_FrictMat_MindlinPhys(en=0.6)],  # collision "physics"
                [Law2_ScGeom_MindlinPhys_Mindlin()]  # contact law -- apply forces
        ),
        NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.0), # No numerical damping of energy in the integrator!!
        TranslationEngine(ids=top_cylinder, 
			velocity=0.1,
			translationAxis=Vector3(0, 0, 1),
			label='transEng',
			dead=True),
	VTKRecorder(fileName='results/angle-of-repose-vtk/AOR-', recorders=['all'], iterPeriod=500),
]


# Function to plot data periodically
def addPlotData():
	zMax=0
	i=0
	for b in O.bodies:
		if isinstance(b.shape,Sphere):
			if b.state.pos[2]>zMax:
				zMax=b.state.pos[2]

	plot.addData(i=O.iter, z=zMax)
	plot.saveDataTxt('results/angle-of-repose.txt')

plot.plots = {'i': ('z')}
plot.plot(noShow=True)

O.engines = O.engines + [PyRunner(command='addPlotData()', iterPeriod=500)]


# Function to pause the simulation if the system is relaxed
def checkUnbalancedForce():
	if unbalancedForce()<1e-2:
		if transEng.dead: # if the translation engine is not activated, then activate it!
			transEng.dead=False # start moving the upper cylinder
		else:
			checkUnb.dead=True
#			O.pause()
			
O.engines=O.engines+[PyRunner(command='checkUnbalancedForce()',iterPeriod=1000,label='checkUnb')]

# Timestep
O.dt = 0.8 * PWaveTimeStep() 

# Save the simulation, so that it can be reloaded later, for experimentation
O.saveTmp()

# Visualisation
v=qt.View()
rndr=qt.Renderer()

v.eyePosition = Vector3(0,-1.75, 0.45)
v.upVector    = Vector3(0, 0, 1)
v.viewDir     = Vector3(0, 1, 0)

colorStyle.setStyle('old',True)
#Gl1_Sphere.quality=2

Tutorial 3 - Granular flow in an hourglass

A random particle packing is generated within an hourglass. This example includes the use of a .stl mesh as a geometry.

# Script to simulate granular flow inside an hourglass
# 2025 © Vasileios Angelidakis <v.angelidakis@qub.ac.uk>
# Developed based on https://gitlab.com/yade-dev/trunk/blob/master/examples/hourglass/hourglass.py with modifications to adopt the viscoelastic Hertz-Mindlin contact law

# Create results directory
import os
import errno
try:
	os.mkdir('./results/')
except OSError as exc:
	if exc.errno != errno.EEXIST:
		raise
	pass

from yade import ymport, qt, pack

# Define materials
O.materials.append(FrictMat(young=5e5, poisson=.25, frictionAngle=0.0, density=2500, label='walls'))
O.materials.append(FrictMat(young=5e5, poisson=.25, frictionAngle=0.0, density=2500, label='spheres'))

# Import facets
facets=ymport.stl('hourglass.stl',material='walls',fixed=True,wire=True)
ids_list=O.bodies.append(facets)

# Packing generation: Three options, to generate an 1. Simple Cubic Packing (SCC) packing, 2. Hexagonal Close Packing packing (HCP), 3. Random Loose Packing
option=3

if option==1:	# Option 1: Generate a SCC sphere packing within a cylinder
	sp=pack.regularOrtho(pack.inCylinder(Vector3(0,0,0.032),Vector3(0,0,0.035),0.01),radius=0.0008,gap=0,material='spheres')
	O.bodies.append(sp)

elif option==2:	# Option 2: Generate a HCP sphere packing within a cylinder
	sp=pack.regularHexa(pack.inCylinder(Vector3(0,0,0.032),Vector3(0,0,0.035),0.01),radius=0.0008,gap=0,material='spheres')
	O.bodies.append(sp)

elif option==3:	# Option 3: Generate sphere packing in random positions within a cuboidal domain
	sp=pack.SpherePack()
	sp.makeCloud(
		minCorner=(-0.01,-0.01,0.020),	# min corner of the domain
		 maxCorner=( 0.01, 0.01,0.035),	# max corner of the domain
		 rMean=0.001,			# average radius
		 rRelFuzz=0.3,			# relative deviation from average radius, e.g. +-30% (optional)
		 seed=5,			# initialise the random number generator engine to get reproducible packings; do not define a seed, to get a different initial packing in each new run
		 num=100,			# number of spheres to be generated (optional)
		 periodic=False)		# whether to generate periodic packing or not 
	sp.toSimulation(material='spheres')	# add spheres to the simulation

# Engines
O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(hertzian=True)],  # collision geometry
                [Ip2_FrictMat_FrictMat_MindlinPhys(en=0.6,es=0.6)],  # collision "physics"
                [Law2_ScGeom_MindlinPhys_Mindlin()]  # contact law -- apply forces
        ),
        # Apply gravity force to particles. damping: numerical dissipation of energy.
        NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.0),
]

# Timestep
O.dt=0.8*PWaveTimeStep()

# Rotation engine: Add it to the O.engines after their initial definition
O.engines=O.engines + [
	RotationEngine(ids=ids_list, 
	angularVelocity=math.pi/0.5, #0.5 is half a period of rotation (in seconds)
	rotateAroundZero=True,
	rotationAxis=Vector3(1, 0, 0),
	zeroPoint=Vector3(0, 0, 0),
	label='rotEng',
	dead=True)
] # dead=True deactivates the engine

# Rotate the hourglass if it static or else stop it if it is alredy rotating
def startStepRotation():
	if (rotEng.dead):			
		rotEng.dead = False		# if the rotation engine is deactivated (i.e. dead=True) then activate it
	else:						
		rotEng.dead = True		# if the rotation engine is activated (i.e. dead=False) then deactivate it
		for b in O.bodies:		# reset the translational and angular velocities of the facets, after deactivating the engine
			if isinstance(b.shape,Facet):
				b.state.vel=[0,0,0]
				b.state.angVel=[0,0,0]

O.engines=O.engines + [PyRunner(iterPeriod=int(0.5 / O.dt), initRun=False, firstIterRun=25000, command='startStepRotation()')]

# Visualisation
v=qt.View()
v.ortho=True # Toggle orthographic view (parallel lines remain parallel) vs perspective view (for ortho=False)

v.eyePosition = Vector3(-0.15,0,0)
v.upVector    = Vector3(0,0,1)
v.viewDir     = Vector3(1,0,0)
colorStyle.setStyle('old',True)
Gl1_Sphere.quality=5
#v.sceneRadius = 0.10 # Radius of visible space. Defined automatically based on initial size of the system

O.saveTmp()

Tutorial 4 - Oedometer (Uniaxial compression) test

A granular packing is generated inside a facetBox() before it is compressed uniaxially.

# Script to simulate a uniaxial compression (oedometer) test with capabilities for batch processing
# 2025 © Vasileios Angelidakis <v.angelidakis@qub.ac.uk>
# Developed based on https://yade-dem.org/doc/tutorial-examples.html with modifications to adopt the viscoelastic Hertz-Mindlin contact law

# The components of the batch are:
# 1. table with parameters, one set of parameters per line (ccc.table)
# 2. readParamsFromTable which reads respective line from the parameter file
# 3. the simulation muse be run using yade-batch, not yade
#
# yade-batch --job-threads=1 oedometer.table oedometer.py
#

# Create results directory
import os
import errno
try:
	os.mkdir('./results/')
except OSError as exc:
	if exc.errno != errno.EEXIST:
		raise
	pass

# load parameters from file if run in batch
# default values are used if not run from batch
readParamsFromTable(rMean=.05, maxLoad=1e5, minLoad=1e4)
# make rMean, maxLoad accessible directly as variables later
from yade.params.table import *

# create box with free top, and ceate loose packing inside the box
from yade import pack, plot
O.bodies.append(geom.facetBox((.5, .5, .5), (.5, .5, .5), wallMask=31))
sp = pack.SpherePack()
sp.makeCloud((0, 0, 0), (1, 1, 1), rMean=rMean, rRelFuzz=0.2)
sp.toSimulation()

O.engines = [
	ForceResetter(),
	# sphere, facet, wall
	InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb(), Bo1_Wall_Aabb()]),
	InteractionLoop(
		# the loading plate is a wall, we need to handle sphere+sphere, sphere+facet, sphere+wall
		[Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom(hertzian=True), Ig2_Wall_Sphere_ScGeom(hertzian=True)],
		[Ip2_FrictMat_FrictMat_MindlinPhys(en=0.6,es=0.6)],
		[Law2_ScGeom_MindlinPhys_Mindlin()]
	),
	NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.0),
	PyRunner(command='checkUnbalanced()', realPeriod=2, label='checker'),
]
O.dt = .5 * PWaveTimeStep()

# the following checkUnbalanced, unloadPlate and stopUnloading functions are all called by the 'checker'
# (the last engine) one after another; this sequence defines progression of different stages of the
# simulation, as each of the functions, when the condition is satisfied, updates 'checker' to call
# the next function when it is run from within the simulation next time


# check whether the gravity deposition has already finished
# if so, add wall on the top of the packing and start the oedometric test
def checkUnbalanced():
	# at the very start, unbalanced force can be low as there is only few contacts, but it does not mean the packing is stable
	if O.iter < 5000:
		return
	# the rest will be run only if unbalanced is < .1 (stabilized packing)
	if unbalancedForce() > .1:
		return
	# add plate at the position on the top of the packing
	# the maximum finds the z-coordinate of the top of the topmost particle
	O.bodies.append(wall(max([b.state.pos[2] + b.shape.radius for b in O.bodies if isinstance(b.shape, Sphere)]), axis=2, sense=-1))
	global plate  # without this line, the plate variable would only exist inside this function
	plate = O.bodies[-1]  # the last particles is the plate
	# Wall objects are "fixed" by default, i.e. not subject to forces
	# prescribing a velocity will therefore make it move at constant velocity (downwards)
	plate.state.vel = (0, 0, -.1)
	# start plotting the data now, it was not interesting before
	O.engines = O.engines + [PyRunner(command='addPlotData()', iterPeriod=200)]
	# next time, do not call this function anymore, but the next one (unloadPlate) instead
	checker.command = 'unloadPlate()'


def unloadPlate():
	# if the force on plate exceeds maximum load, start unloading
	if abs(O.forces.f(plate.id)[2]) > maxLoad:
		plate.state.vel *= -1
		# next time, do not call this function anymore, but the next one (stopUnloading) instead
		checker.command = 'stopUnloading()'


def stopUnloading():
	if abs(O.forces.f(plate.id)[2]) < minLoad:
		# O.tags can be used to retrieve unique identifiers of the simulation
		# if running in batch, subsequent simulation would overwrite each other's output files otherwise
		# d (or description) is simulation description (composed of parameter values)
		# while the id is composed of time and process number
		plot.saveDataTxt('results/oedometer,' + O.tags['d.id'] + '.txt')
		O.pause()


def addPlotData():
	if not isinstance(O.bodies[-1].shape, Wall):
		plot.addData()
		return
	Fz = O.forces.f(plate.id)[2]
	plot.addData(force=Fz, displacement=-(plate.state.pos[2] - plate.state.refPos[2]), unbalanced=unbalancedForce(), i=O.iter, porosity=porosity())


# besides unbalanced force evolution, also plot the displacement-force diagram
plot.plots = {'i': ('unbalanced',), 'displacement': ('force',), 'force ': ('porosity',)}
plot.plot()

try:
	from yade import qt
	v=qt.View()
	colorStyle.setStyle('old',True)
	Gl1_Sphere.quality=2

	# Adjust camera view	
	v.eyePosition=Vector3(0.5,-1.61,0.48)
	v.upVector=Vector3(0,0,1)
	v.viewDir= Vector3(0,1,0)
except: pass

O.run()
# when running with yade-batch, the script must not finish until the simulation is done fully
# this command will wait for that (has no influence in the non-batch mode)
waitIfBatch()

Tutorial 5 - Triaxial compression test

A granular packing is generated inside a periodic space before it is compressed isotropically. Finally, it is compressed deviatorically.

# Script to simulate a periodic triaxial test
# 2025 © Vasileios Angelidakis <v.angelidakis@qub.ac.uk>
# Developed based on https://yade-dem.org/doc/tutorial-examples.html with minimal edits

# The initial packing is either:
#
# 1. random cloud with uniform distribution, or
# 2. cloud with specified granulometry (radii and percentages), or
# 3. cloud of clumps, i.e. rigid aggregates of several particles
#
# The triaxial consists of 2 stages:
#
# 1. isotropic compression, until sigmaIso is reached in all directions;
#    this stage is ended by calling compactionFinished()
# 2. deformation along the z-axis, while maintaining
#    constant stress (sigmaIso) laterally; this stage is ended by calling
#    triaxFinished()
#
# Controlling of strain and stresses is performed via PeriTriaxController,
# of which parameters determine type of control and also stability
# condition (maxUnbalanced) so that the packing is considered stabilized
# and the stage is done.
#

from __future__ import print_function
sigmaIso = -1e5

#import matplotlib
#matplotlib.use('Agg')

# generate loose packing
from yade import pack, qt, plot

O.periodic = True
sp = pack.SpherePack()
if 0:
	## uniform distribution
	sp.makeCloud((0, 0, 0), (2, 2, 2), rMean=.1, rRelFuzz=.3, periodic=True)
else:
	## create packing from clumps
	# configuration of one clump
	c1 = pack.SpherePack([((0, 0, 0), .03333), ((.03, 0, 0), .017), ((0, .03, 0), .017)])
	# make cloud using the configuration c1 (there could c2, c3, ...; selection between them would be random)
	sp.makeClumpCloud((0, 0, 0), (2, 2, 2), [c1], periodic=True, num=500)

# setup periodic boundary, insert the packing
sp.toSimulation()

O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb()]),
        InteractionLoop([Ig2_Sphere_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()]),
        PeriTriaxController(
                label='triax',
                # specify target values and whether they are strains or stresses
                goal=(sigmaIso, sigmaIso, sigmaIso),
                stressMask=7,
                # type of servo-control
                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='addPlotData()', iterPeriod=50),
]
O.dt = .5 * PWaveTimeStep()


def addPlotData():
	s11=min(triax.stress[0],triax.stress[1],triax.stress[2]) # We use min to get the stress of highest magnitude, since stresses are negative here
	s33=max(triax.stress[0],triax.stress[1],triax.stress[2]) # We use max to get the stress of lowest magnutude
	deviatorStress=s11-s33
	meanStress=(triax.stress[0]+triax.stress[1]+triax.stress[2])/3
	evol=triax.strain[0]+triax.strain[1]+triax.strain[2]
	try:
		ph=asin((s11-s33)/(s11+s33))
	except:
		ph=0

	plot.addData(
	        unbalanced=unbalancedForce(),
	        i=O.iter,
	        sxx=triax.stress[0],
	        syy=triax.stress[1],
	        szz=triax.stress[2],
	        exx=triax.strain[0],
	        eyy=triax.strain[1],
	        ezz=triax.strain[2],
	        q=deviatorStress, # deviator stress
		p=meanStress, # mean stress
		ev=evol, # volumetric strain
		phi=ph, # macroscopic friction angle
	        # save all available energy data
	        Etot=O.energy.total(),
	        **O.energy
	)


# enable energy tracking in the code
O.trackEnergy = True

# define what to plot
plot.plots = {
        'i': ('unbalanced',),
        'i ': ('sxx', 'syy', 'szz'),
        ' i': ('exx', 'eyy', 'ezz'),
        'i  ': ('q', 'p'),
        '  i': ('ev'),
        'i    ': ('phi'),
        # energy plot
        ' i ': (O.energy.keys, None, 'Etot'),
}
# show the plot
plot.plot()


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()

from yade import qt
v=qt.View()
colorStyle.setStyle('old',True)
Gl1_Sphere.quality=2

Tutorial 6 - Geogrid pull-out test

A pull-out test on a geogrid embedded in a granular material is conducted. The flexible structures are modelled using spherocylinders and deformable joints, using the GridConnection and GridNode shapes.

# Script to perform a pull out test of a geogrid
# 2025 © Vasileios Angelidakis <v.angelidakis@qub.ac.uk>

# Soil particles are simulated as spheres with visco-elastic Hertz-Mindlin type interactions, and rolling friction.
# The geogrid is simulated using deformable spherocylinders. 

from __future__ import print_function
from yade import pack,geom,plot
from builtins import range
from numpy import *
import os
import shutil

from yade.gridpfacet import *

#------------------------------------------------------------------------------
# Model parameters
particleShape='spheres'; sandType='Ottawa'; N=102; scale=30; gridSegments=4; mu=0.50; subfolder='gridSegments'


if 'description' in O.tags.keys(): 
	O.tags['id']=O.tags['description']
else:
	O.tags['id']=particleShape+'_scale_'+str(scale)+'_gridSegments_'+str(gridSegments)+'_mu_'+str(mu) #'BaselineSimulation' # Name of results subdirectory for simulations NOT ran in batch mode


stage=0 # stage 0 corresponds to the deposition of sand particles
	# stage 1 corresponds to the pulling-out of the geogrid


# ---------------------------------------------------------------------------------------------------------
# Engines
O.engines = [ForceResetter(),
	     InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_GridConnection_Aabb(),Bo1_Facet_Aabb()]),
	     InteractionLoop(
		[
		 Ig2_Sphere_Sphere_ScGeom(),
		 Ig2_Facet_Sphere_ScGeom(),
		 Ig2_Sphere_GridConnection_ScGridCoGeom(),
		 Ig2_GridNode_GridNode_GridNodeGeom6D(),
		 Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
                ], 
                [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True, setCohesionOnNewContacts=True),
		 Ip2_FrictMat_FrictMat_FrictPhys(),
		], 
		[Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),
		 Law2_ScGeom_FrictPhys_CundallStrack(),
		 Law2_ScGridCoGeom_FrictPhys_CundallStrack(),
		 Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),
		]
	),
	NewtonIntegrator(gravity=(0, 0, -9.81), damping=.2, label='newton'),
	TranslationEngine(ids=[0],velocity=-0.2,translationAxis=[1,0,0],label='transEng',dead=True)
]


# ---------------------------------------------------------------------------------------------------------
# Materials
O.materials.append(FrictMat(	young=50e6,  poisson=0.25,density=2650,frictionAngle=atan(mu),label='sandMat')) # Sand
O.materials.append(FrictMat(	young=50e6,  poisson=0.25,density=2650,frictionAngle=atan(0.0),label='containerMat')) # Container -> Same elastic properties as sand, but zero friction
O.materials.append(CohFrictMat( young=1.35e7,poisson=0.40,density= 950,frictionAngle=atan(mu), normalCohesion=5e100,shearCohesion=5e100,momentRotationLaw=True,label='gridMat')) # Geogrid

# ---------------------------------------------------------------------------------------------------------
# Domain - Lateral boundaries made of facets
dimX=0.33
dimY=0.33
dimZ=0.33

O.bodies.append(geom.facetBox((dimX/2, dimY/2, 0.23/2), (dimX/2, dimY/2, 0.23/2), wallMask=1+2+4+8+16, color=(0.2, 1, 0.2), wire=True,material='containerMat'))

seed=5
if particleShape[0:7]=='spheres':
	# ---------------------------------------------------------------------------------------------------------
	# Sphere packing - Below grid
	minCorner=[0,0,0]
	maxCorner=[dimX,dimY,dimZ]
	rMean=0.35e-3 * scale #For the real sands, D50=0.7mm, i.e. rMean=0.35mm
	rRelFuzz=0.15
	num=1000
	periodic=False
	sp = SpherePack()
	sp.makeCloud(minCorner=minCorner, maxCorner=maxCorner, rMean=rMean, rRelFuzz=rRelFuzz, num=num, periodic=periodic,seed=seed)
	sp.toSimulation(material='sandMat')



O.dt=0.8*PWaveTimeStep()

# ---------------------------------------------------------------------------------------------------------
# Run steps until the bottom spheres stabilise, and then create grid:
def checkUnbalancedForce():
	if unbalancedForce()<0.50 and avgNumInteractions(considerClumps=True)>=2:
		print('Iteration',O.iter,'- The initial sample is stable, generating geogrid.')

		# --------------------------------------------------------------------
		# Erase all particles above 0.10m to create a level layer
		for b in O.bodies:
			if isinstance(b.shape,Sphere):
				if b.state.pos[2]+b.shape.radius > 0.10: #this was 0.10
					try:
						if b.isClumpMember:
							O.bodies.deleteClumpBody(O.bodies[b.clumpId])
						else:
							O.bodies.erase(b.id)
					except: pass
		print('Iteration',O.iter,'- Initial packing levelled at a height of 0.10m.')

		generateGeogrid()
		checkUnb.dead=True
O.engines += [PyRunner(command='checkUnbalancedForce()',iterPeriod=1000,label='checkUnb')]

# ---------------------------------------------------------------------------------------------------------
# Run steps until the bottom spheres stabilise, and then create grid:
def generateTopPacking(): 
	if avgNumInteractions(considerClumps=True)>=2: #unbalancedForce()<0.2: # kineticEnergy()<1e-4
		print('Iteration',O.iter,'- The geogrid is installed and the system is stable, generating a second packing of particles.')
		plzMax=findzMax()
		
		# ---------------------------------------------------------------------------------------------------------
		# Find number of particles making the packing below the grid
		num=0
		for b in O.bodies:
			if isinstance(b.shape,Clump) or (isinstance(b.shape,Sphere) and not b.isClumpMember and not (isinstance(b.shape,GridConnection) or isinstance(b.shape,GridNode))):
				num=num+1			
#		print(num)
		
		if particleShape[0:7]=='spheres':
			# ---------------------------------------------------------------------------------------------------------
			# Sphere packing - Above grid
			zMax=findzMax()
			minCorner=[0,0,1.3*zMax]
			maxCorner=[dimX,dimY,dimZ+1.3*zMax]
			sp = SpherePack()
			sp.makeCloud(minCorner=minCorner, maxCorner=maxCorner, rMean=rMean, rRelFuzz=rRelFuzz, num=num, periodic=periodic,seed=seed) #2000 #15000
			sp.toSimulation(material='sandMat')

		else:
		#elif particleShape[0:7]=='clumps':
			# --------------------------------------------------------------------
			# Clump packing - Above grid
			sp = pack.SpherePack()
			sp.makeClumpCloud(	minCorner=[0,0,1.3*zMax],	        	# min corner of the domain
						maxCorner=[dimX,dimY,dimZ+1.3*zMax],	    	# max corner of the domain
						clumps=clumps,
						num=num,
						seed=seed,					# seed used to initialise a random number generator, and achieve repeatable packing conditions
						periodic=False)				# whether to generate periodic packing or not 
			sp.toSimulation(material='sandMat')				# add clumps to the simulation
			
			O.bodies.updateClumpProperties(discretization=25)		# get more accurate clump mass/volume/inertia for the clump particles


		genTopPack.dead=True
		pullGrid.dead=False

O.engines += [PyRunner(command='generateTopPacking()',iterPeriod=1000,label='genTopPack',dead=True)]


# ---------------------------------------------------------------------------------------------------------
# Find zMax
def findzMax():
	zMax=0
	for b in O.bodies:
		if isinstance(b.shape,Sphere):
			if b.state.pos[2]>zMax:
				zMax=b.state.pos[2]
	return zMax


# ---------------------------------------------------------------------------------------------------------
# Generate grid
def generateGeogrid():
	global nodesTG, dx, dy, numX, numY, colorTG
	zMax=findzMax()
	
	dx=0.033 # aperture size
	dy=0.033
	
	numX=11
	numY=11
	
	colorTG = [40. / 255., 102. / 255., 50. / 255.]
	nodesTG = []

	# Top grid (TG):
	# node
	for i in range(numX):
		for j in range(numY):
			nodesTG.append(O.bodies.append(gridNode([i * dx, j * dy, 1.3*zMax], 2.9e-3, wire=False, fixed=False, material='gridMat', color=colorTG)))
			O.bodies[-1].mask=1

	#connection
	for i in range(0, len(nodesTG)):
		for j in range(i + 1, len(nodesTG)):
			dist = (O.bodies[nodesTG[i]].state.pos - O.bodies[nodesTG[j]].state.pos).norm()
			if (dist <= min(dx,dy) * 1.01): # Find the next closest point to form a grid
				if gridSegments==4:
					pos0 = O.bodies[nodesTG[i]].state.pos # Position of start node for this grid cell
					pos4 = O.bodies[nodesTG[j]].state.pos # Position of end node for this grid cell
					
					pos1 = pos0 + (pos4-pos0)* 1/4
					pos2 = pos0 + (pos4-pos0)* 2/4
					pos3 = pos0 + (pos4-pos0)* 3/4
				
					bod1=O.bodies.append(gridNode(pos1, 2e-3, wire=False, fixed=False, material='gridMat', color=colorTG))
					bod2=O.bodies.append(gridNode(pos2, 2e-3, wire=False, fixed=False, material='gridMat', color=colorTG))				
					bod3=O.bodies.append(gridNode(pos3, 2e-3, wire=False, fixed=False, material='gridMat', color=colorTG))			
				
					O.bodies.append(gridConnection(nodesTG[i], bod1, 2e-3, material='gridMat', color=colorTG))
					O.bodies.append(gridConnection(bod1, bod2, 2e-3, material='gridMat', color=colorTG))
					O.bodies.append(gridConnection(bod2, bod3, 2e-3, material='gridMat', color=colorTG))
					O.bodies.append(gridConnection(bod3, nodesTG[j], 2e-3, material='gridMat', color=colorTG))
	genTopPack.dead=False
	
	

# ---------------------------------------------------------------------------------------------------------
# Pull out grid
def pullGeogrid():
	global stage
	if unbalancedForce()<0.5 and avgNumInteractions(considerClumps=True)>=2 and findzMax()<0.2: # and genTopPack.dead kineticEnergy()<1e-4
		print('Iteration',O.iter,'- The top sample is stable, pulling geogrid.')
		transEng.ids=nodesTG[0:numX]
		for i in nodesTG[0:numX]:
			O.bodies[i].state.blockedDOFs='xyz'

		for b in O.bodies:
			b.state.refOri=b.state.ori # Set this orientation (right before pulling) as the reference to measure all rotations
		
		for b in O.bodies:
			if isinstance(b.shape,Sphere):
				for i in range(0,numX-1,2):
					if i*dx <= b.state.pos[0] < (i+1)*dx:
						b.shape.color=[0.1,0.1,0.8]
					elif (i+1)*dx <= b.state.pos[0] < (i+2)*dx:
						b.shape.color=[0.8,0.1,0.1]
			if isinstance(b.shape,GridConnection) or isinstance(b.shape,GridNode):
				 b.shape.color=colorTG
		
		stage=1
		O.resetTime()			
		transEng.dead=False
		pullGrid.dead=True
		pD.dead=False # Start plotting and writing results
O.engines += [PyRunner(command='pullGeogrid()',iterPeriod=1000,label='pullGrid',dead=True)]


# --------------------------------------------------------------------
# Visualisation
try:
	from yade import qt
	v=qt.View()
	v.ortho=True # Toggle orthographic view (parallel lines remain parallel) vs perspective view (for ortho=False)

	v.viewDir=Vector3(0,-1,0)
#	v.eyePosition=Vector3(0.16,0.8,0.16)
	v.eyePosition=Vector3(0,1,0.235)
	v.upVector=Vector3(0,0,1)
	
	colorStyle.setStyle('old',True)

	rndr=qt.Renderer()
	rndr.ghosts=False
	rndr.wire=False
except: pass


##---------------------------------------------------------------------
# Function to plot data
def addPlotData():

	# Calculate average particle rotation
	avgRotation=0
	numParticles=0
	for b in O.bodies:
		if isinstance(b.shape,Sphere):
			numParticles=numParticles+1
#		avgRotation=avgRotation+b.state.rot().norm() # Calculate the RMS rotation of each particle
		avgRotation=avgRotation+abs(b.state.rot()[1]) # Calculate the RMS rotation of each particle
	avgRotation=degrees(avgRotation/numParticles)
	
#	Calculate porosity around the geogrid
	porosity=voxelPorosity(200,start=Vector3(0.03,0.03,0.03),end=Vector3(dimX-0.03,dimY-0.03,0.07))
#	voidRatio=porosity/(1-porosity)
	
	avgForceX=0
	avgDisplX=0
	for i in nodesTG[0:numX]:
		avgForceX=avgForceX+O.forces.f(i)[0]
		avgDisplX=avgDisplX+O.bodies[i].state.displ()[0]
	avgForceX=avgForceX/len(nodesTG[0:numX])
	avgDisplX=avgDisplX/len(nodesTG[0:numX])	
	
	deformation=-avgDisplX
	
	# Plot data
	plot.addData(
		i=O.iter,							# Iteration
		t=O.time,							# Virtual (simulation) time in seconds
		Fx=avgForceX,							# Average horizontal force on pulled grid nodes
		dx=deformation,							# Average horizontal displacement of pulled grid nodes
		unb=unbalancedForce(),						# Unbalanced force ratio
		por=porosity, 							# Local porosity of particles around the geogrid (at z=5+-2cm
#		voidR=voidRatio,						# Void ratio
		Z =avgNumInteractions(skipFree=False,considerClumps=True),	# Average coordination number of sphere-sphere contacts
		Zm=avgNumInteractions(skipFree=True,considerClumps=True),	# Mechanical coordination number of sphere-sphere contacts
#		sli=law.ratioSlidingContacts(),					# Ratio of sliding contacts as percentage of total contacts
		rot=avgRotation,						# Average sphere rotation in degrees
#		ani=anisotropy,							# Scalar quantification of contact force anisotropy (See Barreto et al, Powders and Grains 2009)
	)

	
	# Criterion to stop simulation
	if deformation>0.02:
		print('The simulation is completed')
		O.pause()

plot.plots={'i':('unb',),'i ':('Fx'),' i':('dx'),'i  ':('rot'), '  i':('Z','Zm'),
	'i   ':('por'), #'  i':('voidR') ,  'dZ ':('sli'),
}

plot.plot(noShow=False)

O.engines = O.engines + [PyRunner(command='addPlotData()',iterPeriod=100,label='pD',dead=True)]


O.run(False)

Tutorial 7 - Flow on flexible rockfall protection nets

The impact of a granular flow of debris on a flexible net is simulated. Flexible structures are modelled using spherocylinders and deformable joints, using the GridConnection and GridNode shapes in YADE, respectively. Planar, membrane-like elements also exist, using the PFacet class.

# Script to simulate a granular flow impacting a flexible barrier / rockfall net
# 2025 © Vasileios Angelidakis <v.angelidakis@qub.ac.uk>

from builtins import range
from yade import pack, geom, qt
from yade.gridpfacet import *
from pylab import *

O.engines = [
        ForceResetter(),
        InsertionSortCollider(
		        [   Bo1_Sphere_Aabb(),  
		            Bo1_Facet_Aabb(),
                    Bo1_GridConnection_Aabb(),
                ]),
        InteractionLoop(
                [   Ig2_GridNode_GridNode_GridNodeGeom6D(),
		            Ig2_Sphere_GridConnection_ScGridCoGeom(),
		            Ig2_Sphere_Sphere_ScGeom(),
                    Ig2_Facet_Sphere_ScGeom(hertzian=True)
		        ],
		        [   Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True, setCohesionOnNewContacts=False),
                    Ip2_FrictMat_FrictMat_MindlinPhys(en=0.6)
                ],
                [   Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),
                    Law2_ScGeom_MindlinPhys_Mindlin()
                ]
        ),
        NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.1, label='newton')
]

# Material of the flexible barrier
O.materials.append(
        CohFrictMat(
                young=3e4,
                poisson=0.3,
                density=1e1,
                frictionAngle=radians(10),
                normalCohesion=1e10,
                shearCohesion=1e10,
                momentRotationLaw=True,
                label='gridmat'
        )
)

# Material of the spheres
O.materials.append(
        FrictMat(
                young=3e2,
                poisson=0.3,
                density=1e1,
                frictionAngle=radians(10),
                label='spheremat'
        )
)

colorStyle.setStyle('old',True)

### Parameters of a rectangular grid ###
L = 0.11  #length [m]
l = 0.05  #width	[m]
nbL = 11  #number of nodes for the length	[#]
nbl = 5  #number of nodes for the width	[#]
r = L / 100.  #radius
color = [255. / 255., 102. / 255., 0. / 255.]
nodesIds = []
#Create all nodes first :
for i in range(0, nbL):
	for j in range(0, nbl):
		nodesIds.append(O.bodies.append(gridNode([0.4, i * L / nbL, j * l / nbl], r, wire=False, fixed=False, material='gridmat', color=color)))

#Create connection between the nodes
for i in range(0, len(nodesIds)):
	for j in range(i + 1, len(nodesIds)):
		dist = (O.bodies[i].state.pos - O.bodies[j].state.pos).norm()
		if (dist <= L / nbL * 1.01):
			O.bodies.append(gridConnection(i, j, r, color=color))

#Set a fixed node
O.bodies[0].dynamic = False
O.bodies[4].dynamic = False
O.bodies[50].dynamic = False
O.bodies[54].dynamic = False

extents = (0.25, 0.05, 0.05)
box_1 = O.bodies.append(geom.facetBox((extents[0], extents[1], extents[2]), extents=extents, wallMask=29, color=[0,1,0]))
wall = O.bodies.append(geom.facetBox((0.1, extents[1], extents[2]), (0,extents[1],extents[2]), wallMask=1, color=[1,0,0]))

sp=pack.SpherePack()
sp.makeCloud(minCorner=(0,0, 0),			# min corner of the domain
			 maxCorner=(0.1, 0.1, 0.2),	# max corner of the domain
			 rMean=0.005,			# average radius
			 rRelFuzz=0.2,			# relative deviation from average radius, e.g. +-30% (optional)
			 num=1000,			# number of spheres to be generated (optional)
			 periodic=False)		# whether to generate periodic packing or not 
sp.toSimulation(material='spheremat')			# add spheres to the simulation


def checkUnbalancedForce():
    #if unbalancedForce()<1e-1:
    if O.iter > 10000:
        O.bodies.erase(wall[0])
        O.bodies.erase(wall[1])
        checkUnb.dead=True
        theta = pi/6
        newton.gravity = [sin(theta)*9.81,0, -cos(theta)*9.81] 
			
O.engines=O.engines+[PyRunner(command='checkUnbalancedForce()',iterPeriod=100,label='checkUnb')]


O.dt = 0.8 * PWaveTimeStep()
O.saveTmp()
qt.View()