edem 通过导出颗粒的位置坐标 ID号和接触数获得颗粒团形貌的代码

以下是基于edem导出颗粒位置坐标和接触信息获取颗粒团形貌的示例代码（使用Python语言）：

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

# 导入颗粒位置坐标信息和接触信息
particle_data = np.loadtxt('particle_positions.txt')
contact_data = np.loadtxt('particle_contacts.txt')

# 获取颗粒直径
diameter = 0.05  # 假设颗粒直径为5cm

# 构建颗粒团形貌
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

for i in range(len(particle_data)):
    # 获取颗粒位置坐标
    x, y, z = particle_data[i]
    
    # 根据颗粒直径计算其半径
    radius = diameter / 2
    
    # 构建球体模型
    u, v = np.mgrid[0:2*np.pi:20j, 0:np.pi:10j]
    xx = radius * np.sin(v) * np.cos(u) + x
    yy = radius * np.sin(v) * np.sin(u) + y
    zz = radius * np.cos(v) + z
    
    # 绘制球体
    ax.plot_surface(xx, yy, zz, color='k', alpha=0.2)

# 绘制颗粒接触关系
for i in range(len(contact_data)):
    p1, p2 = contact_data[i]
    x1, y1, z1 = particle_data[int(p1)]
    x2, y2, z2 = particle_data[int(p2)]
    xs = [x1, x2]
    ys = [y1, y2]
    zs = [z1, z2]
    ax.plot(xs, ys, zs, color='k', alpha=0.2)

# 设置坐标轴范围
ax.set_xlim3d([np.min(particle_data[:, 0])-0.1, np.max(particle_data[:, 0])+0.1])
ax.set_ylim3d([np.min(particle_data[:, 1])-0.1, np.max(particle_data[:, 1])+0.1])
ax.set_zlim3d([np.min(particle_data[:, 2])-0.1, np.max(particle_data[:, 2])+0.1])

# 显示图形
plt.show()

请注意，此代码仅是示例代码，并且需要根据实际数据格式和颗粒直径进行适当修改。