From a78a7fd51eb8c89bb95fc0ca1003387f43e97945 Mon Sep 17 00:00:00 2001
From: Jhabriel Varela <jhabriel@gmail.com>
Date: Tue, 8 Nov 2022 23:58:34 +0100
Subject: [PATCH 01/14] ENH: Allow passing `tol` in `domain_boundary_sides()`.

---
 src/porepy/models/abstract_model.py | 42 +++++++++++++++++++----------
 1 file changed, 28 insertions(+), 14 deletions(-)

diff --git a/src/porepy/models/abstract_model.py b/src/porepy/models/abstract_model.py
index 8b1992231..61be3e653 100644
--- a/src/porepy/models/abstract_model.py
+++ b/src/porepy/models/abstract_model.py
@@ -346,8 +346,10 @@ def _export(self) -> None:
         pass
 
     def _domain_boundary_sides(
-        self, g: pp.Grid
-    ) -> Tuple[
+        self,
+        sd: pp.Grid,
+        tol: float = 1e-10,
+    ) -> tuple[
         np.ndarray,
         np.ndarray,
         np.ndarray,
@@ -356,24 +358,36 @@ def _domain_boundary_sides(
         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.
+        """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 correctly tuned.
+
+        Args:
+            sd: Subdomain grid.
+            tol: Tolerance used to determine whether a face center lies on a boundary side.
+
+        Returns:
+            Tuple of boolean arrays, each one of size sd.num_faces. For a given array,
+                a True element indicates that the face lies on the side of the domain.
+                Otherwise, the element is set to False.
+
         """
-        tol = 1e-10
         box = self.box
-        east = g.face_centers[0] > box["xmax"] - tol
-        west = g.face_centers[0] < box["xmin"] + tol
+        east = sd.face_centers[0] > box["xmax"] - tol
+        west = sd.face_centers[0] < box["xmin"] + 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 = sd.face_centers[1] > box["ymax"] - tol
+            south = sd.face_centers[1] < box["ymin"] + 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()
+            top = sd.face_centers[2] > box["zmax"] - tol
+            bottom = sd.face_centers[2] < box["zmin"] + tol
+        all_bf = sd.get_boundary_faces()
         return all_bf, east, west, north, south, top, bottom

From 4de06f8080c95c9f3b853b9a4d8f70c06048ef43 Mon Sep 17 00:00:00 2001
From: Jhabriel Varela <jhabriel@gmail.com>
Date: Wed, 9 Nov 2022 00:13:50 +0100
Subject: [PATCH 02/14] DOC: Improve sentence.

---
 src/porepy/models/abstract_model.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/porepy/models/abstract_model.py b/src/porepy/models/abstract_model.py
index 61be3e653..84a191f7c 100644
--- a/src/porepy/models/abstract_model.py
+++ b/src/porepy/models/abstract_model.py
@@ -362,7 +362,7 @@ def _domain_boundary_sides(
 
         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 correctly tuned.
+        provided `tol` is tuned accordingly.
 
         Args:
             sd: Subdomain grid.

From 0de2d26da9c4bf0655dcb267c5b0a5da02f3de48 Mon Sep 17 00:00:00 2001
From: Eirik Keilegavlen <Eirik.Keilegavlen@uib.no>
Date: Tue, 8 Nov 2022 14:04:59 +0100
Subject: [PATCH 03/14] BLD: Modified mypy.ini

---
 mypy.ini | 8 +-------
 1 file changed, 1 insertion(+), 7 deletions(-)

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

From 208fddf37fc33430e30bc4ef03a13198522844c1 Mon Sep 17 00:00:00 2001
From: Eirik Keilegavlen <Eirik.Keilegavlen@uib.no>
Date: Wed, 9 Nov 2022 08:58:35 +0100
Subject: [PATCH 04/14] STY: Typing after update to Mypy 0.99

---
 src/porepy/fracs/fracture_importer.py         |  5 +-
 src/porepy/fracs/meshing.py                   |  3 +-
 src/porepy/fracs/msh_2_grid.py                | 44 ++++++++--------
 src/porepy/fracs/structured.py                | 16 +++---
 src/porepy/fracs/wells_3d.py                  | 36 +++++++------
 src/porepy/grids/md_grid.py                   | 10 ++--
 src/porepy/grids/mortar_grid.py               |  8 +--
 src/porepy/grids/partition.py                 |  5 +-
 src/porepy/grids/point_grid.py                |  4 +-
 .../grids/standard_grids/md_grids_2d.py       | 43 +++++++++++----
 src/porepy/grids/standard_grids/utils.py      |  4 +-
 src/porepy/numerics/ad/equation_manager.py    |  4 +-
 src/porepy/numerics/ad/operator_functions.py  |  4 +-
 src/porepy/numerics/fv/fvutils.py             | 40 +++++++-------
 src/porepy/numerics/fv/mpsa.py                | 52 +++++++++----------
 .../interface_laws/elliptic_interface_laws.py |  6 ++-
 src/porepy/numerics/linear_solvers.py         |  8 +--
 src/porepy/numerics/mixed_dim/assembler.py    |  6 ++-
 src/porepy/numerics/vem/dual_elliptic.py      |  2 +-
 src/porepy/params/data.py                     |  2 +-
 20 files changed, 176 insertions(+), 126 deletions(-)

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/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..3472d0c6c 100644
--- a/src/porepy/grids/standard_grids/utils.py
+++ b/src/porepy/grids/standard_grids/utils.py
@@ -1,5 +1,3 @@
-from typing import List
-
 import numpy as np
 
 import porepy as pp
@@ -65,7 +63,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/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/operator_functions.py b/src/porepy/numerics/ad/operator_functions.py
index 7adef0ae6..4442e5981 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, Type, Union
+from typing import Callable, List, Optional, Type, Union
 
 import numpy as np
 import scipy.sparse as sps
@@ -466,7 +466,7 @@ def dummy_rel_perm(s):
 
     def __init__(
         self,
-        func: Callable = None,
+        func: Optional[Callable] = None,
         ad_function_type: Type["AbstractFunction"] = Function,
         operator_kwargs: dict = {},
     ) -> None:
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.

From 8cc4209eb02f9e0a151790a5167249e2ebd2818d Mon Sep 17 00:00:00 2001
From: Eirik Keilegavlen <Eirik.Keilegavlen@uib.no>
Date: Wed, 9 Nov 2022 08:59:06 +0100
Subject: [PATCH 05/14] ENH: Cleanup of tags in FractureNetwork-Gmsh interface.

Tags are now systematically taken as numerical in geometry processing
(this is convenient for splitting fractures etc.) and Enums in the
interface to Gmsh.

Also typing.
---
 src/porepy/fracs/fracture_network_2d.py | 188 ++++++++++++++++--------
 src/porepy/fracs/fracture_network_3d.py | 130 +++++++++-------
 src/porepy/fracs/gmsh_interface.py      |  68 +++++----
 3 files changed, 237 insertions(+), 149 deletions(-)

diff --git a/src/porepy/fracs/fracture_network_2d.py b/src/porepy/fracs/fracture_network_2d.py
index ec53b1933..41dc74131 100644
--- a/src/porepy/fracs/fracture_network_2d.py
+++ b/src/porepy/fracs/fracture_network_2d.py
@@ -5,7 +5,7 @@
 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 +14,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,26 +35,37 @@ 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
+        edges (size: (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
+            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)
+            under geometry processing like intersection removal. Additional tags can
+            be assigned
+        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 (double): Tolerance used in geometric computations.
-        tags (dict): Tags for fractures.
-        _decomposition (dict): Decomposition of the fracture network, used for export to
+        tol: Tolerance used in geometric computations.
+        tags: Tags for fractures.
+        _decomposition: Decomposition of the fracture network, used for export to
             gmsh, and available for later processing. Initially empty, is created by
             self.mesh().
 
     """
 
-    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:
@@ -73,20 +85,23 @@ def __init__(self, pts=None, edges=None, domain=None, tol=1e-8):
             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.domain = domain
+
         self.tol = tol
 
-        self.tags = {}
-        self.bounding_box_imposed = False
-        self._decomposition = {}
+        self.tags: dict = {}
+        self.bounding_box_imposed: bool = False
+        self._decomposition: dict = {}
 
         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 +114,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 +200,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 +247,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 +299,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 +338,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 +421,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 +437,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 +507,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 +529,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 +610,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 +636,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 +698,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 +729,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]))
@@ -716,7 +774,7 @@ 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),
                     )
                 )
@@ -1283,7 +1341,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..76def8d2a 100644
--- a/src/porepy/fracs/fracture_network_3d.py
+++ b/src/porepy/fracs/fracture_network_3d.py
@@ -8,7 +8,7 @@
 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 +16,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 +62,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 +100,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 +122,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 +187,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 +385,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 +407,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 +481,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 +510,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 to 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 +582,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 +1287,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 +1503,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 +1635,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 +1845,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 +1871,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 +1888,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 +1905,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 +1998,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 +2068,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 +2081,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 +2400,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..c86cee664 100644
--- a/src/porepy/fracs/gmsh_interface.py
+++ b/src/porepy/fracs/gmsh_interface.py
@@ -15,7 +15,7 @@
 """
 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 +145,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 +169,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 +208,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 +350,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 +379,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 +391,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 +425,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 +468,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 +489,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 +503,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 +526,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 +550,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)
 

From 15c10eceb8e00f72cdf445edee45a22c3c0990f8 Mon Sep 17 00:00:00 2001
From: Eirik Keilegavlen <Eirik.Keilegavlen@uib.no>
Date: Wed, 9 Nov 2022 11:33:47 +0100
Subject: [PATCH 06/14] BLD: Added py.typed file.

Python and mypy seems to want this.
---
 src/porepy/py.typed | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 create mode 100644 src/porepy/py.typed

diff --git a/src/porepy/py.typed b/src/porepy/py.typed
new file mode 100644
index 000000000..e69de29bb

From 68ae26db760125c55d596fa49af924694416a488 Mon Sep 17 00:00:00 2001
From: Eirik Keilegavlen <Eirik.Keilegavlen@uib.no>
Date: Wed, 9 Nov 2022 11:46:29 +0100
Subject: [PATCH 07/14] STY: Import annotations from future.

---
 src/porepy/fracs/fracture_network_2d.py  | 1 +
 src/porepy/fracs/fracture_network_3d.py  | 1 +
 src/porepy/fracs/gmsh_interface.py       | 1 +
 src/porepy/grids/standard_grids/utils.py | 1 +
 4 files changed, 4 insertions(+)

diff --git a/src/porepy/fracs/fracture_network_2d.py b/src/porepy/fracs/fracture_network_2d.py
index 41dc74131..f1d1ebbc6 100644
--- a/src/porepy/fracs/fracture_network_2d.py
+++ b/src/porepy/fracs/fracture_network_2d.py
@@ -1,6 +1,7 @@
 """
 Module contains class for representing a fracture network in a 2d domain.
 """
+from __future__ import annotations
 import copy
 import csv
 import logging
diff --git a/src/porepy/fracs/fracture_network_3d.py b/src/porepy/fracs/fracture_network_3d.py
index 76def8d2a..ca3dadcc8 100644
--- a/src/porepy/fracs/fracture_network_3d.py
+++ b/src/porepy/fracs/fracture_network_3d.py
@@ -4,6 +4,7 @@
 The model relies heavily on functions in the computational geometry library.
 
 """
+from __future__ import annotations
 import copy
 import csv
 import logging
diff --git a/src/porepy/fracs/gmsh_interface.py b/src/porepy/fracs/gmsh_interface.py
index c86cee664..1a6627e6a 100644
--- a/src/porepy/fracs/gmsh_interface.py
+++ b/src/porepy/fracs/gmsh_interface.py
@@ -13,6 +13,7 @@
         gmsh model. Can also mesh.
 
 """
+from __future__ import annotations
 from dataclasses import dataclass
 from enum import Enum
 from typing import Optional, Union
diff --git a/src/porepy/grids/standard_grids/utils.py b/src/porepy/grids/standard_grids/utils.py
index 3472d0c6c..ddedf72cf 100644
--- a/src/porepy/grids/standard_grids/utils.py
+++ b/src/porepy/grids/standard_grids/utils.py
@@ -1,3 +1,4 @@
+from __future__ import annotations
 import numpy as np
 
 import porepy as pp

From 64431b95c904b41d6bce2fe91eacd81deab96b0e Mon Sep 17 00:00:00 2001
From: Eirik Keilegavlen <Eirik.Keilegavlen@uib.no>
Date: Wed, 9 Nov 2022 11:51:51 +0100
Subject: [PATCH 08/14] STY: isort

---
 src/porepy/fracs/fracture_network_2d.py  | 1 +
 src/porepy/fracs/fracture_network_3d.py  | 1 +
 src/porepy/fracs/gmsh_interface.py       | 1 +
 src/porepy/grids/standard_grids/utils.py | 1 +
 4 files changed, 4 insertions(+)

diff --git a/src/porepy/fracs/fracture_network_2d.py b/src/porepy/fracs/fracture_network_2d.py
index f1d1ebbc6..33976d7d6 100644
--- a/src/porepy/fracs/fracture_network_2d.py
+++ b/src/porepy/fracs/fracture_network_2d.py
@@ -2,6 +2,7 @@
 Module contains class for representing a fracture network in a 2d domain.
 """
 from __future__ import annotations
+
 import copy
 import csv
 import logging
diff --git a/src/porepy/fracs/fracture_network_3d.py b/src/porepy/fracs/fracture_network_3d.py
index ca3dadcc8..a09a4987f 100644
--- a/src/porepy/fracs/fracture_network_3d.py
+++ b/src/porepy/fracs/fracture_network_3d.py
@@ -5,6 +5,7 @@
 
 """
 from __future__ import annotations
+
 import copy
 import csv
 import logging
diff --git a/src/porepy/fracs/gmsh_interface.py b/src/porepy/fracs/gmsh_interface.py
index 1a6627e6a..8ad7a3c9d 100644
--- a/src/porepy/fracs/gmsh_interface.py
+++ b/src/porepy/fracs/gmsh_interface.py
@@ -14,6 +14,7 @@
 
 """
 from __future__ import annotations
+
 from dataclasses import dataclass
 from enum import Enum
 from typing import Optional, Union
diff --git a/src/porepy/grids/standard_grids/utils.py b/src/porepy/grids/standard_grids/utils.py
index ddedf72cf..b98030b59 100644
--- a/src/porepy/grids/standard_grids/utils.py
+++ b/src/porepy/grids/standard_grids/utils.py
@@ -1,4 +1,5 @@
 from __future__ import annotations
+
 import numpy as np
 
 import porepy as pp

From 2a2601caaf651e01322bfd72ef47b597791f74a8 Mon Sep 17 00:00:00 2001
From: Eirik Keilegavlen <Eirik.Keilegavlen@uib.no>
Date: Wed, 9 Nov 2022 14:10:31 +0100
Subject: [PATCH 09/14] Update src/porepy/fracs/fracture_network_3d.py

Co-authored-by: Ivar Stefansson <ivar.stefansson@uib.no>
---
 src/porepy/fracs/fracture_network_3d.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/porepy/fracs/fracture_network_3d.py b/src/porepy/fracs/fracture_network_3d.py
index a09a4987f..2c0e0bef9 100644
--- a/src/porepy/fracs/fracture_network_3d.py
+++ b/src/porepy/fracs/fracture_network_3d.py
@@ -536,7 +536,7 @@ def prepare_for_gmsh(self, mesh_args, dfn=False, constraints=None) -> GmshData3d
         polygon_tags: dict[int, GmshInterfaceTags] = dict()
 
         for fi, _ in enumerate(self._fractures):
-            # Translate to from numerical tags to the GmshInterfaceTags system.
+            # Translate from numerical tags to the GmshInterfaceTags system.
             if has_boundary and self.tags["boundary"][fi]:
                 physical_surfaces[fi] = GmshInterfaceTags.DOMAIN_BOUNDARY_SURFACE
                 polygon_tags[fi] = GmshInterfaceTags.DOMAIN_BOUNDARY_SURFACE

From 88f30590eabf42f40c94553845856f9ba5733ebe Mon Sep 17 00:00:00 2001
From: Eirik Keilegavlen <Eirik.Keilegavlen@uib.no>
Date: Wed, 9 Nov 2022 14:55:10 +0100
Subject: [PATCH 10/14] DOC: Improved documentation of FractureNetwork2d.

---
 src/porepy/fracs/fracture_network_2d.py | 81 +++++++++++++++----------
 1 file changed, 48 insertions(+), 33 deletions(-)

diff --git a/src/porepy/fracs/fracture_network_2d.py b/src/porepy/fracs/fracture_network_2d.py
index 33976d7d6..1467d0a90 100644
--- a/src/porepy/fracs/fracture_network_2d.py
+++ b/src/porepy/fracs/fracture_network_2d.py
@@ -40,24 +40,14 @@ class FractureNetwork2d:
     Parameters:
         pts (siz: 2 x num_pts): Start and endpoints of the fractures. Points
             can be shared by fractures.
-        edges (size: (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. 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)
-            under geometry processing like intersection removal. Additional tags can
-            be assigned
-        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.
+        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.
-        tags: Tags for fractures.
-        _decomposition: Decomposition of the fracture network, used for export to
-            gmsh, and available for later processing. Initially empty, is created by
-            self.mesh().
 
     """
 
@@ -74,32 +64,55 @@ def __init__(
             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
+        self.pts = np.zeros((2, 0)) if pts is None else pts
+        """Start and endpoints of the fractures. Points can be shared by fractures."""
 
-        if edges is None:
-            self.edges = np.zeros((2, 0), dtype=int)
-        else:
-            self.edges = edges
+        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.tags: dict = {}
         self.bounding_box_imposed: bool = False
-        self._decomposition: dict = {}
+        """Flag indicating whether the bounding box has been imposed."""
+
+        ## 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")
@@ -764,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]
@@ -782,7 +795,9 @@ def impose_external_boundary(
                 )
 
             # 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(

From eda68622936bce92367850b1ce3e02d2f8ef5a1c Mon Sep 17 00:00:00 2001
From: Eirik Keilegavlen <Eirik.Keilegavlen@uib.no>
Date: Wed, 9 Nov 2022 15:35:29 +0100
Subject: [PATCH 11/14] BLD: Update CODEOWNERS

---
 .github/CODEOWNERS | 3 +++
 1 file changed, 3 insertions(+)

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

From 40fdd7d988180fa20378567c63d99ea8bd586ab2 Mon Sep 17 00:00:00 2001
From: Jhabriel Varela <jhabriel@gmail.com>
Date: Fri, 11 Nov 2022 11:03:42 +0100
Subject: [PATCH 12/14] MAINT: Refactor of domain boundary sides

---
 src/porepy/models/abstract_model.py | 64 ++++++++++++++++++-----------
 1 file changed, 39 insertions(+), 25 deletions(-)

diff --git a/src/porepy/models/abstract_model.py b/src/porepy/models/abstract_model.py
index 84a191f7c..922e7fb78 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,49 +346,62 @@ def _export(self) -> None:
         """
         pass
 
-    def _domain_boundary_sides(
-        self,
-        sd: pp.Grid,
-        tol: float = 1e-10,
-    ) -> tuple[
-        np.ndarray,
-        np.ndarray,
-        np.ndarray,
-        np.ndarray,
-        np.ndarray,
-        np.ndarray,
-        np.ndarray,
-    ]:
+    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.
 
-        Args:
+        Parameters:
             sd: Subdomain grid.
             tol: Tolerance used to determine whether a face center lies on a boundary side.
 
         Returns:
-            Tuple of boolean arrays, each one of size sd.num_faces. For a given array,
-                a True element indicates that the face lies on the side of the domain.
-                Otherwise, the element is set to False.
+            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)
+                print(sides.north)
 
         """
+
+        # Get domain boundary sides
         box = self.box
-        east = sd.face_centers[0] > box["xmax"] - tol
-        west = sd.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(sd.num_faces, dtype=bool)
             south = north.copy()
         else:
-            north = sd.face_centers[1] > box["ymax"] - tol
-            south = sd.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(sd.num_faces, dtype=bool)
             bottom = top.copy()
         else:
-            top = sd.face_centers[2] > box["zmax"] - tol
-            bottom = sd.face_centers[2] < box["zmin"] + tol
+            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()
-        return all_bf, east, west, north, south, top, bottom
+
+        # 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

From 864d4af3de6848c667088555e8c223437683b7a7 Mon Sep 17 00:00:00 2001
From: Jhabriel Varela <jhabriel@gmail.com>
Date: Fri, 11 Nov 2022 13:47:17 +0100
Subject: [PATCH 13/14] DOC: Incorporate @IvarStefansson suggestion.

---
 src/porepy/models/abstract_model.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/porepy/models/abstract_model.py b/src/porepy/models/abstract_model.py
index 922e7fb78..fe057d969 100644
--- a/src/porepy/models/abstract_model.py
+++ b/src/porepy/models/abstract_model.py
@@ -376,7 +376,10 @@ def _domain_boundary_sides(self, sd: pp.Grid, tol: float = 1e-10) -> NamedTuple:
                 model.prepare_simulation()
                 sd = model.mdg.subdomains()[0]
                 sides = model._domain_boundary_sides(sd)
-                print(sides.north)
+                # 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)
 
         """
 

From 958a55d5b6a6985ba4b554951bd345c04ee0312e Mon Sep 17 00:00:00 2001
From: Veljko Lipovac <68282540+vlipovac@users.noreply.github.com>
Date: Mon, 14 Nov 2022 10:25:03 +0100
Subject: [PATCH 14/14] Maint ad docs (#763)

* DOCS: Update docs of some AD objects.

* MAINT: Typing and docstrings of some AD objects.

* Update src/porepy/numerics/ad/operator_functions.py

Co-authored-by: Eirik Keilegavlen <Eirik.Keilegavlen@uib.no>

* STY: Return type in ADMethod.ad_wrapper fixed (mypy).

* STY: ambigous return type of Operator.evaluate

Co-authored-by: Eirik Keilegavlen <Eirik.Keilegavlen@uib.no>
---
 src/porepy/numerics/ad/functions.py          | 131 +++----
 src/porepy/numerics/ad/operator_functions.py | 255 +++++++-------
 src/porepy/numerics/ad/operators.py          | 341 ++++++++++---------
 3 files changed, 354 insertions(+), 373 deletions(-)

diff --git a/src/porepy/numerics/ad/functions.py b/src/porepy/numerics/ad/functions.py
index 45366cb47..f8c149906 100644
--- a/src/porepy/numerics/ad/functions.py
+++ b/src/porepy/numerics/ad/functions.py
@@ -1,22 +1,26 @@
-"""
-This module contains functions to be wrapped in a pp.ad.Function and used as part
-of compound pp.ad.Operators, i.e. as (terms of) equations.
+"""This module contains functions to be wrapped in a
+:class:`~porepy.numerics.ad.operator_functions.Function` and used as part
+of compound :class:`~porepy.numerics.ad.operators.Operator`, i.e. as (terms of) equations.
+
+Some functions depend on non-ad objects. This requires that the function ``f`` be wrapped
+in an ``ad.Function`` using partial evaluation:
 
-Some functions depend on non-ad objects. This requires that the function (f) be wrapped
-in an ad Function using partial evaluation:
+Examples:
+    >>> from functools import partial
+    >>> AdFunction = pp.ad.Function(partial(f, other_parameter), "name")
+    >>> equation: pp.ad.Operator = AdFunction(var) - 2 * var
 
-    from functools import partial
-    AdFunction = pp.ad.Function(partial(f, other_parameter), "name")
-    equation: pp.ad.Operator = AdFunction(var) - 2 * var
+    with ``var`` being some AD variable.
 
-with var being some ad variable.
+    Note that while the argument to ``AdFunction`` is a
+    :class:`~porepy.numerics.ad.operators.Operator, the wrapping in
+    ``pp.ad.Function`` implies that upon parsing,
+    the argument passed to ``f`` will be an Ad_array.
 
-Note that while the argument to AdFunction is a pp.ad.Operator, the wrapping in
-pp.ad.Function implies that upon parsing, the argument passed to f will be an Ad_array.
 """
 from __future__ import annotations
 
-from typing import Callable, Union
+from typing import Callable
 
 import numpy as np
 import scipy.sparse as sps
@@ -90,25 +94,18 @@ def l2_norm(dim: int, var: pp.ad.Ad_array) -> pp.ad.Ad_array:
     """L2 norm of a vector variable.
 
     For the example of dim=3 components and n vectors, the ordering is assumed
-    to be
-        [u0, v0, w0, u1, v1, w1, ..., un, vn, wn]
+    to be ``[u0, v0, w0, u1, v1, w1, ..., un, vn, wn]``
 
     Vectors satisfying ui=vi=wi=0 are assigned zero entries in the jacobi matrix
 
-    Usage note:
+    Note:
         See module level documentation on how to wrap functions like this in ad.Function.
 
-    Parameters
-    ----------
-    dim : int
-        Dimension, i.e. number of vector components.
-    var : pp.ad.Ad_array
-        Ad operator (variable or expression) which is argument of the norm
-        function.
-
-    Returns
-    -------
-    pp.ad.Ad_array
+    Parameters:
+        dim: Dimension, i.e. number of vector components.
+        var: Ad operator (variable or expression) which is argument of the norm function.
+
+    Returns:
         The norm of var with appropriate val and jac attributes.
 
     """
@@ -265,30 +262,24 @@ def heaviside(var, zerovalue: float = 0.5):
 
 
 def heaviside_smooth(var, eps: float = 1e-3):
-    """
-    Smooth (regularized) version of the Heaviside function.
-
-    Parameters
-    ----------
-    var : Ad_array or ndarray
-        Input array.
-    eps : float, optional
-        Regularization parameter. The default is 1E-3. The function will
-        convergence to the Heaviside function in the limit when eps --> 0
-
-    Returns
-    -------
-    Ad_array or ndarray (depending on the input)
-        Regularized heaviside function (and its Jacobian if applicable).
-
-    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
+    """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.
+        eps (optional): Regularization parameter. The function will converge to the
+            Heaviside function in the limit when ``eps --> 0``. The default is ``1e-3``.
+
+    Returns:
+        Regularized heaviside function (and its Jacobian if applicable) in form of a
+        Ad_array or ndarray (depending on the input).
 
     """
     if isinstance(var, Ad_array):
@@ -315,26 +306,19 @@ def __call__(self, var, zerovalue: float = 0.5):
             return np.heaviside(var)  # type: ignore
 
 
-def maximum(
-    var0: pp.ad.Ad_array, var1: Union[pp.ad.Ad_array, np.ndarray]
-) -> pp.ad.Ad_array:
+def maximum(var0: pp.ad.Ad_array, var1: pp.ad.Ad_array | np.ndarray) -> pp.ad.Ad_array:
     """Ad maximum function represented as an Ad_array.
 
     The second argument is allowed to be constant, with a numpy array originally
     wrapped in a pp.ad.Array, whereas the first argument is expected to be an
     Ad_array originating from a pp.ad.Operator.
 
-    Parameters
-    ----------
-    var0 : pp.ad.Ad_array
-        Ad operator (variable or expression).
-    var1 : Union[pp.ad.Ad_array, pp.ad.Array]
-        Ad operator (variable or expression) OR ad Array.
+    Parameters:
+        var0: First argument to the maximum function.
+        var1: Second argument.
 
-    Returns
-    -------
-    pp.ad.Ad_array
-        The maximum of var0 and var1 with appropriate val and jac attributes.
+    Returns:
+        The maximum of ``var0`` and ``var1`` with appropriate val and jac attributes.
 
     """
     vals = [var0.val.copy()]
@@ -357,22 +341,17 @@ def maximum(
 def characteristic_function(tol: float, var: pp.ad.Ad_array):
     """Characteristic function of an ad variable.
 
-    Returns 1 if var.val is within absolute tolerance = tol of zero. The derivative is set to
-    zero independent of var.val.
+    Returns 1 if ``var.val`` is within absolute tolerance = ``tol`` of zero.
+    The derivative is set to zero independent of ``var.val``.
 
-    Usage note:
-        See module level documentation on how to wrap functions like this in ad.Function.
+    Note:
+        See module level documentation on how to wrap functions like this in ``ad.Function``.
 
-    Parameters
-    ----------
-    tol : float
-        Absolute tolerance for comparison with 0 using np.isclose.
-    var : pp.ad.Ad_array
-        Ad operator (variable or expression).
+    Parameters:
+        tol: Absolute tolerance for comparison with 0 using np.isclose.
+        var: Ad operator (variable or expression).
 
-    Returns
-    -------
-    pp.ad.Ad_array
+    Returns:
         The characteristic function of var with appropriate val and jac attributes.
 
     """
diff --git a/src/porepy/numerics/ad/operator_functions.py b/src/porepy/numerics/ad/operator_functions.py
index 4442e5981..3c38ca5b2 100644
--- a/src/porepy/numerics/ad/operator_functions.py
+++ b/src/porepy/numerics/ad/operator_functions.py
@@ -1,14 +1,14 @@
-"""
-Contains callable operators representing functions to be called with other
+"""This module contains callable operators representing functions to be called with other
 operators as input arguments.
 Contains also a decorator class for callables, which transforms them automatically in the
 specified operator function type.
+
 """
 from __future__ import annotations
 
 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
@@ -38,7 +38,7 @@ class AbstractFunction(Operator):
     The abstraction intends to provide means for approximating operators, where values are
     e.g. interpolated and the Jacobian is approximated using FD.
 
-    IMPLEMENTATION NOTE:
+    Note:
         One can flag the operator as ``ad_compatible``. If flagged, the AD framework passes
         AD arrays directly to the callable ``func`` and will **not** call the abstract methods
         for values and the Jacobian during operator parsing.
@@ -47,6 +47,25 @@ class AbstractFunction(Operator):
 
         For now only one child class, porepy.ad.Function, flags itself always as AD compatible.
 
+    Parameters:
+        func: callable Python object representing a (numeric) function.
+            Expected to take numerical information in some form and return numerical
+            information in the same form.
+        name: name of this instance as an AD operator
+        array_compatible (optional): If true, the callable ``func`` will be called
+            using arrays (numpy.typing.ArrayLike). Flagging this true, the user ensures
+            that the callable can work with arrays and return respectively
+            formatted output. If false, the function will be evaluated element-wise
+            (scalar input). Defaults to False.
+        ad_compatible (Optional): If true, the callable ``func`` will be called using
+            the porepy.ad.Ad_array.
+
+            Note that as of now, this will effectively bypass the abstract methods
+            for generating values and the Jacobian, assuming both will be provided
+            correctly by the return value of ``func``.
+
+            Defaults to False.
+
     """
 
     def __init__(
@@ -56,26 +75,6 @@ def __init__(
         array_compatible: bool = False,
         ad_compatible: bool = False,
     ):
-        """Configures this instance as a valid AD operator.
-        The passed callable is expected to take numerical information in some form and
-        return numerical information in the same form form.
-
-        Args:
-            func (Callable): callable Python object representing a (numeric) function
-            name (str): name of this instance as an AD operator
-            array_compatible (default=False): If true, the callable ``func`` will be called
-                using arrays (numpy.typing.ArrayLike). Flagging this true, the user ensures
-                that the callable can work with arrays and return respectively
-                formatted output. If false, the function will be evaluated element-wise
-                (scalar input)
-            ad_compatible (default=False): If true, the callable ``func`` will be called using
-                the porepy.ad.Ad_array.
-
-                Important NOTE: As of now, this will effectively bypass the abstract methods
-                for generating values and the Jacobian, assuming both will be provided
-                correctly by the return value of ``func``
-
-        """
         ### PUBLIC
         # Reference to callable passed at instantiation
         self.func: Callable = func
@@ -89,16 +88,16 @@ def __init__(
         self._operation: Operation = Operation.evaluate
         self._set_tree()
 
-    def __call__(self, *args: Union[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'.
 
-        Args:
+        Parameters:
             *args: AD operators passed as symbolic arguments for the callable passed at
                 instantiation.
 
         Returns:
-            porepy.ad.Operator: Operator with call-arguments as children in the operator tree.
-                The assigned operation is ``evaluate``.
+            Operator with call-arguments as children in the operator tree.
+            The assigned operation is ``evaluate``.
 
         """
         children = [self, *args]
@@ -142,18 +141,18 @@ def parse(self, md: pp.MixedDimensionalGrid):
         The real work will be done by combining the function with arguments, during
         parsing of an operator tree.
 
-        Args:
+        Parameters:
             md: Mixed-dimensional grid.
 
         Returns:
-            The instance itself
+            The instance itself.
+
         """
         return self
 
     @abc.abstractmethod
     def get_values(self, *args: Ad_array) -> np.ndarray:
-        """
-        Abstract method for evaluating the callable passed at instantiation.
+        """Abstract method for evaluating the callable passed at instantiation.
 
         This method will be called during the operator parsing.
         The AD arrays passed as arguments will be in the same order as the operators passed to
@@ -162,12 +161,13 @@ def get_values(self, *args: Ad_array) -> np.ndarray:
         The returned numpy array will be set as 'val' argument for the AD array representing
         this instance.
 
-        Args:
+        Parameters:
             *args: Ad_array representation of the operators passed during the call to this
                 instance
 
         Returns:
-            numpy.ndarray: numeric function values
+            Function values in numerical format.
+
         """
         pass
 
@@ -183,28 +183,35 @@ def get_jacobian(self, *args: Ad_array) -> sps.spmatrix:
         The returned numpy array will be be set as 'jac' argument for the AD array representing
         this instance.
 
-        NOTE: the necessary dimensions for the jacobian can be extracted from the dimensions
-        of the Jacobians of passed Ad_array instances.
+        Note:
+            The necessary dimensions for the jacobian can be extracted from the dimensions
+            of the Jacobians of passed Ad_array instances.
 
-        Args:
+        Parameters:
             *args: Ad_array representation of the operators passed during the call to this
                 instance
 
         Returns:
-            numpy.ndarray: numeric representation of the Jacobian of this function
+            Numeric representation of the Jacobian of this function.
+
         """
         pass
 
 
 class AbstractJacobianFunction(AbstractFunction):
-    """'Half'-abstract base class, providing a call to the callable ``func`` in order to obtain
-    numeric function values.
+    """Partially abstract base class, providing a call to the callable ``func`` in order to
+    obtain numeric function values.
 
     What remains abstract is the Jacobian.
 
     """
 
     def get_values(self, *args: Ad_array) -> np.ndarray:
+        """
+        Returns:
+            The direct evaluation of the callable using ``val`` of passed AD arrays.
+
+        """
         # get values of argument Ad_arrays.
         vals = (arg.val for arg in args)
 
@@ -231,14 +238,15 @@ class Function(AbstractFunction):
     where it is expected that passing Ad_arrays directly to ``func`` will
     return the proper result.
 
-    Here the values AND the Jacobian are obtained exactly by the AD framework.
+    Here the values **and** the Jacobian are obtained exactly by the AD framework.
 
     The intended use is as a wrapper for operations on pp.ad.Ad_array objects,
     in forms which are not directly or easily expressed by the rest of the Ad
     framework.
 
-    NOTE: This is a special case where the abstract methods for getting values and the Jacobian
-    are formally implemented but never used by the AD framework (see ``ad_compatible``)
+    Note:
+        This is a special case where the abstract methods for getting values and the Jacobian
+        are formally implemented but never used by the AD framework (see ``ad_compatible``)
 
     """
 
@@ -259,14 +267,13 @@ class ConstantFunction(AbstractFunction):
     zero Jacobian.
 
     It still has to be called though since it fulfills the notion of a 'function'.
+
+    Parameters:
+        values: constant values per cell.
+
     """
 
     def __init__(self, name: str, values: np.ndarray):
-        """
-        Args:
-            values: constant values per cell
-
-        """
         # dummy function, takes whatever and returns only the pre-set values
         def func(*args):
             return values
@@ -284,13 +291,13 @@ def get_values(self, *args: Ad_array) -> np.ndarray:
 
     def get_jacobian(self, *args: Ad_array) -> sps.spmatrix:
         """
-        IMPLEMENTATION NOTE:
-        The return value is not a sparse matrix as imposed by the parent method signature,
-        but a zero.
-        Numerical operations with a zero always works with any numeric formats in
-        numpy, scipy and PorePy's AD framework.
-        Since the constant function (most likely) gets no arguments passed,
-        we have no way of knowing the necessary SHAPE for a zero matrix. Hence scalar.
+        Note:
+            The return value is not a sparse matrix as imposed by the parent method signature,
+            but a zero.
+            Numerical operations with a zero always works with any numeric formats in
+            numpy, scipy and PorePy's AD framework.
+            Since the constant function (most likely) gets no arguments passed,
+            we have no way of knowing the necessary SHAPE for a zero matrix. Hence scalar.
 
         Returns:
             The trivial derivative of a constant.
@@ -300,25 +307,23 @@ def get_jacobian(self, *args: Ad_array) -> sps.spmatrix:
 
 
 class DiagonalJacobianFunction(AbstractJacobianFunction):
-    """
-    Approximates the Jacobian of the function using identities and scalar
-    multipliers per dependency.
+    """Approximates the Jacobian of the function using identities and scalar multipliers
+    per dependency.
+
+    Parameters:
+        multipliers: scalar multipliers for the identity blocks in the Jacobian,
+            per dependency of ``func``. The order in ``multipliers`` is expected to match
+            the order of AD operators passed to the call of this function.
+
     """
 
     def __init__(
         self,
         func: Callable,
         name: str,
-        multipliers: Union[List[float], float],
+        multipliers: float | list[float],
         array_compatible: bool = False,
     ):
-        """
-        Args:
-            multipliers: scalar multipliers for the identity blocks in the Jacobian,
-                per dependency of ``func``. The order in ``multipliers`` is expected to match
-                the order of AD operators passed to the call of this function.
-
-        """
         super().__init__(func, name, array_compatible)
         # check and format input for further use
         if isinstance(multipliers, list):
@@ -351,6 +356,22 @@ class InterpolatedFunction(AbstractFunction):
     The image of the function is expected to be of dimension 1, while the domain can be
     multidimensional.
 
+    Note:
+        All vector-valued ndarray arguments are assumed to be column vectors.
+        Each row-entry represents a value for an argument of ``func`` in
+        respective order.
+
+    Parameters:
+        min_val: lower bounds for the domain of ``func``.
+        max_val: upper bound for the domain.
+        npt: number of interpolation points per dimension of the domain.
+        order: Order of interpolation. Supports currently only linear order.
+        preval (optional): If True, pre-evaluates the values of the function at
+            the points of interpolation and stores them.
+            If False, evaluates them if necessary.
+            Influences the runtime.
+            Defaults to False.
+
     """
 
     def __init__(
@@ -364,20 +385,6 @@ def __init__(
         preval: bool = False,
         array_compatible: bool = False,
     ):
-        """All vector-valued ndarray arguments are assumed to be column vectors.
-        Each row-entry represents a value for an argument of ``func`` in
-        respective order.
-
-        Args:
-            min_val: lower bounds for the domain of ``func``.
-            max_val: upper bound for the domain.
-            npt: number of interpolation points per dimension of the domain.
-            order: Order of interpolation. Supports currently only linear order.
-            preval (default=False): If True, pre-evaluates the values of the function at
-                the points of interpolation and stores them.
-                If False, evaluates them if necessary.
-                Influences the runtime.
-        """
         super().__init__(func, name, array_compatible)
 
         ### PUBLIC
@@ -432,59 +439,53 @@ def get_jacobian(self, *args: Ad_array) -> sps.spmatrix:
 
 
 class ADmethod:
-    """
-    Automatic-Differentiation method/function.
-
-    (Decorator) Class for methods representing e.g., physical properties.
+    """(Decorator) Class for methods representing e.g., physical properties.
     The decorated function is expected to take scalars/vectors and return a scalar/vector.
-    See example usage below.
 
     The return value will be an AD operator of a type passed to the decorator.
 
-    EXAMPLE USAGE:
-    .. code-block:: python
-        import porepy as pp
+    Examples:
+        .. code:: python
+
+            import porepy as pp
+
+            # decorating class methods
+            class IdealGas:
 
-        # decorating class methods
-        class IdealGas:
+                @ADmethod(ad_operator=pp.ad.DiagonalJacobianFunction,
+                        operators_args={"multipliers"=[1,1]})
+                def density(self, p: float, T: float) -> float:
+                    return p/T
 
-            @ADmethod(ad_operator=pp.ad.DiagonalJacobianFunction,
-                      operators_args={"multipliers"=[1,1]})
-            def density(self, p: float, T: float) -> float:
-                return p/T
+            # decorating function
+            @ADmethod(ad_operator=pp.ad.Function)
+            def dummy_rel_perm(s):
+                return s**2
 
-        # decorating function
-        @ADmethod(ad_operator=pp.ad.Function)
-        def dummy_rel_perm(s):
-            return s**2
+        With above code, the density of an instance of ``IdealGas`` can be called using
+        :class:`~porepy.numerics.ad.operators.MergedVariable` representing
+        pressure and temperature.
+        Analogously, ``dummy_rel_perm`` can be called with one representing the saturation.
+
+    Note:
+        If used as decorator WITHOUT explicit instantiation, the instantiation will be
+        done implicitly with above default arguments (that's how Python decorators work).
+
+    Parameters:
+        func: decorated function object
+        ad_function_type: type reference to an AD operator function to be instantiated.
+            When instantiated, that instance will effectively replace ``func``.
+        operator_kwargs: keyword arguments to be passed when instantiating an operator
+            of type ``ad_function_type``.
 
-    With above code, the density of an instance of 'IdealGas' can be called using
-    :class:`~porepy.numerics.ad.operators.MergedVariable` representing
-    pressure and temperature.
-    Analogously, `dummy_rel_perm` can be called with one representing the saturation.
     """
 
     def __init__(
         self,
         func: Optional[Callable] = None,
-        ad_function_type: Type["AbstractFunction"] = Function,
+        ad_function_type: Type[AbstractFunction] = Function,
         operator_kwargs: dict = {},
     ) -> None:
-        """
-        Decorator class constructor.
-        Saves information about the requested AD Function type and keyword arguments necessary
-        for its instantiation.
-
-        NOTE: If used as decorator WITHOUT explicit instantiation, the instantiation will be
-        done implicitly with above default arguments (that's how Python decorators work).
-
-        Args:
-            func: decorated function object
-            ad_function_type: type reference to an AD operator function to be instantiated.
-                When instantiated, that instance will effectively replace ``func``.
-            operator_kwargs: keyword arguments to be passed when instantiating an operator
-                of type ``ad_function_type``.
-        """
         # reference to decorated function object
         self._func = func
         # mark if decoration without explicit call to constructor
@@ -498,8 +499,7 @@ def __init__(
         self._op_kwargs = operator_kwargs
 
     def __call__(self, *args, **kwargs) -> Union[ADmethod, pp.ad.Operator]:
-        """
-        Wrapper factory.
+        """Wrapper factory.
         The decorated object is wrapped and/or evaluated here.
 
         Dependent on whether the decorated function is a method belonging to a class,
@@ -508,7 +508,7 @@ def __call__(self, *args, **kwargs) -> Union[ADmethod, pp.ad.Operator]:
         If bound to a class instance, the wrapper will include a partial function, where
         the instance of the class was already passed beforehand.
 
-        IMPLEMENTATION NOTE:
+        Note:
             If the decorator was explicitly instantiated during decoration,
             that instance will effectively be replaced by another decorator instance created
             here in the call.
@@ -516,6 +516,7 @@ def __call__(self, *args, **kwargs) -> Union[ADmethod, pp.ad.Operator]:
             used as a decorator, hence properly dereferencing the original instance.
             If used differently or if another reference is saved between explicit instantiation
             and call, this is a potential memory leak.
+
         """
         # If decorated without explicit init, the function is passed during a call to
         # the decorator as first argument.
@@ -572,11 +573,10 @@ def __get__(self, binding_instance: object, binding_type: type) -> Callable:
         are always passed in unbound form to the decorator when the code is
         evaluated
         (i.e. they don't have a reference to the `self` argument, contrary to bound
-         methods)
+        methods)
 
-         Args:
-            binding_instance: instance, whose method/attribute has been decorated
-                by this class.
+        Parameters:
+            binding_instance: instance, whose method has been decorated by this class.
             binding_type: type instance of the decorated method's class/owner.
 
         """
@@ -586,14 +586,15 @@ 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.
 
-        Args:
+        Parameters:
             *args: arguments for the call to the wrapping AD operator function
             **kwargs: keyword argument for the call to the wrapping Ad operator function
+
         """
         # Make sure proper assignment of callable was made
         assert self._func is not None
diff --git a/src/porepy/numerics/ad/operators.py b/src/porepy/numerics/ad/operators.py
index 8b3b04432..3f96b5fbd 100644
--- a/src/porepy/numerics/ad/operators.py
+++ b/src/porepy/numerics/ad/operators.py
@@ -43,18 +43,25 @@ def _get_shape(mat):
 
 
 class Operator:
-    """Superclass for all Ad operators.
+    """Parent class for all AD operators.
 
-    Objects of this class is not meant to be initiated directly; rather the various
+    This class is not meant to be initiated directly, rather the various
     subclasses should be used. Instances of this class will still be created when
     subclasses are combined by operations.
 
-    Attributes:
-        subdomains: List of subdomains (subdomains) on which the operator is defined. Will be
-            empty for operators not associated with specific subdomains.
-        interfaces: List of edges (tuple of subdomains) in the mixed-dimensional grid on which
-            the operator is defined. Will be empty for operators not associated with
-            specific subdomains.
+    Contains a tree structure of child operators for the recursive forward evaluation.
+
+    Provides overload functions for basic arithmetic operations.
+
+    Parameters:
+        name: Name of this operator. Used for string representations
+        subdomains (optional): List of subdomains on which the operator is defined.
+            Will be empty for operators not associated with any subdomains.
+            Defaults to None (converted to empty list).
+        interfaces (optional): List of interfaces in the mixed-dimensional grid on which the
+            operator is defined.
+            Will be empty for operators not associated with any interface.
+            Defaults to None (converted to empty list).
 
     """
 
@@ -68,9 +75,23 @@ def __init__(
         if name is None:
             name = ""
         self._name = name
+
+        self._set_tree(tree)
+
         self.subdomains: list[pp.Grid] = [] if subdomains is None else subdomains
+        """List of subdomains on which the operator is defined, passed at instantiation.
+
+        Will be empty for operators not associated with specific subdomains.
+
+        """
+
         self.interfaces: list[pp.MortarGrid] = [] if interfaces is None else interfaces
-        self._set_tree(tree)
+        """List of mortar grids in the mixed-dimensional grid on which the operator is defined,
+        passed at instantiation.
+
+        Will be empty for operators not associated with specific subdomains.
+
+        """
 
     def _set_tree(self, tree=None):
         if tree is None:
@@ -230,11 +251,8 @@ def _identify_discretizations(
         return _ad_utils.uniquify_discretization_list(all_discr)
 
     def discretize(self, mdg: pp.MixedDimensionalGrid) -> None:
-        """Perform discretization operation on all discretizations identified in
-        the tree of this operator, using data from mdg.
-
-        IMPLEMENTATION NOTE: The discretizations was identified at initialization of
-        Expression - it is now done here to accommodate updates (?) and
+        """Perform the ``discretize`` function of all child operators which are discretizations
+        using data from mdg.
 
         """
         unique_discretizations: dict[
@@ -243,15 +261,23 @@ def discretize(self, mdg: pp.MixedDimensionalGrid) -> None:
         _ad_utils.discretize_from_list(unique_discretizations, mdg)
 
     def is_leaf(self) -> bool:
-        """Check if this operator is a leaf in the tree-representation of an object.
+        """Check if this operator is a leaf in the tree-representation of an expression.
+
+        Note that this implies that the method ``parse()`` is expected to be implemented.
 
         Returns:
-            bool: True if the operator has no children. Note that this implies that the
-                method parse() is expected to be implemented.
+            True if the operator has no children.
+
         """
         return len(self.tree.children) == 0
 
     def set_name(self, name: str) -> None:
+        """Reset this object's name originally passed at instantiation.
+
+        Parameters:
+            name: the new name to be assigned.
+
+        """
         self._name = name
 
     def previous_timestep(self) -> pp.ad.Operator:
@@ -264,7 +290,7 @@ def previous_timestep(self) -> pp.ad.Operator:
 
         Returns:
             A copy of self, with all time dependent operators evaluated at the previous
-                time step.
+            time step.
 
         """
         # Create a copy of the operator tree evaluated at a previous time step. This is done
@@ -335,11 +361,19 @@ def _traverse_tree(op: Operator) -> Operator:
 
     def parse(self, mdg: pp.MixedDimensionalGrid) -> Any:
         """Translate the operator into a numerical expression.
+
         Subclasses that represent atomic operators (leaves in a tree-representation of
         an operator) should override this method to return e.g. a number, an array or a
         matrix.
         This method should not be called on operators that are formed as combinations
-        of atomic operators; such operators should be evaluated by the method evaluate().
+        of atomic operators; such operators should be evaluated by the method :meth:`evaluate`.
+
+        Parameters:
+            mdg: Mixed-dimensional grid on which this operator is to be parsed.
+
+        Returns:
+            A numerical format representing this operator;s values on given domain.
+
         """
         raise NotImplementedError("This type of operator cannot be parsed right away")
 
@@ -650,7 +684,7 @@ def __str__(self) -> str:
         return self._name if self._name is not None else ""
 
     def viz(self):
-        """Give a visualization of the operator tree that has this operator at the top."""
+        """Draws a visualization of the operator tree that has this operator as its root."""
         G = nx.Graph()
 
         def parse_subgraph(node):
@@ -709,16 +743,17 @@ def evaluate(
         """Evaluate the residual and Jacobian matrix for a given state.
 
         Parameters:
-            dof_manager (pp.DofManager): used to represent the problem. Will be used
+            dof_manager: Used to represent the problem. Will be used
                 to parse the sub-operators that combine to form this operator.
-            state (np.ndarray, optional): State vector for which the residual and its
+            state (optional): State vector for which the residual and its
                 derivatives should be formed. If not provided, the state will be pulled from
                 the previous iterate (if this exists), or alternatively from the state
                 at the previous time step.
 
         Returns:
-            An Ad-array representation of the residual and Jacobian. Note that the Jacobian
-                matrix need not be invertible, or ever square; this depends on the operator.
+            A representation of the residual and Jacobian in form of an AD Array.
+            Note that the Jacobian matrix need not be invertible, or even square;
+            this depends on the operator.
 
         """
         # Get the mixed-dimensional grid used for the dof-manager.
@@ -870,27 +905,24 @@ def _parse_other(self, other):
 class Matrix(Operator):
     """Ad representation of a sparse matrix.
 
-    For dense matrices, use an Array instead.
+    For dense matrices, use :class:`Array` instead.
 
     This is a shallow wrapper around the real matrix; it is needed to combine the matrix
     with other types of Ad objects.
 
-    Attributes:
-        shape (Tuple of ints): Shape of the wrapped matrix.
+    Parameters:
+        mat: Sparse matrix to be wrapped as an AD operator.
+        name: Name of this operator
 
     """
 
     def __init__(self, mat: sps.spmatrix, name: Optional[str] = None) -> None:
-        """Construct an Ad representation of a matrix.
-
-        Parameters:
-            mat (sps.spmatrix): Sparse matrix to be represented.
-
-        """
         super().__init__(name=name)
         self._mat = mat
         self._set_tree()
+
         self.shape = mat.shape
+        """Shape of the wrapped matrix."""
 
     def __repr__(self) -> str:
         return f"Matrix with shape {self._mat.shape} and {self._mat.data.size} elements"
@@ -901,32 +933,31 @@ def __str__(self) -> str:
             s += self._name
         return s
 
-    def parse(self, mdg) -> sps.spmatrix:
-        """Convert the Ad matrix into an actual matrix.
-
-        Parameters:
-            mdg (pp.MixedDimensionalGrid): Mixed-dimensional grid. Not used, but it is
-                needed as input to be compatible with parse methods for other operators.
+    def parse(self, mdg: pp.MixedDimensionalGrid) -> sps.spmatrix:
+        """See :meth:`Operator.parse`.
 
         Returns:
-            sps.spmatrix: The wrapped matrix.
+            The wrapped matrix.
 
         """
         return self._mat
 
     def transpose(self) -> "Matrix":
+        """Returns an AD operator representing the transposed matrix."""
         return Matrix(self._mat.transpose())
 
 
 class Array(Operator):
-    """Ad representation of a numpy array.
+    """AD representation of a constant numpy array.
 
-    For sparse matrices, use a Matrix instead.
+    For sparse matrices, use :class:`Matrix` instead.
+    For time-dependent arrays see :class:`TimeDependentArray`.
 
     This is a shallow wrapper around the real array; it is needed to combine the array
-    with other types of Ad objects.
+    with other types of AD operators.
 
-    See also TimeDependentArray.
+    Parameters:
+        values: Numpy array to be represented.
 
     """
 
@@ -934,7 +965,7 @@ def __init__(self, values: np.ndarray, name: Optional[str] = None) -> None:
         """Construct an Ad representation of a numpy array.
 
         Parameters:
-            values (np.ndarray): Numpy array to be represented.
+            values: Numpy array to be represented.
 
         """
         super().__init__(name=name)
@@ -951,14 +982,10 @@ def __str__(self) -> str:
         return s
 
     def parse(self, mdg: pp.MixedDimensionalGrid) -> np.ndarray:
-        """Convert the Ad Array into an actual array.
-
-        Parameters:
-            mdg (pp.MixedDimensionalGrid): Mixed-dimensional grid. Not used, but it is
-                needed as input to be compatible with parse methods for other operators.
+        """See :meth:`Operator.parse`.
 
         Returns:
-            np.ndarray: The wrapped array.
+            The wrapped array.
 
         """
         return self._values
@@ -968,20 +995,27 @@ class TimeDependentArray(Array):
     """An Ad-wrapper around a time-dependent numpy array.
 
     The array is tied to a MixedDimensionalGrid, and is distributed among the data
-    dictionaries associated with subdomains and interfaces. The array values are stored
-    in data[pp.STATE][pp.ITERATE][self._name] for the current time and
-    data[pp.STATE][self._name] for the previous time.
+    dictionaries associated with subdomains and interfaces.
+    The array values are stored
+    in ``data[pp.STATE][pp.ITERATE][self._name]`` for the current time and
+    ``data[pp.STATE][self._name]`` for the previous time.
 
-    The array can be differentiated in time using pp.ad.dt().
+    The array can be differentiated in time using ``pp.ad.dt()``.
 
     The intended use is to represent time-varying quantities in equations, e.g., source
     terms. Future use will also include numerical values of boundary conditions,
     however, this is pending an update to the model classes.
 
-    Attributes:
-        prev_time (boolean): If True, the array will be evaluated using
-            data[pp.STATE] (data being the data dictionaries for subdomains and
-            interfaces), if False, data[pp.STATE][pp.ITERATE] is used.
+    Parameters:
+        name: Name of the variable. Should correspond to items in data[pp.STATE].
+        subdomains: Subdomains on which the array is defined. Defaults to None.
+        interfaces: Interfaces on which the array is defined. Defaults to None.
+            Exactly one of subdomains and interfaces must be non-empty.
+        previous_timestep: Flag indicating if the array should be evaluated at the
+            previous time step.
+
+    Raises:
+        ValueError: If either none of, or both of, subdomains and interfaces are empty.
 
     """
 
@@ -992,23 +1026,6 @@ def __init__(
         interfaces: Optional[list[pp.MortarGrid]] = None,
         previous_timestep: bool = False,
     ):
-        """Initialize a TimeDependentArray.
-
-        Args:
-            name: Name of the variable. Should correspond to items in data[pp.STATE].
-            subdomains: Subdomains on which the array is defined. Defaults to None.
-            interfaces: Interfaces on which the array is defined. Defaults to None.
-
-            Exactly one of subdomains and interfaces must be non-empty.
-
-            previous_timestep: Flag indicating if the array should be evaluated at the
-                previous time step.
-
-        Raises:
-            ValueError: If either none of, or both of, subdomains and interfaces are
-                empty.
-
-        """
 
         self._name: str = name
 
@@ -1047,12 +1064,17 @@ def __init__(
             self._is_interface_array = True
 
         self.prev_time: bool = previous_timestep
+        """If True, the array will be evaluated using ``data[pp.STATE]``
+        (data being the data dictionaries for subdomains and interfaces).
+
+        If False, ``data[pp.STATE][pp.ITERATE]`` is used.
+
+        """
 
         self._set_tree()
 
     def previous_timestep(self) -> TimeDependentArray:
-        """Return a representation of this variable on the previous time step.
-
+        """
         Returns:
             This array represented at the previous time step.
 
@@ -1070,10 +1092,11 @@ def parse(self, mdg: pp.MixedDimensionalGrid) -> np.ndarray:
         """Convert this array into numerical values.
 
         The numerical values will be picked from the representation of the array in
-        data[pp.STATE][pp.ITERATE] (where data is the data dictionary of the subdomains
-        or interfaces of this Array), or, if self.prev_time = True, from data[pp.STATE].
+        ``data[pp.STATE][pp.ITERATE]`` (where data is the data dictionary of the subdomains
+        or interfaces of this Array), or, if ``self.prev_time = True``,
+        from ``data[pp.STATE]``.
 
-        Args:
+        Parameters:
             mdg: Mixed-dimensional grid.
 
         Returns:
@@ -1113,7 +1136,7 @@ def __repr__(self) -> str:
 class Scalar(Operator):
     """Ad representation of a scalar.
 
-    This is a shallow wrapper around the real scalar; it may be useful to combine
+    This is a shallow wrapper around a real scalar. It may be useful to combine
     the scalar with other types of Ad objects.
 
     NOTE: Since this is a wrapper around a Python immutable, copying a Scalar will
@@ -1124,12 +1147,6 @@ class Scalar(Operator):
     """
 
     def __init__(self, value: float, name: Optional[str] = None) -> None:
-        """Construct an Ad representation of a float.
-
-        Parameters:
-            value (float): Number to be represented
-
-        """
         super().__init__(name=name)
         self._value = value
         self._set_tree()
@@ -1144,40 +1161,35 @@ def __str__(self) -> str:
         return s
 
     def parse(self, mdg: pp.MixedDimensionalGrid) -> float:
-        """Convert the Ad Scalar into an actual number.
-
-        Parameters:
-            mdg (pp.MixedDimensionalGrid): Mixed-dimensional grid. Not used, but it is
-                needed as input to be compatible with parse methods for other operators.
+        """See :meth:`Operator.parse`.
 
         Returns:
-            float: The wrapped number.
+            The wrapped number.
 
         """
         return self._value
 
 
 class Variable(Operator):
-    """Ad representation of a variable defined on a single Grid or MortarGrid.
+    """AD operator representing a variable defined on a single grid or mortar grid.
 
-    For combinations of variables on different subdomains, see MergedVariable.
+    For combinations of variables on different subdomains, see :class:`MergedVariable`.
 
-    Conversion of the variable into numerical value should be done with respect to the
-    state of an array; see the method evaluate(). Therefore, the variable does not
-    implement a parse() method.
+    A conversion of the variable into numerical value should be done with respect to the
+    state of an array; see :meth:`Operator.evaluate`. Therefore, the variable does not
+    implement the method :meth:`Operator.parse`.
 
-    Attributes:
-        id (int): Unique identifier of this variable.
-        prev_iter (boolean): Whether the variable represents the state at the
-            previous iteration.
-        prev_time (boolean): Whether the variable represents the state at the
-            previous time step.
-        subdomains: List with one item, giving the single grid on which the operator is
-            defined.
-        interfaces: List with one item, giving the single edge (tuple of subdomains) on
-            which the operator is defined.
+    A variable is associated with either a grid or an interface. Therefore it is assumed that
+    either ``subdomains`` or ``interfaces`` is passed as an argument.
 
-        It is assumed that exactly one of subdomains and interfaces is defined.
+    Parameters:
+        name: Variable name.
+        ndof: Number of dofs per grid element.
+            Valid keys are ``cells``, ``faces`` and ``nodes``.
+        subdomains (length=1): List containing a single grid.
+        interfaces (length=1): List containing a single mortar grid.
+        num_cells: Number of cells in the grid.
+            Only relevant if this is an interface variable.
 
     """
 
@@ -1196,20 +1208,6 @@ def __init__(
         previous_timestep: bool = False,
         previous_iteration: bool = False,
     ):
-        """Initiate an Ad representation of a variable associated with a grid or edge.
-
-        It is assumed that exactly one of subdomains and interfaces is defined.
-        Parameters:
-            name (str): Variable name.
-            ndof (dict): Number of dofs per grid element.
-            subdomains (optional list of pp.Grid ): List with length one containing a grid.
-            interfaces (optional list of pp.MortarGrid): List with length one containing
-                an interface.
-            num_cells (int): Number of cells in the grid. Only sued if the variable
-                is on an interface.
-
-        """
-
         self._name: str = name
         self._cells: int = ndof.get("cells", 0)
         self._faces: int = ndof.get("faces", 0)
@@ -1230,21 +1228,25 @@ def __init__(
             self._g = self.interfaces[0]
             self._is_edge_var = True
 
-        self.prev_time: bool = previous_timestep
-        self.prev_iter: bool = previous_iteration
-
         # The number of cells in the grid. Will only be used if grid_like is a tuple
         # that is, if this is a mortar variable
         self._num_cells = num_cells
 
-        self.id = next(self._ids)
         self._set_tree()
 
-    def size(self) -> int:
-        """Get the number of dofs for this grid.
+        self.prev_time: bool = previous_timestep
+        """Indicates whether the variable represents the state at the previous time step."""
+
+        self.prev_iter: bool = previous_iteration
+        """Indicates whether the variable represents the state at the previous iteration."""
+
+        self.id = next(self._ids)
+        """Unique identifier of this variable."""
 
+    def size(self) -> int:
+        """
         Returns:
-            int: Number of dofs.
+            The number of dofs of this variable on its grid.
 
         """
         if self._is_edge_var:
@@ -1260,10 +1262,10 @@ def size(self) -> int:
             )
 
     def previous_timestep(self) -> Variable:
-        """Return a representation of this variable on the previous time step.
-
+        """
         Returns:
-            Variable: A representation of this variable, with self.prev_time=True.
+            A representation of this variable at the previous time step,
+            with its ``prev_time`` attribute set to ``True``.
 
         """
         ndof = {"cells": self._cells, "faces": self._faces, "nodes": self._nodes}
@@ -1277,10 +1279,10 @@ def previous_timestep(self) -> Variable:
             )
 
     def previous_iteration(self) -> Variable:
-        """Return a representation of this variable on the previous time iteration.
-
+        """
         Returns:
-            Variable: A representation of this variable, with self.prev_iter=True.
+            A representation of this variable on the previous time iteration,
+            with its ``prev_iter`` attribute set to ``True``.
 
         """
         ndof = {"cells": self._cells, "faces": self._faces, "nodes": self._nodes}
@@ -1308,35 +1310,23 @@ def __repr__(self) -> str:
 
 
 class MergedVariable(Variable):
-    """Ad representation of a collection of variables that individually live on separate
-    subdomains of interfaces, but which it is useful to treat jointly.
+    """AD operator representing a collection of variables that individually live on separate
+    subdomains or interfaces, but treated jointly in the mixed-dimensional sense.
 
     Conversion of the variables into numerical value should be done with respect to the
-    state of an array; see the method evaluate().  Therefore, the class does not implement
-    a parse() method.
-
-    Attributes:
-        sub_vars (List of Variable): List of variable on different subdomains or interfaces.
-        id (int): Counter of all variables. Used to identify variables. Usage of this
-            term is not clear, it may change.
-        prev_iter (boolean): Whether the variable represents the state at the
-            previous iteration.
-        prev_time (boolean): Whether the variable represents the state at the
-            previous time step.
+    state of an array; see :meth:`Operator.evaluate`. Therefore, the MergedVariable does not
+    implement the method :meth:`Operator.parse`.
 
-        It is assumed that exactly one of subdomains and interfaces is defined.
+    Parameters:
+        variables: List of variables to be merged. Should all have the same name.
 
     """
 
     def __init__(self, variables: list[Variable]) -> None:
-        """Create a merged representation of variables.
-
-        Parameters:
-            variables (list of Variable): Variables to be merged. Should all have the
-                same name.
 
-        """
         self.sub_vars = variables
+        """List of single-grid variables which are merged into this merged variable,
+        passed at instantiation."""
 
         # Use counter from superclass to ensure unique Variable ids
         self.id = next(Variable._ids)
@@ -1367,19 +1357,18 @@ def __init__(self, variables: list[Variable]) -> None:
         self.prev_iter: bool = False
 
     def size(self) -> int:
-        """Get total size of the merged variable.
-
+        """
         Returns:
-            int: Total size of this merged variable.
+            Total size of this merged variable (sum of sizes of respective sub-variables).
 
         """
         return sum([v.size() for v in self.sub_vars])
 
     def previous_timestep(self) -> MergedVariable:
-        """Return a representation of this merged variable on the previous time step.
-
+        """
         Returns:
-            Variable: A representation of this variable, with self.prev_time=True.
+            A representation of this merged variable on the previous time
+            iteration, with its ``prev_iter`` attribute set to ``True``.
 
         """
         new_subs = [var.previous_timestep() for var in self.sub_vars]
@@ -1388,10 +1377,10 @@ def previous_timestep(self) -> MergedVariable:
         return new_var
 
     def previous_iteration(self) -> MergedVariable:
-        """Return a representation of this merged variable on the previous iteration.
-
+        """
         Returns:
-            Variable: A representation of this variable, with self.prev_iter=True.
+            A representation of this merged variable on the previous
+            iteration, with its ``prev_iter`` attribute set to ``True``
 
         """
         new_subs = [var.previous_iteration() for var in self.sub_vars]
@@ -1400,8 +1389,12 @@ def previous_iteration(self) -> MergedVariable:
         return new_var
 
     def copy(self) -> "MergedVariable":
-        # A shallow copy should be sufficient here; the attributes are not expected to
-        # change.
+        """
+        Returns:
+            A shallow copy should be sufficient here; the attributes are not expected to
+            change.
+
+        """
         return copy.deepcopy(self)
 
     def __repr__(self) -> str:
@@ -1435,9 +1428,17 @@ def __repr__(self) -> str:
 class Tree:
     """Simple implementation of a Tree class. Used to represent combinations of
     Ad operators.
+
+    References:
+        https://stackoverflow.com/questions/2358045/how-can-i-implement-a-tree-in-python
+
+
+    Parameters:
+        operation: See :data:`Operation`
+        children: List of children, either as Ad arrays or other :class:`Operator`.
+
     """
 
-    # https://stackoverflow.com/questions/2358045/how-can-i-implement-a-tree-in-python
     def __init__(
         self,
         operation: Operation,
@@ -1452,5 +1453,5 @@ def __init__(
                 self.add_child(child)
 
     def add_child(self, node: Union[Operator, Ad_array]) -> None:
-        #        assert isinstance(node, (Operator, "pp.ad.Operator"))
+        """Adds a child to this instance."""
         self.children.append(node)