diff --git a/.gitignore b/.gitignore index c31596766..2dfc807bb 100644 --- a/.gitignore +++ b/.gitignore @@ -22,3 +22,4 @@ doc/_build/ .pytest_cache/ .coverage .tox/ +foo.vtk diff --git a/setup.cfg b/setup.cfg index 032cb0125..c2e2af1b1 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,6 +1,6 @@ [metadata] name = meshio -version = 4.4.5 +version = 4.4.6 author = Nico Schlömer et al. author_email = nico.schloemer@gmail.com description = I/O for many mesh formats diff --git a/src/meshio/_common.py b/src/meshio/_common.py index de59b5155..c1e752e19 100644 --- a/src/meshio/_common.py +++ b/src/meshio/_common.py @@ -183,87 +183,3 @@ def _pick_first_int_data(data): other = [] return key, other - - -# https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html -vtk_to_meshio_type = { - 0: "empty", - 1: "vertex", - # 2: 'poly_vertex', - 3: "line", - # 4: 'poly_line', - 5: "triangle", - # 6: 'triangle_strip', - 7: "polygon", - 8: "pixel", - 9: "quad", - 10: "tetra", - # 11: 'voxel', - 12: "hexahedron", - 13: "wedge", - 14: "pyramid", - 15: "penta_prism", - 16: "hexa_prism", - 21: "line3", - 22: "triangle6", - 23: "quad8", - 24: "tetra10", - 25: "hexahedron20", - 26: "wedge15", - 27: "pyramid13", - 28: "quad9", - 29: "hexahedron27", - 30: "quad6", - 31: "wedge12", - 32: "wedge18", - 33: "hexahedron24", - 34: "triangle7", - 35: "line4", - 42: "polyhedron", - # - # 60: VTK_HIGHER_ORDER_EDGE, - # 61: VTK_HIGHER_ORDER_TRIANGLE, - # 62: VTK_HIGHER_ORDER_QUAD, - # 63: VTK_HIGHER_ORDER_POLYGON, - # 64: VTK_HIGHER_ORDER_TETRAHEDRON, - # 65: VTK_HIGHER_ORDER_WEDGE, - # 66: VTK_HIGHER_ORDER_PYRAMID, - # 67: VTK_HIGHER_ORDER_HEXAHEDRON, - # Arbitrary order Lagrange elements - 68: "VTK_LAGRANGE_CURVE", - 69: "VTK_LAGRANGE_TRIANGLE", - 70: "VTK_LAGRANGE_QUADRILATERAL", - 71: "VTK_LAGRANGE_TETRAHEDRON", - 72: "VTK_LAGRANGE_HEXAHEDRON", - 73: "VTK_LAGRANGE_WEDGE", - 74: "VTK_LAGRANGE_PYRAMID", - # Arbitrary order Bezier elements - 75: "VTK_BEZIER_CURVE", - 76: "VTK_BEZIER_TRIANGLE", - 77: "VTK_BEZIER_QUADRILATERAL", - 78: "VTK_BEZIER_TETRAHEDRON", - 79: "VTK_BEZIER_HEXAHEDRON", - 80: "VTK_BEZIER_WEDGE", - 81: "VTK_BEZIER_PYRAMID", -} -meshio_to_vtk_type = {v: k for k, v in vtk_to_meshio_type.items()} - - -def _vtk_to_meshio_order(vtk_type, numnodes, dtype=int): - # meshio uses the same node ordering as VTK for most cell types. However, for the - # linear wedge, the ordering of the gmsh Prism [1] is adopted since this is found in - # most codes (Abaqus, Ansys, Nastran,...). In the vtkWedge [2], the normal of the - # (0,1,2) triangle points outwards, while in gmsh this normal points inwards. - # [1] http://gmsh.info/doc/texinfo/gmsh.html#Node-ordering - # [2] https://vtk.org/doc/nightly/html/classvtkWedge.html - if vtk_type == 13: - return np.array([0, 2, 1, 3, 5, 4], dtype=dtype) - else: - return np.arange(0, numnodes, dtype=dtype) - - -def _meshio_to_vtk_order(meshio_type, numnodes, dtype=int): - if meshio_type == "wedge": - return np.array([0, 2, 1, 3, 5, 4], dtype=dtype) - else: - return np.arange(0, numnodes, dtype=dtype) diff --git a/src/meshio/_vtk_common.py b/src/meshio/_vtk_common.py new file mode 100644 index 000000000..c28dc3e6c --- /dev/null +++ b/src/meshio/_vtk_common.py @@ -0,0 +1,179 @@ +import warnings + +import numpy as np + +from ._common import num_nodes_per_cell +from ._exceptions import ReadError +from ._mesh import CellBlock + +# https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html +vtk_to_meshio_type = { + 0: "empty", + 1: "vertex", + # 2: 'poly_vertex', + 3: "line", + # 4: 'poly_line', + 5: "triangle", + # 6: 'triangle_strip', + 7: "polygon", + 8: "pixel", + 9: "quad", + 10: "tetra", + # 11: 'voxel', + 12: "hexahedron", + 13: "wedge", + 14: "pyramid", + 15: "penta_prism", + 16: "hexa_prism", + 21: "line3", + 22: "triangle6", + 23: "quad8", + 24: "tetra10", + 25: "hexahedron20", + 26: "wedge15", + 27: "pyramid13", + 28: "quad9", + 29: "hexahedron27", + 30: "quad6", + 31: "wedge12", + 32: "wedge18", + 33: "hexahedron24", + 34: "triangle7", + 35: "line4", + 42: "polyhedron", + # + # 60: VTK_HIGHER_ORDER_EDGE, + # 61: VTK_HIGHER_ORDER_TRIANGLE, + # 62: VTK_HIGHER_ORDER_QUAD, + # 63: VTK_HIGHER_ORDER_POLYGON, + # 64: VTK_HIGHER_ORDER_TETRAHEDRON, + # 65: VTK_HIGHER_ORDER_WEDGE, + # 66: VTK_HIGHER_ORDER_PYRAMID, + # 67: VTK_HIGHER_ORDER_HEXAHEDRON, + # Arbitrary order Lagrange elements + 68: "VTK_LAGRANGE_CURVE", + 69: "VTK_LAGRANGE_TRIANGLE", + 70: "VTK_LAGRANGE_QUADRILATERAL", + 71: "VTK_LAGRANGE_TETRAHEDRON", + 72: "VTK_LAGRANGE_HEXAHEDRON", + 73: "VTK_LAGRANGE_WEDGE", + 74: "VTK_LAGRANGE_PYRAMID", + # Arbitrary order Bezier elements + 75: "VTK_BEZIER_CURVE", + 76: "VTK_BEZIER_TRIANGLE", + 77: "VTK_BEZIER_QUADRILATERAL", + 78: "VTK_BEZIER_TETRAHEDRON", + 79: "VTK_BEZIER_HEXAHEDRON", + 80: "VTK_BEZIER_WEDGE", + 81: "VTK_BEZIER_PYRAMID", +} +meshio_to_vtk_type = {v: k for k, v in vtk_to_meshio_type.items()} + + +def vtk_to_meshio_order(vtk_type, numnodes, dtype=int): + # meshio uses the same node ordering as VTK for most cell types. However, for the + # linear wedge, the ordering of the gmsh Prism [1] is adopted since this is found in + # most codes (Abaqus, Ansys, Nastran,...). In the vtkWedge [2], the normal of the + # (0,1,2) triangle points outwards, while in gmsh this normal points inwards. + # [1] http://gmsh.info/doc/texinfo/gmsh.html#Node-ordering + # [2] https://vtk.org/doc/nightly/html/classvtkWedge.html + if vtk_type == 13: + return np.array([0, 2, 1, 3, 5, 4], dtype=dtype) + else: + return np.arange(0, numnodes, dtype=dtype) + + +def meshio_to_vtk_order(meshio_type, numnodes, dtype=int): + if meshio_type == "wedge": + return np.array([0, 2, 1, 3, 5, 4], dtype=dtype) + else: + return np.arange(0, numnodes, dtype=dtype) + + +def vtk_cells_from_data(connectivity, offsets, types, cell_data_raw): + # Translate it into the cells array. + # `connectivity` is a one-dimensional vector with + # (p00, p01, ... ,p0k, p10, p11, ..., p1k, ... + # `offsets` is a pointer array that points to the first position of p0, p1, etc. + if len(offsets) != len(types): + raise ReadError(f"len(offsets) != len(types) ({len(offsets)} != {len(types)})") + + # identify cell blocks + breaks = np.where(types[:-1] != types[1:])[0] + 1 + # all cells with indices between start[k] and end[k] have the same type + start_end = list( + zip( + np.concatenate([[0], breaks]), + np.concatenate([breaks, [len(types)]]), + ) + ) + + cells = [] + cell_data = {} + + for start, end in start_end: + try: + meshio_type = vtk_to_meshio_type[types[start]] + except KeyError: + warnings.warn( + f"File contains cells that meshio cannot handle (type {types[start]})." + ) + continue + + # cells with varying number of points + special_cells = [ + "polygon", + "VTK_LAGRANGE_CURVE", + "VTK_LAGRANGE_TRIANGLE", + "VTK_LAGRANGE_QUADRILATERAL", + "VTK_LAGRANGE_TETRAHEDRON", + "VTK_LAGRANGE_HEXAHEDRON", + "VTK_LAGRANGE_WEDGE", + "VTK_LAGRANGE_PYRAMID", + ] + if meshio_type in special_cells: + # Polygons have unknown and varying number of nodes per cell. + + # Index where the previous block of cells stopped. Needed to know the number + # of nodes for the first cell in the block. + first_node = 0 if start == 0 else offsets[start - 1] + + # Start off the cell-node relation for each cell in this block + start_cn = np.hstack((first_node, offsets[start:end])) + # Find the size of each cell except the last + sizes = np.diff(start_cn) + + # find where the cell blocks start and end + b = np.diff(sizes) + c = np.concatenate([[0], np.where(b != 0)[0] + 1, [len(sizes)]]) + + # Loop over all cell sizes, find all cells with this size, and assign + # connectivity + for cell_block_start, cell_block_end in zip(c, c[1:]): + items = np.arange(cell_block_start, cell_block_end) + sz = sizes[cell_block_start] + indices = np.add.outer( + start_cn[items + 1], + vtk_to_meshio_order(types[start], sz, dtype=offsets.dtype) - sz, + ) + cells.append(CellBlock(meshio_type, connectivity[indices])) + + # Store cell data for this set of cells + for name, d in cell_data_raw.items(): + if name not in cell_data: + cell_data[name] = [] + cell_data[name].append(d[start + items]) + else: + # Non-polygonal cell. Same number of nodes per cell makes everything easier. + n = num_nodes_per_cell[meshio_type] + indices = np.add.outer( + offsets[start:end], + vtk_to_meshio_order(types[start], n, dtype=offsets.dtype) - n, + ) + cells.append(CellBlock(meshio_type, connectivity[indices])) + for name, d in cell_data_raw.items(): + if name not in cell_data: + cell_data[name] = [] + cell_data[name].append(d[start:end]) + + return cells, cell_data diff --git a/src/meshio/vtk/_vtk.py b/src/meshio/vtk/_vtk.py index 5442424f4..bd01a5124 100644 --- a/src/meshio/vtk/_vtk.py +++ b/src/meshio/vtk/_vtk.py @@ -7,16 +7,17 @@ import numpy as np from ..__about__ import __version__ -from .._common import ( - _meshio_to_vtk_order, - _vtk_to_meshio_order, - meshio_to_vtk_type, - vtk_to_meshio_type, -) from .._exceptions import ReadError, WriteError from .._files import open_file from .._helpers import register from .._mesh import CellBlock, Mesh +from .._vtk_common import ( + meshio_to_vtk_order, + meshio_to_vtk_type, + vtk_cells_from_data, + vtk_to_meshio_order, + vtk_to_meshio_type, +) vtk_type_to_numnodes = np.array( [ @@ -130,8 +131,9 @@ def __init__(self): self.cell_data_raw = {} self.point_data = {} self.dataset = {} - self.c = None - self.ct = None + self.connectivity = None + self.offsets = None + self.types = None self.active = None self.is_ascii = False self.split = [] @@ -183,7 +185,15 @@ def read_buffer(f): _read_subsection(f, info) _check_mesh(info) - cells, cell_data = translate_cells(info.c, info.ct, info.cell_data_raw) + + if info.offsets is not None: + cells, cell_data = vtk_cells_from_data( + info.connectivity, info.offsets, info.types, info.cell_data_raw + ) + else: + cells, cell_data = translate_cells( + info.connectivity, info.types, info.cell_data_raw + ) return Mesh( info.points, @@ -225,7 +235,7 @@ def _read_section(f, info): # vtk DataFile Version 5.1 - appearing in Paraview 5.8.1 outputs # No specification found for this file format. # See the question on ParaView Discourse Forum: - # . + # https://discourse.paraview.org/t/specification-of-vtk-datafile-version-5-1/5127 info.num_offsets = int(info.split[1]) info.num_items = int(info.split[2]) dtype = np.dtype(vtk_to_numpy_dtype_name[line.split()[1]]) @@ -234,16 +244,19 @@ def _read_section(f, info): assert "CONNECTIVITY" in line dtype = np.dtype(vtk_to_numpy_dtype_name[line.split()[1]]) connectivity = _read_cells(f, info.is_ascii, info.num_items, dtype) - info.c = (offsets, connectivity) + info.connectivity = connectivity + assert offsets[0] == 0 + assert offsets[-1] == len(connectivity) + info.offsets = offsets[1:] else: f.seek(last_pos) info.num_items = int(info.split[2]) - info.c = _read_cells(f, info.is_ascii, info.num_items) + info.connectivity = _read_cells(f, info.is_ascii, info.num_items) elif info.section == "CELL_TYPES": info.active = "CELL_TYPES" info.num_items = int(info.split[1]) - info.ct = _read_cell_types(f, info.is_ascii, info.num_items) + info.types = _read_cell_types(f, info.is_ascii, info.num_items) elif info.section == "POINT_DATA": info.active = "POINT_DATA" @@ -308,9 +321,9 @@ def _read_subsection(f, info): def _check_mesh(info): if info.dataset["type"] == "UNSTRUCTURED_GRID": - if info.c is None: + if info.connectivity is None: raise ReadError("Required section CELLS not found.") - if info.ct is None: + if info.types is None: raise ReadError("Required section CELL_TYPES not found.") elif info.dataset["type"] == "STRUCTURED_POINTS": dim = info.dataset["DIMENSIONS"] @@ -325,7 +338,7 @@ def _check_mesh(info): for i in range(3) ] info.points = _generate_points(axis) - info.c, info.ct = _generate_cells(dim=info.dataset["DIMENSIONS"]) + info.connectivity, info.types = _generate_cells(dim=info.dataset["DIMENSIONS"]) elif info.dataset["type"] == "RECTILINEAR_GRID": axis = [ info.dataset["X_COORDINATES"], @@ -333,9 +346,9 @@ def _check_mesh(info): info.dataset["Z_COORDINATES"], ] info.points = _generate_points(axis) - info.c, info.ct = _generate_cells(dim=info.dataset["DIMENSIONS"]) + info.connectivity, info.types = _generate_cells(dim=info.dataset["DIMENSIONS"]) elif info.dataset["type"] == "STRUCTURED_GRID": - info.c, info.ct = _generate_cells(dim=info.dataset["DIMENSIONS"]) + info.connectivity, info.types = _generate_cells(dim=info.dataset["DIMENSIONS"]) def _generate_cells(dim): @@ -552,7 +565,7 @@ def _skip_meta(f): break -def translate_cells(data, types, cell_data_raw): +def translate_cells(connectivity, types, cell_data_raw): # https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html # Translate it into the cells array. # `data` is a one-dimensional vector with @@ -569,21 +582,21 @@ def translate_cells(data, types, cell_data_raw): offsets = np.empty(len(types), dtype=int) offsets[0] = 0 for idx in range(numcells - 1): - numnodes[idx] = data[offsets[idx]] + numnodes[idx] = connectivity[offsets[idx]] offsets[idx + 1] = offsets[idx] + numnodes[idx] + 1 idx = numcells - 1 - numnodes[idx] = data[offsets[idx]] - if not np.all(numnodes == data[offsets]): + numnodes[idx] = connectivity[offsets[idx]] + if not np.all(numnodes == connectivity[offsets]): raise ReadError() # TODO: cell_data for idx, vtk_cell_type in enumerate(types): start = offsets[idx] + 1 - cell_idx = start + _vtk_to_meshio_order( + cell_idx = start + vtk_to_meshio_order( vtk_cell_type, numnodes[idx], offsets.dtype ) - cell = data[cell_idx] + cell = connectivity[cell_idx] cell_type = vtk_to_meshio_type[vtk_cell_type] @@ -603,24 +616,21 @@ def translate_cells(data, types, cell_data_raw): for k, c in enumerate(cells): cells[k] = CellBlock(c.type, np.array(c.data)) else: - # Deduct offsets from the cell types. This is much faster than manually going + # Infer offsets from the cell types. This is much faster than manually going # through the data array. Slight disadvantage: This doesn't work for cells with # a custom number of points. + if np.any(types >= len(vtk_type_to_numnodes)): + raise ReadError("File contains cells that meshio cannot handle.") + numnodes = vtk_type_to_numnodes[types] if not np.all(numnodes > 0): raise ReadError("File contains cells that meshio cannot handle.") - if isinstance(data, tuple): - offsets, conn = data - if not np.all(numnodes == np.diff(offsets)): - raise ReadError() - idx0 = 0 - else: - offsets = np.cumsum(numnodes + 1) - (numnodes + 1) - if not np.all(numnodes == data[offsets]): - raise ReadError() - idx0 = 1 - conn = data + offsets = np.cumsum(numnodes + 1) - (numnodes + 1) + + if not np.all(numnodes == connectivity[offsets]): + raise ReadError() + idx0 = 1 b = np.concatenate( [[0], np.where(types[:-1] != types[1:])[0] + 1, [len(types)]] @@ -630,9 +640,9 @@ def translate_cells(data, types, cell_data_raw): continue meshio_type = vtk_to_meshio_type[types[start]] n = numnodes[start] - cell_idx = idx0 + _vtk_to_meshio_order(types[start], n, dtype=offsets.dtype) + cell_idx = idx0 + vtk_to_meshio_order(types[start], n, dtype=offsets.dtype) indices = np.add.outer(offsets[start:end], cell_idx) - cells.append(CellBlock(meshio_type, conn[indices])) + cells.append(CellBlock(meshio_type, connectivity[indices])) for name, d in cell_data_raw.items(): if name not in cell_data: cell_data[name] = [] @@ -728,7 +738,7 @@ def _write_cells(f, cells, binary): if binary: for c in cells: n = c.data.shape[1] - cell_idx = _meshio_to_vtk_order(c.type, n) + cell_idx = meshio_to_vtk_order(c.type, n) dtype = np.dtype(">i4") # One must force endianness here: # @@ -743,7 +753,7 @@ def _write_cells(f, cells, binary): # ascii for c in cells: n = c.data.shape[1] - cell_idx = _meshio_to_vtk_order(c.type, n) + cell_idx = meshio_to_vtk_order(c.type, n) # prepend a column with the value n np.column_stack( [ diff --git a/src/meshio/vtu/_vtu.py b/src/meshio/vtu/_vtu.py index f0d99ba56..e893458da 100644 --- a/src/meshio/vtu/_vtu.py +++ b/src/meshio/vtu/_vtu.py @@ -7,23 +7,16 @@ import logging import re import sys -import warnings import zlib import numpy as np from ..__about__ import __version__ -from .._common import ( - _meshio_to_vtk_order, - _vtk_to_meshio_order, - meshio_to_vtk_type, - num_nodes_per_cell, - raw_from_cell_data, - vtk_to_meshio_type, -) +from .._common import raw_from_cell_data from .._exceptions import ReadError from .._helpers import register from .._mesh import CellBlock, Mesh +from .._vtk_common import meshio_to_vtk_order, meshio_to_vtk_type, vtk_cells_from_data # Paraview 5.8.1's built-in Python doesn't have lzma. try: @@ -38,86 +31,6 @@ def num_bytes_to_num_base64_chars(num_bytes): return -(-num_bytes // 3) * 4 -def _cells_from_data(connectivity, offsets, types, cell_data_raw): - # Translate it into the cells array. - # `connectivity` is a one-dimensional vector with - # (p0, p1, ... ,pk, p10, p11, ..., p1k, ... - if len(offsets) != len(types): - raise ReadError() - - b = np.concatenate([[0], np.where(types[:-1] != types[1:])[0] + 1, [len(types)]]) - - cells = [] - cell_data = {} - - for start, end in zip(b[:-1], b[1:]): - try: - meshio_type = vtk_to_meshio_type[types[start]] - except KeyError: - warnings.warn( - f"File contains cells that meshio cannot handle (type {types[start]})." - ) - continue - - # cells with varying number of points - special_cells = [ - "polygon", - "VTK_LAGRANGE_CURVE", - "VTK_LAGRANGE_TRIANGLE", - "VTK_LAGRANGE_QUADRILATERAL", - "VTK_LAGRANGE_TETRAHEDRON", - "VTK_LAGRANGE_HEXAHEDRON", - "VTK_LAGRANGE_WEDGE", - "VTK_LAGRANGE_PYRAMID", - ] - if meshio_type in special_cells: - # Polygons have unknown and varying number of nodes per cell. - - # Index where the previous block of cells stopped. Needed to know the number - # of nodes for the first cell in the block. - first_node = 0 if start == 0 else offsets[start - 1] - - # Start off the cell-node relation for each cell in this block - start_cn = np.hstack((first_node, offsets[start:end])) - # Find the size of each cell - sizes = np.diff(start_cn) - - # find where the cell blocks start and end - b = np.diff(sizes) - c = np.concatenate([[0], np.where(b != 0)[0] + 1, [len(sizes)]]) - - # Loop over all cell sizes, find all cells with this size, and assign - # connectivity - for cell_block_start, cell_block_end in zip(c, c[1:]): - items = np.arange(cell_block_start, cell_block_end) - sz = sizes[cell_block_start] - indices = np.add.outer( - start_cn[items + 1], - _vtk_to_meshio_order(types[start], sz, dtype=offsets.dtype) - sz, - ) - cells.append(CellBlock(meshio_type, connectivity[indices])) - - # Store cell data for this set of cells - for name, d in cell_data_raw.items(): - if name not in cell_data: - cell_data[name] = [] - cell_data[name].append(d[start + items]) - else: - # Non-polygonal cell. Same number of nodes per cell makes everything easier. - n = num_nodes_per_cell[meshio_type] - indices = np.add.outer( - offsets[start:end], - _vtk_to_meshio_order(types[start], n, dtype=offsets.dtype) - n, - ) - cells.append(CellBlock(meshio_type, connectivity[indices])) - for name, d in cell_data_raw.items(): - if name not in cell_data: - cell_data[name] = [] - cell_data[name].append(d[start:end]) - - return cells, cell_data - - def _polyhedron_cells_from_data(offsets, faces, faceoffsets, cell_data_raw): # In general the number of faces will vary between cells, and the # number of nodes vary between faces for each cell. The information @@ -232,7 +145,7 @@ def _organize_cells(point_offsets, cells, cell_data_raw): else: for offset, cls, cdr in zip(point_offsets, cells, cell_data_raw): - cls, cell_data = _cells_from_data( + cls, cell_data = vtk_cells_from_data( cls["connectivity"].ravel(), cls["offsets"].ravel(), cls["types"].ravel(), @@ -889,7 +802,7 @@ def _polyhedron_face_cells(face_cells): # create connectivity, offset, type arrays connectivity = np.concatenate( [ - v.data[:, _meshio_to_vtk_order(v.type, v.data.shape[1])].reshape(-1) + v.data[:, meshio_to_vtk_order(v.type, v.data.shape[1])].reshape(-1) for v in mesh.cells ] ) diff --git a/src/meshio/xdmf/main.py b/src/meshio/xdmf/main.py index bd506e3e3..4f22135da 100644 --- a/src/meshio/xdmf/main.py +++ b/src/meshio/xdmf/main.py @@ -160,14 +160,19 @@ def read_xdmf2(self, root): # noqa: C901 for c in grid: if c.tag == "Topology": - print(c) data_items = list(c) - print(data_items) if len(data_items) != 1: - raise ReadError( + message = ( "Need exactly 1 data item in , " f"found {len(data_items)}." ) + if len(data_items) == 0: + message += ( + "\nStructured meshes are not supported, see " + "." + ) + raise ReadError(message) + topology_type = c.get("TopologyType") if topology_type == "Mixed": cells = translate_mixed_cells( @@ -183,8 +188,9 @@ def read_xdmf2(self, root): # noqa: C901 cells.append(CellBlock(xdmf_to_meshio_type[topology_type], data)) elif c.tag == "Geometry": - if c.get("GeometryType") not in (None, "XYZ"): - raise ReadError() + geo_type = c.get("GeometryType") + if geo_type not in (None, "XYZ"): + raise ReadError(f'Expected GeometryType "XYZ", not {geo_type}.') data_items = list(c) if len(data_items) != 1: raise ReadError() diff --git a/tests/__init__.py b/tests/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/helpers.py b/tests/helpers.py index e8c56c5ff..5795ae6d9 100644 --- a/tests/helpers.py +++ b/tests/helpers.py @@ -17,6 +17,11 @@ [("line", [[0, 1], [0, 2], [0, 3], [1, 2], [2, 3]])], ) +tri_mesh_one_cell = meshio.Mesh( + [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 1.0, 0.0]], + [("triangle", [[0, 1, 2]])], +) + tri_mesh_2d = meshio.Mesh( [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]], [("triangle", [[0, 1, 2], [0, 2, 3]])], @@ -238,6 +243,19 @@ ], ) +polygon_mesh_one_cell = meshio.Mesh( + [ + [1.0, 0.0, 0.0], + [1.0, 1.0, 0.0], + [1.5, 0.0, 0.0], + [1.7, 0.5, 0.0], + [1.5, 1.2, 0.0], + ], + [ + ("polygon", [[0, 2, 3, 4, 1]]), + ], +) + # Make sure that the polygon cell blocking works. # This mesh is identical with tri_quad_mesh. polygon2_mesh = meshio.Mesh( @@ -318,6 +336,261 @@ ], ) +# From : +lagrange_high_order_mesh = meshio.Mesh( + [ + [0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0], + [0.0, 0.0, 1.0], + [0.14285714924, 0.0, 0.0], + [0.28571429849, 0.0, 0.0], + [0.42857143283, 0.0, 0.0], + [0.57142859697, 0.0, 0.0], + [0.71428573132, 0.0, 0.0], + [0.85714286566, 0.0, 0.0], + [0.0, 0.14285714924, 0.0], + [0.14285714924, 0.14285714924, 0.0], + [0.28571428359, 0.14285714924, 0.0], + [0.42857144773, 0.14285714924, 0.0], + [0.57142858207, 0.14285714924, 0.0], + [0.71428571641, 0.14285714924, 0.0], + [0.85714285076, 0.14285714924, 0.0], + [0.0, 0.28571429849, 0.0], + [0.14285713434, 0.28571429849, 0.0], + [0.28571429849, 0.28571429849, 0.0], + [0.42857143283, 0.28571429849, 0.0], + [0.57142856717, 0.28571429849, 0.0], + [0.71428570151, 0.28571429849, 0.0], + [0.0, 0.42857143283, 0.0], + [0.14285716414, 0.42857143283, 0.0], + [0.28571429849, 0.42857143283, 0.0], + [0.42857143283, 0.42857143283, 0.0], + [0.57142856717, 0.42857143283, 0.0], + [0.0, 0.57142859697, 0.0], + [0.14285713434, 0.57142859697, 0.0], + [0.28571426868, 0.57142859697, 0.0], + [0.42857140303, 0.57142859697, 0.0], + [0.0, 0.71428573132, 0.0], + [0.14285713434, 0.71428573132, 0.0], + [0.28571426868, 0.71428573132, 0.0], + [0.0, 0.85714286566, 0.0], + [0.14285713434, 0.85714286566, 0.0], + [0.0, 0.0, 0.14285714924], + [0.14285714179, 0.0, 0.14285714924], + [0.28571429104, 0.0, 0.14285714924], + [0.42857142538, 0.0, 0.14285714924], + [0.57142855972, 0.0, 0.14285714924], + [0.71428569406, 0.0, 0.14285714924], + [0.8571428284, 0.0, 0.14285714924], + [0.0, 0.14285714179, 0.14285714924], + [0.14285714924, 0.14285714179, 0.14285714924], + [0.28571428359, 0.14285714179, 0.14285714924], + [0.42857141793, 0.14285714179, 0.14285714924], + [0.57142855227, 0.14285714179, 0.14285714924], + [0.71428568661, 0.14285714179, 0.14285714924], + [0.0, 0.28571429104, 0.14285714924], + [0.14285713434, 0.28571429104, 0.14285714924], + [0.28571426868, 0.28571429104, 0.14285714924], + [0.42857140303, 0.28571429104, 0.14285714924], + [0.57142853737, 0.28571429104, 0.14285714924], + [0.0, 0.42857142538, 0.14285714924], + [0.14285713434, 0.42857142538, 0.14285714924], + [0.28571426868, 0.42857142538, 0.14285714924], + [0.42857140303, 0.42857142538, 0.14285714924], + [0.0, 0.57142855972, 0.14285714924], + [0.14285713434, 0.57142855972, 0.14285714924], + [0.28571426868, 0.57142855972, 0.14285714924], + [0.0, 0.71428569406, 0.14285714924], + [0.14285713434, 0.71428569406, 0.14285714924], + [0.0, 0.8571428284, 0.14285714924], + [0.0, 0.0, 0.28571429849], + [0.14285714924, 0.0, 0.28571429849], + [0.28571428359, 0.0, 0.28571429849], + [0.42857144773, 0.0, 0.28571429849], + [0.57142858207, 0.0, 0.28571429849], + [0.71428571641, 0.0, 0.28571429849], + [0.0, 0.14285714924, 0.28571429849], + [0.14285713434, 0.14285714924, 0.28571429849], + [0.28571429849, 0.14285714924, 0.28571429849], + [0.42857143283, 0.14285714924, 0.28571429849], + [0.57142856717, 0.14285714924, 0.28571429849], + [0.0, 0.28571428359, 0.28571429849], + [0.14285716414, 0.28571428359, 0.28571429849], + [0.28571429849, 0.28571428359, 0.28571429849], + [0.42857143283, 0.28571428359, 0.28571429849], + [0.0, 0.42857144773, 0.28571429849], + [0.14285713434, 0.42857144773, 0.28571429849], + [0.28571426868, 0.42857144773, 0.28571429849], + [0.0, 0.57142858207, 0.28571429849], + [0.14285713434, 0.57142858207, 0.28571429849], + [0.0, 0.71428571641, 0.28571429849], + [0.0, 0.0, 0.42857143283], + [0.14285714924, 0.0, 0.42857143283], + [0.28571428359, 0.0, 0.42857143283], + [0.42857141793, 0.0, 0.42857143283], + [0.57142855227, 0.0, 0.42857143283], + [0.0, 0.14285714924, 0.42857143283], + [0.14285713434, 0.14285714924, 0.42857143283], + [0.28571426868, 0.14285714924, 0.42857143283], + [0.42857140303, 0.14285714924, 0.42857143283], + [0.0, 0.28571428359, 0.42857143283], + [0.14285713434, 0.28571428359, 0.42857143283], + [0.28571426868, 0.28571428359, 0.42857143283], + [0.0, 0.42857141793, 0.42857143283], + [0.14285713434, 0.42857141793, 0.42857143283], + [0.0, 0.57142855227, 0.42857143283], + [0.0, 0.0, 0.57142859697], + [0.14285713434, 0.0, 0.57142859697], + [0.28571429849, 0.0, 0.57142859697], + [0.42857143283, 0.0, 0.57142859697], + [0.0, 0.14285713434, 0.57142859697], + [0.14285716414, 0.14285713434, 0.57142859697], + [0.28571429849, 0.14285713434, 0.57142859697], + [0.0, 0.28571429849, 0.57142859697], + [0.14285713434, 0.28571429849, 0.57142859697], + [0.0, 0.42857143283, 0.57142859697], + [0.0, 0.0, 0.71428573132], + [0.14285713434, 0.0, 0.71428573132], + [0.28571426868, 0.0, 0.71428573132], + [0.0, 0.14285713434, 0.71428573132], + [0.14285713434, 0.14285713434, 0.71428573132], + [0.0, 0.28571426868, 0.71428573132], + [0.0, 0.0, 0.85714286566], + [0.14285716414, 0.0, 0.85714286566], + [0.0, 0.14285716414, 0.85714286566], + ], + [ + ( + "VTK_LAGRANGE_TETRAHEDRON", + [ + [ + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 16, + 22, + 27, + 31, + 34, + 36, + 35, + 32, + 28, + 23, + 17, + 10, + 37, + 65, + 86, + 101, + 111, + 117, + 43, + 70, + 90, + 104, + 113, + 118, + 64, + 85, + 100, + 110, + 116, + 119, + 38, + 42, + 112, + 39, + 40, + 41, + 69, + 89, + 103, + 102, + 87, + 66, + 67, + 68, + 88, + 63, + 115, + 49, + 84, + 99, + 109, + 107, + 94, + 75, + 54, + 58, + 61, + 82, + 97, + 79, + 44, + 114, + 62, + 71, + 91, + 105, + 108, + 98, + 83, + 59, + 55, + 50, + 76, + 95, + 80, + 11, + 33, + 15, + 18, + 24, + 29, + 30, + 26, + 21, + 14, + 13, + 12, + 19, + 25, + 20, + 45, + 48, + 60, + 106, + 46, + 47, + 53, + 57, + 56, + 51, + 72, + 92, + 74, + 93, + 81, + 96, + 73, + 78, + 77, + 52, + ] + ], + ) + ], +) + def add_point_data(mesh, dim, num_tags=2, seed=0, dtype=float): np.random.seed(seed) @@ -429,6 +702,8 @@ def cell_sorter(cell): for face_in, face_out in zip(c_in, c_out): assert np.allclose(face_in, face_out, atol=atol, rtol=0.0) else: + print("a", cells0.data) + print("b", cells1.data) assert np.array_equal(cells0.data, cells1.data) for key in input_mesh.point_data.keys(): diff --git a/tests/test_abaqus.py b/tests/test_abaqus.py index 92b325c5a..8c55f49d5 100644 --- a/tests/test_abaqus.py +++ b/tests/test_abaqus.py @@ -1,12 +1,13 @@ import pathlib import tempfile -import helpers import numpy as np import pytest import meshio +from . import helpers + @pytest.mark.parametrize( "mesh", diff --git a/tests/test_ansys.py b/tests/test_ansys.py index d6eced4b3..9be47876d 100644 --- a/tests/test_ansys.py +++ b/tests/test_ansys.py @@ -1,8 +1,9 @@ -import helpers import pytest import meshio +from . import helpers + @pytest.mark.parametrize( "mesh", diff --git a/tests/test_avsucd.py b/tests/test_avsucd.py index 732be8d11..0c76602b0 100644 --- a/tests/test_avsucd.py +++ b/tests/test_avsucd.py @@ -1,8 +1,9 @@ -import helpers import pytest import meshio +from . import helpers + @pytest.mark.parametrize( "mesh", diff --git a/tests/test_cgns.py b/tests/test_cgns.py index ddac5c2e1..38b16dcfe 100644 --- a/tests/test_cgns.py +++ b/tests/test_cgns.py @@ -1,8 +1,9 @@ -import helpers import pytest import meshio +from . import helpers + @pytest.mark.parametrize( "mesh", diff --git a/tests/test_cli.py b/tests/test_cli.py index 8160e1bf9..83cfc72f5 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -1,11 +1,12 @@ import tempfile from pathlib import Path -import helpers import numpy as np import meshio +from . import helpers + def is_same_mesh(mesh0, mesh1, atol): if not np.allclose(mesh0.points, mesh1.points, atol=atol, rtol=0.0): diff --git a/tests/test_dolfin.py b/tests/test_dolfin.py index 8ba13d72a..39f625c36 100644 --- a/tests/test_dolfin.py +++ b/tests/test_dolfin.py @@ -1,9 +1,10 @@ -import helpers import numpy as np import pytest import meshio +from . import helpers + @pytest.mark.parametrize( "mesh", diff --git a/tests/test_exodus.py b/tests/test_exodus.py index 3ab500c6f..fa61db49b 100644 --- a/tests/test_exodus.py +++ b/tests/test_exodus.py @@ -1,8 +1,9 @@ -import helpers import pytest import meshio +from . import helpers + netCDF4 = pytest.importorskip("netCDF4") test_set = [ diff --git a/tests/test_flac3d.py b/tests/test_flac3d.py index 0326b242d..f4a5167f3 100644 --- a/tests/test_flac3d.py +++ b/tests/test_flac3d.py @@ -2,12 +2,13 @@ import pathlib import sys -import helpers import numpy as np import pytest import meshio +from . import helpers + @pytest.mark.parametrize( "mesh, data", diff --git a/tests/test_gmsh.py b/tests/test_gmsh.py index 2132af43f..25a5c9f38 100644 --- a/tests/test_gmsh.py +++ b/tests/test_gmsh.py @@ -2,12 +2,13 @@ import pathlib from functools import partial -import helpers import numpy as np import pytest import meshio +from . import helpers + def gmsh_periodic(): mesh = copy.deepcopy(helpers.quad_mesh) diff --git a/tests/test_hmf.py b/tests/test_hmf.py index ea822e691..b29f5eb0e 100644 --- a/tests/test_hmf.py +++ b/tests/test_hmf.py @@ -1,9 +1,10 @@ -import helpers import numpy as np import pytest import meshio +from . import helpers + test_set_full = [ helpers.empty_mesh, helpers.line_mesh, diff --git a/tests/test_mdpa.py b/tests/test_mdpa.py index 17ed17115..f0a79226f 100644 --- a/tests/test_mdpa.py +++ b/tests/test_mdpa.py @@ -1,8 +1,9 @@ -import helpers import pytest import meshio +from . import helpers + @pytest.mark.parametrize( "mesh", diff --git a/tests/test_med.py b/tests/test_med.py index 7c65d5df7..d51fe3c00 100644 --- a/tests/test_med.py +++ b/tests/test_med.py @@ -1,11 +1,12 @@ import pathlib -import helpers import numpy as np import pytest import meshio +from . import helpers + h5py = pytest.importorskip("h5py") diff --git a/tests/test_medit.py b/tests/test_medit.py index 2d69030f5..f34bafa2d 100644 --- a/tests/test_medit.py +++ b/tests/test_medit.py @@ -1,11 +1,12 @@ import pathlib -import helpers import numpy as np import pytest import meshio +from . import helpers + @pytest.mark.parametrize( "mesh", diff --git a/tests/test_mesh.py b/tests/test_mesh.py index 757289947..b8f550a93 100644 --- a/tests/test_mesh.py +++ b/tests/test_mesh.py @@ -1,11 +1,12 @@ import copy -import helpers import numpy as np import pytest import meshio +from . import helpers + def test_public_attributes(): # Just make sure this is here diff --git a/tests/test_moab.py b/tests/test_moab.py index bc1e4d009..a87f7bfe8 100644 --- a/tests/test_moab.py +++ b/tests/test_moab.py @@ -1,8 +1,9 @@ -import helpers import pytest import meshio +from . import helpers + h5py = pytest.importorskip("h5py") diff --git a/tests/test_nastran.py b/tests/test_nastran.py index 5a4f6ce7f..0b6081e10 100644 --- a/tests/test_nastran.py +++ b/tests/test_nastran.py @@ -1,12 +1,13 @@ import io import pathlib -import helpers import numpy as np import pytest import meshio +from . import helpers + @pytest.mark.parametrize( "mesh", diff --git a/tests/test_netgen.py b/tests/test_netgen.py index f562b10c8..090c27580 100644 --- a/tests/test_netgen.py +++ b/tests/test_netgen.py @@ -1,8 +1,9 @@ -import helpers import pytest import meshio +from . import helpers + test_set = [ helpers.empty_mesh, helpers.line_mesh, diff --git a/tests/test_neuroglancer.py b/tests/test_neuroglancer.py index 07ba06113..25d15ce0d 100644 --- a/tests/test_neuroglancer.py +++ b/tests/test_neuroglancer.py @@ -1,11 +1,12 @@ import pathlib -import helpers import numpy as np import pytest import meshio +from . import helpers + @pytest.mark.parametrize( "mesh", diff --git a/tests/test_obj.py b/tests/test_obj.py index 4aac690b3..f35580675 100644 --- a/tests/test_obj.py +++ b/tests/test_obj.py @@ -1,11 +1,12 @@ import pathlib -import helpers import numpy as np import pytest import meshio +from . import helpers + @pytest.mark.parametrize( "mesh", diff --git a/tests/test_off.py b/tests/test_off.py index a6622347f..df058f67b 100644 --- a/tests/test_off.py +++ b/tests/test_off.py @@ -1,8 +1,9 @@ -import helpers import pytest import meshio +from . import helpers + @pytest.mark.parametrize( "mesh", diff --git a/tests/test_permas.py b/tests/test_permas.py index 1a26f5ccb..bf3ec1bf5 100644 --- a/tests/test_permas.py +++ b/tests/test_permas.py @@ -1,8 +1,9 @@ -import helpers import pytest import meshio +from . import helpers + @pytest.mark.parametrize( "mesh", diff --git a/tests/test_ply.py b/tests/test_ply.py index dde35d78f..f27418e3c 100644 --- a/tests/test_ply.py +++ b/tests/test_ply.py @@ -1,11 +1,12 @@ import pathlib -import helpers import numpy as np import pytest import meshio +from . import helpers + @pytest.mark.parametrize( "mesh", diff --git a/tests/test_stl.py b/tests/test_stl.py index 4c4b0c372..3b5ef44aa 100644 --- a/tests/test_stl.py +++ b/tests/test_stl.py @@ -1,8 +1,9 @@ -import helpers import pytest import meshio +from . import helpers + @pytest.mark.parametrize( "mesh", diff --git a/tests/test_su2.py b/tests/test_su2.py index 077356ba0..eb4e49c64 100644 --- a/tests/test_su2.py +++ b/tests/test_su2.py @@ -1,11 +1,12 @@ import pathlib -import helpers import numpy as np import pytest import meshio +from . import helpers + test_set = [ # helpers.empty_mesh, helpers.tri_mesh_2d, diff --git a/tests/test_svg.py b/tests/test_svg.py index e3888d2e6..fbcfd5ca5 100644 --- a/tests/test_svg.py +++ b/tests/test_svg.py @@ -1,11 +1,12 @@ import pathlib import tempfile -import helpers import pytest import meshio +from . import helpers + test_set = [ helpers.empty_mesh, helpers.line_mesh, diff --git a/tests/test_tecplot.py b/tests/test_tecplot.py index c260ed69a..781ffec50 100644 --- a/tests/test_tecplot.py +++ b/tests/test_tecplot.py @@ -1,12 +1,13 @@ import pathlib from copy import deepcopy -import helpers import numpy as np import pytest import meshio +from . import helpers + @pytest.mark.parametrize( "mesh", diff --git a/tests/test_tetgen.py b/tests/test_tetgen.py index 931d38179..360083709 100644 --- a/tests/test_tetgen.py +++ b/tests/test_tetgen.py @@ -1,10 +1,11 @@ import pathlib -import helpers import pytest import meshio +from . import helpers + test_set = [ # helpers.empty_mesh, helpers.tet_mesh diff --git a/tests/test_ugrid.py b/tests/test_ugrid.py index 59f59b29d..d8a109442 100644 --- a/tests/test_ugrid.py +++ b/tests/test_ugrid.py @@ -1,11 +1,12 @@ import pathlib -import helpers import numpy as np import pytest import meshio +from . import helpers + this_dir = pathlib.Path(__file__).resolve().parent diff --git a/tests/test_vtk.py b/tests/test_vtk.py index f40ad13d6..b7e027b14 100644 --- a/tests/test_vtk.py +++ b/tests/test_vtk.py @@ -1,17 +1,19 @@ import pathlib from functools import partial -import helpers import numpy as np import pytest import meshio +from . import helpers + test_set = [ helpers.empty_mesh, helpers.line_mesh, helpers.tri_mesh_2d, helpers.tri_mesh, + helpers.tri_mesh_one_cell, helpers.triangle6_mesh, helpers.quad_mesh, helpers.quad8_mesh, @@ -23,6 +25,7 @@ helpers.polygon_mesh, helpers.pyramid_mesh, helpers.wedge_mesh, + # helpers.lagrange_high_order_mesh, helpers.add_point_data(helpers.tri_mesh, 1), helpers.add_point_data(helpers.tri_mesh, 2), helpers.add_point_data(helpers.tri_mesh, 3), diff --git a/tests/test_vtu.py b/tests/test_vtu.py index 4d6f078ff..02a247a1c 100644 --- a/tests/test_vtu.py +++ b/tests/test_vtu.py @@ -1,20 +1,23 @@ import pathlib -import helpers import numpy as np import pytest import meshio +from . import helpers + test_set = [ # helpers.empty_mesh, helpers.line_mesh, helpers.tri_mesh, + helpers.tri_mesh_one_cell, helpers.triangle6_mesh, helpers.quad_mesh, helpers.quad8_mesh, helpers.tri_quad_mesh, helpers.polygon_mesh, + helpers.polygon_mesh_one_cell, helpers.polygon2_mesh, helpers.tet_mesh, helpers.tet10_mesh, @@ -23,6 +26,7 @@ helpers.pyramid_mesh, helpers.wedge_mesh, helpers.polyhedron_mesh, + helpers.lagrange_high_order_mesh, helpers.add_point_data(helpers.tri_mesh, 1), helpers.add_point_data(helpers.tri_mesh, 2), helpers.add_point_data(helpers.tri_mesh, 3), diff --git a/tests/test_wkt.py b/tests/test_wkt.py index 16a3d3e09..4e8e85901 100644 --- a/tests/test_wkt.py +++ b/tests/test_wkt.py @@ -1,11 +1,12 @@ import pathlib -import helpers import numpy as np import pytest import meshio +from . import helpers + @pytest.mark.parametrize( "mesh", diff --git a/tests/test_xdmf.py b/tests/test_xdmf.py index 510ea06aa..06ba6763b 100644 --- a/tests/test_xdmf.py +++ b/tests/test_xdmf.py @@ -1,9 +1,10 @@ -import helpers import numpy as np import pytest import meshio +from . import helpers + test_set_full = [ helpers.empty_mesh, helpers.line_mesh,