Repository: MassimoCimmino/pygfunction (language: python)

Issue:
Transition unit tests to pytest
This issue is to replace the `unittest` with `pytest`.

`pytest` has a lot of useful features, such as parametrized tests to simplify the implementation of tests.


Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4da0c751..346d76d6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -20,6 +20,7 @@
 * [Issue 172](https://github.com/MassimoCimmino/pygfunction/issues/172) - Refactored reports of calculation time to use `time.perf_counter()` instead of `time.time()`.
 * [Issue 173](https://github.com/MassimoCimmino/pygfunction/issues/173) - Refactored strings into f-strings instead of using `str.format()`.
 * [Issue 177](https://github.com/MassimoCimmino/pygfunction/issues/177) - Converted `setup.py` script to `setup.cfg` and `pyproject.toml` files. This is motivated by [PEP518](https://www.python.org/dev/peps/pep-0518/) and [PEP621](https://www.python.org/dev/peps/pep-0621/).
+* [Issue 179](https://github.com/MassimoCimmino/pygfunction/issues/179) - Refactored tests to use the `pytest` package instead of `unittests`.
 * [Issue 180](https://github.com/MassimoCimmino/pygfunction/issues/180) - Configured `tox` and github actions for continuous integration.
 
 ### Bug fixes


Test patch:
diff --git a/tests/boreholes_test.py b/tests/boreholes_test.py
index b02ed1af..ebc3c5e9 100644
--- a/tests/boreholes_test.py
+++ b/tests/boreholes_test.py
@@ -1,353 +1,219 @@
 # -*- coding: utf-8 -*-
 """ Test suite for boreholes module.
 """
-import unittest
-
 import numpy as np
-from scipy.constants import pi
-
-
-class TestBorehole(unittest.TestCase):
-    """ Test cases for Borehole() class.
-    """
-
-    def setUp(self):
-        from pygfunction import boreholes
-        self.H = 150.
-        self.D = 2.5
-        self.r_b = 0.075
-        self.tilt = pi/6.
-        self.orientation = pi
-        self.x = 2.
-        self.y = 4.
-        self.bore = boreholes.Borehole(self.H, self.D, self.r_b,
-                                       self.x, self.y, tilt=self.tilt,
-                                       orientation=self.orientation)
-        self.x2 = 5.
-        self.y2 = 8.
-        self.bore2 = boreholes.Borehole(self.H, self.D, self.r_b,
-                                        self.x2, self.y2, tilt=self.tilt,
-                                        orientation=self.orientation)
-
-    def test_inputs(self):
-        """ Tests that all inputs return their value.
-        """
-        self.assertEqual(self.bore.H, self.H,
-                         msg='borehole.H returns incorrect value.')
-        self.assertEqual(self.bore.D, self.D,
-                         msg='borehole.D returns incorrect value.')
-        self.assertEqual(self.bore.r_b, self.r_b,
-                         msg='borehole.r_b returns incorrect value.')
-        self.assertEqual(self.bore.x, self.x,
-                         msg='borehole.x returns incorrect value.')
-        self.assertEqual(self.bore.y, self.y,
-                         msg='borehole.y returns incorrect value.')
-        self.assertEqual(self.bore.tilt, self.tilt,
-                         msg='borehole.tilt returns incorrect value.')
-        self.assertEqual(self.bore.orientation, self.orientation,
-                         msg='borehole.orientation returns incorrect value.')
-
-    def test_position(self):
-        """ Tests position() class method.
-        """
-        pos = self.bore.position()
-        self.assertEqual(pos, (self.x, self.y),
-                         msg='borehole.position() returns incorrect value.')
-
-    def test_distance_with_self(self, abs_tol=1.0e-6):
-        """ Tests distance(borehole) class method returns radius when distance
-            is less than the borehole radius.
-        """
-        dis = self.bore.distance(self.bore)
-        self.assertAlmostEqual(dis, self.r_b, delta=abs_tol,
-                         msg=('borehole.distance() should return the borehole '
-                              'radius when the distance is less than r_b.'))
-
-    def test_distance_with_another(self, abs_tol=1.0e-6):
-        """ Tests distance(borehole2) class method returns separating distance.
-        """
-        dis = self.bore.distance(self.bore2)
-        correct_dis = np.sqrt((self.x2 - self.x)**2 + (self.y2 - self.y)**2)
-        self.assertAlmostEqual(dis, correct_dis, delta=abs_tol,
-                         msg=('borehole.distance() does not evaluate correct '
-                              'separating distance with target borehole.'))
-
-
-class TestBoreFields_RectangleField(unittest.TestCase):
-    """ Test cases for boreholes.rectangle_field().
-    """
-
-    def setUp(self):
-        self.H = 150.
-        self.D = 2.5
-        self.r_b = 0.075
-        self.B_1 = 5.
-        self.B_2 = 6.
-
-    def test_rectangular_1_borehole(self):
-        """ Tests construction of rectangular field with one borehole.
-        """
-        from pygfunction import boreholes
-        N_1 = 1
-        N_2 = 1
-        boreField = boreholes.rectangle_field(N_1, N_2, self.B_1, self.B_2,
-                                              self.H, self.D, self.r_b)
-        self.assertEqual(len(boreField), 1,
-                         msg=('Incorrect number of boreholes in '
-                              'rectangle_field with N1=1 and N2=1.'))
-
-    def test_rectangular_1_row(self):
-        """ Tests construction of rectangular field with one row.
-        """
-        from pygfunction import boreholes
-        N_1 = 4
-        N_2 = 1
-        boreField = boreholes.rectangle_field(N_1, N_2, self.B_1, self.B_2,
-                                              self.H, self.D, self.r_b)
-        self.assertEqual(len(boreField), N_1,
-                         msg=('Incorrect number of boreholes in '
-                              'rectangle_field with N2=1.'))
-
-    def test_rectangular_1_column(self):
-        """ Tests construction of rectangular field with one column.
-        """
-        from pygfunction import boreholes
-        N_1 = 1
-        N_2 = 4
-        boreField = boreholes.rectangle_field(N_1, N_2, self.B_1, self.B_2,
-                                              self.H, self.D, self.r_b)
-        self.assertEqual(len(boreField), N_2,
-                         msg=('Incorrect number of boreholes in '
-                              'rectangle_field with N1=1.'))
-
-    def test_rectangular(self):
-        """ Tests construction of rectangular field with multiple rows/columns.
-        """
-        from pygfunction import boreholes
-        N_1 = 4
-        N_2 = 5
-        boreField = boreholes.rectangle_field(N_1, N_2, self.B_1, self.B_2,
-                                              self.H, self.D, self.r_b)
-        self.assertEqual(len(boreField), N_1*N_2,
-                         msg=('Incorrect number of boreholes in '
-                              'rectangle_field.'))
-
-
-class TestBoreFields_LShapedField(unittest.TestCase):
-    """ Test cases for boreholes.L_shaped_field().
-    """
-
-    def setUp(self):
-        self.H = 150.
-        self.D = 2.5
-        self.r_b = 0.075
-        self.B_1 = 5.
-        self.B_2 = 6.
-
-    def test_L_shaped_1_borehole(self):
-        """ Tests construction of L-shaped field with one borehole.
-        """
-        from pygfunction import boreholes
-        N_1 = 1
-        N_2 = 1
-        boreField = boreholes.L_shaped_field(N_1, N_2, self.B_1, self.B_2,
-                                             self.H, self.D, self.r_b)
-        self.assertEqual(len(boreField), 1,
-                         msg=('Incorrect number of boreholes in '
-                              'L_shaped_field with N1=1 and N2=1.'))
-
-    def test_L_shaped_1_row(self):
-        """ Tests construction of L-shaped field with one borehole.
-        """
-        from pygfunction import boreholes
-        N_1 = 5
-        N_2 = 1
-        boreField = boreholes.L_shaped_field(N_1, N_2, self.B_1, self.B_2,
-                                             self.H, self.D, self.r_b)
-        self.assertEqual(len(boreField), N_1,
-                         msg=('Incorrect number of boreholes in '
-                              'L_shaped_field with N2=1.'))
-
-    def test_L_shaped_1_column(self):
-        """ Tests construction of L-shaped field with one borehole.
-        """
-        from pygfunction import boreholes
-        N_1 = 1
-        N_2 = 5
-        boreField = boreholes.L_shaped_field(N_1, N_2, self.B_1, self.B_2,
-                                             self.H, self.D, self.r_b)
-        self.assertEqual(len(boreField), N_2,
-                         msg=('Incorrect number of boreholes in '
-                              'L_shaped_field with N1=1.'))
-
-    def test_L_shaped(self):
-        """ Tests construction of L-shaped field with with multiple
-            rows/columns.
-        """
-        from pygfunction import boreholes
-        N_1 = 3
-        N_2 = 5
-        boreField = boreholes.L_shaped_field(N_1, N_2, self.B_1, self.B_2,
-                                             self.H, self.D, self.r_b)
-        self.assertEqual(len(boreField), N_1+N_2-1,
-                         msg=('Incorrect number of boreholes in '
-                              'L_shaped_field.'))
-
-
-class TestBoreFields_UShapedField(unittest.TestCase):
-    """ Test cases for boreholes.U_shaped_field().
-    """
-
-    def setUp(self):
-        self.H = 150.
-        self.D = 2.5
-        self.r_b = 0.075
-        self.B_1 = 5.
-        self.B_2 = 6.
-
-    def test_U_shaped_1_borehole(self):
-        """ Tests construction of L-shaped field with one borehole.
-        """
-        from pygfunction import boreholes
-        N_1 = 1
-        N_2 = 1
-        boreField = boreholes.U_shaped_field(N_1, N_2, self.B_1, self.B_2,
-                                             self.H, self.D, self.r_b)
-        self.assertEqual(len(boreField), 1,
-                         msg=('Incorrect number of boreholes in '
-                              'U_shaped_field with N1=1 and N2=1.'))
-
-    def test_U_shaped_1_row(self):
-        """ Tests construction of L-shaped field with one borehole.
-        """
-        from pygfunction import boreholes
-        N_1 = 5
-        N_2 = 1
-        boreField = boreholes.U_shaped_field(N_1, N_2, self.B_1, self.B_2,
-                                             self.H, self.D, self.r_b)
-        self.assertEqual(len(boreField), N_1,
-                         msg=('Incorrect number of boreholes in '
-                              'U_shaped_field with N2=1.'))
-
-    def test_U_shaped_1_column(self):
-        """ Tests construction of U-shaped field with one borehole.
-        """
-        from pygfunction import boreholes
-        N_1 = 1
-        N_2 = 5
-        boreField = boreholes.U_shaped_field(N_1, N_2, self.B_1, self.B_2,
-                                             self.H, self.D, self.r_b)
-        self.assertEqual(len(boreField), N_2,
-                         msg=('Incorrect number of boreholes in '
-                              'U_shaped_field with N1=1.'))
-
-    def test_U_shaped(self):
-        """ Tests construction of U-shaped field with with multiple
-            rows/columns.
-        """
-        from pygfunction import boreholes
-        N_1 = 3
-        N_2 = 5
-        boreField = boreholes.U_shaped_field(N_1, N_2, self.B_1, self.B_2,
-                                             self.H, self.D, self.r_b)
-        self.assertEqual(len(boreField), N_1+2*N_2-2,
-                         msg=('Incorrect number of boreholes in '
-                              'U_shaped_field.'))
-
-
-class TestBoreFields_BoxShapedField(unittest.TestCase):
-    """ Test cases for boreholes.box_shaped_field().
-    """
-
-    def setUp(self):
-        self.H = 150.
-        self.D = 2.5
-        self.r_b = 0.075
-        self.B_1 = 5.
-        self.B_2 = 6.
-
-    def test_box_shaped_1_borehole(self):
-        """ Tests construction of box-shaped field with one borehole.
-        """
-        from pygfunction import boreholes
-        N_1 = 1
-        N_2 = 1
-        boreField = boreholes.box_shaped_field(N_1, N_2, self.B_1, self.B_2,
-                                               self.H, self.D, self.r_b)
-        self.assertEqual(len(boreField), 1,
-                         msg=('Incorrect number of boreholes in '
-                              'box_shaped_field with N1=1 and N2=1.'))
-
-    def test_box_shaped_1_row(self):
-        """ Tests construction of L-shaped field with one borehole.
-        """
-        from pygfunction import boreholes
-        N_1 = 5
-        N_2 = 1
-        boreField = boreholes.box_shaped_field(N_1, N_2, self.B_1, self.B_2,
-                                             self.H, self.D, self.r_b)
-        self.assertEqual(len(boreField), N_1,
-                         msg=('Incorrect number of boreholes in '
-                              'box_shaped_field with N2=1.'))
-
-    def test_box_shaped_1_column(self):
-        """ Tests construction of box-shaped field with one borehole.
-        """
-        from pygfunction import boreholes
-        N_1 = 1
-        N_2 = 5
-        boreField = boreholes.box_shaped_field(N_1, N_2, self.B_1, self.B_2,
-                                               self.H, self.D, self.r_b)
-        self.assertEqual(len(boreField), N_2,
-                         msg=('Incorrect number of boreholes in '
-                              'box_shaped_field with N1=1.'))
-
-    def test_box_shaped(self):
-        """ Tests construction of box-shaped field with with multiple
-            rows/columns.
-        """
-        from pygfunction import boreholes
-        N_1 = 3
-        N_2 = 5
-        boreField = boreholes.box_shaped_field(N_1, N_2, self.B_1, self.B_2,
-                                             self.H, self.D, self.r_b)
-        self.assertEqual(len(boreField), 2*N_1+2*N_2-4,
-                         msg=('Incorrect number of boreholes in '
-                              'box_shaped_field.'))
-
-
-class TestBoreFields_CircleField(unittest.TestCase):
-    """ Test cases for boreholes.circle_field().
-    """
-
-    def setUp(self):
-        self.H = 150.
-        self.D = 2.5
-        self.r_b = 0.075
-        self.R = 5.
-
-    def test_circle_1_borehole(self):
-        """ Tests construction of circle field with one borehole.
-        """
-        from pygfunction import boreholes
-        N = 1
-        boreField = boreholes.circle_field(N, self.R, self.H, self.D, self.r_b)
-        self.assertEqual(len(boreField), 1,
-                         msg=('Incorrect number of boreholes in '
-                              'circle_field with N1=1 and N2=1.'))
-
-    def test_circle(self):
-        """ Tests construction of circle field with with multiple boreholes.
-        """
-        from pygfunction import boreholes
-        N = 8
-        boreField = boreholes.circle_field(N, self.R, self.H, self.D, self.r_b)
-        self.assertEqual(len(boreField), N,
-                         msg=('Incorrect number of boreholes in '
-                              'circle_field.'))
-
-
-if __name__ == '__main__' and __package__ is None:
-    from os import sys, path
-    sys.path.append(path.dirname(path.dirname(path.abspath(__file__))))
-    unittest.main()
+import pytest
+
+import pygfunction as gt
+
+
+# =============================================================================
+# Test Borehole class
+# =============================================================================
+# Test the initialization of the Borehole class
+def test_borehole_init():
+    H = 150.                    # Borehole length [m]
+    D = 4.                      # Borehole buried depth [m]
+    r_b = 0.075                 # Borehole radius [m]
+    x = 0.                      # Borehole x-position [m]
+    y = 0.                      # Borehole y-position [m]
+    tilt = np.pi / 15           # Borehole tilt [rad]
+    orientation = np.pi / 3     # Borehole orientation [rad]
+    # Initialize borehole object
+    borehole = gt.boreholes.Borehole(
+        H, D, r_b, x, y, tilt=tilt, orientation=orientation)
+    assert np.all(
+        [H == borehole.H,
+         D == borehole.D,
+         r_b == borehole.r_b,
+         x == borehole.x,
+         y == borehole.y,
+         orientation == borehole.orientation,
+         tilt == borehole.tilt
+         ])
+
+
+# Test Borehole.distance
+@pytest.mark.parametrize("borehole1, borehole2", [
+        # Same borehole
+        ('single_borehole', 'single_borehole'),
+        # Two boreholes
+        ('single_borehole', 'single_borehole_short'),
+    ])
+def test_borehole_distance(borehole1, borehole2, request):
+    # Extract boreholes from fixtures
+    b1 = request.getfixturevalue(borehole1)[0]
+    b2 = request.getfixturevalue(borehole2)[0]
+    # Expected distance [m]
+    dis = max(b1.r_b, np.sqrt((b1.x - b2.x)**2 + (b1.y - b2.y)**2))
+    assert np.isclose(dis, b1.distance(b2))
+
+
+# Test Borehole.position
+@pytest.mark.parametrize("borehole", [
+        ('single_borehole'),
+        ('single_borehole_short'),
+    ])
+def test_borehole_position(borehole, request):
+    # Extract borehole from fixture
+    b = request.getfixturevalue(borehole)[0]
+    # Evaluate position ([m], [m])
+    (x, y) = b.position()
+    assert x == b.x and y == b.y
+
+
+# =============================================================================
+# Test functions
+# =============================================================================
+# Test rectangular_field
+@pytest.mark.parametrize("N_1, N_2, B_1, B_2", [
+        (1, 1, 5., 5.),     # 1 by 1
+        (2, 1, 5., 5.),     # 2 by 1
+        (1, 2, 5., 5.),     # 1 by 2
+        (2, 2, 5., 7.5),    # 2 by 2 (different x/y spacings)
+        (10, 9, 7.5, 5.),   # 10 by 9 (different x/y spacings)
+    ])
+def test_rectangular_field(N_1, N_2, B_1, B_2):
+    H = 150.        # Borehole length [m]
+    D = 4.          # Borehole buried depth [m]
+    r_b = 0.075     # Borehole radius [m]
+    # Generate the bore field
+    field = gt.boreholes.rectangle_field(N_1, N_2, B_1, B_2, H, D, r_b)
+    # Evaluate the borehole to borehole distances
+    x = np.array([b.x for b in field])
+    y = np.array([b.y for b in field])
+    dis = np.sqrt(
+        np.subtract.outer(x, x)**2 + np.subtract.outer(y, y)**2)[
+            ~np.eye(len(field), dtype=bool)]
+    assert np.all(
+        [len(field) == N_1 * N_2,
+         np.allclose(H, [b.H for b in field]),
+         np.allclose(D, [b.D for b in field]),
+         np.allclose(r_b, [b.r_b for b in field]),
+         len(field) == 1 or np.isclose(np.min(dis), min(B_1, B_2)),
+         ])
+
+
+# Test L_shaped_field
+@pytest.mark.parametrize("N_1, N_2, B_1, B_2", [
+        (1, 1, 5., 5.),     # 1 by 1
+        (2, 1, 5., 5.),     # 2 by 1
+        (1, 2, 5., 5.),     # 1 by 2
+        (2, 2, 5., 7.5),    # 2 by 2 (different x/y spacings)
+        (10, 9, 7.5, 5.),   # 10 by 9 (different x/y spacings)
+    ])
+def test_L_shaped_field(N_1, N_2, B_1, B_2):
+    H = 150.        # Borehole length [m]
+    D = 4.          # Borehole buried depth [m]
+    r_b = 0.075     # Borehole radius [m]
+    # Generate the bore field
+    field = gt.boreholes.L_shaped_field(N_1, N_2, B_1, B_2, H, D, r_b)
+    # Evaluate the borehole to borehole distances
+    x = np.array([b.x for b in field])
+    y = np.array([b.y for b in field])
+    dis = np.sqrt(
+        np.subtract.outer(x, x)**2 + np.subtract.outer(y, y)**2)[
+            ~np.eye(len(field), dtype=bool)]
+    assert np.all(
+        [len(field) == N_1 + N_2 - 1,
+         np.allclose(H, [b.H for b in field]),
+         np.allclose(D, [b.D for b in field]),
+         np.allclose(r_b, [b.r_b for b in field]),
+         len(field) == 1 or np.isclose(np.min(dis), min(B_1, B_2)),
+         ])
+
+
+# Test U_shaped_field
+@pytest.mark.parametrize("N_1, N_2, B_1, B_2", [
+        (1, 1, 5., 5.),     # 1 by 1
+        (2, 1, 5., 5.),     # 2 by 1
+        (1, 2, 5., 5.),     # 1 by 2
+        (2, 2, 5., 7.5),    # 2 by 2 (different x/y spacings)
+        (10, 9, 7.5, 5.),   # 10 by 9 (different x/y spacings)
+    ])
+def test_U_shaped_field(N_1, N_2, B_1, B_2):
+    H = 150.        # Borehole length [m]
+    D = 4.          # Borehole buried depth [m]
+    r_b = 0.075     # Borehole radius [m]
+    # Generate the bore field
+    field = gt.boreholes.U_shaped_field(N_1, N_2, B_1, B_2, H, D, r_b)
+    # Evaluate the borehole to borehole distances
+    x = np.array([b.x for b in field])
+    y = np.array([b.y for b in field])
+    dis = np.sqrt(
+        np.subtract.outer(x, x)**2 + np.subtract.outer(y, y)**2)[
+            ~np.eye(len(field), dtype=bool)]
+    assert np.all(
+        [len(field) == N_1 + 2 * N_2 - 2 if N_1 > 1 else N_2,
+         np.allclose(H, [b.H for b in field]),
+         np.allclose(D, [b.D for b in field]),
+         np.allclose(r_b, [b.r_b for b in field]),
+         len(field) == 1 or np.isclose(np.min(dis), min(B_1, B_2)),
+         ])
+
+
+# Test box_shaped_field
+@pytest.mark.parametrize("N_1, N_2, B_1, B_2", [
+        (1, 1, 5., 5.),     # 1 by 1
+        (2, 1, 5., 5.),     # 2 by 1
+        (1, 2, 5., 5.),     # 1 by 2
+        (2, 2, 5., 7.5),    # 2 by 2 (different x/y spacings)
+        (10, 9, 7.5, 5.),   # 10 by 9 (different x/y spacings)
+    ])
+def test_box_shaped_field(N_1, N_2, B_1, B_2):
+    H = 150.        # Borehole length [m]
+    D = 4.          # Borehole buried depth [m]
+    r_b = 0.075     # Borehole radius [m]
+    # Generate the bore field
+    field = gt.boreholes.box_shaped_field(N_1, N_2, B_1, B_2, H, D, r_b)
+    # Evaluate the borehole to borehole distances
+    x = np.array([b.x for b in field])
+    y = np.array([b.y for b in field])
+    dis = np.sqrt(
+        np.subtract.outer(x, x)**2 + np.subtract.outer(y, y)**2)[
+            ~np.eye(len(field), dtype=bool)]
+    if N_1 == 1 and N_2 == 1:
+        nBoreholes_expected = 1
+    elif N_1 == 1:
+        nBoreholes_expected = N_2
+    elif N_2 == 1:
+        nBoreholes_expected = N_1
+    else:
+        nBoreholes_expected = 2 * (N_1 - 1) + 2 * (N_2 - 1)
+    assert np.all(
+        [len(field) == nBoreholes_expected,
+         np.allclose(H, [b.H for b in field]),
+         np.allclose(D, [b.D for b in field]),
+         np.allclose(r_b, [b.r_b for b in field]),
+         len(field) == 1 or np.isclose(np.min(dis), min(B_1, B_2)),
+         ])
+
+
+# Test circle_field
+@pytest.mark.parametrize("N, R", [
+        (1, 5.),    # 1 borehole
+        (2, 5.),    # 2 boreholes
+        (3, 7.5),   # 3 boreholes
+        (10, 9.),   # 10 boreholes
+    ])
+def test_circle_field(N, R):
+    H = 150.        # Borehole length [m]
+    D = 4.          # Borehole buried depth [m]
+    r_b = 0.075     # Borehole radius [m]
+    # Generate the bore field
+    field = gt.boreholes.circle_field(N, R, H, D, r_b)
+    # Evaluate the borehole to borehole distances
+    x = np.array([b.x for b in field])
+    y = np.array([b.y for b in field])
+    dis = np.sqrt(
+        np.subtract.outer(x, x)**2 + np.subtract.outer(y, y)**2)[
+            ~np.eye(len(field), dtype=bool)]
+    B_min = 2 * R * np.sin(np.pi / N)
+    assert np.all(
+        [len(field) == N,
+         np.allclose(H, [b.H for b in field]),
+         np.allclose(D, [b.D for b in field]),
+         np.allclose(r_b, [b.r_b for b in field]),
+         len(field) == 1 or np.isclose(np.min(dis), B_min),
+         len(field) == 1 or np.max(dis) <= (2 + 1e-6) * R,
+         ])
diff --git a/tests/conftest.py b/tests/conftest.py
new file mode 100644
index 00000000..ea33edef
--- /dev/null
+++ b/tests/conftest.py
@@ -0,0 +1,517 @@
+import numpy as np
+import pytest
+
+import pygfunction as gt
+
+
+# =============================================================================
+# boreholes fixtures
+# =============================================================================
+@pytest.fixture
+def single_borehole():
+    H = 150.        # Borehole length [m]
+    D = 4.          # Borehole buried depth [m]
+    r_b = 0.075     # Borehole radius [m]
+    x = 0.          # Borehole x-position [m]
+    y = 0.          # Borehole y-position [m]
+    return [gt.boreholes.Borehole(H, D, r_b, x, y)]
+
+
+@pytest.fixture
+def single_borehole_short():
+    H = 35.         # Borehole length [m]
+    D = 4.          # Borehole buried depth [m]
+    r_b = 0.075     # Borehole radius [m]
+    x = 3.          # Borehole x-position [m]
+    y = 5.          # Borehole y-position [m]
+    return [gt.boreholes.Borehole(H, D, r_b, x, y)]
+
+
+@pytest.fixture
+def ten_boreholes_rectangular():
+    H = 150.            # Borehole length [m]
+    D = 4.              # Borehole buried depth [m]
+    r_b = 0.075         # Borehole radius [m]
+    B_1 = B_2 = 7.5     # Borehole spacing [m]
+    return gt.boreholes.rectangle_field(5, 2, B_1, B_2, H, D, r_b)
+
+
+@pytest.fixture
+def three_boreholes_unequal():
+    return [gt.boreholes.Borehole(150., 4., 0.075, -1., -2.),
+            gt.boreholes.Borehole(88., 2., 0.065, 5., 3.),
+            gt.boreholes.Borehole(177., 5., 0.085, -1., 7.),]
+
+
+# =============================================================================
+# gfunction fixtures
+# =============================================================================
+@pytest.fixture
+def uniform_segments():
+    options = {'nSegments': 12,
+               'segment_ratios': None,
+               'approximate_FLS': False}
+    return options
+
+
+@pytest.fixture
+def uniform_segments_approx():
+    options = {'nSegments': 12,
+               'segment_ratios': None,
+               'approximate_FLS': True,
+               'nFLS': 10}
+    return options
+
+
+@pytest.fixture
+def unequal_segments():
+    segment_ratios = gt.utilities.segment_ratios(8, end_length_ratio=0.02)
+    options = {'nSegments': 8,
+               'segment_ratios': segment_ratios,
+               'approximate_FLS': False}
+    return options
+
+
+@pytest.fixture
+def unequal_segments_approx():
+    segment_ratios = gt.utilities.segment_ratios(8, end_length_ratio=0.02)
+    options = {'nSegments': 8,
+               'segment_ratios': segment_ratios,
+               'approximate_FLS': True,
+               'nFLS': 10}
+    return options
+
+
+# =============================================================================
+# pipes fixtures
+# =============================================================================
+@pytest.fixture
+def single_Utube(single_borehole):
+    # Extract borehole from fixture
+    borehole = single_borehole[0]
+    # Pipe positions [m]
+    D_s = 0.05
+    pos_pipes = [(-D_s, 0.), (D_s, 0.)]
+    k_s = 2.0               # Ground thermal conductivity [W/m.K]
+    k_g = 1.0               # Grout thermal conductivity [W/m.K]
+    k_p = 0.4               # Pipe thermal conductivity [W/m.K]
+    r_out = 0.02            # Pipe outer radius [m]
+    r_in = 0.015            # Pipe inner radius [m]
+    epsilon = 1.0e-06       # Pipe surface roughness [m]
+    m_flow_borehole = 0.05  # Nominal fluid mass flow rate [kg/s]
+    m_flow_pipe = m_flow_borehole
+    # Fluid is propylene-glycol (20 %) at 20 degC
+    fluid = gt.media.Fluid('MPG', 20.)
+    # Pipe thermal resistance [m.K/W]
+    R_p = gt.pipes.conduction_thermal_resistance_circular_pipe(
+        r_in, r_out, k_p)
+    # Convection heat transfer coefficient [W/m2.K]
+    h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
+        m_flow_pipe, r_in, fluid.mu, fluid.rho, fluid.k, fluid.cp,
+        epsilon)
+    # Film thermal resistance [m.K/W]
+    R_f = 1.0/(h_f*2*np.pi*r_in)
+    # Initialize pipe
+    singleUTube = gt.pipes.SingleUTube(
+            pos_pipes, r_in, r_out, borehole, k_s, k_g, R_f + R_p)
+    return singleUTube
+
+
+@pytest.fixture
+def double_Utube_parallel(single_borehole):
+    # Extract borehole from fixture
+    borehole = single_borehole[0]
+    # Pipe positions [m]
+    D_s = 0.05
+    pos_pipes = [(-D_s, 0.), (0., -D_s), (D_s, 0.), (0., D_s)]
+    k_s = 2.0               # Ground thermal conductivity [W/m.K]
+    k_g = 1.0               # Grout thermal conductivity [W/m.K]
+    k_p = 0.4               # Pipe thermal conductivity [W/m.K]
+    r_out = 0.02            # Pipe outer radius [m]
+    r_in = 0.015            # Pipe inner radius [m]
+    epsilon = 1.0e-06       # Pipe surface roughness [m]
+    m_flow_borehole = 0.05  # Nominal fluid mass flow rate [kg/s]
+    m_flow_pipe = m_flow_borehole
+    # Fluid is propylene-glycol (20 %) at 20 degC
+    fluid = gt.media.Fluid('MPG', 20.)
+    # Pipe thermal resistance [m.K/W]
+    R_p = gt.pipes.conduction_thermal_resistance_circular_pipe(
+        r_in, r_out, k_p)
+    # Convection heat transfer coefficient [W/m2.K]
+    h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
+        m_flow_pipe, r_in, fluid.mu, fluid.rho, fluid.k, fluid.cp,
+        epsilon)
+    # Film thermal resistance [m.K/W]
+    R_f = 1.0 / (h_f * 2 * np.pi * r_in)
+    # Initialize pipe
+    doubleUTube = gt.pipes.MultipleUTube(
+            pos_pipes, r_in, r_out, borehole, k_s, k_g, R_f + R_p, 2,
+            config='parallel')
+    return doubleUTube
+
+
+@pytest.fixture
+def double_Utube_series(single_borehole):
+    # Extract borehole from fixture
+    borehole = single_borehole[0]
+    # Pipe positions [m]
+    D_s = 0.05
+    pos_pipes = [(-D_s, 0.), (0., -D_s), (D_s, 0.), (0., D_s)]
+    k_s = 2.0               # Ground thermal conductivity [W/m.K]
+    k_g = 1.0               # Grout thermal conductivity [W/m.K]
+    k_p = 0.4               # Pipe thermal conductivity [W/m.K]
+    r_out = 0.02            # Pipe outer radius [m]
+    r_in = 0.015            # Pipe inner radius [m]
+    epsilon = 1.0e-06       # Pipe surface roughness [m]
+    m_flow_borehole = 0.05  # Nominal fluid mass flow rate [kg/s]
+    m_flow_pipe = m_flow_borehole
+    # Fluid is propylene-glycol (20 %) at 20 degC
+    fluid = gt.media.Fluid('MPG', 20.)
+    # Pipe thermal resistance [m.K/W]
+    R_p = gt.pipes.conduction_thermal_resistance_circular_pipe(
+        r_in, r_out, k_p)
+    # Convection heat transfer coefficient [W/m2.K]
+    h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
+        m_flow_pipe, r_in, fluid.mu, fluid.rho, fluid.k, fluid.cp,
+        epsilon)
+    # Film thermal resistance [m.K/W]
+    R_f = 1.0 / (h_f * 2 * np.pi * r_in)
+    # Initialize pipe
+    doubleUTube = gt.pipes.MultipleUTube(
+            pos_pipes, r_in, r_out, borehole, k_s, k_g, R_f + R_p, 2,
+            config='series')
+    return doubleUTube
+
+
+@pytest.fixture
+def double_Utube_independent(single_borehole):
+    # Extract borehole from fixture
+    borehole = single_borehole[0]
+    # Pipe positions [m]
+    D_s = 0.05
+    pos_pipes = [(-D_s, 0.), (0., -D_s), (D_s, 0.), (0., D_s)]
+    k_s = 2.0               # Ground thermal conductivity [W/m.K]
+    k_g = 1.0               # Grout thermal conductivity [W/m.K]
+    k_p = 0.4               # Pipe thermal conductivity [W/m.K]
+    r_out = 0.02            # Pipe outer radius [m]
+    r_in = 0.015            # Pipe inner radius [m]
+    epsilon = 1.0e-06       # Pipe surface roughness [m]
+    m_flow_borehole = 0.05  # Nominal fluid mass flow rate [kg/s]
+    m_flow_pipe = m_flow_borehole
+    # Fluid is propylene-glycol (20 %) at 20 degC
+    fluid = gt.media.Fluid('MPG', 20.)
+    # Pipe thermal resistance [m.K/W]
+    R_p = gt.pipes.conduction_thermal_resistance_circular_pipe(
+        r_in, r_out, k_p)
+    # Convection heat transfer coefficient [W/m2.K]
+    h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
+        m_flow_pipe, r_in, fluid.mu, fluid.rho, fluid.k, fluid.cp,
+        epsilon)
+    # Film thermal resistance [m.K/W]
+    R_f = 1.0 / (h_f * 2 * np.pi * r_in)
+    # Initialize pipe
+    doubleUTube = gt.pipes.IndependentMultipleUTube(
+            pos_pipes, r_in, r_out, borehole, k_s, k_g, R_f + R_p, 2)
+    return doubleUTube
+
+
+@pytest.fixture
+def triple_Utube_independent(single_borehole):
+    # Extract borehole from fixture
+    borehole = single_borehole[0]
+    # Pipe positions [m]
+    D_s = 0.05
+    pos_pipes = [(-D_s, 0.),
+                 (-D_s*np.cos(np.pi/3), D_s*np.sin(np.pi/3)),
+                 (-D_s*np.cos(np.pi/3), -D_s*np.sin(np.pi/3)),
+                 (D_s, 0.),
+                 (D_s*np.cos(np.pi/3), -D_s*np.sin(np.pi/3)),
+                 (D_s*np.cos(np.pi/3), D_s*np.sin(np.pi/3)),]
+    k_s = 2.0               # Ground thermal conductivity [W/m.K]
+    k_g = 1.0               # Grout thermal conductivity [W/m.K]
+    k_p = 0.4               # Pipe thermal conductivity [W/m.K]
+    r_out = 0.02            # Pipe outer radius [m]
+    r_in = 0.015            # Pipe inner radius [m]
+    epsilon = 1.0e-06       # Pipe surface roughness [m]
+    m_flow_borehole = 0.05  # Nominal fluid mass flow rate [kg/s]
+    m_flow_pipe = m_flow_borehole
+    # Fluid is propylene-glycol (20 %) at 20 degC
+    fluid = gt.media.Fluid('MPG', 20.)
+    # Pipe thermal resistance [m.K/W]
+    R_p = gt.pipes.conduction_thermal_resistance_circular_pipe(
+        r_in, r_out, k_p)
+    # Convection heat transfer coefficient [W/m2.K]
+    h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
+        m_flow_pipe, r_in, fluid.mu, fluid.rho, fluid.k, fluid.cp,
+        epsilon)
+    # Film thermal resistance [m.K/W]
+    R_f = 1.0 / (h_f * 2 * np.pi * r_in)
+    # Initialize pipe
+    tripleUTube = gt.pipes.IndependentMultipleUTube(
+            pos_pipes, r_in, r_out, borehole, k_s, k_g, R_f + R_p, 3)
+    return tripleUTube
+
+
+@pytest.fixture
+def coaxial_annular_in(single_borehole):
+    # Extract borehole from fixture
+    borehole = single_borehole[0]
+    pos = (0., 0.)      # Pipe position [m]
+    k_s = 2.0           # Ground thermal conductivity [W/m.K]
+    k_g = 1.0           # Grout thermal conductivity [W/m.K]
+    k_p = 0.4           # Pipe thermal conductivity [W/m.K]
+    r_in_in = 0.0221    # Inside pipe inner radius [m]
+    r_in_out = 0.025    # Inside pipe outer radius [m]
+    r_out_in = 0.0487   # Outer pipe inside radius [m]
+    r_out_out = 0.055   # Outer pipe outside radius [m]
+    r_inner = np.array([r_out_in, r_in_in])
+    r_outer = np.array([r_out_out, r_in_out])
+    epsilon = 1.0e-06       # Pipe surface roughness [m]
+    m_flow_borehole = 0.05  # Nominal fluid mass flow rate [kg/s]
+    m_flow_pipe = m_flow_borehole
+    # Fluid is propylene-glycol (20 %) at 20 degC
+    fluid = gt.media.Fluid('MPG', 20.)
+    # Pipe thermal resistances [m.K/W]
+    # Inner pipe
+    R_p_in = gt.pipes.conduction_thermal_resistance_circular_pipe(
+        r_in_in, r_in_out, k_p)
+    # Outer pipe
+    R_p_out = gt.pipes.conduction_thermal_resistance_circular_pipe(
+        r_out_in, r_out_out, k_p)
+    # Fluid-to-fluid thermal resistance [m.K/W]
+    # Inner pipe
+    h_f_in = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
+        m_flow_pipe, r_in_in, fluid.mu, fluid.rho, fluid.k, fluid.cp, epsilon)
+    R_f_in = 1.0 / (h_f_in * 2 * np.pi * r_in_in)
+    # Outer pipe
+    h_f_a_in, h_f_a_out = \
+        gt.pipes.convective_heat_transfer_coefficient_concentric_annulus(
+            m_flow_borehole, r_in_out, r_out_in, fluid.mu, fluid.rho, fluid.k,
+            fluid.cp, epsilon)
+    R_f_out_in = 1.0 / (h_f_a_in * 2 * np.pi * r_in_out)
+    R_ff = R_f_in + R_p_in + R_f_out_in
+    # Coaxial GHE in borehole
+    R_f_out_out = 1.0 / (h_f_a_out * 2 * np.pi * r_out_in)
+    R_fp = R_p_out + R_f_out_out
+    # Initialize pipe
+    pipe = gt.pipes.Coaxial(
+        pos, r_inner, r_outer, borehole, k_s, k_g, R_ff, R_fp)
+    return pipe
+
+
+@pytest.fixture
+def coaxial_annular_out(single_borehole):
+    # Extract borehole from fixture
+    borehole = single_borehole[0]
+    pos = (0., 0.)      # Pipe position [m]
+    k_s = 2.0           # Ground thermal conductivity [W/m.K]
+    k_g = 1.0           # Grout thermal conductivity [W/m.K]
+    k_p = 0.4           # Pipe thermal conductivity [W/m.K]
+    r_in_in = 0.0221    # Inside pipe inner radius [m]
+    r_in_out = 0.025    # Inside pipe outer radius [m]
+    r_out_in = 0.0487   # Outer pipe inside radius [m]
+    r_out_out = 0.055   # Outer pipe outside radius [m]
+    r_inner = np.array([r_in_in, r_out_in])
+    r_outer = np.array([r_in_out, r_out_out])
+    epsilon = 1.0e-06       # Pipe surface roughness [m]
+    m_flow_borehole = 0.05  # Nominal fluid mass flow rate [kg/s]
+    m_flow_pipe = m_flow_borehole
+    # Fluid is propylene-glycol (20 %) at 20 degC
+    fluid = gt.media.Fluid('MPG', 20.)
+    # Pipe thermal resistances [m.K/W]
+    # Inner pipe
+    R_p_in = gt.pipes.conduction_thermal_resistance_circular_pipe(
+        r_in_in, r_in_out, k_p)
+    # Outer pipe
+    R_p_out = gt.pipes.conduction_thermal_resistance_circular_pipe(
+        r_out_in, r_out_out, k_p)
+    # Fluid-to-fluid thermal resistance [m.K/W]
+    # Inner pipe
+    h_f_in = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
+        m_flow_pipe, r_in_in, fluid.mu, fluid.rho, fluid.k, fluid.cp, epsilon)
+    R_f_in = 1.0 / (h_f_in * 2 * np.pi * r_in_in)
+    # Outer pipe
+    h_f_a_in, h_f_a_out = \
+        gt.pipes.convective_heat_transfer_coefficient_concentric_annulus(
+            m_flow_borehole, r_in_out, r_out_in, fluid.mu, fluid.rho, fluid.k,
+            fluid.cp, epsilon)
+    R_f_out_in = 1.0 / (h_f_a_in * 2 * np.pi * r_in_out)
+    R_ff = R_f_in + R_p_in + R_f_out_in
+    # Coaxial GHE in borehole
+    R_f_out_out = 1.0 / (h_f_a_out * 2 * np.pi * r_out_in)
+    R_fp = R_p_out + R_f_out_out
+    # Initialize pipe
+    pipe = gt.pipes.Coaxial(
+        pos, r_inner, r_outer, borehole, k_s, k_g, R_ff, R_fp)
+    return pipe
+
+
+# =============================================================================
+# networks fixtures
+# =============================================================================
+@pytest.fixture
+def single_borehole_network(single_borehole):
+    # Extract borehole from fixture
+    boreField = single_borehole
+    # Pipe positions [m]
+    D_s = 0.05
+    pos_pipes = [(-D_s, 0.), (D_s, 0.)]
+    k_s = 2.0               # Ground thermal conductivity [W/m.K]
+    k_g = 1.0               # Grout thermal conductivity [W/m.K]
+    k_p = 0.4               # Pipe thermal conductivity [W/m.K]
+    r_out = 0.02            # Pipe outer radius [m]
+    r_in = 0.015            # Pipe inner radius [m]
+    epsilon = 1.0e-06       # Pipe surface roughness [m]
+    m_flow_borehole = 0.05  # Nominal fluid mass flow rate [kg/s]
+    m_flow_pipe = m_flow_borehole
+    m_flow_network = m_flow_borehole
+    # Fluid is propylene-glycol (20 %) at 20 degC
+    fluid = gt.media.Fluid('MPG', 20.)
+    # Pipe thermal resistance [m.K/W]
+    R_p = gt.pipes.conduction_thermal_resistance_circular_pipe(
+        r_in, r_out, k_p)
+    # Convection heat transfer coefficient [W/m2.K]
+    h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
+        m_flow_pipe, r_in, fluid.mu, fluid.rho, fluid.k, fluid.cp,
+        epsilon)
+    # Film thermal resistance [m.K/W]
+    R_f = 1.0 / (h_f * 2 * np.pi * r_in)
+
+    # Build network of parallel-connected boreholes
+    bore_connectivity = [-1] * len(boreField)
+    UTubes = []
+    for borehole in boreField:
+        UTube = gt.pipes.SingleUTube(
+            pos_pipes, r_in, r_out, borehole, k_s, k_g, R_f + R_p)
+        UTubes.append(UTube)
+    # Initialize network
+    network = gt.networks.Network(
+        boreField, UTubes, bore_connectivity=bore_connectivity,
+        m_flow_network=m_flow_network, cp_f=fluid.cp)
+    return network
+
+
+@pytest.fixture
+def single_borehole_network_short(single_borehole_short):
+    # Extract borehole from fixture
+    boreField = single_borehole_short
+    # Pipe positions [m]
+    D_s = 0.05
+    pos_pipes = [(-D_s, 0.), (D_s, 0.)]
+    k_s = 2.0               # Ground thermal conductivity [W/m.K]
+    k_g = 1.0               # Grout thermal conductivity [W/m.K]
+    k_p = 0.4               # Pipe thermal conductivity [W/m.K]
+    r_out = 0.02            # Pipe outer radius [m]
+    r_in = 0.015            # Pipe inner radius [m]
+    epsilon = 1.0e-06       # Pipe surface roughness [m]
+    m_flow_borehole = 0.05  # Nominal fluid mass flow rate [kg/s]
+    m_flow_pipe = m_flow_borehole
+    m_flow_network = m_flow_borehole
+    # Fluid is propylene-glycol (20 %) at 20 degC
+    fluid = gt.media.Fluid('MPG', 20.)
+    # Pipe thermal resistance [m.K/W]
+    R_p = gt.pipes.conduction_thermal_resistance_circular_pipe(
+        r_in, r_out, k_p)
+    # Convection heat transfer coefficient [W/m2.K]
+    h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
+        m_flow_pipe, r_in, fluid.mu, fluid.rho, fluid.k, fluid.cp,
+        epsilon)
+    # Film thermal resistance [m.K/W]
+    R_f = 1.0 / (h_f * 2 * np.pi * r_in)
+
+    # Build network of parallel-connected boreholes
+    bore_connectivity = [-1] * len(boreField)
+    UTubes = []
+    for borehole in boreField:
+        UTube = gt.pipes.SingleUTube(
+            pos_pipes, r_in, r_out, borehole, k_s, k_g, R_f + R_p)
+        UTubes.append(UTube)
+    # Initialize network
+    network = gt.networks.Network(
+        boreField, UTubes, bore_connectivity=bore_connectivity,
+        m_flow_network=m_flow_network, cp_f=fluid.cp)
+    return network
+
+
+@pytest.fixture
+def ten_boreholes_network_rectangular(ten_boreholes_rectangular):
+    # Extract bore field from fixture
+    boreField = ten_boreholes_rectangular
+    # Pipe positions [m]
+    D_s = 0.05
+    pos_pipes = [(-D_s, 0.), (D_s, 0.)]
+    k_s = 2.0               # Ground thermal conductivity [W/m.K]
+    k_g = 1.0               # Grout thermal conductivity [W/m.K]
+    k_p = 0.4               # Pipe thermal conductivity [W/m.K]
+    r_out = 0.02            # Pipe outer radius [m]
+    r_in = 0.015            # Pipe inner radius [m]
+    epsilon = 1.0e-06       # Pipe surface roughness [m]
+    m_flow_borehole = 0.05  # Nominal fluid mass flow rate [kg/s]
+    m_flow_pipe = m_flow_borehole
+    m_flow_network = m_flow_borehole * 5
+    # Fluid is propylene-glycol (20 %) at 20 degC
+    fluid = gt.media.Fluid('MPG', 20.)
+    # Pipe thermal resistance [m.K/W]
+    R_p = gt.pipes.conduction_thermal_resistance_circular_pipe(
+        r_in, r_out, k_p)
+    # Convection heat transfer coefficient [W/m2.K]
+    h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
+        m_flow_pipe, r_in, fluid.mu, fluid.rho, fluid.k, fluid.cp,
+        epsilon)
+    # Film thermal resistance [m.K/W]
+    R_f = 1.0 / (h_f * 2 * np.pi * r_in)
+
+    # Build network of parallel-connected boreholes
+    bore_connectivity = [-1] * len(boreField)
+    UTubes = []
+    for borehole in boreField:
+        UTube = gt.pipes.SingleUTube(
+            pos_pipes, r_in, r_out, borehole, k_s, k_g, R_f + R_p)
+        UTubes.append(UTube)
+    # Initialize network
+    network = gt.networks.Network(
+        boreField, UTubes, bore_connectivity=bore_connectivity,
+        m_flow_network=m_flow_network, cp_f=fluid.cp)
+    return network
+
+
+@pytest.fixture
+def ten_boreholes_network_rectangular_series(ten_boreholes_rectangular):
+    # Extract bore field from fixture
+    boreField = ten_boreholes_rectangular
+    # Pipe positions [m]
+    D_s = 0.05
+    pos_pipes = [(-D_s, 0.), (D_s, 0.)]
+    k_s = 2.0               # Ground thermal conductivity [W/m.K]
+    k_g = 1.0               # Grout thermal conductivity [W/m.K]
+    k_p = 0.4               # Pipe thermal conductivity [W/m.K]
+    r_out = 0.02            # Pipe outer radius [m]
+    r_in = 0.015            # Pipe inner radius [m]
+    epsilon = 1.0e-06       # Pipe surface roughness [m]
+    m_flow_borehole = 0.05  # Nominal fluid mass flow rate [kg/s]
+    m_flow_pipe = m_flow_borehole
+    m_flow_network = m_flow_borehole
+    # Fluid is propylene-glycol (20 %) at 20 degC
+    fluid = gt.media.Fluid('MPG', 20.)
+    # Pipe thermal resistance [m.K/W]
+    R_p = gt.pipes.conduction_thermal_resistance_circular_pipe(
+        r_in, r_out, k_p)
+    # Convection heat transfer coefficient [W/m2.K]
+    h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
+        m_flow_pipe, r_in, fluid.mu, fluid.rho, fluid.k, fluid.cp,
+        epsilon)
+    # Film thermal resistance [m.K/W]
+    R_f = 1.0 / (h_f * 2 * np.pi * r_in)
+
+    # Build network of series-connected boreholes
+    bore_connectivity = list(range(-1, len(boreField)-1))
+    UTubes = []
+    for borehole in boreField:
+        UTube = gt.pipes.SingleUTube(
+            pos_pipes, r_in, r_out, borehole, k_s, k_g, R_f + R_p)
+        UTubes.append(UTube)
+    # Initialize network
+    network = gt.networks.Network(
+        boreField, UTubes, bore_connectivity=bore_connectivity,
+        m_flow_network=m_flow_network, cp_f=fluid.cp)
+    return network
diff --git a/tests/gfunction_test.py b/tests/gfunction_test.py
index 10b593c4..a525b569 100644
--- a/tests/gfunction_test.py
+++ b/tests/gfunction_test.py
@@ -1,621 +1,229 @@
 # -*- coding: utf-8 -*-
 """ Test suite for gfunction module.
 """
