Merge pull request #997 from nschloe/hmf

hmf

Author: nschloe

meshio/__about__.py

@@ -2,7 +2,10 @@
     # Python 3.8
     from importlib import metadata
 except ImportError:
-    import importlib_metadata as metadata
+    try:
+        import importlib_metadata as metadata
+    except ImportError:
+        __version__ = "unknown"
 
 try:
     __version__ = metadata.version("meshio")

meshio/__init__.py

@@ -9,6 +9,7 @@
     flac3d,
     gmsh,
     h5m,
+    hmf,
     mdpa,
     med,
     medit,
@@ -44,6 +45,7 @@
     "flac3d",
     "gmsh",
     "h5m",
+    "hmf",
     "mdpa",
     "med",
     "medit",

meshio/hmf/__init__.py

@@ -0,0 +1,3 @@
+from ._hmf import read, write
+
+__all__ = ["read", "write"]

meshio/hmf/_hmf.py

@@ -0,0 +1,158 @@
+import warnings
+
+import meshio
+
+from .._common import cell_data_from_raw, raw_from_cell_data
+from .._helpers import register
+from ..xdmf.common import meshio_to_xdmf_type, xdmf_to_meshio_type
+
+
+def read(filename):
+    import h5py
+
+    with h5py.File(filename, "r") as f:
+        assert f.attrs["type"] == "hmf"
+        assert f.attrs["version"] == "0.1-alpha"
+
+        assert len(f) == 1, "only one domain supported for now"
+        domain = f["domain"]
+
+        assert len(domain) == 1, "only one grid supported for now"
+        grid = domain["grid"]
+
+        points = None
+        cells = {}
+        point_data = {}
+        cell_data_raw = {}
+
+        for key, value in grid.items():
+            if key[:8] == "Topology":
+                cell_type = value.attrs["TopologyType"]
+                cells[xdmf_to_meshio_type[cell_type]] = value[()]
+
+            elif key == "Geometry":
+                # TODO is GeometryType really needed?
+                assert value.attrs["GeometryType"] in ["X", "XY", "XYZ"]
+                points = value[()]
+
+            elif key == "CellAttributes":
+                for name, ca in value.items():
+                    cell_data_raw[name] = ca[()]
+
+            else:
+                assert key == "NodeAttributes"
+                for name, na in value.items():
+                    point_data[name] = na[()]
+
+        cell_data = cell_data_from_raw(cells, cell_data_raw)
+
+        return meshio.Mesh(
+            points,
+            cells,
+            point_data=point_data,
+            cell_data=cell_data,
+        )
+
+
+def write_points_cells(filename, points, cells, **kwargs):
+    write(filename, meshio.Mesh(points, cells), **kwargs)
+
+
+def write(filename, mesh, compression="gzip", compression_opts=4):
+    import h5py
+
+    warnings.warn("Experimental file format. Format can change at any time.")
+    with h5py.File(filename, "w") as h5_file:
+        h5_file.attrs["type"] = "hmf"
+        h5_file.attrs["version"] = "0.1-alpha"
+        domain = h5_file.create_group("domain")
+        grid = domain.create_group("grid")
+
+        _write_points(grid, mesh.points, compression, compression_opts)
+        _write_cells(mesh.cells, grid, compression, compression_opts)
+        _write_point_data(mesh.point_data, grid, compression, compression_opts)
+        _write_cell_data(mesh.cell_data, grid, compression, compression_opts)
+
+
+def _write_points(grid, points, compression, compression_opts):
+    geo = grid.create_dataset(
+        "Geometry",
+        data=points,
+        compression=compression,
+        compression_opts=compression_opts,
+    )
+    geo.attrs["GeometryType"] = "XYZ"[: points.shape[1]]
+
+
+def _write_cells(cell_blocks, grid, compression, compression_opts):
+    for k, cell_block in enumerate(cell_blocks):
+        xdmf_type = meshio_to_xdmf_type[cell_block.type][0]
+        topo = grid.create_dataset(
+            f"Topology{k}",
+            data=cell_block.data,
+            compression=compression,
+            compression_opts=compression_opts,
+        )
+        topo.attrs["TopologyType"] = xdmf_type
+
+
+# In XDMF, the point/cell data are stored as
+#
+# <Attribute Name="phi" AttributeType="Scalar" Center="Node">
+#   <DataItem DataType="Float" Dimensions="4" Format="HDF" Precision="8">
+#     out.h5:/data2
+#   </DataItem>
+# </Attribute>
+#
+# We cannot register multiple entries with the same name in HDF, so instead of
+# "Attribute", use
+# ```
+# NodeAttributes
+#   -> name0 + data0
+#   -> name1 + data0
+#   -> ...
+# CellAttributes
+#   -> ...
+# ```
+# Alternative:
+# ```
+#  NodeAttribute0
+#    -> name
+#    -> data
+#  NodeAttribute1
+#    -> name
+#    -> data
+#  ...
+# ```
+# It's done similarly for Topologies (cells).
+#
+def _write_point_data(point_data, grid, compression, compression_opts):
+    na = grid.create_group("NodeAttributes")
+    for name, data in point_data.items():
+        na.create_dataset(
+            name,
+            data=data,
+            compression=compression,
+            compression_opts=compression_opts,
+        )
+
+
+def _write_cell_data(cell_data, grid, compression, compression_opts):
+    raw = raw_from_cell_data(cell_data)
+    ca = grid.create_group("CellAttributes")
+    for name, data in raw.items():
+        ca.create_dataset(
+            name,
+            data=data,
+            compression=compression,
+            compression_opts=compression_opts,
+        )
+
+
+# TODO register all xdmf except hdf outside this try block
+register(
+    "hmf",
+    [".hmf"],
+    read,
+    {"hmf": write},
+)

