Issue with plotting
I hope someone can help me. I have a well working simulation but I am trying to plot and I keep getting an empty graph. The simulation is of a drill bit (cone) penetrating substrate (packing of spheres in a cylinder) and I need to plot the position of the cone along z axis (penetration depth) against the normal force exerted on the cone (by the substrate/spheres). Here is the code:
from yade import pack, plot, qt, ymport
import numpy as np
cone = yade.geom.facetCone((0,0,7), radiusTop=0.5, radiusBottom=0, height=2, wire = False,color = (1,0,0),
orientation=Quaternion((1, 0, 0), 0), segmentsNumber=10, wallMask=7, angleRange=None, closeGap=False, radiusTopInner=-1, radiusBottomInner=-1)
cone_id = [] # Initialize an empty list to store cone body IDs
for i in range(len(cone)): # Assuming 'cone' is the list of cone bodies
cone_id.append(O.bodies.append(cone[i]))
O.bodies.append(
yade.geom.facetCylinder(center=(0,0,0), radius=1.5, wire=True, color=(0,1,0),
height=4, orientation=Quaternion((1, 0, 0), 0),
segmentsNumber=10, wallMask=6, angleRange=None,
closeGap=False, radiusTopInner=-1.4, radiusBottomInner=-1.4)
)
sp = pack.SpherePack()
# generate randomly spheres with uniform radius distribution
sp.makeCloud((-1, -1, -2), (1, 1, 4), rMean=0.2, rRelFuzz=0)
# add the sphere pack to the simulation
sp.toSimulation()
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
InteractionLoop(
# handle sphere+sphere and facet+sphere collisions
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
# [Ip2_FrictMat_FrictMat_MindlinPhys()],
# [Law2_ScGeom_MindlinPhys_Mindlin()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.04),
# increase velocity to observe if drillbit is moving
TranslationEngine(translationAxis=(0, 0, -1),
velocity=0.1, ids=cone_id, label='translation',
dead=True),
PyRunner(command= 'move()',iterPeriod=1),
# call the checkUnbalanced function (defined below) every 2 seconds
#PyRunner(command='checkUnbalanced()', realPeriod=2),
# call the addPlotData function every 200 steps
#PyRunner(command='addPlotData()', iterPeriod=100),
]
def move():
if O.iter>2040: # at the 500th iteration, start moving the drillbit
translation.dead= False
if O.bodies[2].state.pos[2]<=1:
translation.velocity = 0 # to stop the drillbit
#translation.velocity = -1 # to reverse the drillbit
def addPlotData():
# Get the position of the cone along the z-axis
z_position = [O.bodies[body_id].state.pos[2] for body_id in cone_id]
# Calculate the normal force exerted on the cone
normal_force = []
for i in range(len(cone_id)):
# Accumulate normal forces acting on the cone from all contact interactions
force_sum = Vector3(0, 0, 0)
for ci in O.bodies[cone_id[i]].contacts:
if isinstance(ci.shape1, Facet) and isinstance(ci.shape2, Sphere):
# Check if one of the shapes is the cone and the other is a sphere
force_sum += ci.force
elif isinstance(ci.shape2, Facet) and isinstance(ci.shape1, Sphere):
# Check if the order of shapes is reversed
force_sum -= ci.force
normal_force.append(force_sum.norm())
# Add data to the plot
plot.addData(
z_position=z_position,
normal_force=normal_force
)
O.dt = .8 * PWaveTimeStep()
O.run()
qt.View()
plot.plots = {
'i': ('z_position', 'normal_force'),
}
plot.plot()
Edited by Jan Stránský