yade-users team mailing list archive

[Zhoufeng Shi <question693011@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Answered => Solved

Zhoufeng Shi confirmed that the question is solved:
Thanks Jan for your kind reply!
So I tried the method you mentioned above to debug the problem. After I comment several condition lines that are not suitable for 2d situation and modify some values such as cellSize & cellFill in _getMemoizedPacking function, the code now can work well for saving 2d packing. Lines that I made changes are ended with long '####' symbol for clear inspection. 
However, as I am pretty new to both yade and python, sometimes I may not know what I am doing exactly and whether what I have done is correct or not. so the results can be wrong and I just mistaken it as a correct one. Any critiques and suggestions are welcomed!
the code are attached for future reference
#############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=True,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."), print(X), print(x1),print(scale); 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:
			print('x,y,z:',X,Y,Z)#################################################
			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:
			print('fullDim',fullDim[0],fullDim[1],fullDim[2])####################
			sp.isPeriodic = True
			sp.cellSize=(X,radius*(1+rRelFuzz),Z);################
			if fillPeriodic: sp.cellFill(Vector3(fullDim[0],radius*(1+rRelFuzz),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=True,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.)*3##########################
		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
		print('x1,y1,z1:',x1,y1,z1,'fullDim',fullDim[0],fullDim[1],fullDim[2])##########
		sp=_getMemoizedPacking(memoizeDb,radius,rRelFuzz,x1,y1,z1,fullDim,wantPeri,fillPeriodic=True,spheresInCell=spheresInCell,memoDbg=True)
		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),4.5),radius=.1,rRelFuzz=.5,spheresInCell=500,returnSpherePack=True,memoizeDb='/tmp/bending2DPackCache.sqlite')
sp.toSimulation()
####################### end of the code####################

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