Merge pull request festim-dev#755 from RemDelaporteMathurin/custom-ta…

…gging

Simplify custom tagging of meshes

festim/exports/surface_quantity.py

@@ -12,7 +12,7 @@ class SurfaceQuantity:
 
     Attributes:
         field (festim.Species): species for which the surface flux is computed
-        surface (festim.SurfaceSubdomain1D): surface subdomain
+        surface (festim.SurfaceSubdomain): surface subdomain
         filename (str): name of the file to which the surface flux is exported
         t (list): list of time values
         data (list): list of values of the surface quantity
@@ -47,10 +47,8 @@ def surface(self):
 
     @surface.setter
     def surface(self, value):
-        if not isinstance(value, (int, F.SurfaceSubdomain1D)) or isinstance(
-            value, bool
-        ):
-            raise TypeError("surface should be an int or F.SurfaceSubdomain1D")
+        if not isinstance(value, (int, F.SurfaceSubdomain)) or isinstance(value, bool):
+            raise TypeError("surface should be an int or F.SurfaceSubdomain")
         self._surface = value
 
     @property

festim/heat_transfer_problem.py

@@ -112,26 +112,15 @@ def define_meshtags_and_measures(self):
             self.facet_meshtags = self.mesh.define_surface_meshtags()
             self.volume_meshtags = self.mesh.define_volume_meshtags()
 
-        elif isinstance(self.mesh, F.Mesh1D):
+        elif (
+            isinstance(self.mesh, F.Mesh)
+            and self.facet_meshtags is None
+            and self.volume_meshtags is None
+        ):
             self.facet_meshtags, self.volume_meshtags = self.mesh.define_meshtags(
                 surface_subdomains=self.surface_subdomains,
                 volume_subdomains=self.volume_subdomains,
             )
-        elif isinstance(self.mesh, F.Mesh):
-            # FIXME # refer to issue #647
-            # create empty mesh tags for now
-            facet_indices = np.array([], dtype=np.int32)
-            facet_tags = np.array([], dtype=np.int32)
-            self.facet_meshtags = meshtags(
-                self.mesh.mesh, self.mesh.fdim, facet_indices, facet_tags
-            )
-
-            num_cells = self.mesh.mesh.topology.index_map(self.mesh.vdim).size_local
-            mesh_cell_indices = np.arange(num_cells, dtype=np.int32)
-            tags_volumes = np.full(num_cells, 1, dtype=np.int32)
-            self.volume_meshtags = meshtags(
-                self.mesh.mesh, self.mesh.vdim, mesh_cell_indices, tags_volumes
-            )
 
         # check volume ids are unique
         vol_ids = [vol.id for vol in self.volume_subdomains]

festim/hydrogen_transport_problem.py

@@ -9,8 +9,6 @@
 import tqdm.auto as tqdm
 import festim as F
 
