Python - meshio:如何为纯点数据定义顶点单元?

Python - meshio: How to define vertex cells for point-only data?

问题

我有一个二维时间序列数据,我想使用 meshio 将它保存为 XMDF。问题是,我的网格只是具有关联点数据的点数组,而且我没有定义任何单元格。因此,我尝试使用 "vertex" 单元格类型,这是一个单点单元格,但它不起作用。 Meshio 的文档有点缺乏,所以我卡住了。

代码

根据这两个示例 on their Github page,我做了以下操作。我不确定如何正确定义单元格,因为 meshio 没有正确记录。

# generate some data on a 10x10 mesh with 20 time steps (tested, works)
ts = np.arange(20)
x, y = np.meshgrid(np.arange(10), np.arange(10))
data = np.empty((20, 10, 10))
for i, t in enumerate(ts):
    data[i] = np.sin((x + y) * t)
# data is a 3D NumPy array now with dimensions (20,10,10)

# generate list of points (tested, works)
points = [list(p) for p in zip(*(x.flat, y.flat,))]

# won't use cells, so define vertex cell (1 point per cell) <-- ???
cells = [("vertex", [i,]) for i in range(len(points))]

# as seen in meshio's documentation, write time series data
filename = "test.xdmf"
with meshio.xdmf.TimeSeriesWriter(filename) as writer:
    writer.write_points_cells(points, cells)
    for i, t in enumerate(ts):
        writer.write_data(t, point_data={"sin_city": data[i]})

错误

以上脚本产生以下错误:

Traceback (most recent call last):
  File "/home/ezio/Codes/gfield/_temp.py", line 103, in <module>
    writer.write_points_cells(points, cells)
  File "/home/ezio/anaconda3/envs/radpolpy/lib/python3.8/site-packages/meshio/xdmf/time_series.py", line 284, in write_points_cells
    self.points(grid, points)
  File "/home/ezio/anaconda3/envs/radpolpy/lib/python3.8/site-packages/meshio/xdmf/time_series.py", line 340, in points
    if points.shape[1] == 2:
AttributeError: 'list' object has no attribute 'shape'

我尝试了将一些用于 NumPy 数组的数组的不同组合,但我找不到原因。我请求你的帮助。

更新:

将每个使用的数字数组更改为 NumPy 数组后(归功于评论)- 即,在定义 points 之后直接插入 points = np.array(points),并将元胞生成器行更改为 cells = [("vertex", np.array([i,])) for i in range(len(points))] - 我还有一个不同的错误:

Traceback (most recent call last):
  File "/home/ezio/Codes/gfield/_temp.py", line 105, in <module>
    writer.write_points_cells(points, cells)
  File "/home/ezio/anaconda3/envs/radpolpy/lib/python3.8/site-packages/meshio/xdmf/time_series.py", line 285, in write_points_cells
    self.cells(cells, grid)
  File "/home/ezio/anaconda3/envs/radpolpy/lib/python3.8/site-packages/meshio/xdmf/time_series.py", line 409, in cells
    [
  File "/home/ezio/anaconda3/envs/radpolpy/lib/python3.8/site-packages/meshio/xdmf/time_series.py", line 411, in <listcomp>
    np.insert(
  File "<__array_function__ internals>", line 5, in insert
  File "/home/ezio/anaconda3/envs/radpolpy/lib/python3.8/site-packages/numpy/lib/function_base.py", line 4527, in insert
    axis = normalize_axis_index(axis, ndim)
numpy.AxisError: axis 1 is out of bounds for array of dimension 1

(我还注意到文档在示例中没有使用 NumPy 数组。)

问题是:

  1. 我应该在任何地方都使用 NumPy 数组:即使 meshio 的示例使用 Python 列表,xmdf 模块显然无法处理这些;和
  2. 我也应该将数据展平(很明显)。此外,单元格应该以更好的方式定义,尽管最初的概念也有效。

我的代码的工作版本是:

# generate some data on a 10x10 mesh with 20 time steps (tested, works)
ts = np.arange(20)
x, y = np.meshgrid(np.arange(10), np.arange(10))
data = np.empty((20, 10, 10))
for i, t in enumerate(ts):
    data[i] = np.sin((x + y) * t)
# data is a 3D NumPy array now with dimensions (20,10,10)

# generate list of points (tested, works)
points = [list(p) for p in zip(*(x.flat, y.flat,))]
points = np.array(points)  # add this

# won't use cells, so define vertex cell (1 point per cell)
cells = [("vertex", np.array([[i,] for i in range(len(points)])))]
# instead of cells = [("vertex", [i,]) for i in range(len(points))]

# as seen in meshio's documentation, write time series data
filename = "test.xdmf"
with meshio.xdmf.TimeSeriesWriter(filename) as writer:
    writer.write_points_cells(points, cells)
    for i, t in enumerate(ts):
        # here data[i] also should be flattened
        writer.write_data(t, point_data={"sin_city": data[i].flatten})