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

Problem

I have a 2D time-series data, and I want to save it as XMDF using meshio. The problem is, my mesh is just an array of points with associated point data, and I don’t have any cell defined. As such, I tried to use the "vertex" cell type, which is a single-point cell, but it doesn’t work. Meshio’s documentation is kind of lacking, so I’m stuck.

Code

Following the two examples on their Github page, I did the following. I’m not sure how to define the cells correctly, as meshio doesn’t document this properly.

# 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]})

Error

The above script produces the following error:

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'

I tried different combinations of converting some of the arrays used to NumPy array, but I couldn’t find out the cause. I ask for your help.

Update:

After changing every used number array to NumPy arrays (credit to comments) – that is, inserting points = np.array(points) directly after points is defined, and changing the cell generator line to cells = [("vertex", np.array([i,])) for i in range(len(points))] – I still have a different error:

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

(I also note that the documentation does not use NumPy arrays in the examples.)

Answer

The problem was that:

1. I should’ve use NumPy arrays everywhere: even though meshio’s example use Python lists, the xmdf module apparently can’t handle those; and
2. I should’ve flatten the data as well (obviously). Also the cells should be be defined in a better way, though the original concept works too.

A working version of my code is:

# 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})