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/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 ) 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/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 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] 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_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_reaction.py b/test/test_reaction.py index a34861763..397d82faa 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 @@ -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] diff --git a/test/test_subdomains.py b/test/test_subdomains.py index 9f87ab145..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) @@ -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])