Yade

change O.engies

Asked by xjin

I write a script like this:
-------------------------------------------------------------

from yade import plot,pack,timing
import time, sys, os, copy

runGunplot=False
runInGui=False

QDZ=-0.01
#basic parameter

density=2350
young=3.5e9
poisson=.25
sigmaT=2.2e6
frictionAngle=radians(41.7)
epsCrackOnset=1e-2
relDuctility=55e0
rRelFuzz=0.1

#parameters of face
#parameters of diban
densitydiban=3e9
youngdiban=209e12
podiban=0.29
frictionAnglediban=float(atan(0.7))

#other parameters
intRadius=1.5
dtSafety=.8
damping=0.4
strainRateTension=.04
strainRateCompression=.5
setSpeeds=True

#parameter of the specimen
specimenLength=.07
specimenwidth=.05
specimenwidth2=specimenwidth/2.0
sphereRadius=3.5e-3

# isotropic confinement (should be negative)
isoPrestress=0

concreteId=O.materials.append(CpmMat(density=density,young=young,frictionAngle=frictionAngle,poisson=poisson,sigmaT=sigmaT,relDuctility=relDuctility,epsCrackOnset=epsCrackOnset,isoPrestress=isoPrestress))

dibancailiao=O.materials.append(FrictMat(density=densitydiban,young=youngdiban,poisson=podiban,frictionAngle=frictionAnglediban,label='dibancailiao'))

sps=SpherePack()
pred=pack.inCylinder((0,0,-0.01),(0,0,specimenLength+0.01),specimenwidth2)
sp=pack.randomDensePack(pred,spheresInCell=2000,rRelFuzz=rRelFuzz,radius=sphereRadius,memoizeDb='/home/dj/unix/Uniaxialyuan100.sqlite',returnSpherePack=True)
sps=sp.toSimulation(material=concreteId)

for s in sps:
 O.bodies[s].shape.color=(0.5,0.2,0.5)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb'),Bo1_Facet_Aabb(),],verletDist=.05*sphereRadius),
 InteractionLoop([Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius,label='ss2sc'),Ig2_Facet_Sphere_ScGeom(),],[Ip2_CpmMat_CpmMat_CpmPhys(),Ip2_FrictMat_CpmMat_FrictPhys(),],[Law2_ScGeom_CpmPhys_Cpm(),Law2_ScGeom_FrictPhys_CundallStrack(),],),
 NewtonIntegrator(gravity=(0,0,-9.8),damping=damping,label='damper'),
 CpmStateUpdater(realPeriod=.5),
 #UniaxialStrainer(active=strainactive,strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=setSpeeds,label='strainer'),
 PyRunner(iterPeriod=1,command='gaibian()',label='changestyle'),
 PyRunner(virtPeriod=1e-6/strainRateTension,realPeriod=1,command='addPlotData()',label='plotDataCollector',initRun=True),
 PyRunner(iterPeriod=1000,command='print1()',label='printparameter'),
 PyRunner(realPeriod=200,command='stopIfDamaged()',label='damageChecker'),
]

sl1=len(O.bodies)
L=specimenwidth/2.0

N=36
A=float(360/N)
H=specimenLength+0.02
nw=36
nh=2

#bot facets
for i in range(0,N):
 O.bodies.append(facet([(0,0,QDZ),(L*cos(i*2*pi/N),L*sin(i*2*pi/N),QDZ),(L*cos((i+1)*2*pi/N),L*sin((i+1)* 2*pi/N),QDZ)],fixed=True,material=dibancailiao))

botn=len(O.bodies)
bot=[]

for i in xrange(0,botn):
 bot.append(i)

#cylinder
rCyl2=L*2*0.5/cos(pi/float(nw))
facets=[]
for r in xrange(nw):
 for h in xrange(nh):
  v1=Vector3(rCyl2*cos(2*pi*(r+0)/float(nw)),rCyl2*sin(2*pi*(r+0)/float(nw)),H*(h+0)/float(nh)+QDZ)
  v2=Vector3(rCyl2*cos(2*pi*(r+1)/float(nw)),rCyl2*sin(2*pi*(r+1)/float(nw)),H*(h+0)/float(nh)+QDZ)
  v3=Vector3(rCyl2*cos(2*pi*(r+1)/float(nw)),rCyl2*sin(2*pi*(r+1)/float(nw)),H*(h+1)/float(nh)+QDZ)
  v4=Vector3(rCyl2*cos(2*pi*(r+0)/float(nw)),rCyl2*sin(2*pi*(r+0)/float(nw)),H*(h+1)/float(nh)+QDZ)
  f1=facet((v1,v2,v3),color=(0,0,1),material=dibancailiao)
  f2=facet((v1,v3,v4),color=(0,0,1),material=dibancailiao)
  facets.extend((f1,f2))

