Centre for cylinder & spherical meshes

openmc/mesh.py

@@ -165,7 +165,7 @@ def _grids(self):
         pass
 
     @property
-    def vertices(self):
+    def vertices(self):  # TODO needs changing
         """Return coordinates of mesh vertices.
 
         Returns
@@ -178,7 +178,7 @@ def vertices(self):
         return np.stack(np.meshgrid(*self._grids, indexing='ij'), axis=-1)
 
     @property
-    def centroids(self):
+    def centroids(self):  # TODO needs changing
         """Return coordinates of mesh element centroids.
 
         Returns
@@ -1052,6 +1052,7 @@ def __init__(self, mesh_id=None, name=''):
         self._r_grid = None
         self._phi_grid = [0.0, 2*pi]
         self._z_grid = None
+        self._centre = [0., 0., 0.]
 
     @property
     def dimension(self):
@@ -1063,6 +1064,10 @@ def dimension(self):
     def n_dimension(self):
         return 3
 
+    @property
+    def centre(self):
+        return self._centre
+
     @property
     def r_grid(self):
         return self._r_grid
@@ -1089,6 +1094,11 @@ def indices(self):
                 for p in range(1, np + 1)
                 for r in range(1, nr + 1))
 
+    @centre.setter
+    def centre(self, coords):
+        cv.check_type('mesh r_grid', coords, Iterable, Real)
+        self._centre = np.asarray(coords)
+
     @r_grid.setter
     def r_grid(self, grid):
         cv.check_type('mesh r_grid', grid, Iterable, Real)
@@ -1160,6 +1170,9 @@ def to_xml_element(self):
         subelement = ET.SubElement(element, "z_grid")
         subelement.text = ' '.join(map(str, self.z_grid))
 
+        subelement = ET.SubElement(element, "centre")
+        subelement.text = ' '.join(map(str, self.centre))
+
         return element
 
     @classmethod
@@ -1183,6 +1196,8 @@ def from_xml_element(cls, elem):
         mesh.r_grid = [float(x) for x in get_text(elem, "r_grid").split()]
         mesh.phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()]
         mesh.z_grid = [float(x) for x in get_text(elem, "z_grid").split()]
+        mesh.centre = [float(x) for x in get_text(elem, "centre").split()]
+
         return mesh
 
     @property
@@ -1236,6 +1251,11 @@ def write_data_to_vtk(self, filename, datasets, volume_normalization=True):
         pts_cartesian[:, 0] = r * np.cos(phi)
         pts_cartesian[:, 1] = r * np.sin(phi)
 
+        # offset with centre
+        pts_cartesian[:, 0] = pts_cartesian[:, 0] + self.centre[0]
+        pts_cartesian[:, 1] = pts_cartesian[:, 1] + self.centre[1]
+        pts_cartesian[:, 2] = pts_cartesian[:, 2] + self.centre[2]
+
         return super().write_data_to_vtk(
             points=pts_cartesian,
             filename=filename,
@@ -1284,6 +1304,7 @@ def __init__(self, mesh_id=None, name=''):
         self._r_grid = None
         self._theta_grid = [0, pi]
         self._phi_grid = [0, 2*pi]
+        self._centre = [0., 0., 0.]
 
     @property
     def dimension(self):
@@ -1295,6 +1316,10 @@ def dimension(self):
     def n_dimension(self):
         return 3
 
+    @property
+    def centre(self):
+        return self._centre
+
     @property
     def r_grid(self):
         return self._r_grid
@@ -1321,6 +1346,11 @@ def indices(self):
                 for t in range(1, nt + 1)
                 for r in range(1, nr + 1))
 
+    @centre.setter
+    def centre(self, coords):
+        cv.check_type('mesh r_grid', coords, Iterable, Real)
+        self._centre = np.asarray(coords)
+
     @r_grid.setter
     def r_grid(self, grid):
         cv.check_type('mesh r_grid', grid, Iterable, Real)
@@ -1392,6 +1422,9 @@ def to_xml_element(self):
         subelement = ET.SubElement(element, "phi_grid")
         subelement.text = ' '.join(map(str, self.phi_grid))
 
+        subelement = ET.SubElement(element, "centre")
+        subelement.text = ' '.join(map(str, self.centre))
+
         return element
 
     @classmethod
@@ -1415,6 +1448,8 @@ def from_xml_element(cls, elem):
         mesh.r_grid = [float(x) for x in get_text(elem, "r_grid").split()]
         mesh.theta_grid = [float(x) for x in get_text(elem, "theta_grid").split()]
         mesh.phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()]
+        mesh.centre = [float(x) for x in get_text(elem, "centre").split()]
+
         return mesh
 
     @property
