Move periodic translation from Surface class to BC

include/openmc/boundary_condition.h

@@ -0,0 +1,83 @@
+#ifndef OPENMC_BOUNDARY_CONDITION_H
+#define OPENMC_BOUNDARY_CONDITION_H
+
+#include "openmc/position.h"
+
+namespace openmc {
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+//class BoundaryCondition;
+//
+//namespace model {
+//  extern std::vector<std::unique_ptr<BoundaryCondition>> boundary_conditions;
+//} // namespace model
+
+class Particle;
+class Surface;
+
+class BoundaryCondition
+{
+public:
+  virtual void
+  handle_particle(Particle& p, const Surface& surf) const = 0;
+};
+
+class VacuumBC : public BoundaryCondition
+{
+public:
+  void
+  handle_particle(Particle& p, const Surface& surf) const override;
+};
+
+class ReflectiveBC : public BoundaryCondition
+{
+public:
+  void
+  handle_particle(Particle& p, const Surface& surf) const override;
+};
+
+class WhiteBC : public BoundaryCondition
+{
+public:
+  void
+  handle_particle(Particle& p, const Surface& surf) const override;
+};
+
+class PeriodicBC : public BoundaryCondition
+{
+public:
+  PeriodicBC(int i_surf, int j_surf)
+    : i_surf_(i_surf), j_surf_(j_surf)
+  {};
+
+protected:
+  int i_surf_;
+  int j_surf_;
+};
+
+class TranslationalPeriodicBC : public PeriodicBC
+{
+public:
+  TranslationalPeriodicBC(int i_surf, int j_surf);
+
+  void
+  handle_particle(Particle& p, const Surface& surf) const override;
+
+protected:
+  Position translation_;
+};
+
+class RotationalPeriodicBC : public PeriodicBC
+{
+public:
+  RotationalPeriodicBC(int i_surf, int j_surf);
+
+  void
+  handle_particle(Particle& p, const Surface& surf) const override;
+};
+
+} // namespace openmc
+#endif // OPENMC_BOUNDARY_CONDITION_H

src/boundary_condition.cpp

@@ -1,4 +1,9 @@
 #include "openmc/boundary_condition.h"
+
+#include <exception>
+
+#include <fmt/core.h>
+
 #include "openmc/surface.h"
 
 namespace openmc {
@@ -27,8 +32,82 @@ WhiteBC::handle_particle(Particle& p, const Surface& surf) const
   p.cross_reflective_bc(surf, u);
 }
 
+TranslationalPeriodicBC::TranslationalPeriodicBC(int i_surf, int j_surf)
+  : PeriodicBC(i_surf, j_surf)
+{
+  Surface& surf1 {*model::surfaces[i_surf_]};
+  Surface& surf2 {*model::surfaces[j_surf_]};
+
+  // The following blocks will resolve the type of each surface and compute the
+  // appropriate translation vector
+
+  // Check for a pair of x-planes
+  if (const auto* xplane1 = dynamic_cast<const SurfaceXPlane*>(&surf1)) {
+    if (const auto* xplane2 = dynamic_cast<const SurfaceXPlane*>(&surf2)) {
+      translation_ = {xplane2->x0_ - xplane1->x0_, 0, 0};
+    } else {
+      throw std::invalid_argument(fmt::format("Invalid pair of periodic "
+        "surfaces ({} and {}). For a translational periodic BC, both surfaces "
+        "must be of the same type (e.g. both x-plane).", surf1.id_, surf2.id_));
+    }
+
+  // Check for a pair of y-planes
+  } else if (const auto* yplane1 = dynamic_cast<const SurfaceYPlane*>(&surf1)) {
+    if (const auto* yplane2 = dynamic_cast<const SurfaceYPlane*>(&surf2)) {
+      translation_ = {0, yplane2->y0_ - yplane1->y0_, 0};
+    } else {
+      throw std::invalid_argument(fmt::format("Invalid pair of periodic "
+        "surfaces ({} and {}). For a translational periodic BC, both surfaces "
+        "must be of the same type (e.g. both x-plane).", surf1.id_, surf2.id_));
+    }
+
+  // Check for a pair of z-planes
+  } else if (const auto* zplane1 = dynamic_cast<const SurfaceZPlane*>(&surf1)) {
+    if (const auto* zplane2 = dynamic_cast<const SurfaceZPlane*>(&surf2)) {
+      translation_ = {0, 0, zplane2->z0_ - zplane1->z0_};
+    } else {
+      throw std::invalid_argument(fmt::format("Invalid pair of periodic "
+        "surfaces ({} and {}). For a translational periodic BC, both surfaces "
+        "must be of the same type (e.g. both x-plane).", surf1.id_, surf2.id_));
+    }
+  }
+
+  // TODO: Check for a pair of general planes
+
+}
+
 void