O.bodies.append(facets)

sl2=len(O.bodies)
for i in xrange(sl1,sl2):
 O.bodies[i].dynamic=False

def print1():
  print 'kz=',kz,'ub=',unbalancedForce(),'maxpos=',max([O.bodies[i].state.pos[2] for i in sps]),'cn=',avgNumInteractions()

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_00,coord_25,coord_50,coord_75,coord_100=mm+.0*(mx-mm),mm+.25*(mx-mm),mm+.5*(mx-mm),mm+.75*(mx-mm),mm+1.00*(mx-mm)
area_00,area_25,area_50,area_75,area_100=approxSectionArea(coord_00,axis),approxSectionArea(coord_25,axis),approxSectionArea(coord_50,axis),approxSectionArea(coord_75,axis),approxSectionArea(coord_100,axis)

plot.plots={'eps':('sigma',)}
oter=0
kz=5
def gaibian():
 global kz
 global runInGui
 global bb
 global negIds
 global posIds
 global axis
 global cro
 global mx
 global coord_00
 global coord_25
 global coord_50
 global coord_75
 global coord_100
 global area_00
 global area_25
 global area_50
 global area_75
 global area_100
 global oter
 if kz==5:
   if unbalancedForce()<0.3:
     kz=4
     Oter=O.iter
     runInGui=True
     for i in xrange(sl1,sl2):
      O.bodies.erase(i)
     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_00,coord_25,coord_50,coord_75,coord_100=mm+.0*(mx-mm),mm+.25*(mx-mm),mm+.5*(mx-mm),mm+.75*(mx-mm),mm+1.00*(mx-mm)
     area_00,area_25,area_50,area_75,area_100=approxSectionArea(coord_00,axis),approxSectionArea(coord_25,axis),approxSectionArea(coord_50,axis),approxSectionArea(coord_75,axis),approxSectionArea(coord_100,axis)
     O.engines=O.engines[:5]+[UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=setSpeeds,label='strainer'),]+O.engines[5:]
     O.pause()

def addPlotData():
 if kz==4:
  coord=coord_50
  area=area_50
  plot.addData({'t':O.time,'i':O.iter,'eps':strainer.strain,'sigma':strainer.avgStress+isoPrestress,'sigma':abs(forcesOnCoordPlane(coord,axis)[axis]/area+isoPrestress), })

def stopIfDamaged():
    if kz==4 and O.iter-oter>10:
     if len(plot.data['sigma'])>10:
        maxyinglizhi=max(plot.data['sigma'])
        zuijinm=plot.data['sigma'][-1]
        if 0.8*maxyinglizhi>=zuijinm:
         print O.iter
         O.pause()

#O.step()
ss2sc.interactionDetectionFactor=1.0
is2aabb.aabbEnlargeFactor=1.0

if runInGui:
 plot.plot()
 try:
  from yade import qt
  renderer=qt.Renderer()
 except:
  pass
_______________________________
when the script is running, when the O.iter>216691,I get error:
----------------------------------------------------
NameError Traceback (most recent call last)
/usr/bin/yade in <module>()

/usr/bin/yade in addPlotData()
    177 coord=coord_50
    178 area=area_50
--> 179 plot.addData({'t':O.time,'i':O.iter,'eps':strainer.strain,'sigma':strainer.avgStress+isoPrestress,'sigma':abs(forcesOnCoordPlane(coord,axis)[axis]/area+isoPrestress), })
    180
    181

NameError: global name 'strainer' is not defined
--------------------------------------------------------------
How can I solve it?
Thanks a lot!

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Robert Caulk
Solved:
Last query:
Last reply:
Revision history for this message
Best Robert Caulk (rcaulk) said : 		#1

Hello,

