# yade-users team mailing list archive

## Re: [Question #705741]: Aborted (core dumped) when using VTKRecorder

This is single script of the original posted problem.

It give following error when using "VTKRecorder"

In [1]: python3.10: ./pkg/common/InsertionSortCollider.hpp:125: yade::InsertionSortCollider::Bounds& yade::InsertionSortCollider::VecBounds::operator[](long int): Assertion `idx < long(size()) && idx >= 0' failed.
Aborted (core dumped)

=======

width = 0.05
height = 0.03
margin = 30
# Calculate extra length
dx = width/100/2*margin
dy = width/100/2*margin
dz = height/100/2*margin

# side pannel
p1s = (-width/2,-width/2-dy,-height/2-dz)
p5s = (-width/2,-width/2-dy,height/2+dz)
p6s = (-width/2,width/2+dy,height/2+dz)
p2s = (-width/2,width/2+dy,-height/2-dz)

side1_1 = utils.facet(vertices=[p1s,p5s,p2s], wire=True, highlight=False) #p1 p5 p2
side1_2 = utils.facet(vertices=[p2s,p5s,p6s], wire=True, highlight=False) #p1 p5 p2
O.bodies.append(side1_1)
O.bodies.append(side1_2)

p4s = (width/2,-width/2-dy,-height/2-dz)
p8s = (width/2,-width/2-dy,height/2+dz)
p7s = (width/2,width/2+dy,height/2+dz)
p3s = (width/2,width/2+dy,-height/2-dz)

side2_1 = utils.facet(vertices=[p4s,p8s,p3s], wire=True, highlight=False) #p1 p5 p2
side2_2 = utils.facet(vertices=[p3s,p8s,p7s], wire=True, highlight=False) #p1 p5 p2
O.bodies.append(side2_1)
O.bodies.append(side2_2)

# front pannel
p1f = (-width/2-dx,-width/2,-height/2-dz)
p5f = (-width/2-dx,-width/2,height/2+dz)
p4f = (width/2+dx,-width/2,-height/2-dz)
p8f = (width/2+dx,-width/2,height/2+dz)

front1 = utils.facet(vertices=[p1f,p5f,p8f], wire=True, highlight=False)
front2 = utils.facet(vertices=[p1f,p8f,p4f], wire=True, highlight=False)
O.bodies.append(front1)
O.bodies.append(front2)

# back pannel
p2b = (-width/2-dx,width/2,-height/2-dz)
p6b = (-width/2-dx,width/2,height/2+dz)
p3b = (width/2+dx,width/2,-height/2-dz)
p7b = (width/2+dx,width/2,height/2+dz)

back1 = utils.facet(vertices=[p2b,p6b,p7b], wire=True, highlight=False)
back2 = utils.facet(vertices=[p2b,p7b,p3b], wire=True, highlight=False)
O.bodies.append(back1)
O.bodies.append(back2)

#bottom pannel
p1bt = (-width/2-dx,-width/2-dy,-height/2)
p2bt = (-width/2-dx,width/2+dy,-height/2)
p3bt = (width/2+dx,width/2+dy,-height/2)
p4bt = (width/2+dx,-width/2-dy,-height/2)

bot1 = utils.facet(vertices=[p1bt,p2bt,p3bt], wire=True, highlight=False)
bot2 = utils.facet(vertices=[p1bt,p4bt,p3bt], wire=True, highlight=False)
O.bodies.append(bot1)
O.bodies.append(bot2)

O.materials.append(FrictMat(young=20e6, poisson=0.17, density=2700, frictionAngle=0.523, label='frictmat'))
#O.materials.append(FrictMat(young=70e9, poisson=0.17, density=2300, frictionAngle=0.523, label='frictmat')) # silica https://www.azom.com/properties.aspx?ArticleID=1114

cp1 = (-width/2,-width/2,3*height) #corner point #1
cp2 = (width/2,width/2,2*height+3*height/2) #corner point #2

sp = pack.SpherePack()
sp.makeCloud(cp1, cp2, rMean=radius_mean, rRelFuzz=0.0, num = 10000)
sp.toSimulation(color=(0,0,1)) # pure blue

O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],  # collision geometry
[Ip2_FrictMat_FrictMat_FrictPhys()],  # collision "physics"
[Law2_ScGeom_FrictPhys_CundallStrack()]  # contact law -- apply forces
),
# Apply gravity force to particles. damping: numerical dissipation of energy.
#    GlobalStiffnessTimeStepper(active=True, timestepSafetyCoefficient=0.8,  timeStepUpdateInterval=1000, label = 'timeStepper'),
NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.5),
VTKRecorder(fileName='spheres/3d-vtk-', recorders=['spheres', 'intr', 'facets', 'colors'], iterPeriod=100),
PyRunner(command='checkUnbalanced()', iterPeriod=1000),
#    DomainLimiter(lo=(-width,-width,-height), hi=(width,width,5*height), iterPeriod = 10000, label = 'Domain') # destroy balls outside domain in every 100 steps
]

for b in O.bodies:
b.shape.color=scalarOnColorScale(b.state.vel.norm(),0.0,1e-1)

elapsed_time=0.0
def checkUnbalanced():
#     print(unbalancedForce())
global elapsed_time
if (O.time-elapsed_time) > 1.0:
ball_highest_z = numpy.max([b.state.pos[2] for b in O.bodies if isinstance(b.shape,Sphere)]) # make list of ball z position if body is sphere and find max value
print("The highest ball position is = ", ball_highest_z)
if ball_highest_z < height/2:
sp.makeCloud(cp1, cp2, rMean=radius_mean, rRelFuzz=0.0, num = 10000)
sp.toSimulation(color=(0,0,1)) # pure blue
print("Total number of balls =", len(O.bodies))
elapsed_time=O.time
else:
O.pause(); print("simulation paused")

# plotting