yade-users team mailing list archive

[Seti <question295814@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Answered => Open

Seti is still having a problem:
Hi Jerome/Bruno,


Thanks for reply. Lets make it simple. Would you please help me to understand why below (simple) script does not work?

Thank you so much for your time.

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

utils.readParamsFromTable(rMean=0.02,rRelFuzz=0.1,maxLoad=10e8,minLoad=10)
from yade.params.table import *
from yade import pack, plot

###### Define the material for your facet here
O.bodies.append(utils.geom.facetBox((0.5,0.5,0.5),(0.5,0.5,0.5),wallMask=31))

mat1 = O.materials.append(CohFrictMat(young=30e9,poisson=0.3,frictionAngle=radians(30),momentRotationLaw=True,etaRoll=1,density=2600,isCohesive=True,normalCohesion=1e6, shearCohesion=80e6,label='spheres1'))
mat2 = O.materials.append(CohFrictMat(young=30e9,poisson=0.3,frictionAngle=radians(30),momentRotationLaw=True,etaRoll=1,density=2600,isCohesive=True,normalCohesion=1e6, shearCohesion=80e6,label='spheres2'))

sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1,1,1),rMean=rMean,rRelFuzz=rRelFuzz)

##### add the sphere pack to the simulation
sp.toSimulation()

##### add the sphere pack to the simulation


O.materials.append(FrictMat(young=1e10,poisson=0.1, frictionAngle = radians(0) , label='wallmat'))
wallmat = O.materials[-1]

O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
		 InteractionLoop([Ig2_Sphere_Sphere_L3Geom(),Ig2_Facet_Sphere_L3Geom(),Ig2_Wall_Sphere_L3Geom()],
			[Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(normalCohesion=MatchMaker(matches=((mat1,mat2,50),(mat2,mat2,100))),shearCohesion=MatchMaker(matches=((mat1,mat2,5000),(mat2,mat2,80000))),label="cohesiveIp")],
			[Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(
			label='cohesiveLaw')]

		 ),
	#NewtonIntegrator(damping=0.4,gravity=(0,0,-9.81)),
        NewtonIntegrator(damping=0.4),
	PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'), 
	PyRunner(command='addPlotData()',iterPeriod=100)
	]


O.dt=20*utils.PWaveTimeStep()
O.trackEnergy=True

def checkUnbalanced():
	
	global Forcewall 
	if O.iter<5000: return
	if utils.unbalancedForce()>.1: return
	
	plot.saveDataTxt('Elaheh3.txt.bz2')
	
	
	O.bodies.append(utils.wall((0.5,0.5,1.0),2,sense = -1, material = 'wallmat')) 
	Forcewall = O.bodies[-1]
	
	
	Forcewall.state.vel = (0,0,-0.01) ### applying the velocity to the wall (0,0,-1)=(u,v,w)
	
	O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=200)]
  
	checker.command='unloadPlate()'
def unloadPlate():
   # if the force on Forcewall exceeds maximum load, start unloading
   if abs(O.forces.f(Forcewall.id)[2])>maxLoad:
      Forcewall.state.vel = (0,0,0) ##you hadn't defined Forcewall before. so it gave you an error
      # next time, do not call this function anymore, but the next one (stopUnloading) instead
      checker.command='stopUnloading()'

def stopUnloading():
   if abs(O.forces.f(Forcewall.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')
    

def addPlotData():
    if not isinstance(O.bodies[-1].shape,Wall):
       plot.addData(); return     
    Fz=O.forces.f(Forcewall.id)[2]
    plot.addData(Fz=Fz,w=Forcewall.state.pos[2]-Forcewall.state.refPos[2])#,unbalanced=utils.unbalancedForce(),i=O.iter)

plot.plots={'w':('Fz',)}#,'i':('unbalanced',)
plot.plot()


O.run()
utils.waitIfBatch()
plot.saveDataTxt('m.txt.bz2',vars=('Fz','w'))

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