I am happy to help you here despite the fact that your question requires quite a bit of digging. For example, your script does not produce the error you mention, it looks for a non existent packing. So for documentation's sake please read through [1 https://yade-dem.org/wiki/Howtoask] before posting another question. Take close note of points 2(describe your problem) and 3 (create a minimal working example).

As for your question. Initialize UniaxialStrainer with the rest of your engines using dead=True. Then turn it on when you need it with strainer.dead=False and modify the members as you need (eg. strainer.axis=axis, strainer.crossSectionArea=crossSectionArea).

Cheers,

Robert

[1]https://yade-dem.org/wiki/Howtoask

from yade import plot,pack,timing
import time, sys, os, copy

runGunplot=False
runInGui=False

QDZ=-0.01
#basic parameter

density=2350
young=3.5e9
poisson=.25
sigmaT=2.2e6
frictionAngle=radians(41.7)
epsCrackOnset=1e-2
relDuctility=55e0
rRelFuzz=0.1

#parameters of face
#parameters of diban
densitydiban=3e9
youngdiban=209e12
podiban=0.29
frictionAnglediban=float(atan(0.7))

#other parameters
intRadius=1.5
dtSafety=.8
damping=0.4
strainRateTension=.04
strainRateCompression=.5
setSpeeds=True

#parameter of the specimen
specimenLength=.07
specimenwidth=.05
specimenwidth2=specimenwidth/2.0
sphereRadius=3.5e-3

# isotropic confinement (should be negative)
isoPrestress=0

concreteId=O.materials.append(CpmMat(density=density,young=young,frictionAngle=frictionAngle,poisson=poisson,sigmaT=sigmaT,relDuctility=relDuctility,epsCrackOnset=epsCrackOnset,isoPrestress=isoPrestress))

dibancailiao=O.materials.append(FrictMat(density=densitydiban,young=youngdiban,poisson=podiban,frictionAngle=frictionAnglediban,label='dibancailiao'))

sps=SpherePack()
pred=pack.inCylinder((0,0,-0.01),(0,0,specimenLength+0.01),specimenwidth2)
sp=pack.randomDensePack(pred,spheresInCell=2000,rRelFuzz=rRelFuzz,radius=sphereRadius,returnSpherePack=True) #memoizeDb='/home/dj/unix/Uniaxialyuan100.sqlite',
sps=sp.toSimulation(material=concreteId)

for s in sps:
 O.bodies[s].shape.color=(0.5,0.2,0.5)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb'),Bo1_Facet_Aabb(),],verletDist=.05*sphereRadius),
 InteractionLoop([Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius,label='ss2sc'),Ig2_Facet_Sphere_ScGeom(),],[Ip2_CpmMat_CpmMat_CpmPhys(),Ip2_FrictMat_CpmMat_FrictPhys(),],[Law2_ScGeom_CpmPhys_Cpm(),Law2_ScGeom_FrictPhys_CundallStrack(),],),
 NewtonIntegrator(gravity=(0,0,-9.8),damping=damping,label='damper'),
 CpmStateUpdater(realPeriod=.5),
UniaxialStrainer(blockDisplacements=False,blockRotations=False,setSpeeds=setSpeeds,label='strainer',dead=True), #UniaxialStrainer(active=strainactive,strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=setSpeeds,label='strainer'),
 PyRunner(iterPeriod=1,command='gaibian()',label='changestyle'),
 PyRunner(virtPeriod=1e-6/strainRateTension,realPeriod=1,command='addPlotData()',label='plotDataCollector',initRun=True),
 PyRunner(iterPeriod=1000,command='print1()',label='printparameter'),
 PyRunner(realPeriod=200,command='stopIfDamaged()',label='damageChecker'),
]

sl1=len(O.bodies)
L=specimenwidth/2.0

N=36
A=float(360/N)
H=specimenLength+0.02
nw=36
nh=2

#bot facets
for i in range(0,N):
 O.bodies.append(facet([(0,0,QDZ),(L*cos(i*2*pi/N),L*sin(i*2*pi/N),QDZ),(L*cos((i+1)*2*pi/N),L*sin((i+1)* 2*pi/N),QDZ)],fixed=True,material=dibancailiao))

botn=len(O.bodies)
bot=[]

for i in xrange(0,botn):
 bot.append(i)