-from dolfinx.mesh import meshtags
-
 
 class HydrogenTransportProblem:
     """
@@ -474,30 +472,16 @@ def define_meshtags_and_measures(self):
             self.facet_meshtags = self.mesh.define_surface_meshtags()
             self.volume_meshtags = self.mesh.define_volume_meshtags()
 
-        elif isinstance(self.mesh, F.Mesh1D):
+        elif (
+            isinstance(self.mesh, F.Mesh)
+            and self.facet_meshtags is None
+            and self.volume_meshtags is None
+        ):
             self.facet_meshtags, self.volume_meshtags = self.mesh.define_meshtags(
                 surface_subdomains=self.surface_subdomains,
                 volume_subdomains=self.volume_subdomains,
             )
 
-        elif isinstance(self.mesh, F.Mesh):
-            if not self.facet_meshtags:
-                # create empty facet_meshtags
-                facet_indices = np.array([], dtype=np.int32)
-                facet_tags = np.array([], dtype=np.int32)
-                self.facet_meshtags = meshtags(
-                    self.mesh.mesh, self.mesh.fdim, facet_indices, facet_tags
-                )
-
-            if not self.volume_meshtags:
-                # create meshtags with all cells tagged as 1
-                num_cells = self.mesh.mesh.topology.index_map(self.mesh.vdim).size_local
-                mesh_cell_indices = np.arange(num_cells, dtype=np.int32)
-                tags_volumes = np.full(num_cells, 1, dtype=np.int32)
-                self.volume_meshtags = meshtags(
-                    self.mesh.mesh, self.mesh.vdim, mesh_cell_indices, tags_volumes
-                )
-
         # check volume ids are unique
         vol_ids = [vol.id for vol in self.volume_subdomains]
         if len(vol_ids) != len(np.unique(vol_ids)):

festim/mesh/mesh.py

@@ -1,5 +1,7 @@
 import ufl
 import dolfinx
+import numpy as np
+from dolfinx.mesh import meshtags
 
 
 class Mesh:
@@ -52,3 +54,42 @@ def fdim(self):
     @property
     def n(self):
         return ufl.FacetNormal(self.mesh)
+
+    def define_meshtags(self, surface_subdomains, volume_subdomains):
+        """Defines the facet and volume meshtags of the mesh
+
+        Args:
+            surface_subdomains (list of festim.SufaceSubdomains): the surface subdomains of the model
+            volume_subdomains (list of festim.VolumeSubdomains): the volume subdomains of the model
+
+        Returns:
+            dolfinx.mesh.MeshTags: the facet meshtags
+            dolfinx.mesh.MeshTags: the volume meshtags
+        """
+        # find all cells in domain and mark them as 0
+        num_cells = self.mesh.topology.index_map(self.vdim).size_local
+        mesh_cell_indices = np.arange(num_cells, dtype=np.int32)
+        tags_volumes = np.full(num_cells, 0, dtype=np.int32)
+
+        # find all facets in domain and mark them as 0
+        num_facets = self.mesh.topology.index_map(self.fdim).size_local
+        mesh_facet_indices = np.arange(num_facets, dtype=np.int32)
+        tags_facets = np.full(num_facets, 0, dtype=np.int32)
+
+        for surf in surface_subdomains:
+            # find all facets in subdomain and mark them as surf.id
+            entities = surf.locate_boundary_facet_indices(self.mesh)
+            tags_facets[entities] = surf.id
+
+        for vol in volume_subdomains:
+            # find all cells in subdomain and mark them as vol.id
+            entities = vol.locate_subdomain_entities(self.mesh, self.vdim)
+            tags_volumes[entities] = vol.id
+
+        # define mesh tags
+        facet_meshtags = meshtags(self.mesh, self.fdim, mesh_facet_indices, tags_facets)
+        volume_meshtags = meshtags(
+            self.mesh, self.vdim, mesh_cell_indices, tags_volumes
+        )
+
+        return facet_meshtags, volume_meshtags

festim/mesh/mesh_1d.py

@@ -4,7 +4,6 @@
 import ufl
 import numpy as np
 import festim as F
-from dolfinx.mesh import meshtags
 
 
 class Mesh1D(F.Mesh):
@@ -74,44 +73,6 @@ def check_borders(self, volume_subdomains):
             raise ValueError("borders dont match domain borders")
 
     def define_meshtags(self, surface_subdomains, volume_subdomains):
-        """Defines the facet and volume meshtags of the mesh
-
-        Args:
-            surface_subdomains (list of festim.SufaceSubdomains): the surface subdomains of the model
-            volume_subdomains (list of festim.VolumeSubdomains): the volume subdomains of the model
-
-        Returns:
-            dolfinx.mesh.MeshTags: the facet meshtags
-            dolfinx.mesh.MeshTags: the volume meshtags
-        """
-        facet_indices, tags_facets = [], []
-
-        # find all cells in domain and mark them as 0
-        num_cells = self.mesh.topology.index_map(self.vdim).size_local
-        mesh_cell_indices = np.arange(num_cells, dtype=np.int32)
-        tags_volumes = np.full(num_cells, 0, dtype=np.int32)
-
-        for surf in surface_subdomains:
-            facet_index = surf.locate_boundary_facet_indices(self.mesh, self.fdim)
-            facet_indices.append(facet_index)
-            tags_facets.append(surf.id)
-
-        for vol in volume_subdomains:
-            # find all cells in subdomain and mark them as sub_dom.id
-            entities = vol.locate_subdomain_entities(self.mesh, self.vdim)
-            tags_volumes[entities] = vol.id
-
         # check if all borders are defined
         self.check_borders(volume_subdomains)
-
-        # dofs and tags need to be in np.in32 format for meshtags
-        facet_indices = np.array(facet_indices, dtype=np.int32)
-        tags_facets = np.array(tags_facets, dtype=np.int32)
-
-        # define mesh tags
-        facet_meshtags = meshtags(self.mesh, self.fdim, facet_indices, tags_facets)
-        volume_meshtags = meshtags(
-            self.mesh, self.vdim, mesh_cell_indices, tags_volumes
-        )
-
-        return facet_meshtags, volume_meshtags
+        return super().define_meshtags(surface_subdomains, volume_subdomains)

festim/subdomain/surface_subdomain.py

@@ -1,3 +1,7 @@
+import dolfinx.mesh
+import numpy as np
+
+
 class SurfaceSubdomain:
     """
     Surface subdomain class