-import unittest
-
 import numpy as np
-
-import pygfunction.utilities
-
-
-class TestSolvers(unittest.TestCase):
-    """
-    Test cases for calculation of g-functions using various solver methods
-    available.
-    """
-
-    def setUp(self) -> None:
-        self.H = 150.  # Borehole length [m]
-        self.D = 4.  # Borehole buried depth [m]
-        self.r_b = 0.075  # Borehole radius [m]
-        self.B = 7.5  # Borehole spacing [m]
-        self.alpha = 1.0e-6  # Ground thermal diffusivity [m2/s]
-
-    def test_detailed(self):
-        from pygfunction.gfunction import gFunction
-        from pygfunction.boreholes import rectangle_field
-        from pygfunction.utilities import time_geometric, segment_ratios
-        # Calculation of the g-function at the same time values
-        N_1 = 3
-        N_2 = 2
-        boreField = rectangle_field(N_1, N_2, self.B, self.B,
-                                    self.H, self.D, self.r_b)
-
-        # Thermal properties
-        alpha = 1.0e-6  # Ground thermal diffusivity (m2/s)
-
-        # Geometrically expanding time vector.
-        dt = 100 * 3600.  # Time step
-        tmax = 3000. * 8760. * 3600.  # Maximum time
-        Nt = 10  # Number of time steps
-        time = time_geometric(dt, tmax, Nt)
-
-        # g-Function calculation options
-        options_uniform = {'nSegments': 24,
-                          'segment_ratios': None,
-                          'disp': False}
-        ratios = segment_ratios(8)
-        options_unequal = {'nSegments': 8,
-                          'segment_ratios': ratios,
-                          'disp': False}
-        method = 'detailed'
-
-        gfunc_UBWT_uniform = gFunction(
-            boreField, alpha, time=time, options=options_uniform,
-            method=method).gFunc
-
-        gfunc_UBWT_unequal = gFunction(
-            boreField, alpha, time=time, options=options_unequal,
-            method=method).gFunc
-
-        self.assertTrue(np.allclose(gfunc_UBWT_uniform, gfunc_UBWT_unequal,
-                                    rtol=1e-2, atol=1e-6),
-                        msg='Incorrect g-function for the detailed solver test.'
-                            'A converged solution is compared to a discretized'
-                            'solution using the UBWT boundary condition.')
-
-    def test_similarities(self):
-        from pygfunction.gfunction import gFunction
-        from pygfunction.boreholes import rectangle_field
-        from pygfunction.utilities import time_geometric, segment_ratios
-        # Calculation of the g-function at the same time values
-        N_1 = 3
-        N_2 = 2
-        boreField = rectangle_field(N_1, N_2, self.B, self.B,
-                                    self.H, self.D, self.r_b)
-
-        # Thermal properties
-        alpha = 1.0e-6  # Ground thermal diffusivity (m2/s)
-
-        # Geometrically expanding time vector.
-        dt = 100 * 3600.  # Time step
-        tmax = 3000. * 8760. * 3600.  # Maximum time
-        Nt = 10  # Number of time steps
-        time = time_geometric(dt, tmax, Nt)
-
-        # g-Function calculation options
-        options_uniform = {'nSegments': 24,
-                          'segment_ratios': None,
-                          'disp': False}
-        ratios = segment_ratios(8)
-        options_unequal = {'nSegments': 8,
-                          'segment_ratios': ratios,
-                          'disp': False}
-        method = 'similarities'
-
-        gfunc_UBWT_uniform = gFunction(
-            boreField, alpha, time=time, options=options_uniform,
-            method=method).gFunc
-
-        gfunc_UBWT_unequal = gFunction(
-            boreField, alpha, time=time, options=options_unequal,
-            method=method).gFunc
-
-        self.assertTrue(np.allclose(gfunc_UBWT_uniform, gfunc_UBWT_unequal,
-                                    rtol=1e-2, atol=1e-6),
-                        msg='Incorrect g-function for the similarities solver '
-                            'test.'
-                            'A converged solution is compared to a discretized'
-                            'solution using the UT boundary condition.')
-
-    def test_equivalent(self):
-        from pygfunction.gfunction import gFunction
-        from pygfunction.boreholes import rectangle_field
-        from pygfunction.utilities import time_geometric, segment_ratios
-        # Calculation of the g-function at the same time values
-        N_1 = 3
-        N_2 = 2
-        boreField = rectangle_field(N_1, N_2, self.B, self.B,
-                                    self.H, self.D, self.r_b)
-
-        # Thermal properties
-        alpha = 1.0e-6  # Ground thermal diffusivity (m2/s)
-
-        # Geometrically expanding time vector.
-        dt = 100 * 3600.  # Time step
-        tmax = 3000. * 8760. * 3600.  # Maximum time
-        Nt = 10  # Number of time steps
-        time = time_geometric(dt, tmax, Nt)
-
-        # g-Function calculation options
-        options_uniform = {'nSegments': 24,
-                          'segment_ratios': None,
-                          'disp': False}
-        ratios = segment_ratios(8)
-        options_unequal = {'nSegments': 8,
-                          'segment_ratios': ratios,
-                          'disp': False}
-        method = 'equivalent'
-
-        gfunc_UBWT_uniform = gFunction(
-            boreField, alpha, time=time, options=options_uniform,
-            method=method).gFunc
-
-        gfunc_UBWT_unequal = gFunction(
-            boreField, alpha, time=time, options=options_unequal,
-            method=method).gFunc
-
-        self.assertTrue(np.allclose(gfunc_UBWT_uniform, gfunc_UBWT_unequal,
-                                    rtol=1e-2, atol=1e-6),
-                        msg='Incorrect g-function for the equivalent solver '
-                            'test.'
-                            'A converged solution is compared to a discretized'
-                            'solution using the UT boundary condition.')
-
-
-class TestUniformHeatExtractionRate(unittest.TestCase):
-    """ Test cases for calculation of g-functions using uniform heat extraction
-        rate boundary condition.
-    """
-
-    def setUp(self):
-        self.H = 150.           # Borehole length [m]
-        self.D = 4.             # Borehole buried depth [m]
-        self.r_b = 0.075        # Borehole radius [m]
-        self.B = 7.5            # Borehole spacing [m]
-        self.alpha = 1.0e-6     # Ground thermal diffusivity [m2/s]
-
-    def test_one_borehole(self, rel_tol=1.0e-4):
-        """ Tests the value of the g-function of one borehole.
-        """
-        from pygfunction.gfunction import uniform_heat_extraction
-        from pygfunction.boreholes import rectangle_field
-
-        # Results of Cimmino and Bernier (2014)
-        time = np.array([1070, 33554478])*3600.
-        g_ref = np.array([3.65501492640735,
-                          6.68675948788837])
-
-        # Calculation of the g-function at the same time values
-        N_1 = 1
-        N_2 = 1
-        boreField = rectangle_field(N_1, N_2, self.B, self.B,
-                                    self.H, self.D, self.r_b)
-        g = uniform_heat_extraction(boreField, time, self.alpha)
-
-        self.assertTrue(np.allclose(g, g_ref, rtol=rel_tol, atol=1e-6),
-                        msg='Incorrect values of the g-function of one '
-                            'borehole for uniform heat extraction rate.')
-
-    def test_three_by_two(self, rel_tol=1.0e-4):
-        """ Tests the value of the g-function of a 3 by 2 bore field.
-        """
-        from pygfunction.gfunction import uniform_heat_extraction
-        from pygfunction.boreholes import rectangle_field
-
-        # Results of Cimmino and Bernier (2014)
-        time = np.array([1070, 33554478])*3600.
-        g_ref = np.array([3.66173047944222,
-                          16.0243474955568])
-
-        # Calculation of the g-function at the same time values
-        N_1 = 3
-        N_2 = 2
-        boreField = rectangle_field(N_1, N_2, self.B, self.B,
-                                    self.H, self.D, self.r_b)
-        g = uniform_heat_extraction(boreField, time, self.alpha)
-
-        self.assertTrue(np.allclose(g, g_ref, rtol=rel_tol, atol=1e-6),
-                        msg='Incorrect values of the g-function of three by '
-                            'two field for uniform heat extraction rate.')
-
-
-class TestUniformTemperature(unittest.TestCase):
-    """ Test cases for calculation of g-functions using uniform borehole wall
-        temperature boundary condition.
-    """
-
-    def setUp(self):
-        self.H = 150.           # Borehole length [m]
-        self.D = 4.             # Borehole buried depth [m]
-        self.r_b = 0.075        # Borehole radius [m]
-        self.B = 7.5            # Borehole spacing [m]
-        self.alpha = 1.0e-6     # Ground thermal diffusivity [m2/s]
-
-    def test_one_borehole_one_segment(self, rel_tol=1.0e-4):
-        """ Tests the value of the g-function of one borehole.
-        """
-        from pygfunction.gfunction import uniform_temperature
-        from pygfunction.boreholes import rectangle_field
-
-        # Results of Cimmino and Bernier (2014)
-        time = np.array([1070, 33554478])*3600.
-        g_ref = np.array([3.65502692098609,
-                          6.68675948788828])
-
-        # Calculation of the g-function at the same time values
-        N_1 = 1
-        N_2 = 1
-        boreField = rectangle_field(N_1, N_2, self.B, self.B,
-                                    self.H, self.D, self.r_b)
-        g = uniform_temperature(boreField, time, self.alpha, nSegments=1)
-
-        self.assertTrue(np.allclose(g, g_ref, rtol=rel_tol, atol=1e-6),
-                        msg='Incorrect values of the g-function of one '
-                            'borehole for uniform temperature (1 segment).')
-
-    def test_three_by_two_one_segment(self, rel_tol=1.0e-4):
-        """ Tests the value of the g-function of a 3 by 2 bore field.
-        """
-        from pygfunction.gfunction import uniform_temperature
-        from pygfunction.boreholes import rectangle_field
-
-        # Results of Cimmino and Bernier (2014)
-        time = np.array([1070, 33554478])*3600.
-        g_ref = np.array([3.6617432932388,
-                          15.9710298219412])
-
-        # Calculation of the g-function at the same time values
-        N_1 = 3
-        N_2 = 2
-        boreField = rectangle_field(N_1, N_2, self.B, self.B,
-                                    self.H, self.D, self.r_b)
-        g = uniform_temperature(boreField, time, self.alpha, nSegments=1)
-
-        self.assertTrue(np.allclose(g, g_ref, rtol=rel_tol, atol=1e-6),
-                        msg='Incorrect values of the g-function of three by '
-                            'two field for uniform temperature (1 segment).')
-
-    def test_one_borehole_twelve_segments(self, rel_tol=1.0e-4):
-        """ Tests the value of the g-function of one borehole.
-        """
-        from pygfunction.gfunction import uniform_temperature
-        from pygfunction.boreholes import rectangle_field
-
-        # Results of Cimmino and Bernier (2014)
-        time = np.array([1070, 33554478])*3600.
-        g_ref = np.array([3.65476902065229,
-                          6.6329923757352])
-
-        # Calculation of the g-function at the same time values
-        N_1 = 1
-        N_2 = 1
-        boreField = rectangle_field(N_1, N_2, self.B, self.B,
-                                    self.H, self.D, self.r_b)
-        g = uniform_temperature(
-            boreField, time, self.alpha, nSegments=12, segment_ratios=None)
-
-        self.assertTrue(np.allclose(g, g_ref, rtol=rel_tol, atol=1e-6),
-                        msg='Incorrect values of the g-function of one '
-                            'borehole for uniform temperature (12 segments).')
-
-    def test_three_by_two_twelve_segments(self, rel_tol=1.0e-4):
-        """ Tests the value of the g-function of a 3 by 2 bore field.
-        """
-        from pygfunction.gfunction import uniform_temperature
-        from pygfunction.boreholes import rectangle_field
-
-        # Results of Cimmino and Bernier (2014)
-        time = 33554478*3600.
-        g_ref = 15.1697321426028
-
-        # Calculation of the g-function at the same time values
-        N_1 = 3
-        N_2 = 2
-        boreField = rectangle_field(N_1, N_2, self.B, self.B,
-                                    self.H, self.D, self.r_b)
-        g = uniform_temperature(
-            boreField, time, self.alpha, nSegments=12, segment_ratios=None)
-
-        self.assertTrue(np.allclose(g, g_ref, rtol=rel_tol, atol=1e-6),
-                        msg='Incorrect values of the g-function of three by '
-                            'two field for uniform temperature (12 segments).')
-
-    def test_unequal_segments(self, rel_tol=1.0e-3):
-        from pygfunction.gfunction import gFunction
-        from pygfunction.boreholes import rectangle_field
-        from pygfunction.utilities import time_geometric
-        # Calculation of the g-function at the same time values
-        N_1 = 3
-        N_2 = 2
-        boreField = rectangle_field(N_1, N_2, self.B, self.B,
-                                    self.H, self.D, self.r_b)
-
-        # Geometrically expanding time vector.
-        dt = 100 * 3600.  # Time step
-        tmax = 3000. * 8760. * 3600.  # Maximum time
-        Nt = 10  # Number of time steps
-        time = time_geometric(dt, tmax, Nt)
-
-        nSegments = [12, 11, 13, 12, 11, 13]
-
-        # g-Function calculation option for uniform borehole segment lengths
-        # in the field by defining nSegments as an integer >= 1
-        options = {'nSegments': nSegments[0],
-                   'segment_ratios': None}
-        method = 'similarities'
-        gfunc = gFunction(
-            boreField, self.alpha, time=time, options=options, method=method)
-        g_ref = gfunc.gFunc
-
-        # g-Function calculation with nSegments passed as list, where the list
-        # is of the same length as boreField and each borehole is defined to
-        # have >= 1 segment
-        # Note: nSegments[i] pertains to boreField[i]
-        options = {'nSegments': nSegments,
-                   'segment_ratios': None}
-        method = 'similarities'
-        gfunc = gFunction(
-            boreField, self.alpha, time=time, options=options, method=method)
-        g = gfunc.gFunc
-
-        self.assertTrue(np.allclose(g, g_ref, rtol=rel_tol, atol=1e-6),
-                        msg='Incorrect values of the g-function of six '
-                            'boreholes for uniform temperature and '
-                            'unequal numbers of segments.')
-
-
-class TestMixedInletTemperature(unittest.TestCase):
-    """ Test cases for calculation of g-functions using mixed inlet fluid
-        temperature boundary condition.
-    """
-
-    def setUp(self):
-        self.H = 150.           # Borehole length [m]
-        self.D = 4.             # Borehole buried depth [m]
-        self.r_b = 0.075        # Borehole radius [m]
-        self.B = 7.5            # Borehole spacing [m]
-        self.alpha = 1.0e-6     # Ground thermal diffusivity [m2/s]
-        self.k_s = 2.0          # Ground thermal conductivity [W/m.K]
-        self.k_g = 1.0          # Grout thermal conductivity [W/m.K]
-        self.k_p = 0.4          # Pipe thermal conductivity [W/m.K]
-        self.r_out = 0.02       # Pipe outer radius [m]
-        self.r_in = 0.015       # Pipe inner radius [m]
-        self.D_s = 0.05         # Shank spacing [m]
-        self.epsilon = 1.0e-06  # Pipe roughness [m]
-        self.m_flow_network = 0.25  # Fluid mass flor rate in network [kg/s]
-
-    def test_unequal_segments(self, rel_tol=1.0e-2):
-        """ Tests the value of the g-function of six boreholes in series with
-            unequal numbers of segments.
-        """
-        from pygfunction.gfunction import gFunction
-        from pygfunction.boreholes import rectangle_field
-        from pygfunction.utilities import time_geometric, segment_ratios
-        from pygfunction.pipes import \
-            conduction_thermal_resistance_circular_pipe, \
-            convective_heat_transfer_coefficient_circular_pipe, \
-            SingleUTube
-        from pygfunction.media import Fluid
-        from pygfunction.networks import Network
-        from numpy import pi
-        # Calculation of the g-function at the same time values
-        N_1 = 3
-        N_2 = 2
-        boreField = rectangle_field(N_1, N_2, self.B, self.B,
-                                    self.H, self.D, self.r_b)
-
-        # Geometrically expanding time vector.
-        dt = 100 * 3600.  # Time step
-        tmax = 3000. * 8760. * 3600.  # Maximum time
-        Nt = 10  # Number of time steps
-        time = time_geometric(dt, tmax, Nt)
-
-        nSegments = [12, 11, 13, 12, 11, 13]
-
-        # Fluid is propylene-glycol (20 %) at 20 degC
-        fluid = Fluid('MPG', 20.)
-
-        bore_connectivity = list(range(-1, 5))
-        # Pipe thermal resistance
-        R_p = conduction_thermal_resistance_circular_pipe(
-            self.r_in, self.r_out, self.k_p)
-        m_flow_pipe = self.m_flow_network  # all boreholes in series
-        h_f = convective_heat_transfer_coefficient_circular_pipe(
-            m_flow_pipe, self.r_in, fluid.mu, fluid.rho, fluid.k, fluid.cp,
-            self.epsilon)
-        R_f = 1.0/(h_f*2*pi*self.r_in)
-
-        pos_pipes = [(-self.D_s, 0), (self.D_s, 0)]
-
-        UTubes = []
-        for borebole in boreField:
-            UTube = SingleUTube(pos_pipes, self.r_in, self.r_out, borebole,
-                                self.k_s, self.k_g, R_f + R_p)
-            UTubes.append(UTube)
-        network = Network(
-            boreField, UTubes, bore_connectivity=bore_connectivity,
-            m_flow_network=m_flow_pipe, cp_f=fluid.cp)
-
-        # g-Function calculation option for uniform borehole segment lengths
-        # in the field by defining nSegments as an integer >= 1
-        options = {'nSegments': nSegments[0],
-                   'segment_ratios': None}
-        method = 'similarities'
-        gfunc = gFunction(
-            network, self.alpha, time=time, options=options, method=method)
-        g_reference = gfunc.gFunc
-
-        # g-Function calculation with nSegments passed as list, where the list
-        # is of the same length as boreField and each borehole is defined to
-        # have >= 1 segment
-        # Note: nSegments[i] pertains to boreField[i]
-        options = {'nSegments': nSegments}
-        method = 'similarities'
-        gfunc = gFunction(
-            network, self.alpha, time=time, options=options, method=method)
-        g_unequal_nSegments = gfunc.gFunc
-
-        self.assertTrue(
-            np.allclose(
-                g_unequal_nSegments, g_reference, rtol=rel_tol, atol=1e-6),
-            msg='Incorrect values of the g-function of six '
-                'boreholes for mixed inlet temperature and '
-                'unequal numbers of segments.')
-
-        # Compute g-function with predefined segment lengths
-        nSegments = 8
-        # Define the segment ratios for each borehole in each segment
-        # the segment lengths are defined top to bottom left to right
-        ratios = np.array(
-            [0.05, 0.10, 0.10, 0.25, 0.25, 0.10, 0.10, 0.05])
-        options = {'nSegments': nSegments,
-                   'segment_ratios': ratios,
-                   'disp': False}
-        method = 'similarities'
-
-        gfunc = gFunction(
-            network, self.alpha, time=time, options=options, method=method)
-        g_predefined_ratios = gfunc.gFunc
-
-        self.assertTrue(
-            np.allclose(
-                g_predefined_ratios, g_reference, rtol=rel_tol, atol=1e-6),
-            msg='Incorrect values of the g-function of six '
-                'boreholes for mixed inlet temperature and '
-                'unequal numbers of segments.')
-
-        # Compute g-function with discretized segment lengths
-        ratios = segment_ratios(8)
-        options = {'nSegments': nSegments,
-                   'segment_ratios': ratios,
-                   'disp': False}
-        method = 'similarities'
-
-        gfunc = gFunction(
-            network, self.alpha, time=time, options=options, method=method)
-        g_generated_ratios = gfunc.gFunc
-
-        self.assertTrue(
-            np.allclose(
-                g_generated_ratios, g_reference, rtol=rel_tol, atol=1e-6),
-            msg='Incorrect values of the g-function of six '
-                'boreholes for mixed inlet temperature and '
-                'unequal numbers of segments.')
-
-
-class TestEqualInletTemperature(unittest.TestCase):
-    """ Test cases for calculation of g-functions using equal inlet fluid
-        temperature boundary condition.
-    """
-
-    def setUp(self):
-        self.H = 150.           # Borehole length [m]
-        self.D = 4.             # Borehole buried depth [m]
-        self.r_b = 0.075        # Borehole radius [m]
-        self.B = 7.5            # Borehole spacing [m]
-        self.alpha = 1.0e-6     # Ground thermal diffusivity [m2/s]
-        self.k_s = 2.0          # Ground thermal conductivity [W/m.K]
-        self.m_flow = 0.5       # Fluid mass flow rate (per borehole) [kg/s]
-        self.cp = 4000.         # Fluid specific heat capacity [J/kg.K]
-
-    def test_one_borehole_twenty_four_segments(self, rel_tol=1.0e-4):
-        """ Tests the value of the g-function of one borehole.
-        """
-        from pygfunction.pipes import SingleUTube
-        from pygfunction.gfunction import equal_inlet_temperature
-        from pygfunction.boreholes import rectangle_field
-
-        # Pipe characteristics (will be overwritten by thermal resistances)
-        r_out = 0.01
-        r_in = 0.0075
-        Ds = 0.06
-        pos = [(-Ds, 0.), (Ds, 0.)]
-        k_g = 1.5
-
-        # Reference results
-        # Number of line source segments per borehole
-        nSegments = 24
-        # Very large time values to approximate steady-state
-        time = 1.0e20
-        # Borehole thermal resistances [m.K/W]
-        Rb = np.array([0.1, 0.3])
-        # Borehole short-circuit resistance [m.K/W]
-        Ra = 0.25
-        g_ref = np.array([6.63668744, 6.65139651])
-
-        # Calculation of the g-function at the same time values
-        N_1 = 1
-        N_2 = 1
-        boreField = rectangle_field(N_1, N_2, self.B, self.B,
-                                    self.H, self.D, self.r_b)
-        g = np.zeros_like(g_ref)
-        for i in range(len(Rb)):
-            UTubes = [SingleUTube(
-                    pos, r_in, r_out, borehole, self.k_s, k_g, 1e-30, J=0)
-                      for borehole in boreField]
-            # Overwrite borehole resistances
-            Rd_00 = 2*Rb[i]
-            Rd_01 = 2*Ra*Rd_00/(2*Rd_00 - Ra)
-            for UTube in UTubes:
-                UTube._Rd[0,0] = Rd_00
-                UTube._Rd[1,1] = Rd_00
-                UTube._Rd[0,1] = Rd_01
-                UTube._Rd[1,0] = Rd_01
-                UTube._update_model_variables(self.m_flow, self.cp, nSegments, None)
-            g[i] = equal_inlet_temperature(
-                boreField, UTubes, self.m_flow, self.cp, time, self.alpha,
-                nSegments=nSegments, segment_ratios=None)
-
-        self.assertTrue(np.allclose(g, g_ref, rtol=rel_tol, atol=1e-6),
-                        msg='Incorrect values of the g-function of one '
-                            'borehole for equal inlet fluid temperature '
-                            '(24 segments).')
-
-    def test_four_by_four_twenty_four_segments(self, rel_tol=1.0e-4):
-        """ Tests the value of the g-function of 4 by 4 bore field.
-        """
-        from pygfunction.pipes import SingleUTube
-        from pygfunction.gfunction import equal_inlet_temperature
-        from pygfunction.boreholes import rectangle_field
-
-        # Pipe characteristics (will be overwritten by thermal resistances)
-        r_out = 0.01
-        r_in = 0.0075
-        Ds = 0.06
-        pos = [(-Ds, 0.), (Ds, 0.)]
-        k_g = 1.5
-
-        # Reference results
-        # Number of line source segments per borehole
-        nSegments = 24
-        # Very large time values to approximate steady-state
-        time = 1.0e20
-        # Borehole thermal resistances [m.K/W]
-        Rb = 1.0
-        # Borehole short-circuit resistance [m.K/W]
-        Ra = 0.25
-        g_ref = 28.01257172
-
-        # Calculation of the g-function at the same time values
-        N_1 = 4
-        N_2 = 4
-        boreField = rectangle_field(N_1, N_2, self.B, self.B,
-                                    self.H, self.D, self.r_b)
-        g = np.zeros_like(g_ref)
-        UTubes = [SingleUTube(
-                pos, r_in, r_out, borehole, self.k_s, k_g, 1e-30, J=0)
-                  for borehole in boreField]
-        # Overwrite borehole resistances
-        Rd_00 = 2*Rb
-        Rd_01 = 2*Ra*Rd_00/(2*Rd_00 - Ra)
-        for UTube in UTubes:
-            UTube._Rd[0,0] = Rd_00
-            UTube._Rd[1,1] = Rd_00
-            UTube._Rd[0,1] = Rd_01
-            UTube._Rd[1,0] = Rd_01
-            UTube._update_model_variables(self.m_flow, self.cp, nSegments, None)
-        g = equal_inlet_temperature(
-            boreField, UTubes, self.m_flow, self.cp, time, self.alpha,
-            nSegments=nSegments, segment_ratios=None)
-
-        self.assertTrue(np.allclose(g, g_ref, rtol=rel_tol, atol=1e-6),
-                        msg='Incorrect values of the g-function of four by '
-                            'four field for equal inlet fluid temperature '
-                            '(24 segments).')
-        
-        
-if __name__ == '__main__' and __package__ is None:
-    from os import sys, path
-    sys.path.append(path.dirname(path.dirname(path.abspath(__file__))))
-    unittest.main()
+import pytest
+
+import pygfunction as gt
+
+
+# Test 'UBWT' g-functions for different bore fields using all solvers,
+# unequal/uniform segments, and with/without the FLS approximation
+@pytest.mark.parametrize("field, method, opts, expected", [
+        #  'equivalent' solver - unequal segments
+        ('single_borehole', 'equivalent', 'unequal_segments', np.array([5.59717446, 6.36257605, 6.60517223])),
+        ('single_borehole_short', 'equivalent', 'unequal_segments', np.array([4.15784411, 4.98477603, 5.27975732])),
+        ('ten_boreholes_rectangular', 'equivalent', 'unequal_segments', np.array([10.89935004, 17.09864925, 19.0795435])),
+        # 'equivalent' solver - uniform segments
+        ('single_borehole', 'equivalent', 'uniform_segments', np.array([5.6057331, 6.37369288, 6.61659795])),
+        ('single_borehole_short', 'equivalent', 'uniform_segments', np.array([4.16941861, 4.99989722, 5.29557193])),
+        ('ten_boreholes_rectangular', 'equivalent', 'uniform_segments', np.array([10.96118694, 17.24496533, 19.2536638])),
+        #  'equivalent' solver - unequal segments, FLS approximation
+        ('single_borehole', 'equivalent', 'unequal_segments_approx', np.array([5.59717101, 6.36259907, 6.6050007])),
+        ('single_borehole_short', 'equivalent', 'unequal_segments_approx', np.array([4.15784584, 4.98478735, 5.27961509])),
+        ('ten_boreholes_rectangular', 'equivalent', 'unequal_segments_approx', np.array([10.8993464, 17.09872924, 19.0794071])),
+        #  'equivalent' solver - uniform segments, FLS approximation
+        ('single_borehole', 'equivalent', 'uniform_segments_approx', np.array([5.60572735, 6.37371464, 6.61642409])),
+        ('single_borehole_short', 'equivalent', 'uniform_segments_approx', np.array([4.16941691, 4.99990922, 5.29542863])),
+        ('ten_boreholes_rectangular', 'equivalent', 'uniform_segments_approx', np.array([10.96117468, 17.2450427, 19.25351959])),
+        # 'similarities' solver - unequal segments
+        ('single_borehole', 'similarities', 'unequal_segments', np.array([5.59717446, 6.36257605, 6.60517223])),
+        ('single_borehole_short', 'similarities', 'unequal_segments', np.array([4.15784411, 4.98477603, 5.27975732])),
+        ('ten_boreholes_rectangular', 'similarities', 'unequal_segments', np.array([10.89935004, 17.09864925, 19.0795435])),
+        # 'similarities' solver - uniform segments
+        ('single_borehole', 'similarities', 'uniform_segments', np.array([5.6057331, 6.37369288, 6.61659795])),
+        ('single_borehole_short', 'similarities', 'uniform_segments', np.array([4.16941861, 4.99989722, 5.29557193])),
+        ('ten_boreholes_rectangular', 'similarities', 'uniform_segments', np.array([10.96118694, 17.24496533, 19.2536638])),
+        # 'similarities' solver - unequal segments, FLS approximation
+        ('single_borehole', 'similarities', 'unequal_segments_approx', np.array([5.59717101, 6.36259907, 6.6050007])),
+        ('single_borehole_short', 'similarities', 'unequal_segments_approx', np.array([4.15784584, 4.98478735, 5.27961509])),
+        ('ten_boreholes_rectangular', 'similarities', 'unequal_segments_approx', np.array([10.89852244, 17.09793569, 19.07814962])),
+        # 'similarities' solver - uniform segments, FLS approximation
+        ('single_borehole', 'similarities', 'uniform_segments_approx', np.array([5.60572735, 6.37371464, 6.61642409])),
+        ('single_borehole_short', 'similarities', 'uniform_segments_approx', np.array([4.16941691, 4.99990922, 5.29542863])),
+        ('ten_boreholes_rectangular', 'similarities', 'uniform_segments_approx', np.array([10.96035847, 17.24419784, 19.25220421])),
+        # 'detailed' solver - unequal segments
+        ('single_borehole', 'detailed', 'unequal_segments', np.array([5.59717446, 6.36257605, 6.60517223])),
+        ('single_borehole_short', 'detailed', 'unequal_segments', np.array([4.15784411, 4.98477603, 5.27975732])),
+        ('ten_boreholes_rectangular', 'detailed', 'unequal_segments', np.array([10.89935004, 17.09864925, 19.0795435])),
+        # 'detailed' solver - uniform segments
+        ('single_borehole', 'detailed', 'uniform_segments', np.array([5.6057331, 6.37369288, 6.61659795])),
+        ('single_borehole_short', 'detailed', 'uniform_segments', np.array([4.16941861, 4.99989722, 5.29557193])),
+        ('ten_boreholes_rectangular', 'detailed', 'uniform_segments', np.array([10.96118694, 17.24496533, 19.2536638])),
+        # 'detailed' solver - unequal segments, FLS approximation
+        ('single_borehole', 'detailed', 'unequal_segments_approx', np.array([5.59717101, 6.36259907, 6.6050007])),
+        ('single_borehole_short', 'detailed', 'unequal_segments_approx', np.array([4.15784584, 4.98478735, 5.27961509])),
+        ('ten_boreholes_rectangular', 'detailed', 'unequal_segments_approx', np.array([10.89852244, 17.09793569, 19.07814962])),
+        # 'detailed' solver - uniform segments, FLS approximation
+        ('single_borehole', 'detailed', 'uniform_segments_approx', np.array([5.60572735, 6.37371464, 6.61642409])),
+        ('single_borehole_short', 'detailed', 'uniform_segments_approx', np.array([4.16941691, 4.99990922, 5.29542863])),
+        ('ten_boreholes_rectangular', 'detailed', 'uniform_segments_approx', np.array([10.96035847, 17.24419784, 19.25220421])),
+    ])
+def test_gfunctions_UBWT(field, method, opts, expected, request):
+    # Extract the bore field from the fixture
+    boreholes = request.getfixturevalue(field)
+    # Extract the g-function options from the fixture
+    options = request.getfixturevalue(opts)
+    # Mean borehole length [m]
+    H_mean = np.mean([b.H for b in boreholes])
+    alpha = 1e-6    # Ground thermal diffusivity [m2/s]
+    # Bore field characteristic time [s]
+    ts = H_mean**2 / (9 * alpha)
+    # Times for the g-function [s]
+    time = np.array([0.1, 1., 10.]) * ts
+    # g-Function
+    gFunc = gt.gfunction.gFunction(
+        boreholes, alpha, time=time, method=method, options=options,
+        boundary_condition='UBWT')
+    assert np.allclose(gFunc.gFunc, expected)
+
+
+# Test 'UHTR' g-functions for different bore fields using all solvers,
+# unequal/uniform segments, and with/without the FLS approximation
+@pytest.mark.parametrize("field, method, opts, expected", [
+        #  'equivalent' solver - unequal segments
+        ('single_borehole', 'equivalent', 'unequal_segments', np.array([5.61855789, 6.41336758, 6.66933682])),
+        ('single_borehole_short', 'equivalent', 'unequal_segments', np.array([4.18276733, 5.03671562, 5.34369772])),
+        ('ten_boreholes_rectangular', 'equivalent', 'unequal_segments', np.array([11.27831804, 18.48075762, 21.00669237])),
+        # 'equivalent' solver - uniform segments
+        ('single_borehole', 'equivalent', 'uniform_segments', np.array([5.61855789, 6.41336758, 6.66933682])),
+        ('single_borehole_short', 'equivalent', 'uniform_segments', np.array([4.18276733, 5.03671562, 5.34369772])),
+        ('ten_boreholes_rectangular', 'equivalent', 'uniform_segments', np.array([11.27831804, 18.48075762, 21.00669237])),
+        #  'equivalent' solver - unequal segments, FLS approximation
+        ('single_borehole', 'equivalent', 'unequal_segments_approx', np.array([5.61855411, 6.41337915, 6.66915329])),
+        ('single_borehole_short', 'equivalent', 'unequal_segments_approx', np.array([4.18276637, 5.03673008, 5.34353657])),
+        ('ten_boreholes_rectangular', 'equivalent', 'unequal_segments_approx', np.array([11.27831426, 18.48076919, 21.00650885])),
+        #  'equivalent' solver - uniform segments, FLS approximation
+        ('single_borehole', 'equivalent', 'uniform_segments_approx', np.array([5.61855411, 6.41337915, 6.66915329])),
+        ('single_borehole_short', 'equivalent', 'uniform_segments_approx', np.array([4.18276637, 5.03673008, 5.34353657])),
+        ('ten_boreholes_rectangular', 'equivalent', 'uniform_segments_approx', np.array([11.27831426, 18.48076919, 21.00650885])),
+        # 'similarities' solver - unequal segments
+        ('single_borehole', 'similarities', 'unequal_segments', np.array([5.61855789, 6.41336758, 6.66933682])),
+        ('single_borehole_short', 'similarities', 'unequal_segments', np.array([4.18276733, 5.03671562, 5.34369772])),
+        ('ten_boreholes_rectangular', 'similarities', 'unequal_segments', np.array([11.27831804, 18.48075762, 21.00669237])),
+        # 'similarities' solver - uniform segments
+        ('single_borehole', 'similarities', 'uniform_segments', np.array([5.61855789, 6.41336758, 6.66933682])),
+        ('single_borehole_short', 'similarities', 'uniform_segments', np.array([4.18276733, 5.03671562, 5.34369772])),
+        ('ten_boreholes_rectangular', 'similarities', 'uniform_segments', np.array([11.27831804, 18.48075762, 21.00669237])),
+        # 'similarities' solver - unequal segments, FLS approximation
+        ('single_borehole', 'similarities', 'unequal_segments_approx', np.array([5.61855411, 6.41337915, 6.66915329])),
+        ('single_borehole_short', 'similarities', 'unequal_segments_approx', np.array([4.18276637, 5.03673008, 5.34353657])),
+        ('ten_boreholes_rectangular', 'similarities', 'unequal_segments_approx', np.array([11.27751418, 18.47964006, 21.00475366])),
+        # 'similarities' solver - uniform segments, FLS approximation
+        ('single_borehole', 'similarities', 'uniform_segments_approx', np.array([5.61855411, 6.41337915, 6.66915329])),
+        ('single_borehole_short', 'similarities', 'uniform_segments_approx', np.array([4.18276637, 5.03673008, 5.34353657])),
+        ('ten_boreholes_rectangular', 'similarities', 'uniform_segments_approx', np.array([11.27751418, 18.47964006, 21.00475366])),
+        # 'detailed' solver - unequal segments
+        ('single_borehole', 'detailed', 'unequal_segments', np.array([5.61855789, 6.41336758, 6.66933682])),
+        ('single_borehole_short', 'detailed', 'unequal_segments', np.array([4.18276733, 5.03671562, 5.34369772])),
+        ('ten_boreholes_rectangular', 'detailed', 'unequal_segments', np.array([11.27831804, 18.48075762, 21.00669237])),
+        # 'detailed' solver - uniform segments
+        ('single_borehole', 'detailed', 'uniform_segments', np.array([5.61855789, 6.41336758, 6.66933682])),
+        ('single_borehole_short', 'detailed', 'uniform_segments', np.array([4.18276733, 5.03671562, 5.34369772])),
+        ('ten_boreholes_rectangular', 'detailed', 'uniform_segments', np.array([11.27831804, 18.48075762, 21.00669237])),
+        # 'detailed' solver - unequal segments, FLS approximation
+        ('single_borehole', 'detailed', 'unequal_segments_approx', np.array([5.61855411, 6.41337915, 6.66915329])),
+        ('single_borehole_short', 'detailed', 'unequal_segments_approx', np.array([4.18276637, 5.03673008, 5.34353657])),
+        ('ten_boreholes_rectangular', 'detailed', 'unequal_segments_approx', np.array([11.27751418, 18.47964006, 21.00475366])),
+        # 'detailed' solver - uniform segments, FLS approximation
+        ('single_borehole', 'detailed', 'uniform_segments_approx', np.array([5.61855411, 6.41337915, 6.66915329])),
+        ('single_borehole_short', 'detailed', 'uniform_segments_approx', np.array([4.18276637, 5.03673008, 5.34353657])),
+        ('ten_boreholes_rectangular', 'detailed', 'uniform_segments_approx', np.array([11.27751418, 18.47964006, 21.00475366])),
+    ])
+def test_gfunctions_UHTR(field, method, opts, expected, request):
+    # Extract the bore field from the fixture
+    boreholes = request.getfixturevalue(field)
+    # Extract the g-function options from the fixture
+    options = request.getfixturevalue(opts)
+    # Mean borehole length [m]
+    H_mean = np.mean([b.H for b in boreholes])
+    alpha = 1e-6    # Ground thermal diffusivity [m2/s]
+    # Bore field characteristic time [s]
+    ts = H_mean**2 / (9 * alpha)
+    # Times for the g-function [s]
+    time = np.array([0.1, 1., 10.]) * ts
+    # g-Function
+    gFunc = gt.gfunction.gFunction(
+        boreholes, alpha, time=time, method=method, options=options,
+        boundary_condition='UHTR')
+    assert np.allclose(gFunc.gFunc, expected)
+
+
+# Test 'MIFT' g-functions for different bore fields using all solvers,
+# unequal/uniform segments, and with/without the FLS approximation
+# The 'equivalent' solver is not applied to series-connected boreholes.
+@pytest.mark.parametrize("field, method, opts, expected", [
+        #  'equivalent' solver - unequal segments
+        ('single_borehole_network', 'equivalent', 'unequal_segments', np.array([5.76597302, 6.51058473, 6.73746895])),
+        ('single_borehole_network_short', 'equivalent', 'unequal_segments', np.array([4.17105954, 5.00930075, 5.30832133])),
+        ('ten_boreholes_network_rectangular', 'equivalent', 'unequal_segments', np.array([12.66229998, 18.57852681, 20.33535907])),
+        # 'equivalent' solver - uniform segments
+        ('single_borehole_network', 'equivalent', 'uniform_segments', np.array([5.78644676, 6.5311583, 6.75699875])),
+        ('single_borehole_network_short', 'equivalent', 'uniform_segments', np.array([4.17553236, 5.01476781, 5.31381287])),
+        ('ten_boreholes_network_rectangular', 'equivalent', 'uniform_segments', np.array([12.931553, 18.8892892, 20.63810364])),
+        #  'equivalent' solver - unequal segments, FLS approximation
+        ('single_borehole_network', 'equivalent', 'unequal_segments_approx', np.array([5.76596769, 6.51061169, 6.73731276])),
+        ('single_borehole_network_short', 'equivalent', 'unequal_segments_approx', np.array([4.17105984, 5.00931374, 5.30816983])),
+        ('ten_boreholes_network_rectangular', 'equivalent', 'unequal_segments_approx', np.array([12.66228879, 18.57863253, 20.33526092])),
+        #  'equivalent' solver - uniform segments, FLS approximation
+        ('single_borehole_network', 'equivalent', 'uniform_segments_approx', np.array([5.78644007, 6.53117706, 6.75684456])),
+        ('single_borehole_network_short', 'equivalent', 'uniform_segments_approx', np.array([4.17553118, 5.01478084, 5.31366109])),
+        ('ten_boreholes_network_rectangular', 'equivalent', 'uniform_segments_approx', np.array([12.93153466, 18.88937176, 20.63801163])),
+        # 'similarities' solver - unequal segments
+        ('single_borehole_network', 'similarities', 'unequal_segments', np.array([5.76597302, 6.51058473, 6.73746895])),
+        ('single_borehole_network_short', 'similarities', 'unequal_segments', np.array([4.17105954, 5.00930075, 5.30832133])),
+        ('ten_boreholes_network_rectangular', 'similarities', 'unequal_segments', np.array([12.66229998, 18.57852681, 20.33535907])),
+        ('ten_boreholes_network_rectangular_series', 'similarities', 'unequal_segments', np.array([3.19050169, 8.8595362, 10.84379419])),
+        # 'similarities' solver - uniform segments
+        ('single_borehole_network', 'similarities', 'uniform_segments', np.array([5.78644676, 6.5311583, 6.75699875])),
+        ('single_borehole_network_short', 'similarities', 'uniform_segments', np.array([4.17553236, 5.01476781, 5.31381287])),
+        ('ten_boreholes_network_rectangular', 'similarities', 'uniform_segments', np.array([12.931553, 18.8892892, 20.63810364])),
+        ('ten_boreholes_network_rectangular_series', 'similarities', 'uniform_segments', np.array([3.22186337, 8.92628494, 10.91644607])),
+        # 'similarities' solver - unequal segments, FLS approximation
+        ('single_borehole_network', 'similarities', 'unequal_segments_approx', np.array([5.76596769, 6.51061169, 6.73731276])),
+        ('single_borehole_network_short', 'similarities', 'unequal_segments_approx', np.array([4.17105984, 5.00931374, 5.30816983])),
+        ('ten_boreholes_network_rectangular', 'similarities', 'unequal_segments_approx', np.array([12.66136057, 18.57792276, 20.33429034])),
+        ('ten_boreholes_network_rectangular_series', 'similarities', 'unequal_segments_approx', np.array([3.19002567, 8.85783554, 10.84222402])),
+        # 'similarities' solver - uniform segments, FLS approximation
+        ('single_borehole_network', 'similarities', 'uniform_segments_approx', np.array([5.78644007, 6.53117706, 6.75684456])),
+        ('single_borehole_network_short', 'similarities', 'uniform_segments_approx', np.array([4.17553118, 5.01478084, 5.31366109])),
+        ('ten_boreholes_network_rectangular', 'similarities', 'uniform_segments_approx', np.array([12.93064329, 18.88844718, 20.63710493])),
+        ('ten_boreholes_network_rectangular_series', 'similarities', 'uniform_segments_approx', np.array([3.22139028, 8.92448715, 10.91485947])),
+        # 'detailed' solver - unequal segments
+        ('single_borehole_network', 'detailed', 'unequal_segments', np.array([5.76597302, 6.51058473, 6.73746895])),
+        ('single_borehole_network_short', 'detailed', 'unequal_segments', np.array([4.17105954, 5.00930075, 5.30832133])),
+        ('ten_boreholes_network_rectangular', 'detailed', 'unequal_segments', np.array([12.66229998, 18.57852681, 20.33535907])),
+        ('ten_boreholes_network_rectangular_series', 'detailed', 'unequal_segments', np.array([3.19050169, 8.8595362, 10.84379419])),
+        # 'detailed' solver - uniform segments
+        ('single_borehole_network', 'detailed', 'uniform_segments', np.array([5.78644676, 6.5311583, 6.75699875])),
+        ('single_borehole_network_short', 'detailed', 'uniform_segments', np.array([4.17553236, 5.01476781, 5.31381287])),
+        ('ten_boreholes_network_rectangular', 'detailed', 'uniform_segments', np.array([12.931553, 18.8892892, 20.63810364])),
+        ('ten_boreholes_network_rectangular_series', 'detailed', 'uniform_segments', np.array([3.22186337, 8.92628494, 10.91644607])),
+        # 'detailed' solver - unequal segments, FLS approximation
+        ('single_borehole_network', 'detailed', 'unequal_segments_approx', np.array([5.76596769, 6.51061169, 6.73731276])),
+        ('single_borehole_network_short', 'detailed', 'unequal_segments_approx', np.array([4.17105984, 5.00931374, 5.30816983])),
+        ('ten_boreholes_network_rectangular', 'detailed', 'unequal_segments_approx', np.array([12.66136057, 18.57792276, 20.33429034])),
+        ('ten_boreholes_network_rectangular_series', 'detailed', 'unequal_segments_approx', np.array([3.19002567, 8.85783554, 10.84222402])),
+        # 'detailed' solver - uniform segments, FLS approximation
+        ('single_borehole_network', 'detailed', 'uniform_segments_approx', np.array([5.78644007, 6.53117706, 6.75684456])),
+        ('single_borehole_network_short', 'detailed', 'uniform_segments_approx', np.array([4.17553118, 5.01478084, 5.31366109])),
+        ('ten_boreholes_network_rectangular', 'detailed', 'uniform_segments_approx', np.array([12.93064329, 18.88844718, 20.63710493])),
+        ('ten_boreholes_network_rectangular_series', 'detailed', 'uniform_segments_approx', np.array([3.22139028, 8.92448715, 10.91485947])),
+    ])
+def test_gfunctions_MIFT(field, method, opts, expected, request):
+    # Extract the bore field from the fixture
+    network = request.getfixturevalue(field)
+    # Extract the g-function options from the fixture
+    options = request.getfixturevalue(opts)
+    # Mean borehole length [m]
+    H_mean = np.mean([b.H for b in network.b])
+    alpha = 1e-6    # Ground thermal diffusivity [m2/s]
+    # Bore field characteristic time [s]
+    ts = H_mean**2 / (9 * alpha)
+    # Times for the g-function [s]
+    time = np.array([0.1, 1., 10.]) * ts
+    # g-Function
+    gFunc = gt.gfunction.gFunction(
+        network, alpha, time=time, method=method, options=options,
+        boundary_condition='MIFT')
+    assert np.allclose(gFunc.gFunc, expected)
diff --git a/tests/heat_transfer_test.py b/tests/heat_transfer_test.py
index b6c4ab6f..6fa36b2b 100644
--- a/tests/heat_transfer_test.py
+++ b/tests/heat_transfer_test.py
@@ -1,126 +1,289 @@
 # -*- coding: utf-8 -*-
 """ Test suite for heat_transfer module.
 """
