diff --git a/include/openmc/cell.h b/include/openmc/cell.h index ee12305951..afe40fe433 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -72,9 +72,66 @@ public: //! A geometry primitive that links surfaces, universes, and materials //============================================================================== -class Cell -{ +class Cell { public: + //---------------------------------------------------------------------------- + // Constructors, destructors, factory functions + + explicit Cell(pugi::xml_node cell_node); + Cell() {}; + virtual ~Cell() = default; + + //---------------------------------------------------------------------------- + // Methods + + //! \brief Determine if a cell contains the particle at a given location. + //! + //! The bounds of the cell are detemined by a logical expression involving + //! surface half-spaces. At initialization, the expression was converted + //! to RPN notation. + //! + //! The function is split into two cases, one for simple cells (those + //! involving only the intersection of half-spaces) and one for complex cells. + //! Simple cells can be evaluated with short circuit evaluation, i.e., as soon + //! as we know that one half-space is not satisfied, we can exit. This + //! provides a performance benefit for the common case. In + //! contains_complex, we evaluate the RPN expression using a stack, similar to + //! how a RPN calculator would work. + //! \param r The 3D Cartesian coordinate to check. + //! \param u A direction used to "break ties" the coordinates are very + //! close to a surface. + //! \param on_surface The signed index of a surface that the coordinate is + //! known to be on. This index takes precedence over surface sense + //! calculations. + virtual bool + contains(Position r, Direction u, int32_t on_surface) const = 0; + + //! Find the oncoming boundary of this cell. + virtual std::pair + distance(Position r, Direction u, int32_t on_surface) const = 0; + + //! Write all information needed to reconstruct the cell to an HDF5 group. + //! \param group_id An HDF5 group id. + virtual void to_hdf5(hid_t group_id) const = 0; + + //---------------------------------------------------------------------------- + // Accessors + + //! Get the temperature of a cell instance + //! \param[in] instance Instance index. If -1 is given, the temperature for + //! the first instance is returned. + //! \return Temperature in [K] + double temperature(int32_t instance = -1) const; + + //! Set the temperature of a cell instance + //! \param[in] T Temperature in [K] + //! \param[in] instance Instance index. If -1 is given, the temperature for + //! all instances is set. + void set_temperature(double T, int32_t instance = -1); + + //---------------------------------------------------------------------------- + // Data members + int32_t id_; //!< Unique ID std::string name_; //!< User-defined name int type_; //!< Material, universe, or lattice @@ -116,41 +173,6 @@ public: std::vector rotation_; std::vector offset_; //!< Distribcell offset table - - explicit Cell(pugi::xml_node cell_node); - Cell() {}; - - //! \brief Determine if a cell contains the particle at a given location. - //! - //! The bounds of the cell are detemined by a logical expression involving - //! surface half-spaces. At initialization, the expression was converted - //! to RPN notation. - //! - //! The function is split into two cases, one for simple cells (those - //! involving only the intersection of half-spaces) and one for complex cells. - //! Simple cells can be evaluated with short circuit evaluation, i.e., as soon - //! as we know that one half-space is not satisfied, we can exit. This - //! provides a performance benefit for the common case. In - //! contains_complex, we evaluate the RPN expression using a stack, similar to - //! how a RPN calculator would work. - //! \param r The 3D Cartesian coordinate to check. - //! \param u A direction used to "break ties" the coordinates are very - //! close to a surface. - //! \param on_surface The signed index of a surface that the coordinate is - //! known to be on. This index takes precedence over surface sense - //! calculations. - virtual bool - contains(Position r, Direction u, int32_t on_surface) const = 0; - - //! Find the oncoming boundary of this cell. - virtual std::pair - distance(Position r, Direction u, int32_t on_surface) const = 0; - - //! Write all information needed to reconstruct the cell to an HDF5 group. - //! @param group_id An HDF5 group id. - virtual void to_hdf5(hid_t group_id) const = 0; - - virtual ~Cell() {} }; //============================================================================== diff --git a/src/cell.cpp b/src/cell.cpp index 87b67ef856..04273eded0 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -212,6 +212,39 @@ Universe::to_hdf5(hid_t universes_group) const // Cell implementation //============================================================================== +double +Cell::temperature(int32_t instance) const +{ + if (sqrtkT_.size() < 1) { + throw std::runtime_error{"Cell temperature has not yet been set."}; + } + + if (instance >= 0) { + double sqrtkT = sqrtkT_.size() == 1 ? + sqrtkT_.at(0) : + sqrtkT_.at(instance); + return sqrtkT * sqrtkT / K_BOLTZMANN; + } else { + return sqrtkT_[0] * sqrtkT_[0] / K_BOLTZMANN; + } +} + +void +Cell::set_temperature(double T, int32_t instance) +{ + if (instance >= 0) { + sqrtkT_.at(instance) = std::sqrt(K_BOLTZMANN * T); + } else { + for (auto& T_ : sqrtkT_) { + T_ = std::sqrt(K_BOLTZMANN * T); + } + } +} + +//============================================================================== +// CSGCell implementation +//============================================================================== + CSGCell::CSGCell() {} // empty constructor CSGCell::CSGCell(pugi::xml_node cell_node) @@ -917,27 +950,18 @@ openmc_cell_set_fill(int32_t index, int type, int32_t n, extern "C" int openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance) { - if (index >= 0 && index < model::cells.size()) { - Cell& c {*model::cells[index]}; - - if (instance) { - if (*instance >= 0 && *instance < c.sqrtkT_.size()) { - c.sqrtkT_[*instance] = std::sqrt(K_BOLTZMANN * T); - } else { - strcpy(openmc_err_msg, "Distribcell instance is out of bounds."); - return OPENMC_E_OUT_OF_BOUNDS; - } - } else { - for (auto& T_ : c.sqrtkT_) { - T_ = std::sqrt(K_BOLTZMANN * T); - } - } - - } else { + if (index < 0 || index >= model::cells.size()) { strcpy(openmc_err_msg, "Index in cells array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } + int32_t instance_index = instance ? *instance : -1; + try { + model::cells[index]->set_temperature(T, instance_index); + } catch (const std::exception& e) { + set_errmsg(e.what()); + return OPENMC_E_UNASSIGNED; + } return 0; } @@ -949,25 +973,13 @@ openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T) return OPENMC_E_OUT_OF_BOUNDS; } - Cell& c {*model::cells[index]}; - - if (c.sqrtkT_.size() < 1) { - strcpy(openmc_err_msg, "Cell temperature has not yet been set."); + int32_t instance_index = instance ? *instance : -1; + try { + *T = model::cells[index]->temperature(instance_index); + } catch (const std::exception& e) { + set_errmsg(e.what()); return OPENMC_E_UNASSIGNED; } - - if (instance) { - if (*instance >= 0 && *instance < c.n_instances_) { - double sqrtkT = c.sqrtkT_.size() == 1 ? c.sqrtkT_[0] : c.sqrtkT_[*instance]; - *T = sqrtkT * sqrtkT / K_BOLTZMANN; - } else { - strcpy(openmc_err_msg, "Distribcell instance is out of bounds."); - return OPENMC_E_OUT_OF_BOUNDS; - } - } else { - *T = c.sqrtkT_[0] * c.sqrtkT_[0] / K_BOLTZMANN; - } - return 0; }