diff --git a/include/openmc/boundary_condition.h b/include/openmc/boundary_condition.h index af40131f1c8..0e7cc717679 100644 --- a/include/openmc/boundary_condition.h +++ b/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 diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index 7216ac89649..b3cee5673e8 100644 --- a/src/boundary_condition.cpp +++ b/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(&surf1)) { - surf1_is_xyplane = true; - } else if (const auto* ptr = dynamic_cast(&surf1)) { - surf1_is_xyplane = true; - } else if (const auto* ptr = dynamic_cast(&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(&surf2)) { - surf2_is_xyplane = true; - } else if (const auto* ptr = dynamic_cast(&surf2)) { - surf2_is_xyplane = true; - } else if (const auto* ptr = dynamic_cast(&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);