yade-users team mailing list archive

[JINSUN LEE <question705741@xxxxxxxxxxxxxxxxxxxxx>]

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

JINSUN LEE gave more information on the question:
 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)

=======

from yade import pack, plot


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
radius_mean = 0.002

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='addData()', iterPeriod=1000),
    PyRunner(command='checkUnbalanced()', iterPeriod=1000),
    PyRunner(command='addPlotData()', iterPeriod=1000),
#    DomainLimiter(lo=(-width,-width,-height), hi=(width,width,5*height), iterPeriod = 10000, label = 'Domain') # destroy balls outside domain in every 100 steps
]

def addData():
     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")
             O.save("CDSS_ball_2.yade")

# plotting
def addPlotData():
    plot.addData(i=O.iter, unbalanced=unbalancedForce())

plot.plots={'i':('unbalanced')}    
plot.plot()

O.dt=0.5*PWaveTimeStep()
O.run() #; O.wait()

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