Add support to meshio readers and writers (pyvista#472)

* add support to meshio readers * remove useless tolist * raise error if meshio not installed * fix minor issues * add avs ucd reader * include read_meshio in read * remove meshio import test * make read_meshio public * keep error thrown by meshio * add meshio dependency * use pypi link * fix offset * add save_meshio * fix bug due to wrong input type * use pyvista properties instead of vtk * Update pyvista/utilities/fileio.py Co-Authored-By: Bane Sullivan <banesullivan@gmail.com> * Update pyvista/utilities/fileio.py Co-Authored-By: Bane Sullivan <banesullivan@gmail.com> * Update pyvista/utilities/fileio.py Co-Authored-By: Bane Sullivan <banesullivan@gmail.com> * Update pyvista/utilities/fileio.py Co-Authored-By: Bane Sullivan <banesullivan@gmail.com> * Update pyvista/utilities/fileio.py Co-Authored-By: Bane Sullivan <banesullivan@gmail.com> * convert pixels to quads * add comments and remove duplicate os * add test for meshio reader and writer * handle white spaces in data name * improve meshio test * add uniform to meshio test * add meshio file format test * make sure filename is a string * add save_meshio
hippo91 · Dec 9, 2019 · 272bf7c · 272bf7c
1 parent 4ff55a4
commit 272bf7c
Show file tree

Hide file tree

Showing 6 changed files with 184 additions and 2 deletions.
diff --git a/docs/getting-started/installation.rst b/docs/getting-started/installation.rst
@@ -18,6 +18,7 @@ the following projects are required dependencies of PyVista:
 * `numpy <https://pypi.org/project/numpy/>`_ - NumPy arrays provide a core foundation for PyVista's data array access.
 * `imageio <https://pypi.org/project/imageio/>`_ - This library is used for saving screenshots.
 * `appdirs <https://pypi.org/project/appdirs/>`_ - Data management for our example datasets so users can download tutorials on the fly.
+* `meshio <https://pypi.org/project/meshio/>`_ - Input/Output for many mesh formats.
 
 PyPI
 ~~~~

diff --git a/docs/utilities/utilities.rst b/docs/utilities/utilities.rst
@@ -33,6 +33,8 @@ File IO
 
 .. autofunction:: pyvista.read_legacy
 
+.. autofunction:: pyvista.save_meshio
+
 
 Mesh Creation
 ~~~~~~~~~~~~~

diff --git a/pyvista/utilities/fileio.py b/pyvista/utilities/fileio.py
@@ -8,6 +8,10 @@
 
 import imageio
 
+import numpy
+
+import meshio
+
 READERS = {
     # Standard dataset readers:
     '.vtk': vtk.vtkDataSetReader,
@@ -58,6 +62,7 @@
     #TODO: '.vpc': vtk.vtkVPIC?????,
     # '.bin': vtk.vtkMultiBlockPLOT3DReader,# TODO: non-default routine
     '.tri': vtk.vtkMCubesReader,
+    '.inp': vtk.vtkAVSucdReader,
 }
 
 VTK_MAJOR = vtk.vtkVersion().GetVTKMajorVersion()
@@ -140,7 +145,7 @@ def read_legacy(filename):
     return pyvista.wrap(output)
 
 
-def read(filename, attrs=None):
+def read(filename, attrs=None, file_format=None):
     """Read any VTK file.
 
     It will figure out what reader to use then wrap the VTK object for
@@ -153,13 +158,19 @@ def read(filename, attrs=None):
         the attribute/method names and values are the arguments passed to those
         calls. If you do not have any attributes to call, pass ``None`` as the
         value.
+    file_format : str, optional
+        Format of file to read with meshio.
 
     """
     filename = os.path.abspath(os.path.expanduser(filename))
     if not os.path.isfile(filename):
         raise IOError('File ({}) not found'.format(filename))
     ext = get_ext(filename)
 
+    # Read file using meshio.read if file_format is present
+    if file_format:
+        return read_meshio(filename, file_format)
+
     # From the extension, decide which reader to use
     if attrs is not None:
         reader = get_reader(filename)
@@ -187,7 +198,12 @@ def read(filename, attrs=None):
             reader = get_reader(filename)
             return standard_reader_routine(reader, filename)
         except KeyError:
-            pass
+            # Attempt read with meshio
+            try:
+                return read_meshio(filename)
+            except AssertionError:
+                pass
+
     raise IOError("This file was not able to be automatically read by pyvista.")
 
 
@@ -233,3 +249,126 @@ def read_exodus(filename,
 
     reader.Update()
     return pyvista.wrap(reader.GetOutput())
+
+
+def read_meshio(filename, file_format = None):
+    """Read any mesh file using meshio."""
+    from meshio._vtk import (
+        meshio_to_vtk_type,
+        vtk_type_to_numnodes,
+    )
+
+    # Make sure relative paths will work
+    filename = os.path.abspath(os.path.expanduser(str(filename)))
+
+    # Read mesh file
+    mesh = meshio.read(filename, file_format)
+
+    # Extract cells from meshio.Mesh object
+    offset = []
+    cells = []
+    cell_type = []
+    cell_data = {}
+    next_offset = 0
+    for k, v in mesh.cells.items():
+        vtk_type = meshio_to_vtk_type[k]
+        numnodes = vtk_type_to_numnodes[vtk_type]
+        offset += [next_offset+i*(numnodes+1) for i in range(len(v))]
+        cells.append(numpy.hstack((numpy.full((len(v), 1), numnodes), v)).ravel())
+        cell_type += [vtk_type] * len(v)
+        next_offset = offset[-1] + numnodes + 1
+
+        # Extract cell data
+        if k in mesh.cell_data.keys():
+            for kk, vv in mesh.cell_data[k].items():
+                if kk in cell_data:
+                    cell_data[kk] = numpy.concatenate((cell_data[kk], vv))
+                else:
+                    cell_data[kk] = vv
+
+    # Create pyvista.UnstructuredGrid object
+    grid = pyvista.UnstructuredGrid(
+        numpy.array(offset),
+        numpy.concatenate(cells),
+        numpy.array(cell_type),
+        numpy.array(mesh.points, numpy.float64),
+    )
+
+    # Set point data
+    for k, v in mesh.point_data.items():
+        data = vtk.util.numpy_support.numpy_to_vtk(numpy.array(v, numpy.float64))
+        data.SetName(k)
+        grid.GetPointData().AddArray(data)
+
+    # Set cell data
+    for k, v in cell_data.items():
+        data = vtk.util.numpy_support.numpy_to_vtk(numpy.array(v, numpy.float64))
+        data.SetName(k)
+        grid.GetCellData().AddArray(data)
+    return grid
+
+
+def save_meshio(filename, mesh, file_format = None, **kwargs):
+    """Save mesh to file using meshio.
+
+    Parameters
+    ----------
+    mesh : pyvista.Common
+        Any PyVista mesh/spatial data type.
+    file_format : str
+        File type for meshio to save.
+
+    """
+    from meshio._vtk import vtk_to_meshio_type
+
+    # Make sure relative paths will work
+    filename = os.path.abspath(os.path.expanduser(str(filename)))
+
+    # Cast to pyvista.UnstructuredGrid
+    if not isinstance(mesh, pyvista.UnstructuredGrid):
+        mesh = mesh.cast_to_unstructured_grid()
+
+    # Copy useful arrays to avoid repeated calls to properties
+    vtk_offset = mesh.offset
+    vtk_cells = mesh.cells
+    vtk_cell_type = mesh.celltypes
+
+    # Get cells
+    cells = {k: [] for k in numpy.unique(vtk_cell_type)}
+    if 8 in cells.keys():
+        cells[9] = cells.pop(8)                 # Handle pixels
+    if 11 in cells.keys():
+        cells[12] = cells.pop(11)               # Handle voxels
+    mapper = {k: [] for k in cells.keys()}      # For cell data
+    for i, (offset, cell_type) in enumerate(zip(vtk_offset, vtk_cell_type)):
+        numnodes = vtk_cells[offset]
+        cell = vtk_cells[offset+1:offset+1+numnodes]
+        cell = cell if cell_type not in {8, 11} \
+            else cell[[ 0, 1, 3, 2 ]] if cell_type == 8 \
+            else cell[[ 0, 1, 3, 2, 4, 5, 7, 6 ]]
+        cell_type = cell_type if cell_type not in {8, 11} else cell_type+1
+        cells[cell_type].append(cell)
+        mapper[cell_type].append(i)
+    cells = {vtk_to_meshio_type[k]: numpy.vstack(v) for k, v in cells.items()}
+    mapper = {vtk_to_meshio_type[k]: v for k, v in mapper.items()}
+
+    # Get point data
+    point_data = {k.replace(" ", "_"): v for k, v in mesh.point_arrays.items()}
+
+    # Get cell data
+    vtk_cell_data = mesh.cell_arrays
+    cell_data = {
+        k: {kk.replace(" ", "_"): vv[v] for kk, vv in vtk_cell_data.items()}
+        for k, v in mapper.items()
+    } if vtk_cell_data else {}
+
+    # Save using meshio
+    meshio.write_points_cells(
+        filename = filename,
+        points = numpy.array(mesh.points),
+        cells = cells,
+        point_data = point_data,
+        cell_data = cell_data,
+        file_format = file_format,
+        **kwargs
+    )
diff --git a/requirements.txt b/requirements.txt
@@ -14,3 +14,4 @@ imageio-ffmpeg
 colorcet
 cmocean
 scooby>=0.2.2
+meshio
diff --git a/setup.py b/setup.py
@@ -22,6 +22,7 @@
                     'imageio',
                     'appdirs',
                     'scooby>=0.2.2',
+                    'meshio',
                     ]
 
 # add vtk if not windows and 2.7

diff --git a/tests/test_meshio.py b/tests/test_meshio.py
@@ -0,0 +1,38 @@
+import pytest
+
+import numpy
+
+import pyvista
+
+from pyvista import examples
+
+
+beam = pyvista.UnstructuredGrid(examples.hexbeamfile)
+airplane = examples.load_airplane().cast_to_unstructured_grid()
+uniform = examples.load_uniform().cast_to_unstructured_grid()
+@pytest.mark.parametrize("mesh_in", [beam, airplane, uniform])
+def test_meshio(mesh_in, tmpdir):
+    # Save and read reference mesh using meshio
+    filename = tmpdir.mkdir("tmpdir").join("test_mesh.vtk")
+    pyvista.save_meshio(filename, mesh_in)
+    mesh = pyvista.read_meshio(filename)
+
+    # Assert mesh is still the same
+    assert numpy.allclose(mesh_in.points, mesh.points)
+    if (mesh_in.celltypes == 11).all():
+        cells = mesh_in.cells.reshape((mesh_in.n_cells, 9))[:,[0,1,2,4,3,5,6,8,7]].ravel()
+        assert numpy.allclose(cells, mesh.cells)
+    else:
+        assert numpy.allclose(mesh_in.cells, mesh.cells)
+    for k, v in mesh_in.point_arrays.items():
+        assert numpy.allclose(v, mesh.point_arrays[k.replace(" ", "_")])
+    for k, v in mesh_in.cell_arrays.items():
+        assert numpy.allclose(v, mesh.cell_arrays[k.replace(" ", "_")])
+
+
+def test_file_format():
+    with pytest.raises(AssertionError):
+        _ = pyvista.read_meshio(examples.hexbeamfile, file_format = "bar")
+
+    with pytest.raises(KeyError):
+        pyvista.save_meshio("foo.bar", beam, file_format = "bar")