Hi everyone,
I have particles generated from LEoPart that are advected from my velocity solution. While I know how to output the solution velocity vector (ex: u.vector().get_local()), how can I output the velocity on a per-particle basis? When looking at all of the possible attributes for the particle class, I don’t see anything that can give me the velocity. Since the particles are being advected in the flow, I imagine that there is a way to retrieve this information.
The key here is that particles in LEoPart only carry around the data you tell them to. Particles do not actually have a velocity, but in fact interpolate a velocity field at their respective positions. So if you wanted to find “particle velocities” you need to perform this interpolation into an appropriately sized array associated with each particle. See the following example:
import numpy as np
from dolfin import *
import leopart
# Setup mesh and velocity function
mesh = UnitSquareMesh(4, 4)
V = VectorFunctionSpace(mesh, "CG", 1)
u = Function(V)
u.interpolate(Expression(("-x[1]", "x[0]"), degree=1))
# Setup particles
pres = 100
x = leopart.RegularRectangle(Point(0, 0), Point(1, 1)).generate((10, 10))
u_ptcl_data = np.zeros((len(x), mesh.geometry().dim()), dtype=np.float_)
# Interpolate velocity function onto particles, index slot 1
property_idx = 1
ptcls = leopart.particles(x, [u_ptcl_data], mesh)
ptcls.interpolate(u, property_idx)
# Get the flat data and reshape it to dimensional form
u_ptcl_data = np.array(ptcls.get_property(property_idx)).reshape((-1, mesh.geometry().dim()))
# particle positions are always property zero
positions = np.array(ptcls.get_property(0)).reshape((-1, mesh.geometry().dim()))
# Plot
import matplotlib.pyplot as plt
plot(mesh)
plt.quiver(positions[:,0], positions[:,1], u_ptcl_data[:,0], u_ptcl_data[:,1])
plt.scatter(positions[:,0], positions[:,1])
plt.axis("equal")
plt.show()
Which yields
Thank you Nate, this is exactly what I was looking for! The example really helped as well.