Refactor and generalize the boundary condition implementation

CMakeLists.txt

@@ -225,6 +225,7 @@ target_compile_options(faddeeva PRIVATE ${cxxflags})
 
 list(APPEND libopenmc_SOURCES
   src/bank.cpp
+  src/boundary_condition.cpp
   src/bremsstrahlung.cpp
   src/dagmc.cpp
   src/cell.cpp

docs/source/usersguide/geometry.rst

@@ -162,22 +162,26 @@ surface. To specify a vacuum boundary condition, simply change the
 Reflective and periodic boundary conditions can be set with the strings
 'reflective' and 'periodic'. Vacuum and reflective boundary conditions can be
 applied to any type of surface. Periodic boundary conditions can be applied to
-pairs of planar surfaces. For axis-aligned planes, matching periodic surfaces
-can be determined automatically. For non-axis-aligned planes, it is necessary to
-specify pairs explicitly using the :attr:`Surface.periodic_surface` attribute as
-in the following example::
+pairs of planar surfaces. If there are only two periodic surfaces they will be
+matched automatically. Otherwise it is necessary to specify pairs explicitly
+using the :attr:`Surface.periodic_surface` attribute as in the following
+example::
 
   p1 = openmc.Plane(a=0.3, b=5.0, d=1.0, boundary_type='periodic')
   p2 = openmc.Plane(a=0.3, b=5.0, d=-1.0, boundary_type='periodic')
   p1.periodic_surface = p2
 
-Rotationally-periodic boundary conditions can be specified for a pair of
-:class:`XPlane` and :class:`YPlane`; in that case, the
-:attr:`Surface.periodic_surface` attribute must be specified manually as well.
-
-.. caution:: When using rotationally-periodic boundary conditions, your geometry
-             must be defined in the first quadrant, i.e., above the y-plane and
-             to the right of the x-plane.
+Both rotational and translational periodic boundary conditions are specified in
+the same fashion. If both planes have the same normal vector, a translational
+periodicity is assumed; rotational periodicity is assumed otherwise. Currently,
+only rotations about the :math:`z`-axis are supported.
+
+For a rotational periodic BC, the normal vectors of each surface must point
+inwards---towards the valid geometry. For example, a :class:`XPlane` and
+:class:`YPlane` would be valid for a 90-degree periodic rotation if the geometry
+lies in the first quadrant of the Cartesian grid. If the geometry instead lies
+in the fourth quadrant, the :class:`YPlane` must be replaced by a
+:class:`Plane` with the normal vector pointing in the :math:`-y` direction.
 
 .. _usersguide_cells:
 

include/openmc/boundary_condition.h

@@ -0,0 +1,118 @@
+#ifndef OPENMC_BOUNDARY_CONDITION_H
+#define OPENMC_BOUNDARY_CONDITION_H
+
+#include "openmc/position.h"
+
+namespace openmc {
+
+// Forward declare some types used in function arguments.
+class Particle;
+class Surface;
+
+//==============================================================================
+//! A class that tells particles what to do after they strike an outer boundary.
+//==============================================================================
+
+class BoundaryCondition {
+public:
+  //! Perform tracking operations for a particle that strikes the boundary.
+  //! \param p The particle that struck the boundary.  This class is not meant
+  //!   to directly modify anything about the particle, but it will do so
+  //!   indirectly by calling the particle's appropriate cross_*_bc function.
+  //! \param surf The specific surface on the boundary the particle struck.
+  virtual void
+  handle_particle(Particle& p, const Surface& surf) const = 0;
+
+  //! Return a string classification of this BC.
+  virtual std::string type() const = 0;
+};
+
+//==============================================================================
+//! A BC that kills particles, indicating they left the problem.
+//==============================================================================
+
+class VacuumBC : public BoundaryCondition {
+public:
+  void
+  handle_particle(Particle& p, const Surface& surf) const override;
+
+  std::string type() const override {return "vacuum";}
+};
+
+//==============================================================================
+//! A BC that returns particles via specular reflection.
+//==============================================================================
+
+class ReflectiveBC : public BoundaryCondition {
+public:
+  void
+  handle_particle(Particle& p, const Surface& surf) const override;
+
+  std::string type() const override {return "reflective";}
+};
+
+//==============================================================================
+//! A BC that returns particles via diffuse reflection.
+//==============================================================================
+
+class WhiteBC : public BoundaryCondition {
+public:
+  void
+  handle_particle(Particle& p, const Surface& surf) const override;
+
+  std::string type() const override {return "white";}
+};
+
+//==============================================================================
+//! A BC that moves particles to another part of the problem.
+//==============================================================================
+
+class PeriodicBC : public BoundaryCondition {
+public:
+  PeriodicBC(int i_surf, int j_surf)
+    : i_surf_(i_surf), j_surf_(j_surf)
+  {};
+
+  std::string type() const override {return "periodic";}
+
+protected:
+  int i_surf_;
+  int j_surf_;
+};
+
+//==============================================================================
+//! A BC that moves particles to another part of the problem without rotation.
+//==============================================================================
+
+class TranslationalPeriodicBC : public PeriodicBC {
+public:
+  TranslationalPeriodicBC(int i_surf, int j_surf);
+
+  void
+  handle_particle(Particle& p, const Surface& surf) const override;
+
+protected:
+  //! Vector along which incident particles will be moved
+  Position translation_;
+};
+
+//==============================================================================
+//! A BC that rotates particles about a global axis.
+//
+//! Currently only rotations about the z-axis are supported.
+//==============================================================================
+
+class RotationalPeriodicBC : public PeriodicBC {
+public:
+  RotationalPeriodicBC(int i_surf, int j_surf);
+
+  void
+  handle_particle(Particle& p, const Surface& surf) const override;
+
+protected:
+  //! Angle about the axis by which particle coordinates will be rotated
+  double angle_;
+};
+
+} // namespace openmc
+#endif // OPENMC_BOUNDARY_CONDITION_H

include/openmc/particle.h

@@ -40,6 +40,9 @@ constexpr double CACHE_INVALID {-1.0};
 // Class declarations
 //==============================================================================
 
+// Forward declare the Surface class for use in function arguments.
+class Surface;
+
 class LocalCoord {
 public:
   void rotate(const std::vector<double>& rotation);
@@ -237,6 +240,30 @@ class Particle {
   //! Cross a surface and handle boundary conditions
   void cross_surface();
 
+  //! Cross a vacuum boundary condition.
+  //
+  //! \param surf The surface (with the vacuum boundary condition) that the
+  //!   particle struck.
+  void cross_vacuum_bc(const Surface& surf);
+
+  //! Cross a reflective boundary condition.
+  //
+  //! \param surf The surface (with the reflective boundary condition) that the
+  //!   particle struck.
+  //! \param new_u The direction of the particle after reflection.
+  void cross_reflective_bc(const Surface& surf, Direction new_u);
+
+  //! Cross a periodic boundary condition.
+  //
+  //! \param surf The surface (with the periodic boundary condition) that the
+  //!   particle struck.
+  //! \param new_r The position of the particle after translation/rotation.
+  //! \param new_u The direction of the particle after translation/rotation.
+  //! \param new_surface The signed index of the surface that the particle will
+  //!   reside on after translation/rotation.
+  void cross_periodic_bc(const Surface& surf, Position new_r, Direction new_u,
+                         int new_surface);
+
   //! mark a particle as lost and create a particle restart file
   //! \param message A warning message to display
   void mark_as_lost(const char* message);

include/openmc/surface.h

@@ -10,6 +10,7 @@
 #include "hdf5.h"
 #include "pugixml.hpp"
 
+#include "openmc/boundary_condition.h"
 #include "openmc/constants.h"
 #include "openmc/particle.h"
 #include "openmc/position.h"
@@ -84,18 +85,9 @@ class Surface
 {
 public:
 
-  // Types of available boundary conditions on a surface
-  enum class BoundaryType {
-    TRANSMIT,
-    VACUUM,
-    REFLECT,
-    PERIODIC,
-    WHITE
-  };
-
-  int id_;                    //!< Unique ID
-  BoundaryType bc_;                    //!< Boundary condition
-  std::string name_;          //!< User-defined name
+  int id_; //!< Unique ID
+  std::string name_; //!< User-defined name
+  std::shared_ptr<BoundaryCondition> bc_ {nullptr}; //!< Boundary condition
 
   explicit Surface(pugi::xml_node surf_node);
   Surface();
@@ -141,7 +133,6 @@ class Surface
 
   //! Write all information needed to reconstruct the surface to an HDF5 group.
   //! \param group_id An HDF5 group id.
-  //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.
@@ -180,50 +171,21 @@ class DAGSurface : public Surface
   int32_t dag_index_;      //!< DagMC index of surface
 };
 #endif
-//==============================================================================
-//! A `Surface` that supports periodic boundary conditions.
-//!
-//! Translational periodicity is supported for the `XPlane`, `YPlane`, `ZPlane`,
-//! and `Plane` types.  Rotational periodicity is supported for
-//! `XPlane`-`YPlane` pairs.
-//==============================================================================
-
-class PeriodicSurface : public CSGSurface
-{
-public:
-  int i_periodic_{C_NONE};    //!< Index of corresponding periodic surface
-
-  explicit PeriodicSurface(pugi::xml_node surf_node);
-
-  //! Translate a particle onto this surface from a periodic partner surface.
-  //! \param other A pointer to the partner surface in this periodic BC.
-  //! \param r A point on the partner surface that will be translated onto
-  //!   this surface.
-  //! \param u A direction that will be rotated for systems with rotational
-  //!   periodicity.
-  //! \return true if this surface and its partner make a rotationally-periodic
-  //!   boundary condition.
-  virtual bool periodic_translate(const PeriodicSurface* other, Position& r,
-                                  Direction& u) const = 0;
-
-};
 
 //==============================================================================
 //! A plane perpendicular to the x-axis.
 //
 //! The plane is described by the equation \f$x - x_0 = 0\f$
 //==============================================================================
 
-class SurfaceXPlane : public PeriodicSurface
+class SurfaceXPlane : public CSGSurface
 {
 public:
   explicit SurfaceXPlane(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;
-  bool periodic_translate(const PeriodicSurface* other, Position& r,
-                          Direction& u) const;
   BoundingBox bounding_box(bool pos_side) const;
 
   double x0_;
@@ -235,16 +197,14 @@ class SurfaceXPlane : public PeriodicSurface
 //! The plane is described by the equation \f$y - y_0 = 0\f$
 //==============================================================================
 
-class SurfaceYPlane : public PeriodicSurface
+class SurfaceYPlane : public CSGSurface
 {
 public:
   explicit SurfaceYPlane(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;
-  bool periodic_translate(const PeriodicSurface* other, Position& r,
-                          Direction& u) const;
   BoundingBox bounding_box(bool pos_side) const;
 
   double y0_;
@@ -256,16 +216,14 @@ class SurfaceYPlane : public PeriodicSurface
 //! The plane is described by the equation \f$z - z_0 = 0\f$
 //==============================================================================
 
-class SurfaceZPlane : public PeriodicSurface
+class SurfaceZPlane : public CSGSurface
 {
 public:
   explicit SurfaceZPlane(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;
-  bool periodic_translate(const PeriodicSurface* other, Position& r,
-                          Direction& u) const;
   BoundingBox bounding_box(bool pos_side) const;
 
   double z0_;
@@ -277,16 +235,14 @@ class SurfaceZPlane : public PeriodicSurface
 //! The plane is described by the equation \f$A x + B y + C z - D = 0\f$
 //==============================================================================
 
-class SurfacePlane : public PeriodicSurface
+class SurfacePlane : public CSGSurface
 {
 public:
   explicit SurfacePlane(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;
-  bool periodic_translate(const PeriodicSurface* other, Position& r,
-                          Direction& u) const;
 
   double A_, B_, C_, D_;
 };