离散元后处理与Python实践——ballarray函数
- 2026-04-16 03:37:40
摘要
很多人在PFC里还在用FISH逐颗粒循环，但其实官方早就提供了更高效的Python接口。最近在 Itasca documentation 中偶然发现 ballarray，可以直接将颗粒属性转为NumPy数组，从而实现高效的向量化后处理与筛选分析。
1 ballarray函数解读
ballarray并不是新引入的功能，而是长期被低估的内置能力，在PFC 6.0版本中已经存在。itasca.ballarray 位于 Itasca 官方文档的 Scripting / PFC Python / PFC Python API 模块中，但由于其在文档结构中层级较深且缺乏显式教程说明，该接口在实际使用中并不常见，容易被忽略。然而它提供了对颗粒系统状态变量的批量 NumPy 化访问能力。
import itascaitasca.ballarray# Array interface for Itasca balls下面看几个重要的函数：
球体的半径
itasca.ballarray.radius() '''→ array float{ball}.Get a numpy array of the ball radii.'''球体的位置信息
itasca.ballarray.pos() '''→ array float{ball, 3}.Get a numpy array of the ball centroid location.'''球体的重要物理量或标签
itasca.ballarray.disp() itasca.ballarray.vel()itasca.ballarray.id()'''→ array float{ball, 3}.Get a numpy array of the ball displacement, velocity and ID.'''以前在 PFC 里想拿颗粒数据，往往要写 ball.list() 再配合循环逐个处理，不仅麻烦还不够高效。而现在 ballarray 只需要一行代码就能直接获取全部颗粒信息。更重要的是，它返回的是标准 NumPy 数组，这意味着你可以直接使用 NumPy 的索引、筛选和数学运算能力，把后处理效率提升一个数量级。
2 ballarray函数的应用场景
借助 ballarray 返回的 NumPy 数组，可以直接调用 Python 科学计算生态中的各种工具，从而将原本繁琐的颗粒后处理转化为高效的数组运算。基于这一点，可以很方便地实现但不限于如下几类典型分析场景。
2.1 快速筛选所需颗粒
import itasca as itimport numpy as nppos = itasca.ballarray.pos() # (N,3)rad = itasca.ballarray.radius() # (N,)disp = itasca.ballarray.disp() # (N,3)bid = itasca.ballarray.id() # (N,)vel = itasca.ballarray.vel() # (N,3)# 例1 不同区域的粒径分布region1 = pos[:,0] < 0.5#x<0.5r_region1 = rad[region1]hist, bins = np.histogram(r_region1, bins=20)# 例子2：筛选某一区域的颗粒xmin, xmax = 0.0, 1.0ymin, ymax = 0.0, 1.0zmin, zmax = 0.0, 1.0mask_region = ( (pos[:,0] >= xmin) & (pos[:,0] <= xmax) & (pos[:,1] >= ymin) & (pos[:,1] <= ymax) & (pos[:,2] >= zmin) & (pos[:,2] <= zmax))selected_ids = bid[mask_region]#例子3：筛选位移大于某个临界值的所有球体disp_mag = np.linalg.norm(disp, axis=1)threshold = 1e-4mask_disp = disp_mag > thresholdselected_ids = bid[mask_disp]2.2 动态跟踪
例如在渗流问题中，如果想监测颗粒的移动轨迹，则可以在不同时间内存储pos信息，后续可以根据id信息挑选自己关心的颗粒运移轨迹。
pos_t1 = itasca.ballarray.pos()pos_t2 = itasca.ballarray.pos()2.3 后处理
可视化的vtk写法之前曾经手动写过，可以参照 利用Paraview绘制颗粒信息 。但是如果现在利用ballarray可以更加方便的实现。
VTK输出函数
#球体可视化vtk输出函数,利用pyvista简化import numpy as npimport pyvista as pvimport itasca as itdefWrite_VTK(filename, bin = False):''' filename: string, 文件名 bin： Boolean 控制是否写入二进制文件，False写入Asic II '''# 获取球体数据 points = it.ballarray.pos() balls = pv.PolyData(points) balls["radius"] = it.ballarray.radius() balls["displacement"] = it.ballarray.disp() balls["velocity"] = it.ballarray.vel() balls["ids"] = it.ballarray.ids() balls.save(f"balls_{filename}.vtk", binary=bin) print(f"VTK saved: {filename}.vtk") 调用(PFC IPython环境) PFC IPython环境运行参考为什么在PFC中我总是推荐利用import运行Python，离散元后处理与Python实践——代码运行及勘误，库管理参考PFC中如何管理第三方Python库。
import itasca as itimport Write_VTKimport pyvista as pvit.command('model rest "calm.sav"')Write_VTK('test')可视化(python环境或Paraview) 可视化请参考OpenSeesPy入门（7）- 利用pyvista库可视化vtk文件。
# 可视化import pyvista as pvimport numpy as npmesh = pv.read("test.vtk")# 创建一个单位球作为 glyph 的模板sphere = pv.Sphere(radius=1.0, theta_resolution=16, phi_resolution=16)# 用 glyph 按 radius 缩放glyphs = mesh.glyph( geom=sphere, scale="radius", # 使用 point_data 中的 radius orient=False, # 不需要方向就关掉)plotter = pv.Plotter()plotter.add_mesh( glyphs, color="orange", show_edges=False,)plotter.show()最终可以得到类似的图
以上示例更多是一些基础用法，主要起到抛砖引玉的作用。实际上，一旦将颗粒数据转化为 NumPy 数组，其可扩展性会大大提升。如果读者对 Python，尤其是 NumPy 比较熟悉，建议进一步查阅参考文献中的相关内容，结合自身问题深入探索，往往会有更多意想不到的收获。
此外，颗粒之间的接触信息还可以通过 ballballarray 接口获取，其数据维度相比 ballarray 更为丰富，也为开展更复杂的分析提供了可能。有兴趣的话不妨亲自尝试一下，或许会发现更多有价值的用法。
参考文献
【1】https://docs.itascasoftware.com/itasca900/common/docproject/source/manual/scripting/python/doc/itasca.ballarray.html
