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,

Sorry about all that. I modified my code to make it independant  with a faceted cylinder for penetration and Oedometric-test.py used to make the packing. Here it is below.
The simulation seems to be blocked for more than several minutes. Thee number of particles here is very low
 (example of packing taken from the tutorial examples oedometric-test.py).

All the best,

Vincent

The code:

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

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

#*********************************************************************************************************************
# 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.