@@ -9,6 +13,24 @@ class SurfaceSubdomain:
     def __init__(self, id):
         self.id = id
 
+    def locate_boundary_facet_indices(self, mesh):
+        """Locates the dof of the surface subdomain within the function space
+        and return the index of the dof
+
+        Args:
+            mesh (dolfinx.mesh.Mesh): the mesh of the simulation
+
+        Returns:
+            index (np.array): the first value in the list of surface facet
+                indices of the subdomain
+        """
+        fdim = mesh.topology.dim - 1
+        # By default, all entities are included
+        indices = dolfinx.mesh.locate_entities_boundary(
+            mesh, fdim, lambda x: np.full(x.shape[1], True, dtype=bool)
+        )
+        return indices
+
 
 def find_surface_from_id(id: int, surfaces: list):
     """Returns the correct surface subdomain object from a list of surface ids

festim/subdomain/surface_subdomain_1d.py

@@ -23,19 +23,20 @@ def __init__(self, id, x) -> None:
         super().__init__(id)
         self.x = x
 
-    def locate_boundary_facet_indices(self, mesh, fdim):
+    def locate_boundary_facet_indices(self, mesh):
         """Locates the dof of the surface subdomain within the function space
         and return the index of the dof
 
         Args:
             mesh (dolfinx.mesh.Mesh): the mesh of the simulation
-            fdim (int): the dimension of the model facets
 
         Returns:
             index (np.array): the first value in the list of surface facet
                 indices of the subdomain
         """
+        assert mesh.geometry.dim == 1, "This method is only for 1D meshes"
+        fdim = 0
         indices = dolfinx.mesh.locate_entities_boundary(
             mesh, fdim, lambda x: np.isclose(x[0], self.x)
         )
-        return indices[0]
+        return indices

festim/subdomain/volume_subdomain.py

@@ -1,3 +1,8 @@
+from dolfinx.mesh import locate_entities
+
+import numpy as np
+
+
 class VolumeSubdomain:
     """
     Volume subdomain class
@@ -10,6 +15,22 @@ def __init__(self, id, material):
         self.id = id
         self.material = material
 
+    def locate_subdomain_entities(self, mesh, vdim):
+        """Locates all cells in subdomain borders within domain
+
+        Args:
+            mesh (dolfinx.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
+        """
+        # By default, all entities are included
+        # return array like x full of True
+        entities = locate_entities(mesh, vdim, lambda x: np.full(x.shape[1], True))
+        return entities
+
 
 def find_volume_from_id(id: int, volumes: list):
     """Returns the correct volume subdomain object from a list of volume ids