diff --git a/include/openmc/capi.h b/include/openmc/capi.h index f08ec5c02f..012c7c18a0 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -14,12 +14,16 @@ extern "C" { int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n); int openmc_cell_get_id(int32_t index, int32_t* id); int openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T); + int openmc_cell_get_translation(int32_t index, double xyz[]); + int openmc_cell_get_rotation(int32_t index, double rot[], size_t* n); int openmc_cell_get_name(int32_t index, const char** name); int openmc_cell_get_num_instances(int32_t index, int32_t* num_instances); int openmc_cell_set_name(int32_t index, const char* name); int openmc_cell_set_fill(int32_t index, int type, int32_t n, const int32_t* indices); int openmc_cell_set_id(int32_t index, int32_t id); int openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance, bool set_contained = false); + int openmc_cell_set_translation(int32_t index, const double xyz[]); + int openmc_cell_set_rotation(int32_t index, const double rot[], size_t rot_len); int openmc_energy_filter_get_bins(int32_t index, const double** energies, size_t* n); int openmc_energy_filter_set_bins(int32_t index, size_t n, const double* energies); int openmc_energyfunc_filter_get_energy(int32_t index, size_t* n, const double** energy); diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 634b489de2..8d439b1bbc 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -159,6 +159,10 @@ public: //! temperatures. void set_temperature(double T, int32_t instance = -1, bool set_contained = false); + //! Set the rotation matrix of a cell instance + //! \param[in] rot The rotation matrix of length 3 or 9 + void set_rotation(const vector& rot); + //! Get the name of a cell //! \return Cell name const std::string& name() const { return name_; }; diff --git a/openmc/lib/cell.py b/openmc/lib/cell.py index 7fdc4542db..b53a6be7a7 100644 --- a/openmc/lib/cell.py +++ b/openmc/lib/cell.py @@ -1,7 +1,7 @@ import sys from collections.abc import Mapping, Iterable -from ctypes import c_int, c_int32, c_double, c_char_p, POINTER, c_bool +from ctypes import c_int, c_int32, c_double, c_char_p, POINTER, c_bool, c_size_t from weakref import WeakValueDictionary import numpy as np @@ -35,6 +35,13 @@ _dll.openmc_cell_get_temperature.errcheck = _error_handler _dll.openmc_cell_get_name.argtypes = [c_int32, POINTER(c_char_p)] _dll.openmc_cell_get_name.restype = c_int _dll.openmc_cell_get_name.errcheck = _error_handler +_dll.openmc_cell_get_translation.argtypes = [c_int32, POINTER(c_double)] +_dll.openmc_cell_get_translation.restype = c_int +_dll.openmc_cell_get_translation.errcheck = _error_handler +_dll.openmc_cell_get_rotation.argtypes = [c_int32, POINTER(c_double), + POINTER(c_size_t)] +_dll.openmc_cell_get_rotation.restype = c_int +_dll.openmc_cell_get_rotation.errcheck = _error_handler _dll.openmc_cell_set_name.argtypes = [c_int32, c_char_p] _dll.openmc_cell_set_name.restype = c_int _dll.openmc_cell_set_name.errcheck = _error_handler @@ -49,6 +56,13 @@ _dll.openmc_cell_set_temperature.argtypes = [ c_int32, c_double, POINTER(c_int32), c_bool] _dll.openmc_cell_set_temperature.restype = c_int _dll.openmc_cell_set_temperature.errcheck = _error_handler +_dll.openmc_cell_set_translation.argtypes = [c_int32, POINTER(c_double)] +_dll.openmc_cell_set_translation.restype = c_int +_dll.openmc_cell_set_translation.errcheck = _error_handler +_dll.openmc_cell_set_rotation.argtypes = [ + c_int32, POINTER(c_double), c_size_t] +_dll.openmc_cell_set_rotation.restype = c_int +_dll.openmc_cell_set_rotation.errcheck = _error_handler _dll.openmc_get_cell_index.argtypes = [c_int32, POINTER(c_int32)] _dll.openmc_get_cell_index.restype = c_int _dll.openmc_get_cell_index.errcheck = _error_handler @@ -89,6 +103,13 @@ class Cell(_FortranObjectWithID): Number of unique cell instances bounding_box : 2-tuple of numpy.ndarray Lower-left and upper-right coordinates of bounding box + translation : Iterable of float + 3-D coordinates of the translation vector + rotation : Iterable of float + The rotation matrix or angles of the universe filling the cell. This + can either be a fully specified 3 x 3 rotation matrix or an Iterable + of length 3 with the angles in degrees about the x, y, and z axes, + respectively. """ __instances = WeakValueDictionary() @@ -213,6 +234,48 @@ class Cell(_FortranObjectWithID): _dll.openmc_cell_set_temperature(self._index, T, instance, set_contained) + @property + def translation(self): + translation = np.zeros(3) + _dll.openmc_cell_get_translation( + self._index, translation.ctypes.data_as(POINTER(c_double))) + return translation + + @translation.setter + def translation(self, translation_vec): + vector = np.asarray(translation_vec, dtype=float) + _dll.openmc_cell_set_translation( + self._index, vector.ctypes.data_as(POINTER(c_double))) + + @property + def rotation(self): + rotation_data = np.zeros(12) + rot_size = c_size_t() + + _dll.openmc_cell_get_rotation( + self._index, rotation_data.ctypes.data_as(POINTER(c_double)), + rot_size) + rot_size = rot_size.value + + if rot_size == 9: + return rotation_data[:rot_size].shape(3, 3) + elif rot_size in (0, 12): + # If size is 0, rotation_data[9:] will be zeros. This indicates no + # rotation and is the most straightforward way to always return + # an iterable of floats + return rotation_data[9:] + else: + raise ValueError( + 'Invalid size of rotation matrix: {}'.format(rot_size)) + + @rotation.setter + def rotation(self, rotation_data): + flat_rotation = np.asarray(rotation_data, dtype=float).flatten() + + _dll.openmc_cell_set_rotation( + self._index, flat_rotation.ctypes.data_as(POINTER(c_double)), + c_size_t(len(flat_rotation))) + @property def bounding_box(self): inf = sys.float_info.max diff --git a/src/cell.cpp b/src/cell.cpp index a9e9cc3a5e..97f64e3a1d 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -248,6 +248,48 @@ BoundingBox Universe::bounding_box() const { // Cell implementation //============================================================================== +void +Cell::set_rotation(const vector& rot) { + if (fill_ == C_NONE) { + fatal_error(fmt::format("Cannot apply a rotation to cell {}" + " because it is not filled with another universe", + id_)); + } + + if (rot.size() != 3 && rot.size() != 9) { + fatal_error(fmt::format("Non-3D rotation vector applied to cell {}", id_)); + } + + // Compute and store the rotation matrix. + rotation_.clear(); + rotation_.reserve(rot.size() == 9 ? 9 : 12); + if (rot.size() == 3) { + double phi = -rot[0] * PI / 180.0; + double theta = -rot[1] * PI / 180.0; + double psi = -rot[2] * PI / 180.0; + rotation_.push_back(std::cos(theta) * std::cos(psi)); + rotation_.push_back(-std::cos(phi) * std::sin(psi) + + std::sin(phi) * std::sin(theta) * std::cos(psi)); + rotation_.push_back(std::sin(phi) * std::sin(psi) + + std::cos(phi) * std::sin(theta) * std::cos(psi)); + rotation_.push_back(std::cos(theta) * std::sin(psi)); + rotation_.push_back(std::cos(phi) * std::cos(psi) + + std::sin(phi) * std::sin(theta) * std::sin(psi)); + rotation_.push_back(-std::sin(phi) * std::cos(psi) + + std::cos(phi) * std::sin(theta) * std::sin(psi)); + rotation_.push_back(-std::sin(theta)); + rotation_.push_back(std::sin(phi) * std::cos(theta)); + rotation_.push_back(std::cos(phi) * std::cos(theta)); + + // When user specifies angles, write them at end of vector + rotation_.push_back(rot[0]); + rotation_.push_back(rot[1]); + rotation_.push_back(rot[2]); + } else { + std::copy(rot.begin(), rot.end(), std::back_inserter(rotation_)); + } +} + double Cell::temperature(int32_t instance) const { @@ -572,44 +614,8 @@ CSGCell::CSGCell(pugi::xml_node cell_node) // Read the rotation transform. if (check_for_node(cell_node, "rotation")) { - if (fill_ == C_NONE) { - fatal_error(fmt::format("Cannot apply a rotation to cell {}" - " because it is not filled with another universe", id_)); - } - auto rot {get_node_array(cell_node, "rotation")}; - if (rot.size() != 3 && rot.size() != 9) { - fatal_error(fmt::format( - "Non-3D rotation vector applied to cell {}", id_)); - } - - // Compute and store the rotation matrix. - rotation_.reserve(rot.size() == 9 ? 9 : 12); - if (rot.size() == 3) { - double phi = -rot[0] * PI / 180.0; - double theta = -rot[1] * PI / 180.0; - double psi = -rot[2] * PI / 180.0; - rotation_.push_back(std::cos(theta) * std::cos(psi)); - rotation_.push_back(-std::cos(phi) * std::sin(psi) - + std::sin(phi) * std::sin(theta) * std::cos(psi)); - rotation_.push_back(std::sin(phi) * std::sin(psi) - + std::cos(phi) * std::sin(theta) * std::cos(psi)); - rotation_.push_back(std::cos(theta) * std::sin(psi)); - rotation_.push_back(std::cos(phi) * std::cos(psi) - + std::sin(phi) * std::sin(theta) * std::sin(psi)); - rotation_.push_back(-std::sin(phi) * std::cos(psi) - + std::cos(phi) * std::sin(theta) * std::sin(psi)); - rotation_.push_back(-std::sin(theta)); - rotation_.push_back(std::sin(phi) * std::cos(theta)); - rotation_.push_back(std::cos(phi) * std::cos(theta)); - - // When user specifies angles, write them at end of vector - rotation_.push_back(rot[0]); - rotation_.push_back(rot[1]); - rotation_.push_back(rot[2]); - } else { - std::copy(rot.begin(), rot.end(), std::back_inserter(rotation_)); - } + set_rotation(rot); } } @@ -1299,6 +1305,74 @@ openmc_cell_set_id(int32_t index, int32_t id) } } +//! Return the translation vector of a cell +extern "C" int openmc_cell_get_translation(int32_t index, double xyz[]) +{ + if (index >= 0 && index < model::cells.size()) { + auto& cell = model::cells[index]; + xyz[0] = cell->translation_.x; + xyz[1] = cell->translation_.y; + xyz[2] = cell->translation_.z; + return 0; + } else { + set_errmsg("Index in cells array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } +} + +//! Set the translation vector of a cell +extern "C" int openmc_cell_set_translation(int32_t index, const double xyz[]) +{ + if (index >= 0 && index < model::cells.size()) { + if (model::cells[index]->fill_ == C_NONE) { + set_errmsg(fmt::format("Cannot apply a translation to cell {}" + " because it is not filled with another universe", + index)); + return OPENMC_E_GEOMETRY; + } + model::cells[index]->translation_ = Position(xyz); + return 0; + } else { + set_errmsg("Index in cells array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } +} + +//! Return the rotation matrix of a cell +extern "C" int openmc_cell_get_rotation(int32_t index, double rot[], size_t* n) +{ + if (index >= 0 && index < model::cells.size()) { + auto& cell = model::cells[index]; + *n = cell->rotation_.size(); + std::memcpy(rot, cell->rotation_.data(), *n * sizeof(cell->rotation_[0])); + return 0; + } else { + set_errmsg("Index in cells array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } +} + +//! Set the flattened rotation matrix of a cell +extern "C" int openmc_cell_set_rotation(int32_t index, const double rot[], + size_t rot_len) +{ + if (index >= 0 && index < model::cells.size()) { + if (model::cells[index]->fill_ == C_NONE) { + set_errmsg(fmt::format("Cannot apply a rotation to cell {}" + " because it is not filled with another universe", + index)); + return OPENMC_E_GEOMETRY; + } + std::vector vec_rot(rot, rot + rot_len); + model::cells[index]->set_rotation(vec_rot); + return 0; + } else { + set_errmsg("Index in cells array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } +} + +//! Get the number of instances of the requested cell extern "C" int openmc_cell_get_num_instances(int32_t index, int32_t* num_instances) { diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py index 9a4fc8b9fd..41613832f3 100644 --- a/tests/unit_tests/test_lib.py +++ b/tests/unit_tests/test_lib.py @@ -99,6 +99,21 @@ def lib_run(lib_simulation_init): openmc.lib.run() +@pytest.fixture(scope='module') +def pincell_model_w_univ(): + """Set up a model to test with and delete files when done""" + openmc.reset_auto_ids() + pincell = openmc.examples.pwr_pin_cell() + clad_univ = openmc.Universe(cells=[openmc.Cell(fill=pincell.materials[1])]) + pincell.geometry.root_universe.cells[2].fill = clad_univ + pincell.settings.verbosity = 1 + + # Write XML files in tmpdir + with cdtemp(): + pincell.export_to_xml() + yield + + def test_cell_mapping(lib_init): cells = openmc.lib.cells assert isinstance(cells, Mapping) @@ -693,3 +708,37 @@ def test_trigger_set_n_batches(uo2_trigger_model, mpi_intracomm): # Ensure statepoint was created only at batch 20 when calling set_batches assert not os.path.exists('statepoint.12.h5') assert os.path.exists('statepoint.20.h5') + + +def test_cell_translation(pincell_model_w_univ, mpi_intracomm): + openmc.lib.finalize() + openmc.lib.init(intracomm=mpi_intracomm) + openmc.lib.simulation_init() + # Cell 1 is filled with a material so it has a translation, but we can't + # set it. + cell = openmc.lib.cells[1] + assert cell.translation == pytest.approx([0., 0., 0.]) + with pytest.raises(exc.GeometryError, match='not filled with'): + cell.translation = (1., 0., -1.) + + # Cell 2 was given a universe, so we can assign it a translation vector + cell = openmc.lib.cells[2] + assert cell.translation == pytest.approx([0., 0., 0.]) + # This time we *can* set it + cell.translation = (1., 0., -1.) + assert cell.translation == pytest.approx([1., 0., -1.]) + + +def test_cell_rotation(pincell_model_w_univ): + # Cell 1 is filled with a material so we cannot rotate it, but we can get + # its rotation matrix (which will be the identity matrix) + cell = openmc.lib.cells[1] + assert cell.rotation == pytest.approx([0., 0., 0.]) + with pytest.raises(exc.GeometryError, match='not filled with'): + cell.rotation = (180., 0., 0.) + + # Now repeat with Cell 2 and we will be allowed to do it + cell = openmc.lib.cells[2] + assert cell.rotation == pytest.approx([0., 0., 0.]) + cell.rotation = (180., 0., 0.) + assert cell.rotation == pytest.approx([180., 0., 0.])