Qol auto detect surface to volume

examples/advection_diffusion_multi_material.py

@@ -127,11 +127,6 @@ def velocity_func(x):
 
 # we should be able to automate this
 my_model.interfaces = [F.Interface(6, (bottom_domain, top_domain))]
-my_model.surface_to_volume = {
-    top_surface: top_domain,
-    bottom_surface: bottom_domain,
-    inlet_surf: bottom_domain,
-}
 
 H = F.Species("H", mobile=True)
 

examples/multi_material_1d.py

@@ -55,10 +55,6 @@
     left_surface,
     right_surface,
 ]
-my_model.surface_to_volume = {
-    right_surface: right_domain,
-    left_surface: left_domain,
-}
 
 H = F.Species("H", mobile=True)
 trapped_H = F.Species("H_trapped", mobile=False)

examples/multi_material_2d.py

@@ -85,10 +85,6 @@ def half(x):
 
 # we should be able to automate this
 my_model.interfaces = [F.Interface(5, (bottom_domain, top_domain))]
-my_model.surface_to_volume = {
-    top_surface: top_domain,
-    bottom_surface: bottom_domain,
-}
 
 H = F.Species("H", mobile=True)
 

examples/multi_material_2d_bis.py

@@ -85,10 +85,6 @@ def half(x):
 
 # we should be able to automate this
 my_model.interfaces = [F.Interface(5, (bottom_domain, top_domain))]
-my_model.surface_to_volume = {
-    top_surface: top_domain,
-    bottom_surface: bottom_domain,
-}
 
 H = F.Species("H", mobile=True)
 trapped_H = F.Species("H_trapped", mobile=False)

examples/multi_material_transient.py

@@ -50,10 +50,6 @@
     left_surface,
     right_surface,
 ]
-my_model.surface_to_volume = {
-    right_surface: right_domain,
-    left_surface: left_domain,
-}
 
 H = F.Species("H", mobile=True)
 trapped_H = F.Species("H_trapped", mobile=False)

examples/multi_material_with_one_material.py

@@ -22,7 +22,6 @@
 right_surface = F.SurfaceSubdomain1D(id=2, x=vertices[-1])
 
 my_model.subdomains = [subdomain, left_surface, right_surface]
-my_model.surface_to_volume = {right_surface: subdomain, left_surface: subdomain}
 
 H = F.Species("H", mobile=True)
 trapped_H = F.Species("H_trapped", mobile=False)

src/festim/__init__.py

@@ -79,5 +79,6 @@
     VolumeSubdomain,
     VolumeSubdomain1D,
     find_volume_from_id,
+    map_surface_to_volume_subdomains,
 )
 from .trap import Trap

src/festim/hydrogen_transport_problem.py

