Wrong cell_type in cell_data_dict

[Jesusbill]

Reading the med file attached and looking at cells and cell_data_dict one gets the wrong cell_type keys in the latter for the RES_____FERRAILLAGE field. I list below the length of the data for each cell type:

"cells": {
    "vertex": 2,
    "quad": 45,
    "line": 45,
    "triangle": 4
},
"cell_data": {
    "RES_____FERRAILLAGE": {
        "vertex": 45,
        "quad": 4
    },
    "cell_tags": {
        "vertex": 2,
        "quad": 45,
        "line": 45,
        "triangle": 4
    }
}

Basically, for RES_____FERRAILLAGE in this case vertex should be quad and quad should be triangle.

The issue here is that:

1. This field that applies to each cell (one value per cell - type is ELEM) is defined only for some of the cell types, not all.
2. The construction of cell_data_dict here is based on cell_data which essentially does not contain any information of the cell type and it is assumed that every field contains all cell types by default.

A similar issue (#383), although different in that there were missing values within a given cell type, was tackled in #384 by adding None values. Since cell_data are essentially unaware of the cell type I wonder if a solution could be to add any missing cell type within a field with None values as well. This would allow to fix the issue modifying only the _med.py file, if this issue has not been observed for other formats as well.

output.med.zip

[Jesusbill]

Apologies for any confusion, I was using v4.0.12 to run my code, whereas with the current version of v4.2.1 converting such files throws an assertion error in the initialization of Mesh here. This issue is also apparently related to #876 and #925. @nschloe and @tianyikillua: eager to work on this issue if we can agree on the most sensible solution

[tianyikillua]

Does #927 fix your concerns?

[Jesusbill]

Almost. It gives the following error when trying to access cell_data_dict

File "/home/jesusbill/miniconda3/envs/ifc/lib/python3.8/site-packages/meshio/_mesh.py", line 164, in cell_data_dict
    cell_data_dict[key][cell_type] = numpy.concatenate(val)
  File "<__array_function__ internals>", line 5, in concatenate
ValueError: zero-dimensional arrays cannot be concatenated

Modifying this line to have None in a list itself does the trick and returns the correctly attributed data to cell type.

cell_data[name] = [[None]] * len(cell_types)

Note, however, that trying to convert the med file to another format still returns an error. Does not affect me directly at the moment, though, probably others

 File "/home/jesusbill/miniconda3/envs/ifc/bin/meshio-convert", line 8, in <module>
    sys.exit(convert())
  File "/home/jesusbill/miniconda3/envs/ifc/lib/python3.8/site-packages/meshio/_cli/_convert.py", line 43, in convert
    write(args.outfile, mesh, **kwargs)
  File "/home/jesusbill/miniconda3/envs/ifc/lib/python3.8/site-packages/meshio/_helpers.py", line 144, in write
    return writer(filename, mesh, **kwargs)
  File "/home/jesusbill/miniconda3/envs/ifc/lib/python3.8/site-packages/meshio/vtu/_vtu.py", line 574, in write
    data[k] = dat.astype(dat.dtype.newbyteorder("="), copy=False)
AttributeError: 'list' object has no attribute 'astype'

Thanks @tianyikillua for looking into this and your prompt response

[tianyikillua]

cell_data[name] = [numpy.array([])] * len(cell_types)

should do the trick.

However the problem is that VTU assumes we have a cell data for each cell type.

[tianyikillua]

@nschloe

What's the standard way for cell data not defined for a specific cell data? Use None?

cells = [("quad", ...), ("hexahedron", ...)]
cell_data = [np.array([...]), None]

[nschloe]

@tianyikillua Good question. None seems reasonable, go for it.

[nschloe]

Is this fixed?

[gdmcbain]

No, current master still raises the ValueError of 2020-09-26.

[Jesusbill]

I have been looking at this issue today. Looks like a tough case to resolve ... With the last modification from @tianyikillua back in September the first part of the issue was handled but the functionality is still not there, as these files are not currently readable by meshio. Looks like neither None or np.array([]) are good choices for non existent cell_data. The latter could be a partial resolution as it would allow at least to be able to read the med file and even write it back in med, but it does fail to write in other formats like vtu.

I did a test, however, that I believe can guide to an acceptable solution, a solution that allows to read and write such files in all formats. This solution is similar to what was followed for #384 for non existent data outside a subset. And that would be to construct the field for the given cell type with NaN values.

This is the snippet that I used to impose this solution (first I changed None with np.array([]) in order to be able to read the file)

import meshio
import json
import numpy as np

file = meshio.read("output.med", file_format="med")

# Imposing NaN values (field of 7 components) for RES_____FERRAILLAGE for vertex (2 cells) and line (45 cells)
new_cell_data = {}
for key1 in file.cell_data_dict.keys():
    if key1 == "RES_____FERRAILLAGE":
        cell_data = file.cell_data_dict[key1]
        new_cell_data[key1] = []
        for key2 in cell_data.keys():
            if key2 == "vertex":
                new_cell_data[key1].append(np.array([[np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan] for _ in range(2)]))
            elif key2 == "line":
                new_cell_data[key1].append(np.array([[np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan] for _ in range(45)]))
            else:
                new_cell_data[key1].append(cell_data[key2])

meshio.Mesh(
    file.points,
    file.cells_dict,
    # Optionally provide extra data on points, cells, etc.
    point_data=file.point_data,
    cell_data=new_cell_data,
    field_data=file.field_data
).write(
    "output_new.med",
    file_format="med"
)

meshio.Mesh(
    file.points,
    file.cells_dict,
    # Optionally provide extra data on points, cells, etc.
    point_data=file.point_data,
    cell_data=new_cell_data,
    # field_data=file.field_data
).write(
    "output_new.vtu",
    file_format="vtu"
)

What do you think? Should I work on a possible fix in this direction?

[nschloe]

Looks like neither None or np.array([]) are good choices for non existent cell_data

Why not None? Seems natural enough to me.

[Jesusbill]

Yeah, I only meant that it breaks the read and write functionality. I also did some tests keeping None and debugging the read/write functionality. It would need to modify _mesh.py and every i/o file for each format, whereas I was looking for a solution that would be handled by working only on the med format. In any case, with None and some minor modifications to handle some exceptions, I get the following error when exporting to vtk or vtu

Traceback (most recent call last):
  File "io.py", line 48, in <module>
    meshio.Mesh(
  File "/home/jesusbill/Dev-Projects/github.com/Jesusbill/meshio/meshio/_mesh.py", line 182, in write
    write(path_or_buf, self, file_format, **kwargs)
  File "/home/jesusbill/Dev-Projects/github.com/Jesusbill/meshio/meshio/_helpers.py", line 147, in write
    return writer(filename, mesh, **kwargs)
  File "/home/jesusbill/Dev-Projects/github.com/Jesusbill/meshio/meshio/vtk/_vtk.py", line 695, in write
    _write_field_data(f, mesh.cell_data, binary)
  File "/home/jesusbill/Dev-Projects/github.com/Jesusbill/meshio/meshio/vtk/_vtk.py", line 773, in _write_field_data
    values = np.concatenate(values)
  File "<__array_function__ internals>", line 5, in concatenate
ValueError: all the input arrays must have same number of dimensions, but the array at index 0 has 2 dimension(s) and the array at index 1 has 3 dimension(s)

I agree that None makes sense but I don't know if it can be eventually applicable universally.

But, if you think this is the way to go, I can give it a try, just to see what modifications it would require, if I can make it work of course

[nschloe]

Yeah, I only meant that it breaks the read and write functionality.

I don't get it. Why would that be the case?

It would need to modify _mesh.py and every i/o file for each format, whereas I was looking for a solution that would be handled by working only on the med format.

Ah okay. Well, that can't be expected. You'll probably have to modify most mesh i/o routines, but only one or two lines of code each. I'd consider that doable.

I agree that None makes sense but I don't know if it can be eventually applicable universally.

Just do ahead and try it out, I imagine it won't be that hard.

[Jesusbill]

Just do ahead and try it out, I imagine it won't be that hard.

Ok! Will get onto it. Thanks for the feedback