From e553e7b66908350f3c28752711d0aad6c07cb391 Mon Sep 17 00:00:00 2001 From: Lars Kaiser <71793357+kaiserls@users.noreply.github.com> Date: Fri, 26 Apr 2024 01:41:06 +0200 Subject: [PATCH 1/3] [DGF] Implement the dgf mesh format --- src/meshio/__init__.py | 1 + src/meshio/dgf/__init__.py | 3 + src/meshio/dgf/_dgf.py | 272 +++++++++++++++++++++++++++++++++++++ 3 files changed, 276 insertions(+) create mode 100644 src/meshio/dgf/__init__.py create mode 100644 src/meshio/dgf/_dgf.py diff --git a/src/meshio/__init__.py b/src/meshio/__init__.py index 444e232ab..ea05c8916 100644 --- a/src/meshio/__init__.py +++ b/src/meshio/__init__.py @@ -4,6 +4,7 @@ ansys, avsucd, cgns, + dgf, dolfin, exodus, flac3d, diff --git a/src/meshio/dgf/__init__.py b/src/meshio/dgf/__init__.py new file mode 100644 index 000000000..04645ceda --- /dev/null +++ b/src/meshio/dgf/__init__.py @@ -0,0 +1,3 @@ +from ._dgf import read, write + +__all__ = ["read", "write"] diff --git a/src/meshio/dgf/_dgf.py b/src/meshio/dgf/_dgf.py new file mode 100644 index 000000000..13be441f8 --- /dev/null +++ b/src/meshio/dgf/_dgf.py @@ -0,0 +1,272 @@ +""" +I/O for DUNE grid format, cf. +https://www.dune-project.org/doxygen/2.9.0/group__DuneGridFormatParser.html#details +""" + +import numpy as np + +from ..__about__ import __version__ as version +from .._common import warn +from .._exceptions import ReadError, WriteError +from .._files import open_file +from .._helpers import register_format +from .._mesh import Mesh, CellBlock + + +#: Maps DGF cell types to meshio cell types. +dgf_to_meshio_type = { + ("SIMPLEX",2): "line", + ("SIMPLEX",3): "triangle", + ("CUBE",4): "quad", +} + +#: Maps DGF cell types to meshio cell types. +meshio_to_dgf_type = { + "line": ("SIMPLEX",2), + "triangle": ("SIMPLEX",3), + "quad": ("CUBE",4), +} + +#: The DGF keywords that are supported by meshio. +dgf_keywords = [ + # Directly supported + "VERTEX", + "SIMPLEX", + "CUBE", + # TODO: Add support later! + # Supported by storing the information on the mesh? No face information in meshio? + # "BOUNDARYSEGMENTS", + # "BOUNDARYDOMAIN", + # Unsupported keywords + # "SIMPLEXGENERATOR" + # "GRIDPARAMETER", + # "PERIODICFACETRANSFORMATION", + # ... +] + + +def read(filename): + """Reads a DGF file. + + Args: + filename (PathLike): The path to the file. + + Returns: + meshio.Mesh: The meshio mesh + """ + with open_file(filename, "r") as f: + out = read_buffer(f) + return out + + +def readline(f): + """Reads a line from a dgf file and skips comments. + + Args: + f: The file object. + + Returns: + str: the first non-comment line, stripped of leading/trailing whitespace, and with end-of-line comment removed. + """ + # check file is at end + while True: + line = f.readline() + if not line: + return None + line = line.split("%")[0].strip() + if line: + return line + +def skip_block(f): + """Skips a block in a dgf file. + + Args: + f: The file object. + """ + line = readline(f) + while line and not line.startswith("#"): + line = readline(f) + return + +def append_cell_block(cells, cell_block, cell_data, block_data): + """Appends a cell block to the mesh. + + Args: + cells (list): The current list of cell blocks. + cell_block (CellBlock): The cell block to append. + cell_data (dict): The current cell data. + block_data (dict): The block data to append. + """ + cells.append(cell_block) + + for varname, variable in block_data.items(): + if varname in cell_data: + cell_data[varname].append(variable) + else: + cell_data[varname] = [variable] + +def read_buffer(f): + points = [] + point_data = {} + cells = [] + cell_data = {} + + line = readline(f) + if not line.startswith("DGF"): + raise ReadError("DGF Files must start with 'DGF'") + + line = readline(f) + while line: + if line.upper().startswith("VERTEX"): + points, point_data = _read_vertex(f) + elif line.upper().startswith("SIMPLEX"): + cell_block, block_data = _read_cells(f, "SIMPLEX") + append_cell_block(cells, cell_block, cell_data, block_data) + elif line.upper().startswith("CUBE"): + cell_block, block_data = _read_cells(f, "CUBE") + append_cell_block(cells, cell_block, cell_data, block_data) + elif line.upper().startswith("INTERVAL"): + raise ReadError("INTERVAL blocks are not supported by meshio, please use meshzoo to generate the grid.") + elif line.startswith("#"): + warn("Comment with # found. Skipping line.") + else: + warn(f"Unknown keyword in '{line}'. Skipping block.") + skip_block(f) + line = readline(f) + + return Mesh(points, cells, point_data, cell_data) + +def _read_vertex(f): + """Reads the vertex block of a DGF file. + + Args: + f: The file object, already in the start of the vertex block. + + Returns: + tuple[np.ndarray, dict]: The point coordinates and point data. + """ + data = [] + num_parameters = 0 + points = [] + point_data = {} + + # Determine number of parameters defined on the points + line = readline(f) + if line.startswith("parameters"): + num_parameters = int(line.split()[1]) + line = readline(f) + else: + num_parameters = 0 + + entries_per_line = len(line.split()) + n_dim = entries_per_line - num_parameters + + # Read vertex block including point coordinates and point data + while line: + if line.startswith("#"): + break + data.append([float(x) for x in line.split()]) + line = readline(f) + + # Extract point coordinates and point data + points = np.array([x[:n_dim] for x in data]) + if num_parameters > 0: + point_data = {f"param_{i}": np.array([x[n_dim + i] for x in data],dtype=float) for i in range(num_parameters)} + + return points, point_data + +def _read_cells(f, block_name): + """Reads the simplex block of a DGF file and returns the cells and cell data. + + Args: + f: The file object, already in the start of the simplex block. + block_name (str): The name of the block. + + Returns: + tuple[np.ndarray, dict]: The cells and cell data. + """ + + data = [] + num_parameters = 0 + + # Determine number of parameters defined on the simplex elements + line = readline(f) + if line.startswith("parameters"): + num_parameters = int(line.split()[1]) + line = readline(f) + else: + num_parameters = 0 + + entries_per_cell_definition = len(line.split()) + n_indices = entries_per_cell_definition - num_parameters + + # Read simplex block including cell connectivity and cell data + while line: + if line.startswith("#"): + break + data.append([float(x) for x in line.split()]) + line = readline(f) + + # Determine cell type + if block_name == "SIMPLEX": + dgf_type = ("SIMPLEX", n_indices) + elif block_name == "CUBE": + dgf_type = ("CUBE", n_indices) + + # Extract cell connectivity and cell data + cell_connectivity = np.array([x[:n_indices] for x in data],dtype=int) + cell_block = CellBlock(dgf_to_meshio_type[dgf_type], cell_connectivity) + cell_data = {} + if num_parameters > 0: + cell_data = {f"param_{i}": np.array([x[n_indices + i] for x in data], dtype=float) for i in range(num_parameters)} + + return cell_block, cell_data + +def write(filename, mesh): + """Writes a DGF file.""" + + with open_file(filename, "w") as f: + f.write("DGF\n") + f.write(f'% "Written by meshio v{version}"\n') + + # Write vertex block + f.write("VERTEX\n") + if len(mesh.point_data) > 0: + f.write(f'parameters {len(mesh.point_data)}\n') + for point, data in zip(mesh.points, mesh.point_data.values()): + f.write(" ".join(str(c) for c in point) + " " + " ".join(str(d) for d in data) + "\n") + else: + for point in mesh.points: + f.write(" ".join(str(c) for c in point) + "\n") + f.write("#\n") + + # Write cell blocks + number_of_parameters = len(mesh.cell_data) + for cell_block_index, cell_block in enumerate(mesh.cells): + dgf_cell_type, _ = meshio_to_dgf_type[cell_block.type] + f.write(dgf_cell_type + "\n") + if len(mesh.cell_data) > 0: + f.write(f"parameters {number_of_parameters}\n") + + # grab all param values for the current cell block + # we cant do it directly on the existing list, because it does not support slicing + n_cells = len(cell_block.data) + parameters = np.zeros((number_of_parameters, n_cells)) + for pi, (_, parameter_all_blocks) in enumerate(mesh.cell_data.items()): + parameters[pi,:] = parameter_all_blocks[cell_block_index] + + # write out all cell indices and + for cell_index, cell in enumerate(cell_block.data): + f.write(" ".join(str(id) for id in cell) + " " + " ".join(str(d) for d in parameters[:,cell_index] ) + "\n") + else: + for cell in cell_block.data: + f.write(" ".join(str(id) for id in cell) + "\n") + f.write("#\n") + + # Write default boundary domain + f.write("BOUNDARYDOMAIN\n") + f.write("default 1\n") + f.write("#\n") + + +register_format("dgf", [".dgf"], read, {"dgf": write}) From 8fc6725048e89d295de2b8947606ceae8b8f620d Mon Sep 17 00:00:00 2001 From: Lars Kaiser <71793357+kaiserls@users.noreply.github.com> Date: Fri, 26 Apr 2024 01:42:57 +0200 Subject: [PATCH 2/3] [DGF] Add test for dgf format --- tests/meshes/dgf/cube.dgf | 26 ++++++++++++++++++++++++++ tests/meshes/dgf/line.dgf | 14 ++++++++++++++ tests/meshes/dgf/triangle.dgf | 29 +++++++++++++++++++++++++++++ tests/test_dgf.py | 24 ++++++++++++++++++++++++ 4 files changed, 93 insertions(+) create mode 100644 tests/meshes/dgf/cube.dgf create mode 100644 tests/meshes/dgf/line.dgf create mode 100644 tests/meshes/dgf/triangle.dgf create mode 100644 tests/test_dgf.py diff --git a/tests/meshes/dgf/cube.dgf b/tests/meshes/dgf/cube.dgf new file mode 100644 index 000000000..3eeb58fbf --- /dev/null +++ b/tests/meshes/dgf/cube.dgf @@ -0,0 +1,26 @@ +DGF +% SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root +% SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception +Vertex % the vertices of the grid +-1 -1 % vertex 0 +-0.2 -1 % vertex 1 + 1 -1 % vertex 2 + 1 -0.6 % vertex 3 + 1 1 % vertex 4 +-1 1 % vertex 5 +0.2 -0.2 % vertex 6 +# +CUBE % a cube grid +0 1 5 6 % cube 0, vertices 0,1,5,6 +1 2 6 3 % cube 1 +6 3 5 4 % cube 2 +# +BOUNDARYSEGMENTS +2 1 2 % between vertex 1,2 use id 2 +2 2 3 % between vertex 2,3 use id 2 +4 3 4 % between vertex 3,4 use id 4 +# +BOUNDARYDOMAIN +default 1 % all other boundary segments have id 1 +# +# examplegrid1c.dgf \ No newline at end of file diff --git a/tests/meshes/dgf/line.dgf b/tests/meshes/dgf/line.dgf new file mode 100644 index 000000000..36390e7a7 --- /dev/null +++ b/tests/meshes/dgf/line.dgf @@ -0,0 +1,14 @@ +DGF +Vertex +0.0 -0.002 0.0 +0.0 0.0 0.0 +-5.29908e-05 0.0019993000000000003 0.0 +# +SIMPLEX +parameters 1 +0 1 0.000209184 +1 2 0.000207817 +# +BOUNDARYDOMAIN +default 1 +# \ No newline at end of file diff --git a/tests/meshes/dgf/triangle.dgf b/tests/meshes/dgf/triangle.dgf new file mode 100644 index 000000000..a0f66a87a --- /dev/null +++ b/tests/meshes/dgf/triangle.dgf @@ -0,0 +1,29 @@ +DGF +% SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root +% SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception +Vertex % the vertices of the grid +-1 -1 % vertex 0 +-0.2 -1 % vertex 1 + 1 -1 % vertex 2 + 1 -0.6 % vertex 3 + 1 1 % vertex 4 +-1 1 % vertex 5 +0.2 -0.2 % vertex 6 +# +SIMPLEX % a simplex grid +0 1 5 % triangle 0, vertices 0,1,5 +1 3 6 % triangle 1 +1 2 3 % triangle 2 +6 3 4 % triangle 3 +6 4 5 % triangle 4 +6 5 1 % triangle 5 +# +BOUNDARYSEGMENTS +2 1 2 % between vertex 1,2 use id 2 +2 2 3 % between vertex 2,3 use id 2 +4 3 4 % between vertex 3,4 use id 4 +# +BOUNDARYDOMAIN +default 1 % all other boundary segments have id 1 +# +# examplegrid1s.dgf diff --git a/tests/test_dgf.py b/tests/test_dgf.py new file mode 100644 index 000000000..ef66a1414 --- /dev/null +++ b/tests/test_dgf.py @@ -0,0 +1,24 @@ +import pathlib +import pytest +import meshio + +from . import helpers + +@pytest.mark.parametrize( + "filename, ref_cell_type, ref_num_cells, ref_num_points, ref_dim", + [ + ("line.dgf", "line", 2, 3, 3), + ("triangle.dgf", "triangle", 6, 7, 2), + ("cube.dgf", "quad", 3, 7, 2), + ], +) +def test_dgf(filename, ref_cell_type, ref_num_cells, ref_num_points, ref_dim, tmp_path): + this_dir = pathlib.Path(__file__).resolve().parent + filename = this_dir / "meshes" / "dgf" / filename + mesh = meshio.read(filename) + assert len(mesh.cells) == 1 + assert ref_cell_type == mesh.cells[0].type + assert len(mesh.cells[0].data) == ref_num_cells + assert len(mesh.points) == ref_num_points + assert mesh.points.shape[1] == ref_dim + helpers.write_read(tmp_path, meshio.dgf.write, meshio.dgf.read, mesh, 1.0e-15, extension=".dgf") From 50fdb8c24ad8bff393316b0ba5f0bb277cb91623 Mon Sep 17 00:00:00 2001 From: Lars Kaiser <71793357+kaiserls@users.noreply.github.com> Date: Sun, 28 Apr 2024 16:58:14 +0200 Subject: [PATCH 3/3] [DGF] Formatting and fixing import --- src/meshio/__init__.py | 1 + src/meshio/dgf/_dgf.py | 75 +++++++++++++++++++++++++++--------------- tests/test_dgf.py | 7 +++- 3 files changed, 56 insertions(+), 27 deletions(-) diff --git a/src/meshio/__init__.py b/src/meshio/__init__.py index ea05c8916..ebedd108c 100644 --- a/src/meshio/__init__.py +++ b/src/meshio/__init__.py @@ -49,6 +49,7 @@ "ansys", "avsucd", "cgns", + "dgf", "dolfin", "exodus", "flac3d", diff --git a/src/meshio/dgf/_dgf.py b/src/meshio/dgf/_dgf.py index 13be441f8..b99c0b699 100644 --- a/src/meshio/dgf/_dgf.py +++ b/src/meshio/dgf/_dgf.py @@ -7,24 +7,23 @@ from ..__about__ import __version__ as version from .._common import warn -from .._exceptions import ReadError, WriteError +from .._exceptions import ReadError from .._files import open_file from .._helpers import register_format -from .._mesh import Mesh, CellBlock - +from .._mesh import CellBlock, Mesh #: Maps DGF cell types to meshio cell types. dgf_to_meshio_type = { - ("SIMPLEX",2): "line", - ("SIMPLEX",3): "triangle", - ("CUBE",4): "quad", + ("SIMPLEX", 2): "line", + ("SIMPLEX", 3): "triangle", + ("CUBE", 4): "quad", } #: Maps DGF cell types to meshio cell types. meshio_to_dgf_type = { - "line": ("SIMPLEX",2), - "triangle": ("SIMPLEX",3), - "quad": ("CUBE",4), + "line": ("SIMPLEX", 2), + "triangle": ("SIMPLEX", 3), + "quad": ("CUBE", 4), } #: The DGF keywords that are supported by meshio. @@ -77,6 +76,7 @@ def readline(f): if line: return line + def skip_block(f): """Skips a block in a dgf file. @@ -88,6 +88,7 @@ def skip_block(f): line = readline(f) return + def append_cell_block(cells, cell_block, cell_data, block_data): """Appends a cell block to the mesh. @@ -105,6 +106,7 @@ def append_cell_block(cells, cell_block, cell_data, block_data): else: cell_data[varname] = [variable] + def read_buffer(f): points = [] point_data = {} @@ -126,7 +128,9 @@ def read_buffer(f): cell_block, block_data = _read_cells(f, "CUBE") append_cell_block(cells, cell_block, cell_data, block_data) elif line.upper().startswith("INTERVAL"): - raise ReadError("INTERVAL blocks are not supported by meshio, please use meshzoo to generate the grid.") + raise ReadError( + "INTERVAL blocks are not supported by meshio, please use meshzoo to generate the grid." + ) elif line.startswith("#"): warn("Comment with # found. Skipping line.") else: @@ -136,6 +140,7 @@ def read_buffer(f): return Mesh(points, cells, point_data, cell_data) + def _read_vertex(f): """Reads the vertex block of a DGF file. @@ -157,7 +162,7 @@ def _read_vertex(f): line = readline(f) else: num_parameters = 0 - + entries_per_line = len(line.split()) n_dim = entries_per_line - num_parameters @@ -167,21 +172,25 @@ def _read_vertex(f): break data.append([float(x) for x in line.split()]) line = readline(f) - + # Extract point coordinates and point data points = np.array([x[:n_dim] for x in data]) if num_parameters > 0: - point_data = {f"param_{i}": np.array([x[n_dim + i] for x in data],dtype=float) for i in range(num_parameters)} - + point_data = { + f"param_{i}": np.array([x[n_dim + i] for x in data], dtype=float) + for i in range(num_parameters) + } + return points, point_data + def _read_cells(f, block_name): """Reads the simplex block of a DGF file and returns the cells and cell data. Args: f: The file object, already in the start of the simplex block. block_name (str): The name of the block. - + Returns: tuple[np.ndarray, dict]: The cells and cell data. """ @@ -196,7 +205,7 @@ def _read_cells(f, block_name): line = readline(f) else: num_parameters = 0 - + entries_per_cell_definition = len(line.split()) n_indices = entries_per_cell_definition - num_parameters @@ -214,13 +223,17 @@ def _read_cells(f, block_name): dgf_type = ("CUBE", n_indices) # Extract cell connectivity and cell data - cell_connectivity = np.array([x[:n_indices] for x in data],dtype=int) + cell_connectivity = np.array([x[:n_indices] for x in data], dtype=int) cell_block = CellBlock(dgf_to_meshio_type[dgf_type], cell_connectivity) cell_data = {} if num_parameters > 0: - cell_data = {f"param_{i}": np.array([x[n_indices + i] for x in data], dtype=float) for i in range(num_parameters)} + cell_data = { + f"param_{i}": np.array([x[n_indices + i] for x in data], dtype=float) + for i in range(num_parameters) + } + + return cell_block, cell_data - return cell_block, cell_data def write(filename, mesh): """Writes a DGF file.""" @@ -228,13 +241,18 @@ def write(filename, mesh): with open_file(filename, "w") as f: f.write("DGF\n") f.write(f'% "Written by meshio v{version}"\n') - + # Write vertex block f.write("VERTEX\n") if len(mesh.point_data) > 0: - f.write(f'parameters {len(mesh.point_data)}\n') + f.write(f"parameters {len(mesh.point_data)}\n") for point, data in zip(mesh.points, mesh.point_data.values()): - f.write(" ".join(str(c) for c in point) + " " + " ".join(str(d) for d in data) + "\n") + f.write( + " ".join(str(c) for c in point) + + " " + + " ".join(str(d) for d in data) + + "\n" + ) else: for point in mesh.points: f.write(" ".join(str(c) for c in point) + "\n") @@ -253,16 +271,21 @@ def write(filename, mesh): n_cells = len(cell_block.data) parameters = np.zeros((number_of_parameters, n_cells)) for pi, (_, parameter_all_blocks) in enumerate(mesh.cell_data.items()): - parameters[pi,:] = parameter_all_blocks[cell_block_index] - + parameters[pi, :] = parameter_all_blocks[cell_block_index] + # write out all cell indices and for cell_index, cell in enumerate(cell_block.data): - f.write(" ".join(str(id) for id in cell) + " " + " ".join(str(d) for d in parameters[:,cell_index] ) + "\n") + f.write( + " ".join(str(id) for id in cell) + + " " + + " ".join(str(d) for d in parameters[:, cell_index]) + + "\n" + ) else: for cell in cell_block.data: f.write(" ".join(str(id) for id in cell) + "\n") f.write("#\n") - + # Write default boundary domain f.write("BOUNDARYDOMAIN\n") f.write("default 1\n") diff --git a/tests/test_dgf.py b/tests/test_dgf.py index ef66a1414..b96d7b5a4 100644 --- a/tests/test_dgf.py +++ b/tests/test_dgf.py @@ -1,9 +1,12 @@ import pathlib + import pytest + import meshio from . import helpers + @pytest.mark.parametrize( "filename, ref_cell_type, ref_num_cells, ref_num_points, ref_dim", [ @@ -21,4 +24,6 @@ def test_dgf(filename, ref_cell_type, ref_num_cells, ref_num_points, ref_dim, tm assert len(mesh.cells[0].data) == ref_num_cells assert len(mesh.points) == ref_num_points assert mesh.points.shape[1] == ref_dim - helpers.write_read(tmp_path, meshio.dgf.write, meshio.dgf.read, mesh, 1.0e-15, extension=".dgf") + helpers.write_read( + tmp_path, meshio.dgf.write, meshio.dgf.read, mesh, 1.0e-15, extension=".dgf" + )