crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge pull request #1859 from nelsonag/cell_api

Browse source 
Addition of C-API methods to set/get a Cell's rotation and translation
This commit is contained in:
Paul Romano 2021-08-02 13:39:37 -05:00 • committed by GitHub
commit 1a4cb28211
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
5 changed files with 232 additions and 38 deletions
Download patch file Download diff file Expand all files Collapse all files

4
include/openmc/capi.h
View file

@ -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);

4
include/openmc/cell.h
View file

@ -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<double>& rot);
//! Get the name of a cell
//! \return Cell name
const std::string& name() const { return name_; };

65
openmc/lib/cell.py
View file

@ -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

148
src/cell.cpp
View file

@ -248,6 +248,48 @@ BoundingBox Universe::bounding_box() const {
// Cell implementation
//==============================================================================
void
Cell::set_rotation(const vector<double>& 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<double>(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<double> 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)
{

49
tests/unit_tests/test_lib.py
View file

@ -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.])