-import unittest
+import pytest
 
 import numpy as np
-from scipy.integrate import dblquad
-from scipy.special import erfc
 
+import pygfunction as gt
 
-class TestFiniteLineSource(unittest.TestCase):
-    """ Test cases for finite_line_source function.
-    """
 
-    def setUp(self):
-        self.t = 1. * 8760. * 3600.     # Time is 1 year
-        self.alpha = 1.0e-6             # Thermal diffusivity
-        self.D1 = 4.0                   # Buried depth of source
-        self.D2 = 16.0                  # Buried depth of target
-        self.H1 = 10.0                  # Length of source
-        self.H2 = 7.0                   # Length of target
-        self.dis = 12.0                 # Distance of target
+# =============================================================================
+# Test finite_line_source
+# =============================================================================
+# Test finite_line_source for single borehole to single borehole solution over
+# a single time step, using different combinations of reaSource, imgSource, and
+# with/without the FLS approximation
+@pytest.mark.parametrize("borehole1, borehole2, reaSource, imgSource, approximation, N, expected", [
+        # Same borehole
+        ('single_borehole', 'single_borehole', True, True, False, 10, 6.41336757901249),
+        ('single_borehole', 'single_borehole', True, False, False, 10, 6.531908395192002),
+        ('single_borehole', 'single_borehole', False, True, False, 10, -0.11854082),
+        ('single_borehole', 'single_borehole', False, False, False, 10, 0.),
+        # Same borehole - FLS approximation
+        ('single_borehole', 'single_borehole', True, True, True, 10, 6.41337915),
+        ('single_borehole', 'single_borehole', True, False, True, 10, 6.53191379),
+        ('single_borehole', 'single_borehole', False, True, True, 10, -0.11853464),
+        ('single_borehole', 'single_borehole', False, False, True, 10, 0.),
+        # Two boreholes
+        ('single_borehole', 'single_borehole_short', True, True, False, 10, 1.62660428),
+        ('single_borehole', 'single_borehole_short', True, False, False, 10, 1.99622241),
+        ('single_borehole', 'single_borehole_short', False, True, False, 10, -0.36961813),
+        ('single_borehole', 'single_borehole_short', False, False, False, 10, 0.),
+        # Two boreholes - FLS approximation
+        ('single_borehole', 'single_borehole_short', True, True, True, 10, 1.62655937),
+        ('single_borehole', 'single_borehole_short', True, False, True, 10, 1.99624209),
+        ('single_borehole', 'single_borehole_short', False, True, True, 10, -0.36968273),
+        ('single_borehole', 'single_borehole_short', False, False, True, 10, 0.),
+    ])
+def test_finite_line_source_one_to_one_single_time_step(
+        borehole1, borehole2, reaSource, imgSource, approximation, N, expected,
+        request):
+        # Extract boreholes from fixtures
+        b1 = request.getfixturevalue(borehole1)[0]
+        b2 = request.getfixturevalue(borehole2)[0]
+        alpha = 1e-6                # Ground thermal diffusivity [m2/s]
+        ts = b1.H**2 / (9 * alpha)  # Borehole characteristic time [s]
+        # Time for FLS calculation [s]
+        time = ts
+        # Evaluate FLS
+        h = gt.heat_transfer.finite_line_source(
+            time, alpha, b1, b2, reaSource=reaSource, imgSource=imgSource,
+            approximation=approximation, N=N)
+        assert np.isclose(h, expected)
 
