Qol unmarked cells

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
 

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
         )

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

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

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]
 

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()

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)

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 = [

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]
 

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])