-PeriodicBC::handle_particle(Particle& p, const Surface& surf) const
+TranslationalPeriodicBC::handle_particle(Particle& p, const Surface& surf) const
+{
+  // TODO: off-by-one on surface indices throughout this function.
+  int i_particle_surf = std::abs(p.surface_) - 1;
+
+  // Figure out which of the two BC surfaces were struck then find the
+  // particle's new location and surface.
+  Position new_r;
+  int new_surface;
+  if (i_particle_surf == i_surf_) {
+    new_r = p.r() + translation_;
+    new_surface = p.surface_ > 0 ? j_surf_ + 1 : -(j_surf_ + 1);
+  } else if (i_particle_surf == j_surf_) {
+    new_r = p.r() - translation_;
+    new_surface = p.surface_ > 0 ? i_surf_ + 1 : -(i_surf_ + 1);
+  } else {
+    throw std::runtime_error("Called BoundaryCondition::handle_particle after "
+      "hitting a surface, but that surface is not recognized by the BC.");
+  }
+
+  // Pass the new location and surface to the particle.
+  p.cross_periodic_bc(surf, new_r, p.u(), new_surface);
+}
+
+RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf)
+  : PeriodicBC(i_surf, j_surf)
+{
+}
+
+void
+RotationalPeriodicBC::handle_particle(Particle& p, const Surface& surf) const
 {
   // Get a pointer of the first surface downcast to PeriodicSurface
   auto surf_p = dynamic_cast<const PeriodicSurface*>(&surf);
@@ -38,8 +117,6 @@ PeriodicBC::handle_particle(Particle& p, const Surface& surf) const
     surf.id_ == model::surfaces[i_surf_]->id_ ?
     dynamic_cast<const PeriodicSurface*>(model::surfaces[j_surf_].get()) :
     dynamic_cast<const PeriodicSurface*>(model::surfaces[i_surf_].get());
-  //auto other = dynamic_cast<const PeriodicSurface*>(
-  //  model::surfaces[surf_p->i_periodic_].get());
 
   // Compute the new particle location and direction
   Position r {p.r()};

src/surface.cpp

@@ -1249,8 +1249,19 @@ void read_surfaces(pugi::xml_node node)
     int j_surf = model::surface_map[periodic_pair.second];
     Surface& surf1 {*model::surfaces[i_surf]};
     Surface& surf2 {*model::surfaces[j_surf]};
-    surf1.new_bc_ = std::make_shared<PeriodicBC>(i_surf, j_surf);
-    surf2.new_bc_ = surf1.new_bc_;
+
+    Direction norm1 = surf1.normal({0, 0, 0});
+    Direction norm2 = surf2.normal({0, 0, 0});
+    double dot_prod = norm1.dot(norm2);
+
+    if (std::abs(1.0 - dot_prod) < FP_PRECISION) {
+      surf1.new_bc_ = std::make_shared<TranslationalPeriodicBC>(i_surf, j_surf);
+      surf2.new_bc_ = surf1.new_bc_;
+    } else {
+      surf1.new_bc_ = std::make_shared<RotationalPeriodicBC>(i_surf, j_surf);
+      surf2.new_bc_ = surf1.new_bc_;
+    }
+
     dynamic_cast<PeriodicSurface*>(&surf1)->i_periodic_ = j_surf;
     dynamic_cast<PeriodicSurface*>(&surf2)->i_periodic_ = i_surf;
   }