Merge pull request #7 from jhdark/fenicsx

Fenicsx

festim/hydrogen_transport_problem.py

@@ -138,6 +138,10 @@ def define_markers_and_measures(self):
                 )
                 tags_volumes[entities] = sub_dom.id
 
+        # check if all borders are defined
+        if isinstance(self.mesh, F.Mesh1D):
+            self.mesh.check_borders(self.volume_subdomains)
+
         # dofs and tags need to be in np.in32 format for meshtags
         dofs_facets = np.array(dofs_facets, dtype=np.int32)
         tags_facets = np.array(tags_facets, dtype=np.int32)

festim/mesh/mesh_1d.py

@@ -1,4 +1,4 @@
-from dolfinx import mesh
+import dolfinx.mesh
 from mpi4py import MPI
 import ufl
 import numpy as np
@@ -32,4 +32,32 @@ def generate_mesh(self):
         indexes = np.arange(self.vertices.shape[0])
         cells = np.stack((indexes[:-1], indexes[1:]), axis=-1)
 
-        return mesh.create_mesh(MPI.COMM_WORLD, cells, mesh_points, domain)
+        return dolfinx.mesh.create_mesh(MPI.COMM_WORLD, cells, mesh_points, domain)
+
+    def check_borders(self, volume_subdomains):
+        """Checks that the borders of the subdomain are within the domain
+
+        Args:
+            volume_subdomains (list of festim.VolumeSubdomain1D): the volume subdomains
+
+        Raises:
+            Value error: if borders outside the domain
+        """
+        # check that subdomains are connected
+        all_borders = [border for vol in volume_subdomains for border in vol.borders]
+        sorted_borders = np.sort(all_borders)
+        for start, end in zip(sorted_borders[1:-2:2], sorted_borders[2:-1:2]):
+            if start != end:
+                raise ValueError("Subdomain borders don't match to each other")
+
+        # check volume subdomain is defined
+        # TODO this possible by default
+        if len(all_borders) == 0:
+            raise ValueError("No volume subdomains defined")
+
+        # check that subdomains are within the domain
+        if (
+            sorted_borders[0] != self.vertices[0]
+            or sorted_borders[-1] != self.vertices[-1]
+        ):
+            raise ValueError("borders dont match domain borders")

test/test_subdomains.py

@@ -1,48 +1,94 @@
 import numpy as np
 import festim as F
+import pytest
 
 
 def test_different_surface_ids():
-    my_model = F.HydrogenTransportProblem()
-
-    L = 3e-02
-    my_model.mesh = F.Mesh1D(np.linspace(0, L, num=3))
+    """Checks that different surface ids are correctly set"""
+    my_test_model = F.HydrogenTransportProblem()
+    L = 1e-04
+    my_test_model.mesh = F.Mesh1D(np.linspace(0, L, num=3))
 
-    surfacec_subdomains_ids = [3, 8]
-    surface_subdomain_1 = F.SurfaceSubdomain1D(id=surfacec_subdomains_ids[0], x=0)
-    surface_subdomain_2 = F.SurfaceSubdomain1D(id=surfacec_subdomains_ids[1], x=L)
-    my_model.subdomains = [surface_subdomain_1, surface_subdomain_2]
+    surface_subdomains_ids = [3, 8]
+    surface_subdomain_1 = F.SurfaceSubdomain1D(id=surface_subdomains_ids[0], x=0)
+    surface_subdomain_2 = F.SurfaceSubdomain1D(id=surface_subdomains_ids[1], x=L)
+    volume_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=None)
+    my_test_model.subdomains = [
+        surface_subdomain_1,
+        surface_subdomain_2,
+        volume_subdomain,
+    ]
 
-    my_model.define_function_space()
-    my_model.define_markers_and_measures()
+    my_test_model.define_function_space()
+    my_test_model.define_markers_and_measures()
 
-    for surf_id in surfacec_subdomains_ids:
-        assert surf_id in np.array(my_model.facet_meshtags.values)
+    for surf_id in surface_subdomains_ids:
+        assert surf_id in np.array(my_test_model.facet_meshtags.values)
 
 
 def test_different_volume_ids():
     my_model = F.HydrogenTransportProblem()
