yade-users team mailing list archive

[Rioual <question689131@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Answered => Open

Rioual is still having a problem:
Dear Robert,

You will find in the following the script.
This is inspired by oedometric-test.py (tutorial).
The original code (oedometric-test.py) uses all the cores available, I tested it.

Can you tell me me what could be the origin of the problem in my script ??
I had to add different translationEngines in my function called bt pyrunner:
can it be related ???  

Thanks,
All the best,

V.

#*****************************************************************************************************
#SCRIPT:
#*****************************************************************************************************

# gravity deposition (1) and controlled compression of the packing

# load parameters from file if run in batch
# default values are used if not run from batch
readParamsFromTable(rMean=0.1,rRelFuzz=.04,maxLoad=1e7,minLoad=1e5)
# 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))


## PhysicalParameters 
Density=2400
frictionAngle=radians(25)
#tc = 0.001
#en = 0.05
#et = 0.05
tc = 0.01
en = 0.0001
et = 0.0001

## Import wall's geometry
##params=utils.getViscoelasticFromSpheresInteraction(10e3,tc,en,es)
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
#fctIdsWS=O.bodies.append(ymport.stl('WS-wheelMN:q.stl',color=(1,0,0),material=facetMat))

fctIdscylinder=O.bodies.append(ymport.stl('Pot-
FR-1.stl',color=(1,0,0),material=facetMat))


print 'fctIdscylinder=', fctIdscylinder


###########################################
## Recuperer les dimensions des objets stl
###
minX = 1e99
maxX = -1e99
minY = 1e99
maxY = -1e99
minZ = 1e99
maxZ = -1e99

for facet in fctIdscylinder:
   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 

   minX = min(minX,min(v[0] for v in vs))
   maxX = max(maxX,max(v[0] for v in vs)) ###

   minY = min(minY,min(v[1] for v in vs))
   maxY = max(maxY,max(v[1] for v in vs)) ###

   minZ = min(minZ,min(v[2] for v in vs))
   maxZ = max(maxZ,max(v[2] for v in vs)) ###


print 'minX=',minX,'minY=',minY,'minZ=',minZ
print 'maxX=',maxX,'maxY=',maxY,'maxZ=',maxZ

Rcylext = maxX
Ycylmin = minY
Ycylmax = maxY

global Ycylmax
global Rcylext

###########################################


sp=pack.SpherePack()
sp.makeCloud((-23,10,-23),(23,310,23),rMean=rMean,rRelFuzz=rRelFuzz)
sp.toSimulation(color=(0,0,1),material=sphereMat)

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,-9.81,0),damping=0.1),
	# 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'),
	PyRunner(command='Poros()',iterPeriod=100,label='por'),
]
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<9000: return 
	# the rest will be run only if unbalanced is < .1 (stabilized packing)
	if unbalancedForce()>0.2: return


# Now, we pursue the job

# delete useless particles 1-At the top of the packing 2- On both sides of the packing 
	bodiesToBeDeleted = []
	for b in O.bodies:
	  if ( (b.state.pos[1] > (Ycylmax+rMean) or ((b.state.pos[0])**2+ (b.state.pos[2])**2> (Rcylext-3*rMean)**2) or (b.state.pos[1]<(Ycylmin))) ) and (isinstance(b.shape,Sphere)):
  		  bodiesToBeDeleted.append(b)
	for b in bodiesToBeDeleted:
  		O.bodies.erase(b.id)

 
	# add plate at the position on the top of the packing
	

        fctIdsbouchon=O.bodies.append(ymport.stl('Bouchon-
FR-1.stl',color=(1,0,0),material=facetMat))

        global fctIdsbouchon

        TransEngload=
TranslationEngine(ids=fctIdsbouchon,translationAxis=[0,-1,0],velocity=1,label='load')

        O.engines=O.engines+[TransEngload]

        global TransEngload

	# start plotting the data now, it was not interesting before
	O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=10)]
	# 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

        Fn = sum(O.forces.f(O.bodies[facetid].id)[1] for facetid in
fctIdsbouchon)

        print 'Fn=',Fn,'maxload=',maxLoad

        if abs(Fn)>maxLoad:

#		plate.state.vel*=-1
		TransEngload.Velocity = 0
		TransEngunload = TranslationEngine(ids=fctIdsbouchon,translationAxis=[0,1,0],velocity=0.1,label='unload')
		O.engines=O.engines+[TransEngunload]
		# next time, do not call this function anymore, but the next one (stopUnloading) instead
		checker.command='stopUnloading()'

def stopUnloading():


	Fn = sum(O.forces.f(O.bodies[facetid].id)[1] for facetid in fctIdsbouchon)

        print 'Fn=',Fn,'minload=',minLoad

        if abs(Fn)<minLoad:

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


#	On supprime le bouchon
		for facet in fctIdsbouchon:
			O.bodies.erase(facet)

		O.save('init_Final_state_packing.yade')	# on enregistre le resultats de la simulation au format .yade 
		print '********fin de la construction de l''empilement***********' # impression d'un message de fin
		O.pause()

	
def addPlotData():
#	if not isinstance(O.bodies[-1].shape,Wall):
#		plot.addData(); return

	Fn = sum(O.forces.f(O.bodies[facetid].id)[1] for facetid in fctIdsbouchon)
#	Fz=O.forces.f(plate.id)[2]

# Position du bouchon
	minXb = 1e99
	maxXb = -1e99
	minYb = 1e99
	maxYb = -1e99
	minZb = 1e99
	maxZb = -1e99

	for facet in fctIdsbouchon:
#   		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 

		minXb = min(minXb,min(v[0] for v in vs))
   		maxXb = max(maxXb,max(v[0] for v in vs)) ###

   	        minYb = min(minYb,min(v[1] for v in vs))
  		maxYb = max(maxYb,max(v[1] for v in vs)) ###

 		minZb = min(minZb,min(v[2] for v in vs))
   	        maxZb = max(maxZb,max(v[2] for v in vs)) ###


#	plot.addData(Fz=Fn,w=plate.state.pos[2]-plate.state.refPos[2],unbalanced=unbalancedForce(),i=O.iter)

	
	print 'minYb=',minYb,'maxYb=',maxYb

plot.addData(Fz=Fn,w=(minYb+maxYb)/2,unbalanced=unbalancedForce(),i=O.iter)

#*******************************************************************************************************************************
def Poros():

        dx = (Rcylext-2)/1.414
	
	por = utils.voxelPorosity(200,(-dx ,4,-dx),(dx,Ycylmax,dx))
	print 'tstep=',O.iter,'Porosite=',por

#*******************************************************************************************************************************
# besides unbalanced force evolution, also plot the displacement-force diagram
plot.plots={'i':('unbalanced',),'w':('Fz',)}
plot.plot()

O.run()
# when running with yade-batch, the script must not finish until the simulation is done fully
# this command will wait for that (has no influence in the non-batch mode)
waitIfBatch()


#***********************************************************************************************************************
#END OF SCRIPT
#***********************************************************************************************************************

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