diff --git a/pyproject.toml b/pyproject.toml index f87879b..98787cd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -13,7 +13,7 @@ description = "Convert OpenFOAM files to dolfinx functions" readme = "README/md" requires-python = ">=3.9" license = { file = "LICENSE" } -dependencies = ["fenics-dolfinx>=0.9.0", "pyvista"] +dependencies = ["fenics-dolfinx>=0.9.0", "pyvista", "scipy"] classifiers = [ "Natural Language :: English", "Topic :: Scientific/Engineering", diff --git a/src/foam2dolfinx/helpers.py b/src/foam2dolfinx/helpers.py new file mode 100644 index 0000000..331dcdb --- /dev/null +++ b/src/foam2dolfinx/helpers.py @@ -0,0 +1,42 @@ +import numpy as np +from scipy.spatial import cKDTree + + +def tag_boundary_patch( + dolfinx_mesh, + patch_dataset, + patch_id, + tol=1e-6, + *, + tree: cKDTree, + facet_indices: np.ndarray, + facet_vertices: np.ndarray, +): + """Tags the facets of a dolfinx mesh that belong to a given OpenFOAM boundary patch. + + Args: + dolfinx_mesh: the dolfinx mesh + patch_dataset: the pyvista dataset for the boundary patch + patch_id: integer tag to assign to matched facets + tol: spatial tolerance for matching patch points to mesh vertices + tree: pre-built cKDTree on mesh geometry points + facet_indices: pre-computed exterior facet indices + facet_vertices: pre-computed 2D array (n_facets, n_verts_per_facet) + of vertex indices for each exterior facet + + Returns: + tuple of (matched_facet_indices, tags) as int32 arrays + """ + matched = tree.query_ball_point(patch_dataset.points, tol) + matched_idx = np.unique( + np.fromiter((i for sub in matched for i in sub), dtype=np.intp) + ) + + vertex_matched = np.zeros(len(dolfinx_mesh.geometry.x), dtype=bool) + if len(matched_idx): + vertex_matched[matched_idx] = True + + facet_mask = vertex_matched[facet_vertices].all(axis=1) + matched_facets = facet_indices[facet_mask] + + return matched_facets, np.full(len(matched_facets), patch_id, dtype=np.int32) diff --git a/src/foam2dolfinx/open_foam_reader.py b/src/foam2dolfinx/open_foam_reader.py index cea4dd5..3e57e75 100644 --- a/src/foam2dolfinx/open_foam_reader.py +++ b/src/foam2dolfinx/open_foam_reader.py @@ -5,7 +5,10 @@ import numpy as np import pyvista import ufl -from dolfinx.mesh import create_mesh +from dolfinx.mesh import create_mesh, exterior_facet_indices, meshtags +from scipy.spatial import cKDTree + +from .helpers import tag_boundary_patch __all__ = ["OpenFOAMReader", "find_closest_value"] @@ -34,6 +37,10 @@ class OpenFOAMReader: vertex indices to the original OpenFOAM point indices. Built once per subdomain on the first call to create_dolfinx_function_with_point_data and cached for subsequent calls. + subdomain_cell_offsets: dictionary mapping each subdomain name to a + (start, end) tuple of cell index ranges in the merged pyvista dataset. + Populated by _create_global_dolfinx_mesh and used to assign subdomain + IDs when building cell and interface facet meshtags. Notes: The cell type refers to the VTK cell type, a full list of cells and their @@ -54,12 +61,14 @@ class OpenFOAMReader: connectivities_dict: dict[str, np.ndarray] dolfinx_meshes_dict: dict[str, dolfinx.mesh.Mesh] vertex_maps_dict: dict[str, np.ndarray] + subdomain_cell_offsets: dict[str, tuple[int, int]] def __init__(self, filename, cell_type: int = 12): self.filename = filename self.cell_type = cell_type self.reader = pyvista.POpenFOAMReader(self.filename) + self.OF_multiblock = None self.times = self.reader.time_values self.multidomain = False self.OF_meshes_dict = {} @@ -67,6 +76,8 @@ def __init__(self, filename, cell_type: int = 12): self.connectivities_dict = {} self.dolfinx_meshes_dict = {} self.vertex_maps_dict = {} + self.subdomain_cell_offsets = {} + self._pv_cell_subdomain_ids = None @property def cell_type(self): @@ -78,88 +89,94 @@ def cell_type(self, value): raise TypeError("cell_type value should be an int") self._cell_type = value - def _read_with_pyvista(self, t: float, subdomain: str | None = "default"): - """ - Reads the OpenFOAM data in the filename provided, passes details of the - OpenFOAM mesh to OF_mesh and details of the cells to OF_cells. + def _read_with_pyvista(self, t: float, subdomain: str | None = None): + """Reads the OpenFOAM multiblock data at time t. + + Populates OF_multiblock and OF_meshes_dict for all subdomains. If + subdomain is given, validates it exists and populates OF_cells_dict for + it. In single-domain files, OF_cells_dict["default"] is always populated. + In multidomain files with no subdomain given, cell data is populated + lazily by _create_global_dolfinx_mesh. Args: t: timestamp of the data to read - subdomain: Name of the subdmain in the OpenFOAM file, from which a field is - extracted - + subdomain: if given, validate this subdomain exists and populate its + OF_cells_dict entry. Pass None when reading all subdomains + without targeting a specific one. """ - self.reader.set_active_time_value(t) # Set the time value to read data from - OF_multiblock = self.reader.read() # Read the data from the OpenFOAM file + self.reader.set_active_time_value(t) + self.OF_multiblock = self.reader.read() - # Check if the reader has a multiblock dataset block named "internalMesh" - if "internalMesh" not in OF_multiblock.keys(): + if "internalMesh" not in self.OF_multiblock.keys(): self.multidomain = True - if subdomain not in OF_multiblock.keys(): + named_subdomains = [ + k for k in self.OF_multiblock.keys() if k != "defaultRegion" + ] + if subdomain is not None and subdomain not in named_subdomains: raise ValueError( f"Subdomain {subdomain} not found in the OpenFOAM file. " - f"Available subdomains: {OF_multiblock.keys()}" + f"Available subdomains: {named_subdomains}" ) - - # Extract the internal mesh - if self.multidomain: - for cell_array_name in OF_multiblock.keys(): - self.OF_meshes_dict[cell_array_name] = OF_multiblock[cell_array_name][ - "internalMesh" - ] + for name in named_subdomains: + self.OF_meshes_dict[name] = self.OF_multiblock[name]["internalMesh"] else: - self.OF_meshes_dict[subdomain] = OF_multiblock["internalMesh"] - - # obtain dictionary of cell types in OF_mesh - OF_cell_type_dict = self.OF_meshes_dict[subdomain].cells_dict - - cell_types_in_mesh = [int(k) for k in OF_cell_type_dict.keys()] - - # Raise error if OF_mesh is mixed topology - if len(cell_types_in_mesh) > 1: - raise NotImplementedError("Cannot support mixed-topology meshes") - - self.OF_cells_dict[subdomain] = OF_cell_type_dict.get(self.cell_type) + self.multidomain = False + self.OF_meshes_dict["default"] = self.OF_multiblock["internalMesh"] + subdomain = "default" + + if subdomain is not None: + OF_cell_type_dict = self.OF_meshes_dict[subdomain].cells_dict + cell_types = [int(k) for k in OF_cell_type_dict.keys()] + if len(cell_types) > 1: + raise NotImplementedError("Cannot support mixed-topology meshes") + self.OF_cells_dict[subdomain] = OF_cell_type_dict.get(self.cell_type) + if self.OF_cells_dict[subdomain] is None: + raise ValueError( + f"No cell type {self.cell_type} found in the mesh. " + f"Found {cell_types}" + ) - # Raise error if no cells of the specified type are found in the OF_mesh - if self.OF_cells_dict[subdomain] is None: - raise ValueError( - f"No cell type {self.cell_type} found in the mesh. Found " - f"{cell_types_in_mesh}" - ) + def _get_connectivity(self, cells: np.ndarray) -> tuple[str, np.ndarray]: + """Reorders a raw VTK cell array into DOLFINx vertex ordering. - def _create_dolfinx_mesh(self, subdomain: str | None = "default"): - """Creates a dolfinx.mesh.Mesh based on the elements within the OpenFOAM mesh""" + Args: + cells: 2D array of shape (n_cells, n_verts) from PyVista cells_dict - # Define mesh element and define args conn based on the OF cell type + Returns: + (shape_name, reordered_connectivity) where shape_name is the UFL + cell name ("hexahedron" or "tetrahedron") + """ if self.cell_type == 12: shape = "hexahedron" - args_conn = np.tile( - np.array([0, 1, 3, 2, 4, 5, 7, 6]), - (len(self.OF_cells_dict[subdomain]), 1), - ) - + args_conn = np.tile(np.array([0, 1, 3, 2, 4, 5, 7, 6]), (len(cells), 1)) elif self.cell_type == 10: shape = "tetrahedron" - args_conn = np.argsort( - self.OF_cells_dict[subdomain], axis=1 - ) # Sort the cell connectivity - + args_conn = np.argsort(cells, axis=1) else: raise ValueError( f"Cell type: {self.cell_type}, not supported, please use" " either 12 (hexahedron) or 10 (tetrahedron) cells in OF mesh" ) + rows = np.arange(cells.shape[0])[:, None] + return shape, cells[rows, args_conn] + + def _build_dolfinx_mesh( + self, points: np.ndarray, connectivity: np.ndarray, shape: str + ) -> dolfinx.mesh.Mesh: + """Creates a dolfinx mesh from points and pre-reordered connectivity. + + Also sets mesh_vector_element and mesh_scalar_element on the instance. - # create the connectivity between the OpenFOAM and dolfinx meshes - # Create row indices - rows = np.arange(self.OF_cells_dict[subdomain].shape[0])[:, None] - # Reorder connectivity - self.connectivities_dict[subdomain] = self.OF_cells_dict[subdomain][ - rows, args_conn - ] + Args: + points: (n_points, 3) coordinate array + connectivity: (n_cells, n_verts) reordered connectivity from + _get_connectivity + shape: UFL cell name, e.g. "hexahedron" or "tetrahedron" - degree = 1 # Set polynomial degree + Returns: + the new dolfinx.mesh.Mesh + """ + degree = 1 cell = ufl.Cell(shape) # ufl.Cell.cellname became a property after dolfinx v0.10 cell_name = cell.cellname() if callable(cell.cellname) else cell.cellname @@ -169,16 +186,85 @@ def _create_dolfinx_mesh(self, subdomain: str | None = "default"): self.mesh_scalar_element = basix.ufl.element( "Lagrange", cell_name, degree, shape=() ) - - # Create dolfinx Mesh mesh_ufl = ufl.Mesh( basix.ufl.element("Lagrange", cell_name, degree, shape=(3,)) ) - self.dolfinx_meshes_dict[subdomain] = create_mesh( - comm=MPI.COMM_WORLD, - cells=self.connectivities_dict[subdomain], - x=self.OF_meshes_dict[subdomain].points, - e=mesh_ufl, + return create_mesh( + comm=MPI.COMM_WORLD, cells=connectivity, x=points, e=mesh_ufl + ) + + def _ensure_mesh(self) -> dolfinx.mesh.Mesh: + """Returns the active dolfinx mesh, creating it on first call. + + For multidomain cases returns the merged global mesh; for single-domain + returns the default mesh. + """ + if self.multidomain: + if "_global" not in self.dolfinx_meshes_dict: + self._create_global_dolfinx_mesh() + return self.dolfinx_meshes_dict["_global"] + else: + if "default" not in self.dolfinx_meshes_dict: + self._create_dolfinx_mesh(subdomain="default") + return self.dolfinx_meshes_dict["default"] + + def _create_dolfinx_mesh(self, subdomain: str | None = "default"): + """Creates a dolfinx.mesh.Mesh based on the elements within the OpenFOAM mesh""" + shape, connectivity = self._get_connectivity(self.OF_cells_dict[subdomain]) + self.connectivities_dict[subdomain] = connectivity + self.dolfinx_meshes_dict[subdomain] = self._build_dolfinx_mesh( + self.OF_meshes_dict[subdomain].points, connectivity, shape + ) + + def _create_global_dolfinx_mesh(self): + """Merges all subdomain pyvista meshes into a single global dolfinx mesh. + + Merges all pyvista meshes with clean() to deduplicate shared interface + points, embeds subdomain IDs as cell data so they survive any reordering, + and creates a single dolfinx mesh stored under the key "_global". + """ + subdomain_names = list(self.OF_meshes_dict.keys()) + + # Populate OF_cells_dict for any subdomains not yet read + for name in subdomain_names: + if name not in self.OF_cells_dict: + OF_cell_type_dict = self.OF_meshes_dict[name].cells_dict + cell_types = [int(k) for k in OF_cell_type_dict.keys()] + if len(cell_types) > 1: + raise NotImplementedError("Cannot support mixed-topology meshes") + cells = OF_cell_type_dict.get(self.cell_type) + if cells is None: + raise ValueError( + f"No cell type {self.cell_type} found in subdomain {name}. " + f"Found {cell_types}" + ) + self.OF_cells_dict[name] = cells + + # Record cumulative cell offsets (used for print summaries only) + cumulative = 0 + for name in subdomain_names: + count = len(self.OF_cells_dict[name]) + self.subdomain_cell_offsets[name] = (cumulative, cumulative + count) + cumulative += count + + # Embed subdomain IDs as cell data before merging so the tag survives + # merge+clean regardless of any cell reordering the operation applies + tagged_meshes = [] + for j, name in enumerate(subdomain_names): + pv_mesh = self.OF_meshes_dict[name].copy() + pv_mesh.cell_data["subdomain_id"] = np.full( + pv_mesh.n_cells, j + 1, dtype=np.int32 + ) + tagged_meshes.append(pv_mesh) + + merged = pyvista.merge(tagged_meshes).clean() + self._pv_cell_subdomain_ids = merged.cell_data["subdomain_id"] + + OF_cells = merged.cells_dict.get(self.cell_type) + shape, connectivity = self._get_connectivity(OF_cells) + self.connectivities_dict["_global"] = connectivity + self.dolfinx_meshes_dict["_global"] = self._build_dolfinx_mesh( + merged.points, connectivity, shape ) def _get_mesh( @@ -291,6 +377,171 @@ def create_dolfinx_function_with_point_data( return u + def create_facet_meshtags(self, t: float | None = None) -> dolfinx.mesh.MeshTags: + """Creates a dolfinx.mesh.MeshTags for all tagged facets of the mesh. + + For single-domain meshes, tags external boundary patches (IDs starting at 1). + For multidomain meshes, also tags interface facets between subdomains with + sequential IDs continuing from the last boundary patch ID. + + Args: + t: timestamp of the data to read, defaults to the first available time. + + Returns: + the dolfinx MeshTags for the facets of the mesh + """ + t = t if t is not None else next(tv for tv in self.times if tv != 0) + self._read_with_pyvista(t=t) + + mesh = self._ensure_mesh() + + # Collect boundary patches — in multidomain files each subdomain block + # holds its own "boundary" child; single-domain has one top-level "boundary" + if self.multidomain: + patches = [] + for sd_name in self.OF_meshes_dict.keys(): + sd_block = self.OF_multiblock[sd_name] + if "boundary" in sd_block.keys(): + boundary = sd_block["boundary"] + for patch_name in boundary.keys(): + patches.append((patch_name, boundary[patch_name])) + else: + boundary = self.OF_multiblock["boundary"] + patches = [(name, boundary[name]) for name in boundary.keys()] + + # build shared data once across all patches + fdim = mesh.topology.dim - 1 + mesh.topology.create_connectivity(fdim, 0) + mesh.topology.create_connectivity(0, fdim) + mesh.topology.create_connectivity(fdim, mesh.topology.dim) + facet_indices = exterior_facet_indices(mesh.topology) + c_to_v = mesh.topology.connectivity(fdim, 0) + facet_vertices = np.vstack([c_to_v.links(f) for f in facet_indices]) + tree = cKDTree(mesh.geometry.x) + + all_facets = [] + all_tags = [] + patch_summary = {} + + for i, (patch_name, patch_dataset) in enumerate(patches): + facets, tags = tag_boundary_patch( + mesh, + patch_dataset, + i + 1, + tree=tree, + facet_indices=facet_indices, + facet_vertices=facet_vertices, + ) + all_facets.append(facets) + all_tags.append(tags) + patch_summary[patch_name] = {"id": i + 1, "n_facets": len(facets)} + + print("Boundary patch summary:") + for patch_name, info in patch_summary.items(): + print(f" {patch_name}: id={info['id']}, n_facets={info['n_facets']}") + + next_tag = len(patches) + 1 + + if self.multidomain: + # Map subdomain IDs to dolfinx cell ordering via original_cell_index. + # _pv_cell_subdomain_ids is indexed by the merged pyvista cell index, + # which original_cell_index maps back to for each dolfinx cell. + num_cells = mesh.topology.index_map(mesh.topology.dim).size_local + dolfinx_cell_tags = self._pv_cell_subdomain_ids[ + mesh.topology.original_cell_index[:num_cells] + ] + + # Find internal facets where adjacent cells belong to different subdomains + f_to_c = mesh.topology.connectivity(fdim, mesh.topology.dim) + num_facets = mesh.topology.index_map(fdim).size_local + n_adj = np.array([len(f_to_c.links(f)) for f in range(num_facets)]) + internal_facet_ids = np.where(n_adj == 2)[0].astype(np.int32) + cell_pairs = np.array([f_to_c.links(f) for f in internal_facet_ids]) + tags_a = dolfinx_cell_tags[cell_pairs[:, 0]] + tags_b = dolfinx_cell_tags[cell_pairs[:, 1]] + is_interface = tags_a != tags_b + interface_facet_ids = internal_facet_ids[is_interface] + tags_a_iface = tags_a[is_interface] + tags_b_iface = tags_b[is_interface] + + # Assign one sequential tag per unique subdomain pair + pair_to_tag: dict[tuple[int, int], int] = {} + interface_tags = np.empty(len(interface_facet_ids), dtype=np.int32) + for j in range(len(interface_facet_ids)): + pair_key = ( + min(tags_a_iface[j], tags_b_iface[j]), + max(tags_a_iface[j], tags_b_iface[j]), + ) + if pair_key not in pair_to_tag: + pair_to_tag[pair_key] = next_tag + next_tag += 1 + interface_tags[j] = pair_to_tag[pair_key] + + all_facets.append(interface_facet_ids) + all_tags.append(interface_tags) + + print("Interface summary:") + sd_names = list(self.subdomain_cell_offsets.keys()) + for (id_a, id_b), tag in pair_to_tag.items(): + n_iface = int(np.sum(interface_tags == tag)) + print( + f" {sd_names[id_a - 1]}-{sd_names[id_b - 1]}: " + f"id={tag}, n_facets={n_iface}" + ) + + all_facets_arr = np.concatenate(all_facets) + all_tags_arr = np.concatenate(all_tags) + sort_idx = np.argsort(all_facets_arr) + return meshtags( + mesh, + fdim, + all_facets_arr[sort_idx].astype(np.int32), + all_tags_arr[sort_idx].astype(np.int32), + ) + + def create_cell_meshtags(self, t: float | None = None) -> dolfinx.mesh.MeshTags: + """Creates a dolfinx.mesh.MeshTags for the cells of the mesh. + + For single-domain meshes, all cells are tagged with ID 1. + For multidomain meshes, cells are tagged with their subdomain ID (1-indexed + in the order subdomains appear in the OpenFOAM file). + + Args: + t: timestamp of the data to read, defaults to the first available time. + + Returns: + the dolfinx MeshTags for the cells of the mesh + """ + t = t if t is not None else next(tv for tv in self.times if tv != 0) + self._read_with_pyvista(t=t) + + if self.multidomain: + if "_global" not in self.dolfinx_meshes_dict: + self._create_global_dolfinx_mesh() + mesh = self.dolfinx_meshes_dict["_global"] + + num_cells = mesh.topology.index_map(mesh.topology.dim).size_local + cell_tag_array = self._pv_cell_subdomain_ids[ + mesh.topology.original_cell_index[:num_cells] + ] + + print("Cell zone summary:") + for j, sd_name in enumerate(self.subdomain_cell_offsets.keys()): + n_cells = int(np.sum(cell_tag_array == j + 1)) + print(f" {sd_name}: id={j + 1}, n_cells={n_cells}") + else: + if "default" not in self.dolfinx_meshes_dict: + self._create_dolfinx_mesh(subdomain="default") + mesh = self.dolfinx_meshes_dict["default"] + num_cells = mesh.topology.index_map(mesh.topology.dim).size_local + cell_tag_array = np.ones(num_cells, dtype=np.int32) + + print("Cell zone summary:") + print(f" default: id=1, n_cells={num_cells}") + + cell_indices = np.arange(num_cells, dtype=np.int32) + return meshtags(mesh, mesh.topology.dim, cell_indices, cell_tag_array) + def find_closest_value(values: list[float], target: float) -> float: """ diff --git a/test/conftest.py b/test/conftest.py new file mode 100644 index 0000000..707b296 --- /dev/null +++ b/test/conftest.py @@ -0,0 +1,16 @@ +import zipfile +from pathlib import Path + +import pytest + +from foam2dolfinx import OpenFOAMReader + + +@pytest.fixture +def two_regions_reader(tmp_path): + zip_path = Path("test/data/test_2Regions.zip") + extract_path = tmp_path / "test_2Regions" + with zipfile.ZipFile(zip_path, "r") as zf: + zf.extractall(extract_path) + foam_file = extract_path / "test_2Regions/pv.foam" + return OpenFOAMReader(filename=str(foam_file), cell_type=12) diff --git a/test/test_create_function.py b/test/test_create_function.py index 736a06d..3255a72 100644 --- a/test/test_create_function.py +++ b/test/test_create_function.py @@ -24,7 +24,8 @@ def test_not_finding_function_cell_data(tmpdir): with pytest.raises( ValueError, match=re.escape( - "Function name: coucou not found in the subdomain: solid, in the OpenFOAM file. " + "Function name: coucou not found in the subdomain: solid, " + "in the OpenFOAM file. " "Available functions in subdomain: solid : ['T']" ), ): @@ -50,7 +51,8 @@ def test_not_finding_function_point_data(tmpdir): with pytest.raises( ValueError, match=re.escape( - "Function name: coucou not found in the subdomain: solid, in the OpenFOAM file. " + "Function name: coucou not found in the subdomain: solid, " + "in the OpenFOAM file. " "Available functions in subdomain: solid : ['T']" ), ): diff --git a/test/test_get_connectivity.py b/test/test_get_connectivity.py new file mode 100644 index 0000000..3e0d758 --- /dev/null +++ b/test/test_get_connectivity.py @@ -0,0 +1,82 @@ +import dolfinx +import numpy as np +import pytest +from pyvista import examples + +from foam2dolfinx import OpenFOAMReader + + +@pytest.fixture +def hex_reader(): + return OpenFOAMReader(filename=examples.download_cavity(load=False), cell_type=12) + + +@pytest.fixture +def tet_reader(): + return OpenFOAMReader(filename=examples.download_cavity(load=False), cell_type=10) + + +@pytest.mark.parametrize( + "cell_type, n_verts, expected_shape", + [(12, 8, "hexahedron"), (10, 4, "tetrahedron")], +) +def test_get_connectivity_returns_correct_shape_name( + cell_type, n_verts, expected_shape +): + reader = OpenFOAMReader( + filename=examples.download_cavity(load=False), cell_type=cell_type + ) + shape_name, _ = reader._get_connectivity(np.zeros((3, n_verts), dtype=int)) + assert shape_name == expected_shape + + +@pytest.mark.parametrize("cell_type, n_verts", [(12, 8), (10, 4)]) +def test_get_connectivity_output_shape_matches_input(cell_type, n_verts): + reader = OpenFOAMReader( + filename=examples.download_cavity(load=False), cell_type=cell_type + ) + cells = np.arange(5 * n_verts).reshape(5, n_verts) + _, connectivity = reader._get_connectivity(cells) + assert connectivity.shape == cells.shape + + +def test_get_connectivity_raises_for_unsupported_cell_type(): + reader = OpenFOAMReader(filename=examples.download_cavity(load=False), cell_type=12) + reader._cell_type = 5 + with pytest.raises(ValueError, match="not supported"): + reader._get_connectivity(np.zeros((3, 4), dtype=int)) + + +def test_get_connectivity_hexahedron_applies_vtk_to_dolfinx_permutation(hex_reader): + cells = np.arange(8).reshape(1, 8) + _, connectivity = hex_reader._get_connectivity(cells) + np.testing.assert_array_equal(connectivity[0], [0, 1, 3, 2, 4, 5, 7, 6]) + + +def test_get_connectivity_tetrahedron_sorts_vertices(tet_reader): + cells = np.array([[40, 10, 30, 20]]) + _, connectivity = tet_reader._get_connectivity(cells) + np.testing.assert_array_equal(connectivity[0], [10, 20, 30, 40]) + + +def test_build_dolfinx_mesh_returns_dolfinx_mesh(hex_reader): + hex_reader._read_with_pyvista(t=0) + shape, connectivity = hex_reader._get_connectivity( + hex_reader.OF_cells_dict["default"] + ) + result = hex_reader._build_dolfinx_mesh( + hex_reader.OF_meshes_dict["default"].points, connectivity, shape + ) + assert isinstance(result, dolfinx.mesh.Mesh) + + +@pytest.mark.parametrize("attr", ["mesh_vector_element", "mesh_scalar_element"]) +def test_build_dolfinx_mesh_sets_element_attributes(hex_reader, attr): + hex_reader._read_with_pyvista(t=0) + shape, connectivity = hex_reader._get_connectivity( + hex_reader.OF_cells_dict["default"] + ) + hex_reader._build_dolfinx_mesh( + hex_reader.OF_meshes_dict["default"].points, connectivity, shape + ) + assert getattr(hex_reader, attr) is not None diff --git a/test/test_meshtags.py b/test/test_meshtags.py new file mode 100644 index 0000000..95455d7 --- /dev/null +++ b/test/test_meshtags.py @@ -0,0 +1,83 @@ +import numpy as np +import pytest +from pyvista import examples + +from foam2dolfinx import OpenFOAMReader + + +# --- create_cell_meshtags: multidomain --- + + +def test_create_cell_meshtags_dim_equals_mesh_dim(two_regions_reader): + ct = two_regions_reader.create_cell_meshtags() + mesh = two_regions_reader.dolfinx_meshes_dict["_global"] + assert ct.dim == mesh.topology.dim + + +def test_create_cell_meshtags_multidomain_fluid_count(two_regions_reader): + ct = two_regions_reader.create_cell_meshtags() + assert np.sum(ct.values == 1) == 2100 + + +def test_create_cell_meshtags_multidomain_solid_count(two_regions_reader): + ct = two_regions_reader.create_cell_meshtags() + assert np.sum(ct.values == 2) == 945 + + +def test_create_cell_meshtags_multidomain_no_untagged_cells(two_regions_reader): + ct = two_regions_reader.create_cell_meshtags() + assert np.all(ct.values > 0) + + +# --- create_cell_meshtags: single domain --- + + +def test_create_cell_meshtags_single_domain_dim_equals_mesh_dim(): + reader = OpenFOAMReader(filename=examples.download_cavity(load=False), cell_type=12) + ct = reader.create_cell_meshtags(t=0) + mesh = reader.dolfinx_meshes_dict["default"] + assert ct.dim == mesh.topology.dim + + +def test_create_cell_meshtags_single_domain_all_values_are_one(): + reader = OpenFOAMReader(filename=examples.download_cavity(load=False), cell_type=12) + ct = reader.create_cell_meshtags(t=0) + assert np.all(ct.values == 1) + + +# --- create_facet_meshtags: multidomain --- + + +def test_create_facet_meshtags_dim_equals_fdim(two_regions_reader): + ft = two_regions_reader.create_facet_meshtags() + mesh = two_regions_reader.dolfinx_meshes_dict["_global"] + assert ft.dim == mesh.topology.dim - 1 + + +@pytest.mark.parametrize("tag, expected_count", [(1, 15), (2, 15), (3, 70)]) +def test_create_facet_meshtags_boundary_patch_counts( + two_regions_reader, tag, expected_count +): + ft = two_regions_reader.create_facet_meshtags() + assert np.sum(ft.values == tag) == expected_count + + +def test_create_facet_meshtags_interface_count(two_regions_reader): + ft = two_regions_reader.create_facet_meshtags() + # interface facets receive the highest tag (assigned after all boundary patches) + assert np.sum(ft.values == ft.values.max()) == 60 + + +def test_create_facet_meshtags_all_tag_values_positive(two_regions_reader): + ft = two_regions_reader.create_facet_meshtags() + assert np.all(ft.values > 0) + + +# --- create_facet_meshtags: single domain --- + + +def test_create_facet_meshtags_single_domain_dim_equals_fdim(): + reader = OpenFOAMReader(filename=examples.download_cavity(load=False), cell_type=12) + ft = reader.create_facet_meshtags() + mesh = reader.dolfinx_meshes_dict["default"] + assert ft.dim == mesh.topology.dim - 1 diff --git a/test/test_read_with_pyvista.py b/test/test_read_with_pyvista.py index eb72f96..69944e4 100644 --- a/test/test_read_with_pyvista.py +++ b/test/test_read_with_pyvista.py @@ -28,7 +28,7 @@ def test_error_rasied_when_cells_wanted_are_not_in_file_provided(): my_reader._read_with_pyvista(t=0) -def test_error_rasied_when_subdomain_is_not_given_in_multidomain_case(tmpdir): +def test_error_rasied_when_wrong_subdomain_given_in_multidomain_case(tmpdir): zip_path = Path("test/data/test_2Regions.zip") extract_path = Path(tmpdir) / "test_2Regions" @@ -45,11 +45,53 @@ def test_error_rasied_when_subdomain_is_not_given_in_multidomain_case(tmpdir): with pytest.raises( ValueError, match=( - r"Subdomain None not found in the OpenFOAM file\. " - r"Available subdomains: \['defaultRegion', 'fluid', 'solid']" + r"Subdomain coucou not found in the OpenFOAM file\. " + r"Available subdomains: \['fluid', 'solid']" ), ): - my_of_reader._read_with_pyvista(t=20.0, subdomain=None) + my_of_reader._read_with_pyvista(t=20.0, subdomain="coucou") + + +def test_read_with_pyvista_no_subdomain_single_domain_populates_cells(): + reader = OpenFOAMReader(filename=examples.download_cavity(load=False), cell_type=12) + reader._read_with_pyvista(t=0) + assert "default" in reader.OF_cells_dict + assert reader.OF_cells_dict["default"] is not None + + +def test_read_with_pyvista_no_subdomain_single_domain_sets_multidomain_false(): + reader = OpenFOAMReader(filename=examples.download_cavity(load=False), cell_type=12) + reader._read_with_pyvista(t=0) + assert reader.multidomain is False + + +def test_read_with_pyvista_no_subdomain_multidomain_populates_all_meshes( + two_regions_reader, +): + two_regions_reader._read_with_pyvista(t=20.0) + assert "fluid" in two_regions_reader.OF_meshes_dict + assert "solid" in two_regions_reader.OF_meshes_dict + + +def test_read_with_pyvista_no_subdomain_multidomain_excludes_defaultRegion( + two_regions_reader, +): + two_regions_reader._read_with_pyvista(t=20.0) + assert "defaultRegion" not in two_regions_reader.OF_meshes_dict + + +def test_read_with_pyvista_no_subdomain_multidomain_does_not_populate_cells( + two_regions_reader, +): + two_regions_reader._read_with_pyvista(t=20.0) + assert len(two_regions_reader.OF_cells_dict) == 0 + + +def test_read_with_pyvista_no_subdomain_multidomain_sets_multidomain_true( + two_regions_reader, +): + two_regions_reader._read_with_pyvista(t=20.0) + assert two_regions_reader.multidomain is True @pytest.mark.parametrize("subdomain", ["fluid", "solid"]) diff --git a/test/test_tag_boundary_patch.py b/test/test_tag_boundary_patch.py new file mode 100644 index 0000000..bf1362a --- /dev/null +++ b/test/test_tag_boundary_patch.py @@ -0,0 +1,97 @@ +import dolfinx.mesh +import numpy as np +import pyvista +import pytest +from mpi4py import MPI +from scipy.spatial import cKDTree +from dolfinx.mesh import exterior_facet_indices + +from foam2dolfinx.helpers import tag_boundary_patch + + +@pytest.fixture +def unit_cube_topology(): + mesh = dolfinx.mesh.create_box( + MPI.COMM_WORLD, + [[0, 0, 0], [1, 1, 1]], + [2, 2, 2], + dolfinx.mesh.CellType.hexahedron, + ) + fdim = mesh.topology.dim - 1 + mesh.topology.create_connectivity(fdim, 0) + mesh.topology.create_connectivity(0, fdim) + mesh.topology.create_connectivity(fdim, mesh.topology.dim) + facet_indices = exterior_facet_indices(mesh.topology) + c_to_v = mesh.topology.connectivity(fdim, 0) + facet_vertices = np.vstack([c_to_v.links(f) for f in facet_indices]) + tree = cKDTree(mesh.geometry.x) + return mesh, tree, facet_indices, facet_vertices + + +def test_tag_boundary_patch_returns_empty_when_no_points_match(unit_cube_topology): + mesh, tree, facet_indices, facet_vertices = unit_cube_topology + patch = pyvista.PolyData(np.array([[100.0, 100.0, 100.0]])) + matched, tags = tag_boundary_patch( + mesh, + patch, + 1, + tree=tree, + facet_indices=facet_indices, + facet_vertices=facet_vertices, + ) + assert len(matched) == 0 + assert len(tags) == 0 + + +def test_tag_boundary_patch_tags_equal_patch_id(unit_cube_topology): + mesh, tree, facet_indices, facet_vertices = unit_cube_topology + x0_points = mesh.geometry.x[mesh.geometry.x[:, 0] < 1e-10] + patch = pyvista.PolyData(x0_points) + _, tags = tag_boundary_patch( + mesh, + patch, + 42, + tree=tree, + facet_indices=facet_indices, + facet_vertices=facet_vertices, + ) + assert np.all(tags == 42) + + +def test_tag_boundary_patch_matched_facets_are_subset_of_exterior(unit_cube_topology): + mesh, tree, facet_indices, facet_vertices = unit_cube_topology + x0_points = mesh.geometry.x[mesh.geometry.x[:, 0] < 1e-10] + patch = pyvista.PolyData(x0_points) + matched, _ = tag_boundary_patch( + mesh, + patch, + 1, + tree=tree, + facet_indices=facet_indices, + facet_vertices=facet_vertices, + ) + assert np.all(np.isin(matched, facet_indices)) + + +@pytest.mark.parametrize( + "patch_points", + [ + np.array([[100.0, 100.0, 100.0]]), # no match + None, # match (x=0 face, filled below) + ], +) +def test_tag_boundary_patch_output_arrays_are_int32(unit_cube_topology, patch_points): + mesh, tree, facet_indices, facet_vertices = unit_cube_topology + if patch_points is None: + patch_points = mesh.geometry.x[mesh.geometry.x[:, 0] < 1e-10] + patch = pyvista.PolyData(patch_points) + matched, tags = tag_boundary_patch( + mesh, + patch, + 1, + tree=tree, + facet_indices=facet_indices, + facet_vertices=facet_vertices, + ) + assert matched.dtype == np.int32 + assert tags.dtype == np.int32