-    def test_finite_line_source(self, rel_tol=1.0e-6):
-        """ Tests the value of the FLS solution.
-        """
-        from pygfunction.boreholes import Borehole
-        from pygfunction.heat_transfer import finite_line_source
-        # Evaluate the double integral
-        reference = dblquad(fls_double,
-                            self.D1, self.D1+self.H1,
-                            lambda x: self.D2, lambda x: self.D2+self.H2,
-                            args=(self.t, self.dis, self.alpha))[0]/self.H2
-        # Evaluate using heat_transfer.finite_line_source
-        borehole1 = Borehole(self.H1, self.D1, 0.05, 0., 0.)
-        borehole2 = Borehole(self.H2, self.D2, 0.05, self.dis, 0.)
-        calculated = finite_line_source(self.t, self.alpha,
-                                        borehole1, borehole2)
-        self.assertAlmostEqual(calculated, reference,
-                               delta=rel_tol*reference,
-                               msg='Incorrect value of finite line source '
-                                   'solution.')
 
-    def test_finite_line_source_real_part(self, rel_tol=1.0e-6):
-        """ Tests the value of the real part of the FLS solution.
-        """
-        from pygfunction.boreholes import Borehole
-        from pygfunction.heat_transfer import finite_line_source
-        # Evaluate the double integral
-        reference = dblquad(fls_double,
-                            self.D1, self.D1+self.H1,
-                            lambda x: self.D2, lambda x: self.D2+self.H2,
-                            args=(self.t,
-                                  self.dis,
-                                  self.alpha,
-                                  True,
-                                  False))[0]/self.H2
-        # Evaluate using heat_transfer.finite_line_source
-        borehole1 = Borehole(self.H1, self.D1, 0.05, 0., 0.)
-        borehole2 = Borehole(self.H2, self.D2, 0.05, self.dis, 0.)
-        calculated = finite_line_source(self.t, self.alpha,
-                                        borehole1, borehole2,
-                                        reaSource=True, imgSource=False)
-        self.assertAlmostEqual(calculated, reference,
-                               delta=rel_tol*reference,
-                               msg='Incorrect value of the real part of the '
-                                   'finite line source solution.')
+# Test finite_line_source for single borehole to single borehole solution over
+# a multiple time steps, using different combinations of reaSource, imgSource,
+# and with/without the FLS approximation
+@pytest.mark.parametrize("borehole1, borehole2, reaSource, imgSource, approximation, N, expected", [
+        # Same borehole
+        ('single_borehole', 'single_borehole', True, True, False, 10, np.array([5.61855789, 6.41336758, 6.66933682])),
+        ('single_borehole', 'single_borehole', True, False, False, 10, np.array([5.63709751, 6.5319084, 7.0301884])),
+        ('single_borehole', 'single_borehole', False, True, False, 10, np.array([-0.01853962, -0.11854082, -0.36085158])),
+        ('single_borehole', 'single_borehole', False, False, False, 10, np.array([0., 0., 0.])),
+        # Same borehole - FLS approximation
+        ('single_borehole', 'single_borehole', True, True, True, 10, np.array([5.61855411, 6.41337915, 6.66915329])),
+        ('single_borehole', 'single_borehole', True, False, True, 10, np.array([5.63709751, 6.53191379, 7.0301115])),
+        ('single_borehole', 'single_borehole', False, True, True, 10, np.array([-0.01854339, -0.11853464, -0.3609582])),
+        ('single_borehole', 'single_borehole', False, False, True, 10, np.array([0., 0., 0.])),
+        # Two boreholes
+        ('single_borehole', 'single_borehole_short', True, True, False, 10, np.array([1.1703042, 1.62660428, 1.72510234])),
+        ('single_borehole', 'single_borehole_short', True, False, False, 10, np.array([1.24143841, 1.99622241, 2.46653331])),
+        ('single_borehole', 'single_borehole_short', False, True, False, 10, np.array([-0.07113422, -0.36961813, -0.74143097])),
+        ('single_borehole', 'single_borehole_short', False, False, False, 10, np.array([0., 0., 0.])),
+        # Two boreholes - FLS approximation
+        ('single_borehole', 'single_borehole_short', True, True, True, 10, np.array([1.17008899, 1.62655937, 1.72496239])),
+        ('single_borehole', 'single_borehole_short', True, False, True, 10, np.array([1.24123405, 1.99624209, 2.46638582])),
+        ('single_borehole', 'single_borehole_short', False, True, True, 10, np.array([-0.07114505, -0.36968273, -0.74142343])),
+        ('single_borehole', 'single_borehole_short', False, False, True, 10, np.array([0., 0., 0.])),
+    ])
+def test_finite_line_source_one_to_one_multiple_time_steps(
+        borehole1, borehole2, reaSource, imgSource, approximation, N, expected,
+        request):
+        # Extract boreholes from fixtures
+        b1 = request.getfixturevalue(borehole1)[0]
+        b2 = request.getfixturevalue(borehole2)[0]
+        alpha = 1e-6                # Ground thermal diffusivity [m2/s]
+        ts = b1.H**2 / (9 * alpha)  # Borehole characteristic time [s]
+        # Times for FLS calculation [s]
+        time = np.array([0.1, 1., 10.]) * ts
+        # Evaluate FLS
+        h = gt.heat_transfer.finite_line_source(
+            time, alpha, b1, b2, reaSource=reaSource, imgSource=imgSource,
+            approximation=approximation, N=N)
+        assert np.allclose(h, expected)
 
