Issue210 borefield class

pygfunction/__init__.py

@@ -6,4 +6,5 @@
 from . import media
 from . import networks
 from . import pipes
+from . import solvers
 from . import utilities

pygfunction/borefield.py

@@ -1,13 +1,15 @@
 # -*- coding: utf-8 -*-
+from collections.abc import Callable
 from typing import Union, List, Dict, Tuple
 from typing_extensions import Self     # for compatibility with Python <= 3.10
 
 import matplotlib.pyplot as plt
 from matplotlib.figure import Figure
 import numpy as np
 import numpy.typing as npt
+from scipy.spatial.distance import pdist, squareform
 
-from .boreholes import Borehole
+from .boreholes import Borehole, _EquivalentBorehole
 from .utilities import _initialize_figure, _format_axes, _format_axes_3d
 
 
@@ -45,10 +47,9 @@ def __init__(
             self, H: npt.ArrayLike, D: npt.ArrayLike, r_b: npt.ArrayLike,
             x: npt.ArrayLike, y: npt.ArrayLike, tilt: npt.ArrayLike = 0.,
             orientation: npt.ArrayLike = 0.):
-        # Convert x and y coordinates to arrays
-        x = np.atleast_1d(x)
-        y = np.atleast_1d(y)
-        self.nBoreholes = np.maximum(len(x), len(y))
+        self.nBoreholes = np.max(
+            [len(np.atleast_1d(var))
+             for var in (H, D, r_b, x, y, tilt, orientation)])
 
         # Broadcast all variables to arrays of length `nBoreholes`
         self.H = np.broadcast_to(H, self.nBoreholes)
