From ac8426742abc0a1c8371465e46821526dc7201cc Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 23 Apr 2026 22:21:51 -0400 Subject: [PATCH 01/15] need scipy --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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", From a2ef8fc043545d134730c93871d141e567bfe0dd Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 23 Apr 2026 22:22:08 -0400 Subject: [PATCH 02/15] tag boundary patch --- src/foam2dolfinx/helpers.py | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) create mode 100644 src/foam2dolfinx/helpers.py diff --git a/src/foam2dolfinx/helpers.py b/src/foam2dolfinx/helpers.py new file mode 100644 index 0000000..d3aef71 --- /dev/null +++ b/src/foam2dolfinx/helpers.py @@ -0,0 +1,28 @@ +import numpy as np +from dolfinx.mesh import exterior_facet_indices +from scipy.spatial import cKDTree + + +def tag_boundary_patch(dolfinx_mesh, patch_dataset, patch_id, tol=1e-6): + fdim = dolfinx_mesh.topology.dim - 1 + dolfinx_mesh.topology.create_connectivity(fdim, 0) + dolfinx_mesh.topology.create_connectivity(0, fdim) + dolfinx_mesh.topology.create_connectivity(fdim, dolfinx_mesh.topology.dim) + + facet_indices = exterior_facet_indices(dolfinx_mesh.topology) + x = dolfinx_mesh.geometry.x + patch_points = patch_dataset.points + tree = cKDTree(x) + matched_vertex_indices = tree.query_ball_point(patch_points, tol) + matched_vertex_indices = list(set(i for sub in matched_vertex_indices for i in sub)) + + matched_facets = [] + for facet in facet_indices: + vertices = dolfinx_mesh.topology.connectivity(fdim, 0).links(facet) + if all(v in matched_vertex_indices for v in vertices): + matched_facets.append(facet) + + # print(f"Tagging {len(matched_facets)} facets for patch ID {patch_id}") + return np.array(matched_facets, dtype=np.int32), np.full( + len(matched_facets), patch_id, dtype=np.int32 + ) From dd4441f2277a7ef4138914f186824bc853d1cfb1 Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 23 Apr 2026 22:22:20 -0400 Subject: [PATCH 03/15] create facet meshtags Co-authored-by: Copilot --- src/foam2dolfinx/open_foam_reader.py | 61 +++++++++++++++++++++++----- 1 file changed, 51 insertions(+), 10 deletions(-) diff --git a/src/foam2dolfinx/open_foam_reader.py b/src/foam2dolfinx/open_foam_reader.py index cea4dd5..8f682fb 100644 --- a/src/foam2dolfinx/open_foam_reader.py +++ b/src/foam2dolfinx/open_foam_reader.py @@ -5,7 +5,9 @@ import numpy as np import pyvista import ufl -from dolfinx.mesh import create_mesh +from dolfinx.mesh import create_mesh, meshtags + +from .helpers import tag_boundary_patch __all__ = ["OpenFOAMReader", "find_closest_value"] @@ -60,6 +62,7 @@ def __init__(self, filename, cell_type: int = 12): 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 = {} @@ -90,25 +93,25 @@ def _read_with_pyvista(self, t: float, subdomain: str | None = "default"): """ 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.OF_multiblock = self.reader.read() # Read the data from the OpenFOAM file # 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(): + if subdomain not in self.OF_multiblock.keys(): raise ValueError( f"Subdomain {subdomain} not found in the OpenFOAM file. " - f"Available subdomains: {OF_multiblock.keys()}" + f"Available subdomains: {self.OF_multiblock.keys()}" ) # 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 cell_array_name in self.OF_multiblock.keys(): + self.OF_meshes_dict[cell_array_name] = self.OF_multiblock[ + cell_array_name + ]["internalMesh"] else: - self.OF_meshes_dict[subdomain] = OF_multiblock["internalMesh"] + self.OF_meshes_dict[subdomain] = self.OF_multiblock["internalMesh"] # obtain dictionary of cell types in OF_mesh OF_cell_type_dict = self.OF_meshes_dict[subdomain].cells_dict @@ -291,6 +294,44 @@ def create_dolfinx_function_with_point_data( return u + def create_facet_meshtags( + self, t: float = 0, name: str = "U", subdomain: str | None = "default" + ) -> dolfinx.mesh.MeshTags: + """Creates a dolfinx.mesh.MeshTags for the facets of the mesh based on the + boundary patches in the OpenFOAM file. + + Args: + t: timestamp of the data to read, default to 0. + name: Name of the field in the OpenFOAM file, defaults to "U" for velocity + subdomain: Name of the subdmain in the OpenFOAM file, from which a field is + extracted + + Returns: + the dolfinx MeshTags for the facets of the mesh + """ + mesh = self._get_mesh(t, name, subdomain) + + boundary = self.OF_multiblock["boundary"] + + all_facets = np.array([], dtype=np.int32) + all_tags = np.array([], dtype=np.int32) + + for i, name in enumerate(boundary.keys()): + facets, tags = tag_boundary_patch(mesh, boundary[name], i + 1) + all_facets = np.concatenate([all_facets, facets]) + all_tags = np.concatenate([all_tags, tags]) + + print(f"Tagging {len(facets)} facets for patch {name} with ID {i + 1}") + + facet_meshtags = meshtags( + mesh, + mesh.topology.dim - 1, + all_facets, + all_tags, + ) + + return facet_meshtags + def find_closest_value(values: list[float], target: float) -> float: """ From 92b45860be2bbd31e2f6267e8653a6b6e28f8a4f Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 23 Apr 2026 22:31:35 -0400 Subject: [PATCH 04/15] optimisations by claude --- src/foam2dolfinx/helpers.py | 73 ++++++++++++++++++++-------- src/foam2dolfinx/open_foam_reader.py | 50 +++++++++++++------ 2 files changed, 87 insertions(+), 36 deletions(-) diff --git a/src/foam2dolfinx/helpers.py b/src/foam2dolfinx/helpers.py index d3aef71..5b0963b 100644 --- a/src/foam2dolfinx/helpers.py +++ b/src/foam2dolfinx/helpers.py @@ -3,26 +3,57 @@ from scipy.spatial import cKDTree -def tag_boundary_patch(dolfinx_mesh, patch_dataset, patch_id, tol=1e-6): +def tag_boundary_patch( + dolfinx_mesh, + patch_dataset, + patch_id, + tol=1e-6, + *, + tree=None, + facet_indices=None, + facet_vertices=None, +): + """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: optional pre-built cKDTree on mesh geometry points — pass this when + calling for multiple patches to avoid rebuilding the tree each time + facet_indices: optional pre-computed exterior facet indices — pass together + with facet_vertices when calling for multiple patches + facet_vertices: optional 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 + """ fdim = dolfinx_mesh.topology.dim - 1 - dolfinx_mesh.topology.create_connectivity(fdim, 0) - dolfinx_mesh.topology.create_connectivity(0, fdim) - dolfinx_mesh.topology.create_connectivity(fdim, dolfinx_mesh.topology.dim) - - facet_indices = exterior_facet_indices(dolfinx_mesh.topology) - x = dolfinx_mesh.geometry.x - patch_points = patch_dataset.points - tree = cKDTree(x) - matched_vertex_indices = tree.query_ball_point(patch_points, tol) - matched_vertex_indices = list(set(i for sub in matched_vertex_indices for i in sub)) - - matched_facets = [] - for facet in facet_indices: - vertices = dolfinx_mesh.topology.connectivity(fdim, 0).links(facet) - if all(v in matched_vertex_indices for v in vertices): - matched_facets.append(facet) - - # print(f"Tagging {len(matched_facets)} facets for patch ID {patch_id}") - return np.array(matched_facets, dtype=np.int32), np.full( - len(matched_facets), patch_id, dtype=np.int32 + + if facet_indices is None or facet_vertices is None: + dolfinx_mesh.topology.create_connectivity(fdim, 0) + dolfinx_mesh.topology.create_connectivity(0, fdim) + dolfinx_mesh.topology.create_connectivity(fdim, dolfinx_mesh.topology.dim) + facet_indices = exterior_facet_indices(dolfinx_mesh.topology) + c_to_v = dolfinx_mesh.topology.connectivity(fdim, 0) + facet_vertices = np.vstack([c_to_v.links(f) for f in facet_indices]) + + if tree is None: + tree = cKDTree(dolfinx_mesh.geometry.x) + + 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) ) + + # boolean mask over all vertices, then check all vertices of each facet at once + 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 8f682fb..ee66c04 100644 --- a/src/foam2dolfinx/open_foam_reader.py +++ b/src/foam2dolfinx/open_foam_reader.py @@ -5,7 +5,8 @@ import numpy as np import pyvista import ufl -from dolfinx.mesh import create_mesh, meshtags +from dolfinx.mesh import create_mesh, exterior_facet_indices, meshtags +from scipy.spatial import cKDTree from .helpers import tag_boundary_patch @@ -313,25 +314,44 @@ def create_facet_meshtags( boundary = self.OF_multiblock["boundary"] - all_facets = np.array([], dtype=np.int32) - all_tags = np.array([], dtype=np.int32) - - for i, name in enumerate(boundary.keys()): - facets, tags = tag_boundary_patch(mesh, boundary[name], i + 1) - all_facets = np.concatenate([all_facets, facets]) - all_tags = np.concatenate([all_tags, tags]) + # 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 in enumerate(boundary.keys()): + facets, tags = tag_boundary_patch( + mesh, + boundary[patch_name], + 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(f"Tagging {len(facets)} facets for patch {name} with ID {i + 1}") + print("Boundary patch summary:") + for patch_name, info in patch_summary.items(): + print(f" {patch_name}: id={info['id']}, n_facets={info['n_facets']}") - facet_meshtags = meshtags( + return meshtags( mesh, - mesh.topology.dim - 1, - all_facets, - all_tags, + fdim, + np.concatenate(all_facets), + np.concatenate(all_tags), ) - return facet_meshtags - def find_closest_value(values: list[float], target: float) -> float: """ From 79348b7607f1f42bc2c58435d60590f7f96f2273 Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 23 Apr 2026 23:02:12 -0400 Subject: [PATCH 05/15] create cell meshtags, but all meshtags only work in single domain Co-authored-by: Copilot --- src/foam2dolfinx/open_foam_reader.py | 34 ++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/src/foam2dolfinx/open_foam_reader.py b/src/foam2dolfinx/open_foam_reader.py index ee66c04..cac196a 100644 --- a/src/foam2dolfinx/open_foam_reader.py +++ b/src/foam2dolfinx/open_foam_reader.py @@ -310,6 +310,11 @@ def create_facet_meshtags( Returns: the dolfinx MeshTags for the facets of the mesh """ + if self.multidomain is True: + raise NotImplementedError( + "Facet meshtags cannot be created for multidomain meshes, as the " + "boundary patches are not currently supported for multidomain meshes." + ) mesh = self._get_mesh(t, name, subdomain) boundary = self.OF_multiblock["boundary"] @@ -352,6 +357,35 @@ def create_facet_meshtags( np.concatenate(all_tags), ) + def create_cell_meshtags( + self, t: float = 0, name: str = "U", subdomain: str | None = "default" + ): + """Creates a dolfinx.mesh.MeshTags for the cells of the mesh based on the + cell data in the OpenFOAM file. + + Args: + t: timestamp of the data to read, default to 0. + name: Name of the field in the OpenFOAM file, defaults to "U" for velocity + subdomain: Name of the subdmain in the OpenFOAM file, from which a field is + extracted + + Returns: + the dolfinx MeshTags for the cells of the mesh + """ + if self.multidomain is True: + raise NotImplementedError( + "Facet meshtags cannot be created for multidomain meshes, as the " + "boundary patches are not currently supported for multidomain meshes." + ) + + mesh = self._get_mesh(t, name, subdomain) + + num_cells = mesh.topology.index_map(mesh.topology.dim).size_local + mesh_cell_indices = np.arange(num_cells, dtype=np.int32) + tags_volumes = np.full(num_cells, 1, dtype=np.int32) + + return meshtags(mesh, mesh.topology.dim, mesh_cell_indices, tags_volumes) + def find_closest_value(values: list[float], target: float) -> float: """ From 715e19a737bcc8556288058699e3ac3a8dee4b1a Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 23 Apr 2026 23:29:16 -0400 Subject: [PATCH 06/15] global options for making cell tags and facet tags in multidomain --- src/foam2dolfinx/open_foam_reader.py | 276 +++++++++++++++++++++++---- 1 file changed, 235 insertions(+), 41 deletions(-) diff --git a/src/foam2dolfinx/open_foam_reader.py b/src/foam2dolfinx/open_foam_reader.py index cac196a..f4c1bf4 100644 --- a/src/foam2dolfinx/open_foam_reader.py +++ b/src/foam2dolfinx/open_foam_reader.py @@ -37,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 @@ -57,6 +61,7 @@ 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 @@ -71,6 +76,7 @@ def __init__(self, filename, cell_type: int = 12): self.connectivities_dict = {} self.dolfinx_meshes_dict = {} self.vertex_maps_dict = {} + self.subdomain_cell_offsets = {} @property def cell_type(self): @@ -132,6 +138,38 @@ def _read_with_pyvista(self, t: float, subdomain: str | None = "default"): f"{cell_types_in_mesh}" ) + def _read_with_pyvista_all(self, t: float): + """Reads the OpenFOAM multiblock data at time t for all subdomains. + + Populates OF_multiblock and OF_meshes_dict for every subdomain without + requiring a specific subdomain to be named. For single-domain files, + also populates OF_cells_dict["default"]. For multidomain files, cell + data for each subdomain is read lazily by _create_global_dolfinx_mesh. + + Args: + t: timestamp of the data to read + """ + self.reader.set_active_time_value(t) + self.OF_multiblock = self.reader.read() + + if "internalMesh" not in self.OF_multiblock.keys(): + self.multidomain = True + for name in self.OF_multiblock.keys(): + self.OF_meshes_dict[name] = self.OF_multiblock[name]["internalMesh"] + else: + self.multidomain = False + self.OF_meshes_dict["default"] = self.OF_multiblock["internalMesh"] + OF_cell_type_dict = self.OF_meshes_dict["default"].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["default"] = OF_cell_type_dict.get(self.cell_type) + if self.OF_cells_dict["default"] is None: + raise ValueError( + f"No cell type {self.cell_type} found in the mesh. " + f"Found {cell_types}" + ) + def _create_dolfinx_mesh(self, subdomain: str | None = "default"): """Creates a dolfinx.mesh.Mesh based on the elements within the OpenFOAM mesh""" @@ -185,6 +223,83 @@ def _create_dolfinx_mesh(self, subdomain: str | None = "default"): e=mesh_ufl, ) + def _create_global_dolfinx_mesh(self): + """Merges all subdomain pyvista meshes into a single global dolfinx mesh. + + Reads cell data for any subdomains not yet in OF_cells_dict, records + cumulative cell offsets in subdomain_cell_offsets (used later for cell + and interface facet meshtags), merges all pyvista meshes with clean() + to deduplicate shared interface points, and creates a single dolfinx + mesh stored under the key "_global" in dolfinx_meshes_dict. + """ + 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 before merging (cell order is preserved) + cumulative = 0 + for name in subdomain_names: + count = len(self.OF_cells_dict[name]) + self.subdomain_cell_offsets[name] = (cumulative, cumulative + count) + cumulative += count + + # Merge all pyvista meshes; clean() deduplicates coincident interface points + merged = pyvista.merge( + [self.OF_meshes_dict[name] for name in subdomain_names] + ).clean() + + OF_cells = merged.cells_dict.get(self.cell_type) + + if self.cell_type == 12: + shape = "hexahedron" + args_conn = np.tile(np.array([0, 1, 3, 2, 4, 5, 7, 6]), (len(OF_cells), 1)) + elif self.cell_type == 10: + shape = "tetrahedron" + args_conn = np.argsort(OF_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(OF_cells.shape[0])[:, None] + connectivity = OF_cells[rows, args_conn] + self.connectivities_dict["_global"] = connectivity + + 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 + self.mesh_vector_element = basix.ufl.element( + "Lagrange", cell_name, degree, shape=(3,) + ) + self.mesh_scalar_element = basix.ufl.element( + "Lagrange", cell_name, degree, shape=() + ) + + mesh_ufl = ufl.Mesh( + basix.ufl.element("Lagrange", cell_name, degree, shape=(3,)) + ) + self.dolfinx_meshes_dict["_global"] = create_mesh( + comm=MPI.COMM_WORLD, + cells=connectivity, + x=merged.points, + e=mesh_ufl, + ) + def _get_mesh( self, t: float, name: str, subdomain: str | None ) -> dolfinx.mesh.Mesh: @@ -295,29 +410,44 @@ def create_dolfinx_function_with_point_data( return u - def create_facet_meshtags( - self, t: float = 0, name: str = "U", subdomain: str | None = "default" - ) -> dolfinx.mesh.MeshTags: - """Creates a dolfinx.mesh.MeshTags for the facets of the mesh based on the - boundary patches in the OpenFOAM file. + 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, default to 0. - name: Name of the field in the OpenFOAM file, defaults to "U" for velocity - subdomain: Name of the subdmain in the OpenFOAM file, from which a field is - extracted + t: timestamp of the data to read, defaults to the first available time. Returns: the dolfinx MeshTags for the facets of the mesh """ - if self.multidomain is True: - raise NotImplementedError( - "Facet meshtags cannot be created for multidomain meshes, as the " - "boundary patches are not currently supported for multidomain meshes." - ) - mesh = self._get_mesh(t, name, subdomain) + t = t if t is not None else next(tv for tv in self.times if tv != 0) + self._read_with_pyvista_all(t=t) + + if self.multidomain: + if "_global" not in self.dolfinx_meshes_dict: + self._create_global_dolfinx_mesh() + mesh = self.dolfinx_meshes_dict["_global"] + else: + if "default" not in self.dolfinx_meshes_dict: + self._create_dolfinx_mesh(subdomain="default") + mesh = self.dolfinx_meshes_dict["default"] - boundary = self.OF_multiblock["boundary"] + # 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 @@ -333,10 +463,10 @@ def create_facet_meshtags( all_tags = [] patch_summary = {} - for i, patch_name in enumerate(boundary.keys()): + for i, (patch_name, patch_dataset) in enumerate(patches): facets, tags = tag_boundary_patch( mesh, - boundary[patch_name], + patch_dataset, i + 1, tree=tree, facet_indices=facet_indices, @@ -350,41 +480,105 @@ def create_facet_meshtags( 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: + # Build cell subdomain tag array in dolfinx ordering + total_pv_cells = sum(e for _, e in self.subdomain_cell_offsets.values()) + pv_cell_tags = np.zeros(total_pv_cells, dtype=np.int32) + for j, (sd_name, (start, end)) in enumerate( + self.subdomain_cell_offsets.items() + ): + pv_cell_tags[start:end] = j + 1 + num_cells = mesh.topology.index_map(mesh.topology.dim).size_local + dolfinx_cell_tags = pv_cell_tags[ + 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, - np.concatenate(all_facets), - np.concatenate(all_tags), + all_facets_arr[sort_idx].astype(np.int32), + all_tags_arr[sort_idx].astype(np.int32), ) - def create_cell_meshtags( - self, t: float = 0, name: str = "U", subdomain: str | None = "default" - ): - """Creates a dolfinx.mesh.MeshTags for the cells of the mesh based on the - cell data in the OpenFOAM file. + 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, default to 0. - name: Name of the field in the OpenFOAM file, defaults to "U" for velocity - subdomain: Name of the subdmain in the OpenFOAM file, from which a field is - extracted + t: timestamp of the data to read, defaults to the first available time. Returns: the dolfinx MeshTags for the cells of the mesh """ - if self.multidomain is True: - raise NotImplementedError( - "Facet meshtags cannot be created for multidomain meshes, as the " - "boundary patches are not currently supported for multidomain meshes." - ) - - mesh = self._get_mesh(t, name, subdomain) - - num_cells = mesh.topology.index_map(mesh.topology.dim).size_local - mesh_cell_indices = np.arange(num_cells, dtype=np.int32) - tags_volumes = np.full(num_cells, 1, dtype=np.int32) + t = t if t is not None else next(tv for tv in self.times if tv != 0) + self._read_with_pyvista_all(t=t) - return meshtags(mesh, mesh.topology.dim, mesh_cell_indices, tags_volumes) + if self.multidomain: + if "_global" not in self.dolfinx_meshes_dict: + self._create_global_dolfinx_mesh() + mesh = self.dolfinx_meshes_dict["_global"] + + total_pv_cells = sum(e for _, e in self.subdomain_cell_offsets.values()) + pv_cell_tags = np.zeros(total_pv_cells, dtype=np.int32) + for j, (_, (start, end)) in enumerate(self.subdomain_cell_offsets.items()): + pv_cell_tags[start:end] = j + 1 + num_cells = mesh.topology.index_map(mesh.topology.dim).size_local + cell_tag_array = pv_cell_tags[mesh.topology.original_cell_index[:num_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) + + 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: From 6f1f0c7f80d513f2ac986ec7cdfb2c1ccfef2714 Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 24 Apr 2026 10:18:21 -0400 Subject: [PATCH 07/15] correct tagging in multidomain cases --- src/foam2dolfinx/open_foam_reader.py | 57 ++++++++++++++++++---------- 1 file changed, 38 insertions(+), 19 deletions(-) diff --git a/src/foam2dolfinx/open_foam_reader.py b/src/foam2dolfinx/open_foam_reader.py index f4c1bf4..d76d720 100644 --- a/src/foam2dolfinx/open_foam_reader.py +++ b/src/foam2dolfinx/open_foam_reader.py @@ -77,6 +77,7 @@ def __init__(self, filename, cell_type: int = 12): self.dolfinx_meshes_dict = {} self.vertex_maps_dict = {} self.subdomain_cell_offsets = {} + self._pv_cell_subdomain_ids = None @property def cell_type(self): @@ -105,15 +106,20 @@ def _read_with_pyvista(self, t: float, subdomain: str | None = "default"): # Check if the reader has a multiblock dataset block named "internalMesh" if "internalMesh" not in self.OF_multiblock.keys(): self.multidomain = True - if subdomain not in self.OF_multiblock.keys(): + named_subdomains = [ + k for k in self.OF_multiblock.keys() if k != "defaultRegion" + ] + if subdomain not in named_subdomains: raise ValueError( f"Subdomain {subdomain} not found in the OpenFOAM file. " - f"Available subdomains: {self.OF_multiblock.keys()}" + f"Available subdomains: {named_subdomains}" ) # Extract the internal mesh if self.multidomain: for cell_array_name in self.OF_multiblock.keys(): + if cell_array_name == "defaultRegion": + continue self.OF_meshes_dict[cell_array_name] = self.OF_multiblock[ cell_array_name ]["internalMesh"] @@ -155,6 +161,8 @@ def _read_with_pyvista_all(self, t: float): if "internalMesh" not in self.OF_multiblock.keys(): self.multidomain = True for name in self.OF_multiblock.keys(): + if name == "defaultRegion": + continue self.OF_meshes_dict[name] = self.OF_multiblock[name]["internalMesh"] else: self.multidomain = False @@ -249,17 +257,26 @@ def _create_global_dolfinx_mesh(self): ) self.OF_cells_dict[name] = cells - # Record cumulative cell offsets before merging (cell order is preserved) + # 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) + # Merge all pyvista meshes; clean() deduplicates coincident interface points - merged = pyvista.merge( - [self.OF_meshes_dict[name] for name in subdomain_names] - ).clean() + 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) @@ -483,15 +500,11 @@ def create_facet_meshtags(self, t: float | None = None) -> dolfinx.mesh.MeshTags next_tag = len(patches) + 1 if self.multidomain: - # Build cell subdomain tag array in dolfinx ordering - total_pv_cells = sum(e for _, e in self.subdomain_cell_offsets.values()) - pv_cell_tags = np.zeros(total_pv_cells, dtype=np.int32) - for j, (sd_name, (start, end)) in enumerate( - self.subdomain_cell_offsets.items() - ): - pv_cell_tags[start:end] = j + 1 + # 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 = pv_cell_tags[ + dolfinx_cell_tags = self._pv_cell_subdomain_ids[ mesh.topology.original_cell_index[:num_cells] ] @@ -564,12 +577,15 @@ def create_cell_meshtags(self, t: float | None = None) -> dolfinx.mesh.MeshTags: self._create_global_dolfinx_mesh() mesh = self.dolfinx_meshes_dict["_global"] - total_pv_cells = sum(e for _, e in self.subdomain_cell_offsets.values()) - pv_cell_tags = np.zeros(total_pv_cells, dtype=np.int32) - for j, (_, (start, end)) in enumerate(self.subdomain_cell_offsets.items()): - pv_cell_tags[start:end] = j + 1 num_cells = mesh.topology.index_map(mesh.topology.dim).size_local - cell_tag_array = pv_cell_tags[mesh.topology.original_cell_index[:num_cells]] + 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") @@ -577,6 +593,9 @@ def create_cell_meshtags(self, t: float | None = None) -> dolfinx.mesh.MeshTags: 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) From 0d30a5e4255a9f580913155bd92d1f6a599ca210 Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 24 Apr 2026 10:21:08 -0400 Subject: [PATCH 08/15] fix test --- test/test_read_with_pyvista.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_read_with_pyvista.py b/test/test_read_with_pyvista.py index eb72f96..f1cd8a0 100644 --- a/test/test_read_with_pyvista.py +++ b/test/test_read_with_pyvista.py @@ -46,7 +46,7 @@ def test_error_rasied_when_subdomain_is_not_given_in_multidomain_case(tmpdir): ValueError, match=( r"Subdomain None not found in the OpenFOAM file\. " - r"Available subdomains: \['defaultRegion', 'fluid', 'solid']" + r"Available subdomains: \['fluid', 'solid']" ), ): my_of_reader._read_with_pyvista(t=20.0, subdomain=None) From 6eb54fcbf93831849e4108a3afeccf60ea9d6d25 Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 24 Apr 2026 10:26:30 -0400 Subject: [PATCH 09/15] extract connectivity and mesh methods --- src/foam2dolfinx/open_foam_reader.py | 113 +++++++++++---------------- 1 file changed, 45 insertions(+), 68 deletions(-) diff --git a/src/foam2dolfinx/open_foam_reader.py b/src/foam2dolfinx/open_foam_reader.py index d76d720..5c9dcd5 100644 --- a/src/foam2dolfinx/open_foam_reader.py +++ b/src/foam2dolfinx/open_foam_reader.py @@ -178,38 +178,47 @@ def _read_with_pyvista_all(self, t: float): f"Found {cell_types}" ) - def _create_dolfinx_mesh(self, subdomain: str | None = "default"): - """Creates a dolfinx.mesh.Mesh based on the elements within the OpenFOAM mesh""" + def _get_connectivity(self, cells: np.ndarray) -> tuple[str, np.ndarray]: + """Reorders a raw VTK cell array into DOLFINx vertex ordering. + + 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 @@ -219,26 +228,27 @@ 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 _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. - Reads cell data for any subdomains not yet in OF_cells_dict, records - cumulative cell offsets in subdomain_cell_offsets (used later for cell - and interface facet meshtags), merges all pyvista meshes with clean() - to deduplicate shared interface points, and creates a single dolfinx - mesh stored under the key "_global" in dolfinx_meshes_dict. + 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()) @@ -274,47 +284,14 @@ def _create_global_dolfinx_mesh(self): ) tagged_meshes.append(pv_mesh) - # Merge all pyvista meshes; clean() deduplicates coincident interface points 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) - - if self.cell_type == 12: - shape = "hexahedron" - args_conn = np.tile(np.array([0, 1, 3, 2, 4, 5, 7, 6]), (len(OF_cells), 1)) - elif self.cell_type == 10: - shape = "tetrahedron" - args_conn = np.argsort(OF_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(OF_cells.shape[0])[:, None] - connectivity = OF_cells[rows, args_conn] + shape, connectivity = self._get_connectivity(OF_cells) self.connectivities_dict["_global"] = connectivity - - 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 - self.mesh_vector_element = basix.ufl.element( - "Lagrange", cell_name, degree, shape=(3,) - ) - self.mesh_scalar_element = basix.ufl.element( - "Lagrange", cell_name, degree, shape=() - ) - - mesh_ufl = ufl.Mesh( - basix.ufl.element("Lagrange", cell_name, degree, shape=(3,)) - ) - self.dolfinx_meshes_dict["_global"] = create_mesh( - comm=MPI.COMM_WORLD, - cells=connectivity, - x=merged.points, - e=mesh_ufl, + self.dolfinx_meshes_dict["_global"] = self._build_dolfinx_mesh( + merged.points, connectivity, shape ) def _get_mesh( From bf8f21077c6239c9bec957ceb5f08a06c4c15988 Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 24 Apr 2026 10:32:16 -0400 Subject: [PATCH 10/15] redo read_with_pyvista --- src/foam2dolfinx/open_foam_reader.py | 88 ++++++++-------------------- 1 file changed, 23 insertions(+), 65 deletions(-) diff --git a/src/foam2dolfinx/open_foam_reader.py b/src/foam2dolfinx/open_foam_reader.py index 5c9dcd5..161ed1a 100644 --- a/src/foam2dolfinx/open_foam_reader.py +++ b/src/foam2dolfinx/open_foam_reader.py @@ -89,90 +89,48 @@ 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 - self.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 self.OF_multiblock.keys(): self.multidomain = True named_subdomains = [ k for k in self.OF_multiblock.keys() if k != "defaultRegion" ] - if subdomain not in named_subdomains: + 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: {named_subdomains}" ) - - # Extract the internal mesh - if self.multidomain: - for cell_array_name in self.OF_multiblock.keys(): - if cell_array_name == "defaultRegion": - continue - self.OF_meshes_dict[cell_array_name] = self.OF_multiblock[ - cell_array_name - ]["internalMesh"] - else: - self.OF_meshes_dict[subdomain] = self.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) - - # 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 _read_with_pyvista_all(self, t: float): - """Reads the OpenFOAM multiblock data at time t for all subdomains. - - Populates OF_multiblock and OF_meshes_dict for every subdomain without - requiring a specific subdomain to be named. For single-domain files, - also populates OF_cells_dict["default"]. For multidomain files, cell - data for each subdomain is read lazily by _create_global_dolfinx_mesh. - - Args: - t: timestamp of the data to read - """ - self.reader.set_active_time_value(t) - self.OF_multiblock = self.reader.read() - - if "internalMesh" not in self.OF_multiblock.keys(): - self.multidomain = True - for name in self.OF_multiblock.keys(): - if name == "defaultRegion": - continue + for name in named_subdomains: self.OF_meshes_dict[name] = self.OF_multiblock[name]["internalMesh"] else: self.multidomain = False self.OF_meshes_dict["default"] = self.OF_multiblock["internalMesh"] - OF_cell_type_dict = self.OF_meshes_dict["default"].cells_dict + 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["default"] = OF_cell_type_dict.get(self.cell_type) - if self.OF_cells_dict["default"] is None: + 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}" @@ -418,7 +376,7 @@ def create_facet_meshtags(self, t: float | None = None) -> dolfinx.mesh.MeshTags 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_all(t=t) + self._read_with_pyvista(t=t) if self.multidomain: if "_global" not in self.dolfinx_meshes_dict: @@ -547,7 +505,7 @@ def create_cell_meshtags(self, t: float | None = None) -> dolfinx.mesh.MeshTags: 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_all(t=t) + self._read_with_pyvista(t=t) if self.multidomain: if "_global" not in self.dolfinx_meshes_dict: From a590b42a8ac5c5c376af2b0d418fb0bf3c3f0889 Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 24 Apr 2026 10:35:15 -0400 Subject: [PATCH 11/15] ensure mesh --- src/foam2dolfinx/open_foam_reader.py | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/src/foam2dolfinx/open_foam_reader.py b/src/foam2dolfinx/open_foam_reader.py index 161ed1a..3e57e75 100644 --- a/src/foam2dolfinx/open_foam_reader.py +++ b/src/foam2dolfinx/open_foam_reader.py @@ -193,6 +193,21 @@ def _build_dolfinx_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]) @@ -378,14 +393,7 @@ def create_facet_meshtags(self, t: float | None = None) -> dolfinx.mesh.MeshTags 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"] - else: - if "default" not in self.dolfinx_meshes_dict: - self._create_dolfinx_mesh(subdomain="default") - mesh = self.dolfinx_meshes_dict["default"] + 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" From 949d4a58b794748655dec2f5afc1131a99d3f558 Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 24 Apr 2026 10:39:08 -0400 Subject: [PATCH 12/15] fix test --- test/test_read_with_pyvista.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_read_with_pyvista.py b/test/test_read_with_pyvista.py index f1cd8a0..b6ddbde 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,11 @@ 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"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") @pytest.mark.parametrize("subdomain", ["fluid", "solid"]) From 738d6240265c8d822aea6e98fdef811f9a0d5902 Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 24 Apr 2026 10:39:17 -0400 Subject: [PATCH 13/15] simplify function --- src/foam2dolfinx/helpers.py | 29 ++++++----------------------- 1 file changed, 6 insertions(+), 23 deletions(-) diff --git a/src/foam2dolfinx/helpers.py b/src/foam2dolfinx/helpers.py index 5b0963b..331dcdb 100644 --- a/src/foam2dolfinx/helpers.py +++ b/src/foam2dolfinx/helpers.py @@ -1,5 +1,4 @@ import numpy as np -from dolfinx.mesh import exterior_facet_indices from scipy.spatial import cKDTree @@ -9,9 +8,9 @@ def tag_boundary_patch( patch_id, tol=1e-6, *, - tree=None, - facet_indices=None, - facet_vertices=None, + 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. @@ -20,35 +19,19 @@ def tag_boundary_patch( 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: optional pre-built cKDTree on mesh geometry points — pass this when - calling for multiple patches to avoid rebuilding the tree each time - facet_indices: optional pre-computed exterior facet indices — pass together - with facet_vertices when calling for multiple patches - facet_vertices: optional pre-computed 2D array (n_facets, n_verts_per_facet) + 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 """ - fdim = dolfinx_mesh.topology.dim - 1 - - if facet_indices is None or facet_vertices is None: - dolfinx_mesh.topology.create_connectivity(fdim, 0) - dolfinx_mesh.topology.create_connectivity(0, fdim) - dolfinx_mesh.topology.create_connectivity(fdim, dolfinx_mesh.topology.dim) - facet_indices = exterior_facet_indices(dolfinx_mesh.topology) - c_to_v = dolfinx_mesh.topology.connectivity(fdim, 0) - facet_vertices = np.vstack([c_to_v.links(f) for f in facet_indices]) - - if tree is None: - tree = cKDTree(dolfinx_mesh.geometry.x) - 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) ) - # boolean mask over all vertices, then check all vertices of each facet at once vertex_matched = np.zeros(len(dolfinx_mesh.geometry.x), dtype=bool) if len(matched_idx): vertex_matched[matched_idx] = True From e6b8ce196241fecd0e410ed6d22636abf96b6ccc Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 24 Apr 2026 10:52:09 -0400 Subject: [PATCH 14/15] additional test --- test/conftest.py | 16 +++++++ test/test_get_connectivity.py | 72 +++++++++++++++++++++++++++++ test/test_meshtags.py | 81 +++++++++++++++++++++++++++++++++ test/test_read_with_pyvista.py | 42 +++++++++++++++++ test/test_tag_boundary_patch.py | 71 +++++++++++++++++++++++++++++ 5 files changed, 282 insertions(+) create mode 100644 test/conftest.py create mode 100644 test/test_get_connectivity.py create mode 100644 test/test_meshtags.py create mode 100644 test/test_tag_boundary_patch.py 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_get_connectivity.py b/test/test_get_connectivity.py new file mode 100644 index 0000000..e03fcfa --- /dev/null +++ b/test/test_get_connectivity.py @@ -0,0 +1,72 @@ +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..48add52 --- /dev/null +++ b/test/test_meshtags.py @@ -0,0 +1,81 @@ +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 b6ddbde..69944e4 100644 --- a/test/test_read_with_pyvista.py +++ b/test/test_read_with_pyvista.py @@ -52,6 +52,48 @@ def test_error_rasied_when_wrong_subdomain_given_in_multidomain_case(tmpdir): 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"]) def test_read_with_pyvista_finds_all_mesh_data(tmpdir, subdomain): zip_path = Path("test/data/test_2Regions.zip") diff --git a/test/test_tag_boundary_patch.py b/test/test_tag_boundary_patch.py new file mode 100644 index 0000000..cd37731 --- /dev/null +++ b/test/test_tag_boundary_patch.py @@ -0,0 +1,71 @@ +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 From 6d228eee9689c47bc1c90330fcac663640109900 Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 24 Apr 2026 10:54:05 -0400 Subject: [PATCH 15/15] lint ruff Co-authored-by: Copilot --- test/test_create_function.py | 6 +++-- test/test_get_connectivity.py | 20 +++++++++++---- test/test_meshtags.py | 4 ++- test/test_tag_boundary_patch.py | 44 ++++++++++++++++++++++++++------- 4 files changed, 57 insertions(+), 17 deletions(-) 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 index e03fcfa..3e0d758 100644 --- a/test/test_get_connectivity.py +++ b/test/test_get_connectivity.py @@ -20,15 +20,21 @@ def tet_reader(): "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) +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) + 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 @@ -55,7 +61,9 @@ def test_get_connectivity_tetrahedron_sorts_vertices(tet_reader): 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"]) + 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 ) @@ -65,7 +73,9 @@ def test_build_dolfinx_mesh_returns_dolfinx_mesh(hex_reader): @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"]) + 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 ) diff --git a/test/test_meshtags.py b/test/test_meshtags.py index 48add52..95455d7 100644 --- a/test/test_meshtags.py +++ b/test/test_meshtags.py @@ -55,7 +55,9 @@ def test_create_facet_meshtags_dim_equals_fdim(two_regions_reader): @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): +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 diff --git a/test/test_tag_boundary_patch.py b/test/test_tag_boundary_patch.py index cd37731..bf1362a 100644 --- a/test/test_tag_boundary_patch.py +++ b/test/test_tag_boundary_patch.py @@ -12,7 +12,10 @@ @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 + 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) @@ -29,7 +32,12 @@ def test_tag_boundary_patch_returns_empty_when_no_points_match(unit_cube_topolog 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 + mesh, + patch, + 1, + tree=tree, + facet_indices=facet_indices, + facet_vertices=facet_vertices, ) assert len(matched) == 0 assert len(tags) == 0 @@ -40,7 +48,12 @@ def test_tag_boundary_patch_tags_equal_patch_id(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 + mesh, + patch, + 42, + tree=tree, + facet_indices=facet_indices, + facet_vertices=facet_vertices, ) assert np.all(tags == 42) @@ -50,22 +63,35 @@ def test_tag_boundary_patch_matched_facets_are_subset_of_exterior(unit_cube_topo 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 + 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) -]) +@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 + mesh, + patch, + 1, + tree=tree, + facet_indices=facet_indices, + facet_vertices=facet_vertices, ) assert matched.dtype == np.int32 assert tags.dtype == np.int32