Add temperature interface for Cell

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<double, int32_t>
+  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<double> rotation_;
 
   std::vector<int32_t> 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<double, int32_t>
-  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() {}
 };
 
 //==============================================================================

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;
 }