Added load_mesh_from_file function to support external 3D meshes. (#5162)

* adding load ext mesh functions and few tests

* updated changelog

* adding meshio as optional dependency

* minor doc changes

* adding warnings for non tested file format, docstring instructions & changning name of submesh

* ensuring coverage

---------

Co-authored-by: Robert Timms <[email protected]>

Author: Rishab87

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+- Added helper functions to import external 3D meshes in PyBaMM ([#5162](https://github.com/pybamm-team/PyBaMM/pull/5162))
 - Added support for algebraic and differential surface form in composite models. ([#5165](https://github.com/pybamm-team/PyBaMM/pull/5165))
 - Adds a composite electrode electrode soh model ([#5160](https://github.com/pybamm-team/PyBaMM/pull/5129))
 - Generalises `set_initial_soc` to `set_initial_state` and adds Thevenin initial state setting. ([#5129](https://github.com/pybamm-team/PyBaMM/pull/5129))

pyproject.toml

@@ -118,6 +118,7 @@ jax = [
 # Contains all optional dependencies, except for jax, and dev dependencies
 all = [
     "scikit-fem>=8.1.0",
+    "meshio>=5.3.0",
     "pybamm[examples,plot,cite,bpx,tqdm]",
 ]
 

src/pybamm/__init__.py

@@ -150,6 +150,7 @@
 from .meshes.scikit_fem_submeshes_3d import (
     ScikitFemSubMesh3D,
     ScikitFemGenerator3D,
+    UserSuppliedSubmesh3D,
 )
 
 # Serialisation

src/pybamm/meshes/scikit_fem_submeshes_3d.py

@@ -1,3 +1,5 @@
+import os
+
 import numpy as np
 
 import pybamm
@@ -234,6 +236,87 @@ def __call__(self, lims, npts):
         return submesh
 
 
+class UserSuppliedSubmesh3D(pybamm.MeshGenerator):
+    """
+    A mesh generator that loads a 3D mesh from an external file (e.g., VTK).
+
+    This allows using high-quality meshes from any external tool. The mesh
+    file must contain integer tags identifying the different physical domains
+    and boundaries. This generator works directly with your existing
+    'ScikitFemSubMesh3D' and 'ScikitFiniteElement3D' classes.
+
+    Parameters
+    ----------
+    file_path : str
+        Path to the mesh file.
+    boundary_mapping : dict
+        Maps PyBaMM boundary names to the integer tags in the mesh file.
+    domain_mapping : dict
+        Maps PyBaMM domain names to the integer tags in the mesh file.
+    domain_tag_name : str, optional
+        The name of the data array in the mesh file for domain tags.
+    boundary_tag_name : str, optional
+        The name of the data array in the mesh file for boundary tags.
+    coord_sys : str, optional
+        The coordinate system ("cartesian" or "cylindrical polar").
+
+    Notes
+    -----
+    The external mesh file must meet the following criteria to be compatible:
+
+        1.  **Element Types**: The 3D volume must be meshed with 4-node
+            **tetrahedral** elements. The 2D boundaries must be meshed with
+            3-node **triangular** elements.
+
+        2.  **Physical Tags**: The file must contain integer tags to identify
+            physical groups.
+            - Each tetrahedron in the volume must be tagged with an integer
+            that corresponds to a PyBaMM domain.
+            - Each triangle on a boundary must be tagged with an integer that
+            corresponds to a PyBaMM boundary.
+
+        3.  **File Format**: Any format supported by the `meshio` library is
+            acceptable. For maximum reliability in preserving physical tags,
+            the **Gmsh `.msh` Version 2.2 (ASCII)** format is recommended.
+    """
+
+    def __init__(
+        self,
+        file_path,
+        boundary_mapping,
+        domain_mapping,
+        domain_tag_name=None,
+        boundary_tag_name=None,
+        coord_sys=None,
+    ):
+        super().__init__(ScikitFemSubMesh3D)
+        self.file_path = file_path
+        self.boundary_mapping = boundary_mapping
+        self.domain_mapping = domain_mapping
+        self.domain_tag_name = domain_tag_name
+        self.boundary_tag_name = boundary_tag_name
+
+        if coord_sys:
+            self.coord_sys = coord_sys
+        else:
+            if any("r_min" in k or "r_max" in k for k in boundary_mapping.keys()):
+                self.coord_sys = "cylindrical polar"
+            else:
+                self.coord_sys = "cartesian"
+
+    def __call__(self, lims, npts):
+        skfem_mesh = ScikitFemSubMesh3D.load_mesh_from_file(
+            self.file_path,
+            self.boundary_mapping,
+            self.domain_mapping,
+            self.domain_tag_name,
+            self.boundary_tag_name,
+        )
+        nodes = skfem_mesh.p.T
+        elements = skfem_mesh.t.T
+        return self.submesh_type(skfem_mesh, nodes, elements, self.coord_sys)
+
+
 class ScikitFemSubMesh3D(pybamm.SubMesh):
     """
     A 3D submesh class for unstructured tetrahedral meshes generated by scikit-fem.
@@ -323,3 +406,255 @@ def _from_json(cls, json_dict):
         elements = np.array(json_dict["elements"])
         mesh = skfem.MeshTet(nodes.T, elements.T)
         return cls(mesh, nodes, elements, json_dict["coord_sys"])
+
+    @staticmethod
+    def load_mesh_from_file(
+        file_path,
+        boundary_mapping,
+        domain_mapping,
+        domain_tag_name=None,
+        boundary_tag_name=None,
+    ):
+        """
+        Loads a generic mesh file and converts it to a scikit-fem Mesh object.
+
+        This function is designed to work with meshes from most software packages
+        including GMSH, Ansys, Abaqus, FEniCS, etc. It automatically detects
+        common tag naming conventions or allows manual specification.
+
+        Parameters
+        ----------
+        file_path : str
+            The path to the mesh file (It supports .msh, .xdmf, may also work with others but not tested).
+        boundary_mapping : dict
+            Maps PyBaMM boundary names to integer tags (e.g., {"r_inner": 1}).
+        domain_mapping : dict
+            Maps PyBaMM domain names to integer tags (e.g., {"current collector": 5}).
+        domain_tag_name : str, optional
+            The name of the cell data array that contains domain tags for elements.
+        boundary_tag_name : str, optional
+            The name of the cell data array that contains boundary tags for facets.
+        """
+        meshio = pybamm.util.import_optional_dependency("meshio")
+        skfem = pybamm.util.import_optional_dependency("skfem")
+
+        tested_formats = [".msh", ".xdmf"]
+        _, file_ext = os.path.splitext(file_path)
+        if file_ext not in tested_formats:  # pragma: no cover
+            pybamm.logger.warning(
+                f"File format '{file_ext}' has not been explicitly tested and may not "
+                "work correctly. Recommended formats are .msh and .xdmf."
+            )
+
+        try:
+            m = meshio.read(file_path)
+        except Exception as e:
+            raise pybamm.GeometryError(
+                f"Could not read mesh file '{file_path}': {e}"
+            ) from e
+
+        nodes = m.points
+        if nodes.shape[1] != 3:
+            raise pybamm.GeometryError(
+                f"Mesh must be 3D, got points with shape {nodes.shape}"
+            )
+
+        def find_tag_name(
+            user_provided_name, default_candidates, available_data, context=""
+        ):
+            """
+            Enhanced tag name finder that works with most mesh formats.
+            """
+            # If user explicitly provided a name, use it
+            if user_provided_name:
+                if user_provided_name in available_data:
+                    return user_provided_name
+                else:
+                    available_keys = list(available_data.keys())
+                    raise pybamm.GeometryError(
+                        f"User-specified {context} tag name '{user_provided_name}' not found. "
+                        f"Available arrays: {available_keys}"
+                    )
+
+            for candidate in default_candidates:
+                if candidate in available_data:
+                    return candidate
+
+            integer_arrays = [
+                name
+                for name, data in available_data.items()
+                if np.issubdtype(data.dtype, np.integer)
+            ]
+
+            if len(integer_arrays) == 1:
+                return integer_arrays[0]
+
+            priority_patterns = ["physical", "tag", "id", "group", "material", "region"]
+            for pattern in priority_patterns:
+                for name in integer_arrays:
+                    if pattern.lower() in name.lower():
+                        return name
+
+            raise pybamm.GeometryError(
+                f"Could not automatically detect {context} tag array in '{file_path}'. "
+                f"Available arrays: {list(available_data.keys())}. "
+                f"Integer arrays found: {integer_arrays}. "
+                f"Please specify manually using '{context.replace(' ', '_')}_tag_name'."
+            )
+
+        try:
+            tet_cells = m.get_cells_type("tetra")
+            if len(tet_cells) == 0:
+                raise pybamm.GeometryError(
+                    f"No tetrahedral elements found in '{file_path}'"
+                )
+
+            tet_cell_index = None
+            for i, cell_block in enumerate(m.cells):
+                if cell_block.type == "tetra":
+                    tet_cell_index = i
+                    break
+
+            if tet_cell_index is None:
+                raise pybamm.GeometryError("Could not find tetrahedra in cell list")
+
+            tet_cell_data_dict = {}
+            for name, data_list in m.cell_data.items():
+                if len(data_list) > tet_cell_index:
+                    data = data_list[tet_cell_index]
+                    if len(data) == len(tet_cells):
+                        tet_cell_data_dict[name] = data
+
+            domain_tag_name = find_tag_name(
+                domain_tag_name,
+                [
+                    "gmsh:physical",
+                    "physical",
+                    "PhysicalGroup",
+                    "MaterialId",
+                    "CellEntityIds",
+                    "gmsh:geometrical",
+                ],
+                tet_cell_data_dict,
+                "domain",
+            )
+
+            elements = tet_cells
+            tet_cell_data = tet_cell_data_dict[domain_tag_name]
+
+            subdomains = {}
+            for name, tag in domain_mapping.items():
+                matching_elements = np.where(tet_cell_data == tag)[0]
+                if len(matching_elements) == 0:
+                    pybamm.logger.warning(
+                        f"No elements found for domain '{name}' with tag {tag}"
+                    )
+                else:
+                    subdomains[name] = matching_elements
+                    pybamm.logger.info(
+                        f"Found {len(matching_elements)} elements for domain '{name}' (tag {tag})"
+                    )
+
+        except Exception as e:
+            raise pybamm.GeometryError(
+                f"Failed to extract tetrahedral elements from '{file_path}': {e}"
+            ) from e
+
+        skfem_mesh = skfem.MeshTet(nodes.T, elements.T)
+
+        boundaries = {}
+        try:
+            tri_cells = m.get_cells_type("triangle")
+
+            if len(tri_cells) > 0:
+                tri_cell_indices = []
+                for i, cell_block in enumerate(m.cells):
+                    if cell_block.type == "triangle":
+                        tri_cell_indices.append(i)
+
+                # Combine all triangle cell data
+                all_tri_cells = []
+                all_tri_data = []
+
+                for idx in tri_cell_indices:
+                    block_cells = m.cells[idx].data
+                    all_tri_cells.extend(block_cells)
+
+                    for name, data_list in m.cell_data.items():
+                        if len(data_list) > idx:
+                            if name not in [arr_name for arr_name, _ in all_tri_data]:
+                                all_tri_data.append((name, []))
+                            for arr_name, arr_data in all_tri_data:
+                                if arr_name == name:
+                                    arr_data.extend(data_list[idx])
+                                    break
+
+                all_tri_cells = np.array(all_tri_cells)
+
+                tri_cell_data_dict = {}
+                for name, data in all_tri_data:
+                    tri_cell_data_dict[name] = np.array(data)
+
+                if tri_cell_data_dict:
+                    # Find boundary tag name
+                    boundary_tag_name = find_tag_name(
+                        boundary_tag_name,
+                        [
+                            "gmsh:physical",
+                            "physical",
+                            "PhysicalGroup",
+                            "FaceEntityIds",
+                            "gmsh:geometrical",
+                        ],
+                        tri_cell_data_dict,
+                        "boundary",
+                    )
+
+                    tri_cell_data = tri_cell_data_dict[boundary_tag_name]
+
+                    # Map triangle facets to mesh boundary facets
+                    bnd_facets_indices = skfem_mesh.boundary_facets()
+                    bnd_facets_nodes = skfem_mesh.facets[:, bnd_facets_indices]
+
+                    # Create mapping from node sets to facet indices
+                    facet_map = {}
+                    for i in range(len(bnd_facets_indices)):
+                        node_set = frozenset(bnd_facets_nodes[:, i])
+                        facet_map[node_set] = bnd_facets_indices[i]
+
+                    # Find boundary facets for each tag
+                    for name, tag in boundary_mapping.items():
+                        tagged_triangles = all_tri_cells[tri_cell_data == tag]
+                        facet_indices = []
+
+                        for tri in tagged_triangles:
+                            tri_set = frozenset(tri)
+                            if tri_set in facet_map:
+                                facet_indices.append(facet_map[tri_set])
+
+                        if facet_indices:
+                            boundaries[name] = np.array(facet_indices)
+                            pybamm.logger.info(
+                                f"Found {len(facet_indices)} boundary facets for '{name}' (tag {tag})"
+                            )
+                        else:
+                            pybamm.logger.warning(
+                                f"No boundary facets found for '{name}' with tag {tag}"
+                            )
+
+        except Exception as e:
+            pybamm.logger.warning(
+                f"Could not extract boundary information from '{file_path}': {e}"
+            )
+
+        if boundaries:
+            skfem_mesh = skfem_mesh.with_boundaries(boundaries)
+        if subdomains:
+            skfem_mesh = skfem_mesh.with_subdomains(subdomains)
+
+        pybamm.logger.info(
+            f"Successfully loaded mesh: {len(nodes)} nodes, {len(elements)} elements, "
+            f"{len(boundaries)} boundary groups, {len(subdomains)} subdomains"
+        )
+
+        return skfem_mesh