Segmentation fault(core dumped)

Asked by Cloud on 2019-05-05

Hi,

I’m working on triaxial test simulation of sand particle breaking. I'm trying fragment replace method recently. Refer to Jan's answer to this question: https://answers.launchpad.net/yade/+question/675477. My conditions for breaking are referenced to an article of Ciantia[1]. But segmentation falut occurs when this breaking condition is reached, catchsegv's output is as follow. Can you help me solve this problem? Thank you!

catchsegv's output:

Register dump:

 RAX: 000055ea62954c30 RBX: fffffffffffffff0 RCX: 000055da5e6d5560
 RDX: 000055da60668210 RSI: 000055da60592330 RDI: 000055da60586e00
 RBP: 000055ea62954c30 R8 : 000055da5e6cffe0 R9 : 000055da5e2c31c0
 R10: 000055da5e7349e0 R11: 00007f956a1923a0 R12: 0000000000000000
 R13: 0000000000000000 R14: 000055da60668210 R15: 00007f95b6ead0b0
 RSP: 00007f956bffe610

 RIP: 00007f95b8cfe55d EFLAGS: 00010206

 CS: 0033 FS: 0000 GS: 0000

 Trap: 0000000e Error: 00000004 OldMask: 00000000 CR2: 62954c30

 FPUCW: 0000037f FPUSW: 00000020 TAG: 00000000
 RIP: a1eee2c4 RDP: 00000000

 ST(0) ffff 8000000000000000 ST(1) 0000 0000000000000000
 ST(2) 0000 0000000000000000 ST(3) ffff f800000000000000
 ST(4) ffff 81ceb32c4b43fcf5 ST(5) ffff a000000000000000
 ST(6) ffff 81ceb32c4b43fcf5 ST(7) 81ce 81ceb32c4b43fcf5
 mxcsr: 1fa0
 XMM0: 000000000000000000000000ffffffff XMM1: 000000000000000000000000ffffffff
 XMM2: 000000000000000000000000ffffffff XMM3: 000000000000000000000000ffffffff
 XMM4: 000000000000000000000000ffffffff XMM5: 000000000000000000000000ffffffff
 XMM6: 000000000000000000000000ffffffff XMM7: 000000000000000000000000ffffffff
 XMM8: 000000000000000000000000ffffffff XMM9: 000000000000000000000000ffffffff
 XMM10: 000000000000000000000000ffffffff XMM11: 000000000000000000000000ffffffff
 XMM12: 000000000000000000000000ffffffff XMM13: 000000000000000000000000ffffffff
 XMM14: 000000000000000000000000ffffffff XMM15: 000000000000000000000000ffffffff