+    my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 10, 11))
 
-    sub_dom_1 = np.linspace(0, 1e-04, num=3)
-    sub_dom_2 = np.linspace(1e-04, 2e-04, num=4)
-    sub_dom_3 = np.linspace(2e-04, 3e-04, num=5)
-    my_model.mesh = F.Mesh1D(
-        np.unique(np.concatenate([sub_dom_1, sub_dom_2, sub_dom_3]))
-    )
-
-    vol_subdomains_ids = [2, 16, 7]
+    vol_subdom_ids = [2, 16, 7]
     vol_subdomain_1 = F.VolumeSubdomain1D(
-        id=vol_subdomains_ids[0], borders=[sub_dom_1[0], sub_dom_1[-1]], material=None
+        id=vol_subdom_ids[0], borders=[0, 2], material=None
     )
     vol_subdomain_2 = F.VolumeSubdomain1D(
-        id=vol_subdomains_ids[1], borders=[sub_dom_2[0], sub_dom_2[-1]], material=None
+        id=vol_subdom_ids[1], borders=[2, 6], material=None
     )
     vol_subdomain_3 = F.VolumeSubdomain1D(
-        id=vol_subdomains_ids[2], borders=[sub_dom_3[0], sub_dom_3[-1]], material=None
+        id=vol_subdom_ids[2], borders=[6, 10], material=None
     )
     my_model.subdomains = [vol_subdomain_1, vol_subdomain_2, vol_subdomain_3]
 
     my_model.define_markers_and_measures()
 
-    for vol_id in vol_subdomains_ids:
+    for vol_id in vol_subdom_ids:
         assert vol_id in np.array(my_model.volume_meshtags.values)
+
+
+def test_non_matching_volume_borders():
+    """Checks that non-matching borders raise an error"""
+    mesh = F.Mesh1D(vertices=np.linspace(0, 5, 6))
+    vol_subdomain_1 = F.VolumeSubdomain1D(id=1, borders=[0, 2], material=None)
+    vol_subdomain_2 = F.VolumeSubdomain1D(id=1, borders=[3, 5], material=None)
+    subdomains = [vol_subdomain_1, vol_subdomain_2]
+
+    with pytest.raises(ValueError, match="Subdomain borders don't match to each other"):
+        mesh.check_borders(subdomains)
+
+
+def test_matching_volume_borders_non_ascending_order():
+    """Checks that subdomain placed in non ascending order still passes"""
+    mesh = F.Mesh1D(vertices=np.linspace(0, 8, 9))
+    vol_subdomain_1 = F.VolumeSubdomain1D(id=1, borders=[0, 2], material=None)
+    vol_subdomain_2 = F.VolumeSubdomain1D(id=2, borders=[4, 8], material=None)
+    vol_subdomain_3 = F.VolumeSubdomain1D(id=3, borders=[2, 4], material=None)
+    subdomains = [vol_subdomain_1, vol_subdomain_2, vol_subdomain_3]
+
+    mesh.check_borders(subdomains)
+
+
+def test_borders_out_of_domain():
+    """Checks that borders outside of the domain raise an error"""
+    mesh = F.Mesh1D(vertices=np.linspace(0, 2))
+    subdomains = [F.VolumeSubdomain1D(id=1, borders=[1, 15], material=None)]
+    with pytest.raises(ValueError, match="borders dont match domain borders"):
+        mesh.check_borders(subdomains)
+
+
+def test_borders_inside_domain():
+    """Checks that borders inside of the domain raise an error"""
+    mesh = F.Mesh1D(vertices=np.linspace(0, 20))
+    subdomains = [F.VolumeSubdomain1D(id=1, borders=[1, 6], material=None)]
+    with pytest.raises(ValueError, match="borders dont match domain borders"):
+        mesh.check_borders(subdomains)
+
+
+def test_raise_error_with_no_volume_subdomain():
+    """Checks that error is rasied when no volume subdomain is defined"""
+    mesh = F.Mesh1D(vertices=np.linspace(0, 20))
+
+    with pytest.raises(ValueError, match="No volume subdomains defined"):
+        mesh.check_borders([])