double free or corruption (out)

hi,

I have these errors "double free or corruption (out)/segment error/bus error" a few seconds after I click the run button of my script. It is basicly about generating a cylindrical specimen consisted of spheres and cylinders. My yade version is 2018.02b and I use ubuntu18.04. Below is my code.

#====================================================================
from yade.gridpfacet import *
from yade import pack, plot
#from random import random as r
import numpy as np
from numpy import *
import math

#parameters
# target size after equilibrium (before shearing): cylinder with hight=80e-3, radius=20e-3
# target size before equilibrium (initial size): cylinder with hight=800e-3, radius=20e-3
height = 500e-3
width = 40e-3
#rParticle = 0.11*1e-3/2
#rRelFuzz=.0005
mi,ma = (-20e-3,-20e-3,0),(20e-3,20e-3,500e-3)
nCyls,nSphs = 80,140
frictionAngleSph=30
frictionAngleCyl=30

#=============================meterials========================================
O.materials.append(FrictMat(young=4.0e6,poisson=.3,frictionAngle=30,density=2630e0,label='sphereMat'))#for sphere

O.materials.append(FrictMat(young=3e9,poisson=.15,frictionAngle=20,density=910e+0,label='extcylMat'))#for sphere-cylinder
O.materials.append(CohFrictMat(young=3e9,poisson=.15,density=910e0,frictionAngle=20,normalCohesion=1e40,shearCohesion=1e40,momentRotationLaw=True,label='intcylMat'))#for cylinder-cylinder

O.materials.append( CohFrictMat( young=3e8,poisson=0.15,density=910e0,frictionAngle=20,normalCohesion=3e100,shearCohesion=3e100,momentRotationLaw=True,label='gridNodeMat' ) )#for gridNodes
#O.materials.append(CohFrictMat(young=3e9,poisson=.15,density=910e6,frictionAngle=20,normalCohesion=1e40,shearCohesion=1e40,momentRotationLaw=True,label='gridNodeMat'))#for gridNodes
O.materials.append(FrictMat(young=4e6,poisson=0.3,density=1000,frictionAngle=20,label='pFacetMat')) #for pfacet

#==============================Engines=========================================
O.engines=[
 ForceResetter(),
 InsertionSortCollider([
  Bo1_Sphere_Aabb(),
  Bo1_GridConnection_Aabb(),
  Bo1_Facet_Aabb()
 ]),
 InteractionLoop([
  Ig2_Sphere_Sphere_ScGeom(),
  Ig2_Box_Sphere_ScGeom(),
  Ig2_GridNode_GridNode_GridNodeGeom6D(),
  Ig2_Sphere_GridConnection_ScGridCoGeom(),
  Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
  Ig2_Facet_Sphere_ScGeom(),
 ],
 [
  Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False), # internal cylinder physics
  Ip2_FrictMat_FrictMat_FrictPhys(), # physics for external interactions, i.e., cylinder-cylinder, sphere-sphere, cylinder-sphere

  Ip2_FrictMat_CpmMat_FrictPhys(),
 ],
 [
  Law2_ScGeom_FrictPhys_CundallStrack(), # contact law for sphere-sphere
  Law2_ScGridCoGeom_FrictPhys_CundallStrack(), # contact law for cylinder-sphere
  Law2_ScGeom6D_CohFrictPhys_CohesionMoment(), # contact law for "internal" cylinder forces
  Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(), # contact law for cylinder-cylinder interaction
  #Law2_GridCoGridCoGeom_CohFrictPhys_CundallStrack(),
  Law2_ScGeom_CpmPhys_Cpm(),
 ]
 ),
 GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.5,label='ts'),
 NewtonIntegrator(gravity=(0,0,-9.81),damping=0.5,label='newton'),
]

#=======================================sphere=================================================
sp = yade.pack.SpherePack()
sp.makeCloud(mi,ma,psdSizes=[0.5e-3,0.75e-3,0.8e-3,0.9e-3,1.2e-3,2e-3],psdCumm=[0.,0.2,0.4,0.6,0.8,1.0],num=1080)
#sp.makeCloud(mi,ma,porosity=0.66/1.66,psdSizes=[0.5e-3,0.75e-3,0.8e-3,0.9e-3,1.2e-3,2e-3],psdCumm=[0.,0.2,0.4,0.6,0.8,1.0],num=108000)
# e min =0.56, and e max =0.89, e=0.66
spheres=sp.toSimulation(color=(0,0.5,0.7),material='sphereMat')

