yade-users team mailing list archive

[Rioual <question688203@xxxxxxxxxxxxxxxxxxxxx>]

Question #688203 on Yade changed:
https://answers.launchpad.net/yade/+question/688203

    Status: Needs information => Open

Rioual gave more information on the question:
Hello Jan,

Here is the code again.

Best

V.


******************************************************************

# gravity deposition, continuing with oedometric test after stabilization
# shows also how to run parametric studies with yade-batch

# 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 03-oedometric-test.table 
03-oedometric-test.py
#

# load parameters from file if run in batch
# default values are used if not run from batch
readParamsFromTable(rMean=.05,rRelFuzz=.3,maxLoad=1e6,minLoad=1e4)
# make rMean, rRelFuzz, 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=rRelFuzz)
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(),Ig2_Wall_Sphere_ScGeom()],
		[Ip2_FrictMat_FrictMat_FrictPhys()],
		[Law2_ScGeom_FrictPhys_CundallStrack()]
	),
	NewtonIntegrator(gravity=(0,0,-9.81),damping=0.5),
	# the label creates an automatic variable referring to this engine
	# we use it below to change its attributes from the functions called
	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(O.tags['d.id']+'.txt')
		print 'End building packing'
		O.pause()

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

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

O.run();O.wait()
#*********************************************************************************************************************
# Let the object penetrate the packing


#Start penetration

import time
from yade import qt
from yade import ymport
from yade.gridpfacet import *

###########################################
## PhysicalParameters
Density=2400
frictionAngle=radians(35)
tc = 0.001
en = 0.05
et = 0.05

## Import wall's geometry
facetMat=O.materials.append(ViscElMat(frictionAngle=frictionAngle,tc=tc, 
en=en, et=et)) # **params sets kn, cn, ks, cs
sphereMat=O.materials.append(ViscElMat(density=Density,frictionAngle=frictionAngle,tc=tc,en=en,et=et))
from yade import ymport
#******************************************
global TransEngload2
#******************************************
# Importation de la roue

global fctIdsWS

fctIdsWS=O.bodies.append( facetCylinder(center=(0.5,0.5,1.2), radius=0.3, 
height=0.2, orientation=Quaternion((1, 0, 0), 0), segmentsNumber=10, 
wallMask=7, angleRange=None, closeGap=False, radiusTopInner=-1, 
radiusBottomInner=-1, material=facetMat) )


## Timestep
O.dt=.2*tc

## Engines
O.engines=[
## Resets forces and momenta the act on bodies
ForceResetter(),
## Using bounding boxes find possible body collisions.
#InsertionSortCollider(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_GridConnection_Aabb()]),
## Interactions
InteractionLoop(
		# the loading plate is a facet, we need to handle sphere+sphere, 
sphere+facet
		[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom(),Ig2_Sphere_GridConnection_ScGridCoGeom()],
		[Ip2_FrictMat_FrictMat_FrictPhys()],
		[Law2_ScGeom_FrictPhys_CundallStrack()]
		),
## Apply gravity
GravityEngine(gravity=[0,0,-9.81]),
## Cundall damping must been disabled!
NewtonIntegrator(damping=0),


## Apply kinematics to wheel
PyRunner(command='kinematics_WS()',realPeriod=1,label='kine'),
]

from yade import qt
qt.View()


def kinematics_WS():


	h_WS = calc_h()[0]
	hmaxSpheres = calc_h()[1]


	print 'test0','h_WS=',h_WS,'hmaxSpheres=',hmaxSpheres

	TransEngload2 = 
TranslationEngine(ids=fctIdsWS,translationAxis=[0,0,-1],velocity=10)

        O.engines=O.engines+[TransEngload2]


	while h_WS > hmaxSpheres:

                print 'test01'

                TransEngload2.dead = False

                print 'test02'


######################### Here the simulation seems to be 
blocked#################################################################

                O.run(1,True)

#               O.engines=O.engines+[PyRunner(command='calc_h()')]


		print 'test03'

		h_WS = calc_h()[0]
		hmaxSpheres = calc_h()[1]

                print 'test1','h_WS=',h_WS,'hmaxSpheres=',hmaxSpheres

                TransEngload2.dead = True

        else:

                amTOT = sum((O.bodies[facetid].state.angMom)[1] for
facetid in fctIdsWS)


		while (amTOT< 1e-10):

O.engines=O.engines+[PyRunner(command='addTorque()',iterPeriod=1)]

                        amTOT = sum((O.bodies[facetid].state.angMom)[1]
for facetid in fctIdsWS)

                else:

# Stop simulation et measurement of the shear torque and cohesion			C_WS = 3* (imposed_T)/(2*Pi*(R0^3-R1^3))
			print 'C_WS=',C_WS

			print '********End of the calculation for Cohesion***********' #
			O.pause()

#*********************************************************************************************************************************
#This function calculate the height of the boundary h_WS and the maximum 
height of the packing hmax_Spheres
def calc_h():

###########################################
#Calculate h_WS

	minZ = 1e99
	maxZ = -1e99

	for facet in fctIdsWS:
#		print 'facet=', facet
   		facet= O.bodies[facet]		vs 
= [facet.state.pos + facet.state.ori * v for v in facet.shape.vertices] # 
vertices in global coord system
#		print 'vs=',vs

		minZ = min(minZ,min(v[2] for v in vs))
		maxZ = max(maxZ,max(v[2] for v in vs)) ####	print 'maxZ=',maxZ,'minZ=',minZ

        h_WS = maxZ

###########################################
#Calculate hmax_Spheres

	idHMax=0				# on definit une variable pour identifier quel corps a la 
hauteur la plus haute
	hMax=0.0				# initialisation de la hauteur a zero
	for i in O.bodies:			# on parcours tout les corps du systeme
		h=i.state.pos[2]		# on extrait la position selon l'axe y
		if (type(i.shape).__name__ == 'Sphere' and h>hMax):	#si le solide est une 
sphere et sa position est plus haut que hmax
			hMax=h			# on enregistre sa hauteur
			idHMax=i.id		# on garde en memoire de quel corps il s'agit
#	if (idHMax == 0): return rMax		# valeur de retour par default, rMAX
#	else: return (hMax+O.bodies[idHMax].shape.radius)

        hmaxSpheres = hMax+O.bodies[idHMax].shape.radius

        return h_WS,hmaxSpheres

#*****************************************************************************************************************************
#*****************************************************************************************************************************
# Add incremental torque on the WS boundary

def addTorque():


	increment_T= 0.001

	for Idfacets in fctIdsWS:
		O.forces.addT(Idfacets,(0,increment_T,0))
		imposed_T = O.forces.permT(Idfacets)
	return imposed_T

#*****************************************************************************************************************************

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.