Merge branch 'develop' into refactor_write_source_bank

include/openmc/constants.h

@@ -25,7 +25,7 @@ using double_4dvec = vector<vector<vector<vector<double>>>>;
 constexpr int HDF5_VERSION[] {3, 0};
 
 // Version numbers for binary files
-constexpr array<int, 2> VERSION_STATEPOINT {18, 0};
+constexpr array<int, 2> VERSION_STATEPOINT {18, 1};
 constexpr array<int, 2> VERSION_PARTICLE_RESTART {2, 0};
 constexpr array<int, 2> VERSION_TRACK {3, 0};
 constexpr array<int, 2> VERSION_SUMMARY {6, 0};
@@ -51,6 +51,10 @@ constexpr double FP_PRECISION {1e-14};
 constexpr double FP_REL_PRECISION {1e-5};
 constexpr double FP_COINCIDENT {1e-12};
 
+// Coincidence tolerances
+constexpr double TORUS_TOL {1e-10};
+constexpr double RADIAL_MESH_TOL {1e-10};
+
 // Maximum number of random samples per history
 constexpr int MAX_SAMPLE {100000};
 

include/openmc/distribution.h

@@ -31,6 +31,39 @@ using UPtrDist = unique_ptr<Distribution>;
 //! \return Unique pointer to distribution
 UPtrDist distribution_from_xml(pugi::xml_node node);
 
+//==============================================================================
+//! A discrete distribution index (probability mass function)
+//==============================================================================
+
+class DiscreteIndex {
+public:
+  DiscreteIndex() {};
+  DiscreteIndex(pugi::xml_node node);
+  DiscreteIndex(const double* p, int n);
+
+  void assign(const double* p, int n);
+
+  //! Sample a value from the distribution
+  //! \param seed Pseudorandom number seed pointer
+  //! \return Sampled value
+  size_t sample(uint64_t* seed) const;
+
+  // Properties
+  const vector<double>& prob() const { return prob_; }
+  const vector<size_t>& alias() const { return alias_; }
+
+private:
+  vector<double> prob_; //!< Probability of accepting the uniformly sampled bin,
+                        //!< mapped to alias method table
+  vector<size_t> alias_; //!< Alias table
+
+  //! Normalize distribution so that probabilities sum to unity
+  void normalize();
+
+  //! Initialize alias tables for distribution
+  void init_alias();
+};
+
 //==============================================================================
 //! A discrete distribution (probability mass function)
 //==============================================================================
@@ -47,20 +80,13 @@ public:
 
   // Properties
   const vector<double>& x() const { return x_; }
-  const vector<double>& prob() const { return prob_; }
-  const vector<size_t>& alias() const { return alias_; }
+  const vector<double>& prob() const { return di_.prob(); }
+  const vector<size_t>& alias() const { return di_.alias(); }
 
 private:
-  vector<double> x_;    //!< Possible outcomes
-  vector<double> prob_; //!< Probability of accepting the uniformly sampled bin,
-                        //!< mapped to alias method table
-  vector<size_t> alias_; //!< Alias table
-
-  //! Normalize distribution so that probabilities sum to unity
-  void normalize();
-
-  //! Initialize alias tables for distribution
-  void init_alias(vector<double>& x, vector<double>& p);
+  vector<double> x_; //!< Possible outcomes
+  DiscreteIndex di_; //!< discrete probability distribution of
+                     //!< outcome indices
 };
 
 //==============================================================================

include/openmc/distribution_spatial.h

@@ -114,11 +114,8 @@ public:
 
 private:
   int32_t mesh_idx_ {C_NONE};
-  double total_strength_ {0.0};
-  // TODO: move to an independent class in the future that's similar
-  // to a discrete distribution without outcomes
-  std::vector<double> mesh_CDF_;
-  std::vector<double> mesh_strengths_;
+  DiscreteIndex elem_idx_dist_; //!< Distribution of
+                                //!< mesh element indices
 };
 
 //==============================================================================

include/openmc/material.h

@@ -90,6 +90,12 @@ public:
   void set_densities(
     const vector<std::string>& name, const vector<double>& density);
 
+  //! Clone the material by deep-copying all members, except for the ID,
+  //  which will get auto-assigned to the next available ID. After creating
+  //  the new material, it is added to openmc::model::materials.
+  //! \return reference to the cloned material
+  Material & clone();
+
   //----------------------------------------------------------------------------
   // Accessors
 

include/openmc/mesh.h

@@ -75,6 +75,12 @@ public:
 
   // Methods
 
+  //! Update a position to the local coordinates of the mesh
+  virtual void local_coords(Position& r) const {};
+
+  //! Return a position in the local coordinates of the mesh
+  virtual Position local_coords(const Position& r) const { return r; };
+
   //! Sample a mesh volume using a certain seed
   //
   //! \param[in] seed Seed to use for random sampling
@@ -262,6 +268,23 @@ public:
 protected:
 };
 
+class PeriodicStructuredMesh : public StructuredMesh {
+
+public:
+  PeriodicStructuredMesh() = default;
+  PeriodicStructuredMesh(pugi::xml_node node) : StructuredMesh {node} {};
+
+  void local_coords(Position& r) const override { r -= origin_; };
+
+  Position local_coords(const Position& r) const override
+  {
+    return r - origin_;
+  };
+
+  // Data members
+  Position origin_ {0.0, 0.0, 0.0}; //!< Origin of the mesh
+};
+
 //==============================================================================
 //! Tessellation of n-dimensional Euclidean space by congruent squares or cubes
 //==============================================================================
@@ -352,7 +375,7 @@ public:
   int set_grid();
 };
 
-class CylindricalMesh : public StructuredMesh {
+class CylindricalMesh : public PeriodicStructuredMesh {
 public:
   // Constructors
   CylindricalMesh() = default;
@@ -406,7 +429,7 @@ private:
   }
 };
 
-class SphericalMesh : public StructuredMesh {
+class SphericalMesh : public PeriodicStructuredMesh {
 public:
   // Constructors
   SphericalMesh() = default;

include/openmc/xml_interface.h

@@ -9,6 +9,7 @@
 #include "xtensor/xadapt.hpp"
 #include "xtensor/xarray.hpp"
 
+#include "openmc/position.h"
 #include "openmc/vector.h"
 
 namespace openmc {
@@ -49,5 +50,8 @@ xt::xarray<T> get_node_xarray(
   return xt::adapt(v, shape);
 }
 
+Position get_node_position(
+  pugi::xml_node node, const char* name, bool lowercase = false);
+
 } // namespace openmc
 #endif // OPENMC_XML_INTERFACE_H