num=0
vol=0
mass=0
for b in O.bodies:
    r=b.shape.radius
    x=b.state.pos[0]
    y=b.state.pos[1]
    z=b.state.pos[2]
    ID=b.id
    if x**2+y**2>=(20e-3-r)**2:
        O.bodies.erase(b.id)
    else:
        num+=1
        vol+=4./3.*np.pi*r**3
        mass+=2.69*4./3.*np.pi*r**3
print 'num, e, voidRatio =',num, 0.66, ((80e-3*(40e-3*0.5)**2*np.pi)-vol)/vol

#=======================================fiber=================================================
length=25.4e-3
diameter=0.3e-3
Xw=0.25e-2 # fiber content
fiber_vol=(80e-3*(40e-3*0.5)**2*np.pi)*Xw
Vfiber=length*np.pi*(0.5*diameter)**2
numFiber=fiber_vol/Vfiber
print "the numer of fibres =",int(numFiber)
#numFiber=156
Ne=int(length/(0.9e-3))
print 'Ne=',Ne
nodesIds=[]
cylIds=[]
numFibre=0
target_num=int(numFiber)
target_num=int(8)
####fiberfilewrite
fiberwr=open('fibers.txt','w')
for n in range(500):
    if numFibre<target_num:
  random.seed(n*5)
  z0=random.uniform(0, 80e-3 )
  random.seed( n+2)
  sita0=random.uniform(-np.pi,np.pi)
  random.seed( n+9)
  l_y=random.uniform(0,20e-3)
  x0=l_y*sin(sita0)
  y0=l_y*cos(sita0)

  random.seed(n+8)
  sita=random.uniform(-0.5*np.pi, 0.5*np.pi)
  random.seed(n+9)
  phi=random.uniform(0, 2*np.pi)

  hz=length*sin(sita)
  hx=length*cos(sita)*cos(phi)
  hy=length*cos(sita)*sin(phi)

  x1=x0+hx
  y1=y0+hy
  z1=z0+hz

  if (x1**2+y1**2)<(40e-3*0.5)**2 and 0<(z1)<80e-3 and (x0**2+y0**2)<(40e-3*0.5)**2 and 0<z0<80e-3 :
      fiberwr.write(str(numFibre)+'\t'+str(x0)+'\t'+str(y0)+'\t'+str(z0)+'\t'+str(x1)+'\t'+str(y1)+'\t'+str(z1)+'\n')
      numFibre+=1
      vertices=[]
      for i in range(0, Ne+1):
          px=float(i)*hx/float(Ne)+x0; py=float(i)*hy/float(Ne)+y0; pz=float(i)*hz/float(Ne)+z0;
   vertices.append([px,py,pz])
      cylinderConnection(vertices,0.15e-3,nodesIds,cylIds,color=[0.5,0.5,0],fixed=True,highlight=False,intMaterial='intcylMat',extMaterial='extcylMat')

#=======================================facet=================================================
nw = 24
nh = 60
R=0.00002
color=[255./255.,102./255.,0./255.]
rCyl2 = .5*width / cos(pi/float(nw))
# facets
facets = []
rCyl2 = .5*width / cos(pi/float(nw))
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)), height*(h+0)/float(nh) )
  v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), height*(h+0)/float(nh) )
  v3 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), height*(h+1)/float(nh) )
  v4 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), rCyl2*sin(2*pi*(r+0)/float(nw)), height*(h+1)/float(nh) )
  f1 = facet((v1,v2,v3),color=[255./255.,102./255.,0./255.],material='gridNodeMat')
  f2 = facet((v1,v3,v4),color=[255./255.,102./255.,0./255.],material='gridNodeMat')
  facets.extend((f1,f2))
O.bodies.append(facets)
#===============================================================================

Thanks for any help!!!!!!!