-    def test_finite_line_source_image_part(self, rel_tol=1.0e-6):
-        """ Tests the value of the image part of the FLS solution.
-        """
-        from pygfunction.boreholes import Borehole
-        from pygfunction.heat_transfer import finite_line_source
-        # Evaluate the double integral
-        reference = dblquad(fls_double,
-                            self.D1, self.D1+self.H1,
-                            lambda x: self.D2, lambda x: self.D2+self.H2,
-                            args=(self.t,
-                                  self.dis,
-                                  self.alpha,
-                                  False,
-                                  True))[0]/self.H2
-        # Evaluate using heat_transfer.finite_line_source
-        borehole1 = Borehole(self.H1, self.D1, 0.05, 0., 0.)
-        borehole2 = Borehole(self.H2, self.D2, 0.05, self.dis, 0.)
-        calculated = finite_line_source(self.t, self.alpha,
-                                        borehole1, borehole2,
-                                        reaSource=False, imgSource=True)
-        self.assertAlmostEqual(calculated, reference,
-                               delta=np.abs(rel_tol*reference),
-                               msg='Incorrect value of the image part of the '
-                                   'finite line source solution.')
 
-    def test_finite_line_source_no_part(self, rel_tol=1.0e-6):
-        """ Tests the value of the FLS solution when considering no source.
-        """
-        from pygfunction.boreholes import Borehole
-        from pygfunction.heat_transfer import finite_line_source
-        # Evaluate using heat_transfer.finite_line_source
-        borehole1 = Borehole(self.H1, self.D1, 0.05, 0., 0.)
-        borehole2 = Borehole(self.H2, self.D2, 0.05, self.dis, 0.)
-        calculated = finite_line_source(self.t, self.alpha,
-                                        borehole1, borehole2,
-                                        reaSource=False, imgSource=False)
-        self.assertEqual(calculated, 0.,
-                               msg='Incorrect value of no part of the '
-                                   'finite line source solution.')
+# Test finite_line_source for single borehole to single borehole solution at
+# steady-state, using different combinations of reaSource, imgSource
+@pytest.mark.parametrize("borehole1, borehole2, reaSource, imgSource, expected", [
+        # Same borehole
+        ('single_borehole', 'single_borehole', True, True, 6.68879600214517),
+        ('single_borehole', 'single_borehole', True, False, 7.29454957717299),
+        ('single_borehole', 'single_borehole', False, True, -0.6057535750278188),
+        ('single_borehole', 'single_borehole', False, False, 0.),
+        # Two boreholes
+        ('single_borehole', 'single_borehole_short', True, True, 1.7307294312317285),
+        ('single_borehole', 'single_borehole_short', True, False, 2.729558982105579),
+        ('single_borehole', 'single_borehole_short', False, True, -0.9988295508738515),
+        ('single_borehole', 'single_borehole_short', False, False, 0.),
+    ])
+def test_finite_line_source_one_to_one_steady_state(
+        borehole1, borehole2, reaSource, imgSource, expected, request):
+        b1 = request.getfixturevalue(borehole1)[0]
+        b2 = request.getfixturevalue(borehole2)[0]
+        alpha = 1e-6     # Ground thermal diffusivity [m2/s]
+        # Time for FLS calculation [s]
+        time = np.inf
+        # Evaluate FLS
+        h = gt.heat_transfer.finite_line_source(
+            time, alpha, b1, b2, reaSource=reaSource, imgSource=imgSource)
+        assert np.isclose(h, expected)
 
 
-def fls_double(z2, z1, t, dis, alpha, reaSource=True, imgSource=True):
-    """ FLS expression for double integral solution.
-    """
-    r_pos = np.sqrt(dis**2 + (z2 - z1)**2)
-    r_neg = np.sqrt(dis**2 + (z2 + z1)**2)
-    fls = 0.
-    if reaSource:
-        fls += 0.5*erfc(r_pos/np.sqrt(4*alpha*t))/r_pos
-    if imgSource:
-        fls += -0.5*erfc(r_neg/np.sqrt(4*alpha*t))/r_neg
-    return fls
+# Test finite_line_source for multiple boreholes to multiple boreholes solution
+# over a single time step, using different combinations of reaSource,
+# imgSource, and with/without the FLS approximation
+@pytest.mark.parametrize("reaSource, imgSource, approximation, N, expected", [
+        # Without FLS approximation
+        (True, True, False, 10, np.array(
+            [[6.37361733, 1.00243165, 1.75179834],
+             [1.70869031, 6.17773914, 1.78839136],
+             [1.48457487, 0.88914372, 6.32534129]])),
+        (True, False, False, 10, np.array(
+            [[6.47984031, 1.11115916, 1.8478299],
+             [1.8940213, 6.39079975, 1.96672666],
+             [1.56595754, 0.9778076, 6.40728418]])),
+        (False, True, False, 10, np.array(
+            [[-0.10622298, -0.10872751, -0.09603156],
+             [-0.18533099, -0.21306061, -0.1783353],
+             [-0.08138268, -0.08866388, -0.0819429]])),
+        (False, False, False, 10, np.zeros((3, 3))),
+        # With FLS approximation
+        (True, True, True, 10, np.array(
+            [[6.37363678, 1.00243611, 1.75176654],
+             [1.70869791, 6.17777608, 1.7884002],
+             [1.48454791, 0.88914812, 6.32536189]])),
+        (True, False, True, 10, np.array(
+            [[6.4798501, 1.11115569, 1.84778942],
+             [1.89401537, 6.3908279, 1.96672028],
+             [1.56592324, 0.97780443, 6.40729653]])),
+        (False, True, True, 10, np.array(
+            [[-0.10621332, -0.10871958, -0.09602288],
+             [-0.18531747, -0.21305182, -0.17832008],
+             [-0.08137532, -0.08865631, -0.08193464]])),
+        (False, False, True, 10, np.zeros((3, 3))),
+    ])
+def test_finite_line_source_multiple_boreholes_single_time_step(
+        three_boreholes_unequal, reaSource, imgSource, approximation, N,
+        expected):
+        # Extract boreholes from fixture
+        boreholes = three_boreholes_unequal
+        alpha = 1e-6    # Ground thermal diffusivity [m2/s]
+        # Bore field characteristic time [s]
+        ts = np.mean([b.H for b in boreholes])**2 / (9 * alpha)
+        # Time for FLS calculation [s]
+        time = ts
+        # Evaluate FLS
+        h = gt.heat_transfer.finite_line_source(
+            time, alpha, boreholes, boreholes, reaSource=reaSource,
+            imgSource=imgSource, approximation=approximation, N=N)
+        assert np.allclose(h, expected)
 
 
-if __name__ == '__main__' and __package__ is None:
-    from os import sys, path
-    sys.path.append(path.dirname(path.dirname(path.abspath(__file__))))
-    unittest.main()
+# Test finite_line_source for multiple boreholes to multiple boreholes solution
+# over multiple time steps, using different combinations of reaSource,
+# imgSource, and with/without the FLS approximation
+@pytest.mark.parametrize("reaSource, imgSource, approximation, N, expected", [
+        # Without FLS approximation
+        (True, True, False, 10, np.array(
+            [[[5.54962392, 6.37361733, 6.66453291],
+              [0.56736312, 1.00243165, 1.13057208],
+              [0.88843196, 1.75179834, 2.10314226]],
+             [[0.96709622, 1.70869031, 1.9271115],
+              [5.58527306, 6.17773914, 6.28409069],
+              [1.03296293, 1.78839136, 2.04396068]],
+             [[0.75290844, 1.48457487, 1.78232395],
+              [0.51356349, 0.88914372, 1.01620644],
+              [5.44663175, 6.32534129, 6.69487339]]])),
+        (True, False, False, 10, np.array(
+            [[[5.5653799, 6.47984031, 7.00849476],
+              [0.58387558, 1.11115916, 1.42217343],
+              [0.89939015, 1.8478299, 2.44818701]],
+             [[0.99524247, 1.8940213, 2.42415926],
+              [5.63143639, 6.39079975, 6.73519825],
+              [1.05867673, 1.96672666, 2.55106101]],
+             [[0.76219504, 1.56595754, 2.07473475],
+              [0.52634775, 0.9778076, 1.26832412],
+              [5.45693032, 6.40728418, 6.99916997]]])),
+        (False, True, False, 10, np.array(
+            [[[-0.01575598, -0.10622298, -0.34396185],
+              [-0.01651247, -0.10872751, -0.29160135],
+              [-0.01095819, -0.09603156, -0.34504475]],
+             [[-0.02814625, -0.18533099, -0.49704776],
+              [-0.04616333, -0.21306061, -0.45110757],
+              [-0.02571379, -0.1783353, -0.50710034]],
+             [[-0.0092866, -0.08138268, -0.2924108],
+              [-0.01278426, -0.08866388, -0.25211768],
+              [-0.01029857, -0.0819429, -0.30429658]]])),
+        (False, False, False, 10, np.zeros((3, 3, 3))),
+        # With FLS approximation
+        (True, True, True, 10, np.array(
+            [[[5.54962044, 6.37363678, 6.66437927],
+              [0.56717846, 1.00243611, 1.13049774],
+              [0.88813623, 1.75176654, 2.10302321]],
+             [[0.96678146, 1.70869791, 1.92698479],
+              [5.58527022, 6.17777608, 6.28433246],
+              [1.03266465, 1.7884002, 2.04390268]],
+             [[0.75265782, 1.48454791, 1.78222306],
+              [0.51341519, 0.88914812, 1.0161776],
+              [5.44663013, 6.32536189, 6.69481429]]])),
+        (True, False, True, 10, np.array(
+            [[[5.5653799, 6.4798501, 7.00842175],
+              [0.58369262, 1.11115569, 1.42213585],
+              [0.8990927, 1.84778942, 2.44811819]],
+             [[0.9949306, 1.89401537, 2.42409521],
+              [5.63143661, 6.3908279, 6.73530821],
+              [1.05838072, 1.96672028, 2.55103547]],
+             [[0.76194297, 1.56592324, 2.07467643],
+              [0.52620058, 0.97780443, 1.26831142],
+              [5.45693032, 6.40729653, 6.99913332]]])),
+        (False, True, True, 10, np.array(
+            [[[-0.01575946, -0.10621332, -0.34404249],
+              [-0.01651416, -0.10871958, -0.29163811],
+              [-0.01095647, -0.09602288, -0.34509498]],
+             [[-0.02814914, -0.18531747, -0.49711042],
+              [-0.0461664, -0.21305182, -0.45097575],
+              [-0.02571607, -0.17832008, -0.50713279]],
+             [[-0.00928515, -0.08137532, -0.29245337],
+              [-0.01278539, -0.08865631, -0.25213381],
+              [-0.01030019, -0.08193464, -0.30431902]]])),
+        (False, False, True, 10, np.zeros((3, 3, 3))),
+    ])
+def test_finite_line_source_multiple_boreholes_multiple_time_steps(
+        three_boreholes_unequal, reaSource, imgSource, approximation, N,
+        expected):
+        # Extract boreholes from fixture
+        boreholes = three_boreholes_unequal
+        alpha = 1e-6    # Ground thermal diffusivity [m2/s]
+        # Bore field characteristic time [s]
+        ts = np.mean([b.H for b in boreholes])**2 / (9 * alpha)
+        # Times for FLS calculation [s]
+        time = np.array([0.1, 1., 10.]) * ts
+        # Evaluate FLS
+        h = gt.heat_transfer.finite_line_source(
+            time, alpha, boreholes, boreholes, reaSource=reaSource,
+            imgSource=imgSource, approximation=approximation, N=N)
+        assert np.allclose(h, expected)
+
+
+# Test finite_line_source for multiple boreholes to multiple boreholes solution
+# at steady-state, using different combinations of reaSource, imgSource
+@pytest.mark.parametrize("reaSource, imgSource, expected", [
+        (True, True, np.array(
+            [[6.688796, 1.13927583, 2.13603998],
+             [1.94194744, 6.28943104, 2.06405893],
+             [1.81020337, 1.02619879, 6.73274814]])),
+        (True, False, np.array(
+            [[7.29454958, 1.59004241, 2.78444571],
+             [2.71029956, 6.9045905, 2.8863451],
+             [2.35969976, 1.43501903, 7.3348811]])),
+        (False, True, np.array(
+            [[-0.60575358, -0.45076657, -0.64840574],
+             [-0.76835211, -0.61515946, -0.82228616],
+             [-0.54949639, -0.40882024, -0.60213297]])),
+        (False, False, np.zeros((3, 3))),
+    ])
+def test_finite_line_source_multiple_boreholes_steady_state(
+        three_boreholes_unequal, reaSource, imgSource, expected):
+        # Extract boreholes from fixture
+        boreholes = three_boreholes_unequal
+        alpha = 1e-6    # Ground thermal diffusivity [m2/s]
+        # Time for FLS calculation [s]
+        time = np.inf
+        # Evaluate FLS
+        h = gt.heat_transfer.finite_line_source(
+            time, alpha, boreholes, boreholes, reaSource=reaSource,
+            imgSource=imgSource)
+        assert np.allclose(h, expected)
diff --git a/tests/pipes_test.py b/tests/pipes_test.py
index eef2ae57..2dd903c1 100644
--- a/tests/pipes_test.py
+++ b/tests/pipes_test.py
@@ -1,662 +1,538 @@
 # -*- coding: utf-8 -*-
 """ Test suite for pipes module.
 """
