Skip to content

Qol auto detect surface to volume #1122

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
RemDelaporteMathurin merged 15 commits into from
Apr 28, 2026
Merged

Qol auto detect surface to volume #1122

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
61a28c4
added tests
RemDelaporteMathurin Apr 22, 2026
cbf5064
added check in define_meshtags
RemDelaporteMathurin Apr 22, 2026
7e64f63
added an assertion to avoid id=0
RemDelaporteMathurin Apr 22, 2026
6632f99
fixed test
RemDelaporteMathurin Apr 22, 2026
649e10b
fixed more tests
RemDelaporteMathurin Apr 22, 2026
630e281
ruff fixes
RemDelaporteMathurin Apr 22, 2026
cb2026b
fixed buidl
RemDelaporteMathurin Apr 22, 2026
0480c8b
added test
RemDelaporteMathurin Apr 23, 2026
ad0c94e
added function
RemDelaporteMathurin Apr 23, 2026
a681a6a
documentation
RemDelaporteMathurin Apr 23, 2026
d3d9411
surface_to_volume no longer explicit
RemDelaporteMathurin Apr 23, 2026
a6306eb
INFO instead of DEBUG
RemDelaporteMathurin Apr 23, 2026
a9503b1
parallel safe
RemDelaporteMathurin Apr 27, 2026
333a554
fix indent
RemDelaporteMathurin Apr 27, 2026
fd21218
Merge remote-tracking branch 'origin/main' into qol-auto-detect-subdo…
RemDelaporteMathurin Apr 27, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 0 additions & 5 deletions examples/advection_diffusion_multi_material.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
4 changes: 0 additions & 4 deletions examples/multi_material_1d.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
4 changes: 0 additions & 4 deletions examples/multi_material_2d.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
4 changes: 0 additions & 4 deletions examples/multi_material_2d_bis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
4 changes: 0 additions & 4 deletions examples/multi_material_transient.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
1 change: 0 additions & 1 deletion examples/multi_material_with_one_material.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
1 change: 1 addition & 0 deletions src/festim/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -79,5 +79,6 @@
VolumeSubdomain,
VolumeSubdomain1D,
find_volume_from_id,
map_surface_to_volume_subdomains,
)
from .trap import Trap
24 changes: 22 additions & 2 deletions src/festim/hydrogen_transport_problem.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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 = {}
Copy link
Copy Markdown
Collaborator

@jhdark jhdark Apr 27, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should we add some kind of "deprecation" warning only for the next release, just noting that this is no longer required to be given by the user if a non-empty dict is given? Even though its not a deprecation

Copy link
Copy Markdown
Collaborator Author

@RemDelaporteMathurin RemDelaporteMathurin Apr 27, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes there is a a DeprecationWarning in HydrogenTransportProblemDiscontinuous.initialise()

https://github.com/festim-dev/FESTIM/pull/1122/changes#diff-586fa3ee395f75eb89f7d327f2f051a08b9d4f0cb520918f7e112d3ac7271af8R1142-R1150

self.subdomain_to_species = {} # maps subdomain to species defined in it

@property
Expand Down Expand Up @@ -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,
Comment thread
RemDelaporteMathurin marked this conversation as resolved.
comm=self.mesh.mesh.comm,
)

# create submeshes and transfer meshtags to subdomains
for subdomain in self.volume_subdomains:
Expand Down
7 changes: 6 additions & 1 deletion src/festim/subdomain/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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",
Expand All @@ -9,4 +13,5 @@
"SurfaceSubdomain1D",
"VolumeSubdomain",
"VolumeSubdomain1D",
"map_surface_to_volume_subdomains",
]
87 changes: 87 additions & 0 deletions src/festim/subdomain/volume_subdomain.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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]:
Comment thread
RemDelaporteMathurin marked this conversation as resolved.
"""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)

Comment thread
RemDelaporteMathurin marked this conversation as resolved.
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)
Comment thread
RemDelaporteMathurin marked this conversation as resolved.
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)

Comment thread
RemDelaporteMathurin marked this conversation as resolved.
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
5 changes: 0 additions & 5 deletions test/system_tests/test_advection_term_integration.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
4 changes: 0 additions & 4 deletions test/system_tests/test_cylindrical_coordinates.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down
5 changes: 0 additions & 5 deletions test/system_tests/test_misc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down
13 changes: 0 additions & 13 deletions test/system_tests/test_multi_mat_penalty.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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),
Expand Down
20 changes: 0 additions & 20 deletions test/system_tests/test_multi_material.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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))
Expand Down
4 changes: 0 additions & 4 deletions test/system_tests/test_spherical_coordinates.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down
4 changes: 0 additions & 4 deletions test/test_complex_reaction.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down
5 changes: 0 additions & 5 deletions test/test_h_transport_problem.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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),
Expand Down
Loading
Loading