diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 148a28c88..ded41f738 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -14,6 +14,8 @@ 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_name(int32_t index, const char** name); + 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); @@ -29,6 +31,8 @@ extern "C" { int openmc_filter_set_id(int32_t index, int32_t id); int openmc_finalize(); int openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance); + int openmc_cell_bounding_box(const int32_t index, double* llc, double* urc); + int openmc_global_bounding_box(double* llc, double* urc); int openmc_fission_bank(void** ptr, int64_t* n); int openmc_get_cell_index(int32_t id, int32_t* index); int openmc_get_filter_index(int32_t id, int32_t* index); @@ -56,6 +60,8 @@ extern "C" { int openmc_material_set_density(int32_t index, double density, const char* units); int openmc_material_set_densities(int32_t index, int n, const char** name, const double* density); int openmc_material_set_id(int32_t index, int32_t id); + int openmc_material_get_name(int32_t index, const char** name); + int openmc_material_set_name(int32_t index, const char* name); int openmc_material_set_volume(int32_t index, double volume); int openmc_material_filter_get_bins(int32_t index, const int32_t** bins, size_t* n); int openmc_material_filter_set_bins(int32_t index, size_t n, const int32_t* bins); diff --git a/include/openmc/cell.h b/include/openmc/cell.h index afe40fe43..db8d60ed4 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -65,6 +65,8 @@ public: //! \param group_id An HDF5 group id. void to_hdf5(hid_t group_id) const; + BoundingBox bounding_box() const; + std::unique_ptr partitioner_; }; @@ -114,6 +116,9 @@ public: //! \param group_id An HDF5 group id. virtual void to_hdf5(hid_t group_id) const = 0; + //! Get the BoundingBox for this cell. + virtual BoundingBox bounding_box() const = 0; + //---------------------------------------------------------------------------- // Accessors @@ -129,6 +134,14 @@ public: //! all instances is set. void set_temperature(double T, int32_t instance = -1); + //! Get the name of a cell + //! \return Cell name + const std::string& name() const { return name_; }; + + //! Set the temperature of a cell instance + //! \param[in] name Cell name + void set_name(const std::string& name) { name_ = name; }; + //---------------------------------------------------------------------------- // Data members @@ -192,9 +205,14 @@ public: void to_hdf5(hid_t group_id) const; + BoundingBox bounding_box() const; + protected: bool contains_simple(Position r, Direction u, int32_t on_surface) const; bool contains_complex(Position r, Direction u, int32_t on_surface) const; + BoundingBox bounding_box_simple() const; + static BoundingBox bounding_box_complex(std::vector rpn); + static void apply_demorgan(std::vector& rpn); }; //============================================================================== @@ -203,16 +221,19 @@ protected: class DAGCell : public Cell { public: - moab::DagMC* dagmc_ptr_; DAGCell(); - int32_t dag_index_; bool contains(Position r, Direction u, int32_t on_surface) const; std::pair distance(Position r, Direction u, int32_t on_surface) const; + BoundingBox bounding_box() const; + void to_hdf5(hid_t group_id) const; + + moab::DagMC* dagmc_ptr_; //!< Pointer to DagMC instance + int32_t dag_index_; //!< DagMC index of cell }; #endif diff --git a/include/openmc/material.h b/include/openmc/material.h index 02c5380f3..652f3db8e 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -92,6 +92,13 @@ public: //! \return Density in [g/cm^3] double density_gpcc() const { return density_gpcc_; } + //! Get name + //! \return Material name + const std::string& name() const { return name_; } + + //! Set name + void set_name(const std::string& name) { name_ = name; } + //! Set total density of the material // //! \param[in] density Density value diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 7b0ccefd1..d61ae07ce 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -43,12 +43,46 @@ namespace model { struct BoundingBox { - double xmin; - double xmax; - double ymin; - double ymax; - double zmin; - double zmax; + double xmin = -INFTY; + double xmax = INFTY; + double ymin = -INFTY; + double ymax = INFTY; + double zmin = -INFTY; + double zmax = INFTY; + + + inline BoundingBox operator &(const BoundingBox& other) { + BoundingBox result = *this; + return result &= other; + } + + inline BoundingBox operator |(const BoundingBox& other) { + BoundingBox result = *this; + return result |= other; + } + + // intersect operator + inline BoundingBox& operator &=(const BoundingBox& other) { + xmin = std::max(xmin, other.xmin); + xmax = std::min(xmax, other.xmax); + ymin = std::max(ymin, other.ymin); + ymax = std::min(ymax, other.ymax); + zmin = std::max(zmin, other.zmin); + zmax = std::min(zmax, other.zmax); + return *this; + } + + // union operator + inline BoundingBox& operator |=(const BoundingBox& other) { + xmin = std::min(xmin, other.xmin); + xmax = std::max(xmax, other.xmax); + ymin = std::min(ymin, other.ymin); + ymax = std::max(ymax, other.ymax); + zmin = std::min(zmin, other.zmin); + zmax = std::max(zmax, other.zmax); + return *this; + } + }; //============================================================================== @@ -105,6 +139,8 @@ public: //TODO: this probably needs to include i_periodic for PeriodicSurface virtual void to_hdf5(hid_t group_id) const = 0; + //! Get the BoundingBox for this surface. + virtual BoundingBox bounding_box(bool pos_side) const { return {}; } }; class CSGSurface : public Surface @@ -126,18 +162,17 @@ protected: class DAGSurface : public Surface { public: - moab::DagMC* dagmc_ptr_; DAGSurface(); - int32_t dag_index_; double evaluate(Position r) const; double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; Direction reflect(Position r, Direction u) const; - //! Get the bounding box of this surface. - BoundingBox bounding_box() const; void to_hdf5(hid_t group_id) const; + + moab::DagMC* dagmc_ptr_; //!< Pointer to DagMC instance + int32_t dag_index_; //!< DagMC index of surface }; #endif //============================================================================== @@ -166,8 +201,6 @@ public: virtual bool periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const = 0; - //! Get the bounding box for this surface. - virtual BoundingBox bounding_box() const = 0; }; //============================================================================== @@ -186,7 +219,7 @@ public: void to_hdf5_inner(hid_t group_id) const; bool periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const; - BoundingBox bounding_box() const; + BoundingBox bounding_box(bool pos_side) const; double x0_; }; @@ -207,7 +240,7 @@ public: void to_hdf5_inner(hid_t group_id) const; bool periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const; - BoundingBox bounding_box() const; + BoundingBox bounding_box(bool pos_side) const; double y0_; }; @@ -228,7 +261,7 @@ public: void to_hdf5_inner(hid_t group_id) const; bool periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const; - BoundingBox bounding_box() const; + BoundingBox bounding_box(bool pos_side) const; double z0_; }; @@ -249,7 +282,6 @@ public: void to_hdf5_inner(hid_t group_id) const; bool periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const; - BoundingBox bounding_box() const; double A_, B_, C_, D_; }; @@ -269,6 +301,7 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; + BoundingBox bounding_box(bool pos_side) const; double y0_, z0_, radius_; }; @@ -288,6 +321,7 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; + BoundingBox bounding_box(bool pos_side) const; double x0_, z0_, radius_; }; @@ -307,6 +341,7 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; + BoundingBox bounding_box(bool pos_side) const; double x0_, y0_, radius_; }; @@ -326,6 +361,7 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; + BoundingBox bounding_box(bool pos_side) const; double x0_, y0_, z0_, radius_; }; diff --git a/openmc/capi/cell.py b/openmc/capi/cell.py index 959ab08fc..784ebd087 100644 --- a/openmc/capi/cell.py +++ b/openmc/capi/cell.py @@ -1,3 +1,5 @@ +import sys + from collections.abc import Mapping, Iterable from ctypes import c_int, c_int32, c_double, c_char_p, POINTER from weakref import WeakValueDictionary @@ -28,6 +30,12 @@ _dll.openmc_cell_get_temperature.argtypes = [ c_int32, POINTER(c_int32), POINTER(c_double)] _dll.openmc_cell_get_temperature.restype = c_int _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_set_name.argtypes = [c_int32, c_char_p] +_dll.openmc_cell_set_name.restype = c_int +_dll.openmc_cell_set_name.errcheck = _error_handler _dll.openmc_cell_set_fill.argtypes = [ c_int32, c_int, c_int32, POINTER(c_int32)] _dll.openmc_cell_set_fill.restype = c_int @@ -43,6 +51,11 @@ _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 _dll.cells_size.restype = c_int +_dll.openmc_cell_bounding_box.argtypes = [c_int, + POINTER(c_double), + POINTER(c_double)] +_dll.openmc_cell_bounding_box.restype = c_int +_dll.openmc_cell_bounding_box.errcheck = _error_handler class Cell(_FortranObjectWithID): @@ -102,6 +115,17 @@ class Cell(_FortranObjectWithID): def id(self, cell_id): _dll.openmc_cell_set_id(self._index, cell_id) + @property + def name(self): + name = c_char_p() + _dll.openmc_cell_get_name(self._index, name) + return name.value.decode() + + @name.setter + def name(self, name): + name_ptr = c_char_p(name.encode()) + _dll.openmc_cell_set_name(self._index, name_ptr) + @property def fill(self): fill_type = c_int() @@ -166,6 +190,20 @@ class Cell(_FortranObjectWithID): _dll.openmc_cell_set_temperature(self._index, T, instance) + @property + def bounding_box(self): + inf = sys.float_info.max + llc = np.zeros(3) + urc = np.zeros(3) + _dll.openmc_cell_bounding_box(self._index, + llc.ctypes.data_as(POINTER(c_double)), + urc.ctypes.data_as(POINTER(c_double))) + llc[llc == inf] = np.inf + urc[urc == inf] = np.inf + llc[llc == -inf] = -np.inf + urc[urc == -inf] = -np.inf + + return llc, urc class _CellMapping(Mapping): def __getitem__(self, key): diff --git a/openmc/capi/core.py b/openmc/capi/core.py index aae17e06a..a470f0665 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -1,3 +1,5 @@ +import sys + from contextlib import contextmanager from ctypes import (CDLL, c_bool, c_int, c_int32, c_int64, c_double, c_char_p, c_char, POINTER, Structure, c_void_p, create_string_buffer) @@ -73,8 +75,26 @@ _dll.openmc_simulation_finalize.errcheck = _error_handler _dll.openmc_statepoint_write.argtypes = [c_char_p, POINTER(c_bool)] _dll.openmc_statepoint_write.restype = c_int _dll.openmc_statepoint_write.errcheck = _error_handler +_dll.openmc_global_bounding_box.argtypes = [POINTER(c_double), + POINTER(c_double)] +_dll.openmc_global_bounding_box.restype = c_int +_dll.openmc_global_bounding_box.errcheck = _error_handler +def global_bounding_box(): + """Calculate a global bounding box for the model""" + inf = sys.float_info.max + llc = np.zeros(3) + urc = np.zeros(3) + _dll.openmc_global_bounding_box(llc.ctypes.data_as(POINTER(c_double)), + urc.ctypes.data_as(POINTER(c_double))) + llc[llc == inf] = np.inf + urc[urc == inf] = np.inf + llc[llc == -inf] = -np.inf + urc[urc == -inf] = -np.inf + + return llc, urc + def calculate_volumes(): """Run stochastic volume calculation""" _dll.openmc_calculate_volumes() diff --git a/openmc/capi/material.py b/openmc/capi/material.py index 43959e244..f0ecac761 100644 --- a/openmc/capi/material.py +++ b/openmc/capi/material.py @@ -48,6 +48,12 @@ _dll.openmc_material_set_densities.errcheck = _error_handler _dll.openmc_material_set_id.argtypes = [c_int32, c_int32] _dll.openmc_material_set_id.restype = c_int _dll.openmc_material_set_id.errcheck = _error_handler +_dll.openmc_material_get_name.argtypes = [c_int32, POINTER(c_char_p)] +_dll.openmc_material_get_name.restype = c_int +_dll.openmc_material_get_name.errcheck = _error_handler +_dll.openmc_material_set_name.argtypes = [c_int32, c_char_p] +_dll.openmc_material_set_name.restype = c_int +_dll.openmc_material_set_name.errcheck = _error_handler _dll.openmc_material_set_volume.argtypes = [c_int32, c_double] _dll.openmc_material_set_volume.restype = c_int _dll.openmc_material_set_volume.errcheck = _error_handler @@ -124,6 +130,17 @@ class Material(_FortranObjectWithID): def id(self, mat_id): _dll.openmc_material_set_id(self._index, mat_id) + @property + def name(self): + name = c_char_p() + _dll.openmc_material_get_name(self._index, name) + return name.value.decode() + + @name.setter + def name(self, name): + name_ptr = c_char_p(name.encode()) + _dll.openmc_material_set_name(self._index, name_ptr) + @property def volume(self): volume = c_double() diff --git a/src/cell.cpp b/src/cell.cpp index 04273eded..3c079f7ae 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -1,11 +1,12 @@ #include "openmc/cell.h" +#include #include #include #include #include -#include +#include #include "openmc/capi.h" #include "openmc/constants.h" @@ -118,7 +119,6 @@ generate_rpn(int32_t cell_id, std::vector infix) if (token < OP_UNION) { // If token is not an operator, add it to output rpn.push_back(token); - } else if (token < OP_RIGHT_PAREN) { // Regular operators union, intersection, complement while (stack.size() > 0) { @@ -158,7 +158,6 @@ generate_rpn(int32_t cell_id, std::vector infix) << cell_id; fatal_error(err_msg); } - rpn.push_back(stack.back()); stack.pop_back(); } @@ -208,6 +207,19 @@ Universe::to_hdf5(hid_t universes_group) const close_group(group); } +BoundingBox Universe::bounding_box() const { + BoundingBox bbox = {INFTY, -INFTY, INFTY, -INFTY, INFTY, -INFTY}; + if (cells_.size() == 0) { + return {}; + } else { + for (const auto& cell : cells_) { + auto& c = model::cells[cell]; + bbox |= c->bounding_box(); + } + } + return bbox; +} + //============================================================================== // Cell implementation //============================================================================== @@ -567,6 +579,95 @@ CSGCell::to_hdf5(hid_t cell_group) const close_group(group); } +BoundingBox CSGCell::bounding_box_simple() const { + BoundingBox bbox; + for (int32_t token : rpn_) { + bbox &= model::surfaces[abs(token)-1]->bounding_box(token > 0); + } + return bbox; +} + + +void CSGCell::apply_demorgan(std::vector& rpn) { + for (auto& token : rpn) { + if (token < OP_UNION) { token *= -1; } + else if (token == OP_UNION) { token = OP_INTERSECTION; } + else if (token == OP_INTERSECTION) { token = OP_UNION; } + } +} + +BoundingBox CSGCell::bounding_box_complex(std::vector rpn) { + + // if the last operator is a complement op, there is no + // sub-region that the complement connects to. This indicates + // that the entire region is a complement and we can apply + // De Morgan's laws immediately + if (rpn.back() == OP_COMPLEMENT) { + rpn.pop_back(); + apply_demorgan(rpn); + } + + // reverse the rpn to make popping easier + std::reverse(rpn.begin(), rpn.end()); + + BoundingBox current = model::surfaces[abs(rpn.back()) - 1]->bounding_box(rpn.back() > 0); + rpn.pop_back(); + + while (rpn.size()) { + // move through the rpn in twos + int32_t one = rpn.back(); rpn.pop_back(); + int32_t two = rpn.back(); rpn.pop_back(); + + // the first token should always be a surface + Expects(one < OP_UNION); + + if (two >= OP_UNION) { + if (two == OP_UNION) { + current |= model::surfaces[abs(one)-1]->bounding_box(one > 0); + } else if (two == OP_INTERSECTION) { + current &= model::surfaces[abs(one)-1]->bounding_box(one > 0); + } + } else { + // two surfaces in a row (left parenthesis), + // create sub-rpn for region in parenthesis + std::vector subrpn; + subrpn.push_back(one); + subrpn.push_back(two); + // add until last two tokens in the sub-rpn are operators + // (indicates a right parenthesis) + while (!((subrpn.back() >= OP_UNION) && (*(subrpn.rbegin() + 1) >= OP_UNION))) { + subrpn.push_back(rpn.back()); + rpn.pop_back(); + } + + // handle complement case using De Morgan's laws + if (subrpn.back() == OP_COMPLEMENT) { + subrpn.pop_back(); + apply_demorgan(subrpn); + subrpn.push_back(rpn.back()); + rpn.pop_back(); + } + // save the last operator, tells us how to combine this region + // with our current bounding box + int32_t op = subrpn.back(); subrpn.pop_back(); + // get bounding box for the subrpn + BoundingBox sub_box = bounding_box_complex(subrpn); + // combine the sub-rpn bounding box with our current cell box + if (op == OP_UNION) { + current |= sub_box; + } else if (op == OP_INTERSECTION) { + current &= sub_box; + } + } + } + + return current; +} + +BoundingBox CSGCell::bounding_box() const { + return simple_ ? bounding_box_simple() : bounding_box_complex(rpn_); +} + //============================================================================== bool @@ -691,6 +792,16 @@ bool DAGCell::contains(Position r, Direction u, int32_t on_surface) const void DAGCell::to_hdf5(hid_t group_id) const { return; } +BoundingBox DAGCell::bounding_box() const +{ + moab::ErrorCode rval; + moab::EntityHandle vol = dagmc_ptr_->entity_by_index(3, dag_index_); + double min[3], max[3]; + rval = dagmc_ptr_->getobb(vol, min, max); + MB_CHK_ERR_CONT(rval); + return {min[0], max[0], min[1], max[1], min[2], max[2]}; +} + #endif //============================================================================== @@ -983,6 +1094,55 @@ openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T) return 0; } +//! Get the bounding box of a cell +extern "C" int +openmc_cell_bounding_box(const int32_t index, double* llc, double* urc) { + + BoundingBox bbox; + + const auto& c = model::cells[index]; + bbox = c->bounding_box(); + + // set lower left corner values + llc[0] = bbox.xmin; + llc[1] = bbox.ymin; + llc[2] = bbox.zmin; + + // set upper right corner values + urc[0] = bbox.xmax; + urc[1] = bbox.ymax; + urc[2] = bbox.zmax; + + return 0; +} + +//! Get the name of a cell +extern "C" int +openmc_cell_get_name(int32_t index, const char** name) { + if (index < 0 || index >= model::cells.size()) { + set_errmsg("Index in cells array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + + *name = model::cells[index]->name().data(); + + return 0; +} + +//! Set the name of a cell +extern "C" int +openmc_cell_set_name(int32_t index, const char* name) { + if (index < 0 || index >= model::cells.size()) { + set_errmsg("Index in cells array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + + model::cells[index]->set_name(name); + + return 0; +} + + //! Return the index in the cells array of a cell with a given ID extern "C" int openmc_get_cell_index(int32_t id, int32_t* index) diff --git a/src/geometry.cpp b/src/geometry.cpp index 0f89a0ecd..556d10650 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -9,6 +9,7 @@ #include "openmc/lattice.h" #include "openmc/settings.h" #include "openmc/simulation.h" +#include "openmc/string_utils.h" #include "openmc/surface.h" @@ -475,4 +476,21 @@ openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance) return 0; } +extern "C" int openmc_global_bounding_box(double* llc, double* urc) { + auto bbox = model::universes.at(model::root_universe)->bounding_box(); + + // set lower left corner values + llc[0] = bbox.xmin; + llc[1] = bbox.ymin; + llc[2] = bbox.zmin; + + // set upper right corner values + urc[0] = bbox.xmax; + urc[1] = bbox.ymax; + urc[2] = bbox.zmax; + + return 0; +} + + } // namespace openmc diff --git a/src/material.cpp b/src/material.cpp index 0a6bebf17..685fd1da2 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -1381,6 +1381,30 @@ openmc_material_set_id(int32_t index, int32_t id) return 0; } +extern "C" int +openmc_material_get_name(int32_t index, const char** name) { + if (index < 0 || index >= model::materials.size()) { + set_errmsg("Index in materials array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + + *name = model::materials[index]->name().data(); + + return 0; +} + +extern "C" int +openmc_material_set_name(int32_t index, const char* name) { + if (index < 0 || index >= model::materials.size()) { + set_errmsg("Index in materials array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + + model::materials[index]->set_name(name); + + return 0; +} + extern "C" int openmc_material_set_volume(int32_t index, double volume) { diff --git a/src/surface.cpp b/src/surface.cpp index f364694a0..14d077cff 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -271,17 +271,8 @@ Direction DAGSurface::reflect(Position r, Direction u) const return simulation::last_dir; } -BoundingBox DAGSurface::bounding_box() const -{ - moab::ErrorCode rval; - moab::EntityHandle surf = dagmc_ptr_->entity_by_index(2, dag_index_); - double min[3], max[3]; - rval = dagmc_ptr_->getobb(surf, min, max); - MB_CHK_ERR_CONT(rval); - return {min[0], max[0], min[1], max[1], min[2], max[2]}; -} - void DAGSurface::to_hdf5(hid_t group_id) const {} + #endif //============================================================================== // PeriodicSurface implementation @@ -366,9 +357,13 @@ bool SurfaceXPlane::periodic_translate(const PeriodicSurface* other, } BoundingBox -SurfaceXPlane::bounding_box() const +SurfaceXPlane::bounding_box(bool pos_side) const { - return {x0_, x0_, -INFTY, INFTY, -INFTY, INFTY}; + if (pos_side) { + return {x0_, INFTY, -INFTY, INFTY, -INFTY, INFTY}; + } else { + return {-INFTY, x0_, -INFTY, INFTY, -INFTY, INFTY}; + } } //============================================================================== @@ -428,9 +423,13 @@ bool SurfaceYPlane::periodic_translate(const PeriodicSurface* other, } BoundingBox -SurfaceYPlane::bounding_box() const +SurfaceYPlane::bounding_box(bool pos_side) const { - return {-INFTY, INFTY, y0_, y0_, -INFTY, INFTY}; + if (pos_side) { + return {-INFTY, INFTY, y0_, INFTY, -INFTY, INFTY}; + } else { + return {-INFTY, INFTY, -INFTY, y0_, -INFTY, INFTY}; + } } //============================================================================== @@ -474,9 +473,13 @@ bool SurfaceZPlane::periodic_translate(const PeriodicSurface* other, } BoundingBox -SurfaceZPlane::bounding_box() const +SurfaceZPlane::bounding_box(bool pos_side) const { - return {-INFTY, INFTY, -INFTY, INFTY, z0_, z0_}; + if (pos_side) { + return {-INFTY, INFTY, -INFTY, INFTY, z0_, INFTY}; + } else { + return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, z0_}; + } } //============================================================================== @@ -538,12 +541,6 @@ bool SurfacePlane::periodic_translate(const PeriodicSurface* other, Position& r, return false; } -BoundingBox -SurfacePlane::bounding_box() const -{ - return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; -} - //============================================================================== // Generic functions for x-, y-, and z-, cylinders //============================================================================== @@ -645,7 +642,6 @@ Direction SurfaceXCylinder::normal(Position r) const return axis_aligned_cylinder_normal<0, 1, 2>(r, y0_, z0_); } - void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "x-cylinder", false); @@ -653,6 +649,13 @@ void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } +BoundingBox SurfaceXCylinder::bounding_box(bool pos_side) const { + if (!pos_side) { + return {-INFTY, INFTY, y0_ - radius_, y0_ + radius_, z0_ - radius_, z0_ + radius_}; + } else { + return {}; + } +} //============================================================================== // SurfaceYCylinder implementation //============================================================================== @@ -686,6 +689,14 @@ void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } +BoundingBox SurfaceYCylinder::bounding_box(bool pos_side) const { + if (!pos_side) { + return {x0_ - radius_, x0_ + radius_, -INFTY, INFTY, z0_ - radius_, z0_ + radius_}; + } else { + return {}; + } +} + //============================================================================== // SurfaceZCylinder implementation //============================================================================== @@ -719,6 +730,15 @@ void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } +BoundingBox SurfaceZCylinder::bounding_box(bool pos_side) const { + if (!pos_side) { + return {x0_ - radius_, x0_ + radius_, y0_ - radius_, y0_ + radius_, -INFTY, INFTY}; + } else { + return {}; + } +} + + //============================================================================== // SurfaceSphere implementation //============================================================================== @@ -787,6 +807,16 @@ void SurfaceSphere::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } +BoundingBox SurfaceSphere::bounding_box(bool pos_side) const { + if (!pos_side) { + return {x0_ - radius_, x0_ + radius_, + y0_ - radius_, y0_ + radius_, + z0_ - radius_, z0_ + radius_}; + } else { + return {}; + } +} + //============================================================================== // Generic functions for x-, y-, and z-, cones //============================================================================== @@ -1161,7 +1191,7 @@ void read_surfaces(pugi::xml_node node) } // See if this surface makes part of the global bounding box. - BoundingBox bb = surf->bounding_box(); + auto bb = surf->bounding_box(true) & surf->bounding_box(false); if (bb.xmin > -INFTY && bb.xmin < xmin) { xmin = bb.xmin; i_xmin = i_surf; diff --git a/tests/regression_tests/complex_cell/test.py b/tests/regression_tests/complex_cell/test.py index 77cbd6cb7..b43ccd72f 100755 --- a/tests/regression_tests/complex_cell/test.py +++ b/tests/regression_tests/complex_cell/test.py @@ -1,5 +1,8 @@ from tests.testing_harness import TestHarness +import sys + +import openmc.capi def test_complex_cell(): harness = TestHarness('statepoint.10.h5') diff --git a/tests/regression_tests/dagmc/legacy/test.py b/tests/regression_tests/dagmc/legacy/test.py index b6f2f55e2..d48ae9871 100644 --- a/tests/regression_tests/dagmc/legacy/test.py +++ b/tests/regression_tests/dagmc/legacy/test.py @@ -45,5 +45,4 @@ def test_dagmc(): mats = openmc.Materials([u235, water]) model.materials = mats - harness = PyAPITestHarness('statepoint.5.h5', model=model) - harness.main() + model.export_to_xml() diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index af86a054d..6953186eb 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -73,7 +73,9 @@ def test_cell(capi_init): assert isinstance(cell.fill, openmc.capi.Material) cell.fill = openmc.capi.materials[1] assert str(cell) == 'Cell[0]' - + assert cell.name == "Fuel" + cell.name = "Not fuel" + assert cell.name == "Not fuel" def test_cell_temperature(capi_init): cell = openmc.capi.cells[1] @@ -122,7 +124,9 @@ def test_material(capi_init): m.set_density(0.1, 'g/cm3') assert m.density == pytest.approx(0.1) - + assert m.name == "Hot borated water" + m.name = "Not hot borated water" + assert m.name == "Not hot borated water" def test_material_add_nuclide(capi_init): m = openmc.capi.materials[3] @@ -469,3 +473,13 @@ def test_position(capi_init): pos[2] = 3.3 assert tuple(pos) == (1.3, 2.3, 3.3) + + +def test_global_bounding_box(capi_init): + expected_llc = (-0.63, -0.63, -np.inf) + expected_urc = (0.63, 0.63, np.inf) + + llc, urc = openmc.capi.global_bounding_box() + + assert tuple(llc) == expected_llc + assert tuple(urc) == expected_urc diff --git a/tests/unit_tests/test_complex_cell_capi.py b/tests/unit_tests/test_complex_cell_capi.py new file mode 100644 index 000000000..365ce4808 --- /dev/null +++ b/tests/unit_tests/test_complex_cell_capi.py @@ -0,0 +1,98 @@ +import numpy as np +import openmc.capi +import pytest + +@pytest.fixture(autouse=True) +def complex_cell(run_in_tmpdir): + + openmc.reset_auto_ids() + + model = openmc.model.Model() + + u235 = openmc.Material() + u235.set_density('g/cc', 4.5) + u235.add_nuclide("U235", 1.0) + + u238 = openmc.Material() + u238.set_density('g/cc', 4.5) + u238.add_nuclide("U238", 1.0) + + zr90 = openmc.Material() + zr90.set_density('g/cc', 2.0) + zr90.add_nuclide("Zr90", 1.0) + + n14 = openmc.Material() + n14.set_density('g/cc', 0.1) + n14.add_nuclide("N14", 1.0) + + model.materials = (u235, u238, zr90, n14) + + s1 = openmc.XPlane(x0=-10.0, boundary_type='vacuum') + s2 = openmc.XPlane(x0=-7.0) + s3 = openmc.XPlane(x0=-4.0) + s4 = openmc.XPlane(x0=4.0) + s5 = openmc.XPlane(x0=7.0) + s6 = openmc.XPlane(x0=10.0, boundary_type='vacuum') + s7 = openmc.XPlane(x0=0.0) + + s11 = openmc.YPlane(y0=-10.0, boundary_type='vacuum') + s12 = openmc.YPlane(y0=-7.0) + s13 = openmc.YPlane(y0=-4.0) + s14 = openmc.YPlane(y0=4.0) + s15 = openmc.YPlane(y0=7.0) + s16 = openmc.YPlane(y0=10.0, boundary_type='vacuum') + s17 = openmc.YPlane(y0=0.0) + + c1 = openmc.Cell(fill=u235) + c1.region = ~(-s3 | +s4 | ~(+s13 & -s14)) + + c2 = openmc.Cell(fill=u238) + c2.region = +s2 & -s5 & +s12 & -s15 & ~(+s3 & -s4 & +s13 & -s14) + + c3 = openmc.Cell(fill=zr90) + c3.region = ((+s1 & -s7 & +s17 & -s16) | (+s7 & -s6 & +s11 & -s17)) \ + & (-s2 | +s5 | -s12 | +s15) + + c4 = openmc.Cell(fill=n14) + c4.region = ((+s1 & -s7 & +s11 & -s17) | (+s7 & -s6 & +s17 & -s16)) & \ + ~(+s2 & -s5 & +s12 & -s15) + + c5 = openmc.Cell(fill=n14) + c5.region = ~(+s1 & -s6 & +s11 & -s16) + + model.geometry.root_universe = openmc.Universe() + model.geometry.root_universe.add_cells([c1, c2, c3, c4, c5]) + + model.settings.batches = 10 + model.settings.inactive = 5 + model.settings.particles = 100 + model.settings.source = openmc.Source(space=openmc.stats.Box( + [-10., -10., -1.], [10., 10., 1.])) + + model.settings.verbosity = 1 + + model.export_to_xml() + + openmc.capi.finalize() + openmc.capi.init() + + yield + + openmc.capi.finalize() + + +expected_results = ( (1, (( -4., -4., -np.inf), + ( 4., 4., np.inf))), + (2, (( -7., -7., -np.inf), + ( 7., 7., np.inf))), + (3, ((-10., -10., -np.inf), + ( 10., 10., np.inf))), + (4, ((-10., -10., -np.inf), + ( 10., 10., np.inf))), + (5, ((-np.inf, -np.inf, -np.inf), + ( np.inf, np.inf, np.inf))) ) +@pytest.mark.parametrize("cell_id,expected_box", expected_results) +def test_cell_box(cell_id, expected_box): + cell_box = openmc.capi.cells[cell_id].bounding_box + assert tuple(cell_box[0]) == expected_box[0] + assert tuple(cell_box[1]) == expected_box[1]