From 61a28c49368dcbf33678890cbab0909b2f59c0bc Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 22 Apr 2026 14:15:35 -0400 Subject: [PATCH 01/14] added tests --- test/system_tests/test_multi_material.py | 49 ++++++++++++++++++++++++ test/test_subdomains.py | 20 ++++++++++ 2 files changed, 69 insertions(+) diff --git a/test/system_tests/test_multi_material.py b/test/system_tests/test_multi_material.py index 86526d7b4..3a7c1551a 100644 --- a/test/system_tests/test_multi_material.py +++ b/test/system_tests/test_multi_material.py @@ -3,6 +3,7 @@ import dolfinx import dolfinx.fem.petsc import numpy as np +import pytest import ufl import festim as F @@ -396,3 +397,51 @@ def test_2_mats_particle_flux_bc(tmpdir): my_model.initialise() my_model.run() + + +def test_all_cells_are_not_tagged(): + """ + Checks that an error is raised when not all cells are tagged with a non-zero value. + This can be caused by a volume subdomain not being defined correctly, or by a mesh + that is too coarse to capture the geometry of the volume subdomains. + """ + + my_model = F.HydrogenTransportProblemDiscontinuous() + + mat = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0) + + vol1 = F.VolumeSubdomain1D(id=1, borders=[0, 0.5], material=mat) + vol2 = F.VolumeSubdomain1D(id=2, borders=[0.5, 1], material=mat) + + surf1 = F.SurfaceSubdomain1D(id=1, x=0) + surf2 = F.SurfaceSubdomain1D(id=2, x=1) + + 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)) + + my_model.species = [F.Species("H", subdomains=[vol1, vol2])] + + my_model.boundary_conditions = [ + F.FixedConcentrationBC(species=my_model.species[0], subdomain=surf1, value=1), + ] + + my_model.temperature = 300 + + my_model.settings = F.Settings(transient=False, atol=1e-9, rtol=1e-9) + + my_model.exports = [ + F.AverageVolume(field=my_model.species[0], volume=vol1), + F.AverageVolume(field=my_model.species[0], volume=vol2), + ] + + with pytest.raises( + AssertionError, match="All cells must be tagged with a non-zero value" + ): + my_model.initialise() diff --git a/test/test_subdomains.py b/test/test_subdomains.py index 9f87ab145..757a9af01 100644 --- a/test/test_subdomains.py +++ b/test/test_subdomains.py @@ -214,3 +214,23 @@ def test_name_setter(): with pytest.raises(TypeError, match="Name must be a string"): F.VolumeSubdomain(id=1, material=None, name=1) + + +def test_all_cells_are_not_tagged(): + """ + Checks that an error is raised when not all cells are tagged with a non-zero value. + This can be caused by a volume subdomain not being defined correctly, or by a mesh + that is too coarse to capture the geometry of the volume subdomains. + """ + + mat = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0) + + vol1 = F.VolumeSubdomain1D(id=1, borders=[0, 0.5], material=mat) + vol2 = F.VolumeSubdomain1D(id=2, borders=[0.5, 1], material=mat) + + mesh = F.Mesh1D(np.linspace(0, 1, 10)) + + with pytest.raises( + AssertionError, match="All cells must be tagged with a non-zero value" + ): + mesh.define_meshtags(surface_subdomains=[], volume_subdomains=[vol1, vol2]) From cbf5064b34b17ec0e639e5e6a269e854d1266b2a Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 22 Apr 2026 14:16:22 -0400 Subject: [PATCH 02/14] added check in define_meshtags --- src/festim/mesh/mesh.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/festim/mesh/mesh.py b/src/festim/mesh/mesh.py index 032897c7a..d08d4cab6 100644 --- a/src/festim/mesh/mesh.py +++ b/src/festim/mesh/mesh.py @@ -180,6 +180,14 @@ def define_meshtags(self, surface_subdomains, volume_subdomains, interfaces=None ) tags_volumes[:] = vol.id + # make sure there are no zeros in the tags_volumes + assert np.all(tags_volumes != 0), ( + "All cells must be tagged with a non-zero value" + "This error can be caused by a volume subdomain not being defined" + "correctly, or by a mesh that is too coarse to capture the geometry" + "of the volume subdomains" + ) + volume_meshtags = meshtags( self._mesh, self.vdim, mesh_cell_indices, tags_volumes ) From 7e64f6375e144a2dfbf674ab83bbdcd5a06838e1 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 22 Apr 2026 14:17:11 -0400 Subject: [PATCH 03/14] added an assertion to avoid id=0 --- src/festim/subdomain/volume_subdomain.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/festim/subdomain/volume_subdomain.py b/src/festim/subdomain/volume_subdomain.py index 56bdbd806..22a735b00 100644 --- a/src/festim/subdomain/volume_subdomain.py +++ b/src/festim/subdomain/volume_subdomain.py @@ -24,7 +24,7 @@ class VolumeSubdomain: Volume subdomain class Args: - id: the id of the volume subdomain + id: the id of the volume subdomain (> 0) submesh: the submesh of the volume subdomain cell_map: the cell map of the volume subdomain parent_mesh: the parent mesh of the volume subdomain @@ -58,6 +58,7 @@ class VolumeSubdomain: def __init__( self, id, material, locator: Callable | None = None, name: str | None = None ): + assert id != 0, "Volume subdomain id cannot be 0" self.id = id self.material = material self.locator = locator From 6632f995b5dc1ec4a7561fe51b70c98b5cd5934a Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 22 Apr 2026 14:21:17 -0400 Subject: [PATCH 04/14] fixed test --- test/system_tests/test_io.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/system_tests/test_io.py b/test/system_tests/test_io.py index ffae401a8..f0eceb58a 100644 --- a/test/system_tests/test_io.py +++ b/test/system_tests/test_io.py @@ -19,7 +19,7 @@ def test_writing_and_reading_of_species_function_using_checkpoints(tmpdir): my_model.mesh = F.Mesh(mesh) my_mat = F.Material(name="mat", D_0=1, E_D=0) - vol = F.VolumeSubdomain(id=0, material=my_mat) + vol = F.VolumeSubdomain(id=1, material=my_mat) surf = F.SurfaceSubdomain(id=1) my_model.subdomains = [vol, surf] From 649e10bd5311985fc7b492b7818aab6544774d60 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 22 Apr 2026 14:24:13 -0400 Subject: [PATCH 05/14] fixed more tests --- test/test_heat_transfer_problem.py | 2 +- test/test_initial_condition.py | 4 ++-- test/test_subdomains.py | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/test/test_heat_transfer_problem.py b/test/test_heat_transfer_problem.py index 39210a320..c4e3dd5f1 100644 --- a/test/test_heat_transfer_problem.py +++ b/test/test_heat_transfer_problem.py @@ -484,7 +484,7 @@ def test_meshtags_from_xdmf(tmp_path, mesh): def test_raise_error_non_unique_vol_ids(): """Test that an error is raised if the volume ids are not unique""" my_problem = F.HeatTransferProblem() - my_problem.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 2000)) + my_problem.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 11)) left = F.SurfaceSubdomain1D(id=1, x=0) right = F.SurfaceSubdomain1D(id=2, x=1) mat = F.Material(D_0=1, E_D=None) diff --git a/test/test_initial_condition.py b/test/test_initial_condition.py index 6ced0b703..588785445 100644 --- a/test/test_initial_condition.py +++ b/test/test_initial_condition.py @@ -169,7 +169,7 @@ def f(x): my_problem.mesh = F.Mesh(mesh) mat = F.Material(D_0=1, E_D=0.1, name="dummy_mat") - vol = F.VolumeSubdomain(id=0, material=mat) + vol = F.VolumeSubdomain(id=1, material=mat) my_problem.subdomains = [vol] function_initial_value = F.read_function_from_file( @@ -227,7 +227,7 @@ def f2(x): my_problem.mesh = F.Mesh(mesh) mat = F.Material(D_0=1, E_D=0.1, name="dummy_mat") - vol = F.VolumeSubdomain(id=0, material=mat) + vol = F.VolumeSubdomain(id=1, material=mat) my_problem.subdomains = [vol] my_problem.initial_conditions = [ diff --git a/test/test_subdomains.py b/test/test_subdomains.py index 757a9af01..351ac0e42 100644 --- a/test/test_subdomains.py +++ b/test/test_subdomains.py @@ -125,7 +125,7 @@ def test_ValueError_raised_when_id_not_found_in_volumes_subdomains(): def test_ValueError_rasied_when_volume_ids_are_not_unique(): """Checks""" my_test_model = F.HydrogenTransportProblem( - mesh=F.Mesh1D(np.linspace(0, 2, num=10)), species=[F.Species("H")] + mesh=F.Mesh1D(np.linspace(0, 2, num=11)), species=[F.Species("H")] ) vol_1 = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=None) From 630e28106de0dd59131cd6e6ed2ef4b7d88276a7 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 22 Apr 2026 14:25:07 -0400 Subject: [PATCH 06/14] ruff fixes --- src/festim/hydrogen_transport_problem.py | 1 - src/festim/problem.py | 2 +- test/test_reaction.py | 2 +- 3 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 14d9c35b4..58e196c62 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -41,7 +41,6 @@ is_it_time_to_export, nmm_interpolate, ) -from festim.material import SolubilityLaw from .mesh import CoordinateSystem, Mesh diff --git a/src/festim/problem.py b/src/festim/problem.py index 15901c797..5694d23f0 100644 --- a/src/festim/problem.py +++ b/src/festim/problem.py @@ -1,3 +1,4 @@ +import warnings from typing import Any from mpi4py import MPI @@ -10,7 +11,6 @@ from dolfinx import fem from dolfinx.nls.petsc import NewtonSolver from packaging.version import Version -import warnings import festim as F from festim.mesh.mesh import Mesh as _Mesh diff --git a/test/test_reaction.py b/test/test_reaction.py index a34861763..0dd24e331 100644 --- a/test/test_reaction.py +++ b/test/test_reaction.py @@ -1,10 +1,10 @@ from mpi4py import MPI +import numpy as np import pytest from dolfinx.fem import Function, functionspace from dolfinx.mesh import create_unit_cube from ufl import exp -import numpy as np import festim as F From cb2026bb30918833f4923bca1809559327698d03 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 22 Apr 2026 14:44:58 -0400 Subject: [PATCH 07/14] fixed buidl --- test/test_reaction.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_reaction.py b/test/test_reaction.py index 0dd24e331..397d82faa 100644 --- a/test/test_reaction.py +++ b/test/test_reaction.py @@ -401,7 +401,7 @@ def test_product_setter_raise_error_E_p_no_product(): vol1 = F.VolumeSubdomain1D(id=1, borders=[0, 0.5], material=mat) vol2 = F.VolumeSubdomain1D(id=2, borders=[0.5, 1], material=mat) my_model = F.HydrogenTransportProblemDiscontinuous() -my_model.mesh = F.Mesh1D(np.linspace(0, 1, 100)) +my_model.mesh = F.Mesh1D(np.linspace(0, 1, 101)) my_model.subdomains = [vol1, vol2] From 0480c8bbb2a9c3f7fd79c9857df7ee03fc57f03b Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 23 Apr 2026 10:03:00 -0400 Subject: [PATCH 08/14] added test --- test/test_subdomains.py | 46 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/test/test_subdomains.py b/test/test_subdomains.py index 351ac0e42..844b0fd38 100644 --- a/test/test_subdomains.py +++ b/test/test_subdomains.py @@ -234,3 +234,49 @@ def test_all_cells_are_not_tagged(): AssertionError, match="All cells must be tagged with a non-zero value" ): mesh.define_meshtags(surface_subdomains=[], volume_subdomains=[vol1, vol2]) + + +def test_map_surface_to_volume_subdomains(): + """ + Tests that the function map_surface_to_volume_subdomains correctly maps surface + subdomains to volume subdomains based on the facet to cell connectivity + """ + + n = 20 + mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, n, n) + + festim_mesh = F.Mesh(mesh) + + surface_1 = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[0], 0)) + surface_2 = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[0], 1)) + + material_1 = F.Material(D_0=1, E_D=0, name="material_1") + volume_1 = F.VolumeSubdomain( + id=2, material=material_1, locator=lambda x: x[0] <= 0.5 + ) + volume_2 = F.VolumeSubdomain( + id=3, material=material_1, locator=lambda x: x[0] >= 0.5 + ) + + ft, ct = festim_mesh.define_meshtags( + surface_subdomains=[surface_1, surface_2], + volume_subdomains=[volume_1, volume_2], + ) + + facet_to_cell = mesh.topology.connectivity(mesh.topology.dim - 1, mesh.topology.dim) + + surface_to_subdomain = F.map_surface_to_volume_subdomains( + ft, ct, facet_to_cell, [volume_1, volume_2], [surface_1, surface_2] + ) + print(surface_to_subdomain) + for surface, volume in surface_to_subdomain.items(): + print(f"Surface {surface.id} is connected to volume {volume.id}") + + assert surface_to_subdomain[surface_1] == volume_1, ( + "Expected surface 1 to be connected to volume 1, " + f"got {surface_to_subdomain[surface_1].id}" + ) + assert surface_to_subdomain[surface_2] == volume_2, ( + "Expected surface 2 to be connected to volume 2, " + f"got {surface_to_subdomain[surface_2].id}" + ) From ad0c94ee69e1164e05c618d7b6bc74da55979429 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 23 Apr 2026 10:03:24 -0400 Subject: [PATCH 09/14] added function --- src/festim/__init__.py | 1 + src/festim/subdomain/volume_subdomain.py | 59 ++++++++++++++++++++++++ 2 files changed, 60 insertions(+) diff --git a/src/festim/__init__.py b/src/festim/__init__.py index 823d519d9..b6133ed86 100644 --- a/src/festim/__init__.py +++ b/src/festim/__init__.py @@ -79,5 +79,6 @@ VolumeSubdomain, VolumeSubdomain1D, find_volume_from_id, + map_surface_to_volume_subdomains, ) from .trap import Trap diff --git a/src/festim/subdomain/volume_subdomain.py b/src/festim/subdomain/volume_subdomain.py index 22a735b00..29bcbc2b1 100644 --- a/src/festim/subdomain/volume_subdomain.py +++ b/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,61 @@ 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], +): + + # 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) + + 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: + print(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 From a681a6ae3946de29384798541b9a6baeb6cc7641 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 23 Apr 2026 10:09:50 -0400 Subject: [PATCH 10/14] documentation --- src/festim/subdomain/volume_subdomain.py | 25 ++++++++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/src/festim/subdomain/volume_subdomain.py b/src/festim/subdomain/volume_subdomain.py index 29bcbc2b1..e111f4b70 100644 --- a/src/festim/subdomain/volume_subdomain.py +++ b/src/festim/subdomain/volume_subdomain.py @@ -194,7 +194,25 @@ def map_surface_to_volume_subdomains( facet_to_cell: dolfinx.cpp.graph.AdjacencyList_int32, volume_subdomains: list[VolumeSubdomain], surface_subdomains: list[SurfaceSubdomain], -): +) -> 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 + + 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] @@ -231,7 +249,10 @@ def map_surface_to_volume_subdomains( surface_to_subdomain = {} for s_tag, v_tag in unique_pairs: - print(f"Facet tag {s_tag} is connected to cell tag {v_tag}") + dolfinx.log.log( + dolfinx.log.LogLevel.DEBUG, + 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) From d3d94113067806f176bd4c7012cc820f918d0a77 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 23 Apr 2026 10:26:39 -0400 Subject: [PATCH 11/14] surface_to_volume no longer explicit --- .../advection_diffusion_multi_material.py | 5 ---- examples/multi_material_1d.py | 4 ---- examples/multi_material_2d.py | 4 ---- examples/multi_material_2d_bis.py | 4 ---- examples/multi_material_transient.py | 4 ---- examples/multi_material_with_one_material.py | 1 - src/festim/hydrogen_transport_problem.py | 23 +++++++++++++++++-- src/festim/subdomain/__init__.py | 7 +++++- .../test_advection_term_integration.py | 5 ---- .../test_cylindrical_coordinates.py | 4 ---- test/system_tests/test_misc.py | 5 ---- test/system_tests/test_multi_mat_penalty.py | 13 ----------- test/system_tests/test_multi_material.py | 20 ---------------- .../test_spherical_coordinates.py | 4 ---- test/test_complex_reaction.py | 4 ---- test/test_h_transport_problem.py | 5 ---- test/test_initial_condition.py | 4 ---- test/test_profile.py | 10 -------- test/test_vtx.py | 1 - 19 files changed, 27 insertions(+), 100 deletions(-) diff --git a/examples/advection_diffusion_multi_material.py b/examples/advection_diffusion_multi_material.py index a1869335a..b28d53a17 100644 --- a/examples/advection_diffusion_multi_material.py +++ b/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) diff --git a/examples/multi_material_1d.py b/examples/multi_material_1d.py index 2fb2a815c..6fbd9ec31 100644 --- a/examples/multi_material_1d.py +++ b/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) diff --git a/examples/multi_material_2d.py b/examples/multi_material_2d.py index 7a64cbe8c..eb7e49e38 100644 --- a/examples/multi_material_2d.py +++ b/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) diff --git a/examples/multi_material_2d_bis.py b/examples/multi_material_2d_bis.py index cc6ad33da..384052a94 100644 --- a/examples/multi_material_2d_bis.py +++ b/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) diff --git a/examples/multi_material_transient.py b/examples/multi_material_transient.py index 8495336fa..a98c5e274 100644 --- a/examples/multi_material_transient.py +++ b/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) diff --git a/examples/multi_material_with_one_material.py b/examples/multi_material_with_one_material.py index ca8c80422..d58f3ffc2 100644 --- a/examples/multi_material_with_one_material.py +++ b/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) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 58e196c62..29349bc8c 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1054,7 +1054,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 @@ -1092,7 +1091,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 @@ -1140,6 +1139,26 @@ 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, + ) # create submeshes and transfer meshtags to subdomains for subdomain in self.volume_subdomains: diff --git a/src/festim/subdomain/__init__.py b/src/festim/subdomain/__init__.py index d3fb7ec57..7679593f2 100644 --- a/src/festim/subdomain/__init__.py +++ b/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", ] diff --git a/test/system_tests/test_advection_term_integration.py b/test/system_tests/test_advection_term_integration.py index 3d3e728a0..3a11b1fea 100644 --- a/test/system_tests/test_advection_term_integration.py +++ b/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) diff --git a/test/system_tests/test_cylindrical_coordinates.py b/test/system_tests/test_cylindrical_coordinates.py index eaa313e98..460fa1249 100644 --- a/test/system_tests/test_cylindrical_coordinates.py +++ b/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] diff --git a/test/system_tests/test_misc.py b/test/system_tests/test_misc.py index 47d855942..5b9152017 100644 --- a/test/system_tests/test_misc.py +++ b/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] diff --git a/test/system_tests/test_multi_mat_penalty.py b/test/system_tests/test_multi_mat_penalty.py index 116da435e..f5ef5c132 100644 --- a/test/system_tests/test_multi_mat_penalty.py +++ b/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), diff --git a/test/system_tests/test_multi_material.py b/test/system_tests/test_multi_material.py index 3a7c1551a..41a31ee19 100644 --- a/test/system_tests/test_multi_material.py +++ b/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)) diff --git a/test/system_tests/test_spherical_coordinates.py b/test/system_tests/test_spherical_coordinates.py index d477f484d..2b2887405 100644 --- a/test/system_tests/test_spherical_coordinates.py +++ b/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] diff --git a/test/test_complex_reaction.py b/test/test_complex_reaction.py index bcbbb7923..e6fbabf18 100644 --- a/test/test_complex_reaction.py +++ b/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] diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index d33e002a2..bc41336de 100644 --- a/test/test_h_transport_problem.py +++ b/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), diff --git a/test/test_initial_condition.py b/test/test_initial_condition.py index 588785445..6214fe278 100644 --- a/test/test_initial_condition.py +++ b/test/test_initial_condition.py @@ -347,10 +347,6 @@ def test_initial_condition_discontinuous(): my_model.interfaces = [ F.Interface(id=3, subdomains=[vol1, vol2], penalty_term=1000) ] - my_model.surface_to_volume = { - left_surf: vol1, - right_surf: vol2, - } my_model.temperature = 300 diff --git a/test/test_profile.py b/test/test_profile.py index c52de4544..bd36291fd 100644 --- a/test/test_profile.py +++ b/test/test_profile.py @@ -144,11 +144,6 @@ def test_profile_discontinuous(): 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), @@ -215,11 +210,6 @@ def test_profile_discontinuous_single_species(): 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 = [ F.FixedConcentrationBC(left_surf, value=10, species=protium), F.FixedConcentrationBC(right_surf, value=0, species=protium), diff --git a/test/test_vtx.py b/test/test_vtx.py index edbb2e531..33ce7bf69 100644 --- a/test/test_vtx.py +++ b/test/test_vtx.py @@ -68,7 +68,6 @@ def test_vtx_DG(tmpdir): my_model.temperature = 55 my_model.subdomains = [s0, s1, l0, l1] - my_model.surface_to_volume = {l0: s0, l1: s1} # NOTE: Ask Remi why `H` has to live in both s0 and s1 my_model.species = [ F.Species("H", subdomains=[s0, s1]), From a6306eb577f51009d26ac63223ee6aa000e3b72c Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 23 Apr 2026 10:48:33 -0400 Subject: [PATCH 12/14] INFO instead of DEBUG --- src/festim/subdomain/volume_subdomain.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/festim/subdomain/volume_subdomain.py b/src/festim/subdomain/volume_subdomain.py index e111f4b70..5b018c4bc 100644 --- a/src/festim/subdomain/volume_subdomain.py +++ b/src/festim/subdomain/volume_subdomain.py @@ -250,7 +250,7 @@ def map_surface_to_volume_subdomains( for s_tag, v_tag in unique_pairs: dolfinx.log.log( - dolfinx.log.LogLevel.DEBUG, + 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) From a9503b1e6708fda5e03087dda9e04c26f6762c27 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Mon, 27 Apr 2026 19:23:20 -0400 Subject: [PATCH 13/14] parallel safe Co-authored-by: James Dark <65899899+jhdark@users.noreply.github.com> --- src/festim/hydrogen_transport_problem.py | 1 + src/festim/subdomain/volume_subdomain.py | 7 +++++++ 2 files changed, 8 insertions(+) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 29349bc8c..f172b0ed0 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1158,6 +1158,7 @@ def initialise(self): 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 diff --git a/src/festim/subdomain/volume_subdomain.py b/src/festim/subdomain/volume_subdomain.py index 5b018c4bc..f3d6c1a96 100644 --- a/src/festim/subdomain/volume_subdomain.py +++ b/src/festim/subdomain/volume_subdomain.py @@ -194,6 +194,7 @@ def map_surface_to_volume_subdomains( 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. @@ -208,6 +209,7 @@ def map_surface_to_volume_subdomains( 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 @@ -242,6 +244,11 @@ def map_surface_to_volume_subdomains( 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} From 333a554fba7b51503d27bfa3b68f8d2b4f26890c Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Mon, 27 Apr 2026 19:25:24 -0400 Subject: [PATCH 14/14] fix indent --- src/festim/subdomain/volume_subdomain.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/festim/subdomain/volume_subdomain.py b/src/festim/subdomain/volume_subdomain.py index f3d6c1a96..5fcc40269 100644 --- a/src/festim/subdomain/volume_subdomain.py +++ b/src/festim/subdomain/volume_subdomain.py @@ -244,7 +244,7 @@ def map_surface_to_volume_subdomains( 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: + 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: