Merge pull request #593 from jhdark/subdomains

Subdomain classes

festim/__init__.py

@@ -21,4 +21,7 @@
 
 from .species import Species, Trap
 
+from .subdomain.surface_subdomain import SurfaceSubdomain1D
+from .subdomain.volume_subdomain import VolumeSubdomain1D
+
 from .hydrogen_transport_problem import HydrogenTransportProblem

festim/hydrogen_transport_problem.py

@@ -3,12 +3,9 @@
 import ufl
 from mpi4py import MPI
 from dolfinx.fem import Function
-from ufl import (
-    TestFunction,
-    dot,
-    grad,
-    exp,
-)
+from dolfinx.mesh import meshtags, locate_entities
+from ufl import TestFunction, dot, grad, exp, Measure
+import numpy as np
 
 import festim as F
 
@@ -38,8 +35,8 @@ class HydrogenTransportProblem:
         dx (dolfinx.fem.dx): the volume measure of the model
         ds (dolfinx.fem.ds): the surface measure of the model
         function_space (dolfinx.fem.FunctionSpace): the function space of the model
-        facet_tags (dolfinx.cpp.mesh.MeshTags): the facet tags of the model
-        volume_tags (dolfinx.cpp.mesh.MeshTags): the volume tags of the model
+        facet_meshtags (dolfinx.cpp.mesh.MeshTags): the facet tags of the model
+        volume_meshtags (dolfinx.cpp.mesh.MeshTags): the volume tags of the model
         formulation (ufl.form.Form): the formulation of the model
         solver (dolfinx.nls.newton.NewtonSolver): the solver of the model
 
