diff --git a/src/io4dolfinx/backends/pyvista/backend.py b/src/io4dolfinx/backends/pyvista/backend.py index 261c47a..77372b6 100644 --- a/src/io4dolfinx/backends/pyvista/backend.py +++ b/src/io4dolfinx/backends/pyvista/backend.py @@ -22,28 +22,10 @@ from io4dolfinx.utils import check_file_exists from .. import FileMode, ReadMode +from ..vtk_cell_types import _arbitrary_lagrange_vtk, _cell_degree, _first_order_vtk # Cell types can be found at # https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html -_first_order_vtk = { - 1: "point", - 3: "interval", - 5: "triangle", - 9: "quadrilateral", - 10: "tetrahedron", - 12: "hexahedron", -} - -_arbitrary_lagrange_vtk = { - 68: "interval", - 69: "triangle", - 70: "quadrilateral", - 71: "tetrahedron", - 72: "hexahedron", - 73: "prism", - 74: "pyramid", -} - _quadratric_vtk = { 21: "interval", 22: "triangle", @@ -56,52 +38,6 @@ read_mode = ReadMode.serial -def _cell_degree(ct: str, num_nodes: int) -> int: - if ct == "point": - return 1 - elif ct == "interval": - return int(num_nodes - 1) - elif ct == "triangle": - n = (np.sqrt(1 + 8 * num_nodes) - 1) / 2 - if 2 * num_nodes != n * (n + 1): - raise ValueError(f"Unknown triangle layout. Number of nodes: {num_nodes}") - return int(n - 1) - elif ct == "tetrahedron": - n = 0 - while n * (n + 1) * (n + 2) < 6 * num_nodes: - n += 1 - if n * (n + 1) * (n + 2) != 6 * num_nodes: - raise ValueError(f"Unknown tetrahedron layout. Number of nodes: {num_nodes}") - return int(n - 1) - - elif ct == "quadrilateral": - n = np.sqrt(num_nodes) - if num_nodes != n * n: - raise ValueError(f"Unknown quadrilateral layout. Number of nodes: {num_nodes}") - return int(n - 1) - elif ct == "hexahedron": - n = np.cbrt(num_nodes) - if num_nodes != n * n * n: - raise ValueError(f"Unknown hexahedron layout. Number of nodes: {num_nodes}") - return int(n - 1) - elif ct == "prism": - if num_nodes == 6: - return 1 - elif num_nodes == 15: - return 2 - else: - raise ValueError(f"Unknown prism layout. Number of nodes: {num_nodes}") - elif ct == "pyramid": - if num_nodes == 5: - return 1 - elif num_nodes == 13: - return 2 - else: - raise ValueError(f"Unknown pyramid layout. Number of nodes: {num_nodes}") - else: - raise ValueError(f"Unknown cell type {ct} with {num_nodes=}.") - - def get_default_backend_args(arguments: dict[str, Any] | None) -> dict[str, Any]: """Get default backend arguments given a set of input arguments. diff --git a/src/io4dolfinx/backends/vtk_cell_types.py b/src/io4dolfinx/backends/vtk_cell_types.py new file mode 100644 index 0000000..aa0bf45 --- /dev/null +++ b/src/io4dolfinx/backends/vtk_cell_types.py @@ -0,0 +1,65 @@ +import numpy as np + +# Cell types can be found at +# https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html + +_first_order_vtk = { + 1: "point", + 3: "interval", + 5: "triangle", + 9: "quadrilateral", + 10: "tetrahedron", + 12: "hexahedron", +} + +_arbitrary_lagrange_vtk = { + 68: "interval", + 69: "triangle", + 70: "quadrilateral", + 71: "tetrahedron", + 72: "hexahedron", + 73: "prism", + 74: "pyramid", +} + + +def _cell_degree(cell_type: str, num_nodes: int) -> int: + if cell_type == "point": + return 1 + elif cell_type == "interval": + return int(num_nodes - 1) + elif cell_type == "triangle": + degree = (np.sqrt(1 + 8 * num_nodes) - 1) / 2 + if 2 * num_nodes != degree * (degree + 1): + raise ValueError(f"Unknown triangle layout. Number of nodes: {num_nodes}") + return int(degree - 1) + elif cell_type == "tetrahedron": + degree = 0 + while degree * (degree + 1) * (degree + 2) < 6 * num_nodes: + degree += 1 + if degree * (degree + 1) * (degree + 2) != 6 * num_nodes: + raise ValueError(f"Unknown tetrahedron layout. Number of nodes: {num_nodes}") + return int(degree - 1) + elif cell_type == "quadrilateral": + degree = np.sqrt(num_nodes) + if num_nodes != degree * degree: + raise ValueError(f"Unknown quadrilateral layout. Number of nodes: {num_nodes}") + return int(degree - 1) + elif cell_type == "hexahedron": + degree = np.cbrt(num_nodes) + if num_nodes != degree * degree * degree: + raise ValueError(f"Unknown hexahedron layout. Number of nodes: {num_nodes}") + return int(degree - 1) + elif cell_type == "prism": + if num_nodes == 6: + return 1 + elif num_nodes == 15: + return 2 + raise ValueError(f"Unknown prism layout. Number of nodes: {num_nodes}") + elif cell_type == "pyramid": + if num_nodes == 5: + return 1 + elif num_nodes == 13: + return 2 + raise ValueError(f"Unknown pyramid layout. Number of nodes: {num_nodes}") + raise ValueError(f"Unknown cell type {cell_type} with {num_nodes=}.") diff --git a/src/io4dolfinx/backends/vtkhdf/backend.py b/src/io4dolfinx/backends/vtkhdf/backend.py index e4f1737..f3297d8 100644 --- a/src/io4dolfinx/backends/vtkhdf/backend.py +++ b/src/io4dolfinx/backends/vtkhdf/backend.py @@ -18,7 +18,7 @@ from .. import FileMode, ReadMode from ..h5py.backend import convert_file_mode, h5pyfile -from ..pyvista.backend import _arbitrary_lagrange_vtk, _cell_degree, _first_order_vtk +from ..vtk_cell_types import _arbitrary_lagrange_vtk, _cell_degree, _first_order_vtk read_mode = ReadMode.parallel @@ -349,7 +349,7 @@ def write_attributes( attributes: Dictionary of attributes to write backend_args: Arguments to backend """ - raise NotImplementedError("The Pyvista backend cannot write attributes.") + raise NotImplementedError("The VTKHDF backend cannot write attributes.") def read_attributes( @@ -369,7 +369,7 @@ def read_attributes( Returns: Dictionary of attributes read from file """ - raise NotImplementedError("The Pyvista backend cannot read attributes.") + raise NotImplementedError("The VTKHDF backend cannot read attributes.") def read_timestamps( @@ -1000,7 +1000,7 @@ def read_dofmap( Returns: Dofmap as an AdjacencyList """ - raise NotImplementedError("The Pyvista backend cannot make checkpoints.") + raise NotImplementedError("The VTKHDF backend cannot make checkpoints.") def read_dofs( @@ -1024,7 +1024,7 @@ def read_dofs( and the global starting point on the process. Process 0 has [0, M), process 1 [M, N), process 2 [N, O) etc. """ - raise NotImplementedError("The Pyvista backend cannot make checkpoints.") + raise NotImplementedError("The VTKHDF backend cannot make checkpoints.") def read_cell_perms(