problems in uniax.py

Asked by Chen Lifan on 2020-09-19

Hello:
I have some problems in the example uniax.py.
1. I found the number (by the command 'print len(O.bodies)') in the cell isn't the same as the value appointed by sphereInCell. When I change the shape of the packing (i.e. the values in pack.inHyperboloid), the number of speres in the packing changes as well. I wonder the reason for it.
2. How to prepare a 2D packing using pack.randomDensePack? I tried to give the 'dim' in it a value as 2, but error occured.
Any kind of help is appreciated. Thanks.

Here is the scripts and error:

scripts:

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division

from yade import plot,pack,timing
import time, sys, os, copy
readParamsFromTable(noTableOk=True, # unknownOk=True,
 young=24e9,
 poisson=.2,

 sigmaT=3.5e6,
 frictionAngle=atan(0.8),
 epsCrackOnset=1e-4,
 relDuctility=30,

 intRadius=1.5,
 dtSafety=.8,
 damping=0.4,
 strainRateTension=.05,
 strainRateCompression=.5,
 setSpeeds=True,
 # 1=tension, 2=compression (ANDed; 3=both)
 doModes=3,

 specimenLength=.15,
 sphereRadius=3.5e-3,

 # isotropic confinement (should be negative)
 isoPrestress=0,
)

from yade.params.table import *

if 'description' in O.tags.keys(): O.tags['id']=O.tags['id']+O.tags['description']

# make geom; the dimensions are hard-coded here; could be in param table if desired
# z-oriented hyperboloid, length 20cm, diameter 10cm, skirt 8cm
# using spheres 7mm of diameter
concreteId=O.materials.append(CpmMat(young=young,frictionAngle=frictionAngle,poisson=poisson,density=4800,sigmaT=sigmaT,relDuctility=relDuctility,epsCrackOnset=epsCrackOnset,isoPrestress=isoPrestress))

sps=SpherePack()
sp=pack.randomDensePack(pack.inHyperboloid((0,0,-.5*specimenLength),(0,0,.5*specimenLength),.5*specimenLength,.17*specimenLength),dim=2,spheresInCell=2000,radius=sphereRadius,memoizeDb='/tmp/triaxPackCache.sqlite',returnSpherePack=True)
#sp=pack.randomDensePack(pack.inAlignedBox((-.25*specimenLength,-.25*specimenLength,-.5*specimenLength),(.25*specimenLength,.25*specimenLength,.5*specimenLength)),spheresInCell=2000,radius=sphereRadius,memoizeDb='/tmp/triaxPackCache.sqlite',returnSpherePack=True)
sp.toSimulation(material=concreteId)
bb=uniaxialTestFeatures()
negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
O.dt=dtSafety*PWaveTimeStep()
print 'Timestep',O.dt

mm,mx=[pt[axis] for pt in aabbExtrema()]
coord_25,coord_50,coord_75=mm+.25*(mx-mm),mm+.5*(mx-mm),mm+.75*(mx-mm)
area_25,area_50,area_75=approxSectionArea(coord_25,axis),approxSectionArea(coord_50,axis),approxSectionArea(coord_75,axis)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb'),],verletDist=.05*sphereRadius),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius,label='ss2sc')],
  [Ip2_CpmMat_CpmMat_CpmPhys()],
  [Law2_ScGeom_CpmPhys_Cpm()],
 ),
 NewtonIntegrator(damping=damping,label='damper'),
 CpmStateUpdater(realPeriod=.5),
 UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=setSpeeds,label='strainer'),
 PyRunner(virtPeriod=1e-6/strainRateTension,realPeriod=1,command='addPlotData()',label='plotDataCollector',initRun=True),
 PyRunner(realPeriod=4,command='stopIfDamaged()',label='damageChecker'),
]
#O.miscParams=[Gl1_CpmPhys(dmgLabel=False,colorStrain=False,epsNLabel=False,epsT=False,epsTAxes=False,normal=False,contactLine=True)]

# plot stresses in ¼, ½ and ¾ if desired as well; too crowded in the graph that includes confinement, though
plot.plots={'eps':('sigma',)} #,'sigma.50')},'t':('eps')} #'sigma.25','sigma.50','sigma.75')}

O.saveTmp('initial');

O.timingEnabled=False

global mode
mode='tension' if doModes & 1 else 'compression'

def initTest():
 global mode
 print "init"
 if O.iter>0:
  O.wait();
  O.loadTmp('initial')
  print "Reversing plot data"; plot.reverseData()
 else: plot.plot(subPlots=False)
 strainer.strainRate=abs(strainRateTension) if mode=='tension' else -abs(strainRateCompression)
 try:
  from yade import qt
  renderer=qt.Renderer()
  renderer.dispScale=(1000,1000,1000) if mode=='tension' else (100,100,100)
 except ImportError: pass
 print "init done, will now run."
 O.step(); # to create initial contacts
 # now reset the interaction radius and go ahead
 ss2sc.interactionDetectionFactor=1.
 is2aabb.aabbEnlargeFactor=1.

 O.run()

