diff --git a/.github/CODEOWNERS b/.github/CODEOWNERS
index a7ff711ca..76e1d6dad 100644
--- a/.github/CODEOWNERS
+++ b/.github/CODEOWNERS
@@ -1,2 +1,5 @@
 # Set Ivar Stefansson and Eirik Keilegavlen as global owners
 * @keileg @IvarStefansson
+
+# Make vlipovac code owner of the docs folder
+docs/* @vlipovac
diff --git a/mypy.ini b/mypy.ini
index b716e3f9a..c2089e799 100644
--- a/mypy.ini
+++ b/mypy.ini
@@ -6,15 +6,9 @@
 mypy_path = src/porepy
 
 
-# EK: It is not at all clear to me what to do with various libraries such as numpy etc.
+# EK: It is not at all clear to me what to do with various libraries such as scipy etc.
 # For the moment, ignore some of them (seems to be done as below), and hope my understanding
 # improves in the future.
-[mypy-numpy]
-ignore_missing_imports = True
-[mypy-numpy.matlib]
-ignore_missing_imports = True
-[mypy-numpy.linalg]
-ignore_missing_imports = True
 
 [mypy-scipy]
 ignore_missing_imports = True
diff --git a/src/porepy/fracs/fracture_importer.py b/src/porepy/fracs/fracture_importer.py
index 82465e4c0..59183c2e2 100644
--- a/src/porepy/fracs/fracture_importer.py
+++ b/src/porepy/fracs/fracture_importer.py
@@ -1,4 +1,5 @@
 import csv
+from typing import Optional
 
 import gmsh
 import numpy as np
@@ -403,7 +404,9 @@ def dfm_3d_from_fab(
         return mdg
 
 
-def network_3d_from_fab(f_name: str, return_all: bool = False, tol: float = None):
+def network_3d_from_fab(
+    f_name: str, return_all: bool = False, tol: Optional[float] = None
+):
     """Read fractures from a .fab file, as specified by FracMan.
 
     The filter is based on the .fab-files available at the time of writing, and
diff --git a/src/porepy/fracs/fracture_network_2d.py b/src/porepy/fracs/fracture_network_2d.py
index ec53b1933..1467d0a90 100644
--- a/src/porepy/fracs/fracture_network_2d.py
+++ b/src/porepy/fracs/fracture_network_2d.py
@@ -1,11 +1,13 @@
 """
 Module contains class for representing a fracture network in a 2d domain.
 """
+from __future__ import annotations
+
 import copy
 import csv
 import logging
 import time
-from typing import Dict, List, Optional, Tuple, Union
+from typing import Optional, Tuple, Union
 
 import meshio
 import numpy as np
@@ -14,12 +16,13 @@
 import porepy.fracs.simplex
 from porepy.fracs import tools
 
-from .gmsh_interface import GmshData2d, GmshWriter, Tags
+from .gmsh_interface import GmshData2d, GmshWriter
+from .gmsh_interface import Tags as GmshInterfaceTags
 
 logger = logging.getLogger(__name__)
 
 
-class FractureNetwork2d(object):
+class FractureNetwork2d:
     """Class representation of a set of fractures in a 2D domain.
 
     The fractures are represented by their endpoints. Poly-line fractures are
@@ -34,59 +37,86 @@ class FractureNetwork2d(object):
     between these components may change in the future, specifically, utility functions
     may be removed.
 
-    Attributes:
-        pts (np.array, 2 x num_pts): Start and endpoints of the fractures. Points
+    Parameters:
+        pts (siz: 2 x num_pts): Start and endpoints of the fractures. Points
             can be shared by fractures.
-        edges (np.array, (2 + num_tags) x num_fracs): The first two rows represent
-            indices, referring to pts, of the start and end points of the fractures.
-            Additional rows are optional tags of the fractures.
-        domain (dictionary or np.ndarray): The domain in which the fracture set is
-            defined. If dictionary, it should contain keys 'xmin', 'xmax', 'ymin',
-            'ymax', each of which maps to a double giving the range of the domain.
-            If np.array, it should be of size 2 x n, and given the vertexes of the
-            domain. The fractures need not lay inside the domain.
-        tol (double): Tolerance used in geometric computations.
-        tags (dict): Tags for fractures.
-        _decomposition (dict): Decomposition of the fracture network, used for export to
-            gmsh, and available for later processing. Initially empty, is created by
-            self.mesh().
+        edges (size: (2 + num_tags) x num_fracs): Fractures, defined as connections
+            between the points.
+        domain: The domain in which the fracture set is defined. If dictionary, it
+            should contain keys 'xmin', 'xmax', 'ymin', 'ymax', each of which maps to a
+            double giving the range of the domain. If np.array, it should be of size
+            2 x n, and given the vertexes of the domain. The fractures need not lay
+            inside the domain.
+        tol: Tolerance used in geometric computations.
 
     """
 
-    def __init__(self, pts=None, edges=None, domain=None, tol=1e-8):
+    def __init__(
+        self,
+        pts: Optional[np.ndarray] = None,
+        edges: Optional[np.ndarray] = None,
+        domain: Optional[dict | np.ndarray] = None,
+        tol: float = 1e-8,
+    ) -> None:
         """Define the fracture set.
 
         Parameters:
             pts (np.array, 2 x n): Start and endpoints of the fractures. Points
             can be shared by fractures.
         edges (np.array, (2 + num_tags) x num_fracs): The first two rows represent
-            indices, referring to pts, of the start and end points of the fractures.
+            indices, refering to pts, of the start and end points of the fractures.
             Additional rows are optional tags of the fractures.
-        domain (dictionary or set of points): The domain in which the fracture set is
-
-             defined. See self.attributes for description.
         tol (double, optional): Tolerance used in geometric computations. Defaults to
             1e-8.
 
         """
-        if pts is None:
-            self.pts = np.zeros((2, 0))
-        else:
-            self.pts = pts
-        if edges is None:
-            self.edges = np.zeros((2, 0), dtype=int)
-        else:
-            self.edges = edges
+        self.pts = np.zeros((2, 0)) if pts is None else pts
+        """Start and endpoints of the fractures. Points can be shared by fractures."""
+
+        self.edges = np.zeros((2, 0), dtype=int) if edges is None else edges
+        """The fractures as an array of start and end points, referring to ``pts``
+
+        Additional rows are optional tags of the fractures. In the standard form, the
+        third row (first row of tags) identifies the type of edges, referring to the
+        numbering system in GmshInterfaceTags. The second row of tags keeps track of the
+        numbering of the edges (referring to the original order of the edges) in
+        geometry processing like intersection removal. Additional tags can be assigned
+        by the user.
+
+        """
+
         self.domain = domain
+        """The domain for this fracture network.
+
+        The domain is defined by a dictionary with keys 'xmin', 'xmax', 'ymin', 'ymax'.
+        If not specified, the domain will be set to the bounding box of the fractures.
+
+        """
+
         self.tol = tol
+        """Tolerance used in geometric computations."""
+
+        self.tags: dict[int | str, np.ndarray] = dict()
+        """Tags for the fractures."""
+        # TODO: The current system of tags is a bit confusing, there is both self.tags
+        # and the tags located in self.edges. The latter is used for the gmsh interface,
+        # and there may be inconsistencies in the transfer of information between the
+        # two systems.
+
+        self.bounding_box_imposed: bool = False
+        """Flag indicating whether the bounding box has been imposed."""
 
-        self.tags = {}
-        self.bounding_box_imposed = False
-        self._decomposition = {}
+        ## PRIVATE
+
+        self._decomposition: dict = dict()
+        """Dictionary of geometric information obtained from the meshing process.
+
+        This will include intersection points identified.
+        """
 
         if pts is None and edges is None:
             logger.info("Generated empty fracture set")
-        else:
+        elif pts is not None and edges is not None:
             logger.info("Generated a fracture set with %i fractures", self.num_frac())
             if pts.size > 0:
                 logger.info(
@@ -99,6 +129,10 @@ def __init__(self, pts=None, edges=None, domain=None, tol=1e-8):
                     pts[0].max(),
                     pts[1].max(),
                 )
+        else:
+            raise ValueError(
+                "Specify both points and connections for a 2d fracture network."
+            )
         if domain is not None:
             logger.info("Domain specification :" + str(domain))
 
@@ -181,18 +215,18 @@ def add_network(self, fs):
     def mesh(
         self,
         mesh_args,
-        tol=None,
-        do_snap=True,
-        constraints=None,
-        file_name=None,
-        dfn=False,
-        tags_to_transfer: Optional[List[str]] = None,
-        remove_small_fractures=False,
+        tol: Optional[float] = None,
+        do_snap: bool = True,
+        constraints: Optional[np.ndarray] = None,
+        file_name: Optional[str] = None,
+        dfn: bool = False,
+        tags_to_transfer: Optional[list[str]] = None,
+        remove_small_fractures: bool = False,
         write_geo: bool = True,
         finalize_gmsh: bool = True,
         clear_gmsh: bool = False,
         **kwargs,
-    ):
+    ) -> pp.MixedDimensionalGrid:
         """Create MixedDimensionalGrid (mixed-dimensional grid) for this fracture network.
 
         Parameters:
@@ -228,6 +262,12 @@ def mesh(
         """
         if file_name is None:
             file_name = "gmsh_frac_file.msh"
+        # No constraints if not available.
+        if constraints is None:
+            constraints = np.empty(0, dtype=int)
+        else:
+            constraints = np.atleast_1d(constraints)
+        assert isinstance(constraints, np.ndarray)
 
         gmsh_repr = self.prepare_for_gmsh(
             mesh_args, tol, do_snap, constraints, dfn, remove_small_fractures
@@ -274,13 +314,13 @@ def mesh(
 
     def prepare_for_gmsh(
         self,
-        mesh_args,
-        tol=None,
-        do_snap=True,
-        constraints=None,
-        dfn=False,
-        remove_small_fractures=False,
-    ):
+        mesh_args: dict,
+        tol: Optional[float] = None,
+        do_snap: bool = True,
+        constraints: Optional[np.ndarray] = None,
+        dfn: bool = False,
+        remove_small_fractures: bool = False,
+    ) -> GmshData2d:
         """Process network intersections and write a gmsh .geo configuration file,
         ready to be processed by gmsh.
 
@@ -313,6 +353,7 @@ def prepare_for_gmsh(
             constraints = np.empty(0, dtype=int)
         else:
             constraints = np.atleast_1d(constraints)
+        assert isinstance(constraints, np.ndarray)
         constraints = np.sort(constraints)
 
         p = self.pts
@@ -395,13 +436,13 @@ def prepare_for_gmsh(
 
         # map the constraint index
         index_map = np.where(np.logical_not(to_delete))[0]
-        constraints = np.arange(index_map.size)[np.in1d(index_map, constraints)]
+        mapped_constraints = np.arange(index_map.size)[np.in1d(index_map, constraints)]
 
         # update the tags
         for key, value in self.tags.items():
             self.tags[key] = np.delete(value, to_delete)
 
-        self._find_and_split_intersections(constraints)
+        self._find_and_split_intersections(mapped_constraints)
         # Insert auxiliary points and determine mesh size.
         # _insert_auxiliary_points(..) does both.
         # _set_mesh_size_without_auxiliary_points() sets the mesh size
@@ -411,39 +452,65 @@ def prepare_for_gmsh(
             self._insert_auxiliary_points(**mesh_args)
         else:
             self._set_mesh_size_without_auxiliary_points(**mesh_args)
-        # Transfer data to the format expected by the gmsh interface
+
+        # Transfer data to the format expected by the gmsh interface.
+        # This requires some information processing and translation between data
+        # formats: In the geometry processing undertaken up to this point, it has been
+        # convenient to use numerical values for identifying the different line types
+        # (fracture, constraint, boundary). For the Gmsh processing, a string-based
+        # system is used, as this is more readable and closer to the system employed in
+        # Gmsh. In practice, this requires translating from GmshInterfaceTags values to
+        # names.
+        # In addition to this translation, the below code also does some interpretation
+        # of the information obtained during geometry processing.
+
         decomp = self._decomposition
 
         edges = decomp["edges"]
-
-        # Only lines specified as fractures are defined as physical
-        phys_line_tags = {}
+        # Information about line types is found in the third row of edges
         edge_types = edges[2]
+
+        # Process information about lines that should be tagged as physical by Gmsh.
+        # These are fractures, domain boundaries and auxiliary (constraints).
+        # phys_line_tags is a mapping from line index to the Tag.
+        phys_line_tags: dict[int, GmshInterfaceTags] = {}
+
         for ei, tag in enumerate(edge_types):
             if tag in (
-                Tags.FRACTURE.value,
-                Tags.DOMAIN_BOUNDARY_LINE.value,
-                Tags.AUXILIARY_LINE.value,
+                GmshInterfaceTags.FRACTURE.value,
+                GmshInterfaceTags.DOMAIN_BOUNDARY_LINE.value,
+                GmshInterfaceTags.AUXILIARY_LINE.value,
             ):
-                phys_line_tags[ei] = tag
-
-        phys_point_tags = {
-            i: Tags.FRACTURE_INTERSECTION_POINT.value for i in decomp["intersections"]
+                # Note: phys_line_tags contains the GmshInterfaceTags instead of
+                # the numbers in edges[2].
+                phys_line_tags[ei] = GmshInterfaceTags(tag)
+
+        # Tag all points that have been defined as intersections between fractures.
+        # phys_point_tags is a mapping from the point index to the tag.
+        phys_point_tags: dict[int, GmshInterfaceTags] = {
+            i: GmshInterfaceTags.FRACTURE_INTERSECTION_POINT
+            for i in decomp["intersections"]
         }
 
-        point_on_fracture = edges[:2, edge_types == Tags.FRACTURE.value].ravel()
+        # Find points on the boundary, and mark these as physical points.
         point_on_boundary = edges[
-            :2, edge_types == Tags.DOMAIN_BOUNDARY_LINE.value
+            :2, edge_types == GmshInterfaceTags.DOMAIN_BOUNDARY_LINE.value
         ].ravel()
-        fracture_boundary_points = np.intersect1d(point_on_fracture, point_on_boundary)
-
         phys_point_tags.update(
-            {pi: Tags.DOMAIN_BOUNDARY_POINT.value for pi in point_on_boundary}
+            {pi: GmshInterfaceTags.DOMAIN_BOUNDARY_POINT for pi in point_on_boundary}
         )
 
-        # register fracture boundary points, override previously added boundary points
+        # Find points that are both on the boundary and on a fracture. These have
+        # a special tag, thus override the values set for normal boundary points.
+        point_on_fracture = edges[
+            :2, edge_types == GmshInterfaceTags.FRACTURE.value
+        ].ravel()
+        fracture_boundary_points = np.intersect1d(point_on_fracture, point_on_boundary)
         phys_point_tags.update(
-            {pi: Tags.FRACTURE_BOUNDARY_POINT.value for pi in fracture_boundary_points}
+            {
+                pi: GmshInterfaceTags.FRACTURE_BOUNDARY_POINT
+                for pi in fracture_boundary_points
+            }
         )
 
         data = GmshData2d(
@@ -455,7 +522,7 @@ def prepare_for_gmsh(
         )
         return data
 
-    def _find_and_split_intersections(self, constraints: np.ndarray):
+    def _find_and_split_intersections(self, constraints: np.ndarray) -> None:
         """Unified description of points and lines for domain and fractures.
 
         FIXME: update documentation
@@ -477,9 +544,11 @@ def _find_and_split_intersections(self, constraints: np.ndarray):
 
         tags = np.zeros((2, edges.shape[1]), dtype=int)
 
-        tags[0][np.logical_not(self.tags["boundary"])] = Tags.FRACTURE.value
-        tags[0][self.tags["boundary"]] = Tags.DOMAIN_BOUNDARY_LINE.value
-        tags[0][constraints] = Tags.AUXILIARY_LINE.value
+        tags[0][
+            np.logical_not(self.tags["boundary"])
+        ] = GmshInterfaceTags.FRACTURE.value
+        tags[0][self.tags["boundary"]] = GmshInterfaceTags.DOMAIN_BOUNDARY_LINE.value
+        tags[0][constraints] = GmshInterfaceTags.AUXILIARY_LINE.value
 
         tags[1] = np.arange(edges.shape[1])
 
@@ -556,17 +625,17 @@ def _find_and_split_intersections(self, constraints: np.ndarray):
             }
         )
 
-    def _find_intersection_points(self, lines):
+    def _find_intersection_points(self, lines: np.ndarray) -> np.ndarray:
 
-        frac_id = np.ravel(lines[:2, lines[2] == Tags.FRACTURE.value])
+        frac_id = np.ravel(lines[:2, lines[2] == GmshInterfaceTags.FRACTURE.value])
         _, frac_ia, frac_count = np.unique(frac_id, True, False, True)
 
         # In the case we have auxiliary points remove do not create a 0d point in
         # case one intersects a single fracture. In the case of multiple fractures intersection
         # with an auxiliary point do consider the 0d.
         aux_id = np.logical_or(
-            lines[2] == Tags.AUXILIARY_LINE.value,
-            lines[2] == Tags.DOMAIN_BOUNDARY_LINE.value,
+            lines[2] == GmshInterfaceTags.AUXILIARY_LINE.value,
+            lines[2] == GmshInterfaceTags.DOMAIN_BOUNDARY_LINE.value,
         )
         if np.any(aux_id):
             aux_id = np.ravel(lines[:2, aux_id])
@@ -582,8 +651,11 @@ def _find_intersection_points(self, lines):
         return frac_id[frac_ia[frac_count > 1]]
 
     def _insert_auxiliary_points(
-        self, mesh_size_frac=None, mesh_size_bound=None, mesh_size_min=None
-    ):
+        self,
+        mesh_size_frac: Optional[float] = None,
+        mesh_size_bound: Optional[float] = None,
+        mesh_size_min: Optional[float] = None,
+    ) -> None:
         # Mesh size
         # Tag points at the domain corners
         logger.info("Determine mesh size")
@@ -641,7 +713,7 @@ def _set_mesh_size_without_auxiliary_points(
 
     def impose_external_boundary(
         self,
-        domain: Optional[Union[Dict, np.ndarray]] = None,
+        domain: Optional[Union[dict, np.ndarray]] = None,
         add_domain_edges: bool = True,
     ) -> Tuple[np.ndarray, np.ndarray]:
         """
@@ -672,11 +744,12 @@ def impose_external_boundary(
             x_max = domain["xmax"]
             y_min = domain["ymin"]
             y_max = domain["ymax"]
-            dom_p = np.array(
+            dom_p: np.ndarray = np.array(
                 [[x_min, x_max, x_max, x_min], [y_min, y_min, y_max, y_max]]
             )
             dom_lines = np.array([[0, 1], [1, 2], [2, 3], [3, 0]]).T
         else:
+            assert isinstance(domain, np.ndarray)
             dom_p = domain
             tmp = np.arange(dom_p.shape[1])
             dom_lines = np.vstack((tmp, (tmp + 1) % dom_p.shape[1]))
@@ -704,7 +777,7 @@ def impose_external_boundary(
 
         # Define boundary tags. Set False to all existing edges (after cutting those
         # outside the boundary).
-        boundary_tags = self.tags.get("boundary", [False] * e.shape[1])
+        boundary_tags = self.tags.get("boundary", np.zeros(e.shape[1], dtype=bool))
 
         if add_domain_edges:
             num_p = p.shape[1]
@@ -716,13 +789,15 @@ def impose_external_boundary(
                 self.tags[key] = np.hstack(
                     (
                         value[edges_kept],
-                        Tags.DOMAIN_BOUNDARY_LINE.value
+                        GmshInterfaceTags.DOMAIN_BOUNDARY_LINE.value
                         * np.ones(dom_lines.shape[1], dtype=int),
                     )
                 )
 
             # Define the new boundary tags
-            new_boundary_tags = boundary_tags + dom_lines.shape[1] * [True]
+            new_boundary_tags = np.hstack(
+                [boundary_tags, np.ones(dom_lines.shape[1], bool)]
+            )
             self.tags["boundary"] = np.array(new_boundary_tags)
 
             self._decomposition["domain_boundary_points"] = num_p + np.arange(
@@ -1283,7 +1358,7 @@ def to_csv(self, file_name: str, with_header=True):
                 csv_writer.writerow(data)
 
     def to_file(
-        self, file_name: str, data: Dict[str, np.ndarray] = None, **kwargs
+        self, file_name: str, data: Optional[dict[str, np.ndarray]] = None, **kwargs
     ) -> None:
         """
         Export the fracture network to file.
diff --git a/src/porepy/fracs/fracture_network_3d.py b/src/porepy/fracs/fracture_network_3d.py
index 089560ebc..2c0e0bef9 100644
--- a/src/porepy/fracs/fracture_network_3d.py
+++ b/src/porepy/fracs/fracture_network_3d.py
@@ -4,11 +4,13 @@
 The model relies heavily on functions in the computational geometry library.
 
 """
+from __future__ import annotations
+
 import copy
 import csv
 import logging
 import time
-from typing import Dict, List, Optional, Tuple, Union
+from typing import Optional, Union
 
 import meshio
 import networkx as nx
@@ -16,10 +18,12 @@
 from scipy.spatial import ConvexHull
 
 import porepy as pp
-from porepy.fracs.gmsh_interface import GmshData3d, GmshWriter, Tags
+from porepy.fracs.gmsh_interface import GmshData3d, GmshWriter
 from porepy.fracs.plane_fracture import PlaneFracture
 from porepy.utils import setmembership, sort_points
 
+from .gmsh_interface import Tags as GmshInterfaceTags
+
 # Module-wide logger
 logger = logging.getLogger(__name__)
 
@@ -60,8 +64,8 @@ class FractureNetwork3d(object):
 
     def __init__(
         self,
-        fractures: Optional[List[PlaneFracture]] = None,
-        domain: Optional[Union[Dict[str, float], List[np.ndarray]]] = None,
+        fractures: Optional[list[PlaneFracture]] = None,
+        domain: Optional[Union[dict[str, float], list[np.ndarray]]] = None,
         tol: float = 1e-8,
         run_checks: bool = False,
     ) -> None:
@@ -98,7 +102,7 @@ def __init__(
         # and whether the intersection is on the boundary of the fractures.
         # Note that a Y-type intersection between three fractures is represented
         # as three intersections.
-        self.intersections: Dict[str, np.ndarray] = {
+        self.intersections: dict[str, np.ndarray] = {
             "first": np.array([], dtype=object),
             "second": np.array([], dtype=object),
             "start": np.zeros((3, 0)),
@@ -120,7 +124,7 @@ def __init__(
         self.mesh_size_frac = None
         self.mesh_size_bound = None
         # Assign an empty tag dictionary
-        self.tags: Dict[str, List[bool]] = {}
+        self.tags: dict[str, list[bool]] = {}
 
         # No auxiliary points have been added
         self.auxiliary_points_added = False
@@ -185,12 +189,12 @@ def num_frac(self) -> int:
 
     def mesh(
         self,
-        mesh_args: Dict[str, float],
+        mesh_args: dict[str, float],
         dfn: bool = False,
         file_name: Optional[str] = None,
         constraints: Optional[np.ndarray] = None,
         write_geo: bool = True,
-        tags_to_transfer: Optional[List[str]] = None,
+        tags_to_transfer: Optional[list[str]] = None,
         finalize_gmsh: bool = True,
         clear_gmsh: bool = False,
         **kwargs,
@@ -383,7 +387,7 @@ def prepare_for_gmsh(self, mesh_args, dfn=False, constraints=None) -> GmshData3d
         auxiliary_line[some_boundary_edge] = False
         # The edge tags for internal lines were set accordingly in self._classify_edges.
         # Update to auxiliary line if this was really what we had.
-        edge_tags[auxiliary_line] = Tags.AUXILIARY_LINE.value
+        edge_tags[auxiliary_line] = GmshInterfaceTags.AUXILIARY_LINE.value
 
         # .. and we're done with edges (until someone defines a new special case)
         # Next, find intersection points.
@@ -405,8 +409,8 @@ def prepare_for_gmsh(self, mesh_args, dfn=False, constraints=None) -> GmshData3d
 
         # We should also consider points that are both on a fracture tip and a domain
         # boundary line - these will be added to the fracture and boundary points below.
-        domain_boundary_line = edge_tags == Tags.DOMAIN_BOUNDARY_LINE.value
-        fracture_tip_line = edge_tags == Tags.FRACTURE_TIP.value
+        domain_boundary_line = edge_tags == GmshInterfaceTags.DOMAIN_BOUNDARY_LINE.value
+        fracture_tip_line = edge_tags == GmshInterfaceTags.FRACTURE_TIP.value
         domain_boundary_point = edges[:2, domain_boundary_line]
         fracture_tip_point = edges[:2, fracture_tip_line]
         domain_and_tip_point = np.intersect1d(domain_boundary_point, fracture_tip_point)
@@ -479,17 +483,21 @@ def prepare_for_gmsh(self, mesh_args, dfn=False, constraints=None) -> GmshData3d
         # but it can be useful to mark them for other purposes (EK: DFM upscaling)
         point_tags[
             fracture_constraint_intersection
-        ] = Tags.FRACTURE_CONSTRAINT_INTERSECTION_POINT.value
+        ] = GmshInterfaceTags.FRACTURE_CONSTRAINT_INTERSECTION_POINT.value
 
-        point_tags[fracture_and_boundary_points] = Tags.FRACTURE_BOUNDARY_POINT.value
+        point_tags[
+            fracture_and_boundary_points
+        ] = GmshInterfaceTags.FRACTURE_BOUNDARY_POINT.value
 
         # We're done! Hurrah!
 
         # Find points tagged as on the domain boundary
-        boundary_points = np.where(point_tags == Tags.DOMAIN_BOUNDARY_POINT.value)[0]
+        boundary_points = np.where(
+            point_tags == GmshInterfaceTags.DOMAIN_BOUNDARY_POINT.value
+        )[0]
 
         fracture_boundary_points = np.where(
-            point_tags == Tags.FRACTURE_BOUNDARY_POINT.value
+            point_tags == GmshInterfaceTags.FRACTURE_BOUNDARY_POINT.value
         )[0]
 
         # Intersections on the boundary should not have a 0d grid assigned
@@ -504,48 +512,53 @@ def prepare_for_gmsh(self, mesh_args, dfn=False, constraints=None) -> GmshData3d
 
         line_in_poly = self.decomposition["line_in_frac"]
 
-        physical_points: Dict[int, Tags] = {}
+        physical_points: dict[int, GmshInterfaceTags] = {}
         for pi in fracture_boundary_points:
-            physical_points[pi] = Tags.FRACTURE_BOUNDARY_POINT
+            physical_points[pi] = GmshInterfaceTags.FRACTURE_BOUNDARY_POINT
 
         for pi in fracture_constraint_intersection:
-            physical_points[pi] = Tags.FRACTURE_CONSTRAINT_INTERSECTION_POINT
+            physical_points[
+                pi
+            ] = GmshInterfaceTags.FRACTURE_CONSTRAINT_INTERSECTION_POINT
 
         for pi in boundary_points:
-            physical_points[pi] = Tags.DOMAIN_BOUNDARY_POINT
+            physical_points[pi] = GmshInterfaceTags.DOMAIN_BOUNDARY_POINT
 
         for pi in true_intersection_points:
-            physical_points[pi] = Tags.FRACTURE_INTERSECTION_POINT
+            physical_points[pi] = GmshInterfaceTags.FRACTURE_INTERSECTION_POINT
 
         # Use separate structures to store tags and physical names for the polygons.
         # The former is used for feeding information into gmsh, while the latter is
         # used to tag information in the output from gmsh. They are kept as separate
         # variables since a user may want to modify the physical surfaces before
         # generating the .msh file.
-        physical_surfaces: Dict[int, Tags] = {}
-        polygon_tags = np.zeros(len(self._fractures), dtype=int)
+        physical_surfaces: dict[int, GmshInterfaceTags] = dict()
+        polygon_tags: dict[int, GmshInterfaceTags] = dict()
 
         for fi, _ in enumerate(self._fractures):
+            # Translate from numerical tags to the GmshInterfaceTags system.
             if has_boundary and self.tags["boundary"][fi]:
-                physical_surfaces[fi] = Tags.DOMAIN_BOUNDARY_SURFACE
-                polygon_tags[fi] = Tags.DOMAIN_BOUNDARY_SURFACE.value
+                physical_surfaces[fi] = GmshInterfaceTags.DOMAIN_BOUNDARY_SURFACE
+                polygon_tags[fi] = GmshInterfaceTags.DOMAIN_BOUNDARY_SURFACE
             elif fi in constraints:
-                physical_surfaces[fi] = Tags.AUXILIARY_PLANE
-                polygon_tags[fi] = Tags.AUXILIARY_PLANE.value
+                physical_surfaces[fi] = GmshInterfaceTags.AUXILIARY_PLANE
+                polygon_tags[fi] = GmshInterfaceTags.AUXILIARY_PLANE
             else:
-                physical_surfaces[fi] = Tags.FRACTURE
-                polygon_tags[fi] = Tags.FRACTURE.value
-
-        physical_lines = {}
-        for ei in range(edges.shape[1]):
-            if edges[2, ei] in (
-                Tags.FRACTURE_TIP.value,
-                Tags.FRACTURE_INTERSECTION_LINE.value,
-                Tags.DOMAIN_BOUNDARY_LINE.value,
-                Tags.FRACTURE_BOUNDARY_LINE.value,
-                Tags.AUXILIARY_LINE.value,
+                physical_surfaces[fi] = GmshInterfaceTags.FRACTURE
+                polygon_tags[fi] = GmshInterfaceTags.FRACTURE
+
+        physical_lines: dict[int, GmshInterfaceTags] = dict()
+        for ei, tag in enumerate(edges[2]):
+            if tag in (
+                GmshInterfaceTags.FRACTURE_TIP.value,
+                GmshInterfaceTags.FRACTURE_INTERSECTION_LINE.value,
+                GmshInterfaceTags.DOMAIN_BOUNDARY_LINE.value,
+                GmshInterfaceTags.FRACTURE_BOUNDARY_LINE.value,
+                GmshInterfaceTags.AUXILIARY_LINE.value,
             ):
-                physical_lines[ei] = edges[2, ei]
+                # Translate from the numerical value of the line to the tag in the
+                # Gmsh interface.
+                physical_lines[ei] = GmshInterfaceTags(tag)
 
         gmsh_repr = GmshData3d(
             dim=3,
@@ -571,7 +584,7 @@ def __getitem__(self, position):
 
     def intersections_of_fracture(
         self, frac: Union[int, PlaneFracture]
-    ) -> Tuple[List[int], List[bool]]:
+    ) -> tuple[list[int], list[bool]]:
         """Get all known intersections for a fracture.
 
         If called before find_intersections(), the returned list will be empty.
@@ -1276,7 +1289,7 @@ def bounding_box(self):
 
     def impose_external_boundary(
         self,
-        domain: Optional[Union[Dict[str, float], List[np.ndarray]]] = None,
+        domain: Optional[Union[dict[str, float], list[np.ndarray]]] = None,
         keep_box: bool = True,
         area_threshold: float = 1e-4,
         clear_gmsh: bool = True,
@@ -1492,7 +1505,7 @@ def impose_external_boundary(
             # List of nodes removed (isolated by an introduced diagonal)
             removed = []
             # List of triangles introduced
-            tris: List[List[int]] = []
+            tris: list[list[int]] = []
 
             # Helper function to get i2 move to the next available node
             def _next_i2(i2):
@@ -1624,7 +1637,7 @@ def _next_i2(i2):
             # Find the pairs looping over all vertexes, find all its
             # neighboring vertexes.
             main_vertex_list = []
-            other_vertex_list: List[int] = []
+            other_vertex_list: list[int] = []
             indptr_list = [0]
 
             for i in range(num_p):
@@ -1834,9 +1847,9 @@ def _classify_edges(self, polygon_edges, constraints):
         # auxiliary type. These will have tag zero; and treated in a special
         # manner by the interface to gmsh.
         not_boundary_ind = np.setdiff1d(np.arange(num_is_bound), bound_ind)
-        tag[bound_ind] = Tags.FRACTURE_TIP.value
+        tag[bound_ind] = GmshInterfaceTags.FRACTURE_TIP.value
 
-        tag[not_boundary_ind] = Tags.FRACTURE_INTERSECTION_LINE.value
+        tag[not_boundary_ind] = GmshInterfaceTags.FRACTURE_INTERSECTION_LINE.value
 
         return tag, np.logical_not(all_bound), some_bound
 
@@ -1860,7 +1873,7 @@ def _on_domain_boundary(self, edges, edge_tags):
         )[0]
 
         # ... on the points...
-        point_tags = Tags.NEUTRAL.value * np.ones(
+        point_tags = GmshInterfaceTags.NEUTRAL.value * np.ones(
             self.decomposition["points"].shape[1], dtype=int
         )
         # and the mapping between fractures and edges.
@@ -1877,9 +1890,11 @@ def _on_domain_boundary(self, edges, edge_tags):
                 if all(edge_of_domain_boundary):
                     # The point is not associated with a fracture extending to the
                     # boundary
-                    edge_tags[e] = Tags.DOMAIN_BOUNDARY_LINE.value
+                    edge_tags[e] = GmshInterfaceTags.DOMAIN_BOUNDARY_LINE.value
                     # The points of this edge are also associated with the boundary
-                    point_tags[edges[:, e]] = Tags.DOMAIN_BOUNDARY_POINT.value
+                    point_tags[
+                        edges[:, e]
+                    ] = GmshInterfaceTags.DOMAIN_BOUNDARY_POINT.value
                 else:
                     # The edge is associated with at least one fracture. Still, if it is
                     # also the edge of at least one boundary point, we will consider it
@@ -1892,13 +1907,17 @@ def _on_domain_boundary(self, edges, edge_tags):
 
                     # The line is on the boundary
                     if on_one_domain_edge:
-                        edge_tags[e] = Tags.DOMAIN_BOUNDARY_LINE.value
-                        point_tags[edges[:, e]] = Tags.DOMAIN_BOUNDARY_POINT.value
+                        edge_tags[e] = GmshInterfaceTags.DOMAIN_BOUNDARY_LINE.value
+                        point_tags[
+                            edges[:, e]
+                        ] = GmshInterfaceTags.DOMAIN_BOUNDARY_POINT.value
                     else:
                         # The edge is an intersection between a fracture and a boundary
                         # polygon
-                        edge_tags[e] = Tags.FRACTURE_BOUNDARY_LINE.value
-                        point_tags[edges[:, e]] = Tags.FRACTURE_BOUNDARY_POINT.value
+                        edge_tags[e] = GmshInterfaceTags.FRACTURE_BOUNDARY_LINE.value
+                        point_tags[
+                            edges[:, e]
+                        ] = GmshInterfaceTags.FRACTURE_BOUNDARY_POINT.value
             else:
                 # This is not an edge on the domain boundary, and the tag assigned in
                 # in self._classify_edges() is still valid: It is either a fracture tip
@@ -1981,7 +2000,7 @@ def _determine_mesh_size(self, point_tags, **kwargs):
         mesh_size = np.minimum(mesh_size_min, self.mesh_size_frac * np.ones(num_pts))
 
         if self.mesh_size_bound is not None:
-            on_boundary = point_tags == Tags.DOMAIN_BOUNDARY_POINT.value
+            on_boundary = point_tags == GmshInterfaceTags.DOMAIN_BOUNDARY_POINT.value
             mesh_size_bound = np.minimum(
                 mesh_size_min, self.mesh_size_bound * np.ones(num_pts)
             )
@@ -2051,7 +2070,7 @@ def dist_p(a, b):
             return np.sqrt(np.sum(np.power(b - a, 2), axis=0))
 
         # Dictionary that for each fracture maps the index of all other fractures.
-        intersecting_fracs: Dict[List[int]] = {}
+        intersecting_fracs: dict[list[int]] = {}
         # Loop over all fractures
         for fi, f in enumerate(self._fractures):
 
@@ -2064,7 +2083,7 @@ def dist_p(a, b):
 
             # Keep track of which other fractures are intersecting - will be
             # needed later on
-            intersections_this_fracture: List[int] = []
+            intersections_this_fracture: list[int] = []
 
             # Loop over all intersections of the fracture
             isects, is_first_isect = self.intersections_of_fracture(f)
@@ -2383,7 +2402,10 @@ def _add_intersection(
         )
 
     def to_file(
-        self, file_name: str, data: Dict[str, Union[np.ndarray, List]] = None, **kwargs
+        self,
+        file_name: str,
+        data: Optional[dict[str, Union[np.ndarray, list]]] = None,
+        **kwargs,
     ) -> None:
         """
         Export the fracture network to file.
diff --git a/src/porepy/fracs/gmsh_interface.py b/src/porepy/fracs/gmsh_interface.py
index 41ee85e19..8ad7a3c9d 100644
--- a/src/porepy/fracs/gmsh_interface.py
+++ b/src/porepy/fracs/gmsh_interface.py
@@ -13,9 +13,11 @@
         gmsh model. Can also mesh.
 
 """
+from __future__ import annotations
+
 from dataclasses import dataclass
 from enum import Enum
-from typing import Dict, List, Optional, Tuple, Union
+from typing import Optional, Union
 
 import gmsh
 import numpy as np
@@ -145,10 +147,10 @@ class _GmshData:
     lines: np.ndarray
     # Mapping from indices of points (referring to the ordering in points) to
     # numerical tags for the points.
-    physical_points: Dict[int, Tags]
+    physical_points: dict[int, Tags]
     # Mapping from indices of lines (referring to the ordering in lines) to numerical tags
     # for the lines.
-    physical_lines: Dict[int, Tags]
+    physical_lines: dict[int, Tags]
 
 
 @dataclass
@@ -169,15 +171,15 @@ class GmshData3d(_GmshData):
 
     # Polygons (both boundary surfaces, fractures and auxiliary lines).
     # See FractureNetwork3d for examples of usage.
-    polygons: List[np.ndarray]
+    polygons: list[np.ndarray]
     # Tags used to identify the type of polygons.
-    polygon_tags: np.ndarray
+    polygon_tags: dict[int, Tags]
     # Physical name information for polygons. Used to set gmsh tags for objects to be
     # represented in hte output mesh.
-    physical_surfaces: Dict[int, Tags]
+    physical_surfaces: dict[int, Tags]
     # Lines to be embedded in surfaces. Outer list has one item per polygon, inner has
     # index of lines to embed.
-    lines_in_surface: List[List[int]]
+    lines_in_surface: list[list[int]]
 
     dim: int = 3
 
@@ -208,7 +210,7 @@ def __init__(self, data: Union[GmshData2d, GmshData3d]) -> None:
 
         self.define_geometry()
 
-    def set_gmsh_options(self, options: Optional[Dict] = None) -> None:
+    def set_gmsh_options(self, options: Optional[dict] = None) -> None:
         """Set Gmsh options. See Gmsh documentation for choices available.
 
         Parameters:
@@ -350,7 +352,7 @@ def generate(
             gmsh.finalize()
             GmshWriter.gmsh_initialized = False
 
-    def _add_points(self) -> List[int]:
+    def _add_points(self) -> list[int]:
         """Add points to gmsh. Also set physical names for points if required."""
         point_tags = []
         for pi, sz in enumerate(self._data.mesh_size):
@@ -379,6 +381,8 @@ def _add_fractures_2d(self) -> None:
         """Add all fracture lines to gmsh."""
         # Only add the lines that correspond to fractures or auxiliary lines.
         # Boundary lines are added elsewhere.
+        # NOTE: Here we operate on the numerical tagging of lines (the attribute edges
+        # in FractureNetwork2d), thus we must compare with the values of the Tags.
         inds = np.argwhere(
             np.logical_or.reduce(
                 (
@@ -389,19 +393,23 @@ def _add_fractures_2d(self) -> None:
         ).ravel()
         self._frac_tags = self._add_lines(inds, embed_in_domain=True)
 
-    def _add_fractures_3d(self):
+    def _add_fractures_3d(self) -> None:
         """Add fracture polygons to gmsh"""
-        inds = np.argwhere(
-            np.logical_or.reduce(
-                (
-                    self._data.polygon_tags == Tags.FRACTURE.value,
-                    self._data.polygon_tags == Tags.AUXILIARY_PLANE.value,
-                )
-            )
-        ).ravel()
-        self._frac_tags = self._add_polygons_3d(inds, embed_in_domain=True)
 
-    def _add_polygons_3d(self, inds: np.ndarray, embed_in_domain: bool) -> List[int]:
+        if not isinstance(self._data, GmshData3d):
+            raise ValueError("Need 3d geometry to write fractures")
+
+        # Loop over the polygons, find those tagged as fractures or auxiliary planes
+        # (which could be constraints in the meshing). Add these to the information
+        # to be written.
+        indices_to_write = []
+        for ind, tag in self._data.polygon_tags.items():
+            if tag in (Tags.FRACTURE, Tags.AUXILIARY_PLANE):
+                indices_to_write.append(ind)
+
+        self._frac_tags = self._add_polygons_3d(indices_to_write, embed_in_domain=True)
+
+    def _add_polygons_3d(self, inds: list, embed_in_domain: bool) -> list[int]:
         """Generic method to add polygons to a 3d geometry"""
 
         if not isinstance(self._data, GmshData3d):
@@ -419,7 +427,7 @@ def _add_polygons_3d(self, inds: np.ndarray, embed_in_domain: bool) -> List[int]
 
         surf_tag = []
 
-        to_phys_tags: List[Tuple[int, str, List[int]]] = []
+        to_phys_tags: list[tuple[int, str, list[int]]] = []
 
         for pi in inds:
             line_tags = []
@@ -462,9 +470,9 @@ def _add_polygons_3d(self, inds: np.ndarray, embed_in_domain: bool) -> List[int]
 
         return surf_tag
 
-    def _add_lines(self, ind: np.ndarray, embed_in_domain: bool) -> List[int]:
+    def _add_lines(self, ind: np.ndarray, embed_in_domain: bool) -> list[int]:
         # Helper function to write lines.
-        line_tags: List[int] = []
+        line_tags: list[int] = []
 
         if ind.size == 0:
             return line_tags
@@ -483,7 +491,7 @@ def _add_lines(self, ind: np.ndarray, embed_in_domain: bool) -> List[int]:
 
         # Temporary storage of the lines that are to be assigned physical
         # groups
-        to_physical_group: List[Tuple[int, str, List[int]]] = []
+        to_physical_group: list[tuple[int, str, list[int]]] = []
 
         for i in range_id:
             loc_tags = []
@@ -497,9 +505,9 @@ def _add_lines(self, ind: np.ndarray, embed_in_domain: bool) -> List[int]:
             # Store local tags
             line_tags += loc_tags
 
-            # Add this line to the set of physical groups to be assigned
+            # Add this line to the set of physical groups to be assigned.
             # We do not assign physical groupings inside this for-loop, as this would
-            # require multiple costly synchronizations (gmsh style)
+            # require multiple costly synchronizations (gmsh style).
             if has_physical_lines and i in self._data.physical_lines:
                 to_physical_group.append((i, physical_name, loc_tags))
 
@@ -520,6 +528,9 @@ def _add_lines(self, ind: np.ndarray, embed_in_domain: bool) -> List[int]:
 
     def _add_domain_2d(self) -> int:
         """Write boundary lines, and tie them together to a closed loop. Define domain."""
+
+        # Here we operate on ``FractureNetwork2d.edges``, thus we should use the
+        # numerical values of the tag for comparison.
         bound_line_ind = np.argwhere(
             self._data.lines[2] == Tags.DOMAIN_BOUNDARY_LINE.value
         ).ravel()
@@ -541,9 +552,10 @@ def _add_domain_3d(self) -> int:
         if not isinstance(self._data, GmshData3d):
             raise ValueError("Need 3d geometry to specify 3d domain")
 
-        inds = np.argwhere(
-            self._data.polygon_tags == Tags.DOMAIN_BOUNDARY_SURFACE.value
-        ).ravel()
+        inds = []
+        for i, tag in self._data.polygon_tags.items():
+            if tag == Tags.DOMAIN_BOUNDARY_SURFACE:
+                inds.append(i)
 
         bound_surf_tags = self._add_polygons_3d(inds, embed_in_domain=False)
 
diff --git a/src/porepy/fracs/meshing.py b/src/porepy/fracs/meshing.py
index 2f6945991..8bacf2326 100644
--- a/src/porepy/fracs/meshing.py
+++ b/src/porepy/fracs/meshing.py
@@ -11,6 +11,7 @@
 
 import logging
 import time
+from typing import Optional
 
 import numpy as np
 import scipy.sparse as sps
@@ -28,7 +29,7 @@
 
 def subdomains_to_mdg(
     subdomains: list[list[pp.Grid]],
-    time_tot: float = None,
+    time_tot: Optional[float] = None,
     **kwargs,
 ) -> pp.MixedDimensionalGrid:
     """Convert a list of grids to a full MixedDimensionalGrid.
diff --git a/src/porepy/fracs/msh_2_grid.py b/src/porepy/fracs/msh_2_grid.py
index d76ec3b66..482fe40cc 100644
--- a/src/porepy/fracs/msh_2_grid.py
+++ b/src/porepy/fracs/msh_2_grid.py
@@ -1,7 +1,9 @@
 """
 Module for converting gmsh output file to our grid structure.
 """
-from typing import Dict, List, Tuple, Union
+from __future__ import annotations
+
+from typing import Optional, Union
 
 import numpy as np
 
@@ -10,7 +12,7 @@
 from .gmsh_interface import PhysicalNames
 
 
-def create_3d_grids(pts: np.ndarray, cells: Dict[str, np.ndarray]) -> List[pp.Grid]:
+def create_3d_grids(pts: np.ndarray, cells: dict[str, np.ndarray]) -> list[pp.Grid]:
     """Create a tetrahedral grid from a gmsh tessalation.
 
     Parameters:
@@ -37,13 +39,13 @@ def create_3d_grids(pts: np.ndarray, cells: Dict[str, np.ndarray]) -> List[pp.Gr
 
 def create_2d_grids(
     pts: np.ndarray,
-    cells: Dict[str, np.ndarray],
-    phys_names: Dict[str, str],
-    cell_info: Dict,
+    cells: dict[str, np.ndarray],
+    phys_names: dict[str, str],
+    cell_info: dict,
     is_embedded: bool = False,
-    surface_tag: str = None,
-    constraints: np.ndarray = None,
-) -> List[pp.Grid]:
+    surface_tag: Optional[str] = None,
+    constraints: Optional[np.ndarray] = None,
+) -> list[pp.Grid]:
     """Create 2d grids for lines of a specified type from a gmsh tessalation.
 
     Only surfaces that were defined as 'physical' in the gmsh sense may have a grid
@@ -85,7 +87,7 @@ def create_2d_grids(
 
     if is_embedded:
         # List of 2D grids, one for each surface
-        g_list: List[pp.Grid] = []
+        g_list: list[pp.Grid] = []
 
         # Special treatment of the case with no fractures
         if "triangle" not in cells:
@@ -199,14 +201,14 @@ def create_2d_grids(
 
 def create_1d_grids(
     pts: np.ndarray,
-    cells: Dict[str, np.ndarray],
-    phys_names: Dict,
-    cell_info: Dict,
-    line_tag: str = None,
+    cells: dict[str, np.ndarray],
+    phys_names: dict,
+    cell_info: dict,
+    line_tag: Optional[str] = None,
     tol: float = 1e-4,
-    constraints: np.ndarray = None,
+    constraints: Optional[np.ndarray] = None,
     return_fracture_tips: bool = True,
-) -> Union[List[pp.Grid], Tuple[List[pp.Grid], np.ndarray]]:
+) -> Union[list[pp.Grid], tuple[list[pp.Grid], np.ndarray]]:
     """Create 1d grids for lines of a specified type from a gmsh tessalation.
 
     Only lines that were defined as 'physical' in the gmsh sense may have a grid
@@ -253,7 +255,7 @@ def create_1d_grids(
     # fractures), fracture tips, and auxiliary lines (to be disregarded)
 
     # Data structure for the point grids
-    g_1d: List[pp.Grid] = []
+    g_1d: list[pp.Grid] = []
 
     # If there are no fracture intersections, we return empty lists
     if "line" not in cells:
@@ -318,11 +320,11 @@ def create_1d_grids(
 
 def create_0d_grids(
     pts: np.ndarray,
-    cells: Dict[str, np.ndarray],
-    phys_names: Dict[int, str],
-    cell_info: Dict[str, np.ndarray],
-    target_tag_stem: str = None,
-) -> List[pp.PointGrid]:
+    cells: dict[str, np.ndarray],
+    phys_names: dict[int, str],
+    cell_info: dict[str, np.ndarray],
+    target_tag_stem: Optional[str] = None,
+) -> list[pp.PointGrid]:
     """Create 0d grids for points of a specified type from a gmsh tessalation.
 
     Only points that were defined as 'physical' in the gmsh sense may have a grid
diff --git a/src/porepy/fracs/structured.py b/src/porepy/fracs/structured.py
index d92c554b4..571f8f166 100644
--- a/src/porepy/fracs/structured.py
+++ b/src/porepy/fracs/structured.py
@@ -3,7 +3,9 @@
 
 The functions in this module can be accessed through the meshing wrapper module.
 """
-from typing import List
+from __future__ import annotations
+
+from typing import Optional
 
 import numpy as np
 import scipy.sparse as sps
@@ -15,8 +17,8 @@
 
 
 def _cart_grid_3d(
-    fracs: List[np.ndarray], nx: np.ndarray, physdims: np.ndarray = None
-) -> List[pp.Grid]:
+    fracs: list[np.ndarray], nx: np.ndarray, physdims: Optional[np.ndarray] = None
+) -> list[pp.Grid]:
     """
     Create grids for a domain with possibly intersecting fractures in 3d.
 
@@ -61,8 +63,8 @@ def _cart_grid_3d(
 
 
 def _tensor_grid_3d(
-    fracs: List[np.ndarray], x: np.ndarray, y: np.ndarray, z: np.ndarray
-) -> List[List[pp.Grid]]:
+    fracs: list[np.ndarray], x: np.ndarray, y: np.ndarray, z: np.ndarray
+) -> list[list[pp.Grid]]:
     """
     Create a grids for a domain with possibly intersecting fractures in 3d.
 
@@ -140,8 +142,8 @@ def _cart_grid_2d(fracs, nx, physdims=None):
 
 
 def _tensor_grid_2d(
-    fracs: List[np.ndarray], x: np.ndarray, y: np.ndarray
-) -> List[List[pp.Grid]]:
+    fracs: list[np.ndarray], x: np.ndarray, y: np.ndarray
+) -> list[list[pp.Grid]]:
     """
     Create a grids for a domain with possibly intersecting fractures in 2d.
 
diff --git a/src/porepy/fracs/wells_3d.py b/src/porepy/fracs/wells_3d.py
index c4f877a41..59e616de2 100644
--- a/src/porepy/fracs/wells_3d.py
+++ b/src/porepy/fracs/wells_3d.py
@@ -9,8 +9,10 @@
     compute_well_fracture_intersections(well_network, fracture_network)
     well_network.mesh(mdg)
 """
+from __future__ import annotations
+
 import logging
-from typing import Dict, Iterator, List, Optional, Tuple
+from typing import Iterator, Optional
 
 import numpy as np
 import scipy.sparse as sps
@@ -40,7 +42,7 @@ def __init__(
         self,
         points: np.ndarray,
         index: Optional[int] = None,
-        tags: Optional[Dict] = None,
+        tags: Optional[dict] = None,
     ) -> None:
         """Initialize fractures.
 
@@ -76,7 +78,7 @@ def set_index(self, i: int) -> None:
         """
         self.index = i
 
-    def segments(self) -> Iterator[Tuple[List[int], np.ndarray]]:
+    def segments(self) -> Iterator[tuple[list[int], np.ndarray]]:
         """
         Iterate over the segments defined through segment indices and endpoints.
         """
@@ -91,7 +93,7 @@ def num_points(self) -> int:
     def num_segments(self) -> int:
         return self.num_points() - 1
 
-    def add_point(self, point: np.ndarray, ind: int = None) -> None:
+    def add_point(self, point: np.ndarray, ind: Optional[int] = None) -> None:
         """Add new pts (3 x 1) to self.pts.
         If ind is not specified, the point is appended at the end of
         self.pts. Otherwise, it is inserted between the old points number ind
@@ -102,7 +104,7 @@ def add_point(self, point: np.ndarray, ind: int = None) -> None:
         else:
             self.pts = np.hstack((self.pts[:, :ind], point, self.pts[:, ind:]))
 
-    def _mesh_size(self, segment_ind: Tuple[int] = None) -> Optional[float]:
+    def _mesh_size(self, segment_ind: Optional[tuple[int]] = None) -> Optional[float]:
         """Return the mesh size for a well or one of its segments.
 
         Args:
@@ -146,10 +148,10 @@ class WellNetwork3d:
 
     def __init__(
         self,
-        wells: Optional[List[Well]] = None,
-        domain: Optional[Dict[str, float]] = None,
+        wells: Optional[list[Well]] = None,
+        domain: Optional[dict[str, float]] = None,
         tol: float = 1e-8,
-        parameters: Optional[Dict] = None,
+        parameters: Optional[dict] = None,
     ) -> None:
         """Initialize well network.
 
@@ -183,12 +185,12 @@ def __init__(
             self.parameters = parameters
 
         if domain is None:
-            domain = {}
-        self.domain: Dict = domain
+            domain = dict()
+        self.domain: dict = domain
         self.tol = tol
 
         # Assign an empty tag dictionary
-        self.tags: Dict[str, List[bool]] = {}
+        self.tags: dict[str, list[bool]] = dict()
 
     def add(self, w: Well) -> None:
         """Add a well to the network.
@@ -208,7 +210,7 @@ def add(self, w: Well) -> None:
             w.set_index(0)
         self._wells.append(w)
 
-    def _mesh_size(self, w: Well, segment_ind: Tuple[int] = None) -> float:
+    def _mesh_size(self, w: Well, segment_ind: Optional[tuple[int]] = None) -> float:
         """Return the mesh size for a well or one of its segments.
 
         Args:
@@ -442,7 +444,7 @@ def compute_well_fracture_intersections(
 
 def compute_well_rock_matrix_intersections(
     mdg: pp.MixedDimensionalGrid,
-    cells: np.ndarray = None,
+    cells: Optional[np.ndarray] = None,
     min_length: float = 1e-10,
     tol: float = 1e-5,
 ) -> None:
@@ -470,7 +472,7 @@ def compute_well_rock_matrix_intersections(
     tree.from_grid(sd_max, cells)
 
     # Extract the grids of the wells of co-dimension 2
-    well_subdomains: List[pp.Grid] = [
+    well_subdomains: list[pp.Grid] = [
         g for g in mdg.subdomains(dim=dim_max - 2) if hasattr(g, "well_num")
     ]
 
@@ -562,7 +564,7 @@ def compute_well_rock_matrix_intersections(
 
 def _argsort_points_along_line_segment(
     seg: np.ndarray,
-) -> Tuple[np.ndarray, np.ndarray]:
+) -> tuple[np.ndarray, np.ndarray]:
     """Sort point lying along a segment.
 
     Args:
@@ -594,10 +596,10 @@ def _argsort_points_along_line_segment(
 def _intersection_segment_fracture(
     segment_points: np.ndarray,
     fracture: pp.PlaneFracture,
-    tags: List[np.ndarray],
+    tags: list[np.ndarray],
     ignore_endpoint_tag: bool,
     tol: float = 1e-8,
-) -> Tuple[np.ndarray, List[np.ndarray]]:
+) -> tuple[np.ndarray, list[np.ndarray]]:
     """Compute intersection between a single line segment and fracture.
 
     Computes intersection point. If no intersection exists (distance > 0), no updates are
diff --git a/src/porepy/grids/md_grid.py b/src/porepy/grids/md_grid.py
index 501703e74..14a78f4e5 100644
--- a/src/porepy/grids/md_grid.py
+++ b/src/porepy/grids/md_grid.py
@@ -589,7 +589,7 @@ def replace_subdomains_and_interfaces(
     # ---- Methods for getting information on the bucket, or its components ----
 
     def diameter(
-        self, cond: Callable[[Union[pp.Grid, pp.MortarGrid]], bool] = None
+        self, cond: Optional[Callable[[Union[pp.Grid, pp.MortarGrid]], bool]] = None
     ) -> float:
         """Compute the cell diameter (mesh size) of the mixed-dimensional grid.
 
@@ -632,7 +632,9 @@ def dim_max(self) -> int:
         # lead to a recursion error.
         return np.max([sd.dim for sd in self._subdomain_data.keys()])
 
-    def num_subdomain_cells(self, cond: Callable[[pp.Grid], bool] = None) -> int:
+    def num_subdomain_cells(
+        self, cond: Optional[Callable[[pp.Grid], bool]] = None
+    ) -> int:
         """Compute the total number of cells of the mixed-dimensional grid.
 
         A function can be passed to filter subdomains and/or interfaces.
@@ -650,7 +652,9 @@ def num_subdomain_cells(self, cond: Callable[[pp.Grid], bool] = None) -> int:
             [grid.num_cells for grid in self.subdomains() if cond(grid)], dtype=int
         )
 
-    def num_interface_cells(self, cond: Callable[[pp.MortarGrid], bool] = None) -> int:
+    def num_interface_cells(
+        self, cond: Optional[Callable[[pp.MortarGrid], bool]] = None
+    ) -> int:
         """
         Compute the total number of mortar cells of the mixed-dimensional grid.
 
diff --git a/src/porepy/grids/mortar_grid.py b/src/porepy/grids/mortar_grid.py
index ca87f259e..34853f49a 100644
--- a/src/porepy/grids/mortar_grid.py
+++ b/src/porepy/grids/mortar_grid.py
@@ -212,7 +212,7 @@ def compute_geometry(self) -> None:
     ### Methods to update the mortar grid, or the neighboring grids.
 
     def update_mortar(
-        self, new_side_grids: dict[MortarSides, pp.Grid], tol: float = None
+        self, new_side_grids: dict[MortarSides, pp.Grid], tol: Optional[float] = None
     ) -> None:
         """
         Update the low_to_mortar_int and high_to_mortar_int maps when the mortar grids
@@ -319,7 +319,7 @@ def update_mortar(
 
         self._check_mappings()
 
-    def update_secondary(self, new_g: pp.Grid, tol: float = None) -> None:
+    def update_secondary(self, new_g: pp.Grid, tol: Optional[float] = None) -> None:
         """Update the mappings between Mortar and secondary grid when the latter is changed.
 
         NOTE: This function assumes that the secondary grid is only updated once: A change
@@ -395,7 +395,9 @@ def update_secondary(self, new_g: pp.Grid, tol: float = None) -> None:
 
         self._check_mappings()
 
-    def update_primary(self, g_new: pp.Grid, g_old: pp.Grid, tol: float = None) -> None:
+    def update_primary(
+        self, g_new: pp.Grid, g_old: pp.Grid, tol: Optional[float] = None
+    ) -> None:
         """
 
         Update the _primary_to_mortar_int map when the primary (higher-dimensional) grid is
diff --git a/src/porepy/grids/partition.py b/src/porepy/grids/partition.py
index 2e4fa71c1..4bd6db830 100644
--- a/src/porepy/grids/partition.py
+++ b/src/porepy/grids/partition.py
@@ -7,6 +7,7 @@
 from __future__ import annotations
 
 import warnings
+from typing import Optional
 
 import numpy as np
 import scipy.sparse as sps
@@ -62,7 +63,7 @@ def partition_metis(g: pp.Grid, num_part: int) -> np.ndarray:
 
 
 def partition_structured(
-    g: pp.TensorGrid, num_part: int = 1, coarse_dims: np.ndarray = None
+    g: pp.TensorGrid, num_part: int = 1, coarse_dims: Optional[np.ndarray] = None
 ) -> np.ndarray:
     """Define a partitioning of a grid based on logical Cartesian indexing.
 
@@ -767,7 +768,7 @@ def overlap(
 
 
 def grid_is_connected(
-    g: pp.Grid, cell_ind: np.ndarray = None
+    g: pp.Grid, cell_ind: Optional[np.ndarray] = None
 ) -> tuple[bool, list[np.ndarray]]:
     """
     Check if a grid is fully connected, as defined by its cell_connection_map().
diff --git a/src/porepy/grids/point_grid.py b/src/porepy/grids/point_grid.py
index 15196d7c4..d8935cd0b 100644
--- a/src/porepy/grids/point_grid.py
+++ b/src/porepy/grids/point_grid.py
@@ -1,3 +1,5 @@
+from typing import Optional
+
 import numpy as np
 import scipy.sparse as sps
 
@@ -5,7 +7,7 @@
 
 
 class PointGrid(Grid):
-    def __init__(self, pt: np.ndarray, name: str = None) -> None:
+    def __init__(self, pt: np.ndarray, name: Optional[str] = None) -> None:
         """Constructor for 0D grid.
 
         Args:
diff --git a/src/porepy/grids/standard_grids/md_grids_2d.py b/src/porepy/grids/standard_grids/md_grids_2d.py
index 39fb33bef..c1c540999 100644
--- a/src/porepy/grids/standard_grids/md_grids_2d.py
+++ b/src/porepy/grids/standard_grids/md_grids_2d.py
@@ -12,7 +12,7 @@
 """
 from __future__ import annotations
 
-from typing import Union
+from typing import Optional, Union
 
 import numpy as np
 
@@ -31,7 +31,11 @@ def _n_cells(mesh_args: Union[np.ndarray, dict, None]) -> np.ndarray:
         return mesh_args
 
 
-def single_horizontal(mesh_args=None, x_endpoints=None, simplex=True):
+def single_horizontal(
+    mesh_args: Optional[dict | np.ndarray] = None,
+    x_endpoints: Optional[np.ndarray] = None,
+    simplex: bool = True,
+):
     """
     Create a grid bucket for a domain containing a single horizontal fracture at y=0.5.
 
@@ -39,7 +43,7 @@ def single_horizontal(mesh_args=None, x_endpoints=None, simplex=True):
         mesh_args:  For triangular grids: Dictionary containing at least "mesh_size_frac". If
                         the optional values of "mesh_size_bound" and "mesh_size_min" are
                         not provided, these are set by utils.set_mesh_sizes.
-                    For cartesian grids: List containing number of cells in x and y
+                    For cartesian grids: np.array containing number of cells in x and y
                         direction.
         x_endpoints (list): Contains the x coordinates of the two endpoints. If not
             provided, the endpoints will be set to [0, 1]
@@ -55,6 +59,11 @@ def single_horizontal(mesh_args=None, x_endpoints=None, simplex=True):
     if simplex:
         if mesh_args is None:
             mesh_args = {"mesh_size_frac": 0.2}
+
+        if not isinstance(mesh_args, dict):
+            # The numpy-array format for mesh_args should not be used for simplex grids
+            raise ValueError("Mesh arguments should be a dictionary for simplex grids")
+
         points = np.array([x_endpoints, [0.5, 0.5]])
         edges = np.array([[0], [1]])
         mdg = utils.make_mdg_2d_simplex(mesh_args, points, edges, domain)
@@ -66,9 +75,9 @@ def single_horizontal(mesh_args=None, x_endpoints=None, simplex=True):
 
 
 def single_vertical(
-    mesh_args: Union[np.ndarray, dict] = None,
-    y_endpoints: np.ndarray = None,
-    simplex: bool = True,
+    mesh_args: Optional[Union[np.ndarray, dict]] = None,
+    y_endpoints: Optional[np.ndarray] = None,
+    simplex: Optional[bool] = True,
 ):
     """
     Create a grid bucket for a domain containing a single vertical fracture at x=0.5.
@@ -97,7 +106,11 @@ def single_vertical(
     if simplex:
         if mesh_args is None:
             mesh_args = {"mesh_size_frac": 0.2}
-        assert type(mesh_args) == dict
+
+        if not isinstance(mesh_args, dict):
+            # The numpy-array format for mesh_args should not be used for simplex grids
+            raise ValueError("Mesh arguments should be a dictionary for simplex grids")
+
         points = np.array([[0.5, 0.5], y_endpoints])
         edges = np.array([[0], [1]])
         mdg = utils.make_mdg_2d_simplex(mesh_args, points, edges, domain)
@@ -108,7 +121,12 @@ def single_vertical(
     return mdg, domain
 
 
-def two_intersecting(mesh_args=None, x_endpoints=None, y_endpoints=None, simplex=True):
+def two_intersecting(
+    mesh_args: Optional[dict | np.ndarray] = None,
+    x_endpoints: Optional[list] = None,
+    y_endpoints: Optional[list] = None,
+    simplex: bool = True,
+):
     """
     Create a grid bucket for a domain containing fractures, one horizontal and one vertical
     at y=0.5 and x=0.5 respectively.
@@ -138,6 +156,11 @@ def two_intersecting(mesh_args=None, x_endpoints=None, y_endpoints=None, simplex
     if simplex:
         if mesh_args is None:
             mesh_args = {"mesh_size_frac": 0.2}
+
+        if not isinstance(mesh_args, dict):
+            # The numpy-array format for mesh_args should not be used for simplex grids
+            raise ValueError("Mesh arguments should be a dictionary for simplex grids")
+
         points = np.array(
             [
                 [x_endpoints[0], x_endpoints[1], 0.5, 0.5],
@@ -159,7 +182,7 @@ def two_intersecting(mesh_args=None, x_endpoints=None, y_endpoints=None, simplex
     return mdg, domain
 
 
-def seven_fractures_one_L_intersection(mesh_args=None):
+def seven_fractures_one_L_intersection(mesh_args: dict):
     """
     Create a grid bucket for a domain containing the network introduced as example 1 of
     Berge et al. 2019: Finite volume discretization for poroelastic media with fractures
@@ -197,7 +220,7 @@ def seven_fractures_one_L_intersection(mesh_args=None):
     return mdg, domain
 
 
-def benchmark_regular(mesh_args, is_coarse=False):
+def benchmark_regular(mesh_args: dict, is_coarse: bool = False):
     """
     Create a grid bucket for a domain containing the network introduced as example 2 of
     Berre et al. 2018: Benchmarks for single-phase flow in fractured porous media.
diff --git a/src/porepy/grids/standard_grids/utils.py b/src/porepy/grids/standard_grids/utils.py
index c903ef06c..b98030b59 100644
--- a/src/porepy/grids/standard_grids/utils.py
+++ b/src/porepy/grids/standard_grids/utils.py
@@ -1,4 +1,4 @@
-from typing import List
+from __future__ import annotations
 
 import numpy as np
 
@@ -65,7 +65,7 @@ def make_mdg_2d_simplex(
 
 
 def make_mdg_2d_cartesian(
-    n_cells: np.ndarray, fractures: List[np.ndarray], domain: dict
+    n_cells: np.ndarray, fractures: list[np.ndarray], domain: dict
 ):
     """
     Construct a grid bucket with Cartesian grid in the 2d domain.
diff --git a/src/porepy/models/abstract_model.py b/src/porepy/models/abstract_model.py
index 8b1992231..fe057d969 100644
--- a/src/porepy/models/abstract_model.py
+++ b/src/porepy/models/abstract_model.py
@@ -18,7 +18,8 @@
 import logging
 import time
 import warnings
-from typing import Any, Dict, Optional, Tuple
+from collections import namedtuple
+from typing import Any, Dict, NamedTuple, Optional, Tuple
 
 import numpy as np
 import scipy.sparse as sps
@@ -345,35 +346,65 @@ def _export(self) -> None:
         """
         pass
 
-    def _domain_boundary_sides(
-        self, g: pp.Grid
-    ) -> Tuple[
-        np.ndarray,
-        np.ndarray,
-        np.ndarray,
-        np.ndarray,
-        np.ndarray,
-        np.ndarray,
-        np.ndarray,
-    ]:
-        """Obtain indices of the faces of a grid that lie on each side of the domain
-        boundaries. It is assumed the domain is box shaped.
+    def _domain_boundary_sides(self, sd: pp.Grid, tol: float = 1e-10) -> NamedTuple:
+        """Obtain indices of the faces lying on the sides of the domain boundaries.
+
+        The method is primarily intended for box-shaped domains. However, it can also be
+        applied to non-box-shaped domains (e.g., domains with perturbed boundary nodes)
+        provided `tol` is tuned accordingly.
+
+        Parameters:
+            sd: Subdomain grid.
+            tol: Tolerance used to determine whether a face center lies on a boundary side.
+
+        Returns:
+            NamedTuple containing the domain boundary sides. Available attributes are:
+
+                - all_bf (np.ndarray of int): indices of the boundary faces.
+                - east (np.ndarray of bool): flags of the faces lying on the East side.
+                - west (np.ndarray of bool): flags of the faces lying on the West side.
+                - north (np.ndarray of bool): flags of the faces lying on the North side.
+                - south (np.ndarray of bool): flags of the faces lying on the South side.
+                - top (np.ndarray of bool): flags of the faces lying on the Top side.
+                - bottom (np.ndarray of bool): flags of the faces lying on Bottom side.
+
+        Examples:
+
+            .. code:: python
+
+                model = pp.IncompressibleFlow({})
+                model.prepare_simulation()
+                sd = model.mdg.subdomains()[0]
+                sides = model._domain_boundary_sides(sd)
+                # Access north faces using index or name is equivalent:
+                north_by_index = sides[3]
+                north_by_name = sides.north
+                assert all(north_by_index == north_by_name)
+
         """
-        tol = 1e-10
+
+        # Get domain boundary sides
         box = self.box
-        east = g.face_centers[0] > box["xmax"] - tol
-        west = g.face_centers[0] < box["xmin"] + tol
+        east = np.abs(box["xmax"] - sd.face_centers[0]) <= tol
+        west = np.abs(box["xmin"] - sd.face_centers[0]) <= tol
         if self.mdg.dim_max() == 1:
-            north = np.zeros(g.num_faces, dtype=bool)
+            north = np.zeros(sd.num_faces, dtype=bool)
             south = north.copy()
         else:
-            north = g.face_centers[1] > box["ymax"] - tol
-            south = g.face_centers[1] < box["ymin"] + tol
+            north = np.abs(box["ymax"] - sd.face_centers[1]) <= tol
+            south = np.abs(box["ymin"] - sd.face_centers[1]) <= tol
         if self.mdg.dim_max() < 3:
-            top = np.zeros(g.num_faces, dtype=bool)
+            top = np.zeros(sd.num_faces, dtype=bool)
             bottom = top.copy()
         else:
-            top = g.face_centers[2] > box["zmax"] - tol
-            bottom = g.face_centers[2] < box["zmin"] + tol
-        all_bf = g.get_boundary_faces()
-        return all_bf, east, west, north, south, top, bottom
+            top = np.abs(box["zmax"] - sd.face_centers[2]) <= tol
+            bottom = np.abs(box["zmin"] - sd.face_centers[2]) <= tol
+        all_bf = sd.get_boundary_faces()
+
+        # Create a namedtuple to store the arrays
+        DomainSides = namedtuple(
+            "DomainSides", "all_bf east west north south top bottom"
+        )
+        domain_sides = DomainSides(all_bf, east, west, north, south, top, bottom)
+
+        return domain_sides
diff --git a/src/porepy/numerics/ad/equation_manager.py b/src/porepy/numerics/ad/equation_manager.py
index e84f08aa4..1c2d43761 100644
--- a/src/porepy/numerics/ad/equation_manager.py
+++ b/src/porepy/numerics/ad/equation_manager.py
@@ -527,7 +527,9 @@ def _column_projection(self, variables: Sequence["pp.ad.Variable"]) -> sps.spmat
 
     def _variable_set_complement(
         self,
-        variables: Sequence[Union["pp.ad.Variable", "pp.ad.MergedVariable"]] = None,
+        variables: Optional[
+            Sequence[Union["pp.ad.Variable", "pp.ad.MergedVariable"]]
+        ] = None,
     ) -> List["pp.ad.Variable"]:
         """
         Take the complement of a set of variables, with respect to the full set of
diff --git a/src/porepy/numerics/ad/functions.py b/src/porepy/numerics/ad/functions.py
index d2c6b0163..f8c149906 100644
--- a/src/porepy/numerics/ad/functions.py
+++ b/src/porepy/numerics/ad/functions.py
@@ -262,8 +262,15 @@ def heaviside(var, zerovalue: float = 0.5):
 
 
 def heaviside_smooth(var, eps: float = 1e-3):
-    """
-    Smooth (regularized) version of the Heaviside function.
+    """Smooth (regularized) version of the Heaviside function.
+
+    Note:
+        The analytical expression for the smooth version Heaviside function reads:
+            ``H_eps(x) = (1/2) * (1 + (2/pi) * arctan(x/eps))``,
+        with its derivative smoothly approximating the Dirac delta function:
+            ``d(H(x))/dx = delta_eps = (1/pi) * (eps / (eps^2 + x^2))``.
+
+        Reference: https://ieeexplore.ieee.org/document/902291
 
     Parameters:
         var: Input array.
@@ -274,14 +281,6 @@ def heaviside_smooth(var, eps: float = 1e-3):
         Regularized heaviside function (and its Jacobian if applicable) in form of a
         Ad_array or ndarray (depending on the input).
 
-    Note:
-        The analytical expression for the smooth version Heaviside function reads:
-            ``H_eps(x) = (1/2) * (1 + (2/pi) * arctan(x/eps))``,
-        with its derivative smoothly approximating the Dirac delta function:
-            ``d(H(x))/dx = delta_eps = (1/pi) * (eps / (eps^2 + x^2))``.
-
-        Reference: https://ieeexplore.ieee.org/document/902291
-
     """
     if isinstance(var, Ad_array):
         val = 0.5 * (1 + 2 * np.pi ** (-1) * np.arctan(var.val * eps ** (-1)))
@@ -314,7 +313,7 @@ def maximum(var0: pp.ad.Ad_array, var1: pp.ad.Ad_array | np.ndarray) -> pp.ad.Ad
     wrapped in a pp.ad.Array, whereas the first argument is expected to be an
     Ad_array originating from a pp.ad.Operator.
 
-    Args:
+    Parameters:
         var0: First argument to the maximum function.
         var1: Second argument.
 
@@ -345,10 +344,10 @@ def characteristic_function(tol: float, var: pp.ad.Ad_array):
     Returns 1 if ``var.val`` is within absolute tolerance = ``tol`` of zero.
     The derivative is set to zero independent of ``var.val``.
 
-    Notes:
+    Note:
         See module level documentation on how to wrap functions like this in ``ad.Function``.
 
-    Args:
+    Parameters:
         tol: Absolute tolerance for comparison with 0 using np.isclose.
         var: Ad operator (variable or expression).
 
diff --git a/src/porepy/numerics/ad/operator_functions.py b/src/porepy/numerics/ad/operator_functions.py
index 41da85cca..3c38ca5b2 100644
--- a/src/porepy/numerics/ad/operator_functions.py
+++ b/src/porepy/numerics/ad/operator_functions.py
@@ -8,7 +8,7 @@
 
 import abc
 from functools import partial
-from typing import Callable, List, Optional, Type, Union
+from typing import Callable, Optional, Type, Union
 
 import numpy as np
 import scipy.sparse as sps
@@ -88,7 +88,7 @@ def __init__(
         self._operation: Operation = Operation.evaluate
         self._set_tree()
 
-    def __call__(self, *args: pp.ad.Operator | Ad_array):
+    def __call__(self, *args: pp.ad.Operator | Ad_array) -> pp.ad.Operator:
         """Renders this function operator callable, fulfilling its notion as 'function'.
 
         Parameters:
@@ -321,7 +321,7 @@ def __init__(
         self,
         func: Callable,
         name: str,
-        multipliers: Union[List[float], float],
+        multipliers: float | list[float],
         array_compatible: bool = False,
     ):
         super().__init__(func, name, array_compatible)
@@ -586,7 +586,7 @@ def __get__(self, binding_instance: object, binding_type: type) -> Callable:
         # This will trigger the function evaluation.
         return partial(self.__call__, binding_instance)
 
-    def ad_wrapper(self, *args, **kwargs) -> AbstractFunction:
+    def ad_wrapper(self, *args, **kwargs) -> Operator:
         """Actual wrapper function.
         Constructs the necessary AD-Operator class wrapping the decorated callable
         and performs the evaluation/call.
diff --git a/src/porepy/numerics/ad/operators.py b/src/porepy/numerics/ad/operators.py
index 7bfda11cd..3f96b5fbd 100644
--- a/src/porepy/numerics/ad/operators.py
+++ b/src/porepy/numerics/ad/operators.py
@@ -254,9 +254,6 @@ def discretize(self, mdg: pp.MixedDimensionalGrid) -> None:
         """Perform the ``discretize`` function of all child operators which are discretizations
         using data from mdg.
 
-        IMPLEMENTATION NOTE: The discretizations was identified at initialization of
-        Expression - it is now done here to accommodate updates (?) and
-
         """
         unique_discretizations: dict[
             _ad_utils.MergedOperator, GridLike
@@ -742,7 +739,7 @@ def evaluate(
         self,
         dof_manager: pp.DofManager,
         state: Optional[np.ndarray] = None,
-    ) -> Ad_array:
+    ):
         """Evaluate the residual and Jacobian matrix for a given state.
 
         Parameters:
diff --git a/src/porepy/numerics/fv/fvutils.py b/src/porepy/numerics/fv/fvutils.py
index 815a207d4..4bd0ba587 100644
--- a/src/porepy/numerics/fv/fvutils.py
+++ b/src/porepy/numerics/fv/fvutils.py
@@ -6,7 +6,9 @@
 between these methods, the current structure with multiple auxiliary methods emerged.
 
 """
-from typing import Any, Callable, Dict, Generator, List, Optional, Tuple
+from __future__ import annotations
+
+from typing import Any, Callable, Generator, Optional
 
 import numpy as np
 import scipy.sparse as sps
@@ -299,8 +301,8 @@ def determine_eta(sd: pp.Grid) -> float:
 
 
 def find_active_indices(
-    parameter_dictionary: Dict[str, Any], sd: pp.Grid
-) -> Tuple[np.ndarray, np.ndarray]:
+    parameter_dictionary: dict[str, Any], sd: pp.Grid
+) -> tuple[np.ndarray, np.ndarray]:
     """Process information in parameter dictionary on whether the discretization
     should consider a subgrid. Look for fields in the parameter dictionary called
     specified_cells, specified_faces or specified_nodes. These are then processed by
@@ -351,7 +353,7 @@ def find_active_indices(
 def subproblems(
     sd: pp.Grid, max_memory: int, peak_memory_estimate: int
 ) -> Generator[
-    Tuple[pp.Grid, np.ndarray, np.ndarray, np.ndarray, np.ndarray], None, None
+    tuple[pp.Grid, np.ndarray, np.ndarray, np.ndarray, np.ndarray], None, None
 ]:
 
     if sd.dim == 0:
@@ -614,7 +616,7 @@ def vector_divergence(sd: pp.Grid) -> sps.csr_matrix:
 
 def scalar_tensor_vector_prod(
     sd: pp.Grid, k: pp.SecondOrderTensor, subcell_topology: SubcellTopology
-) -> Tuple[sps.csr_matrix, np.ndarray, np.ndarray]:
+) -> tuple[sps.csr_matrix, np.ndarray, np.ndarray]:
     """
     Compute product of normal vectors and tensors on a sub-cell level.
     This is essentially defining Darcy's law for each sub-face in terms of
@@ -1002,18 +1004,18 @@ def keep_neumann(self, other, transform=True):
 
 def partial_update_discretization(
     sd: pp.Grid,  # Grid
-    data: Dict,  # full data dictionary for this grid
+    data: dict,  # full data dictionary for this grid
     keyword: str,  # keyword for the target discretization
     discretize: Callable,  # Discretization operation
     dim: Optional[int] = None,  # dimension. Used to expand vector quantities
-    scalar_cell_right: Optional[List[str]] = None,  # See method documentation
-    vector_cell_right: Optional[List[str]] = None,
-    scalar_face_right: Optional[List[str]] = None,
-    vector_face_right: Optional[List[str]] = None,
-    scalar_cell_left: Optional[List[str]] = None,
-    vector_cell_left: Optional[List[str]] = None,
-    scalar_face_left: Optional[List[str]] = None,
-    vector_face_left: Optional[List[str]] = None,
+    scalar_cell_right: Optional[list[str]] = None,  # See method documentation
+    vector_cell_right: Optional[list[str]] = None,
+    scalar_face_right: Optional[list[str]] = None,
+    vector_face_right: Optional[list[str]] = None,
+    scalar_cell_left: Optional[list[str]] = None,
+    vector_cell_left: Optional[list[str]] = None,
+    scalar_face_left: Optional[list[str]] = None,
+    vector_face_left: Optional[list[str]] = None,
     second_keyword: Optional[str] = None,  # Used for biot discertization
 ) -> None:
     """Do partial update of discretization scheme.
@@ -1191,10 +1193,10 @@ def partial_update_discretization(
 
 def cell_ind_for_partial_update(
     sd: pp.Grid,
-    cells: np.ndarray = None,
-    faces: np.ndarray = None,
-    nodes: np.ndarray = None,
-) -> Tuple[np.ndarray, np.ndarray]:
+    cells: Optional[np.ndarray] = None,
+    faces: Optional[np.ndarray] = None,
+    nodes: Optional[np.ndarray] = None,
+) -> tuple[np.ndarray, np.ndarray]:
     """Obtain indices of cells and faces needed for a partial update of the
     discretization stencil.
 
@@ -1400,7 +1402,7 @@ def map_subgrid_to_grid(
     loc_cells: np.ndarray,
     is_vector: bool,
     nd: Optional[int] = None,
-) -> Tuple[np.ndarray, np.ndarray]:
+) -> tuple[np.ndarray, np.ndarray]:
     """Obtain mappings from the cells and faces of a subgrid back to a larger grid.
 
     Args:
diff --git a/src/porepy/numerics/fv/mpsa.py b/src/porepy/numerics/fv/mpsa.py
index f6e3a6987..0d3d497d5 100644
--- a/src/porepy/numerics/fv/mpsa.py
+++ b/src/porepy/numerics/fv/mpsa.py
@@ -11,7 +11,7 @@
 
 import logging
 from time import time
-from typing import Any, Dict, Tuple
+from typing import Any, Optional
 
 import numpy as np
 import scipy.sparse as sps
@@ -84,7 +84,7 @@ def ndof(self, sd: pp.Grid) -> int:
         return sd.dim * sd.num_cells
 
     def extract_displacement(
-        self, sd: pp.Grid, solution_array: np.ndarray, d: Dict
+        self, sd: pp.Grid, solution_array: np.ndarray, d: dict
     ) -> np.ndarray:
         """Extract the displacement part of a solution.
 
@@ -102,7 +102,7 @@ def extract_displacement(
         return solution_array
 
     def extract_stress(
-        self, sd: pp.Grid, solution_array: np.ndarray, d: Dict
+        self, sd: pp.Grid, solution_array: np.ndarray, d: dict
     ) -> np.ndarray:
         """Extract the stress corresponding to a solution
 
@@ -128,7 +128,7 @@ def extract_stress(
 
         return stress * solution_array + bound_stress * bc_val
 
-    def discretize(self, sd: pp.Grid, data: Dict) -> None:
+    def discretize(self, sd: pp.Grid, data: dict) -> None:
         """
         Discretize the second order vector elliptic equation using multi-point
         stress approximation.
@@ -175,8 +175,8 @@ def discretize(self, sd: pp.Grid, data: Dict) -> None:
         data (dict): For entries, see above.
 
         """
-        parameter_dictionary: Dict[str, Any] = data[pp.PARAMETERS][self.keyword]
-        matrix_dictionary: Dict[str, sps.spmatrix] = data[pp.DISCRETIZATION_MATRICES][
+        parameter_dictionary: dict[str, Any] = data[pp.PARAMETERS][self.keyword]
+        matrix_dictionary: dict[str, sps.spmatrix] = data[pp.DISCRETIZATION_MATRICES][
             self.keyword
         ]
         constit: pp.FourthOrderTensor = parameter_dictionary["fourth_order_tensor"]
@@ -415,8 +415,8 @@ def update_discretization(self, sd: pp.Grid, data: dict) -> None:
         )
 
     def assemble_matrix_rhs(
-        self, sd: pp.Grid, data: Dict
-    ) -> Tuple[sps.spmatrix, np.ndarray]:
+        self, sd: pp.Grid, data: dict
+    ) -> tuple[sps.spmatrix, np.ndarray]:
         """
         Return the matrix and right-hand side for a discretization of a second
         order elliptic equation using a FV method with a multi-point stress
@@ -441,7 +441,7 @@ def assemble_matrix_rhs(
         """
         return self.assemble_matrix(sd, data), self.assemble_rhs(sd, data)
 
-    def assemble_matrix(self, sd: pp.Grid, data: Dict) -> sps.spmatrix:
+    def assemble_matrix(self, sd: pp.Grid, data: dict) -> sps.spmatrix:
         """
         Return the matrix for a discretization of a second order elliptic vector
         equation using a FV method.
@@ -470,7 +470,7 @@ def assemble_matrix(self, sd: pp.Grid, data: Dict) -> sps.spmatrix:
 
         return M
 
-    def assemble_rhs(self, sd: pp.Grid, data: Dict) -> np.ndarray:
+    def assemble_rhs(self, sd: pp.Grid, data: dict) -> np.ndarray:
         """Return the right-hand side for a discretization of a second
         order elliptic equation using a finite volume method.
 
@@ -506,11 +506,11 @@ def _stress_disrcetization(
         sd: pp.Grid,
         constit: pp.FourthOrderTensor,
         bound: pp.BoundaryConditionVectorial,
-        eta: float = None,
+        eta: Optional[float] = None,
         inverter: str = "numba",
-        hf_disp: bool = False,
-        hf_eta: float = None,
-    ) -> Tuple[sps.spmatrix, sps.spmatrix, sps.spmatrix, sps.spmatrix]:
+        hf_disp: Optional[bool] = False,
+        hf_eta: Optional[float] = None,
+    ) -> tuple[sps.spmatrix, sps.spmatrix, sps.spmatrix, sps.spmatrix]:
         """
         Actual implementation of the MPSA W-method. To calculate the MPSA
         discretization on a grid, either call this method, or, to respect the
@@ -621,7 +621,7 @@ def _stress_disrcetization(
             k = pp.SecondOrderTensor(2 * constit.mu + constit.lmbda)
             params["second_order_tensor"] = k
 
-            d: Dict = {
+            d: dict = {
                 pp.PARAMETERS: {tpfa_key: params},
                 pp.DISCRETIZATION_MATRICES: {tpfa_key: {}},
             }
@@ -728,7 +728,7 @@ def _create_inverse_gradient_matrix(
         bound_exclusion: pp.fvutils.ExcludeBoundaries,
         eta: float,
         inverter: str,
-    ) -> Tuple[sps.spmatrix, sps.spmatrix, np.ndarray]:
+    ) -> tuple[sps.spmatrix, sps.spmatrix, np.ndarray]:
         """
         This is the function where the real discretization takes place. It contains
         the parts that are common for elasticity and poro-elasticity, and was thus
@@ -1084,8 +1084,8 @@ def _reconstruct_displacement(
         self,
         sd: pp.Grid,
         subcell_topology: pp.fvutils.SubcellTopology,
-        eta: float = None,
-    ) -> Tuple[sps.csr_matrix, sps.csr_matrix]:
+        eta: Optional[float] = None,
+    ) -> tuple[sps.csr_matrix, sps.csr_matrix]:
         """
         Function for reconstructing the displacement at the half faces given the
         local gradients. For a subcell Ks associated with cell K and node s, the
@@ -1212,7 +1212,7 @@ def _get_displacement_submatrices(
         eta: float,
         num_sub_cells: int,
         bound_exclusion: pp.fvutils.ExcludeBoundaries,
-    ) -> Tuple[sps.spmatrix, sps.spmatrix]:
+    ) -> tuple[sps.spmatrix, sps.spmatrix]:
         nd = sd.dim
         # Distance from cell centers to face centers, this will be the
         # contribution from gradient unknown to equations for displacement
@@ -1247,7 +1247,7 @@ def _get_displacement_submatrices_rob(
         eta: float,
         num_sub_cells: int,
         bound_exclusion: pp.fvutils.ExcludeBoundaries,
-    ) -> Tuple[sps.spmatrix, sps.spmatrix]:
+    ) -> tuple[sps.spmatrix, sps.spmatrix]:
         nd = sd.dim
         # Distance from cell centers to face centers, this will be the
         # contribution from gradient unknown to equations for displacement
@@ -1297,7 +1297,7 @@ def _get_displacement_submatrices_rob(
 
     def _split_stiffness_matrix(
         self, constit: pp.FourthOrderTensor
-    ) -> Tuple[np.ndarray, np.ndarray]:
+    ) -> tuple[np.ndarray, np.ndarray]:
         """
         Split the stiffness matrix into symmetric and asymetric part
 
@@ -1361,7 +1361,7 @@ def _tensor_vector_prod(
         sd: pp.Grid,
         constit: pp.FourthOrderTensor,
         subcell_topology: pp.fvutils.SubcellTopology,
-    ) -> Tuple[sps.spmatrix, sps.spmatrix, np.ndarray, np.ndarray]:
+    ) -> tuple[sps.spmatrix, sps.spmatrix, np.ndarray, np.ndarray]:
         """Compute product between stiffness tensor and face normals.
 
         The method splits the stiffness matrix into a symmetric and asymmetric
@@ -1540,7 +1540,7 @@ def _block_diagonal_structure(
         nno: np.ndarray,
         bound_exclusion: pp.fvutils.ExcludeBoundaries,
         nd: int,
-    ) -> Tuple[sps.spmatrix, sps.spmatrix, np.ndarray]:
+    ) -> tuple[sps.spmatrix, sps.spmatrix, np.ndarray]:
         """
         Define matrices to turn linear system into block-diagonal form.
 
@@ -1681,7 +1681,7 @@ def _rearange_columns_displacement_eqs(
         d_cont_cell: np.ndarray,
         num_sub_cells: int,
         nd: int,
-    ) -> Tuple[np.ndarray, np.ndarray]:
+    ) -> tuple[np.ndarray, np.ndarray]:
         """Transform columns of displacement balance from increasing cell
         ordering (first x-variables of all cells, then y) to increasing
         variables (first all variables of the first cells, then...)
@@ -1720,7 +1720,7 @@ def _rearange_columns_displacement_eqs(
         d_cont_cell = d_cont_cell[:, d_cont_cell_map]
         return d_cont_grad, d_cont_cell
 
-    def _row_major_to_col_major(self, shape: Tuple, nd: int, axis: int) -> sps.spmatrix:
+    def _row_major_to_col_major(self, shape: tuple, nd: int, axis: int) -> sps.spmatrix:
         """Transform columns of displacement balance from increasing cell
         ordering (first x-variables of all cells, then y) to increasing
         variables (first all variables of the first cells, then...)
@@ -1821,7 +1821,7 @@ def _eliminate_ncasym_neumann(
 
     def _reduce_grid_constit_2d(
         self, sd: pp.Grid, constit: pp.FourthOrderTensor
-    ) -> Tuple[pp.Grid, pp.FourthOrderTensor]:
+    ) -> tuple[pp.Grid, pp.FourthOrderTensor]:
         sd = sd.copy()
 
         (
diff --git a/src/porepy/numerics/interface_laws/elliptic_interface_laws.py b/src/porepy/numerics/interface_laws/elliptic_interface_laws.py
index 2f9e46ab5..a194be522 100644
--- a/src/porepy/numerics/interface_laws/elliptic_interface_laws.py
+++ b/src/porepy/numerics/interface_laws/elliptic_interface_laws.py
@@ -1104,8 +1104,10 @@ def assemble_matrix_rhs(
         sd_data_secondary: Dict,
         intf_data: Dict,
         matrix: sps.spmatrix,
-    ) -> Tuple[np.ndarray, np.ndarray]:
-        pass
+    ) -> Tuple[np.ndarray, np.ndarray]:  # type:ignore
+        raise NotImplementedError(
+            "Wells are not compatible with the Assembler framework"
+        )
 
     def __repr__(self) -> str:
         return f"Interface coupling of Well type, with keyword {self.keyword}"
diff --git a/src/porepy/numerics/linear_solvers.py b/src/porepy/numerics/linear_solvers.py
index 57acdbbf8..0a6704e55 100644
--- a/src/porepy/numerics/linear_solvers.py
+++ b/src/porepy/numerics/linear_solvers.py
@@ -5,7 +5,9 @@
 is just a wrapper around that, mostly for compliance with the nonlinear
 case, see numerics.nonlinear.nonlinear_solvers.
 """
-from typing import Dict, Tuple
+from __future__ import annotations
+
+from typing import Optional
 
 from porepy.models.abstract_model import AbstractModel
 
@@ -15,7 +17,7 @@ class LinearSolver:
     model class before and after solving the problem.
     """
 
-    def __init__(self, params: Dict = None) -> None:
+    def __init__(self, params: Optional[dict] = None) -> None:
         """Define linear solver.
 
         Parameters:
@@ -29,7 +31,7 @@ def __init__(self, params: Dict = None) -> None:
         # default_options.update(params)
         self.params = params  # default_options
 
-    def solve(self, setup: AbstractModel) -> Tuple[float, bool]:
+    def solve(self, setup: AbstractModel) -> tuple[float, bool]:
         """Solve a linear problem defined by the current state of the model.
 
         Parameters:
diff --git a/src/porepy/numerics/mixed_dim/assembler.py b/src/porepy/numerics/mixed_dim/assembler.py
index d45157674..378a37780 100644
--- a/src/porepy/numerics/mixed_dim/assembler.py
+++ b/src/porepy/numerics/mixed_dim/assembler.py
@@ -88,14 +88,16 @@ def __init__(
         self._identify_variable_combinations()
 
     @staticmethod
-    def _discretization_key(row: str, col: str = None) -> str:
+    def _discretization_key(row: str, col: Optional[str] = None) -> str:
         if col is None or row == col:
             return row
         else:
             return row + "_" + col
 
     @staticmethod
-    def _variable_term_key(term: str, key_1: str, key_2: str, key_3: str = None) -> str:
+    def _variable_term_key(
+        term: str, key_1: str, key_2: str, key_3: Optional[str] = None
+    ) -> str:
         """Get the key-variable combination used to identify a specific term in the
         equation.
 
diff --git a/src/porepy/numerics/vem/dual_elliptic.py b/src/porepy/numerics/vem/dual_elliptic.py
index a05759a7e..e0c77826d 100644
--- a/src/porepy/numerics/vem/dual_elliptic.py
+++ b/src/porepy/numerics/vem/dual_elliptic.py
@@ -322,7 +322,7 @@ def _assemble_neumann_common(
         data: dict,
         M: sps.csr_matrix,
         mass: sps.csr_matrix,
-        bc_weight: float = None,
+        bc_weight: Optional[float] = None,
     ) -> tuple[sps.csr_matrix, np.ndarray]:
         """Impose Neumann boundary discretization on an already assembled
         system matrix.
diff --git a/src/porepy/params/data.py b/src/porepy/params/data.py
index c304a49b8..a9815e6b6 100644
--- a/src/porepy/params/data.py
+++ b/src/porepy/params/data.py
@@ -91,7 +91,7 @@ class Parameters(Dict):
     def __init__(
         self,
         grid: Optional[Union[pp.Grid, pp.MortarGrid]] = None,
-        keywords: list[str] = None,
+        keywords: Optional[list[str]] = None,
         dictionaries: Optional[list[dict]] = None,
     ):
         """Initialize Data object.
diff --git a/src/porepy/py.typed b/src/porepy/py.typed
new file mode 100644
index 000000000..e69de29bb