Merge branch '112-vtk-fails-to-open-xdmf-for-large-topologies' into 'master'

Resolve "VTK fails to open XDMF for large topologies"

Closes #112

See merge request mechanics-and-materials/ae108!101

tools/generate_xdmf.py

@@ -30,6 +30,7 @@ import re
 import typing
 import sys
 import xml.etree.ElementTree as ET
+import itertools
 
 import h5py
 import numpy
@@ -311,32 +312,31 @@ class MeshInformationMissing(Exception):
     """
 
 
-def read_element_vertices(hdf_file: h5py.File) -> typing.Iterator[typing.List[int]]:
+def read_element_vertices(hdf_file: h5py.File) -> typing.Iterator[numpy.ndarray]:
     """
     Returns the vertex indices of the elements.
     """
     if not "/topology" in hdf_file:
         raise MeshInformationMissing()
 
-    cells_data = hdf_file["/topology/cells"]
     cones_data = hdf_file["/topology/cones"][()]
     number_of_elements = numpy.sum(cones_data > 0)
+    cells_data = hdf_file["/topology/cells"][()] - number_of_elements
 
     offset = 0
-    for cone_size in numpy.nditer(cones_data):
-        if cone_size == 0:
-            continue
-        yield list(
-            numpy.nditer(cells_data[offset : offset + cone_size] - number_of_elements)
-        )
+    for cone_size in filter(lambda x: x > 0, numpy.nditer(cones_data)):
+        yield cells_data[offset : offset + cone_size]
         offset += cone_size
 
 
 def add_topology(
-    parent: ET.Element, hdf_file: h5py.File, prefer_planar: bool
+    parent: ET.Element,
+    hdf_file: h5py.File,
+    prefer_planar: bool,
 ) -> ET.Element:
     """
-    Adds the topology element to the parent.
+    Writes topology data compatible with XDMF to "/topology/mixed"
+    in `hdf_file`, and adds the topology element to the parent.
     """
     if not "/topology" in hdf_file:
         raise MeshInformationMissing()
@@ -346,41 +346,36 @@ def add_topology(
     number_of_elements = numpy.sum(cones_data > 0)
     coordinate_dimension = read_coordinate_dimension(hdf_file)
 
-    dataitem_dimension = sum(
-        cone_size
-        + len(
-            number_of_corners_to_type(
-                int(cone_size), prefer_planar or coordinate_dimension <= 2
-            )
-        )
-        for cone_size in numpy.nditer(cones_data)
-        if cone_size > 0
-    )
-
     topology = ET.SubElement(
         parent,
         "Topology",
         {"TopologyType": "Mixed", "NumberOfElements": str(number_of_elements)},
     )
 
-    dataitem = ET.SubElement(
-        topology,
-        "DataItem",
-        {"Format": "XML", "Dimensions": str(dataitem_dimension)},
+    topology_data = numpy.fromiter(
+        itertools.chain.from_iterable(
+            (
+                itertools.chain(
+                    iter(
+                        number_of_corners_to_type(
+                            vertices.shape[0],
+                            prefer_planar or coordinate_dimension <= 2,
+                        )
+                    ),
+                    numpy.nditer(vertices),
+                )
+                for vertices in read_element_vertices(hdf_file)
+            )
+        ),
+        dtype="int32",
     )
 
-    dataitem.text = "  ".join(
-        f"""{" ".join(str(type_)
-          for type_ in number_of_corners_to_type(
-            len(vertices),
-            prefer_planar or coordinate_dimension <= 2,
-          )
-        )
-        } {
-          " ".join(str(vertex) for vertex in vertices)
-        }"""
-        for vertices in read_element_vertices(hdf_file)
+    hdf_file.require_dataset(
+        "/topology/mixed", shape=topology_data.shape, dtype=topology_data.dtype
     )
+    hdf_file["/topology/mixed"].write_direct(topology_data)
+
+    add_hdf_dataitem(topology, hdf_file, "/topology/mixed")
 
     return topology
 
@@ -532,7 +527,12 @@ def main() -> None:
     """
     filename, prefer_planar = parse_command_line_arguments()
     with open(filename.stem + ".xdmf", "w", encoding="utf-8") as xdmf_file:
-        xdmf_file.write(hdf_to_xdmf_string(h5py.File(filename, "r"), prefer_planar))
+        xdmf_file.write(
+            hdf_to_xdmf_string(
+                h5py.File(filename, "r+"),
+                prefer_planar,
+            )
+        )
 
 
 class ErrorCode(IntEnum):