New question #693011 on Yade:
https://answers.launchpad.net/yade/+question/693011
Dear all,
I am trying to save the previously generated packing into file and load in the
next simulation just like what was done in uniax.py (using memoizeDb). however,
when I try to do the same for 2d randomDensePack, it always says that ' No
suitable packing in database found, running PERIODIC compression' and 'Packing
saved to the database /tmp/bending2DPackCache.sqlite'. So I wonder if there is
anything I need to do with '_memoizePacking' and '_getMemoizedPacking'.A sample
code is attached below and most of the code are copy from pack.py. Thanks in
advance!
###############code begin###################
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from yade import pack,plot,timing,qt
import time, sys, os, copy
def
_memoizePacking(memoizeDb,sp,radius,rRelFuzz,wantPeri,fullDim,noPrint=False):
import sys
if not memoizeDb: return
import pickle,sqlite3,time,os
if os.path.exists(memoizeDb):
conn=sqlite3.connect(memoizeDb)
else:
conn=sqlite3.connect(memoizeDb)
c=conn.cursor()
c.execute('create table packings (radius real, rRelFuzz real,
dimx real, dimy real, dimz real, N integer, timestamp real, periodic integer,
pack blob)')
c=conn.cursor()
if(sys.version_info[0]<3):
packBlob=buffer(pickle.dumps(sp.toList(),pickle.HIGHEST_PROTOCOL))
else:
packBlob=memoryview(pickle.dumps(sp.toList(),pickle.HIGHEST_PROTOCOL))
packDim=sp.cellSize if wantPeri else fullDim
c.execute('insert into packings values
(?,?,?,?,?,?,?,?,?)',(radius,rRelFuzz,packDim[0],0,packDim[2],len(sp),time.time(),wantPeri,packBlob,))
c.close()
conn.commit()
if not noPrint: print("Packing saved to the database",memoizeDb)
def
_getMemoizedPacking(memoizeDb,radius,rRelFuzz,x1,y1,z1,fullDim,wantPeri,fillPeriodic,spheresInCell,memoDbg=False,noPrint=False):
"""Return suitable SpherePack read from *memoizeDb* if found, None
otherwise.
:param fillPeriodic: whether to fill fullDim by repeating
periodic packing
:param wantPeri: only consider periodic packings
"""
import os,os.path,sqlite3,time,pickle,sys
if memoDbg and not noPrint:
def memoDbgMsg(s): print(s)
else:
def memoDbgMsg(s): pass
if not memoizeDb or not os.path.exists(memoizeDb):
if memoizeDb: memoDbgMsg("Database %s does not
exist."%memoizeDb)
return None
# find suitable packing and return it directly
conn=sqlite3.connect(memoizeDb); c=conn.cursor();
try:
c.execute('select
radius,rRelFuzz,dimx,dimy,dimz,N,timestamp,periodic from packings order by N')
except sqlite3.OperationalError:
raise RuntimeError("ERROR: database `"+memoizeDb+"' not
compatible with randomDensePack (corrupt, deprecated format or not a db created
by randomDensePack)")
for row in c:
R,rDev,X,Y,Z,NN,timestamp,isPeri=row[0:8]; scale=radius/R
rDev*=scale; X*=scale; Y*=scale; Z*=scale
memoDbgMsg("Considering packing
(radius=%g±%g,N=%g,dim=%g×%g×%g,%s,scale=%g), created
%s"%(R,.5*rDev,NN,X,Y,Z,"periodic" if isPeri else
"non-periodic",scale,time.asctime(time.gmtime(timestamp))))
if not isPeri and wantPeri: memoDbgMsg("REJECT: is not
periodic, which is requested."); continue
if wantPeri and (X/x1>0.9 or X/x1<0.6): memoDbgMsg("REJECT:
initSize differs too much from scaled packing size."); continue
if (rRelFuzz==0 and rDev!=0) or (rRelFuzz!=0 and rDev==0) or
(rRelFuzz!=0 and abs((rDev-rRelFuzz)/rRelFuzz)>1e-2): memoDbgMsg("REJECT:
radius fuzz differs too much (%g, %g desired)"%(rDev,rRelFuzz)); continue #
radius fuzz differs too much
if isPeri and wantPeri:
if spheresInCell>NN and spheresInCell>0:
memoDbgMsg("REJECT: Number of spheres in the packing too small"); continue
if abs((y1/x1)/(Y/X)-1)>0.3 or
abs((z1/x1)/(Z/X)-1)>0.3: memoDbgMsg("REJECT: proportions (y/x=%g, z/x=%g)
differ too much from what is desired (%g, %g)."%(Y/X,Z/X,y1/x1,z1/x1)); continue
else:
if (X<fullDim[0] or Y<fullDim[1] or Z<fullDim[2]):
memoDbgMsg("REJECT: not large enough"); continue # not large enough
memoDbgMsg("ACCEPTED");
if not noPrint: print("Found suitable packing in %s
(radius=%g±%g,N=%g,dim=%g×%g×%g,%s,scale=%g), created
%s"%(memoizeDb,R,rDev,NN,X,Y,Z,"periodic" if isPeri else
"non-periodic",scale,time.asctime(time.gmtime(timestamp))))
c.execute('select pack from packings where
timestamp=?',(timestamp,))
sp=SpherePack(pickle.loads(c.fetchone()[0]))
sp.scale(scale);
if isPeri and wantPeri:
sp.isPeriodic = True
sp.cellSize=(X,Y,Z);
if fillPeriodic:
sp.cellFill(Vector3(fullDim[0],0,fullDim[2]));###
#sp.cellSize=(0,0,0) # resetting cellSize avoids warning when
rotating
return sp
#if orientation: sp.rotate(*orientation.toAxisAngle())
#return filterSpherePack(predicate,sp,material=material)
#print "No suitable packing in database found, running",'PERIODIC
compression' if wantPeri else 'triaxial'
#sys.stdout.flush()
def
randomDensePack2d(predicate,radius,material=-1,dim=None,cropLayers=0,rRelFuzz=0.,spheresInCell=0,memoizeDb=None,useOBB=False,memoDbg=False,color=None,returnSpherePack=None,seed=0):
import sqlite3, os.path, pickle, time, sys, numpy
from math import pi
from yade import _packPredicates
wantPeri=(spheresInCell>0)
if 'inGtsSurface' in dir(_packPredicates) and
type(predicate)==inGtsSurface and useOBB:
center,dim,orientation=gtsSurfaceBestFitOBB(predicate.surf)
print("Best-fit oriented-bounding-box computed for GTS surface,
orientation is",orientation)
dim*=2 # gtsSurfaceBestFitOBB returns halfSize
else:
if not dim: dim=predicate.dim()
if max(dim)==float('inf'): raise RuntimeError("Infinite
predicate and no dimension of packing requested.")
center=predicate.center()
orientation=None
if not wantPeri: fullDim=tuple([dim[i]+4*cropLayers*radius for i in
(0,1,2)])
else:
# compute cell dimensions now, as they will be compared to ones
stored in the db
# they have to be adjusted to not make the cell to small WRT
particle radius
fullDim=dim
cloudPorosity=0.25 # assume this number for the initial cloud
(can be underestimated)
beta,gamma=fullDim[1]/fullDim[0],fullDim[2]/fullDim[0] # ratios
β=y₀/x₀, γ=z₀/x₀
N100=spheresInCell/cloudPorosity # number of spheres for cell
being filled by spheres without porosity
x1=radius*(1/(beta*gamma)*N100*(4/3.)*pi)**(1/3.)*5
y1,z1=beta*x1,gamma*x1; vol0=x1*y1*z1
maxR=radius*(1+rRelFuzz)
x1=max(x1,8*maxR); y1=max(y1,8*maxR); z1=max(z1,8*maxR);
vol1=x1*y1*z1
N100*=vol1/vol0 # volume might have been increased, increase
number of spheres to keep porosity the same
sp=_getMemoizedPacking(memoizeDb,radius,rRelFuzz,x1,y1,z1,fullDim,wantPeri,fillPeriodic=True,spheresInCell=spheresInCell,memoDbg=False)
if sp:
if orientation:
sp.cellSize=(0,0,0) # resetting cellSize avoids
warning when rotating
sp.rotate(*orientation.toAxisAngle())
return
filterSpherePack(predicate,sp,material=material,returnSpherePack=returnSpherePack)
else: print("No suitable packing in database found,
running",'PERIODIC compression' if wantPeri else 'triaxial')
sys.stdout.flush()
O.switchScene(); O.resetThisScene() ### !!
if wantPeri:
# x1,y1,z1 already computed above
sp=SpherePack()
O.periodic=True
#O.cell.refSize=(x1,y1,z1)
O.cell.setBox((x1,y1,z1))
#print cloudPorosity,beta,gamma,N100,x1,y1,z1,O.cell.refSize
#print x1,y1,z1,radius,rRelFuzz
O.materials.append(FrictMat(young=3e10,density=2400))
vmin,vmax=Vector3().Zero,O.cell.refSize
vmax[1]=0 ## diff
num=sp.makeCloud(vmin,vmax,radius,rRelFuzz,spheresInCell,True)
O.engines=[ForceResetter(),InsertionSortCollider([Bo1_Sphere_Aabb()],verletDist=.05*radius),InteractionLoop([Ig2_Sphere_Sphere_ScGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_ScGeom_FrictPhys_CundallStrack()]),PeriIsoCompressor(charLen=2*radius,stresses=[-100e9,-1e8],maxUnbalanced=1e-2,doneHook='O.pause();',globalUpdateInt=5,keepProportions=True),NewtonIntegrator(damping=.6)]
O.materials.append(FrictMat(young=30e9,frictionAngle=.5,poisson=.3,density=1e3))
for s in sp:
b=utils.sphere(s[0],s[1])
b.state.blockDOFs='yXZ'##diff
O.bodies.append(b)
O.dt=utils.PWaveTimeStep()
O.run(); O.wait()
sp=SpherePack(); sp.fromSimulation()
#print 'Resulting
cellSize',sp.cellSize,'proportions',sp.cellSize[1]/sp.cellSize[0],sp.cellSize[2]/sp.cellSize[0]
# repetition to the required cell size will be done below,
after memoizing the result
else:
assumedFinalDensity=0.6
V=(4.0/3.0)*pi*radius**3.0;
N=assumedFinalDensity*fullDim[0]*fullDim[1]*fullDim[2]/V;
TriaxialTest(
numberOfGrains=int(N),radiusMean=radius,radiusStdDev=rRelFuzz,
# upperCorner is just size ratio, if radiusMean is
specified
upperCorner=fullDim,
seed=seed,
## no need to touch any the following
noFiles=True,lowerCorner=[0,0,0],sigmaIsoCompaction=-4e4,sigmaLateralConfinement=-5e2,compactionFrictionDeg=1,StabilityCriterion=.02,strainRate=.2,thickness=0,maxWallVelocity=.1,wallOversizeFactor=1.5,autoUnload=True,autoCompressionActivation=False,
internalCompaction=True).load()
while ( numpy.isnan(utils.unbalancedForce()) or
utils.unbalancedForce()>0.005 ) :
O.run(500,True)
sp=SpherePack(); sp.fromSimulation()
O.switchScene() ### !!
_memoizePacking(memoizeDb,sp,radius,rRelFuzz,wantPeri,fullDim)
if wantPeri: sp.cellFill(Vector3(fullDim[0],maxR,fullDim[2]))##diff
if orientation:
sp.cellSize=(0,0,0); # reset periodicity to avoid warning when
rotating periodic packing
sp.rotate(*orientation.toAxisAngle())
return
filterSpherePack(predicate,sp,material=material,color=color,returnSpherePack=returnSpherePack)
###test
# tests
from yade import pack
#sp =
randomDensePack2d(pack.inAlignedBox((-1,-1,-1),(2,2,2)),radius=.1,rRelFuzz=.5,spheresInCell=200)
sp =
randomDensePack2d(pack.inSphere((0,0,0),3),radius=.1,rRelFuzz=.5,spheresInCell=200,returnSpherePack=True,memoizeDb='/tmp/bending2DPackCache.sqlite')
sp.toSimulation()
##################ends#######################
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