@@ -88,8 +85,8 @@ def __init__(
         self.dx = None
         self.ds = None
         self.function_space = None
-        self.facet_tags = None
-        self.volume_tags = None
+        self.facet_meshtags = None
+        self.volume_meshtags = None
         self.formulation = None
 
     def initialise(self):
@@ -98,19 +95,54 @@ def initialise(self):
         """
 
         self.define_function_space()
-        (
-            self.facet_tags,
-            self.volume_tags,
-            self.dx,
-            self.ds,
-        ) = self.mesh.create_measures_and_tags(self.function_space)
+        self.define_markers_and_measures()
         self.assign_functions_to_species()
         self.create_formulation()
 
     def define_function_space(self):
         elements = ufl.FiniteElement("CG", self.mesh.mesh.ufl_cell(), 1)
         self.function_space = fem.FunctionSpace(self.mesh.mesh, elements)
 
+    def define_markers_and_measures(self):
+        dofs_facets, tags_facets = [], []
+
+        # TODO this should be a property of mesh
+        fdim = self.mesh.mesh.topology.dim - 1
+        vdim = self.mesh.mesh.topology.dim
+
+        # find all cells in domain and mark them as 0
+        num_cells = self.mesh.mesh.topology.index_map(vdim).size_local
+        mesh_cell_indices = np.arange(num_cells, dtype=np.int32)
+        tags_volumes = np.full(num_cells, 0, dtype=np.int32)
+
+        for sub_dom in self.subdomains:
+            if isinstance(sub_dom, F.SurfaceSubdomain1D):
+                dof = sub_dom.locate_dof(self.function_space)
+                dofs_facets.append(dof)
+                tags_facets.append(sub_dom.id)
+            if isinstance(sub_dom, F.VolumeSubdomain1D):
+                # find all cells in subdomain and mark them as sub_dom.id
+                entities = sub_dom.locate_subdomain_entities(self.mesh.mesh, vdim)
+                tags_volumes[entities] = sub_dom.id
+
+        # 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)
+
+        # define mesh tags
+        self.facet_meshtags = meshtags(self.mesh.mesh, fdim, dofs_facets, tags_facets)
+        self.volume_meshtags = meshtags(
+            self.mesh.mesh, vdim, mesh_cell_indices, tags_volumes
+        )
+
+        # define measures
+        self.ds = Measure(
+            "ds", domain=self.mesh.mesh, subdomain_data=self.facet_meshtags
+        )
+        self.dx = Measure(
+            "dx", domain=self.mesh.mesh, subdomain_data=self.volume_meshtags
+        )
+
     def assign_functions_to_species(self):
         """Creates for each species the solution, prev solution and test function"""
         if len(self.species) > 1:
@@ -124,8 +156,6 @@ def create_formulation(self):
         """Creates the formulation of the model"""
         if len(self.sources) > 1:
             raise NotImplementedError("Sources not implemented yet")
-        if len(self.subdomains) > 1:
-            raise NotImplementedError("Multiple subdomains not implemented yet")
         if len(self.species) > 1:
             raise NotImplementedError("Multiple species not implemented yet")
 

festim/mesh/mesh.py

@@ -26,17 +26,3 @@ def __init__(self, mesh=None):
             self.mesh.topology.create_connectivity(
                 self.mesh.topology.dim - 1, self.mesh.topology.dim
             )
-
-    def create_measures_and_tags(self, function_space):
-        """Creates the ufl.measure.Measure objects for self.ds and
-        self.dx, also passes the facet and volume tags
-        """
-        facet_tags, volume_tags = self.create_meshtags(function_space)
-        dx = ufl.Measure("dx", domain=self.mesh, subdomain_data=volume_tags)
-        ds = ufl.Measure("ds", domain=self.mesh, subdomain_data=facet_tags)
-        return (
-            facet_tags,
-            volume_tags,
-            dx,
-            ds,
-        )

festim/mesh/mesh_1d.py

@@ -36,34 +36,3 @@ def generate_mesh(self):
         cells = np.stack((indexes[:-1], indexes[1:]), axis=-1)
 
         return mesh.create_mesh(MPI.COMM_WORLD, cells, mesh_points, domain)
-
-    def create_meshtags(self, function_space):
-        """Creates the meshtags for a given function space
-        Args:
-            function_space (dolfinx.fem.function.FunctionSpace): the function
-                space of the model
-
-        Returns:
-            dolfinx.mesh.MeshTagsMetaClass: the tags containing the facet
-                and volume tags
-        """
-
-        dofs_L = fem.locate_dofs_geometrical(
-            function_space, lambda x: np.isclose(x[0], min(self.vertices))
-        )
-        dofs_R = fem.locate_dofs_geometrical(
-            function_space, lambda x: np.isclose(x[0], max(self.vertices))
-        )
-
-        dofs_facets = np.array([dofs_L[0], dofs_R[0]], dtype=np.int32)
-        tags_facets = np.array([1, 2], dtype=np.int32)
-
-        facet_dimension = self.mesh.topology.dim - 1
-        mesh_tags_facets = mesh.meshtags(
-            self.mesh, facet_dimension, dofs_facets, tags_facets
-        )
-
-        # TODO implement this
-        mesh_tags_volumes = None
-
-        return mesh_tags_facets, mesh_tags_volumes

festim/subdomain/surface_subdomain.py

@@ -0,0 +1,39 @@
+from dolfinx import fem
+import numpy as np
+
+
+class SurfaceSubdomain1D:
+    """
+    Surface subdomain class for 1D cases
+
+    Args:
+        id (int): the id of the surface subdomain
+        x (float): the x coordinate of the surface subdomain
+
+    Attributes:
+        id (int): the id of the surface subdomain
+        x (float): the x coordinate of the surface subdomain
+
+    Usage:
+        >>> surf_subdomain = F.SurfaceSubdomain1D(id=1, x=1)
+    """
+
+    def __init__(self, id, x) -> None:
+        self.id = id
+        self.x = x
+
+    def locate_dof(self, function_space):
+        """Locates the dof of the surface subdomain within the function space
+
+        Args:
+            function_space (dolfinx.fem.FunctionSpace): the function space of
+                the model
+
+        Returns:
+            dof (np.array): the first value in the list of dofs of the surface
+                subdomain
+        """
+        dofs = fem.locate_dofs_geometrical(
+            function_space, lambda x: np.isclose(x[0], self.x)
+        )
+        return dofs[0]

festim/subdomain/volume_subdomain.py

@@ -0,0 +1,45 @@
+from dolfinx.mesh import locate_entities
+import numpy as np
+
+
+class VolumeSubdomain1D:
+    """
+    Volume subdomain class for 1D cases
+
+    Args:
+        id (int): the id of the volume subdomain
+        borders (list of float): the borders of the volume subdomain
+        material (festim.Material): the material of the volume subdomain
+
+    Attributes:
+        id (int): the id of the volume subdomain
+        borders (list of float): the borders of the volume subdomain
+        material (festim.Material): the material of the volume subdomain
+
+    Usage:
+        >>> vol_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, 1],
+        ...     material=F.Material(...))
+    """
+
+    def __init__(self, id, borders, material) -> None:
+        self.borders = borders
+        self.material = material
+        self.id = id
+
+    def locate_subdomain_entities(self, mesh, vdim):
+        """Locates all cells in subdomain borders within domain
+
+        Args:
+            mesh (dolfinx.cpp.mesh.Mesh): the mesh of the model
+            vdim (int): the dimension of the volumes of the mesh,
+                for 1D this is always 1
+
+        Returns:
+            entities (np.array): the entities of the subdomain
+        """
+        entities = locate_entities(
+            mesh,
+            vdim,
+            lambda x: np.logical_and(x[0] >= self.borders[0], x[0] <= self.borders[1]),
+        )
+        return entities

test/test_permeation_problem.py

@@ -26,6 +26,11 @@ def test_permeation_problem():
     my_model = F.HydrogenTransportProblem()
     my_model.mesh = my_mesh
 
+    my_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=None)
+    left_surface = F.SurfaceSubdomain1D(id=1, x=0)
+    right_surface = F.SurfaceSubdomain1D(id=2, x=L)
+    my_model.subdomains = [my_subdomain, left_surface, right_surface]
+
     mobile_H = F.Species("H")
     my_model.species = [mobile_H]
 
@@ -45,9 +50,9 @@ def siverts_law(T, S_0, E_S, pressure):
         return S * pressure**0.5
 
     fdim = my_mesh.mesh.topology.dim - 1
-    left_facets = my_model.facet_tags.find(1)
+    left_facets = my_model.facet_meshtags.find(1)
     left_dofs = locate_dofs_topological(V, fdim, left_facets)
-    right_facets = my_model.facet_tags.find(2)
+    right_facets = my_model.facet_meshtags.find(2)
     right_dofs = locate_dofs_topological(V, fdim, right_facets)
 
     S_0 = 4.02e21
@@ -116,10 +121,9 @@ def siverts_law(T, S_0, E_S, pressure):
     analytical_flux = P_up**0.5 * permeability / L * (2 * summation + 1)
 
     analytical_flux = np.abs(analytical_flux)
-
     flux_values = np.array(np.abs(flux_values))
 
-    relative_error = (flux_values - analytical_flux) / analytical_flux
+    relative_error = np.abs((flux_values - analytical_flux) / analytical_flux)
 
     relative_error = relative_error[
         np.where(analytical_flux > 0.01 * np.max(analytical_flux))

test/test_subdomains.py

@@ -0,0 +1,48 @@
+import numpy as np
+import festim as F
+
+
+def test_different_surface_ids():
+    my_model = F.HydrogenTransportProblem()
+
+    L = 3e-02
+    my_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]
+
+    my_model.define_function_space()
+    my_model.define_markers_and_measures()
+
+    for surf_id in surfacec_subdomains_ids:
+        assert surf_id in np.array(my_model.facet_meshtags.values)
+
+
+def test_different_volume_ids():
+    my_model = F.HydrogenTransportProblem()
+
+    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_subdomain_1 = F.VolumeSubdomain1D(
+        id=vol_subdomains_ids[0], borders=[sub_dom_1[0], sub_dom_1[-1]], material=None
+    )
+    vol_subdomain_2 = F.VolumeSubdomain1D(
+        id=vol_subdomains_ids[1], borders=[sub_dom_2[0], sub_dom_2[-1]], material=None
+    )
+    vol_subdomain_3 = F.VolumeSubdomain1D(
+        id=vol_subdomains_ids[2], borders=[sub_dom_3[0], sub_dom_3[-1]], 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:
+        assert vol_id in np.array(my_model.volume_meshtags.values)