From 7479a9022a056d697f79ba2caeccb067eb7aee67 Mon Sep 17 00:00:00 2001
From: Patrick Shriwise <pshriwise@gmail.com>
Date: Tue, 6 Jun 2023 11:01:59 -0500
Subject: [PATCH] Fix spherical coords again. (#2538)

---
 openmc/mesh.py                       |  8 +++---
 tests/unit_tests/test_mesh_to_vtk.py | 42 ++++++++++++++++++++++++++++
 2 files changed, 46 insertions(+), 4 deletions(-)

diff --git a/openmc/mesh.py b/openmc/mesh.py
index c59f5aaf6..fde9f302d 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -1338,7 +1338,7 @@ class CylindricalMesh(StructuredMesh):
         pts_cartesian[:, 0] = r * np.cos(phi) + self.origin[0]
         pts_cartesian[:, 1] = r * np.sin(phi) + self.origin[1]
         pts_cartesian[:, 2] += self.origin[2]
-        
+
         return super().write_data_to_vtk(
             points=pts_cartesian,
             filename=filename,
@@ -1585,9 +1585,9 @@ class SphericalMesh(StructuredMesh):
 
         r, theta, phi = pts_spherical[:, 0], pts_spherical[:, 1], pts_spherical[:, 2]
 
-        pts_cartesian[:, 0] = r * np.sin(phi) * np.cos(theta) + self.origin[0]
-        pts_cartesian[:, 1] = r * np.sin(phi) * np.sin(theta) + self.origin[1]
-        pts_cartesian[:, 2] = r * np.cos(phi) + self.origin[2]
+        pts_cartesian[:, 0] = r * np.sin(theta) * np.cos(phi) + self.origin[0]
+        pts_cartesian[:, 1] = r * np.sin(theta) * np.sin(phi) + self.origin[1]
+        pts_cartesian[:, 2] = r * np.cos(theta) + self.origin[2]
 
         return super().write_data_to_vtk(
             points=pts_cartesian,
diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py
index 4c0a43877..2e2790e0a 100644
--- a/tests/unit_tests/test_mesh_to_vtk.py
+++ b/tests/unit_tests/test_mesh_to_vtk.py
@@ -273,3 +273,45 @@ def test_vtk_write_ordering(model, mesh, surface):
             # need to get flat index with axes reversed due to ordering passed into the VTK file
             flat_idx = np.ravel_multi_index(tuple(ijk[::-1]), mesh.dimension[::-1])
             assert vtk_data[flat_idx] == 0.0, err_msg
+
+
+def test_sphere_mesh_coordinates(run_in_tmpdir):
+    mesh = openmc.SphericalMesh()
+    mesh.r_grid = np.linspace(0.1, 10, 30)
+    mesh.phi_grid = np.linspace(0, 1.5*np.pi, 25)
+    mesh.theta_grid = np.linspace(0, np.pi / 2, 15)
+
+    # write the data to a VTK file (no data)
+    vtk_filename = 'test.vtk'
+    mesh.write_data_to_vtk(vtk_filename, {})
+
+    # read file
+    reader = vtk.vtkStructuredGridReader()
+    reader.SetFileName(str(vtk_filename))
+    reader.Update()
+
+    vtk_grid = reader.GetOutput()
+
+    # create a region that matches the spherical mesh description
+    x = openmc.XPlane()
+    z = openmc.ZPlane()
+    y = openmc.YPlane()
+    s = openmc.Sphere(r=10.0)
+
+    region = +z & +y & -s | -x & -y & +z & -s
+
+    # the VTK interface will update this list when GetCentroid is called
+    centroid = np.zeros(3)
+
+    # ensure all centroids of the sphere mesh are inside the cell region
+    for i in range(vtk_grid.GetNumberOfCells()):
+        # get the cell from the stuctured mesh object
+        cell = vtk_grid.GetCell(i)
+        cell.GetCentroid(centroid)
+
+        # if the coordinate conversion is happening correctly,
+        # every one of the cell centroids should be in the CSG region
+        assert centroid in region, \
+            f'Cell centroid {centroid} not in equivalent ' \
+            f'CSG region for spherical mesh {mesh}'
+