Add unstructured mesh infrastructure

CHANGELOG.md

@@ -1,5 +1,9 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
+## Features
+
+- Added unstructured mesh support (`UnstructuredSubMesh`, generators, and interface coupling) for arbitrary 2D/3D domains. ([#5687](https://github.com/pybamm-team/PyBaMM/pull/5687))
+
 ## Bug fixes
 
 - `BatchStudy.solve` no longer ignores its `solver` argument: previously the loop over study inputs shadowed it, so a caller-supplied solver was silently dropped. A solver from `BatchStudy(solvers=...)` still takes precedence. ([#5677](https://github.com/pybamm-team/PyBaMM/pull/5677))

docs/source/api/meshes/index.rst

@@ -8,3 +8,4 @@ Meshes
   one_dimensional_submeshes
   two_dimensional_submeshes
   three_dimensional_submeshes
+  unstructured_submeshes

docs/source/api/meshes/unstructured_submeshes.rst

@@ -0,0 +1,16 @@
+Unstructured Sub Meshes
+=======================
+
+.. autoclass:: pybamm.UnstructuredSubMesh
+    :members:
+
+.. autoclass:: pybamm.UnstructuredMeshGenerator
+    :members:
+
+.. autoclass:: pybamm.UserSuppliedUnstructuredMesh
+    :members:
+
+.. autoclass:: pybamm.TaggedSubMeshGenerator
+    :members:
+
+.. autofunction:: pybamm.compute_interface_data

packages/pybamm/src/pybamm/__init__.py

@@ -168,6 +168,15 @@
     UserSuppliedSubmesh3D,
 )
 
+
+from .meshes.unstructured_submesh import (
+    UnstructuredSubMesh,
+    UnstructuredMeshGenerator,
+    UserSuppliedUnstructuredMesh,
+    TaggedSubMeshGenerator,
+    compute_interface_data,
+)
+
 # Serialisation
 from .models.base_model import load_model
 

packages/pybamm/src/pybamm/meshes/__init__.py

@@ -1,2 +1,3 @@
 __all__ = ['meshes', 'one_dimensional_submeshes', 'scikit_fem_submeshes',
-           'zero_dimensional_submesh', 'scikit_fem_submeshes_3d']
+           'zero_dimensional_submesh', 'scikit_fem_submeshes_3d',
+           'unstructured_submesh']

packages/pybamm/src/pybamm/meshes/meshes.py

@@ -170,6 +170,9 @@ def __init__(self, geometry, submesh_types, var_pts):
             self[domain] = submesh_types[domain](geometry[domain], submesh_pts[domain])
             self.base_domains.append(domain)
 
+        # compute interface data for unstructured meshes
+        self._compute_unstructured_interfaces()
+
         # add ghost meshes
         self.add_ghost_meshes()
 
@@ -214,6 +217,8 @@ def combine_submeshes(self, *submeshnames):
                 raise pybamm.GeometryError(
                     "Cannot combine submeshes of different dimensions"
                 )
+            elif isinstance(self[submeshnames[i]], pybamm.UnstructuredSubMesh):
+                pass
             elif self[submeshnames[i]].dimension == 2:
                 if "left" in submeshnames[i] or "right" in submeshnames[i + 1]:
                     # Make sure that the lr edges are aligned
@@ -264,7 +269,11 @@ def combine_submeshes(self, *submeshnames):
                 )
 
         coord_sys = self[submeshnames[0]].coord_sys
-        if self[submeshnames[0]].dimension == 1:
+        if isinstance(self[submeshnames[0]], pybamm.UnstructuredSubMesh):
+            return pybamm.UnstructuredSubMesh.combine(
+                [self[name] for name in submeshnames]
+            )
+        elif self[submeshnames[0]].dimension == 1:
             combined_submesh_edges = np.concatenate(
                 [self[submeshnames[0]].edges]
                 + [self[submeshname].edges[1:] for submeshname in submeshnames[1:]]
@@ -349,6 +358,35 @@ def combine_submeshes(self, *submeshnames):
                 submesh.internal_boundaries.append(self[submeshname].edges_lr[0] + min)
         return submesh
 
+    def _compute_unstructured_interfaces(self):
+        """
+        For adjacent domains backed by :class:`UnstructuredSubMesh`, compute
+        and store interface coupling data.
+        """
+        unstructured_domains = [
+            d
+            for d in self.base_domains
+            if isinstance(self[d], pybamm.UnstructuredSubMesh)
+        ]
+        for i in range(len(unstructured_domains) - 1):
+            left_name = unstructured_domains[i]
+            right_name = unstructured_domains[i + 1]
+            left_mesh = self[left_name]
+            right_mesh = self[right_name]
+            if (
+                "right" in left_mesh.boundary_faces
+                and "left" in right_mesh.boundary_faces
+            ):
+                try:
+                    pybamm.compute_interface_data(
+                        left_mesh,
+                        right_mesh,
+                        left_name=left_name,
+                        right_name=right_name,
+                    )
+                except pybamm.GeometryError:
+                    pass
+
     def add_ghost_meshes(self):
         """
         Create meshes for potential ghost nodes on either side of each submesh, using
@@ -365,7 +403,8 @@ def add_ghost_meshes(self):
                     submesh,
                     pybamm.SubMesh0D
                     | pybamm.ScikitSubMesh2D
-                    | pybamm.ScikitFemSubMesh3D,
+                    | pybamm.ScikitFemSubMesh3D
+                    | pybamm.UnstructuredSubMesh,
                 )
             )
         ]