From f8bc1a37a5985b5d12fb7cba83ab3b01b1ed3313 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Mon, 12 Oct 2026 10:33:24 -0400 Subject: [PATCH] feat(equations): added rigid rotation helper This is a class which adds some helper utilities for computing the effective potential and the hessian of a body under rigid rotation. Note this is not an equation in and of itself --- .../serif/equations/rigid_rotation.hpp | 236 ++++++++++++++++++ 1 file changed, 236 insertions(+) create mode 100644 src/include/serif/equations/rigid_rotation.hpp diff --git a/src/include/serif/equations/rigid_rotation.hpp b/src/include/serif/equations/rigid_rotation.hpp new file mode 100644 index 00000000..3d2cbbad --- /dev/null +++ b/src/include/serif/equations/rigid_rotation.hpp @@ -0,0 +1,236 @@ +#pragma once + +/** + * @file rigid_rotation.hpp + * @author Emily M. Boudreaux + * + * @brief An equation implementing rigid rotation + */ + +#include +#include + +#include "serif/discretization/forms/value.hpp" +#include "serif/utils/error/errors.hpp" + +namespace serif::equations +{ + /// Rigid rotation: ω = Ωn̂, λ = r − rₐ, v = ω × λ, and Ψ = ½|v|² = ½Ω²s². + /// Here n̂ is the unit axis, rₐ a point on it, and s = |n̂ × λ| the cylindrical radius. + /// Ψ is the positive centrifugal potential: ∇Ψ points outward and Φ_eff = Φ_grav − Ψ. + + /** + * @class RigidRotation + * + * @brief A class that implements rigid rotation, providing methods to compute the centrifugal potential, its gradient, and its Hessian. + * + * @note This is not an equation in the same way that other parts of this module are; rather, this is a utility class used in other + * equations. + * + * @todo This should probably be refactored into its own submodule to clean up equations a bit + * + * Broadly this class handles a number of the details of rigid rotation. Note that this is really setup to be static physics; that is to say, + * there are a series of things you setup with the constructor (the angular velocity vector and the axis origin), then those + * are used to evaluate a series of relations. There is not however any way to update and re-evaluate the same + * object. + * + * The relations implemented are + * - Centrifugal potential: \f$\Psi = \frac{1}{2}|\omega \times \lambda|^2\f$ where \f$\lambda = r - r_a\f$ is the position relative to the axis origin. This can be further simplified to \f$\Psi = \frac{1}{2}\Omega^{2}s^{2}\f$ where \f$s\f$ is the + * cylindrical radius from the axis of rotation. + * - Centrifugal potential gradient: \f$\nabla \Psi = H\ lambda\f$ where \f$H = |\omega|^2 I - \omega \otimes \omega\f$ is the Hessian of the centrifugal potential. + */ + class RigidRotation final + { + public: + using Vector = discretization::forms::VectorValue; ///< Alias for some vector (rank 1 tensor) type + using Matrix = discretization::forms::MatrixValue; ///< Alias for some matrix (rank 2 tensor) type + + /** + * @brief Construct a RigidRotation object with the given angular velocity vector and axis origin. + * + * @param angular_velocity_vector The angular velocity vector, which must be finite. + * @param axis_origin The origin of the rotation axis, which must be finite. (default to the origin, (0, 0, 0)) + * + * @throws std::invalid_argument If the angular velocity vector or axis origin is not finite. + * + * The constructor also evaluates the Hessian of the centrifugal potential at construction time. This is defined as + * + * \f[ + * H = \nabla\nabla \Psi = |\omega|^2 I - \omega \otimes \omega + * \f] + * + * @todo We make an assumption here about the dimensionality of the space (e.g that the hessian is a 3x3 matrix). This + * is likely fine; however, if we ever want to be able to model 1D or 2D systems we will need to refactor this to be more general. + */ + explicit RigidRotation( + const Vector &angular_velocity_vector, + const Vector &axis_origin = {} + ) : + m_angular_velocity_vector(angular_velocity_vector), + m_axis_origin(axis_origin) + { + using namespace utils::errors; + + double angular_speed_squared = 0.0; + bool finite_axis_origin = true; + + // |ω|² = ωₓ² + ωᵧ² + ω_z². + for (int component = 0; component < 3; ++component) + { + finite_axis_origin = finite_axis_origin && std::isfinite(axis_origin[component]); + angular_speed_squared += angular_velocity_vector[component] * angular_velocity_vector[component]; + } + + serif_ensure(finite_axis_origin && std::isfinite(angular_speed_squared), "Rigid rotation requires finite center coordinates and finite squared angular speed; reduce overflowing parameters."); + + // H = ∇∇Ψ = |ω|²I − ω⊗ω. + for (int row = 0; row < 3; ++row) + { + for (int column = 0; column < 3; ++column) + { + m_potential_hessian[row][column] = (row == column ? angular_speed_squared : 0.0) - angular_velocity_vector[row] * angular_velocity_vector[column]; + } + } + } + + /** + * @brief Construct a RigidRotation object with a unit rotation axis and an axis origin. + * + * @param rotation_axis The rotation axis, which must be finite and nonzero. (default to the z-axis, (0, 0, 1)) + * @param axis_origin The origin of the rotation axis, which must be finite. (default to the origin, (0, 0, 0)) + * + * @throws std::invalid_argument If the norm of the rotation axis is not finite and positive. + * + * This method normalizes the rotation axis to create a unit vector. It then constructs a RigidRotation object using this unit vector and the provided axis origin. + * + * @note This constructs a new RigidRotation object does not modify the original object in place. + * + * @todo This may make sense to turn into a free-function rather than a class method. + */ + [[nodiscard]] static RigidRotation unit_rotation(Vector rotation_axis = {0.0, 0.0, 1.0}, const Vector &axis_origin = {}) + { + using namespace utils::errors; + + const double axis_norm = std::hypot(rotation_axis[0], rotation_axis[1], rotation_axis[2]); + serif_ensure(std::isfinite(axis_norm) && axis_norm > 0.0, "A rotation axis must be finite and nonzero."); + + for (double& component : rotation_axis) + { + component /= axis_norm; + } + return RigidRotation(rotation_axis, axis_origin); + } + + /// Ψ = ½|ω × λ|²; using the cross product avoids subtracting nearly equal squared lengths near the axis. + + /** + * @brief Compute the centrifugal potential at a given position. This is computed as + * \f$ \Psi = \frac{1}{2}|\omega \times \lambda|^2\f$ where \f$\lambda = r - r_a\f$ is the position relative to the axis origin. + * + * @param position The position at which to compute the centrifugal potential, which must be finite. + * @return The centrifugal potential at the given position. + * + * @todo We should add some validation for position into this + */ + [[nodiscard]] double centrifugal_potential(const Vector& position) const + { + double linear_speed_squared = 0.0; + for (int component = 0; component < 3; ++component) + { + // These are allowed to cycle through the components of the vector. + // For some component c this will look like. Note this is + // just a trick to get the appropriate indices for + // the cross product. + + /* + * | component | next | last | + * |===========|======|======| + * | 0 | 1 | 2 | + * | 1 | 2 | 0 | + * | 2 | 0 | 1 | + */ + const int next = (component + 1) % 3; + const int last = (component + 2) % 3; + + + // vᵢ = ωⱼλₖ − ωₖλⱼ, with (i, j, k) cyclic. + + // Get to this by remembering + // v = ω x λ + // so then, trivially from the cross product, v has the form + // v_x = ω_y λ_z - ω_z λ_y + // v_y = ω_z λ_x - ω_x λ_z + // v_z = ω_x λ_y - ω_y λ_x + + // The equation for the potential is + // φ = 1/2 |v|^2 + + // Since we can calculate each component of v using the cyclic permutation of the + // indices it is trivially to compute the squared magnitude of v by summing the squares + // of each component. Then we simply multiply by 1/2 (outside of the loop) + // to get the potential. + const double velocity_component = m_angular_velocity_vector[next] * (position[last] - m_axis_origin[last]) - m_angular_velocity_vector[last] * (position[next] - m_axis_origin[next]); + linear_speed_squared += velocity_component * velocity_component; + } + return 0.5 * linear_speed_squared; + } + + /** + * @brief Compute the gradient of the centrifugal potential at a given position. This is computed as + * \f$ \nabla \Psi = H \lambda\f$ where \f$H = |\omega|^2 I - \omega \otimes \omega\f$ is the Hessian of the centrifugal potential. + * + * @param position The position at which to compute the gradient of the centrifugal potential, which must be finite. + * @return The gradient of the centrifugal potential at the given position. + * + * @todo We should add some validation for position into this + */ + [[nodiscard]] Vector centrifugal_potential_gradient(const Vector& position) const + { + Vector result{}; + for (int row = 0; row < 3; ++row) + { + for (int column = 0; column < 3; ++column) + { + result[row] += m_potential_hessian[row][column] * (position[column] - m_axis_origin[column]); + } + } + return result; + } + + /** + * @brief Get the Hessian of the centrifugal potential, which is independent of position for rigid rotation. This is computed as + * \f$ H = |\omega|^2 I - \omega \otimes \omega\f$. + * + * @return The Hessian of the centrifugal potential. + */ + [[nodiscard]] const Matrix& centrifugal_potential_hessian() const + { + return m_potential_hessian; + } + + /** + * @brief Get the angular velocity vector of the rigid rotation. + * + * @return The angular velocity vector. + */ + [[nodiscard]] const Vector& angular_velocity_vector() const + { + return m_angular_velocity_vector; + } + + /** + * @brief Get the origin of the rotation axis. + * + * @return The origin of the rotation axis. + */ + [[nodiscard]] const Vector& axis_origin() const + { + return m_axis_origin; + } + + private: + Vector m_angular_velocity_vector; ///< The angular velocity vector of the rigid rotation, which must be finite. + Vector m_axis_origin; ///< The origin of the rotation axis, which must be finite. + Matrix m_potential_hessian{}; ///< The Hessian of the centrifugal potential, which is independent of position for rigid rotation and is computed as \f$ H = |\omega|^2 I - \omega \otimes \omega\f$. + }; +} // namespace serif::equations