Backtrace:
/usr/lib/x86_64-linux-gnu/yade/libyade.so(_ZN13BodyContainer5eraseEib+0x4d)[0x7f95b8cfe55d]
/usr/lib/x86_64-linux-gnu/yade/py/yade/wrapper.so(_ZN5boost6python7objects23caller_py_function_implINS0_6detail6callerIM15pyBodyContainerFbibENS0_21default_call_policiesENS_3mpl7vector4IbRS5_ibEEEEEclEP7_objectSG_+0xf7)[0x7f95a6160017]
/usr/lib/x86_64-linux-gnu/libboost_python-py27.so.1.65.1(_ZNK5boost6python7objects8function4callEP7_objectS4_+0x285)[0x7f95b6eace55]
/usr/lib/x86_64-linux-gnu/libboost_python-py27.so.1.65.1(+0x22008)[0x7f95b6ead008]
/usr/lib/x86_64-linux-gnu/libboost_python-py27.so.1.65.1(_ZN5boost6python21handle_exception_implENS_9function0IvEE+0x73)[0x7f95b6eb39a3]
/usr/lib/x86_64-linux-gnu/libboost_python-py27.so.1.65.1(+0x1f6c9)[0x7f95b6eaa6c9]
/usr/bin/python(PyEval_EvalFrameEx+0x54a0)[0x55da5e2b0420]
/usr/bin/python(PyEval_EvalFrameEx+0x52b2)[0x55da5e2b0232]
/usr/bin/python(PyEval_EvalCodeEx+0x6da)[0x55da5e2a8d0a]
/usr/bin/python(PyEval_EvalCode+0x19)[0x55da5e2a8629]
/usr/bin/python(+0x12461f)[0x55da5e2d961f]
/usr/bin/python(PyRun_StringFlags+0x66)[0x55da5e3063d6]
/usr/bin/python(PyRun_SimpleStringFlags+0x3c)[0x55da5e306adc]
/usr/lib/x86_64-linux-gnu/yade/libyade.so(_Z11pyRunStringRKNSt7__cxx1112basic_stringIcSt11char_traitsIcESaIcEEE+0x1b)[0x7f95b987bd2b]
/usr/lib/x86_64-linux-gnu/yade/libyade.so(_ZN5Scene18moveToNextTimeStepEv+0x4b2)[0x7f95b8d91882]
/usr/lib/x86_64-linux-gnu/yade/libyade.so(_ZN14SimulationFlow12singleActionEv+0x33)[0x7f95b8d94d83]
/usr/lib/x86_64-linux-gnu/yade/libyade.so(_ZN12ThreadWorker16callSingleActionEv+0x4f)[0x7f95b8d9879f]
/usr/lib/x86_64-linux-gnu/yade/libyade.so(_ZN12ThreadRunner4callEv+0x47)[0x7f95b8d959b7]
/usr/lib/x86_64-linux-gnu/yade/libyade.so(_ZN12ThreadRunner3runEv+0x50)[0x7f95b8d95c90]
/usr/lib/x86_64-linux-gnu/yade/libyade.so(_ZN5boost6detail11thread_dataINS_9function0IvEEE3runEv+0x32)[0x7f95b8d97d52]
/usr/lib/x86_64-linux-gnu/libboost_thread.so.1.65.1(+0x11bcd)[0x7f95b5f54bcd]
/lib/x86_64-linux-gnu/libpthread.so.0(+0x76db)[0x7f95bc16c6db]
/lib/x86_64-linux-gnu/libc.so.6(clone+0x3f)[0x7f95bc4a588f]

breaking criterion:

def selectBodiesToCrush():
    breakparticle = []
    for b in O.bodies:
        if not isinstance(b.shape,Sphere):continue # exclude walls
        sigma_lim0 = 500e4
        diamter0 = 2.0
        diamter = 2.0 * b.shape.radius
        m = 10.0
        sigma_lim = sigma_lim0 * math.pow((diamter/diamter0), (-3.0/m))

        Fn_max = 0
        for i in b.intrs():
            Fn_abs = math.sqrt(i.phys.normalForce[0]*i.phys.normalForce[0]+i.phys.normalForce[1]*i.phys.normalForce[1]+i.phys.normalForce[2]*i.phys.normalForce[2])
            if Fn_abs > Fn_max:Fn_max = Fn_abs
            if Fn_abs == Fn_max:
                b1,b2=[O.bodies[ii] for ii in (i.id1,i.id2)]
                if not isinstance(b1.shape,Sphere):continue
                if not isinstance(b2.shape,Sphere):continue
                r_apostrophe = math.pow(1/b1.shape.radius + 1/b2.shape.radius,-1.0)
                E_apostrophe = math.pow(((1-b1.material.poisson*b1.material.poisson)/b1.material.young)+((1-b2.material.poisson*b2.material.poisson)/b2.material.young),-1.0)
                rH = math.pow((3*r_apostrophe)/(4*E_apostrophe),1.0/3.0)
                F_lim = math.pow((sigma_lim * math.pi * rH * rH),3.0)
                if Fn_max >= F_lim:
                    breakparticle.append(b)
    return breakparticle