@@ -1470,6 +1505,11 @@ def write_data_to_vtk(self, filename, datasets, volume_normalization=True):
         pts_cartesian[:, 1] = r * np.sin(phi) * np.sin(theta)
         pts_cartesian[:, 2] = r * np.cos(phi)
 
+        # offset with centre
+        pts_cartesian[:, 0] = pts_cartesian[:, 0] + self.centre[0]
+        pts_cartesian[:, 1] = pts_cartesian[:, 1] + self.centre[1]
+        pts_cartesian[:, 2] = pts_cartesian[:, 2] + self.centre[2]
+
         return super().write_data_to_vtk(
             points=pts_cartesian,
             filename=filename,

src/mesh.cpp

@@ -955,6 +955,7 @@ CylindricalMesh::CylindricalMesh(pugi::xml_node node) : StructuredMesh {node}
   grid_[0] = get_node_array<double>(node, "r_grid");
   grid_[1] = get_node_array<double>(node, "phi_grid");
   grid_[2] = get_node_array<double>(node, "z_grid");
+  centre_ = get_node_array<double>(node, "centre");
 
   if (int err = set_grid()) {
     fatal_error(openmc_err_msg);
@@ -983,6 +984,11 @@ StructuredMesh::MeshIndex CylindricalMesh::get_indices(
       mapped_r[1] += 2 * M_PI;
   }
 
+  // TODO: pass centre as argument 
+  mapped_r[0] += centre_[0]
+  mapped_r[1] += centre_[1]
+  mapped_r[2] += centre_[2]
+
   MeshIndex idx = StructuredMesh::get_indices(mapped_r, in_mesh);
 
   idx[1] = sanitize_phi(idx[1]);
@@ -1187,6 +1193,7 @@ SphericalMesh::SphericalMesh(pugi::xml_node node) : StructuredMesh {node}
   grid_[0] = get_node_array<double>(node, "r_grid");
   grid_[1] = get_node_array<double>(node, "theta_grid");
   grid_[2] = get_node_array<double>(node, "phi_grid");
+  centre_ = get_node_array<double>(node, "centre");
 
   if (int err = set_grid()) {
     fatal_error(openmc_err_msg);
@@ -1216,6 +1223,14 @@ StructuredMesh::MeshIndex SphericalMesh::get_indices(
       mapped_r[2] += 2 * M_PI;
   }
 
+  // test
+  coucou = grid_[0]
+
+  // TODO: pass centre as argument 
+  mapped_r[0] += centre_[0]
+  mapped_r[1] += centre_[1]
+  mapped_r[2] += centre_[2]
+
   MeshIndex idx = StructuredMesh::get_indices(mapped_r, in_mesh);
 
   idx[1] = sanitize_theta(idx[1]);

tests/unit_tests/test_cylindrical_mesh.py

@@ -0,0 +1,30 @@
+import openmc
+import numpy as np
+
+def test_centre_read_write_to_xml():
+    """Tests that the centre attribute can be written and read back to XML
+    """
+    # build
+    mesh = openmc.CylindricalMesh()
+    mesh.phi_grid = [1, 2, 3]
+    mesh.z_grid = [1, 2, 3]
+    mesh.r_grid = [1, 2, 3]
+
+    mesh.centre = [0.1, 0.2, 0.3]
+
+    tally = openmc.Tally()
+
+    mesh_filter = openmc.MeshFilter(mesh)
+    tally.filters.append(mesh_filter)
+
+    tally.scores.append("heating")
+
+    tallies = openmc.Tallies([tally])
+
+    tallies.export_to_xml()
+
+    # read back
+    new_tallies = openmc.Tallies.from_xml()
+    new_tally = new_tallies[0]
+    new_mesh = new_tally.filters[0].mesh
+    assert np.allclose(new_mesh.centre, mesh.centre)

tests/unit_tests/test_spherical_mesh.py

@@ -0,0 +1,30 @@
+import openmc
+import numpy as np
+
+def test_centre_read_write_to_xml():
+    """Tests that the centre attribute can be written and read back to XML
+    """
+    # build
+    mesh = openmc.SphericalMesh()
+    mesh.phi_grid = [1, 2, 3]
+    mesh.theta_grid = [1, 2, 3]
+    mesh.r_grid = [1, 2, 3]
+
+    mesh.centre = [0.1, 0.2, 0.3]
+
+    tally = openmc.Tally()
+
+    mesh_filter = openmc.MeshFilter(mesh)
+    tally.filters.append(mesh_filter)
+
+    tally.scores.append("heating")
+
+    tallies = openmc.Tallies([tally])
+
+    tallies.export_to_xml()
+
+    # read back
+    new_tallies = openmc.Tallies.from_xml()
+    new_tally = new_tallies[0]
+    new_mesh = new_tally.filters[0].mesh
+    assert np.allclose(new_mesh.centre, mesh.centre)