meshio/vtu/_vtu.py

@@ -5,7 +5,6 @@
 """
 import base64
 import logging
-import lzma
 import re
 import sys
 import zlib
@@ -25,6 +24,12 @@
 from .._helpers import register
 from .._mesh import CellBlock, Mesh
 
+# Paraview 5.8.1's built-in Python doesn't have lzma.
+try:
+    import lzma
+except ModuleNotFoundError:
+    lzma = None
+
 
 def num_bytes_to_num_base64_chars(num_bytes):
     # Rounding up in integer division works by double negation since Python

meshio/xdmf/main.py

@@ -351,11 +351,11 @@ def __init__(
         #     grid, "Information", Name="Information", Value=str(len(mesh.field_data))
         # )
 
-        self.points(grid, mesh.points)
+        self.write_points(grid, mesh.points)
         # self.field_data(mesh.field_data, information)
-        self.cells(mesh.cells, grid)
-        self.point_data(mesh.point_data, grid)
-        self.cell_data(mesh.cell_data, grid)
+        self.write_cells(mesh.cells, grid)
+        self.write_point_data(mesh.point_data, grid)
+        self.write_cell_data(mesh.cell_data, grid)
 
         ET.register_namespace("xi", "https://www.w3.org/2001/XInclude/")
 
@@ -388,15 +388,11 @@ def numpy_to_xml_string(self, data):
         )
         return os.path.basename(self.h5_filename) + ":/" + name
 
-    def points(self, grid, points):
-        if points.shape[1] == 1:
-            geometry_type = "X"
-        elif points.shape[1] == 2:
-            geometry_type = "XY"
-        else:
-            if points.shape[1] != 3:
-                raise WriteError()
-            geometry_type = "XYZ"
+    def write_points(self, grid, points):
+        if points.shape[1] > 3:
+            raise WriteError("Can only write points up to dimension 3.")
+
+        geometry_type = "XYZ"[: points.shape[1]]
 
         geo = ET.SubElement(grid, "Geometry", GeometryType=geometry_type)
         dt, prec = numpy_to_xdmf_dtype[points.dtype.name]
@@ -411,7 +407,7 @@ def points(self, grid, points):
         )
         data_item.text = self.numpy_to_xml_string(points)
 
-    def cells(self, cells, grid):
+    def write_cells(self, cells, grid):
         if len(cells) == 1:
             meshio_type = cells[0].type
             num_cells = len(cells[0].data)
@@ -474,7 +470,7 @@ def cells(self, cells, grid):
             )
             data_item.text = self.numpy_to_xml_string(cd)
 
-    def point_data(self, point_data, grid):
+    def write_point_data(self, point_data, grid):
         for name, data in point_data.items():
             att = ET.SubElement(
                 grid,
@@ -495,7 +491,7 @@ def point_data(self, point_data, grid):
             )
             data_item.text = self.numpy_to_xml_string(data)
 
-    def cell_data(self, cell_data, grid):
+    def write_cell_data(self, cell_data, grid):
         raw = raw_from_cell_data(cell_data)
         for name, data in raw.items():
             att = ET.SubElement(

setup.cfg

@@ -1,6 +1,6 @@
 [metadata]
 name = meshio
-version = 4.3.7
+version = 4.3.8
 author = Nico Schlömer et al.
 author_email = [email protected]
 description = I/O for many mesh formats

test/test_hmf.py

@@ -0,0 +1,37 @@
+import helpers
+import numpy
+import pytest
+
+import meshio
+
+test_set_full = [
+    helpers.line_mesh,
+    helpers.tri_mesh,
+    helpers.line_tri_mesh,
+    helpers.tri_mesh_2d,
+    helpers.triangle6_mesh,
+    helpers.quad_mesh,
+    helpers.quad8_mesh,
+    helpers.tri_quad_mesh,
+    helpers.tet_mesh,
+    helpers.tet10_mesh,
+    helpers.hex_mesh,
+    helpers.hex20_mesh,
+    helpers.add_point_data(helpers.tri_mesh, 1),
+    helpers.add_cell_data(helpers.tri_mesh, [("a", (), numpy.float64)]),
+]
+
+
+@pytest.mark.parametrize("mesh", test_set_full)
+@pytest.mark.parametrize("compression", [None, "gzip"])
+def test_xdmf3(mesh, compression):
+    def write(*args, **kwargs):
+        return meshio.xdmf.write(*args, compression=compression, **kwargs)
+
+    helpers.write_read(write, meshio.xdmf.read, mesh, 1.0e-14)
+
+
+def test_generic_io():
+    helpers.generic_io("test.hmf")
+    # With additional, insignificant suffix:
+    helpers.generic_io("test.0.hmf")

test/test_xdmf.py

@@ -34,19 +34,18 @@
 
 
 @pytest.mark.parametrize("mesh", test_set_full)
-@pytest.mark.parametrize("data_format", ["XML", "Binary", "HDF"])
-def test_xdmf3(mesh, data_format):
+@pytest.mark.parametrize(
+    "kwargs0",
+    [
+        {"data_format": "XML"},
+        {"data_format": "Binary"},
+        {"data_format": "HDF", "compression": None},
+        {"data_format": "HDF", "compression": "gzip"},
+    ],
+)
+def test_xdmf3(mesh, kwargs0):
     def write(*args, **kwargs):
-        return meshio.xdmf.write(*args, data_format=data_format, **kwargs)
-
-    helpers.write_read(write, meshio.xdmf.read, mesh, 1.0e-14)
-
-
-# HDF5 compressed I/O
-@pytest.mark.parametrize("mesh", test_set_full)
-def test_compression(mesh):
-    def write(*args, **kwargs):
-        return meshio.xdmf.write(*args, data_format="HDF", compression="gzip", **kwargs)
+        return meshio.xdmf.write(*args, **{**kwargs0, **kwargs})
 
     helpers.write_read(write, meshio.xdmf.read, mesh, 1.0e-14)