Fix Exodus HEX20 node ordering on read and write

src/meshio/exodus/_exodus.py

@@ -65,6 +65,10 @@
 }
 meshio_to_exodus_type = {v: k for k, v in exodus_to_meshio_type.items()}
 
+# HEX20 orders the vertical midside nodes before the top-face midside nodes,
+# unlike meshio's VTK ordering. This permutation is its own inverse.
+_hex20_order = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 12, 13, 14, 15]
+
 
 def read(filename):  # noqa: C901
     import netCDF4
@@ -105,7 +109,10 @@ def read(filename):  # noqa: C901
                     info += [b"".join(c).decode("UTF-8") for c in val[:]]
             elif key[:7] == "connect":
                 meshio_type = exodus_to_meshio_type[value.elem_type.upper()]
-                cells.append((meshio_type, value[:] - 1))
+                connectivity = value[:] - 1
+                if meshio_type == "hexahedron20":
+                    connectivity = connectivity[:, _hex20_order]
+                cells.append((meshio_type, connectivity))
             elif key == "coord":
                 points = nc.variables["coord"][:].T
             elif key == "coordx":
@@ -316,8 +323,11 @@ def write(filename, mesh):
             dtype = numpy_to_exodus_dtype[cell_block.data.dtype.name]
             data = rootgrp.createVariable(f"connect{k + 1}", dtype, (dim1, dim2))
             data.elem_type = meshio_to_exodus_type[cell_block.type]
+            connectivity = cell_block.data
+            if cell_block.type == "hexahedron20":
+                connectivity = connectivity[:, _hex20_order]
             # Exodus is 1-based
-            data[:] = cell_block.data + 1
+            data[:] = connectivity + 1
 
         # point data
         # The variable `name_nod_var` holds the names and indices of the node variables, the

tests/test_exodus.py

@@ -1,3 +1,4 @@
+import numpy as np
 import pytest
 
 import meshio
@@ -35,3 +36,84 @@ def test_generic_io(tmp_path):
     helpers.generic_io(tmp_path / "test.e")
     # With additional, insignificant suffix:
     helpers.generic_io(tmp_path / "test.0.e")
+
+
+@pytest.mark.parametrize("mixed", [False, True])
+@pytest.mark.parametrize("elem_type", ["HEX20", "hex20"])
+@pytest.mark.parametrize("dtype", [np.int32, np.int64])
+def test_read_hex20_connectivity(tmp_path, mixed, elem_type, dtype):
+    # Exodus numbers the vertical midside nodes before the top-face midside
+    # nodes (Exodus manual, Figure 4.13); meshio uses VTK's opposite order.
+    connectivity = np.arange(1, 41, dtype=dtype).reshape(2, 20)
+    expected = np.array(
+        [[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 12, 13, 14, 15]]
+    ) + [[0], [20]]
+    blocks = [(elem_type, "hexahedron20", connectivity, expected)]
+    if mixed:
+        blocks = [
+            ("HEX8", "hexahedron", connectivity[:, :8], connectivity[:, :8] - 1),
+            *blocks,
+            ("QUAD8", "quad8", connectivity[:, :8], connectivity[:, :8] - 1),
+            *blocks,
+        ]
+    points = np.arange(120, dtype=float).reshape(40, 3)
+    filename = tmp_path / "hex20.e"
+    # Construct the input directly, independently of meshio's writer.
+    with netCDF4.Dataset(filename, "w") as nc:
+        nc.createDimension("num_nodes", len(points))
+        nc.createDimension("num_dim", 3)
+        nc.createVariable("coord", "f8", ("num_dim", "num_nodes"))[:] = points.T
+        for k, (exodus_type, _, data, _) in enumerate(blocks, start=1):
+            dims = (f"num_el_in_blk{k}", f"num_nod_per_el{k}")
+            nc.createDimension(dims[0], data.shape[0])
+            nc.createDimension(dims[1], data.shape[1])
+            variable = nc.createVariable(f"connect{k}", data.dtype, dims)
+            variable.elem_type = exodus_type
+            variable[:] = data
+
+    mesh = meshio.exodus.read(filename)
+    np.testing.assert_array_equal(mesh.points, points)
+    assert len(mesh.cells) == len(blocks)
+    for cell, (_, meshio_type, data, expected) in zip(mesh.cells, blocks):
+        assert cell.type == meshio_type
+        # Preserve the existing masked-array subtraction's dtype, which may
+        # differ from the on-disk dtype between NumPy versions.
+        assert cell.data.dtype == (np.ma.array(data) - 1).dtype
+        np.testing.assert_array_equal(cell.data, expected)
+
+
+@pytest.mark.parametrize("mixed", [False, True])
+@pytest.mark.parametrize("dtype", [np.int32, np.int64])
+def test_write_hex20_connectivity(tmp_path, mixed, dtype):
+    connectivity = np.arange(40, dtype=dtype).reshape(2, 20)
+    expected = np.array(
+        [[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 17, 18, 19, 20, 13, 14, 15, 16]]
+    ) + [[0], [20]]
+    blocks = [("hexahedron20", "HEX20", connectivity, expected)]
+    if mixed:
+        blocks = [
+            ("hexahedron", "HEX8", connectivity[:, :8], connectivity[:, :8] + 1),
+            *blocks,
+            ("quad8", "QUAD8", connectivity[:, :8], connectivity[:, :8] + 1),
+            *blocks,
+        ]
+    mesh = meshio.Mesh(
+        np.arange(120, dtype=float).reshape(40, 3),
+        [(cell_type, data) for cell_type, _, data, _ in blocks],
+    )
+    original = [cell.data.copy() for cell in mesh.cells]
+    for cell in mesh.cells:
+        cell.data.setflags(write=False)
+    filename = tmp_path / "hex20.e"
+    meshio.exodus.write(filename, mesh)
+
+    # Inspect raw Exodus connectivity, independently of meshio's reader.
+    with netCDF4.Dataset(filename) as nc:
+        assert len(nc.dimensions["num_el_blk"]) == len(blocks)
+        for k, (_, exodus_type, _, expected) in enumerate(blocks, start=1):
+            variable = nc.variables[f"connect{k}"]
+            assert variable.elem_type == exodus_type
+            assert variable.dtype == np.dtype(dtype)
+            np.testing.assert_array_equal(variable[:], expected)
+    for cell, data in zip(mesh.cells, original):
+        np.testing.assert_array_equal(cell.data, data)