yade-users team mailing list archive

[Cloud <question676370@xxxxxxxxxxxxxxxxxxxxx>]

New question #676370 on Yade:
https://answers.launchpad.net/yade/+question/676370

Hi, 

I'm trying to simulation triaxial test by using https://github.com/yade/trunk/blob/master/examples/triax-tutorial/script-session1.py wiht some minor modifications. It works well with young = 100e6, but it encountered some problems when I set young = 500e6. Most of this simulation works normally, but then there will be the following prompts in the terminal and the simulation stopped.

In [1]: Unhandled exception in thread started by <function liveUpdate at 0x7f0388a37320>
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
/usr/lib/x86_64-linux-gnu/yade/py/yade/plot.py in liveUpdate(timestamp)
    506                         for ax in axes:
    507                                 try:
--> 508                                         ax.relim() # recompute axes limits
    509                                         ax.autoscale_view()
    510                                 except RuntimeError: pass # happens if data are being updated and have not the same dimension at the very moment

/usr/lib/python2.7/dist-packages/matplotlib/axes/_base.pyc in relim(self, visible_only)
   1936         for line in self.lines:
   1937             if not visible_only or line.get_visible():
-> 1938                 self._update_line_limits(line)
   1939 
   1940         for p in self.patches:

/usr/lib/python2.7/dist-packages/matplotlib/axes/_base.pyc in _update_line_limits(self, line)
   1799         Figures out the data limit of the given line, updating self.dataLim.
   1800         """
-> 1801         path = line.get_path()
   1802         if path.vertices.size == 0:
   1803             return

/usr/lib/python2.7/dist-packages/matplotlib/lines.pyc in get_path(self)
    955         """
    956         if self._invalidy or self._invalidx:
--> 957             self.recache()
    958         return self._path
    959 

/usr/lib/python2.7/dist-packages/matplotlib/lines.pyc in recache(self, always)
    665             y = self._y
    666 
--> 667         self._xy = np.column_stack(np.broadcast_arrays(x, y)).astype(float)
    668         self._x, self._y = self._xy.T  # views
    669 

/usr/lib/python2.7/dist-packages/numpy/lib/stride_tricks.pyc in broadcast_arrays(*args, **kwargs)
    247     args = [np.array(_m, copy=False, subok=subok) for _m in args]
    248 
--> 249     shape = _broadcast_shape(*args)
    250 
    251     if all(array.shape == shape for array in args):

/usr/lib/python2.7/dist-packages/numpy/lib/stride_tricks.pyc in _broadcast_shape(*args)
    182     # use the old-iterator because np.nditer does not handle size 0 arrays
    183     # consistently
--> 184     b = np.broadcast(*args[:32])
    185     # unfortunately, it cannot handle 32 or more arguments directly
    186     for pos in range(32, len(args), 31):

ValueError: shape mismatch: objects cannot be broadcast to a single shape
axial strain reach Limitation, stop loading!

Has anybody met the same problem? How can I solve it ?

regards

Cloud

The whole code :
from yade import pack

num_spheres=5000
key='_My_triax_'
targetPorosity=0.40
compFricDegree=20
finalFricDegree=30
rate=-0.02
damp=0.2
stabilityThreshold=0.01
young=500e6
mn,mx=Vector3(0,0,0),Vector3(1,1,2)
confiningStress=100000

O.materials.append(FrictMat(young=young,poisson=0.1,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.1,frictionAngle=0,density=0,label='walls'))

walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0.1,num_spheres,False, 0.95,seed=1)
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

triax=TriaxialStressController(
    maxMultiplier=1.+2e4/young,
    finalMaxMultiplier=1.+2e3/young,
    thickness = 0,
    stressMask = 7,
    internalCompaction=True,
)

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
    triax,
    TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+key),
    NewtonIntegrator(damping=damp)
]

Gl1_Sphere.stripes=0

yade.qt.Controller(),yade.qt.View()

triax.goal1=triax.goal2=triax.goal3=-confiningStress

while 1:
    O.run(1000, True)
    unb=unbalancedForce()
    print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
    if unb < stabilityThreshold and abs(-confiningStress-triax.meanStress)/confiningStress<0.001:
        break

O.save('confinedState'+key+'.yade.gz')
print "###   Isotropic state saved   ###"

import sys
while triax.porosity > targetPorosity:
    compFricDegree = 0.95*compFricDegree
    setContactFriction(radians(compFricDegree))
    print "\r Friction: ", compFricDegree,"porosity:", triax.porosity,
    sys.stdout.flush()
    O.run(500,1)

O.save('compactedState'+key+'.yade.gz')
print "###   Compacted state saved   ###"

triax.internalCompaction=False

setContactFriction(radians(finalFricDegree))

triax.stressMask = 3
triax.goal3 = rate
triax.goal1=-confiningStress
triax.goal2=-confiningStress

newton.damping=0.1

O.saveTmp()

from yade import plot

def history():
    plot.addData(
            e11=-triax.strain[0],
            e22=-triax.strain[1],
            e33=-triax.strain[2],
            ev=triax.strain[0]+triax.strain[1]+triax.strain[2],
            s11=-triax.stress(triax.wall_right_id)[0],
            s22=-triax.stress(triax.wall_top_id)[1],
            s33=-triax.stress(triax.wall_front_id)[2],
            devi = -triax.stress(triax.wall_front_id)[2] + triax.stress(triax.wall_top_id)[1],
            i=O.iter)

def stoploading():
    axial_strain = abs(triax.strain[2])
    if axial_strain > 0.2:
        print 'axial strain reach Limitation, stop loading!'
        O.pause()

if 1:
    O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]
    O.engines=O.engines+[PyRunner(command='stoploading()',iterPeriod=10)]
else:
    O.engines[4]=PyRunner(iterPeriod=20,command='history()',label='recorder')

O.run(100,True)

plot.plots={'e33':('ev',),' e33 ':('devi')
}
plot.plot()

plot.saveDataTxt('results'+key)
plot.saveGnuplot('plotScript'+key)

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