Add test for new circumcenter()

Author: xylar

conda_package/tests/test_mesh_creation_util.py

@@ -0,0 +1,113 @@
+import collections
+
+import numpy as np
+
+from mpas_tools.mesh.creation.util import circumcenter
+
+point = collections.namedtuple('Point', ['x', 'y', 'z'])
+
+
+def test_circumcenter_planar():
+    # Use provided cell coordinates and vertex as expected output
+    x1, y1, z1 = 0.0, 17320.50807569, 0.0
+    x2, y2, z2 = 10000.0, 17320.50807569, 0.0
+    x3, y3, z3 = 5000.0, 25980.76211353, 0.0
+    expected = (5000.0, 20207.25942164, 0.0)
+
+    xv, yv, zv = circumcenter(False, x1, y1, z1, x2, y2, z2, x3, y3, z3)
+    assert np.allclose([xv, yv, zv], expected, atol=1e-8), (
+        f'Planar circumcenter failed: got {[xv, yv, zv]}, expected {expected}'
+    )
+
+    xv_old, yv_old, zv_old = _old_circumcenter(
+        False, x1, y1, z1, x2, y2, z2, x3, y3, z3
+    )
+    # Print both sets and their difference
+    print('Planar circumcenter:')
+    print('  circumcenter:      ', (xv, yv, zv))
+    print('  _old_circumcenter: ', (xv_old, yv_old, zv_old))
+    print('  difference:        ', (xv - xv_old, yv - yv_old, zv - zv_old))
+    # Check for exact equality
+    assert np.array_equal([xv, yv, zv], [xv_old, yv_old, zv_old]), (
+        f'circumcenter and _old_circumcenter results differ: '
+        f'({xv}, {yv}, {zv}) != ({xv_old}, {yv_old}, {zv_old})'
+    )
+
+
+def test_circumcenter_sphere():
+    # Use provided cell coordinates and vertex as expected output
+    x1, y1, z1 = -3590.95704026, -684.52527287, -6371218.95125671
+    x2, y2, z2 = -9926.94769501, 18734.26792821, -6371184.72274307
+    x3, y3, z3 = 12921.92507847, 11340.69768405, -6371196.8028643
+    expected = (296.83807146, 11327.05862805, -6371209.92418473)
+
+    xv, yv, zv = circumcenter(True, x1, y1, z1, x2, y2, z2, x3, y3, z3)
+    assert np.allclose([xv, yv, zv], expected, atol=1e-8), (
+        f'Spherical circumcenter failed: got {[xv, yv, zv]}, expected '
+        f'{expected}'
+    )
+
+    xv_old, yv_old, zv_old = _old_circumcenter(
+        True, x1, y1, z1, x2, y2, z2, x3, y3, z3
+    )
+    # Print both sets and their difference
+    print('Spherical circumcenter:')
+    print('  circumcenter:      ', (xv, yv, zv))
+    print('  _old_circumcenter: ', (xv_old, yv_old, zv_old))
+    print('  difference:        ', (xv - xv_old, yv - yv_old, zv - zv_old))
+    # Check for exact equality
+    assert np.array_equal([xv, yv, zv], [xv_old, yv_old, zv_old]), (
+        f'circumcenter and _old_circumcenter results differ: '
+        f'({xv}, {yv}, {zv}) != ({xv_old}, {yv_old}, {zv_old})'
+    )
+
+
+def _old_circumcenter(on_sphere, x1, y1, z1, x2, y2, z2, x3, y3, z3):
+    """
+    Compute the circumcenter of the triangle (possibly on a sphere)
+    with the three given vertices in Cartesian coordinates.
+
+    Returns
+    -------
+    center : point
+        The circumcenter of the triangle with x, y and z attributes
+
+    """
+    p1 = point(x1, y1, z1)
+    p2 = point(x2, y2, z2)
+    p3 = point(x3, y3, z3)
+    if on_sphere:
+        a = (p2.x - p3.x) ** 2 + (p2.y - p3.y) ** 2 + (p2.z - p3.z) ** 2
+        b = (p3.x - p1.x) ** 2 + (p3.y - p1.y) ** 2 + (p3.z - p1.z) ** 2
+        c = (p1.x - p2.x) ** 2 + (p1.y - p2.y) ** 2 + (p1.z - p2.z) ** 2
+
+        pbc = a * (-a + b + c)
+        apc = b * (a - b + c)
+        abp = c * (a + b - c)
+
+        xv = (pbc * p1.x + apc * p2.x + abp * p3.x) / (pbc + apc + abp)
+        yv = (pbc * p1.y + apc * p2.y + abp * p3.y) / (pbc + apc + abp)
+        zv = (pbc * p1.z + apc * p2.z + abp * p3.z) / (pbc + apc + abp)
+    else:
+        d = 2 * (
+            p1.x * (p2.y - p3.y) + p2.x * (p3.y - p1.y) + p3.x * (p1.y - p2.y)
+        )
+
+        xv = (
+            (p1.x**2 + p1.y**2) * (p2.y - p3.y)
+            + (p2.x**2 + p2.y**2) * (p3.y - p1.y)
+            + (p3.x**2 + p3.y**2) * (p1.y - p2.y)
+        ) / d
+        yv = (
+            (p1.x**2 + p1.y**2) * (p3.x - p2.x)
+            + (p2.x**2 + p2.y**2) * (p1.x - p3.x)
+            + (p3.x**2 + p3.y**2) * (p2.x - p1.x)
+        ) / d
+        zv = 0.0
+
+        # Optional method to use barycenter instead.
+        # xv = p1.x + p2.x + p3.x
+        # xv = xv / 3.0
+        # yv = p1.y + p2.y + p3.y
+        # yv = yv / 3.0
+    return point(xv, yv, zv)