-import unittest
+import pytest
 
 import numpy as np
-from scipy.constants import pi
-
-
-class TestFluidFrictionFactor(unittest.TestCase):
-    """ Test cases for calculation of friction factor in circular pipes.
-    """
-
-    def setUp(self):
-        self.r_in = 0.01            # Inner radius [m]
-        self.epsilon = 1.5e-6       # Pipe surface roughness [m]
-        self.visc = 0.00203008      # Fluid viscosity [kg/m.s]
-        self.den = 1014.78          # Fluid density
-
-    def test_laminar(self, rel_tol=1.0e-3):
-        """ Tests the value of the Darcy friction factor in laminar flow.
-        """
-        from pygfunction.pipes import fluid_friction_factor_circular_pipe
-        m_flow = 0.05               # Fluid mass flow rate
-
-        # Result from EES
-        reference = 0.04082
-        # Calculate using pipes.fluid_friction_factor_circular_pipe()
-        f = fluid_friction_factor_circular_pipe(m_flow, self.r_in, self.visc,
-                                                self.den, self.epsilon)
-        self.assertAlmostEqual(f, reference, delta=rel_tol*reference,
-                               msg='Incorrect value of Darcy friction '
-                                   'factor for laminar flow.')
-
-    def test_turbulent_rough(self, rel_tol=1.0e-3):
-        """ Tests the value of the Darcy friction factor in turbulent flow
-            in rough pipes.
-        """
-        from pygfunction.pipes import fluid_friction_factor_circular_pipe
-        m_flow = 0.50               # Fluid mass flow rate
-
-        # Result from EES
-        reference = 0.02766
-        # Calculate using pipes.fluid_friction_factor_circular_pipe()
-        f = fluid_friction_factor_circular_pipe(m_flow, self.r_in, self.visc,
-                                                self.den, self.epsilon,
-                                                tol=1.0e-8)
-        self.assertAlmostEqual(f, reference, delta=rel_tol*reference,
-                               msg='Incorrect value of Darcy friction '
-                                   'factor for turbulent flow in rough pipes.')
-
-
-class TestConvectiveHeatTransferCoefficient(unittest.TestCase):
-    """ Test cases for calculation of convection coefficient in circular pipes.
-    """
-
-    def setUp(self):
-        self.r_in = 0.01            # Inner radius [m]
-        self.epsilon = 1.5e-6       # Pipe surface roughness [m]
-        self.visc = 0.00203008      # Fluid viscosity [kg/m.s]
-        self.den = 1014.78          # Fluid density [kg/m3]
-        self.cp = 3977.             # Fluid specific heat capacity [J/kg.K]
-        self.k = 0.4922             # Fluid thermal conductivity [W/m.K]
-
-    def test_turbulent(self, rel_tol=1.0e-3):
-        """ Tests the value of the convection coefficient in laminar flow.
-        """
-        from pygfunction.pipes import convective_heat_transfer_coefficient_circular_pipe
-        m_flow = 0.50               # Fluid mass flow rate
-
-        # Result from EES
-        reference = 4037.58769
-        # Calculate using pipes.convective_heat_transfer_coefficient_circular_pipe()
-        h = convective_heat_transfer_coefficient_circular_pipe(m_flow,
-                                                               self.r_in,
-                                                               self.visc,
-                                                               self.den,
-                                                               self.k,
-                                                               self.cp,
-                                                               self.epsilon)
-        self.assertAlmostEqual(h, reference, delta=rel_tol*reference,
-                               msg='Incorrect value of the convection coeff. '
-                                   'for turbulent flow.')
-
-
-class TestPipeThermalResistance(unittest.TestCase):
-    """ Test cases for calculation of the thermal resistance through the pipe
-        wall.
-    """
-
-    def setUp(self):
-        self.r_in = 0.01            # Inner radius [m]
-        self.r_out = 0.02           # Outer radius [m]
-        self.k = 0.6                # Fluid thermal conductivity [W/m.K]
-
-    def test_turbulent(self, rel_tol=1.0e-6):
-        """ Tests the value of the convection coefficient in laminar flow.
-        """
-        from pygfunction.pipes import conduction_thermal_resistance_circular_pipe
-
-        # Exact solution
-        reference = np.log(self.r_out/self.r_in)/(2*pi*self.k)
-        # Calculate using pipes.conduction_thermal_resistance_circular_pipe()
-        R = conduction_thermal_resistance_circular_pipe(self.r_in, self.r_out,
-                                                        self.k)
-        self.assertAlmostEqual(R, reference, delta=rel_tol*reference,
-                               msg='Incorrect value of the thermal resistance '
-                                   'through pipe wall.')
-
-
-class TestBoreholeThermalResistances(unittest.TestCase):
-    """ Test cases for calculation of the internal thermal resistances of a
-        borehole.
-    """
-
-    def setUp(self):
-        # Pipe positions [m]
-        self.pos_2pipes = [(0.03, 0.00), (-0.03, 0.02)]
-        self.r_out = 0.02       # Pipe outer radius [m]
-        self.r_b = 0.07         # Borehole radius [m]
-        self.k_s = 2.5          # Ground thermal conductivity [W/m.K]
-        self.k_g = 1.5          # Grout thermal conductivity [W/m.K]
-        beta = 1.2
-        # Fluid to outer pipe wall thermal resistance [m.K/W]
-        self.Rfp = beta/(2*pi*self.k_g)
-
-    def test_line_source_approximation(self):
-        """ Tests the value of the internal resistances and conductances.
-        """
-        from pygfunction.pipes import thermal_resistances
-
-        # Reference solution (Claesson and Hellstrom, 2011)
-        R_ref = np.array([[25.486e-2, 01.538e-2],
-                          [01.538e-2, 25.207e-2]])
-        S_ref = np.array([[3.698, 0.240],
-                          [0.240, 3.742]])
-        # Calculate using pipes.fluid_friction_factor_circular_pipe()
-        R, Rd = thermal_resistances(self.pos_2pipes, self.r_out, self.r_b,
-                                    self.k_s, self.k_g, self.Rfp, J=0)
-        S = 1/Rd
-        self.assertTrue(np.allclose(R, R_ref, rtol=1e-8, atol=1e-5),
-                        msg='Incorrect value of the internal thermal '
-                            'resistances.')
-        self.assertTrue(np.allclose(S, S_ref, rtol=1e-8, atol=1e-3),
-                        msg='Incorrect value of the delta-circuit thermal '
-                            'resistances.')
-
-    def test_first_order_multipole(self):
-        """ Tests the value of the internal resistances and conductances.
-        """
-        from pygfunction.pipes import thermal_resistances
-
-        # Reference solution (Claesson and Hellstrom, 2011)
-        R_ref = np.array([[25.569e-2, 01.562e-2],
-                          [01.562e-2, 25.288e-2]])
-        S_ref = np.array([[3.683, 0.243],
-                          [0.243, 3.727]])
-        # Calculate using pipes.fluid_friction_factor_circular_pipe()
-        R, Rd = thermal_resistances(self.pos_2pipes, self.r_out, self.r_b,
-                                    self.k_s, self.k_g, self.Rfp, J=1)
-        S = 1/Rd
-        self.assertTrue(np.allclose(R, R_ref, rtol=1e-8, atol=1e-5),
-                        msg='Incorrect value of the internal thermal '
-                            'resistances.')
-        self.assertTrue(np.allclose(S, S_ref, rtol=1e-8, atol=1e-3),
-                        msg='Incorrect value of the delta-circuit thermal '
-                            'resistances.')
-
-    def test_second_order_multipole(self):
-        """ Tests the value of the internal resistances and conductances.
-        """
-        from pygfunction.pipes import thermal_resistances
-
-        # Reference solution (Claesson and Hellstrom, 2011)
-        R_ref = np.array([[25.590e-2, 01.561e-2],
-                          [01.561e-2, 25.309e-2]])
-        S_ref = np.array([[3.681, 0.242],
-                          [0.242, 3.724]])
-        # Calculate using pipes.fluid_friction_factor_circular_pipe()
-        R, Rd = thermal_resistances(self.pos_2pipes, self.r_out, self.r_b,
-                                    self.k_s, self.k_g, self.Rfp, J=2)
-        S = 1/Rd
-        self.assertTrue(np.allclose(R, R_ref, rtol=1e-8, atol=1e-5),
-                        msg='Incorrect value of the internal thermal '
-                            'resistances.')
-        self.assertTrue(np.allclose(S, S_ref, rtol=1e-8, atol=1e-3),
-                        msg='Incorrect value of the delta-circuit thermal '
-                            'resistances.')
-
-
-class TestMultipole(unittest.TestCase):
-    """ Test cases for multipole function.
-    """
-
-    def setUp(self):
-        # Pipe positions [m]
-        self.pos_4pipes = [(0.03, 0.03), (-0.03, 0.03),
-                           (-0.03, -0.03), (0.03, -0.03)]
-        self.r_out = 0.02       # Pipe outer radius [m]
-        self.r_b = 0.07         # Borehole radius [m]
-        self.k_s = 2.5          # Ground thermal conductivity [W/m.K]
-        self.k_g = 1.5          # Grout thermal conductivity [W/m.K]
-        self.T_b = 0.0          # Borehole wall temperature [degC]
-        # Pipe heat transfer rates [W/m]
-        self.Q_p = np.array([10., 9., 8., 7.])
-        self.J = 3              # Number of multipole per pipe
-        beta = 1.2
-        # Fluid to outer pipe wall thermal resistance [m.K/W]
-        self.Rfp = beta/(2*pi*self.k_g)
-
-    def test_third_order_fluid_temperatures_four_pipes(self):
-        """ Tests the value of the fluid temperatures.
-        """
-        from pygfunction.pipes import multipole
-
-        # Reference solution (Claesson, 2012)
-        reference = np.array([2.719532, 2.519227, 2.193383, 1.993078])
-        # Calculate using pipes.multipole()
-        Tf = multipole(self.pos_4pipes, self.r_out, self.r_b,
-                       self.k_s, self.k_g, self.Rfp,
-                       self.T_b, self.Q_p, self.J)[0]
-        self.assertTrue(np.allclose(Tf, reference, rtol=1e-4, atol=1e-10),
-                        msg='Incorrect value of fluid temperatures.')
-
-
-class TestSingleUTube(unittest.TestCase):
-    """ Test cases for SingleUTube class.
-    """
-
-    def setUp(self):
-        self.r_b = 0.075        # Borehole radius [m]
-        self.H = 100.0          # Borehole length [m]
-        self.r_out = 0.010      # Pipe outer radius [m]
-        self.D_s = 0.060        # Shank spacing [m]
-        self.k_s = 2.0          # Ground thermal conductivity [W/m.K]
-        self.k_g = 1.0          # Grout thermal conductivity [W/m.K]
-        self.cp = 4000          # Fluid specific heat capacity [J/kg.K]
-        # Fluid to outer pipe wall thermal resistance [m.K/W]
-        self.Rfp = 1e-30
-
-    def test_outlet_fluid_temperature(self):
-        """ Tests the value of the outlet fluid temperature for a single U-tube
-            borehole.
-        """
-        from pygfunction.pipes import SingleUTube
-        from pygfunction.boreholes import Borehole
-
-        # Reference solution (Cimmino, 2016)
-        # Fluid mass flow rate [kg/s]
-        m_flow = np.array([0.2, 0.3])
-        T_b = 1.0               # Borehole wall temperature [degC]
-        Tf_in = 5.0             # Inlet fluid temperature [degC]
-        borehole = Borehole(self.H, 0., self.r_b, 0., 0.)
-
-        Tf_out_ref = np.array([2.2161, 2.8394])
-        # Calculate using pipes.SingleUTube()
-        Tf_out = np.zeros_like(Tf_out_ref)
-        for i in range(len(m_flow)):
-            pos = [(-self.D_s, 0.), (self.D_s, 0.)]
-            UTube = SingleUTube(pos, 1e-30, self.r_out, borehole,
-                                self.k_s, self.k_g, self.Rfp, J=0)
-            # Obtain outlet fluid temperature
-            Tf_out[i] = UTube.get_outlet_temperature(Tf_in, T_b,
-                                                     m_flow[i], self.cp)
-        self.assertTrue(np.allclose(Tf_out, Tf_out_ref, rtol=1e-4, atol=1e-10),
-                        msg='Incorrect value of outlet fluid temperatures.')
-
-    def test_inlet_fluid_temperature(self):
-        """ Tests the value of the inlet fluid temperature for a single U-tube
-            borehole.
-        """
-        from pygfunction.pipes import SingleUTube
-        from pygfunction.boreholes import Borehole
-
-        # Reference solution (Cimmino, 2016)
-        # Fluid mass flow rate [kg/s]
-        m_flow = np.array([0.2, 0.3])
-        T_b = 1.0               # Borehole wall temperature [degC]
-        # Outlet fluid temperature [degC]
-        Tf_out = np.array([2.2161, 2.8394])
-        # Inlet fluid temperature [degC]
-        Tf_in_ref = np.array([5.0, 5.0])
-        # Total heat transfer rate [W]
-        Qf = m_flow*self.cp*(Tf_out - Tf_in_ref)
-        borehole = Borehole(self.H, 0., self.r_b, 0., 0.)
-
-        # Calculate using pipes.SingleUTube()
-        Tf_in = np.zeros_like(Tf_in_ref)
-        for i in range(len(m_flow)):
-            pos = [(-self.D_s, 0.), (self.D_s, 0.)]
-            UTube = SingleUTube(pos, 1e-30, self.r_out, borehole,
-                                self.k_s, self.k_g, self.Rfp, J=0)
-            # Obtain inlet fluid temperature
-            Tf_in[i] = UTube.get_inlet_temperature(Qf[i], T_b,
-                                                    m_flow[i], self.cp)
-        self.assertTrue(np.allclose(Tf_in, Tf_in_ref, rtol=1e-4, atol=1e-10),
-                        msg='Incorrect value of inlet fluid temperatures.')
-
-    def test_fluid_heat_transfer_rate(self):
-        """ Tests the value of the fluid heat transfer rate for a single U-tube
-            borehole.
-        """
-        from pygfunction.pipes import SingleUTube
-        from pygfunction.boreholes import Borehole
-
-        # Reference solution (Cimmino, 2016)
-        # Fluid mass flow rate [kg/s]
-        m_flow = np.array([0.2, 0.3])
-        T_b = 1.0               # Borehole wall temperature [degC]
-        # Outlet fluid temperature [degC]
-        Tf_out = np.array([2.2161, 2.8394])
-        # Inlet fluid temperature [degC]
-        Tf_in = np.array([5.0, 5.0])
-        # Total heat transfer rate [W]
-        Qf_ref = m_flow*self.cp*(Tf_out - Tf_in)
-        borehole = Borehole(self.H, 0., self.r_b, 0., 0.)
-
-        # Calculate using pipes.SingleUTube()
-        Qf = np.zeros_like(Qf_ref)
-        for i in range(len(m_flow)):
-            pos = [(-self.D_s, 0.), (self.D_s, 0.)]
-            UTube = SingleUTube(pos, 1e-30, self.r_out, borehole,
-                                self.k_s, self.k_g, self.Rfp, J=0)
-            # Obtain fluid heat transfer rate
-            Qf[i] = UTube.get_fluid_heat_extraction_rate(Tf_in[i], T_b,
-                                                        m_flow[i], self.cp)
-        self.assertTrue(np.allclose(Qf, Qf_ref, rtol=1e-4, atol=1e-10),
-                        msg='Incorrect value of fluid heat extraction rate.')
-
-
-class TestMultipleUTube(unittest.TestCase):
-    """ Test cases for MultipleUTube class.
-    """
-
-    def setUp(self):
-        self.r_b = 0.075        # Borehole radius [m]
-        self.H = 100.0          # Borehole length [m]
-        self.r_out = 0.010      # Pipe outer radius [m]
-        self.D_s = 0.060        # Shank spacing [m]
-        self.k_s = 2.0          # Ground thermal conductivity [W/m.K]
-        self.k_g = 1.0          # Grout thermal conductivity [W/m.K]
-        self.cp = 4000          # Fluid specific heat capacity [J/kg.K]
-        # Fluid to outer pipe wall thermal resistance [m.K/W]
-        self.Rfp = 1e-30
-
-    def test_outlet_fluid_temperature_parallel(self):
-        """ Tests the value of the outlet fluid temperature for a borehole with
-            multiple U-tubes in parallel.
-        """
-        from pygfunction.pipes import MultipleUTube
-        from pygfunction.boreholes import Borehole
-
-        # Reference solution (Cimmino, 2016)
-        nPipes = [2, 3, 4]      # Number of U-tubes
-        # Fluid mass flow rate [kg/s]
-        m_flow = np.array([0.2*n for n in nPipes])
-        T_b = 1.0               # Borehole wall temperature [degC]
-        Tf_in = 5.0             # Inlet fluid temperature [degC]
-        borehole = Borehole(self.H, 0., self.r_b, 0., 0.)
-
-        Tf_out_ref = np.array([2.0860, 2.2279, 2.4627])
-        # Calculate using pipes.MultipleUTube()
-        Tf_out = np.zeros_like(Tf_out_ref)
-        for i in range(len(nPipes)):
-            pos = self._pipePositions(nPipes[i])
-            UTube = MultipleUTube(pos, 1e-30, self.r_out, borehole,
-                                  self.k_s, self.k_g, self.Rfp, nPipes[i], J=0)
-            # Obtain outlet fluid temperature
-            Tf_out[i] = UTube.get_outlet_temperature(Tf_in, T_b,
-                                                     m_flow[i], self.cp)
-        self.assertTrue(np.allclose(Tf_out, Tf_out_ref, rtol=1e-4, atol=1e-10),
-                        msg='Incorrect value of outlet fluid temperatures for '
-                            'U-tubes in parallel.')
-
-    def test_outlet_fluid_temperature_series(self):
-        """ Tests the value of the outlet fluid temperature for a borehole with
-            multiple U-tubes in series.
-        """
-        from pygfunction.pipes import MultipleUTube
-        from pygfunction.boreholes import Borehole
-
-        # Reference solution (Cimmino, 2016)
-        nPipes = [2, 3, 4]      # Number of U-tubes
-        # Fluid mass flow rate [kg/s]
-        m_flow = 0.3
-        T_b = 1.0               # Borehole wall temperature [degC]
-        Tf_in = 5.0             # Inlet fluid temperature [degC]
-        borehole = Borehole(self.H, 0., self.r_b, 0., 0.)
-
-        Tf_out_ref = np.array([1.6215, 1.2051, 1.1065])
-        # Calculate using pipes.MultipleUTube()
-        Tf_out = np.zeros_like(Tf_out_ref)
-        for i in range(len(nPipes)):
-            pos = self._pipePositions(nPipes[i])
-            UTube = MultipleUTube(pos, 1e-30, self.r_out, borehole,
-                                  self.k_s, self.k_g, self.Rfp, nPipes[i],
-                                  config='series', J=0)
-            # Obtain outlet fluid temperature
-            Tf_out[i] = UTube.get_outlet_temperature(Tf_in, T_b,
-                                                     m_flow, self.cp)
-        self.assertTrue(np.allclose(Tf_out, Tf_out_ref, rtol=1e-4, atol=1e-10),
-                        msg='Incorrect value of outlet fluid temperatures for '
-                            'U-tubes in series.')
-
-    def test_inlet_fluid_temperature_parallel(self):
-        """ Tests the value of the inlet fluid temperature for a borehole with
-            multiple U-tubes in parallel.
-        """
-        from pygfunction.pipes import MultipleUTube
-        from pygfunction.boreholes import Borehole
-
-        # Reference solution (Cimmino, 2016)
-        nPipes = [2, 3, 4]      # Number of U-tubes
-        # Fluid mass flow rate [kg/s]
-        m_flow = np.array([0.2*n for n in nPipes])
-        T_b = 1.0               # Borehole wall temperature [degC]
-        # Outlet fluid temperature [degC]
-        Tf_out = np.array([2.0860, 2.2279, 2.4627])
-        # Inlet fluid temperature [degC]
-        Tf_in_ref = np.array([5.0, 5.0, 5.0])
-        # Total heat transfer rate [W]
-        Qf = m_flow*self.cp*(Tf_out - Tf_in_ref)
-        borehole = Borehole(self.H, 0., self.r_b, 0., 0.)
-
-        # Calculate using pipes.MultipleUTube()
-        Tf_in = np.zeros_like(Tf_in_ref)
-        for i in range(len(nPipes)):
-            pos = self._pipePositions(nPipes[i])
-            UTube = MultipleUTube(pos, 1e-30, self.r_out, borehole,
-                                  self.k_s, self.k_g, self.Rfp, nPipes[i], J=0)
-            # Obtain inlet fluid temperature
-            Tf_in[i] = UTube.get_inlet_temperature(Qf[i], T_b,
-                                                   m_flow[i], self.cp)
-        self.assertTrue(np.allclose(Tf_in, Tf_in_ref, rtol=1e-4, atol=1e-10),
-                        msg='Incorrect value of inlet fluid temperatures for '
-                            'U-tubes in parallel.')
-
-    def test_inlet_fluid_temperature_series(self):
-        """ Tests the value of the inlet fluid temperature for a borehole with
-            multiple U-tubes in series.
-        """
-        from pygfunction.pipes import MultipleUTube
-        from pygfunction.boreholes import Borehole
-
-        # Reference solution (Cimmino, 2016)
-        nPipes = [2, 3, 4]      # Number of U-tubes
-        # Fluid mass flow rate [kg/s]
-        m_flow = 0.3
-        T_b = 1.0               # Borehole wall temperature [degC]
-        # Outlet fluid temperature [degC]
-        Tf_out = np.array([1.6215, 1.2051, 1.1065])
-        # Inlet fluid temperature [degC]
-        Tf_in_ref = np.array([5.0, 5.0, 5.0])
-        # Total heat transfer rate [W]
-        Qf = m_flow*self.cp*(Tf_out - Tf_in_ref)
-        borehole = Borehole(self.H, 0., self.r_b, 0., 0.)
-
-        # Calculate using pipes.MultipleUTube()
-        Tf_in = np.zeros_like(Tf_in_ref)
-        for i in range(len(nPipes)):
-            pos = self._pipePositions(nPipes[i])
-            UTube = MultipleUTube(pos, 1e-30, self.r_out, borehole,
-                                  self.k_s, self.k_g, self.Rfp, nPipes[i],
-                                  config='series', J=0)
-            # Obtain inlet fluid temperature
-            Tf_in[i] = UTube.get_inlet_temperature(Qf[i], T_b,
-                                                   m_flow, self.cp)
-        self.assertTrue(np.allclose(Tf_in, Tf_in_ref, rtol=1e-4, atol=1e-10),
-                        msg='Incorrect value of inlet fluid temperatures for '
-                            'U-tubes in series.')
-
-    def test_fluid_heat_transfer_rate_parallel(self):
-        """ Tests the value of the fluid heat transfer rate for a borehole with
-            multiple U-tubes in parallel.
-        """
-        from pygfunction.pipes import MultipleUTube
-        from pygfunction.boreholes import Borehole
-
-        # Reference solution (Cimmino, 2016)
-        nPipes = [2, 3, 4]      # Number of U-tubes
-        # Fluid mass flow rate [kg/s]
-        m_flow = np.array([0.2*n for n in nPipes])
-        T_b = 1.0               # Borehole wall temperature [degC]
-        # Outlet fluid temperature [degC]
-        Tf_out = np.array([2.0860, 2.2279, 2.4627])
-        # Inlet fluid temperature [degC]
-        Tf_in = np.array([5.0, 5.0, 5.0])
-        # Total heat transfer rate [W]
-        Qf_ref = m_flow*self.cp*(Tf_out - Tf_in)
-        borehole = Borehole(self.H, 0., self.r_b, 0., 0.)
-
-        # Calculate using pipes.MultipleUTube()
-        Qf = np.zeros_like(Qf_ref)
-        for i in range(len(nPipes)):
-            pos = self._pipePositions(nPipes[i])
-            UTube = MultipleUTube(pos, 1e-30, self.r_out, borehole,
-                                  self.k_s, self.k_g, self.Rfp, nPipes[i], J=0)
-            # Obtain fluid heat extraction rate
-            Qf[i] = UTube.get_fluid_heat_extraction_rate(Tf_in[i], T_b,
-                                                         m_flow[i], self.cp)
-        self.assertTrue(np.allclose(Qf, Qf_ref, rtol=1e-4, atol=1e-10),
-                        msg='Incorrect value of fluid heat extraction rate '
-                            'for U-tubes in parallel.')
-
-    def test_fluid_heat_transfer_rate_series(self):
-        """ Tests the value of the fluid heat transfer rate for a borehole with
-            multiple U-tubes in series.
-        """
-        from pygfunction.pipes import MultipleUTube
-        from pygfunction.boreholes import Borehole
-
-        # Reference solution (Cimmino, 2016)
-        nPipes = [2, 3, 4]      # Number of U-tubes
-        # Fluid mass flow rate [kg/s]
-        m_flow = 0.3
-        T_b = 1.0               # Borehole wall temperature [degC]
-        # Outlet fluid temperature [degC]
-        Tf_out = np.array([1.6215, 1.2051, 1.1065])
-        # Inlet fluid temperature [degC]
-        Tf_in = np.array([5.0, 5.0, 5.0])
-        # Total heat transfer rate [W]
-        Qf_ref = m_flow*self.cp*(Tf_out - Tf_in)
-        borehole = Borehole(self.H, 0., self.r_b, 0., 0.)
-
-        # Calculate using pipes.MultipleUTube()
-        Qf = np.zeros_like(Qf_ref)
-        for i in range(len(nPipes)):
-            pos = self._pipePositions(nPipes[i])
-            UTube = MultipleUTube(pos, 1e-30, self.r_out, borehole,
-                                  self.k_s, self.k_g, self.Rfp, nPipes[i],
-                                  config='series', J=0)
-            # Obtain fluid heat extraction rate
-            Qf[i] = UTube.get_fluid_heat_extraction_rate(Tf_in[i], T_b,
-                                                         m_flow, self.cp)
-        self.assertTrue(np.allclose(Qf, Qf_ref, rtol=1e-4, atol=1e-10),
-                        msg='Incorrect value of fluid heat extraction rate '
-                            'for U-tubes in series.')
-
-    def _pipePositions(self, nPipes):
-        """ Positions pipes in an axisymetric configuration.
-        """
-        dt = pi / float(nPipes)
-        pos = [(0., 0.) for i in range(2*nPipes)]
-        for i in range(nPipes):
-            pos[i] = (self.D_s*np.cos(2.0*i*dt+pi),
-                      self.D_s*np.sin(2.0*i*dt+pi))
-            pos[i+nPipes] = (self.D_s*np.cos(2.0*i*dt+pi+dt),
-                             self.D_s*np.sin(2.0*i*dt+pi+dt))
-        return pos
-
-
-class TestIndependentMultipleUTube(unittest.TestCase):
-    """ Test cases for MultipleUTube class.
-    """
-
-    def setUp(self):
-        self.r_b = 0.075        # Borehole radius [m]
-        self.H = 100.0          # Borehole length [m]
-        self.T_b = 2.0          # Borehole wall temperature [degC]
-        self.r_out = 0.010      # Pipe outer radius [m]
-        self.D_s = 0.060        # Shank spacing [m]
-        self.nPipes = 4         # Number of U-tubes
-        self.k_s = 2.0          # Ground thermal conductivity [W/m.K]
-        self.k_g = 1.0          # Grout thermal conductivity [W/m.K]
-        # Fluid specific heat capacity [J/kg.K]
-        self.cp = 4000
-        # Fluid mass flow rate in each U-tube [kg/s]
-        self.m_flow = np.array([0.40, 0.35, 0.30, 0.25])
-        # Fluid specific heat capacity [J/kg.K]
-        self.cp = 4000.*np.ones_like(self.m_flow)
-        # Fluid to outer pipe wall thermal resistance [m.K/W]
-        self.Rfp = 1e-30
-        # Inlet fluid temperature in each U-tube [degC]
-        self.Tf_in_ref = np.array([6.0, -6.0, 5.0, -5.0])
-        # outlet fluid temperature in each U-tube [degC]
-        self.Tf_out_ref = np.array([4.4434, -2.5261, 3.3923, -0.89299])
-        # Total heat transfer rate [W]
-        self.Qf_ref = self.m_flow*self.cp*(self.Tf_out_ref - self.Tf_in_ref)
-
-    def test_outlet_fluid_temperature(self):
-        """ Tests the value of the outlet fluid temperature for a borehole with
-            multiple independent U-tubes.
-        """
-        from pygfunction.pipes import IndependentMultipleUTube
-        from pygfunction.boreholes import Borehole
-
-        borehole = Borehole(self.H, 0., self.r_b, 0., 0.)
-
-        # Calculate using pipes.IndependentMultipleUTube()
-        pos = self._pipePositions(self.nPipes)
-        UTube = IndependentMultipleUTube(pos, 1e-30, self.r_out, borehole,
-                                         self.k_s, self.k_g, self.Rfp,
-                                         self.nPipes, J=0)
-        # Obtain outlet fluid temperature
-        Tf_out = UTube.get_outlet_temperature(self.Tf_in_ref, self.T_b,
-                                              self.m_flow, self.cp)
-        self.assertTrue(np.allclose(Tf_out, self.Tf_out_ref,
-                                    rtol=1e-4, atol=1e-10),
-                        msg='Incorrect value of outlet fluid temperatures for '
-                            'independent U-tubes.')
-
-    def test_inlet_fluid_temperature(self):
-        """ Tests the value of the inlet fluid temperature for a borehole with
-            multiple independent U-tubes.
-        """
-        from pygfunction.pipes import IndependentMultipleUTube
-        from pygfunction.boreholes import Borehole
-
-        borehole = Borehole(self.H, 0., self.r_b, 0., 0.)
-
-        # Calculate using pipes.IndependentMultipleUTube()
-        pos = self._pipePositions(self.nPipes)
-        UTube = IndependentMultipleUTube(pos, 1e-30, self.r_out, borehole,
-                                         self.k_s, self.k_g, self.Rfp,
-                                         self.nPipes, J=0)
-        # Obtain inlet fluid temperature
-        Tf_in = UTube.get_inlet_temperature(self.Qf_ref, self.T_b,
-                                            self.m_flow, self.cp)
-        self.assertTrue(np.allclose(Tf_in, self.Tf_in_ref,
-                                    rtol=1e-4, atol=1e-10),
-                        msg='Incorrect value of inlet fluid temperatures for '
-                            'independent U-tubes.')
-
-    def test_fluid_heat_transfer_rate(self):
-        """ Tests the value of the fluid heat transfer rates for a borehole
-            with multiple independent U-tubes.
-        """
-        from pygfunction.pipes import IndependentMultipleUTube
-        from pygfunction.boreholes import Borehole
-
-        borehole = Borehole(self.H, 0., self.r_b, 0., 0.)
-
-        # Calculate using pipes.IndependentMultipleUTube()
-        pos = self._pipePositions(self.nPipes)
-        UTube = IndependentMultipleUTube(pos, 1e-30, self.r_out, borehole,
-                                         self.k_s, self.k_g, self.Rfp,
-                                         self.nPipes, J=0)
-        # Obtain outlet fluid temperature
-        Qf = UTube.get_fluid_heat_extraction_rate(self.Tf_in_ref, self.T_b,
-                                                  self.m_flow, self.cp)
-        self.assertTrue(np.allclose(Qf, self.Qf_ref,
-                                    rtol=1e-4, atol=1e-10),
-                        msg='Incorrect value of fluid heat extraction rates '
-                            'for independent U-tubes.')
-
-    def _pipePositions(self, nPipes):
-        """ Positions pipes in an axisymetric configuration.
-        """
-        dt = pi / float(nPipes)
-        pos = [(0., 0.) for i in range(2*nPipes)]
-        for i in range(nPipes):
-            pos[i] = (self.D_s*np.cos(2.0*i*dt+pi),
-                      self.D_s*np.sin(2.0*i*dt+pi))
-            pos[i+nPipes] = (self.D_s*np.cos(2.0*i*dt+pi+dt),
-                             self.D_s*np.sin(2.0*i*dt+pi+dt))
-        return pos
-
-
-if __name__ == '__main__' and __package__ is None:
-    from os import sys, path
-    sys.path.append(path.dirname(path.dirname(path.abspath(__file__))))
-    unittest.main()
+
+import pygfunction as gt
+
+
+# =============================================================================
+# Test functions
+# =============================================================================
+# Test convective_heat_transfer_coefficient_circular_pipe
+@pytest.mark.parametrize("m_flow, expected", [
+    (0.05, 90.07260000000001),  # Laminar flow
+    (0.10, 572.256944167273),   # Transition flow
+    (0.50, 4036.196217814895)   # Turbulent flow
+    ])
+def test_convective_heat_transfer_coefficient_circular_pipe(m_flow, expected):
+    r_in = 0.01         # Inner radius [m]
+    epsilon = 1.5e-6    # Pipe surface roughness [m]
+    visc = 0.00203008   # Fluid viscosity [kg/m.s]
+    den = 1014.78       # Fluid density [kg/m3]
+    cp = 3977.          # Fluid specific heat capacity [J/kg.K]
+    k = 0.4922          # Fluid thermal conductivity [W/m.K]
+    # Convective heat transfer coefficient [W/m2.K]
+    h = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
+        m_flow, r_in, visc, den, k, cp, epsilon)
+    assert np.isclose(h, expected)
+
+
+# Test convective_heat_transfer_coefficient_concentric_annulus
+@pytest.mark.parametrize("m_flow, expected", [
+    (0.05, (141.4907984705223, 110.95487746200112)),    # Laminar flow
+    (0.40, (904.4869811625874, 904.4869811625874)),     # Transition flow
+    (0.60, (1411.2063074288633, 1411.2063074288633))    # Turbulent flow
+    ])
+def test_convective_heat_transfer_coefficient_concentric_annulus(
+        m_flow, expected):
+    r_in = 0.01         # Inner radius [m]
+    r_out = 0.02        # Outer radius [m]
+    epsilon = 1.5e-6    # Pipe surface roughness [m]
+    visc = 0.00203008   # Fluid viscosity [kg/m.s]
+    den = 1014.78       # Fluid density [kg/m3]
+    cp = 3977.          # Fluid specific heat capacity [J/kg.K]
+    k = 0.4922          # Fluid thermal conductivity [W/m.K]
+    # Convective heat transfer coefficients [W/m2.K]
+    h = gt.pipes.convective_heat_transfer_coefficient_concentric_annulus(
+        m_flow, r_in, r_out, visc, den, k, cp, epsilon)
+    assert np.allclose(h, expected)
+
+
+# Test conduction_thermal_resistance_circular_pipe
+def test_conduction_thermal_resistance_circular_pipe():
+    r_in = 0.01     # Inner radius [m]
+    r_out = 0.02    # Outer radius [m]
+    k = 0.6         # Fluid thermal conductivity [W/m.K]
+    # Conduction thermal resistance [m.K/W]
+    expected = 0.18386300012720966
+    R = gt.pipes.conduction_thermal_resistance_circular_pipe(r_in, r_out, k)
+    assert np.isclose(R, expected)
+
+
+# Test fluid_friction_factor_circular_pipe
+@pytest.mark.parametrize("m_flow, expected", [
+    (0.05, 0.04081718025087723),    # Laminar flow
+    (0.50, 0.027641340780182006)    # Turbulent flow
+    ])
+def test_fluid_friction_factor_circular_pipe(m_flow, expected):
+    r_in = 0.01         # Inner radius [m]
+    epsilon = 1.5e-6    # Pipe surface roughness [m]
+    visc = 0.00203008   # Fluid viscosity [kg/m.s]
+    den = 1014.78       # Fluid density [kg/m3]
+    # Fluid Darcy friction factor [-]
+    f = gt.pipes.fluid_friction_factor_circular_pipe(
+        m_flow, r_in, visc, den, epsilon)
+    assert np.isclose(f, expected)
+
+
+# Test thermal_resistances
+@pytest.mark.parametrize("J, R_expected, Rd_expected", [
+    (0,     # Zero-th order multipoles
+     np.array([[0.25486306, 0.01538038],
+               [0.01538038, 0.25206829]]),
+     np.array([[0.27042505, 4.16155713],
+               [4.16155713, 0.26726918]])),
+    (1,     # First order multipoles
+     np.array([[0.25569372, 0.01562313],
+               [0.01562313, 0.25288076]]),
+     np.array([[0.27150208, 4.12311447],
+               [4.12311447, 0.26832082]])),
+    (2,     # Second order multipoles
+     np.array([[0.25590404, 0.01560503],
+               [0.01560503, 0.25308681]]),
+     np.array([[0.27169419, 4.13472001],
+               [4.13472001, 0.26850888]])),
+    (3,     # Third order multipoles
+     np.array([[0.25592405, 0.01560826],
+               [0.01560826, 0.25310667]]),
+     np.array([[0.27171747, 4.1345064 ],
+               [4.1345064, 0.26853194]])),
+    ])
+def test_thermal_resistances(J, R_expected, Rd_expected):
+    # Pipe positions [m]
+    pos_2pipes = [(0.03, 0.00), (-0.03, 0.02)]
+    r_out = 0.02    # Pipe outer radius [m]
+    r_b = 0.07      # Borehole radius [m]
+    k_s = 2.5       # Ground thermal conductivity [W/m.K]
+    k_g = 1.5       # Grout thermal conductivity [W/m.K]
+    # Fluid to outer pipe wall thermal resistance [m.K/W]
+    beta = 1.2
+    Rfp = beta / (2 * np.pi * k_g)
+    # Thermal resistances [m.k/W]
+    R, Rd = gt.pipes.thermal_resistances(
+        pos_2pipes, r_out, r_b, k_s, k_g, Rfp, J=J)
+    assert np.allclose(R, R_expected) and np.allclose(Rd, Rd_expected)
+
+
+# Test multipole
+@pytest.mark.parametrize("J, expected", [
+    (0, np.array([2.67263436, 2.47271955, 2.15219567, 1.95228086])),
+    (1, np.array([2.71748588, 2.51729508, 2.19369127, 1.99350047])),
+    (2, np.array([2.71914947, 2.51894712, 2.19312216, 1.99291981])),
+    (3, np.array([2.71942944, 2.51913631, 2.19328373, 1.9929906])),
+    ])
+def test_multipole(J, expected):
+    # Pipe positions [m]
+    pos_4pipes = [
+        (0.03, 0.03), (-0.03, 0.03), (-0.03, -0.03), (0.03, -0.03)]
+    r_out = 0.02    # Pipe outer radius [m]
+    r_b = 0.07      # Borehole radius [m]
+    k_s = 2.5       # Ground thermal conductivity [W/m.K]
+    k_g = 1.5       # Grout thermal conductivity [W/m.K]
+    T_b = 0.0       # Borehole wall temperature [degC]
+    # Pipe heat transfer rates [W/m]
+    q_p = np.array([10., 9., 8., 7.])
+    # Fluid to outer pipe wall thermal resistance [m.K/W]
+    beta = 1.2
+    R_fp = beta / (2 * np.pi * k_g)
+    # Fluid temperatures [degC]
+    T_f = gt.pipes.multipole(
+        pos_4pipes, r_out, r_b, k_s, k_g, R_fp, T_b, q_p, J)[0]
+    assert np.allclose(T_f, expected)
+
+
+# =============================================================================
+# Test pipe classes
+# =============================================================================
+# Test get_temperature
+@pytest.mark.parametrize("pipe_fixture, segment_ratios, T_b, z, expected", [
+    # Single U-tube
+    ('single_Utube', None, 1., 65., np.array([4.34676755, 3.07354134])),
+    ('single_Utube', None, np.array([1., 2., 3., 1.]), 65., np.array([4.41754093, 3.49949295])),
+    ('single_Utube', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([4.47310416, 3.66490249])),
+    ('single_Utube', None, 1., np.array([65., 75.]), np.array([[4.34676755, 3.07354134], [4.25566624, 3.13435325]])),
+    ('single_Utube', None, np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.41754093, 3.49949295], [4.35173147, 3.54346564]])),
+    ('single_Utube', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.47310416, 3.66490249], [4.42647994, 3.69214797]])),
+    # Double U-tube (Parallel)
+    ('double_Utube_parallel', None, 1., 65., np.array([3.87525104, 3.87525104, 2.20313908, 2.20313908])),
+    ('double_Utube_parallel', None, np.array([1., 2., 3., 1.]), 65., np.array([4.00464852, 4.00464852, 2.84788608, 2.84788608])),
+    ('double_Utube_parallel', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([4.09862044, 4.09862044, 3.07439258, 3.07439258])),
+    ('double_Utube_parallel', None, 1., np.array([65., 75.]), np.array([[3.87525104, 3.87525104, 2.20313908, 2.20313908], [3.73265141, 3.73265141, 2.26719823, 2.26719823]])),
+    ('double_Utube_parallel', None, np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.00464852, 4.00464852, 2.84788608, 2.84788608], [3.90522192, 3.90522192, 2.89301847, 2.89301847]])),
+    ('double_Utube_parallel', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.09862044, 4.09862044, 3.07439258, 3.07439258], [4.03210009, 4.03210009, 3.09222735, 3.09222735]])),
+    # Double U-tube (Series)
+    ('double_Utube_series', None, 1., 65., np.array([4.36908096, 2.53231146, 3.13441957, 2.03763963])),
+    ('double_Utube_series', None, np.array([1., 2., 3., 1.]), 65., np.array([4.44022419, 2.94528677, 3.54323578, 2.65057213])),
+    ('double_Utube_series', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([4.49394782, 3.15796672, 3.7037625, 2.90601555])),
+    ('double_Utube_series', None, 1., np.array([65., 75.]), np.array([[4.36908096, 2.53231146, 3.13441957, 2.03763963], [4.28094228, 2.49706752, 3.19348974, 2.0612353]])),
+    ('double_Utube_series', None, np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.44022419, 2.94528677, 3.54323578, 2.65057213], [4.37608674, 2.92420008, 3.58625745, 2.66472128]])),
+    ('double_Utube_series', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.49394782, 3.15796672, 3.7037625, 2.90601555], [4.44772426, 3.15021634, 3.73124558, 2.90769612]])),
+    # Coaxial (Annular pipe is inlet pipe)
+    ('coaxial_annular_in', None, 1., 65., np.array([3.15203088, 2.18408362])),
+    ('coaxial_annular_in', None, np.array([1., 2., 3., 1.]), 65., np.array([3.4176666 , 2.73205968])),
+    ('coaxial_annular_in', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([3.57428289, 3.01759194])),
+    ('coaxial_annular_in', None, 1., np.array([65., 75.]), np.array([[3.15203088, 2.18408362], [2.96401382, 2.15051705]])),
+    ('coaxial_annular_in', None, np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[3.4176666, 2.73205968], [3.2920645, 2.7081367]])),
+    ('coaxial_annular_in', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[3.57428289, 3.01759194], [3.4962546, 2.99796987]])),
+    # Coaxial (Annular pipe is outlet pipe)
+    ('coaxial_annular_out', None, 1., 65., np.array([4.50649998, 2.92933532])),
+    ('coaxial_annular_out', None, np.array([1., 2., 3., 1.]), 65., np.array([4.62416027, 3.50307539])),
+    ('coaxial_annular_out', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([4.64884211, 3.58227143])),
+    ('coaxial_annular_out', None, 1., np.array([65., 75.]), np.array([[4.50649998, 2.92933532], [4.44976224, 3.02086677]])),
+    ('coaxial_annular_out', None, np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.62416027, 3.50307539], [4.58402116, 3.57860389]])),
+    ('coaxial_annular_out', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.64884211, 3.58227143], [4.61006377, 3.60745651]])),
+    ])
+def test_temperature(
+        pipe_fixture, segment_ratios, T_b, z, expected, request):
+    # Extract pipe from fixture
+    pipe = request.getfixturevalue(pipe_fixture)
+    # Fluid is propylene-glycol 20%
+    fluid = gt.media.Fluid('MPG', 20.)
+    m_flow_borehole = 0.2   # Fluid mass flow rate [kg/s]
+    T_f_in = 5.0            # Inlet fluid temperature [degC]
+    # Fluid temperatures [degC]
+    T_f = pipe.get_temperature(
+        z, T_f_in, T_b, m_flow_borehole, fluid.cp,
+        segment_ratios=segment_ratios)
+    assert np.allclose(T_f, expected)
+
+
+# Test get_outlet_temperature
+@pytest.mark.parametrize("pipe_fixture, segment_ratios, T_b, expected", [
+    # Single U-tube
+    ('single_Utube', None, 1., 2.712371852688313),
+    ('single_Utube', None, np.array([1., 2., 3., 1.]), 3.1377635748663573),
+    ('single_Utube', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 3.335501080169065),
+    # Double U-tube (Parallel)
+    ('double_Utube_parallel', None, 1., 1.8553031331306218),
+    ('double_Utube_parallel', None, np.array([1., 2., 3., 1.]), 2.4278457017624655),
+    ('double_Utube_parallel', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 2.691532668379643),
+    # Double U-tube (Series)
+    ('double_Utube_series', None, 1., 1.8983711735742064),
+    ('double_Utube_series', None, np.array([1., 2., 3., 1.]), 2.4755999700741573),
+    ('double_Utube_series', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 2.744018786582172),
+    # Coaxial (Annular pipe is inlet pipe)
+    ('coaxial_annular_in', None, 1., 2.581130521333567),
+    ('coaxial_annular_in', None, np.array([1., 2., 3., 1.]), 3.0276625795763357),
+    ('coaxial_annular_in', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 3.2733787998105672),
+    # Coaxial (Annular pipe is outlet pipe)
+    ('coaxial_annular_out', None, 1., 2.5811305213335674),
+    ('coaxial_annular_out', None, np.array([1., 2., 3., 1.]), 2.981638747649938),
+    ('coaxial_annular_out', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 3.148206616090016),
+    ])
+def test_outlet_temperature(
+        pipe_fixture, segment_ratios, T_b, expected, request):
+    # Extract pipe from fixture
+    pipe = request.getfixturevalue(pipe_fixture)
+    # Fluid is propylene-glycol 20%
+    fluid = gt.media.Fluid('MPG', 20.)
+    m_flow_borehole = 0.2   # Fluid mass flow rate [kg/s]
+    T_f_in = 5.0            # Inlet fluid temperature [degC]
+    # Outlet fluid temperature [degC]
+    T_f_out = pipe.get_outlet_temperature(
+        T_f_in, T_b, m_flow_borehole, fluid.cp, segment_ratios=segment_ratios)
+    assert np.isclose(T_f_out, expected)
+
+
+# Test get_inlet_temperature
+@pytest.mark.parametrize("pipe_fixture, segment_ratios, T_b, expected", [
+    # Single U-tube
+    ('single_Utube', None, 1., 7.595314034714041),
+    ('single_Utube', None, np.array([1., 2., 3., 1.]), 8.33912674339739),
+    ('single_Utube', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 8.68487787525871),
+    # Double U-tube (Parallel)
+    ('double_Utube_parallel', None, 1., 5.7977998086638305),
+    ('double_Utube_parallel', None, np.array([1., 2., 3., 1.]), 6.526064048901171),
+    ('double_Utube_parallel', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 6.861469307697622),
+    # Double U-tube (Series)
+    ('double_Utube_series', None, 1., 5.8644202354664365),
+    ('double_Utube_series', None, np.array([1., 2., 3., 1.]), 6.60884044665609),
+    ('double_Utube_series', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 6.955005421937585),
+    # Coaxial (Annular pipe is inlet pipe)
+    ('coaxial_annular_in', None, 1., 7.237470090568812),
+    ('coaxial_annular_in', None, np.array([1., 2., 3., 1.]), 7.97588456424095),
+    ('coaxial_annular_in', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 8.382216898252953),
+    # Coaxial (Annular pipe is outlet pipe)
+    ('coaxial_annular_out', None, 1., 7.237470090568813),
+    ('coaxial_annular_out', None, np.array([1., 2., 3., 1.]), 7.899776560345228),
+    ('coaxial_annular_out', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 8.175224028526785),
+    ])
+def test_inlet_temperature(pipe_fixture, segment_ratios, T_b, expected, request):
+    # Extract pipe from fixture
+    pipe = request.getfixturevalue(pipe_fixture)
+    # Fluid is propylene-glycol 20%
+    fluid = gt.media.Fluid('MPG', 20.)
+    m_flow_borehole = 0.2   # Fluid mass flow rate [kg/s]
+    Q_f = -3000.0           # Total heat transfer rate [W]
+    # Inlet fluid temperature [degC]
+    T_f_in = pipe.get_inlet_temperature(
+        Q_f, T_b, m_flow_borehole, fluid.cp, segment_ratios=segment_ratios)
+    assert np.isclose(T_f_in, expected)
+
+
+# Test get_borehole_heat_extraction_rate
+@pytest.mark.parametrize("pipe_fixture, segment_ratios, T_b, expected", [
+    # Single U-tube
+    ('single_Utube', None, 1., -1819.4736348927008),
+    ('single_Utube', None, np.array([1., 2., 3., 1.]), np.array([-507.98022943, -330.29924271, -155.92399643, -486.93326314])),
+    ('single_Utube', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-211.44166366, -492.72964222, -279.36944922, -303.45149861])),
+    # Double U-tube (Parallel)
+    ('double_Utube_parallel', None, 1., -2501.14645849),
+    ('double_Utube_parallel', None, np.array([1., 2., 3., 1.]), np.array([-796.48662356, -444.22614316, -108.02227066, -697.03753979])),
+    ('double_Utube_parallel', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-342.52693756, -711.1206948, -251.40681559, -447.57582464])),
+    # Double U-tube (Series)
+    ('double_Utube_series', None, 1., -2466.89213085),
+    ('double_Utube_series', None, np.array([1., 2., 3., 1.]), np.array([-745.16518357, -428.05472293, -114.8035859, -719.7675482])),
+    ('double_Utube_series', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-319.92216781, -677.35565178, -267.23882139, -463.25689612])),
+    # Coaxial (Annular pipe is inlet pipe)
+    ('coaxial_annular_in', None, 1., -1923.85692048),
+    ('coaxial_annular_in', None, np.array([1., 2., 3., 1.]), np.array([-757.51176437, -346.76503548, -48.92829119, -415.50088061])),
+    ('coaxial_annular_in', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-335.11374414, -618.64814833, -139.77685268, -279.11910948])),
+    # Coaxial (Annular pipe is outlet pipe)
+    ('coaxial_annular_out', None, 1., -1923.85692048),
+    ('coaxial_annular_out', None, np.array([1., 2., 3., 1.]), np.array([-480.81667849, -324.83211948, -133.10520419, -666.55719699])),
+    ('coaxial_annular_out', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-199.09314899, -449.21544246, -227.76223639, -389.6723144])),
+    ])
+def test_borehole_heat_extraction_rate(
+        pipe_fixture, segment_ratios, T_b, expected, request):
+    # Extract pipe from fixture
+    pipe = request.getfixturevalue(pipe_fixture)
+    # Fluid is propylene-glycol 20%
+    fluid = gt.media.Fluid('MPG', 20.)
+    m_flow_borehole = 0.2   # Fluid mass flow rate [kg/s]
+    T_f_in = 5.0            # Inlet fluid temperature [degC]
+    # Borehole heat extraction rates [W]
+    Q_b = pipe.get_borehole_heat_extraction_rate(
+        T_f_in, T_b, m_flow_borehole, fluid.cp, segment_ratios=segment_ratios)
+    assert np.allclose(Q_b, expected)
+
+
+# Test get_fluid_heat_extraction_rate
+@pytest.mark.parametrize("pipe_fixture, segment_ratios, T_b, expected", [
+    # Single U-tube
+    ('single_Utube', None, 1., -1819.4736348927008),
+    ('single_Utube', None, np.array([1., 2., 3., 1.]), -1481.1367317058312),
+    ('single_Utube', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1323.8654645418355),
+    # Double U-tube (Parallel)
+    ('double_Utube_parallel', None, 1., -2501.14645849),
+    ('double_Utube_parallel', None, np.array([1., 2., 3., 1.]), -2045.7725771641726),
+    ('double_Utube_parallel', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1836.0481583644691),
+    # Double U-tube (Series)
+    ('double_Utube_series', None, 1., -2466.89213085),
+    ('double_Utube_series', None, np.array([1., 2., 3., 1.]), -2007.7910405893485),
+    ('double_Utube_series', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1794.303127215247),
+    # Coaxial (Annular pipe is inlet pipe)
+    ('coaxial_annular_in', None, 1., -1923.85692048),
+    ('coaxial_annular_in', None, np.array([1., 2., 3., 1.]), -1568.705971637178),
+    ('coaxial_annular_in', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1373.2746534366859),
+    # Coaxial (Annular pipe is outlet pipe)
+    ('coaxial_annular_out', None, 1., -1923.85692048),
+    ('coaxial_annular_out', None, np.array([1., 2., 3., 1.]), -1605.3111991367698),
+    ('coaxial_annular_out', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1472.830819664631),
+    ])
+def test_fluid_heat_extraction_rate(
+        pipe_fixture, segment_ratios, T_b, expected, request):
+    # Extract pipe from fixture
+    pipe = request.getfixturevalue(pipe_fixture)
+    # Fluid is propylene-glycol 20%
+    fluid = gt.media.Fluid('MPG', 20.)
+    m_flow_borehole = 0.2   # Fluid mass flow rate [kg/s]
+    T_f_in = 5.0            # Inlet fluid temperature [degC]
+    # Fluid heat extraction rate [W]
+    Q_f = pipe.get_fluid_heat_extraction_rate(
+        T_f_in, T_b, m_flow_borehole, fluid.cp, segment_ratios=segment_ratios)
+    assert np.isclose(Q_f, expected)
+
+
+# Test get_total_heat_extraction_rate
+@pytest.mark.parametrize("pipe_fixture, segment_ratios, T_b, expected", [
+    # Single U-tube
+    ('single_Utube', None, 1., -1819.4736348927008),
+    ('single_Utube', None, np.array([1., 2., 3., 1.]), -1481.1367317058312),
+    ('single_Utube', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1323.8654645418355),
+    # Double U-tube (Parallel)
+    ('double_Utube_parallel', None, 1., -2501.14645849),
+    ('double_Utube_parallel', None, np.array([1., 2., 3., 1.]), -2045.7725771641726),
+    ('double_Utube_parallel', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1836.0481583644691),
+    # Double U-tube (Series)
+    ('double_Utube_series', None, 1., -2466.89213085),
+    ('double_Utube_series', None, np.array([1., 2., 3., 1.]), -2007.7910405893485),
+    ('double_Utube_series', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1794.303127215247),
+    # Coaxial (Annular pipe is inlet pipe)
+    ('coaxial_annular_in', None, 1., -1923.85692048),
+    ('coaxial_annular_in', None, np.array([1., 2., 3., 1.]), -1568.705971637178),
+    ('coaxial_annular_in', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1373.2746534366859),
+    # Coaxial (Annular pipe is outlet pipe)
+    ('coaxial_annular_out', None, 1., -1923.85692048),
+    ('coaxial_annular_out', None, np.array([1., 2., 3., 1.]), -1605.3111991367698),
+    ('coaxial_annular_out', np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1472.830819664631),
+    ])
+def test_total_heat_extraction_rate(
+        pipe_fixture, segment_ratios, T_b, expected, request):
+    # Extract pipe from fixture
+    pipe = request.getfixturevalue(pipe_fixture)
+    # Fluid is propylene-glycol 20%
+    fluid = gt.media.Fluid('MPG', 20.)
+    m_flow_borehole = 0.2   # Fluid mass flow rate [kg/s]
+    T_f_in = 5.0            # Inlet fluid temperature [degC]
+    # Total heat extraction rate [W]
+    Q_t = pipe.get_total_heat_extraction_rate(
+        T_f_in, T_b, m_flow_borehole, fluid.cp, segment_ratios=segment_ratios)
+    assert np.isclose(Q_t, expected)
+
+
+# =============================================================================
+# Test IndependentMultipleUTube class
+# =============================================================================
+# Test get_temperature
+@pytest.mark.parametrize(
+    "pipe_fixture, m_flow, T_f_in, segment_ratios, T_b, z, expected", [
+    # Double U-tube
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, 1., 65., np.array([4.33561246, -0.53401739,  3.03985865, 0.28974217])),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, np.array([1., 2., 3., 1.]), 65., np.array([4.40351925, -0.44632268,  3.43990994, 0.77984857])),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([4.45571533, -0.3797212, 3.59641415, 0.96525637])),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, 1., np.array([65., 75.]), np.array([[4.33561246, -0.53401739, 3.03985865, 0.28974217], [4.2430049, -0.47133142, 3.10196228, 0.25289182]])),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.40351925, -0.44632268, 3.43990994, 0.77984857], [4.33450267, -0.35334337, 3.48624558, 0.72533714]])),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.45571533, -0.3797212, 3.59641415, 0.96525637], [4.40441296, -0.26381686, 3.62735271, 0.89129191]])),
+    # Triple U-tube
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, 1., 65., np.array([4.42233734, -0.45836746, 5.35867251, 3.21820626, 0.54732768, 2.80410243])),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, np.array([1., 2., 3., 1.]), 65., np.array([4.47430559, -0.3726503, 5.47798812, 3.56864252, 0.98391184, 3.36084568])),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), 65., np.array([4.51735491, -0.31260448, 5.56070713, 3.7101146, 1.15077001, 3.55786056])),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, 1., np.array([65., 75.]), np.array([[4.42233734, -0.45836746, 5.35867251, 3.21820626, 0.54732768, 2.80410243], [4.33960106, -0.38478969, 5.148416, 3.27807117, 0.49886893, 2.90891002]])),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.47430559, -0.3726503, 5.47798812, 3.56864252, 0.98391184, 3.36084568], [4.41115886, -0.27214867, 5.30491836, 3.61292754, 0.92062819, 3.45008461]])),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([65., 75.]), np.array([[4.51735491, -0.31260448, 5.56070713, 3.7101146, 1.15077001, 3.55786056], [4.4690976, -0.1924684, 5.41547435, 3.74095158, 1.07113754, 3.62431085]])),
+    ])
+def test_temperature_independent(
+        pipe_fixture, m_flow, T_f_in, segment_ratios, T_b, z, expected,
+        request):
+    # Extract pipe from fixture
+    pipe = request.getfixturevalue(pipe_fixture)
+    # Fluid is propylene-glycol 20%
+    fluid = gt.media.Fluid('MPG', 20.)
+    # Fluid temperatures [degC]
+    T_f = pipe.get_temperature(
+        z, T_f_in, T_b, m_flow, fluid.cp, segment_ratios=segment_ratios)
+    assert np.allclose(T_f, expected)
+
+
+# Test get_outlet_temperature
+@pytest.mark.parametrize(
+    "pipe_fixture, m_flow, T_f_in, segment_ratios, T_b, expected", [
+    # Double U-tube
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, 1., np.array([2.66975268, 0.50433911])),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, np.array([1., 2., 3., 1.]), np.array([3.07152772, 0.9760115])),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([3.25838257, 1.19473512])),
+    # Triple U-tube
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, 1., np.array([2.85759687, 0.84655363, 2.20811155])),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, np.array([1., 2., 3., 1.]), np.array([3.22566425, 1.26691386, 2.70350023])),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([3.39731878, 1.46195771, 2.93172173])),
+    ])
+def test_outlet_temperature_independent(
+        pipe_fixture, m_flow, T_f_in, segment_ratios, T_b, expected, request):
+    # Extract pipe from fixture
+    pipe = request.getfixturevalue(pipe_fixture)
+    # Fluid is propylene-glycol 20%
+    fluid = gt.media.Fluid('MPG', 20.)
+    # Outlet fluid temperatures [degC]
+    T_f_out = pipe.get_outlet_temperature(
+        T_f_in, T_b, m_flow, fluid.cp, segment_ratios=segment_ratios)
+    assert np.allclose(T_f_out, expected)
+
+
+# Test get_inlet_temperature
+@pytest.mark.parametrize(
+    "pipe_fixture, m_flow, Q_f, segment_ratios, T_b, expected", [
+    # Double U-tube
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([-3000., 2000.]), None, 1., np.array([7.40595748, -3.59946781])),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([-3000., 2000.]), None, np.array([1., 2., 3., 1.]), np.array([8.15037424, -2.85931237])),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([-3000., 2000.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([8.49653478, -2.51603365])),
+    # Triple U-tube
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([-3000., 2000., -2500.]), None, 1., np.array([7.87321014, -2.88443189, 8.87646527])),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([-3000., 2000., -2500.]), None, np.array([1., 2., 3., 1.]), np.array([8.62102769, -2.14384786, 9.60745667])),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([-3000., 2000., -2500.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([8.96946474, -1.80024805, 9.94472697])),
+    ])
+def test_inlet_temperature_independent(
+        pipe_fixture, m_flow, Q_f, segment_ratios, T_b, expected, request):
+    # Extract pipe from fixture
+    pipe = request.getfixturevalue(pipe_fixture)
+    # Fluid is propylene-glycol 20%
+    fluid = gt.media.Fluid('MPG', 20.)
+    # Inlet fluid temperatures [degC]
+    T_f_in = pipe.get_inlet_temperature(
+        Q_f, T_b, m_flow, fluid.cp, segment_ratios=segment_ratios)
+    assert np.allclose(T_f_in, expected)
+
+
+# Test get_borehole_heat_extraction_rate
+@pytest.mark.parametrize(
+    "pipe_fixture, m_flow, T_f_in, segment_ratios, T_b, expected", [
+    # Double U-tube
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, 1., -956.00963184),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, np.array([1., 2., 3., 1.]), np.array([-311.65264343, -12.22888682, 289.43868666, -320.65369703])),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-136.63705522, -61.48218337, 405.20011405, -214.40320201])),
+    # Triple U-tube
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, 1., -2508.09408199),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, np.array([1., 2., 3., 1.]), np.array([-754.8446691, -351.75861521, 43.23971775, -704.23081162])),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-321.07744291, -563.32294806, -1.60256316, -449.20675261])),
+    ])
+def test_borehole_heat_extraction_rate_independent(
+        pipe_fixture, m_flow, T_f_in, segment_ratios, T_b, expected, request):
+    # Extract pipe from fixture
+    pipe = request.getfixturevalue(pipe_fixture)
+    # Fluid is propylene-glycol 20%
+    fluid = gt.media.Fluid('MPG', 20.)
+    # Borehole heat extraction rates [W]
+    Q_b = pipe.get_borehole_heat_extraction_rate(
+        T_f_in, T_b, m_flow, fluid.cp, segment_ratios=segment_ratios)
+    assert np.allclose(Q_b, expected)
+
+
+# Test get_fluid_heat_extraction_rate
+@pytest.mark.parametrize(
+    "pipe_fixture, m_flow, T_f_in, segment_ratios, T_b, expected", [
+    # Double U-tube
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, 1., np.array([-1853.37094997, 897.36131814])),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, np.array([1., 2., 3., 1.]), np.array([-1533.81766537, 1178.72112475])),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-1385.20195806, 1309.19311228])),
+    # Triple U-tube
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, 1., np.array([-1703.96836927, 1101.4975216, -1905.62323433])),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, np.array([1., 2., 3., 1.]), np.array([-1411.22460372, 1352.24883733, -1708.61861179])),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), np.array([-1274.69852553, 1468.59548314, -1617.86019738])),
+    ])
+def test_fluid_heat_extraction_rate_independent(
+        pipe_fixture, m_flow, T_f_in, segment_ratios, T_b, expected, request):
+    # Extract pipe from fixture
+    pipe = request.getfixturevalue(pipe_fixture)
+    # Fluid is propylene-glycol 20%
+    fluid = gt.media.Fluid('MPG', 20.)
+    # Fluid heat extraction rates [W]
+    Q_f = pipe.get_fluid_heat_extraction_rate(
+        T_f_in, T_b, m_flow, fluid.cp, segment_ratios=segment_ratios)
+    assert np.allclose(Q_f, expected)
+
+
+# Test get_total_heat_extraction_rate
+@pytest.mark.parametrize(
+    "pipe_fixture, m_flow, T_f_in, segment_ratios, T_b, expected", [
+    # Double U-tube
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, 1., -956.0096318353349),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), None, np.array([1., 2., 3., 1.]), -355.0965406202829),
+    ('double_Utube_independent', np.array([0.2, 0.15]), np.array([5., -1.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -76.00884577783381),
+    # Triple U-tube
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, 1., -2508.094081991294),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), None, np.array([1., 2., 3., 1.]), -1767.5943781810597),
+    ('triple_Utube_independent', np.array([0.2, 0.15, 0.10]), np.array([5., -1., 7.]), np.array([0.1, 0.35, 0.40, 0.15]), np.array([1., 2., 3., 1.]), -1423.9632397638252),
+    ])
+def test_total_heat_extraction_rate_independent(
+        pipe_fixture, m_flow, T_f_in, segment_ratios, T_b, expected, request):
+    # Extract pipe from fixture
+    pipe = request.getfixturevalue(pipe_fixture)
+    # Fluid is propylene-glycol 20%
+    fluid = gt.media.Fluid('MPG', 20.)
+    # Total heat extraction rate [W]
+    Q_t = pipe.get_total_heat_extraction_rate(
+        T_f_in, T_b, m_flow, fluid.cp, segment_ratios=segment_ratios)
+    assert np.isclose(Q_t, expected)