#cylinder
rCyl2=L*2*0.5/cos(pi/float(nw))
facets=[]
for r in xrange(nw):
 for h in xrange(nh):
  v1=Vector3(rCyl2*cos(2*pi*(r+0)/float(nw)),rCyl2*sin(2*pi*(r+0)/float(nw)),H*(h+0)/float(nh)+QDZ)
  v2=Vector3(rCyl2*cos(2*pi*(r+1)/float(nw)),rCyl2*sin(2*pi*(r+1)/float(nw)),H*(h+0)/float(nh)+QDZ)
  v3=Vector3(rCyl2*cos(2*pi*(r+1)/float(nw)),rCyl2*sin(2*pi*(r+1)/float(nw)),H*(h+1)/float(nh)+QDZ)
  v4=Vector3(rCyl2*cos(2*pi*(r+0)/float(nw)),rCyl2*sin(2*pi*(r+0)/float(nw)),H*(h+1)/float(nh)+QDZ)
  f1=facet((v1,v2,v3),color=(0,0,1),material=dibancailiao)
  f2=facet((v1,v3,v4),color=(0,0,1),material=dibancailiao)
  facets.extend((f1,f2))

O.bodies.append(facets)

sl2=len(O.bodies)
for i in xrange(sl1,sl2):
 O.bodies[i].dynamic=False

def print1():
  print 'kz=',kz,'ub=',unbalancedForce(),'maxpos=',max([O.bodies[i].state.pos[2] for i in sps]),'cn=',avgNumInteractions()

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_00,coord_25,coord_50,coord_75,coord_100=mm+.0*(mx-mm),mm+.25*(mx-mm),mm+.5*(mx-mm),mm+.75*(mx-mm),mm+1.00*(mx-mm)
area_00,area_25,area_50,area_75,area_100=approxSectionArea(coord_00,axis),approxSectionArea(coord_25,axis),approxSectionArea(coord_50,axis),approxSectionArea(coord_75,axis),approxSectionArea(coord_100,axis)

plot.plots={'eps':('sigma',)}
oter=0
kz=5
def gaibian():
 global kz
 global runInGui
 global bb
 global negIds
 global posIds
 global axis
 global cro
 global mx
 global coord_00
 global coord_25
 global coord_50
 global coord_75
 global coord_100
 global area_00
 global area_25
 global area_50
 global area_75
 global area_100
 global oter
 if kz==5:
   if unbalancedForce()<0.3:
     kz=4
     Oter=O.iter
     runInGui=True
     for i in xrange(sl1,sl2):
      O.bodies.erase(i)
     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_00,coord_25,coord_50,coord_75,coord_100=mm+.0*(mx-mm),mm+.25*(mx-mm),mm+.5*(mx-mm),mm+.75*(mx-mm),mm+1.00*(mx-mm)
     area_00,area_25,area_50,area_75,area_100=approxSectionArea(coord_00,axis),approxSectionArea(coord_25,axis),approxSectionArea(coord_50,axis),approxSectionArea(coord_75,axis),approxSectionArea(coord_100,axis)
strainer.axis=axis
strainer.asymmetry=0
strainer.posIds=posIds
strainer.negIds=negIds
strainer.dead=False
strainer.strainRate=strainRateTension
strainer.crossSectionArea=crossSectionArea
O.pause()
O
def addPlotData():
 if kz==4:
  coord=coord_50
  area=area_50
  plot.addData({'t':O.time,'i':O.iter,'eps':strainer.strain,'sigma':strainer.avgStress+isoPrestress,'sigma':abs(forcesOnCoordPlane(coord,axis)[axis]/area+isoPrestress), })

def stopIfDamaged():
    if kz==4 and O.iter-oter>10:
     if len(plot.data['sigma'])>10:
        maxyinglizhi=max(plot.data['sigma'])
        zuijinm=plot.data['sigma'][-1]
        if 0.8*maxyinglizhi>=zuijinm:
         print O.iter
         O.pause()

#O.step()
ss2sc.interactionDetectionFactor=1.0
is2aabb.aabbEnlargeFactor=1.0

if runInGui:
 plot.plot()
 try:
  from yade import qt
  renderer=qt.Renderer()
 except:
  pass

Revision history for this message
xjin (jpeng22) said : 		#2

Thanks Robert Caulk, that solved my question.