yade-users team mailing list archive

[Jérôme Duriez <question254747@xxxxxxxxxxxxxxxxxxxxx>]

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

Jérôme Duriez gave more information on the question:
Just in case, and for completeness, I'm pasting below the script if someone wants to try to reproduce this feature / bug (I could reproduce with this script).
I apologize, this time I did not make any effort to shorten it, or even to translate it in english. It includes an internal compaction of a cloud, then unloading. The point is that it generates several .xml saves, for different mean pressures reached in descending order. 

You might try to laucnh, after that :
O.load('ech1conf5kPa.xml') # the first one generated
O.engines=[NewtonIntegrator(label='newt')] # the choice of the engine does not matter
print newt # => OK

Or 
O.load('ech1conf1000Pa.xml') # or e.g. ech1conf500Pa.xml
O.engines=[NewtonIntegrator(label='newt')] # the choice of the engine does not matter
print newt # => error....

So, here is the script :

# Les spheres
sp=SpherePack()
sp.makeCloud(maxCorner=Vector3(0.001,0.001,0.001),num=10000,porosity=0.74)
sphMat=O.materials.append( FrictMat(young=1e6,poisson=0.3,frictionAngle=radians(0.5)) )# cf Table 5.2 p.103 Khosravani, sauf frictionAngle = mise en place, cf p.103
sp.toSimulation(material=O.materials[sphMat])

# Les box
boxMat=O.materials.append( FrictMat(young=1e6,poisson=0.3,frictionAngle=radians(0)) ) # 'frictionless' 
rSph=O.bodies[0].shape.radius
[wallX2id,wallX1id,wallY2id,wallY1id,wallZ2id,wallZ1id]=O.bodies.append(aabbWalls(thickness=rSph,material=O.materials[boxMat])) # 1 = max, 2 = min

from yade import plot
nBod = O.bodies.__len__()
def computeThings():
    global volume,pBox,pLW
    dim = aabbDim()
    lX = dim[0]
    lY = dim[1]
    lZ = dim[2]
    volume = lX * lY * lZ
    matStr = - getStress(volume)
    pLW = (matStr[0,0] + matStr[1,1] + matStr[2,2]) /3.0
    signMaxZ = O.forces.f(wallZ1id)[2]/(lX*lY)
    signMinZ = -O.forces.f(wallZ2id)[2]/(lX*lY)
    signMaxY = O.forces.f(wallY1id)[1]/(lX*lZ)
    signMinY = -O.forces.f(wallY2id)[1]/(lX*lZ)
    signMaxX = O.forces.f(wallX1id)[0]/(lY*lZ)
    signMinX = -O.forces.f(wallX2id)[0]/(lY*lZ)
    pBox = (signMaxZ + signMinZ + signMaxY + signMinY + signMaxX + signMinY)/6.0
    uF = unbalancedForce()

    plot.addData(snZ1 = signMaxZ , snZ2 = signMinZ,
                 snY1 = signMaxY , snY2 = signMinY,
                 snX1 = signMaxX , snX2 = signMinX,
                 radi = O.bodies[0].shape.radius,uF=uF, it = O.iter, pLW=pLW)

O.engines=[ForceResetter()
           ,InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()])
           ,InteractionLoop(
               [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
               [Ip2_FrictMat_FrictMat_CapillaryPhys()],
               [Law2_ScGeom_FrictPhys_CundallStrack()]
           )
           ,NewtonIntegrator(damping=.2)
           ,PyRunner(command='computeThings()',iterPeriod=25)
          ]
O.dt=8.e-7 # 8e-7 d apres tests GSTS

# Grossissement des particules
multiplier=1.02
growParticles(multiplier)
O.run(100,True)
yade.qt.Controller(), yade.qt.View();
plot.plots={'it':'radi',' it': 'uF','it ':('snZ1','snZ2','snY1','snY2','snX1','snX2','pLW')}
plot.plot()

global volume,pBox,pLW
reached=0

while abs ((pBox+pLW)/2.0 < 0.95*5000):
    growParticles(multiplier)
    print 'Radius = ',O.bodies[0].shape.radius
    O.step()
    matStr=-getStress(volume) # la valeur de computeThings n est calculee que tous les 25 iterations, j ai besoin d une tres precise ici
    pLW = (matStr[0,0] + matStr[1,1] + matStr[2,2]) /3.0
    if pLW >= 0.95*5000:
        print '95% 5 kPa atteint'
        reached=1
    if reached:
        multiplier=1.001
        while unbalancedForce() > 0.005:
            print 'Recherche Equilibre depuis pLW ~', pLW
            O.run(100,True)
    O.run(400,True)
    
from yade import export
export.text('ech1conf5kPa.spheres')
O.save('ech1conf5kPa.xml')

print ''
print 'Stabilise a 5 kPa, debut decharge'
maxStrainRate = 1
vitesse = maxStrainRate * numpy.mean(aabbDim())
O.bodies[wallX2id].state.vel=Vector3(-vitesse,0,0)
O.bodies[wallX1id].state.vel=Vector3(vitesse,0,0)
O.bodies[wallY2id].state.vel=Vector3(0,-vitesse,0)
O.bodies[wallY1id].state.vel=Vector3(0,vitesse,0)
O.bodies[wallZ2id].state.vel=Vector3(0,0,-vitesse)
O.bodies[wallZ1id].state.vel=Vector3(0,0,vitesse)

stressThrPassed = [0,0,0,0]
stressThr = [1000,750,500,250]
def checkUnload():
    global volume
    matStr=-getStress(volume)
    pLW = (matStr[0,0] + matStr[1,1] + matStr[2,2]) /3.0 # checkUnload et computeThings ont meme iterPeriod, ils peuvent ttefois etre desynchronises en fonction de leurs iterations de depart
    for i in range(0,4):
        if pLW < stressThr[i] and not stressThrPassed[i]:
            export.text('ech1conf'+str(stressThr[i])+'Pa.spheres')
            O.save('ech1conf'+str(stressThr[i])+'Pa.xml')
            print 'p =', str(stressThr[i]), 'atteint a iter =', O.iter
            stressThrPassed[i]=1
            if i==3:
                O.pause()
        
O.engines=O.engines+[PyRunner(command='checkUnload()',iterPeriod=25)]
plot.plots={'it': 'uF','it ':('snZ1','snZ2','snY1','snY2','snX1','snX2','pLW')}
plot.plot()
O.run(wait=1)
plot.saveDataTxt('data_genEch1.txt')

-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.