From 101cda5ac38d822b1bcb31b2e763ef61ae3887f5 Mon Sep 17 00:00:00 2001 From: chavlin Date: Tue, 10 Dec 2024 15:01:04 -0500 Subject: [PATCH] add cartesian bbox calculation from spherical bounding box edges --- yt_idv/coordinate_utilities.pyx | 131 +++++++++++++++++++++--- yt_idv/scene_data/block_collection.py | 36 +++---- yt_idv/tests/test_spherical_bboxes.py | 138 ++++++++++++++++---------- 3 files changed, 217 insertions(+), 88 deletions(-) diff --git a/yt_idv/coordinate_utilities.pyx b/yt_idv/coordinate_utilities.pyx index cac0c91..0f7df45 100644 --- a/yt_idv/coordinate_utilities.pyx +++ b/yt_idv/coordinate_utilities.pyx @@ -377,28 +377,43 @@ def cartesian_bboxes(MixedCoordBBox bbox_handler, np.float64_t[:] dpos0, np.float64_t[:] dpos1, np.float64_t[:] dpos2, - np.float64_t[:] x, - np.float64_t[:] y, - np.float64_t[:] z, - np.float64_t[:] dx, - np.float64_t[:] dy, - np.float64_t[:] dz, ): - # calculates the cartesian bounding boxes around non-cartesian - # volume elements - # - # bbox_handler : a MixedCoordBBox child instance - # pos0, pos1, pos2: native coordinates of element centers - # dpos0, dpos1, dpos2: element widths in native coordinates - # x, y, z: cartesian centers of bounding boxes, modified in place - # dx, dy, dz : full-widths of the cartesian bounding boxes, modified in place + """ + Calculate the cartesian bounding boxes around non-cartesian volume elements + + Parameters + ---------- + bbox_handler + a MixedCoordBBox child instance + pos0, pos1, pos2 + native coordinates of element centers + dpos0, dpos1, dpos2 + element widths in native coordinates + + Returns + ------- + x_y_z + a 3-element tuples containing the cartesian bounding box centers + dx_y_z + a 3-element tuples containing the cartesian bounding box widths + + """ cdef int i, n_pos, i_result cdef np.float64_t[3] xyz_i cdef np.float64_t[3] dxyz_i cdef int failure = 0 + cdef np.float64_t[:] x, y, z # bbox centers + cdef np.float64_t[:] dx, dy, dz # bbox full widths n_pos = pos0.size + x = np.empty_like(pos0) + y = np.empty_like(pos0) + z = np.empty_like(pos0) + dx = np.empty_like(pos0) + dy = np.empty_like(pos0) + dz = np.empty_like(pos0) + with nogil: for i in range(n_pos): @@ -421,3 +436,91 @@ def cartesian_bboxes(MixedCoordBBox bbox_handler, if failure > 0: raise RuntimeError("Unexpected error in get_cartesian_bbox.") + + x_y_z = (np.asarray(x), np.asarray(y), np.asarray(z)) + dx_y_z = (np.asarray(dx), np.asarray(dy), np.asarray(dz)) + return x_y_z, dx_y_z + + +@cython.cdivision(True) +@cython.boundscheck(False) +@cython.wraparound(False) +def cartesian_bboxes_edges(MixedCoordBBox bbox_handler, + np.float64_t[:] le0, + np.float64_t[:] le1, + np.float64_t[:] le2, + np.float64_t[:] re0, + np.float64_t[:] re1, + np.float64_t[:] re2, + ): + """ + Calculate the cartesian bounding boxes around non-cartesian volume elements + + Same as cartesian_bboxes, but supplying and returning element left/right edges. + + Parameters + ---------- + bbox_handler + a MixedCoordBBox child instance + le0, le1, le2 + native coordinates of element left edges + re0, re1, re2 + native coordinates of element right edges + + Returns + ------- + xyz_le + a 3-element tuples containing the cartesian bounding box left edges + xyz_re + a 3-element tuples containing the cartesian bounding box right edges + + """ + cdef int i, n_pos, i_result + cdef np.float64_t[3] xyz_i + cdef np.float64_t[3] dxyz_i + cdef np.float64_t pos0, pos1, pos2, dpos0, dpos1, dpos2 + cdef int failure = 0 + cdef np.float64_t[:] x_le, y_le, z_le # bbox left edges + cdef np.float64_t[:] x_re, y_re, z_re # bbox right edges + + n_pos = le0.size + x_le = np.empty_like(le0) + y_le = np.empty_like(le0) + z_le = np.empty_like(le0) + x_re = np.empty_like(le0) + y_re = np.empty_like(le0) + z_re = np.empty_like(le0) + + with nogil: + for i in range(n_pos): + pos0 = (le0[i] + re0[i]) / 2. + pos1 = (le1[i] + re1[i]) / 2. + pos2 = (le2[i] + re2[i]) / 2. + dpos0 = re0[i] - le0[i] + dpos1 = re1[i] - le1[i] + dpos2 = re2[i] - le2[i] + i_result = bbox_handler.get_cartesian_bbox(pos0, + pos1, + pos2, + dpos0, + dpos1, + dpos2, + xyz_i, + dxyz_i) + failure += i_result + + x_le[i] = xyz_i[0] - dxyz_i[0] / 2. + x_re[i] = xyz_i[0] + dxyz_i[0] / 2. + y_le[i] = xyz_i[1] - dxyz_i[1] / 2. + y_re[i] = xyz_i[1] + dxyz_i[1] / 2. + z_le[i] = xyz_i[2] - dxyz_i[2] / 2. + z_re[i] = xyz_i[2] + dxyz_i[2] / 2. + + + if failure > 0: + raise RuntimeError("Unexpected error in get_cartesian_bbox.") + + xyz_le = (np.asarray(x_le), np.asarray(y_le), np.asarray(z_le)) + xyz_re = (np.asarray(x_re), np.asarray(y_re), np.asarray(z_re)) + + return xyz_le, xyz_re diff --git a/yt_idv/scene_data/block_collection.py b/yt_idv/scene_data/block_collection.py index 1e7873b..280ee94 100644 --- a/yt_idv/scene_data/block_collection.py +++ b/yt_idv/scene_data/block_collection.py @@ -123,32 +123,26 @@ def _set_geometry_attributes(self, le, re): elif self._yt_geom_str == "spherical": from yt_idv.coordinate_utilities import ( SphericalMixedCoordBBox, - cartesian_bboxes, + cartesian_bboxes_edges, ) axis_id = self.data_source.ds.coordinates.axis_id - # cartesian bbox calcualtions - # TODO: clean this up by rewriting the cython a bit... - widths = re - le - centers = (le + re) / 2.0 + # cartesian bbox calculations bbox_handler = SphericalMixedCoordBBox() - r = centers[:, axis_id["r"]].astype("float64") - theta = centers[:, axis_id["theta"]].astype("float64") - phi = centers[:, axis_id["phi"]].astype("float64") - dr = widths[:, axis_id["r"]].astype("float64") - dtheta = widths[:, axis_id["theta"]].astype("float64") - dphi = widths[:, axis_id["phi"]].astype("float64") - x = np.full(r.shape, np.nan, dtype="float64") - y = np.full(r.shape, np.nan, dtype="float64") - z = np.full(r.shape, np.nan, dtype="float64") - dx = np.full(r.shape, np.nan, dtype="float64") - dy = np.full(r.shape, np.nan, dtype="float64") - dz = np.full(r.shape, np.nan, dtype="float64") - cartesian_bboxes( - bbox_handler, r, theta, phi, dr, dtheta, dphi, x, y, z, dx, dy, dz + + r_le = le[:, axis_id["r"]].astype("float64") + r_re = re[:, axis_id["r"]].astype("float64") + theta_le = le[:, axis_id["theta"]].astype("float64") + theta_re = re[:, axis_id["theta"]].astype("float64") + phi_le = le[:, axis_id["phi"]].astype("float64") + phi_re = re[:, axis_id["phi"]].astype("float64") + + xyz_le, xyz_re = cartesian_bboxes_edges( + bbox_handler, r_le, theta_le, phi_le, r_re, theta_re, phi_re ) - le_cart = np.column_stack([x - dx / 2.0, y - dy / 2.0, z - dz / 2.0]) - re_cart = np.column_stack([x + dx / 2.0, y + dy / 2.0, z + dz / 2.0]) + + le_cart = np.column_stack(xyz_le) + re_cart = np.column_stack(xyz_re) # cartesian le, re, width of whole domain domain_le = [] diff --git a/yt_idv/tests/test_spherical_bboxes.py b/yt_idv/tests/test_spherical_bboxes.py index 088c664..81734fe 100644 --- a/yt_idv/tests/test_spherical_bboxes.py +++ b/yt_idv/tests/test_spherical_bboxes.py @@ -4,6 +4,7 @@ from yt_idv.coordinate_utilities import ( SphericalMixedCoordBBox, cartesian_bboxes, + cartesian_bboxes_edges, cartesian_to_spherical, spherical_to_cartesian, ) @@ -23,66 +24,59 @@ def test_cartesian_bboxes_for_spherical(): phi = np.array([0.05]) dphi = np.array([0.05]) - x = np.full(r.shape, np.nan, dtype="float64") - y = np.full(r.shape, np.nan, dtype="float64") - z = np.full(r.shape, np.nan, dtype="float64") - dx = np.full(r.shape, np.nan, dtype="float64") - dy = np.full(r.shape, np.nan, dtype="float64") - dz = np.full(r.shape, np.nan, dtype="float64") - bbox_handler = SphericalMixedCoordBBox() - cartesian_bboxes(bbox_handler, r, theta, phi, dr, dtheta, dphi, x, y, z, dx, dy, dz) - assert z + dz / 2 == 1.0 - assert np.allclose(x - dx / 2, 0.0) - assert np.allclose(y - dy / 2, 0.0) + xyz, dxyz = cartesian_bboxes(bbox_handler, r, theta, phi, dr, dtheta, dphi) + assert xyz[2] + dxyz[2] / 2 == 1.0 + assert np.allclose(xyz[0] - dxyz[0] / 2, 0.0) + assert np.allclose(xyz[1] - dxyz[1] / 2, 0.0) # now theta = np.pi theta = np.array([np.pi - dtheta[0] / 2]) - cartesian_bboxes(bbox_handler, r, theta, phi, dr, dtheta, dphi, x, y, z, dx, dy, dz) - assert z - dz / 2 == -1.0 - assert np.allclose(x - dx / 2, 0.0) - assert np.allclose(y - dy / 2, 0.0) + xyz, dxyz = cartesian_bboxes(bbox_handler, r, theta, phi, dr, dtheta, dphi) + assert xyz[2] - dxyz[2] / 2 == -1.0 + assert np.allclose(xyz[0] - dxyz[0] / 2, 0.0) + assert np.allclose(xyz[1] - dxyz[1] / 2, 0.0) # element at equator, overlapping the +y axis theta = np.array([np.pi / 2]) phi = np.array([np.pi / 2]) - cartesian_bboxes(bbox_handler, r, theta, phi, dr, dtheta, dphi, x, y, z, dx, dy, dz) + xyz, dxyz = cartesian_bboxes(bbox_handler, r, theta, phi, dr, dtheta, dphi) - assert y + dy / 2 == 1.0 - assert np.allclose(x, 0.0) - assert np.allclose(z, 0.0) + assert xyz[1] + dxyz[1] / 2 == 1.0 + assert np.allclose(xyz[0], 0.0) + assert np.allclose(xyz[2], 0.0) # element at equator, overlapping the -x axis phi = np.array([np.pi]) - cartesian_bboxes(bbox_handler, r, theta, phi, dr, dtheta, dphi, x, y, z, dx, dy, dz) + xyz, dxyz = cartesian_bboxes(bbox_handler, r, theta, phi, dr, dtheta, dphi) - assert x - dx / 2 == -1.0 - assert np.allclose(y, 0.0) - assert np.allclose(z, 0.0) + assert xyz[0] - dxyz[0] / 2 == -1.0 + assert np.allclose(xyz[1], 0.0) + assert np.allclose(xyz[2], 0.0) # element at equator, overlapping the -y axis phi = np.array([3 * np.pi / 2]) - cartesian_bboxes(bbox_handler, r, theta, phi, dr, dtheta, dphi, x, y, z, dx, dy, dz) + xyz, dxyz = cartesian_bboxes(bbox_handler, r, theta, phi, dr, dtheta, dphi) - assert y - dy / 2 == -1.0 - assert np.allclose(x, 0.0) - assert np.allclose(z, 0.0) + assert xyz[1] - dxyz[1] / 2 == -1.0 + assert np.allclose(xyz[0], 0.0) + assert np.allclose(xyz[2], 0.0) # element at equator, overlapping +x axis phi = dphi / 2.0 - cartesian_bboxes(bbox_handler, r, theta, phi, dr, dtheta, dphi, x, y, z, dx, dy, dz) - assert x + dx / 2 == 1.0 + xyz, dxyz = cartesian_bboxes(bbox_handler, r, theta, phi, dr, dtheta, dphi) + assert xyz[0] + dxyz[0] / 2 == 1.0 # element with edge on +x axis in -theta direction theta = np.pi / 2 - dtheta / 2 - cartesian_bboxes(bbox_handler, r, theta, phi, dr, dtheta, dphi, x, y, z, dx, dy, dz) - assert x + dx / 2 == 1.0 + xyz, dxyz = cartesian_bboxes(bbox_handler, r, theta, phi, dr, dtheta, dphi) + assert xyz[0] + dxyz[0] / 2 == 1.0 # element with edge on +x axis in +theta direction theta = np.pi / 2 + dtheta / 2 - cartesian_bboxes(bbox_handler, r, theta, phi, dr, dtheta, dphi, x, y, z, dx, dy, dz) - assert x + dx / 2 == 1.0 + xyz, dxyz = cartesian_bboxes(bbox_handler, r, theta, phi, dr, dtheta, dphi) + assert xyz[0] + dxyz[0] / 2 == 1.0 # finally, check that things work OK with a wide range of # angles @@ -104,10 +98,7 @@ def test_cartesian_bboxes_for_spherical(): r_th_ph = [r_th_ph[i].ravel() for i in range(3)] d_r_th_ph = [d_r_th_ph[i].ravel() for i in range(3)] - x_y_z = [np.full(r_th_ph[0].shape, np.nan, dtype="float64") for _ in range(3)] - d_x_y_z = [np.full(r_th_ph[0].shape, np.nan, dtype="float64") for _ in range(3)] - - cartesian_bboxes( + x_y_z, d_x_y_z = cartesian_bboxes( bbox_handler, r_th_ph[0], r_th_ph[1], @@ -115,12 +106,6 @@ def test_cartesian_bboxes_for_spherical(): d_r_th_ph[0], d_r_th_ph[1], d_r_th_ph[2], - x_y_z[0], - x_y_z[1], - x_y_z[2], - d_x_y_z[0], - d_x_y_z[1], - d_x_y_z[2], ) assert np.all(np.isfinite(x_y_z)) @@ -171,10 +156,7 @@ def test_large_elements(n_angles): r_th_ph = [r_th_ph[i].ravel() for i in range(3)] d_r_th_ph = [d_r_th_ph[i].ravel() for i in range(3)] - x_y_z = [np.full(r_th_ph[0].shape, np.nan, dtype="float64") for _ in range(3)] - d_x_y_z = [np.full(r_th_ph[0].shape, np.nan, dtype="float64") for _ in range(3)] - - cartesian_bboxes( + x_y_z, d_x_y_z = cartesian_bboxes( bbox_handler, r_th_ph[0], r_th_ph[1], @@ -182,12 +164,6 @@ def test_large_elements(n_angles): d_r_th_ph[0], d_r_th_ph[1], d_r_th_ph[2], - x_y_z[0], - x_y_z[1], - x_y_z[2], - d_x_y_z[0], - d_x_y_z[1], - d_x_y_z[2], ) x_y_z = np.column_stack(x_y_z) @@ -201,3 +177,59 @@ def test_large_elements(n_angles): assert np.all(le == -1.0) assert np.all(re == 1.0) + + +def test_spherical_boxes_edges(): + # check that you get the same result from supplying centers + # vs edges. + + bbox_handler = SphericalMixedCoordBBox() + + r_edges = np.linspace(0.0, 10.0, 20, dtype="float64") + theta_edges = np.linspace(0, np.pi, 20, dtype="float64") + phi_edges = np.linspace(0.0, 2 * np.pi, 20, dtype="float64") + + r = (r_edges[0:-1] + r_edges[1:]) / 2.0 + theta = (theta_edges[0:-1] + theta_edges[1:]) / 2.0 + phi = (phi_edges[0:-1] + phi_edges[1:]) / 2.0 + + dr = r_edges[1:] - r_edges[:-1] + dtheta = theta_edges[1:] - theta_edges[:-1] + dphi = phi_edges[1:] - phi_edges[:-1] + + r_th_ph = np.meshgrid(r, theta, phi) + d_r_th_ph = np.meshgrid(dr, dtheta, dphi) + r_th_ph = [r_th_ph[i].ravel() for i in range(3)] + d_r_th_ph = [d_r_th_ph[i].ravel() for i in range(3)] + + x_y_z, d_x_y_z = cartesian_bboxes( + bbox_handler, + r_th_ph[0], + r_th_ph[1], + r_th_ph[2], + d_r_th_ph[0], + d_r_th_ph[1], + d_r_th_ph[2], + ) + x_y_z = np.column_stack(x_y_z) + d_x_y_z = np.column_stack(d_x_y_z) + + le = [r_th_ph[i] - d_r_th_ph[i] / 2.0 for i in range(3)] + re = [r_th_ph[i] + d_r_th_ph[i] / 2.0 for i in range(3)] + xyz_le, xyz_re = cartesian_bboxes_edges( + bbox_handler, + le[0], + le[1], + le[2], + re[0], + re[1], + re[2], + ) + xyz_le = np.column_stack(xyz_le) + xyz_re = np.column_stack(xyz_re) + + centers = (xyz_le + xyz_re) / 2.0 + widths = xyz_re - xyz_le + + assert np.allclose(centers, x_y_z) + assert np.allclose(widths, d_x_y_z)