allow RotationalPeriodicBC around x,y, or z axis

Author: lewisgross1296

include/openmc/boundary_condition.h

@@ -138,18 +138,23 @@ class TranslationalPeriodicBC : public PeriodicBC {
 //==============================================================================
 //! A BC that rotates particles about a global axis.
 //
-//! Currently only rotations about the z-axis are supported.
+//! Only rotations about the x, y, and z axes are supported.
 //==============================================================================
 
 class RotationalPeriodicBC : public PeriodicBC {
 public:
   RotationalPeriodicBC(int i_surf, int j_surf);
-
+  float compute_periodic_rotation(float rise_1, float run_1, float rise_2, float run_2) const;
   void handle_particle(Particle& p, const Surface& surf) const override;
 
 protected:
   //! Angle about the axis by which particle coordinates will be rotated
   double angle_;
+  enum PeriodicAxis {x, y, z};
+  //! Ensure that choice of axes is right handed. axis_1_idx_ corresponds to the independent axis and axis_2_idx_ corresponds to the dependent axis in the 2D plane  perpendicular to the planes' axis of rotation
+  int zero_axis_idx;
+  int axis_1_idx_; 
+  int axis_2_idx_;
 };
 
 } // namespace openmc

src/boundary_condition.cpp

@@ -158,56 +158,46 @@ void TranslationalPeriodicBC::handle_particle(
 // RotationalPeriodicBC implementation
 //==============================================================================
 
-RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf)
+RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf, PeriodicAxis axis)
   : PeriodicBC(i_surf, j_surf)
 {
   Surface& surf1 {*model::surfaces[i_surf_]};
   Surface& surf2 {*model::surfaces[j_surf_]};
 
-  // Check the type of the first surface
-  bool surf1_is_xyplane;
-  if (const auto* ptr = dynamic_cast<const SurfaceXPlane*>(&surf1)) {
-    surf1_is_xyplane = true;
-  } else if (const auto* ptr = dynamic_cast<const SurfaceYPlane*>(&surf1)) {
-    surf1_is_xyplane = true;
-  } else if (const auto* ptr = dynamic_cast<const SurfacePlane*>(&surf1)) {
-    surf1_is_xyplane = false;
-  } else {
-    throw std::invalid_argument(fmt::format(
-      "Surface {} is an invalid type for "
-      "rotational periodic BCs. Only x-planes, y-planes, or general planes "
-      "(that are perpendicular to z) are supported for these BCs.",
-      surf1.id_));
+  // below convention for right handed coordinate system
+  switch(axis) {
+    case x:
+      zero_axis_idx = 0; // x component of plane must be zero
+      axis_1_idx_ = 1;   // y component independent
+      axis_2_idx_ = 2;   // z component dependent
+      break;
+    case y:
+      zero_axis_idx = 1; // y component of plane must be zero
+      axis_1_idx_ = 2;   // z component independent
+      axis_2_idx_ = 0;   // x component dependent
+      break;
+    case z:
+      zero_axis_idx = 2; // z component of plane must be zero
+      axis_1_idx_ = 0;   // x component independent
+      axis_2_idx_ = 1;   // y component dependent
+      break;
+    default:
+      throw std::invalid_argument(fmt::format("You've specified an axis {} that is not x, y, or z."),axis)
   }
 
-  // Check the type of the second surface
-  bool surf2_is_xyplane;
-  if (const auto* ptr = dynamic_cast<const SurfaceXPlane*>(&surf2)) {
-    surf2_is_xyplane = true;
-  } else if (const auto* ptr = dynamic_cast<const SurfaceYPlane*>(&surf2)) {
-    surf2_is_xyplane = true;
-  } else if (const auto* ptr = dynamic_cast<const SurfacePlane*>(&surf2)) {
-    surf2_is_xyplane = false;
-  } else {
-    throw std::invalid_argument(fmt::format(
-      "Surface {} is an invalid type for "
-      "rotational periodic BCs. Only x-planes, y-planes, or general planes "
-      "(that are perpendicular to z) are supported for these BCs.",
-      surf2.id_));
-  }
 
   // Compute the surface normal vectors and make sure they are perpendicular
-  // to the z-axis
+  // to the correct axis
   Direction norm1 = surf1.normal({0, 0, 0});
   Direction norm2 = surf2.normal({0, 0, 0});
-  if (std::abs(norm1.z) > FP_PRECISION) {
+  if (std::abs(norm1[zero_axis_idx]) > FP_PRECISION) {
     throw std::invalid_argument(fmt::format(
       "Rotational periodic BCs are only "
       "supported for rotations about the z-axis, but surface {} is not "
       "perpendicular to the z-axis.",
       surf1.id_));
   }
-  if (std::abs(norm2.z) > FP_PRECISION) {
+  if (std::abs(norm2[zero_axis_idx]) > FP_PRECISION) {
     throw std::invalid_argument(fmt::format(
       "Rotational periodic BCs are only "
       "supported for rotations about the z-axis, but surface {} is not "
@@ -231,15 +221,7 @@ RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf)
       surf2.id_));
   }
 
-  // Compute the BC rotation angle.  Here it is assumed that both surface
-  // normal vectors point inwards---towards the valid geometry region.
-  // Consequently, the rotation angle is not the difference between the two
-  // normals, but is instead the difference between one normal and one
-  // anti-normal.  (An incident ray on one surface must be an outgoing ray on
-  // the other surface after rotation hence the anti-normal.)
-  double theta1 = std::atan2(norm1.y, norm1.x);
-  double theta2 = std::atan2(norm2.y, norm2.x) + PI;
-  angle_ = theta2 - theta1;
+  angle_ = compute_periodic_rotation(norm1[axis_2_idx_], norm1[axis_1_idx_], norm2[axis_2_idx_], norm2[axis_1_idx_])
 
   // Warn the user if the angle does not evenly divide a circle
   double rem = std::abs(std::remainder((2 * PI / angle_), 1.0));
@@ -251,6 +233,20 @@ RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf)
   }
 }
 
