Merge branch 'develop' into heat-photon

include/openmc/cell.h

@@ -43,6 +43,7 @@ constexpr int32_t OP_UNION        {std::numeric_limits<int32_t>::max() - 4};
 
 class Cell;
 class Universe;
+class UniversePartitioner;
 
 namespace model {
   extern std::vector<std::unique_ptr<Cell>> cells;
@@ -65,6 +66,8 @@ public:
   //! \brief Write universe information to an HDF5 group.
   //! \param group_id An HDF5 group id.
   void to_hdf5(hid_t group_id) const;
+
+  std::unique_ptr<UniversePartitioner> partitioner_;
 };
 
 //==============================================================================
@@ -193,6 +196,36 @@ public:
 };
 #endif
 
+//==============================================================================
+//! Speeds up geometry searches by grouping cells in a search tree.
+//
+//! Currently this object only works with universes that are divided up by a
+//! bunch of z-planes.  It could be generalized to other planes, cylinders,
+//! and spheres.
+//==============================================================================
+
+class UniversePartitioner
+{
+public:
+  explicit UniversePartitioner(const Universe& univ);
+
+  //! Return the list of cells that could contain the given coordinates.
+  const std::vector<int32_t>& get_cells(Position r, Direction u) const;
+
+private:
+  //! A sorted vector of indices to surfaces that partition the universe
+  std::vector<int32_t> surfs_;
+
+  //! Vectors listing the indices of the cells that lie within each partition
+  //
+  //! There are n+1 partitions with n surfaces.  `partitions_.front()` gives the
+  //! cells that lie on the negative side of `surfs_.front()`.
+  //! `partitions_.back()` gives the cells that lie on the positive side of
+  //! `surfs_.back()`.  Otherwise, `partitions_[i]` gives cells sandwiched
+  //! between `surfs_[i-1]` and `surfs_[i]`.
+  std::vector<std::vector<int32_t>> partitions_;
+};
+
 //==============================================================================
 // Non-member functions
 //==============================================================================

include/openmc/geometry.h

@@ -35,6 +35,14 @@ struct BoundaryInfo {
   std::array<int, 3> lattice_translation {}; //!< which way lattice indices will change
 };
 
+//==============================================================================
+//! Check two distances by coincidence tolerance
+//==============================================================================
+
+inline bool coincident(double d1, double d2) {
+  return std::abs(d1 - d2) < FP_COINCIDENT;
+}
+
 //==============================================================================
 //! Check for overlapping cells at a particle's position.
 //==============================================================================

include/openmc/surface.h

@@ -180,7 +180,6 @@ public:
 
 class SurfaceXPlane : public PeriodicSurface
 {
-  double x0_;
 public:
   explicit SurfaceXPlane(pugi::xml_node surf_node);
   double evaluate(Position r) const;
@@ -190,6 +189,8 @@ public:
   bool periodic_translate(const PeriodicSurface* other, Position& r,
                           Direction& u) const;
   BoundingBox bounding_box() const;
+
+  double x0_;
 };
 
 //==============================================================================
@@ -200,7 +201,6 @@ public:
 
 class SurfaceYPlane : public PeriodicSurface
 {
-  double y0_;
 public:
   explicit SurfaceYPlane(pugi::xml_node surf_node);
   double evaluate(Position r) const;
@@ -210,6 +210,8 @@ public:
   bool periodic_translate(const PeriodicSurface* other, Position& r,
                           Direction& u) const;
   BoundingBox bounding_box() const;
+
+  double y0_;
 };
 
 //==============================================================================
@@ -220,7 +222,6 @@ public:
 
 class SurfaceZPlane : public PeriodicSurface
 {
-  double z0_;
 public:
   explicit SurfaceZPlane(pugi::xml_node surf_node);
   double evaluate(Position r) const;
@@ -230,6 +231,8 @@ public:
   bool periodic_translate(const PeriodicSurface* other, Position& r,
                           Direction& u) const;
   BoundingBox bounding_box() const;
+
+  double z0_;
 };
 
 //==============================================================================
@@ -240,7 +243,6 @@ public:
 
 class SurfacePlane : public PeriodicSurface
 {
-  double A_, B_, C_, D_;
 public:
   explicit SurfacePlane(pugi::xml_node surf_node);
   double evaluate(Position r) const;
@@ -250,6 +252,8 @@ public:
   bool periodic_translate(const PeriodicSurface* other, Position& r,
                           Direction& u) const;
   BoundingBox bounding_box() const;
+
+  double A_, B_, C_, D_;
 };
 
 //==============================================================================
@@ -261,13 +265,14 @@ public:
 
 class SurfaceXCylinder : public CSGSurface
 {
-  double y0_, z0_, radius_;
 public:
   explicit SurfaceXCylinder(pugi::xml_node surf_node);
   double evaluate(Position r) const;
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
+
+  double y0_, z0_, radius_;
 };
 
 //==============================================================================
@@ -279,13 +284,14 @@ public:
 
 class SurfaceYCylinder : public CSGSurface
 {
-  double x0_, z0_, radius_;
 public:
   explicit SurfaceYCylinder(pugi::xml_node surf_node);
   double evaluate(Position r) const;
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
+
+  double x0_, z0_, radius_;
 };
 
 //==============================================================================
@@ -297,13 +303,14 @@ public:
 
 class SurfaceZCylinder : public CSGSurface
 {
-  double x0_, y0_, radius_;
 public:
   explicit SurfaceZCylinder(pugi::xml_node surf_node);
   double evaluate(Position r) const;
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
+
+  double x0_, y0_, radius_;
 };
 
 //==============================================================================
@@ -315,13 +322,14 @@ public:
 
 class SurfaceSphere : public CSGSurface
 {
-  double x0_, y0_, z0_, radius_;
 public:
   explicit SurfaceSphere(pugi::xml_node surf_node);
   double evaluate(Position r) const;
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
+
+  double x0_, y0_, z0_, radius_;
 };
 
 //==============================================================================
@@ -333,13 +341,14 @@ public:
 
 class SurfaceXCone : public CSGSurface
 {
-  double x0_, y0_, z0_, radius_sq_;
 public:
   explicit SurfaceXCone(pugi::xml_node surf_node);
   double evaluate(Position r) const;
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
+
+  double x0_, y0_, z0_, radius_sq_;
 };
 
 //==============================================================================
@@ -351,13 +360,14 @@ public:
 
 class SurfaceYCone : public CSGSurface
 {
-  double x0_, y0_, z0_, radius_sq_;
 public:
   explicit SurfaceYCone(pugi::xml_node surf_node);
   double evaluate(Position r) const;
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
+
+  double x0_, y0_, z0_, radius_sq_;
 };
 
 //==============================================================================
@@ -369,13 +379,14 @@ public:
 
 class SurfaceZCone : public CSGSurface
 {
-  double x0_, y0_, z0_, radius_sq_;
 public:
   explicit SurfaceZCone(pugi::xml_node surf_node);
   double evaluate(Position r) const;
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
+
+  double x0_, y0_, z0_, radius_sq_;
 };
 
 //==============================================================================
@@ -386,14 +397,15 @@ public:
 
 class SurfaceQuadric : public CSGSurface
 {
-  // Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0
-  double A_, B_, C_, D_, E_, F_, G_, H_, J_, K_;
 public:
   explicit SurfaceQuadric(pugi::xml_node surf_node);
   double evaluate(Position r) const;
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
+
+  // Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0
+  double A_, B_, C_, D_, E_, F_, G_, H_, J_, K_;
 };
 
 //==============================================================================