@@ -62,6 +63,8 @@ def __init__(
         self._is_tilted = np.broadcast_to(
             np.greater(np.abs(tilt), 1e-6),
             self.nBoreholes)
+        self._all_tilted = np.all(self._is_tilted)
+        self._all_vertical = not np.any(self._is_tilted)
         # Vertical boreholes default to an orientation of zero
         if not np.any(self._is_tilted):
             self.orientation = np.broadcast_to(0., self.nBoreholes)
@@ -109,6 +112,99 @@ def __ne__(
         check = not self == other_field
         return check
 
+    @property
+    def all_tilted(self) -> np.ndarray:
+        """Return a bool. True if the all boreholes are inclined."""
+        return self._all_tilted
+
+    @property
+    def all_vertical(self) -> np.ndarray:
+        """Return an bool. True if all boreholes are vertical."""
+        return self._all_vertical
+
+    @property
+    def is_tilted(self) -> np.ndarray:
+        """Return an boolean array. True if the corresponding borehole is inclined."""
+        return self._is_tilted
+
+    @property
+    def is_vertical(self) -> np.ndarray:
+        """Return an boolean array. True if the corresponding borehole is vertical."""
+        return np.logical_not(self._is_tilted)
+
+    def distance(
+            self, other_field: Union[Borehole, Self, List[Borehole]],
+            outer: bool = True
+            ) -> np.ndarray:
+        """
+        Evaluate horizontal distances to boreholes of another borefield.
+
+        Parameters
+        ----------
+        other_field : Borehole or Borefield object
+            Borefield with which to evaluate distances.
+        outer : , optional
+            True if the horizontal distance is to be evaluated for all
+            boreholes of the borefield onto all boreholes of the other
+            borefield to return a (nBoreholes_other_field, nBoreholes,)
+            array. If false, the horizontal distance is evaluated pairwise
+            between boreholes of the borefield and boreholes of the other
+            borefield. The numbers of boreholes should be the same
+            (i.e. nBoreholes == nBoreholes_other_field) and a (nBoreholes,)
+            array is returned.
+            Default is True.
+
+        Returns
+        -------
+        dis : array, shape (nBoreholes_other_field, nBoreholes,) or (nBoreholes,)
+            Horizontal distances between boreholes (in meters).
+
+        """
+        if outer:
+            dis = np.maximum(
+                np.sqrt(
+                    np.subtract.outer(other_field.x, self.x)**2
+                    + np.subtract.outer(other_field.y, self.y)**2
+                    ),
+                self.r_b)
+        else:
+            dis = np.maximum(
+                np.sqrt(
+                    (other_field.x - self.x)**2 + (other_field.y - self.y)**2),
+                self.r_b)
+        return dis
+
+    def distance_to_self(
+            self, outer: bool = True) -> np.ndarray:
+        """
+        Evaluate horizontal distances between boreholes within the borefield.
+
+        Parameters
+        ----------
+        outer : , optional
+            True if the horizontal distance is to be evaluated for all
+            boreholes of the borefield onto all of the boreholesto return a
+            (nBoreholes, nBoreholes,) array. If false, a condensed distance
+            vector corresponding to the upper triangle of the distance matrix
+            is return (excluding the diagonal elements).
+            Default is True.
+
+        Returns
+        -------
+        dis : array, shape (nBoreholes, nBoreholes,) or (nBoreholes,)
+            Horizontal distances between boreholes (in meters).
+
+        """
+        # Condensed vector of the distances between boreholes
+        dis = pdist(np.stack((self.x, self.y), axis=1))
+        if outer:
+            # Convert to a square matrix
+            dis = squareform(dis)
+            # Overwrite the diagonal with the borehole radii
+            i, j = np.diag_indices(len(self))
+            dis[i, j] = self.r_b
+        return dis
+
     def evaluate_g_function(
             self,
             alpha: float,
@@ -306,6 +402,77 @@ def evaluate_g_function(
 
         return gfunc.gFunc
 
+    def segments(
+            self, nSegments: Union[int, npt.ArrayLike],
+            segment_ratios: Union[
+                None, npt.ArrayLike, List[npt.ArrayLike], Callable[[int], npt.ArrayLike]
+                ] = None) -> Self:
+        """
+        Split boreholes in the bore field into segments.
+
+        Parameters
+        ----------
+        nSegments : int or (nBoreholes,) array
+            Number of segments per borehole.
+        segment_ratios : array or list of arrays, optional
+            Ratio of the borehole length represented by each segment. The sum
+            of ratios must be equal to 1. The shape of the array is of
+            (nSegments,). If all boreholes have the same number of segments
+            and the same ratios, a single (nSegments,) array can be provided.
+            Otherwise, a list of (nSegments_i,) arrays with the ratios
+            associated with each borehole must be provided. If
+            segment_ratios==None, segments of equal lengths are considered.
+            Default is None.
+
+        Returns
+        -------
+        segments : Borefield object
+            Segments in the bore field.
+
+        Examples
+        --------
+        >>> borefield = gt.borefield.Borefield.rectangle_field(
+            N_1=2, N_2=3, B_1 = 7.5, B_2=7.5, H=150., D=4., r_b=0.075)
+        >>> borefield.segments(5)
+
+        """
+        nSegments = np.broadcast_to(nSegments, self.nBoreholes)
+        if callable(segment_ratios):
+            segment_ratios = [segment_ratios(nSeg) for nSeg in nSegments]
+        if segment_ratios is None:
+            # Local coordinates of the top-edges of the segments
+            xi = np.concatenate(
+                [np.arange(nSeg_i) / nSeg_i for nSeg_i in nSegments])
+            # Segment ratios
+            segment_ratios = np.concatenate(
+                [np.full(nSeg_i, 1. / nSeg_i) for nSeg_i in nSegments])
+        elif isinstance(segment_ratios, list):
+            # Local coordinates of the top-edges of the segments
+            xi = np.concatenate(
+                [np.cumsum(np.concatenate([[0.], seg_rat_i[:-1]]))
+                 for seg_rat_i in segment_ratios])
+            # Segment ratios
+            segment_ratios = np.concatenate(segment_ratios)
+        else:
+            # Local coordinates of the top-edges of the segments
+            xi = np.tile(
+                np.cumsum(np.concatenate([[0.], segment_ratios[:-1]])),
+                self.nBoreholes)
+            # Segment ratios
+            segment_ratios = np.tile(segment_ratios, self.nBoreholes)
+        xi = xi * np.repeat(self.H, nSegments)
+
+        # Geometrical parameters of the segments
+        H = np.repeat(self.H, nSegments) * segment_ratios
+        r_b = np.repeat(self.r_b, nSegments)
+        tilt = np.repeat(self.tilt, nSegments)
+        orientation = np.repeat(self.orientation, nSegments)
+        D = np.repeat(self.D, nSegments) + xi * np.cos(tilt)
+        x = np.repeat(self.x, nSegments) + xi * np.sin(tilt) * np.cos(orientation)
+        y = np.repeat(self.y, nSegments) + xi * np.sin(tilt) * np.sin(orientation)
+
+        return Borefield(H, D, r_b, x, y, tilt=tilt, orientation=orientation)
+
     def visualize_field(
             self, viewTop: bool = True, view3D: bool = True,
             labels: bool = True, showTilt: bool = True) -> Figure:
@@ -1098,3 +1265,97 @@ def circle_field(
         # Create the bore field
         borefield = cls(H, D, r_b, x, y, tilt=tilt, orientation=orientation)
         return borefield
+
+
+class _EquivalentBorefield(object):
+
+    def __init__(
+            self, H: npt.ArrayLike, D: npt.ArrayLike, r_b: npt.ArrayLike,
+            x: List[npt.ArrayLike], y: List[npt.ArrayLike],
+            tilt: npt.ArrayLike, orientation: List[npt.ArrayLike]):
+        self.H = H
+        self.D = D
+        self.r_b = r_b
+        self.x = x
+        self.y = y
+        self.tilt = tilt
+        self.orientation = orientation
+        self.nBoreholes = len(x)
+
+        # Identify tilted boreholes
+        self._is_tilted = np.broadcast_to(
+            np.greater(np.abs(tilt), 1e-6),
+            self.nBoreholes)
+        # Vertical boreholes default to an orientation of zero
+        if not np.any(self._is_tilted):
+            self.orientation = [np.broadcast_to(0., len(x_i)) for x_i in x]
+        else:
+            self.orientation = [
+                np.multiply(orientation_i, _is_tilted_i)
+                for (orientation_i, _is_tilted_i)
+                in zip(self.orientation, self._is_tilted)]
+
+    def __getitem__(self, key):
+        if isinstance(key, (int, np.integer)):
+            # Returns a borehole object if only one borehole is indexed
+            output_class = _EquivalentBorehole
+        else:
+            # Returns a _EquivalentBorehole object for slices and lists of
+            # indexes
+            output_class = _EquivalentBorefield
+        return output_class(
+            self.H[key], self.D[key], self.r_b[key], self.x[key], self.y[key],
+            tilt=self.tilt[key], orientation=self.orientation[key])
+
+    def __len__(self) -> int:
+        """Return the number of boreholes."""
+        return self.nBoreholes
+
+    def segments(
+            self, nSegments: int,
+            segment_ratios: Union[
+                None, npt.ArrayLike, Callable[[int], npt.ArrayLike]
+                ] = None) -> Self:
+        """
+        Split boreholes in the bore field into segments.
+
+        Parameters
+        ----------
+        nSegments : int or (nBoreholes,) array
+            Number of segments per borehole.
+        segment_ratios : array or list of arrays, optional
+            Ratio of the borehole length represented by each segment. The sum
+            of ratios must be equal to 1. The shape of the array is of
+            (nSegments,). If all boreholes have the same number of segments
+            and the same ratios, a single (nSegments,) array can be provided.
+            Otherwise, a list of (nSegments_i,) arrays with the ratios
+            associated with each borehole must be provided. If
+            segment_ratios==None, segments of equal lengths are considered.
+            Default is None.
+
+        Returns
+        -------
+        segments : Borefield object
+            Segments in the bore field.
+
+        Examples
+        --------
+        >>> borefield = gt.borefield.Borefield.rectangle_field(
+            N_1=2, N_2=3, B_1 = 7.5, B_2=7.5, H=150., D=4., r_b=0.075)
+        >>> borefield.segments(5)
+
+        """
+        return _EquivalentBorefield.from_equivalent_boreholes(
+            [seg for b in self for seg in b.segments(nSegments, segment_ratios=segment_ratios)
+             ])
+
+    @classmethod
+    def from_equivalent_boreholes(cls, equivalent_boreholes):
+        H = np.array([b.H for b in equivalent_boreholes])
+        D = np.array([b.D for b in equivalent_boreholes])
+        r_b = np.array([b.r_b for b in equivalent_boreholes])
+        x = [b.x for b in equivalent_boreholes]
+        y = [b.y for b in equivalent_boreholes]
+        tilt = np.array([b.tilt for b in equivalent_boreholes])
+        orientation = [b.orientation for b in equivalent_boreholes]
+        return cls(H, D, r_b, x, y, tilt, orientation)