yade-users team mailing list archive

[Jan Stránský <question688203@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Open => Answered

Jan Stránský proposed the following answer:
I have tried once more the code using Yade 2018.02b on Ubuntu 18.04 with
the same result - the problematic code is not reached at all..

Ibelow is the exact code I used (I had some problem with formatting the
copy-pasted code and got syntax errors)

If the code is problematic for you, then the yade version contributes to
the problem (and personally I don't want to search for problems with
2016 version..)

cheers
Jan

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