def createChildren(bb):
    p = bb.state.pos
    r = bb.shape.radius
    # Particle replacement mode
    s1 = sphere(p+(0,0.5359*r,0),0.4641*r,material='spheres',color=(0,0,1))
    s2 = sphere(p+(-0.4641*r,-0.2679*r,0),0.4641*r,material='spheres',color=(0,0,1))
    s3 = sphere(p+(0.4641*r,-0.2679*r,0),0.4641*r,material='spheres',color=(0,0,1))
    s4 = sphere(p+(-0.6720*r,0.3880*r,0),0.2240*r,material='spheres',color=(0,0,1))
    s5 = sphere(p+(0.6720*r,0.3880*r,0),0.2240*r,material='spheres',color=(0,0,1))
    s6 = sphere(p+(0,-0.7760*r,0),0.2240*r,material='spheres',color=(0,0,1))
    s7 = sphere(p+(0,0,0.6339*r),0.3659*r,material='spheres',color=(0,0,1))
    s8 = sphere(p+(0,0,-0.6339*r),0.3659*r,material='spheres',color=(0,0,1))
    s9 = sphere(p+(0,-0.5942*r,0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s10 = sphere(p+(0.5146*r,0.2671*r,0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s11 = sphere(p+(-0.5146*r,0.2971*r,0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s12 = sphere(p+(0,-0.5942*r,-0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s13 = sphere(p+(0.5146*r,0.2971*r,-0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s14 = sphere(p+(-0.5146*r,0.2971*r,-0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    return s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,s14

def replace():
    bodiesToDelete = selectBodiesToCrush()
    for bb in bodiesToDelete:
        if isinstance(bb.shape,Sphere) and bb.shape.radius > minr/4:
            children = createChildren(bb)
            O.bodies.erase(bb.id)
            O.bodies.append(children)

[1]https://www.researchgate.net/publication/272794711_An_approach_to_enhance_efficiency_of_DEM_modelling_of_soils_with_crushable_grains

Regards,
Cloud

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Cloud
Solved:
2019-05-08
Last query:
2019-05-08
Last reply:
2019-05-06
Janek Kozicki (cosurgi) said : #1

We need a MWE to reproduce this. Not part of it. https://yade-dev.gitlab.io/trunk/user.html#questions-and-answers

Cloud (strife-cloud) said : #2

Hi, Janek,

Sorry, my MWE is as follow:

# -*- coding: utf-8 -*-
from yade import pack
import math

mn,mx=Vector3(0,0,0),Vector3(1,2,1)
young = 0.208e9
poisson = 0.3
compFricDegree = 15.0
finalFricDegree = 30.0
porosity = 0.95
num_spheres = 10000
confiningS = 50e3
isoconfining = 5e3
damp = 0.2
rate = -0.02
AxialStrainLimit = 0.3
key="triax_"
cd = 1 # cd=1 - consolidate drained cd=0- undrained
stabilityThreshold = 0.001
thick = 0.01

O.materials.append(FrictMat(
    young=young,
    poisson=poisson,
    density=2650,
    frictionAngle=radians(compFricDegree),
    label='spheres')
    )

O.materials.append(FrictMat(
    young=15e6,
    poisson=0.4,
    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,1.0/7.0,num_spheres,False, 0.95,seed=1)
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

O.dt=.1*utils.PWaveTimeStep()
O.usesTimeStepper=True

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

newton = NewtonIntegrator(damping=damp)

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Wall_Aabb(),Bo1_Facet_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_MindlinPhys()],
        [Law2_ScGeom_MindlinPhys_Mindlin()]
    ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8, defaultDt=4*utils.PWaveTimeStep()),
    triax,
    TriaxialStateRecorder(iterPeriod=1000,file='TriaxialRecorder'+key,truncate=1),
    newton,
]

triax.goal1=triax.goal2=triax.goal3=-isoconfining #(5kpa)

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

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

####### mini particle radius#####
minr=0
for b0 in O.bodies:
    if not isinstance(b0.shape,Sphere):continue
    if (b0.shape.radius)<=minr or (minr==0.0):minr=(b0.shape.radius)

#################################
##### Particle breakage #####
#################################
import math
def selectBodiesToCrush():
    breakparticle = []
    for b in O.bodies:
        if not isinstance(b.shape,Sphere):continue # exclude walls
        sigma_lim0 = 400e4
        diamter0 = 2.0
        diamter = 2.0 * b.shape.radius
        m = 10.0
        sigma_lim = sigma_lim0 * math.pow((diamter/diamter0), (-3.0/m))

        Fn_max = 0
        for i in b.intrs():
            Fn_abs = math.sqrt(i.phys.normalForce[0]*i.phys.normalForce[0]+i.phys.normalForce[1]*i.phys.normalForce[1]+i.phys.normalForce[2]*i.phys.normalForce[2])
            if Fn_abs > Fn_max:Fn_max = Fn_abs
            if Fn_abs == Fn_max:
                b1,b2=[O.bodies[ii] for ii in (i.id1,i.id2)]
                if not isinstance(b1.shape,Sphere):continue
                if not isinstance(b2.shape,Sphere):continue
                r_apostrophe = math.pow(1/b1.shape.radius + 1/b2.shape.radius,-1.0)
                E_apostrophe = math.pow(((1-b1.material.poisson*b1.material.poisson)/b1.material.young)+((1-b2.material.poisson*b2.material.poisson)/b2.material.young),-1.0)
                rH = math.pow((3*r_apostrophe)/(4*E_apostrophe),1.0/3.0)
                F_lim = math.pow((sigma_lim * math.pi * rH * rH),3.0)
                if Fn_max >= F_lim:
                    breakparticle.append(b)
    return breakparticle

def createChildren(bb):
    p = bb.state.pos
    r = bb.shape.radius
    s1 = sphere(p+(0,0.5359*r,0),0.4641*r,material='spheres',color=(0,0,1))
    s2 = sphere(p+(-0.4641*r,-0.2679*r,0),0.4641*r,material='spheres',color=(0,0,1))
    s3 = sphere(p+(0.4641*r,-0.2679*r,0),0.4641*r,material='spheres',color=(0,0,1))
    s4 = sphere(p+(-0.6720*r,0.3880*r,0),0.2240*r,material='spheres',color=(0,0,1))
    s5 = sphere(p+(0.6720*r,0.3880*r,0),0.2240*r,material='spheres',color=(0,0,1))
    s6 = sphere(p+(0,-0.7760*r,0),0.2240*r,material='spheres',color=(0,0,1))
    s7 = sphere(p+(0,0,0.6339*r),0.3659*r,material='spheres',color=(0,0,1))
    s8 = sphere(p+(0,0,-0.6339*r),0.3659*r,material='spheres',color=(0,0,1))
    s9 = sphere(p+(0,-0.5942*r,0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s10 = sphere(p+(0.5146*r,0.2671*r,0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s11 = sphere(p+(-0.5146*r,0.2971*r,0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s12 = sphere(p+(0,-0.5942*r,-0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s13 = sphere(p+(0.5146*r,0.2971*r,-0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s14 = sphere(p+(-0.5146*r,0.2971*r,-0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    return s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,s14

def replace():
    bodiesToDelete = selectBodiesToCrush()
    for bb in bodiesToDelete:
        if isinstance(bb.shape,Sphere) and bb.shape.radius > minr/4:
            children = createChildren(bb)
            O.bodies.erase(bb.id)
            O.bodies.append(children)

triax.internalCompaction=False

setContactFriction(radians(finalFricDegree))

if cd:
    triax.stressMask = 5
    triax.goal2 = rate
    triax.goal1 = -confiningS
    triax.goal3 = -confiningS
else:
    triax.stressMask = 0
    triax.goal2 = rate
    triax.goal1 = triax.goal3 = -rate/2

newton.damping=0.1

O.saveTmp()

from yade import plot

def history():
    plot.addData(
            e22=-triax.strain[1],
            e11=-triax.strain[0],
            e33=-triax.strain[2],
       ev=triax.volumetricStrain,
      s11=-triax.stress(triax.wall_right_id)[0],
      s22=-triax.stress(triax.wall_top_id)[1],
      s33=-triax.stress(triax.wall_front_id)[2],
            args=(-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_top_id)[1]-triax.stress(triax.wall_front_id)[2])/3,
            eta = 3*(-triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_right_id)[0])/(-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_top_id)[1]-triax.stress(triax.wall_front_id)[2]),
            devi = -triax.stress(triax.wall_top_id# -*- coding: utf-8 -*-
from yade import pack
import math

nRead=utils.readParamsFromTable(
        young=0.208e9,
        poisson = 0.3,
        compFricDegree=15.0,
        alphaKr = 0.2,
        Eta = 1.0,
        finalFricDegree=30.0,
        porosity = 0.95,
        num_spheres=10000,
        confiningS=50e3,
        isoconfining = 5e3,
        AxialStrainLimit=0.3,
        rate=-0.02,
        damp=0.2,
        unknownOk = True
)
from yade.params import table

mn,mx=Vector3(0,0,0),Vector3(1,2,1)
young = table.young
poisson = table.poisson
compFricDegree = table.compFricDegree
alphaKr = table.alphaKr
Eta = table.Eta
finalFricDegree = table.finalFricDegree
porosity = table.porosity
num_spheres = table.num_spheres
confiningS = table.confiningS
isoconfining = table.isoconfining
damp = table.damp
rate = table.rate
AxialStrainLimit=table.AxialStrainLimit
key="triax_"
cd = 1 # cd=1 - consolidate drained cd=0- undrained
stabilityThreshold = 0.001
thick = 0.01

O.materials.append(FrictMat(
    young=young,
    poisson=poisson,
    density=2650,
    frictionAngle=radians(compFricDegree),
    label='spheres')
    )

O.materials.append(FrictMat(
    young=15e6,
    poisson=0.4,
    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,1.0/7.0,num_spheres,False, 0.95,seed=1)
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

O.dt=.1*utils.PWaveTimeStep()
O.usesTimeStepper=True

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

newton = NewtonIntegrator(damping=damp)

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Wall_Aabb(),Bo1_Facet_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_MindlinPhys()],
        [Law2_ScGeom_MindlinPhys_Mindlin()]
    ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8, defaultDt=4*utils.PWaveTimeStep()),
    triax,
    TriaxialStateRecorder(iterPeriod=1000,file='TriaxialRecorder'+key,truncate=1),
    newton,
]

triax.goal1=triax.goal2=triax.goal3=-isoconfining #(5kpa)

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

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

####### mini particle radius#####
minr=0
for b0 in O.bodies:
    if not isinstance(b0.shape,Sphere):continue
    if (b0.shape.radius)<=minr or (minr==0.0):minr=(b0.shape.radius)

#################################
##### Particle breakage #####
#################################
import math
def selectBodiesToCrush():
    breakparticle = []
    for b in O.bodies:
        if not isinstance(b.shape,Sphere):continue # exclude walls
        sigma_lim0 = 400e4
        diamter0 = 2.0
        diamter = 2.0 * b.shape.radius
        m = 10.0
        sigma_lim = sigma_lim0 * math.pow((diamter/diamter0), (-3.0/m))

        Fn_max = 0
        for i in b.intrs():
            Fn_abs = math.sqrt(i.phys.normalForce[0]*i.phys.normalForce[0]+i.phys.normalForce[1]*i.phys.normalForce[1]+i.phys.normalForce[2]*i.phys.normalForce[2])
            if Fn_abs > Fn_max:Fn_max = Fn_abs
            if Fn_abs == Fn_max:
                b1,b2=[O.bodies[ii] for ii in (i.id1,i.id2)]
                if not isinstance(b1.shape,Sphere):continue
                if not isinstance(b2.shape,Sphere):continue
                r_apostrophe = math.pow(1/b1.shape.radius + 1/b2.shape.radius,-1.0)
                E_apostrophe = math.pow(((1-b1.material.poisson*b1.material.poisson)/b1.material.young)+((1-b2.material.poisson*b2.material.poisson)/b2.material.young),-1.0)
                rH = math.pow((3*r_apostrophe)/(4*E_apostrophe),1.0/3.0)
                F_lim = math.pow((sigma_lim * math.pi * rH * rH),3.0)
                if Fn_max >= F_lim:
                    breakparticle.append(b)
    return breakparticle

def createChildren(bb):
    p = bb.state.pos
    r = bb.shape.radius
    s1 = sphere(p+(0,0.5359*r,0),0.4641*r,material='spheres',color=(0,0,1))
    s2 = sphere(p+(-0.4641*r,-0.2679*r,0),0.4641*r,material='spheres',color=(0,0,1))
    s3 = sphere(p+(0.4641*r,-0.2679*r,0),0.4641*r,material='spheres',color=(0,0,1))
    s4 = sphere(p+(-0.6720*r,0.3880*r,0),0.2240*r,material='spheres',color=(0,0,1))
    s5 = sphere(p+(0.6720*r,0.3880*r,0),0.2240*r,material='spheres',color=(0,0,1))
    s6 = sphere(p+(0,-0.7760*r,0),0.2240*r,material='spheres',color=(0,0,1))
    s7 = sphere(p+(0,0,0.6339*r),0.3659*r,material='spheres',color=(0,0,1))
    s8 = sphere(p+(0,0,-0.6339*r),0.3659*r,material='spheres',color=(0,0,1))
    s9 = sphere(p+(0,-0.5942*r,0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s10 = sphere(p+(0.5146*r,0.2671*r,0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s11 = sphere(p+(-0.5146*r,0.2971*r,0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s12 = sphere(p+(0,-0.5942*r,-0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s13 = sphere(p+(0.5146*r,0.2971*r,-0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    s14 = sphere(p+(-0.5146*r,0.2971*r,-0.4456*r),0.2573*r,material='spheres',color=(0,0,1))
    return s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,s14

def replace():
    bodiesToDelete = selectBodiesToCrush()
    for bb in bodiesToDelete:
        if isinstance(bb.shape,Sphere) and bb.shape.radius > minr/4:
            children = createChildren(bb)
            O.bodies.erase(bb.id)
            O.bodies.append(children)

triax.internalCompaction=False

setContactFriction(radians(finalFricDegree))

if cd:
    triax.stressMask = 5
    triax.goal2 = rate
    triax.goal1 = -confiningS
    triax.goal3 = -confiningS
else:
    triax.stressMask = 0
    triax.goal2 = rate
    triax.goal1 = triax.goal3 = -rate/2

newton.damping=0.1

O.saveTmp()

from yade import plot

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

if 1:
    O.engines=O.engines[0:5]+[PyRunner(iterPeriod=500,command='history()',label='recorder')]+O.engines[5:7]
    O.engines=O.engines+[PyRunner(command='replace()',iterPeriod=100,label ='replace')]

############################
####### Start loading ##########
############################
while 1:
 O.run(1000,True)
 if O.iter%100 == 0:
    print "step = ", O.iter,"number of particles:",len(O.bodies)
         print "axialstrain = ",abs(triax.strain[1])
        if abs(triax.strain[1]) > AxialStrainLimit:
  break

plot.plots={'e22':('ev',),' e22 ':('eta')}
plot.plot()

plot.saveDataTxt('results'+key)
utils.waitIfBatch()
)[1]+triax.stress(triax.wall_right_id)[0],
            e=triax.porosity/(1-triax.porosity),
            i=O.iter,
            )

if 1:
    O.engines=O.engines[0:5]+[PyRunner(iterPeriod=500,command='history()',label='recorder')]+O.engines[5:7]
    O.engines=O.engines+[PyRunner(command='replace()',iterPeriod=100,label ='replace')]

###################################
####### Start loading ##########
###################################
while 1:
 O.run(1000,True)
 if O.iter%100 == 0:
    print "step = ", O.iter,"number of particles:",len(O.bodies)
         print "axialstrain = ",abs(triax.strain[1])
        if abs(triax.strain[1]) > AxialStrainLimit:
  break

plot.plots={'e22':('ev',),' e22 ':('eta')}
plot.plot()

plot.saveDataTxt('results'+key)

Janek Kozicki (cosurgi) said : #3

please make sure that yuor MWE reproduces the issue.

```
  File "./test.py", line 180
    from yade import pack
```
SyntaxError: invalid syntax

Janek Kozicki (cosurgi) said : #4

also make sure that it has no mixed indentation (use either tabs or spaces) and use brackets ( ) for your print function.

Cloud (strife-cloud) said : #5

It may be a time-step problem. After making some modification to time-step, it runs normally.