+float RotationalPeriodicBC::compute_periodic_rotation(
+  float rise_1, float run_1, float rise_2, float run_2) const
+{
+  // Compute the BC rotation angle.  Here it is assumed that both surface
+  // normal vectors point inwards---towards the valid geometry region.
+  // Consequently, the rotation angle is not the difference between the two
+  // normals, but is instead the difference between one normal and one
+  // anti-normal.  (An incident ray on one surface must be an outgoing ray on
+  // the other surface after rotation hence the anti-normal.)
+  double theta1 = std::atan2(rise_1, run_1);
+  double theta2 = std::atan2(rise_2, run_2) + PI;
+  return theta2 - theta1;
+}
+
 void RotationalPeriodicBC::handle_particle(
   Particle& p, const Surface& surf) const
 {
@@ -279,9 +275,9 @@ void RotationalPeriodicBC::handle_particle(
   double cos_theta = std::cos(theta);
   double sin_theta = std::sin(theta);
   Position new_r = {
-    cos_theta * r.x - sin_theta * r.y, sin_theta * r.x + cos_theta * r.y, r.z};
+    cos_theta * r[axis_1_idx_] - sin_theta * r[axis_2_idx_], sin_theta * r[axis_1_idx_] + cos_theta * r[axis_2_idx_], r[zero_axis_idx]};
   Direction new_u = {
-    cos_theta * u.x - sin_theta * u.y, sin_theta * u.x + cos_theta * u.y, u.z};
+    cos_theta * u[axis_1_idx_] - sin_theta * u[axis_2_idx_], sin_theta * u[axis_1_idx_] + cos_theta * u[axis_2_idx_], u[zero_axis_idx]};
 
   // Handle the effects of the surface albedo on the particle's weight.
   BoundaryCondition::handle_albedo(p, surf);