@@ -1055,7 +1055,6 @@ def __init__(
         exports=None,
         traps=None,
         interfaces: list[_subdomain.Interface] | None = None,
-        surface_to_volume: dict | None = None,
         petsc_options: dict | None = None,
     ):
         """Class for a multi-material hydrogen transport problem
@@ -1093,7 +1092,7 @@ def __init__(
             petsc_options=petsc_options,
         )
         self.interfaces = interfaces or []
-        self.surface_to_volume = surface_to_volume or {}
+        self.surface_to_volume = {}
         self.subdomain_to_species = {}  # maps subdomain to species defined in it
 
     @property
@@ -1141,6 +1140,27 @@ def initialise(self):
                 raise TypeError("subdomains attribute in Species should be list")
 
         self.define_meshtags_and_measures()
+        if self.surface_to_volume:
+            # tell users that this is no longer required
+            warnings.warn(
+                f"The surface_to_volume attribute of the {self.__class__.__name__}"
+                " class is no longer required and can be removed."
+                "The mapping between surface and volume subdomains is now done"
+                "automatically based on the connectivity of the mesh and the meshtags",
+                DeprecationWarning,
+            )
+        else:
+            facet_to_cell = self.mesh.mesh.topology.connectivity(
+                self.mesh.mesh.topology.dim - 1, self.mesh.mesh.topology.dim
+            )
+            self.surface_to_volume = _subdomain.map_surface_to_volume_subdomains(
+                ft=self.facet_meshtags,
+                ct=self.volume_meshtags,
+                facet_to_cell=facet_to_cell,
+                volume_subdomains=self.volume_subdomains,
+                surface_subdomains=self.surface_subdomains,
+                comm=self.mesh.mesh.comm,
+            )
 
         # create submeshes and transfer meshtags to subdomains
         for subdomain in self.volume_subdomains:

src/festim/subdomain/__init__.py

@@ -1,6 +1,10 @@
 from .interface import Interface
 from .surface_subdomain import SurfaceSubdomain, SurfaceSubdomain1D
-from .volume_subdomain import VolumeSubdomain, VolumeSubdomain1D
+from .volume_subdomain import (
+    VolumeSubdomain,
+    VolumeSubdomain1D,
+    map_surface_to_volume_subdomains,
+)
 
 __all__ = [
     "Interface",
@@ -9,4 +13,5 @@
     "SurfaceSubdomain1D",
     "VolumeSubdomain",
     "VolumeSubdomain1D",
+    "map_surface_to_volume_subdomains",
 ]

src/festim/subdomain/volume_subdomain.py

@@ -9,6 +9,7 @@
 from scifem.mesh import transfer_meshtags_to_submesh
 
 from festim.material import Material
+from festim.subdomain.surface_subdomain import SurfaceSubdomain
 
 # Define the appropriate method based on the version
 try:
@@ -185,3 +186,89 @@ def find_volume_from_id(id: int, volumes: list):
         if vol.id == id:
             return vol
     raise ValueError(f"id {id} not found in list of volumes")
+
+
+def map_surface_to_volume_subdomains(
+    ft: dolfinx.mesh.MeshTags,
+    ct: dolfinx.mesh.MeshTags,
+    facet_to_cell: dolfinx.cpp.graph.AdjacencyList_int32,
+    volume_subdomains: list[VolumeSubdomain],
+    surface_subdomains: list[SurfaceSubdomain],
+    comm=None,
+) -> dict[SurfaceSubdomain, VolumeSubdomain]:
+    """Maps surface subdomains to volume subdomains based on the facet and cell meshtags
+    and the facet to cell connectivity.
+
+
+    Raises:
+        AssertionError: if a surface subdomain is connected to multiple volume subdomains
+
+    Args:
+        ft: the facet meshtags of the parent mesh
+        ct: the cell meshtags of the parent mesh
+        facet_to_cell: the facet to cell connectivity of the parent mesh
+        volume_subdomains: the list of volume subdomains
+        surface_subdomains: the list of surface subdomains
+        comm: MPI communicator (required for parallel runs)
+
+    Returns:
+        dict[SurfaceSubdomain, VolumeSubdomain]: a dictionary mapping surface subdomains
+            to volume subdomains
+    """
+
+    # get connected cells for tagged facets
+    start_indices = facet_to_cell.offsets[ft.indices]
+    end_indices = facet_to_cell.offsets[ft.indices + 1]
+    num_connections = end_indices - start_indices
+
+    # A facet is connected to at most 2 cells (boundary = 1, interior = 2)
+    cell_ids_0 = facet_to_cell.array[start_indices]
+    has_second_cell = num_connections == 2
+    cell_ids_1 = facet_to_cell.array[start_indices[has_second_cell] + 1]
+
+    connected_cells = np.concatenate([cell_ids_0, cell_ids_1])
+    connected_facet_tags = np.concatenate([ft.values, ft.values[has_second_cell]])
+
+    # map connected cells to their cell tags
+    sort_idx = np.argsort(ct.indices)
+    sorted_ct_indices = ct.indices[sort_idx]
+    sorted_ct_values = ct.values[sort_idx]
+
+    idx = np.searchsorted(sorted_ct_indices, connected_cells)
+    # mask out-of-bounds
+    valid = idx < len(sorted_ct_indices)
+    # of those in bounds, check if they actually match
+    valid[valid] = sorted_ct_indices[idx[valid]] == connected_cells[valid]
+
+    valid_cell_tags = sorted_ct_values[idx[valid]]
+    valid_facet_tags = connected_facet_tags[valid]
+
+    unique_pairs = np.unique(np.vstack((valid_facet_tags, valid_cell_tags)).T, axis=0)
+    if comm is not None and comm.size > 1:
+        all_pairs = comm.allgather(unique_pairs)
+        non_empty = [p for p in all_pairs if len(p) > 0]
+        if non_empty:
+            unique_pairs = np.unique(np.vstack(non_empty), axis=0)
+
+    surface_tag_to_subdomain = {s.id: s for s in surface_subdomains}
+    volume_tag_to_subdomain = {v.id: v for v in volume_subdomains}
+
+    surface_to_subdomain = {}
+
+    for s_tag, v_tag in unique_pairs:
+        dolfinx.log.log(
+            dolfinx.log.LogLevel.INFO,
+            f"Facet tag {s_tag} is connected to cell tag {v_tag}",
+        )
+        s_subdomain = surface_tag_to_subdomain.get(s_tag)
+        v_subdomain = volume_tag_to_subdomain.get(v_tag)
+
+        if s_subdomain and v_subdomain:
+            if s_subdomain in surface_to_subdomain:
+                assert surface_to_subdomain[s_subdomain] == v_subdomain, (
+                    f"Surface subdomain {s_subdomain.id} is connected to multiple volume subdomains: "
+                    f"{surface_to_subdomain[s_subdomain].id} and {v_subdomain.id}"
+                )
+            else:
+                surface_to_subdomain[s_subdomain] = v_subdomain
+    return surface_to_subdomain

test/system_tests/test_advection_term_integration.py

@@ -250,11 +250,6 @@ def velocity_func(x):
 
     # we should be able to automate this
     my_model.interfaces = [F.Interface(6, (bottom_domain, top_domain))]
-    my_model.surface_to_volume = {
-        top_surface: top_domain,
-        bottom_surface: bottom_domain,
-        inlet_surf: bottom_domain,
-    }
 
     H = F.Species("H", mobile=True)
 

test/system_tests/test_cylindrical_coordinates.py

@@ -115,10 +115,6 @@ def c_exact_right(x):
     my_model.subdomains = [vol_1, vol_2, left, right]
 
     my_model.interfaces = [F.Interface(5, (vol_1, vol_2), penalty_term=100)]
-    my_model.surface_to_volume = {
-        left: vol_1,
-        right: vol_2,
-    }
 
     H = F.Species("H", mobile=True, subdomains=[vol_1, vol_2])
     my_model.species = [H]

test/system_tests/test_misc.py

@@ -210,11 +210,6 @@ def test_temp_dependent_bc_mixed_domain_temperature_as_function():
         bottom_surface,
     ]
 
-    my_model.surface_to_volume = {
-        top_surface: top_domain,
-        bottom_surface: bottom_domain,
-    }
-
     H = F.Species("H", mobile=True, subdomains=[bottom_domain, top_domain])
 
     my_model.species = [H]

test/system_tests/test_multi_mat_penalty.py

@@ -112,10 +112,6 @@ def c_exact_bot_np(x):
     my_model.interfaces = [
         F.Interface(5, (bottom_domain, top_domain), penalty_term=1),
     ]
-    my_model.surface_to_volume = {
-        top_surface: top_domain,
-        bottom_surface: bottom_domain,
-    }
 
     H = F.Species("H", mobile=True)
 
@@ -195,10 +191,6 @@ def test_derived_quantities_multi_mat():
     my_model.interfaces = [
         F.Interface(5, (bottom_domain, top_domain), penalty_term=1),
     ]
-    my_model.surface_to_volume = {
-        top_surface: top_domain,
-        bottom_surface: bottom_domain,
-    }
 
     H = F.Species("H", mobile=True)
 
@@ -270,11 +262,6 @@ def test_penalty_multispecies():
     for spe in my_model.species:
         spe.subdomains = [vol1, vol2]
 
-    my_model.surface_to_volume = {
-        left_surf: vol1,
-        right_surf: vol2,
-    }
-
     my_model.boundary_conditions = [
         # Protium BCs
         F.FixedConcentrationBC(left_surf, value=10, species=protium),

test/system_tests/test_multi_material.py

@@ -110,10 +110,6 @@ def c_exact_bot_np(x):
     ]
 
     my_model.interfaces = [F.Interface(5, (bottom_domain, top_domain))]
-    my_model.surface_to_volume = {
-        top_surface: top_domain,
-        bottom_surface: bottom_domain,
-    }
 
     H = F.Species("H", mobile=True)
 
@@ -182,10 +178,6 @@ def test_1_material_discontinuous_version(tmpdir):
     right_surface = F.SurfaceSubdomain1D(id=2, x=vertices[-1])
 
     my_model.subdomains = [subdomain, left_surface, right_surface]
-    my_model.surface_to_volume = {
-        right_surface: subdomain,
-        left_surface: subdomain,
-    }
 
     H = F.Species("H", mobile=True)
     trapped_H = F.Species("H_trapped", mobile=False)
@@ -278,10 +270,6 @@ def test_3_materials_transient(tmpdir):
         left_surface,
         right_surface,
     ]
-    my_model.surface_to_volume = {
-        right_surface: right_domain,
-        left_surface: left_domain,
-    }
 
     H = F.Species("H", mobile=True)
     trapped_H = F.Species("H_trapped", mobile=False)
@@ -367,10 +355,6 @@ def test_2_mats_particle_flux_bc(tmpdir):
 
     # we should be able to automate this
     my_model.interfaces = [F.Interface(5, (bottom_domain, top_domain))]
-    my_model.surface_to_volume = {
-        top_surface: top_domain,
-        bottom_surface: bottom_domain,
-    }
 
     H = F.Species("H", mobile=True)
 
@@ -418,10 +402,6 @@ def test_all_cells_are_not_tagged():
 
     my_model.subdomains = [vol1, vol2, surf1, surf2]
 
-    my_model.surface_to_volume = {
-        surf1: vol1,
-        surf2: vol2,
-    }
     my_model.interfaces = [F.Interface(id=3, subdomains=[vol1, vol2])]
 
     my_model.mesh = F.Mesh1D(np.linspace(0, 1, 10))

test/system_tests/test_spherical_coordinates.py

@@ -116,10 +116,6 @@ def c_exact_right(x):
     my_model.subdomains = [vol_1, vol_2, left, right]
 
     my_model.interfaces = [F.Interface(5, (vol_1, vol_2), penalty_term=1e5)]
-    my_model.surface_to_volume = {
-        left: vol_1,
-        right: vol_2,
-    }
 
     H = F.Species("H", mobile=True, subdomains=[vol_1, vol_2])
     my_model.species = [H]

test/test_complex_reaction.py

@@ -188,10 +188,6 @@ def test_surface_reaction_in_discontinuous_case():
     right_surface = F.SurfaceSubdomain1D(id=2, x=1)
 
     my_model.subdomains = [subdomain, left_surface, right_surface]
-    my_model.surface_to_volume = {
-        right_surface: subdomain,
-        left_surface: subdomain,
-    }
 
     H = F.Species("H", mobile=True, subdomains=[subdomain])
     my_model.species = [H]

test/test_h_transport_problem.py

@@ -1411,11 +1411,6 @@ def test_surface_reaction_BC_discontinuous():
         outlet,
     ]
 
-    my_model.surface_to_volume = {
-        inlet: left_vol,
-        outlet: right_vol,
-    }
-
     my_model.method_interface = F.InterfaceMethod.nitsche
     my_model.interfaces = [
         F.Interface(id=5, subdomains=[left_vol, right_vol], penalty_term=10),