def stopIfDamaged():
 global mode
 if O.iter<2 or not plot.data.has_key('sigma'): return # do nothing at the very beginning
 sigma,eps=plot.data['sigma'],plot.data['eps']
 extremum=max(sigma) if (strainer.strainRate>0) else min(sigma)
 minMaxRatio=0.5 if mode=='tension' else 0.5
 if extremum==0: return
 import sys; sys.stdout.flush()
 if abs(sigma[-1]/extremum)<minMaxRatio or abs(strainer.strain)>(5e-3 if isoPrestress==0 else 5e-2):
  if mode=='tension' and doModes & 2: # only if compression is enabled
   mode='compression'
   O.save('/tmp/uniax-tension.yade.gz')
   print "Saved /tmp/uniax-tension.yade.gz (for use with interaction-histogram.py and uniax-post.py)"
   print "Damaged, switching to compression... "; O.pause()
   # important! initTest must be launched in a separate thread;
   # otherwise O.load would wait for the iteration to finish,
   # but it would wait for initTest to return and deadlock would result
   import thread; thread.start_new_thread(initTest,())
   return
  else:
   print "Damaged, stopping."
   ft,fc=max(sigma),min(sigma)
   print 'Strengths fc=%g, ft=%g, |fc/ft|=%g'%(fc,ft,abs(fc/ft))
   title=O.tags['description'] if 'description' in O.tags.keys() else O.tags['params']
   print 'gnuplot',plot.saveGnuplot(O.tags['id'],title=title)
   print 'Bye.'
   O.pause()
   #sys.exit(0) # results in some threading exception

def addPlotData():
 yade.plot.addData({'t':O.time,'i':O.iter,'eps':strainer.strain,'sigma':strainer.avgStress+isoPrestress,
  'sigma.25':forcesOnCoordPlane(coord_25,axis)[axis]/area_25+isoPrestress,
  'sigma.50':forcesOnCoordPlane(coord_50,axis)[axis]/area_50+isoPrestress,
  'sigma.75':forcesOnCoordPlane(coord_75,axis)[axis]/area_75+isoPrestress,
  })
plot.plot(subPlots=False)
#O.run()
initTest()
waitIfBatch()

error:

Welcome to Yade 2018.02b
TCP python prompt on localhost:9000, auth cookie `uaekcd'
XMLRPC info provider on http://localhost:21000
Running script uniax.py
Traceback (most recent call last):
  File "/usr/bin/yade", line 182, in runScript
    execfile(script,globals())
  File "uniax.py", line 78, in <module>
    sp=pack.randomDensePack(pack.inHyperboloid((0,0,-.5*specimenLength),(0,0,.5*specimenLength),.5*specimenLength,.17*specimenLength),dim=2,spheresInCell=2000,radius=sphereRadius,memoizeDb='/tmp/triaxPackCache.sqlite',returnSpherePack=True)
  File "/usr/lib/x86_64-linux-gnu/yade/py/yade/pack.py", line 491, in randomDensePack
    if max(dim)==float('inf'): raise RuntimeError("Infinite predicate and no dimension of packing requested.")
TypeError: 'int' object is not iterable

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Zhoufeng Shi
Solved:
2020-09-20
Last query:
2020-09-20
Last reply:
2020-09-20
Best Zhoufeng Shi (zhoufeng-s) said : #1

hi lifan,
for Q2, please check #660243 for 2d random dense pack. it provides some instruction on how to modify the randomDensePack code to make it 2d.
As for Q1, you can also try increasing the periodic cell size in the randomDensePack code and the spheresincell will increase to the value as requested.

Chen Lifan (chenlifan) said : #2

Thanks Zhoufeng Shi, that solved my question.

Jan Stránský (honzik) said : #3

Hello,

1)
if spheresInCell > 0, then periodic compaction is created, containing spheresInCell number of particles.
This periodic packing is then copied wherever needed and then cropped to form the desired packing.
So, spheresInCell is independent of number of particles in the resulting packing.
It is especially advantageous for large initial packing (but is usable universally), e.g. 1M particles, where randomDensePack with 1M particles would take much time. This way, you can generate a periodic cube with e.g. 1k particles and quickly copy it in the space.
With this approach, the packing generation has "almost" O(1) complexity.

**Personally** I do not use randomDensePack without spheresInCell parameter

cheers
Jan

Zhoufeng Shi (zhoufeng-s) said : #4

hi Jan
if the spheresInCell is not approximately equal to the number of particles in resulting packing(eg. sphereincel/final particles in packing =100), there seems to be clear pattern which is a kind of regularity that I don't desire. So what is the best ratio between sphereInCell and final particle number to make sure its randomness and save computation cost at the same time?
regards,
zhoufeng

Jan Stránský (honzik) said : #5

> there seems to be clear pattern which is a kind of regularity

there **is** a regular pattern, by design (one periodic packing is copied wherever needed to fill the final packing)

> So what is the best ratio ...

You have to choose. According to your definition of "best", "randomness" etc.

cheers
Jan

Zhoufeng Shi (zhoufeng-s) said : #6

thanks Jan, that is very clear