Merge branch 'develop' into fix_tests_in_vscode

Author: IvarStefansson

src/porepy/examples/geothermal_reservoir.py

@@ -81,15 +81,14 @@ def bc_values_pressure(self, bg: pp.BoundaryGrid) -> np.ndarray:
         if self.is_well_grid(sd):
             well = self.well_network.wells[sd.tags["parent_well_index"]]
             well_tag = well.tags["well_name"]
-            protocol = self.well_protocols()[well_tag]
             # Find indices of the well boundary sides.
             domain_sides = self.domain_boundary_sides(bg)
             # The top of the domain is '.top' in 3d, '.north' in 2d.
             inds = domain_sides.top if self.nd == 3 else domain_sides.north
             # Set pressure values according to the well protocol.
             values[inds] = self.units.convert_units(
                 self.get_well_value(
-                    protocol["pressures"],
+                    self.well_protocols(well_tag, "pressures"),
                     self.time_manager.schedule,
                     self.time_manager.time,
                 ),
@@ -109,16 +108,14 @@ def bc_values_temperature(self, bg: pp.BoundaryGrid) -> np.ndarray:
         sd = bg.parent
         values = super().bc_values_temperature(bg)  # type: ignore[misc]
         if self.is_well_grid(sd):
-            # Retrieve well protocol.
             well_tag = self.well_names[sd.tags["parent_well_index"]]
-            protocol = self.well_protocols()[well_tag]
             # Find indices of the well boundary sides.
             domain_sides = self.domain_boundary_sides(bg)
             inds = domain_sides.top if self.nd == 3 else domain_sides.north
             # Set temperature values according to the well protocol.
             values[inds] = self.units.convert_units(
                 self.get_well_value(
-                    protocol["temperatures"],
+                    self.well_protocols(well_tag, "temperatures"),
                     self.time_manager.schedule,
                     self.time_manager.time,
                 ),
@@ -153,45 +150,32 @@ def get_well_value(
         else:
             return float(np.interp(current_time, times, values))
 
-    def well_protocols(self) -> dict[str, dict[str, NDArray[np.float64]]]:
-        """Dictionary mapping well tags to well protocols.
+    def well_protocols(self, well_tag: str, variable: str) -> NDArray[np.float64]:
+        """Return the time-dependent protocol array for a given well and variable.
+
+        The value is read from ``self.params`` under the key
+        ``"{well_tag}_{variable}"``. A scalar is broadcast to all schedule times; an
+        array must match schedule length.
+
+        Parameters:
+            well_tag: Name of the well (e.g. ``"injection_well"``).
+            variable: Protocol variable name (e.g. ``"pressures"``, ``"temperatures"``,
+                ``"mass_rates"``).
 
         Returns:
-            Dictionary with well protocols, each containing a dictionary with
-            time-dependent temperatures and pressures, with each value being an array of
-            size equal to the number of scheduled times in the time manager.
+            Array of protocol values, one entry per scheduled time point.
         """
         num_times = self.time_manager.schedule.size
-        protocols: dict[str, dict[str, NDArray[np.float64]]] = {}
-        # Construct protocols for each well.
-        for well_tag in self.well_names:
-            # Initialize protocol dictionary for the well.
-            protocols[well_tag] = {}
-            # Set values for temperatures and pressures.
-            for variable in ["temperatures", "pressures"]:
-                input_values = self.params.get(f"{well_tag}_{variable}", 0.0)
-                if isinstance(input_values, (float, int)):
-                    # Broadcast single value to all time steps for convenient user
-                    # definition of well protocols.
-                    values = np.full(num_times, input_values, dtype=float)
-
-                elif isinstance(input_values, (list, np.ndarray)):
-                    # Enforce array of float values.
-                    values = np.array(input_values, dtype=float)
-                    if values.size != num_times:
-                        raise ValueError(
-                            f"Well protocol for {well_tag} {variable} has size "
-                            f"{values.size}, expected {num_times}."
-                        )
-                else:
-                    raise TypeError(
-                        f"Well protocol for {well_tag} {variable} has unsupported "
-                        f"type {type(input_values)}."
-                    )
-                # Populate well dictionary for the current variable.
-                protocols[well_tag][variable] = values
-
-        return protocols
+        raw = self.params.get(f"{well_tag}_{variable}", 0.0)
+        if isinstance(raw, (int, float)):
+            return np.full(num_times, float(raw))
+        values = np.asarray(raw, dtype=float)
+        if values.size != num_times:
+            raise ValueError(
+                f"Protocol '{well_tag}_{variable}' has {values.size} entries, "
+                f"expected {num_times} (one per schedule point)."
+            )
+        return values
 
 
 class NeumannWellBCsFirstTimeInterval(pp.PorePyModel):

src/porepy/fracs/elliptic_fracture.py

@@ -25,10 +25,28 @@ class EllipticFracture(Fracture):
     The fracture is defined by its center position, major and minor axes, and its
     spatial orientation given by three rotation angles in radians.
 
+    Coordinate convention:
+        The global coordinate system is right-handed in 3D with +x = east, +y = north,
+        and +z = up.
+
+    Orientation convention:
+        The orientation of the fracture is described using a strike-dip representation.
+
+        The strike angle is defined as the azimuth, i.e., the angle measured clockwise
+        from the geographic north (y-axis) in the horizontal plane. This definition is
+        widely adopted in geosciences, although alternative notation conventions also
+        exist.
+
+        The strike-dip convention follows the right-hand rule in geology: when moving
+        along the strike direction, the fracture plane dips downward to the right.
+
+    References:
+        - [1] Fossen, H.: Structural Geology (2016).
+
     Example:
-        Fracture centered at ``[0, 1, 0]``, with a ratio of lengths of 2, rotation in
-        xy-plane of 45 degrees, and an incline of 30 degrees rotated around the x-axis,
-        due to the strike angle of 0 radians:
+        Fracture centered at ``[0, 1, 0]``, with a ratio of lengths of 2, a rotation
+        in the xy-plane of 45 degrees, and an incline of 30 degrees rotated around the
+        y-axis (north), due to the strike angle of 0 radians:
 
 
         >>> import numpy as np
@@ -43,19 +61,16 @@ class EllipticFracture(Fracture):
 
     Parameters:
         center: ``shape=(3, 1)``
-
             Center coordinates of fracture.
         major_axis: Length of major axis (radius-like, not diameter).
-        minor_axis: Length of minor axis.
-
-            There are no checks on whether the minor axis is less or equal to the major.
+        minor_axis: Length of minor axis. There are no checks on whether the minor axis
+            is less than or equal to the major axis.
         major_axis_angle: Rotation of the major axis from the x-axis in radians.
-            Measured before strike-dip rotation, see below.
-        strike_angle: Line of rotation for the dip. Given as angle in radians from the
-            x-direction.
+            Measured before the strike-dip rotation, see below.
+        strike_angle: Geological strike angle in radians, measured clockwise from the
+            y-axis (north) in the horizontal plane.
         dip_angle: Dip angle in radians, i.e., rotation around the strike direction.
         index: ``default=None``
-
             Index to be assigned to the fracture.
 
     """
@@ -97,9 +112,9 @@ def fracture_to_gmsh(self) -> int:
             dimTags, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, self.major_axis_angle
         )
 
-        # 3) Rotate around the strike direction by the dip angle.
-        strike_x = math.cos(self.strike_angle)
-        strike_y = math.sin(self.strike_angle)
+        # 3) Rotate around the geological strike direction by the dip angle.
+        strike_x = math.sin(self.strike_angle)
+        strike_y = math.cos(self.strike_angle)
         strike_z = 0.0
 
         gmsh.model.occ.rotate(

src/porepy/fracs/fracture_network.py

@@ -349,14 +349,16 @@ def _entity_on_domain_boundary(self, target_dim: int, ind: list[int]) -> bool:
         # having been split into multiple parts.
         domain_entities = gmsh.model.get_entities(self.nd)
         boundary_entities = gmsh.model.get_boundary(
-            [(self.nd, tag) for _, tag in domain_entities]
+            [(self.nd, tag) for _, tag in domain_entities], oriented=False
         )
         # Get hold of the boundary points of the entity to check.
         if target_dim == 0:
             boundary_points = [(target_dim, i) for i in ind]
         else:
             assert len(ind) == 1, "Only single entity indices are supported."
-            boundary_points = gmsh.model.get_boundary([(target_dim, ind[0])])
+            boundary_points = gmsh.model.get_boundary(
+                [(target_dim, ind[0])], oriented=False
+            )
 
         # For each boundary surface of the domain, compute the distance between the
         # entity and all boundary points to check if they are all zero.
@@ -405,7 +407,7 @@ def _insert_mesh_size_control_points(
         ### Get hold of entities representing fractures and boundaries.
         domain_entities = gmsh.model.get_entities(self.nd)
         boundaries = gmsh.model.get_boundary(
-            [(self.nd, tag) for _, tag in domain_entities]
+            [(self.nd, tag) for _, tag in domain_entities], oriented=False
         )
         fractures = [
             f for f in gmsh.model.get_entities(self.nd - 1) if f not in boundaries
@@ -433,7 +435,9 @@ def _insert_mesh_size_control_points(
         # Take note of the boundary points of all entities, to avoid inserting points
         # there (doing so may confuse Gmsh).
         for ent in entities:
-            bp = gmsh.model.get_boundary([(self.nd - 1, ent)], recursive=True)
+            bp = gmsh.model.get_boundary(
+                [(self.nd - 1, ent)], recursive=True, oriented=False
+            )
             for b in bp:
                 if b[0] != 0:
                     continue

src/porepy/fracs/fracture_network_2d.py

@@ -384,7 +384,10 @@ def _impose_boundary_process_intersections(
             boundary_tags = []
         else:
             boundary_tags = [
-                t for _, t in gmsh.model.get_boundary([(self.nd, domain_tag)])
+                t
+                for _, t in gmsh.model.get_boundary(
+                    [(self.nd, domain_tag)], oriented=False
+                )
             ]
 
         # Mapping from the new fracture tags (gmsh assigned) to the input fractures.
@@ -441,7 +444,7 @@ def _impose_boundary_process_intersections(
                     updated_mesh_size_points[segment[1]] = mesh_size_points[
                         old_gmsh_tag
                     ]
-                pt_index = gmsh.model.get_boundary([segment])
+                pt_index = gmsh.model.get_boundary([segment], oriented=False)
 
                 if fi not in constraints:
                     # If this is not a constraint, collect the boundary points for
@@ -507,7 +510,7 @@ def _set_mesh_size_fields(
         ### Get hold of lines representing fractures and boundaries.
         domain_entities = gmsh.model.get_entities(2)
         boundaries = gmsh.model.get_boundary(
-            [(self.nd, tag) for _, tag in domain_entities]
+            [(self.nd, tag) for _, tag in domain_entities], oriented=False
         )
 
         line_tags = set(tag for _, tag in gmsh.model.getEntities(self.nd - 1))
@@ -540,7 +543,9 @@ def _set_mesh_size_fields(
             end_points = np.array(
                 [
                     gmsh.model.occ.get_bounding_box(0, p[1])[:3]
-                    for p in gmsh.model.get_boundary([(1, line)], combined=False)
+                    for p in gmsh.model.get_boundary(
+                        [(1, line)], combined=False, oriented=False
+                    )
                 ]
             ).T
             length = np.linalg.norm(end_points[:, 1] - end_points[:, 0])

src/porepy/fracs/fracture_network_3d.py

@@ -419,10 +419,10 @@ def _impose_boundary_process_intersections(
             num_orig_subfrac.append(len(frac))
 
             for sfi, sub_frac in enumerate(frac):
-                bounding_lines = gmsh.model.get_boundary([sub_frac])
+                bounding_lines = gmsh.model.get_boundary([sub_frac], oriented=False)
                 bounding_points = []
                 for line in bounding_lines:
-                    bounding_points += gmsh.model.get_boundary([line])
+                    bounding_points += gmsh.model.get_boundary([line], oriented=False)
 
                 if len(bounding_points) == 0:
                     # This is most likely a disc fracture, which has no bounding
@@ -636,7 +636,7 @@ def _impose_boundary_process_intersections(
                 # At most one of the parents was not a constraint. This line should not
                 # produce a point.
                 continue
-            for bp in gmsh.model.get_boundary([(1, line)]):
+            for bp in gmsh.model.get_boundary([(1, line)], oriented=False):
                 points_of_intersection_lines.append(bp[1])
 
         num_point_occ = np.bincount(points_of_intersection_lines)
@@ -736,7 +736,9 @@ def _set_mesh_size_fields(
         ### Get hold of lines representing fractures and boundaries.
         domain_entities = gmsh.model.get_entities(nd)
         # Get the boundaries.
-        boundaries = gmsh.model.get_boundary([(nd, tag) for _, tag in domain_entities])
+        boundaries = gmsh.model.get_boundary(
+            [(nd, tag) for _, tag in domain_entities], oriented=False
+        )
 
         surface_tags = set(tag for _, tag in gmsh.model.get_entities(nd - 1))
         boundary_tags = set(tag for _, tag in boundaries)