Merge branch 'main' into spack-build-fixes

src/CMakeLists.txt

@@ -56,8 +56,8 @@ target_compile_features(openturbine_library PUBLIC cxx_std_17)
 set_target_properties(openturbine_library PROPERTIES PUBLIC_HEADER "types.hpp")
 
 # Add subdirectories for additional components of the library
-add_subdirectory(beams)
 add_subdirectory(constraints)
+add_subdirectory(elements)
 add_subdirectory(math)
 add_subdirectory(model)
 add_subdirectory(solver)

src/constraints/constraint_type.hpp

@@ -2,28 +2,43 @@
 
 namespace openturbine {
 
-/// @brief Enum class to define the type of the constraint
+/**
+ * @brief Enumeration to define the type of the constraint
+ */
 enum class ConstraintType : std::uint8_t {
-    kNone = 0,           // No constraint (default)
-    kFixedBC = 1,        // Fixed boundary condition/clamped constraint (all 6 DOFs fixed)
-    kPrescribedBC = 2,   // Prescribed boundary condition (displacement and rotation are specified
+    kNone = 0,           //< No constraint (default)
+    kFixedBC = 1,        //< Fixed boundary condition/clamped constraint (all 6 DOFs fixed)
+    kPrescribedBC = 2,   //< Prescribed boundary condition (displacement and rotation are specified
                          // i.e. all 6 DOFs are fixed)
     kRigidJoint = 3,     // Rigid constraint between two nodes (no relative motion between nodes i.e.
                          // all 6 DOFs of target node constrained)
-    kRevoluteJoint = 4,  // Target node rotates freely around a specified axis (5 DOFs constrained)
-    kRotationControl = 5,  // A rotation is specified about a given axis (1 DOF allowed about the
+    kRevoluteJoint = 4,  //< Target node rotates freely around a specified axis (5 DOFs constrained)
+    kRotationControl = 5,  //< A rotation is specified about a given axis (1 DOF allowed about the
                            // axis specified by the user - the other 5 DOFs are constrained)
 };
 
+/**
+ * @brief Returns the number of nodes used by the constraint type
+ *
+ * @param t Constraint type
+ * @return Number of nodes
+ */
 KOKKOS_INLINE_FUNCTION
 constexpr size_t GetNumberOfNodes(ConstraintType t) {
+    // Rigid joint, revolute joint, and rotation control have two nodes
     const auto has_two_nodes = t == ConstraintType::kRigidJoint ||
                                t == ConstraintType::kRevoluteJoint ||
                                t == ConstraintType::kRotationControl;
+    // Default is one node (fixed and prescribed BCs)
     return 1U + static_cast<size_t>(has_two_nodes);
 }
 
-/// @brief Returns the number of degrees of freedom used/fixed by the constraint type
+/**
+ * @brief Returns the number of degrees of freedom used/fixed by the constraint type
+ *
+ * @param type Constraint type
+ * @return Number of degrees of freedom
+ */
 KOKKOS_INLINE_FUNCTION
 constexpr size_t NumDOFsForConstraint(ConstraintType type) {
     switch (type) {
@@ -34,4 +49,5 @@ constexpr size_t NumDOFsForConstraint(ConstraintType type) {
             return static_cast<size_t>(6U);  // Default: Fixes 6 DOFs
     }
 }
+
 }  // namespace openturbine

src/dof_management/assemble_node_freedom_allocation_table.hpp

@@ -4,8 +4,8 @@
 
 #include "freedom_signature.hpp"
 
-#include "src/beams/beams.hpp"
 #include "src/constraints/constraints.hpp"
+#include "src/elements/elements.hpp"
 #include "src/state/state.hpp"
 
 namespace openturbine {
@@ -28,6 +28,18 @@ struct AssembleNodeFreedomMapTable_Beams {
     }
 };
 
+struct AssembleNodeFreedomMapTable_Masses {
+    Kokkos::View<size_t*>::const_type node_state_indices;
+    Kokkos::View<FreedomSignature*>::const_type element_freedom_signature;
+    Kokkos::View<FreedomSignature*> node_freedom_allocation_table;
+
+    KOKKOS_FUNCTION
+    void operator()(size_t i) const {
+        const auto node_index = node_state_indices(i);
+        Kokkos::atomic_or(&node_freedom_allocation_table(node_index), element_freedom_signature(i));
+    }
+};
+
 struct AssembleNodeFreedomMapTable_Constraints {
     Kokkos::View<ConstraintType*>::const_type type;
     Kokkos::View<size_t*>::const_type target_node_index;
@@ -55,12 +67,20 @@ struct AssembleNodeFreedomMapTable_Constraints {
 };
 
 inline void assemble_node_freedom_allocation_table(
-    State& state, const Beams& beams, const Constraints& constraints
+    State& state, const Elements& elements, const Constraints& constraints
 ) {
     Kokkos::parallel_for(
-        "AssembleNodeFreedomMapTable_Beams", beams.num_elems,
+        "AssembleNodeFreedomMapTable_Beams", elements.beams.num_elems,
         AssembleNodeFreedomMapTable_Beams{
-            beams.num_nodes_per_element, beams.node_state_indices, beams.element_freedom_signature,
+            elements.beams.num_nodes_per_element, elements.beams.node_state_indices,
+            elements.beams.element_freedom_signature, state.node_freedom_allocation_table
+        }
+    );
+
+    Kokkos::parallel_for(
+        "AssembleNodeFreedomMapTable_Masses", elements.masses.num_elems,
+        AssembleNodeFreedomMapTable_Masses{
+            elements.masses.state_indices, elements.masses.element_freedom_signature,
             state.node_freedom_allocation_table
         }
     );

src/dof_management/compute_node_freedom_map_table.hpp

@@ -4,7 +4,6 @@
 
 #include "freedom_signature.hpp"
 
-#include "src/beams/beams.hpp"
 #include "src/state/state.hpp"
 
 namespace openturbine {

src/dof_management/create_element_freedom_table.hpp

@@ -4,12 +4,12 @@
 
 #include "freedom_signature.hpp"
 
-#include "src/beams/beams.hpp"
+#include "src/elements/elements.hpp"
 #include "src/state/state.hpp"
 
 namespace openturbine {
 
-struct CreateElementFreedomTable {
+struct CreateElementFreedomTable_Beams {
     Kokkos::View<size_t*>::const_type num_nodes_per_element;
     Kokkos::View<size_t**>::const_type node_state_indices;
     Kokkos::View<size_t*>::const_type node_freedom_map_table;
@@ -26,12 +26,34 @@ struct CreateElementFreedomTable {
     }
 };
 
-inline void create_element_freedom_table(Beams& beams, const State& state) {
+struct CreateElementFreedomTable_Masses {
+    Kokkos::View<size_t*>::const_type node_state_indices;
+    Kokkos::View<size_t*>::const_type node_freedom_map_table;
+    Kokkos::View<size_t**> element_freedom_table;
+
+    KOKKOS_FUNCTION
+    void operator()(size_t i) const {
+        const auto node_index = node_state_indices(i);
+        for (auto k = 0U; k < 6U; ++k) {
+            element_freedom_table(i, k) = node_freedom_map_table(node_index) + k;
+        }
+    }
+};
+
+inline void create_element_freedom_table(Elements& elements, const State& state) {
+    Kokkos::parallel_for(
+        "CreateElementFreedomTable_Beams", elements.beams.num_elems,
+        CreateElementFreedomTable_Beams{
+            elements.beams.num_nodes_per_element, elements.beams.node_state_indices,
+            state.node_freedom_map_table, elements.beams.element_freedom_table
+        }
+    );
+
     Kokkos::parallel_for(
-        "Create Element Freedom Table", beams.num_elems,
-        CreateElementFreedomTable{
-            beams.num_nodes_per_element, beams.node_state_indices, state.node_freedom_map_table,
-            beams.element_freedom_table
+        "CreateElementFreedomTable_Masses", elements.masses.num_elems,
+        CreateElementFreedomTable_Masses{
+            elements.masses.state_indices, state.node_freedom_map_table,
+            elements.masses.element_freedom_table
         }
     );
 }

src/elements/CMakeLists.txt

@@ -0,0 +1,4 @@
+target_sources(openturbine_library PRIVATE)
+
+add_subdirectory(beams)
+add_subdirectory(masses)

src/beams/beam_element.hpp → src/elements/beams/beam_element.hpp

@@ -10,6 +10,13 @@
 
 namespace openturbine {
 
+/**
+ * @brief Beam element constitutes flexible beams material behavior in openturbine.
+ *
+ * @details A beam element is defined by a set of nodes and sections. Each section is defined by a
+ * 6x6 mass matrix and a 6x6 stiffness matrix. The element quadrature is used to integrate the
+ * mass and stiffness matrices along the length of the beam.
+ */
 struct BeamElement {
     std::vector<BeamNode> nodes;        // Element node positions/rotations in material frame
     std::vector<BeamSection> sections;  // Element mass/stiffness in material frame

src/beams/beams_input.hpp → src/elements/beams/beams_input.hpp

@@ -10,12 +10,13 @@
 namespace openturbine {
 
 /**
- * @brief Represents the input data for beam simulations
+ * @brief Represents the input data for creating flexible beams
  *
- * This struct encapsulates all necessary input parameters for beam simulations,
- * including the beam elements and environmental factors such as gravity.
- * It also provides utility methods for accessing and computing various
- * properties of the beam structure.
+ * This struct encapsulates all necessary input parameters for instantiating flex beams in
+ * openturbine i.e. the beam elements and environmental factors such as gravity.
+ * It also provides some utilities for computing properties such as
+ * - total number of nodes/quadrature points
+ * - maximum number of nodes/quadrature points per element etc.
  */
 struct BeamsInput {
     std::vector<BeamElement> elements;  //< Elements in the beam

src/beams/calculate_QP_position.hpp → src/elements/beams/calculate_QP_position.hpp

@@ -6,6 +6,22 @@
 
 namespace openturbine {
 
+/**
+ * @brief Functor to calculate current position and orientation at quadrature points
+ *
+ * This functor updates the current position and orientation (qp_x_) of quadrature points
+ * by combining:
+ * - Initial position (qp_x0_) with displacement (qp_u_) for translational components
+ * - Initial orientation (qp_r0_) with rotation (qp_r_) for rotational components using quaternions
+ *
+ * @param i_elem Element index
+ * @param qp_x0_ Initial quadrature point positions (num_elems x num_qps x 3)
+ * @param qp_u_ Displacements at quadrature points (num_elems x num_qps x 3)
+ * @param qp_r0_ Initial orientations as quaternions (num_elems x num_qps x 4)
+ * @param qp_r_ Rotations as quaternions (num_elems x num_qps x 4)
+ * @param qp_x_ Output current positions and orientations (num_elems x num_qps x 7)
+ *              where [0:3] = position, [3:7] = orientation quaternion
+ */
 struct CalculateQPPosition {
     size_t i_elem;
     Kokkos::View<double** [3]>::const_type qp_x0_;
@@ -33,4 +49,5 @@ struct CalculateQPPosition {
         qp_x_(i_elem, i_qp, 6) = RR0(3);
     }
 };
+
 }  // namespace openturbine

src/beams/create_beams.hpp → src/elements/beams/create_beams.hpp

@@ -1,6 +1,7 @@
 #pragma once
 
 #include "beams.hpp"
+#include "beams_input.hpp"
 #include "calculate_jacobian.hpp"
 #include "interpolate_QP_position.hpp"
 #include "interpolate_QP_rotation.hpp"

src/beams/interpolate_QP_position.hpp → src/elements/beams/interpolate_QP_position.hpp

@@ -4,18 +4,27 @@
 
 namespace openturbine {
 
+/**
+ * @brief Interpolates quadrature point positions from nodal positions using shape functions
+ *
+ * This functor computes the global positions of quadrature points by interpolating
+ * from the nodal positions using shape functions:
+ * - Takes nodal positions and rotations for each element
+ * - Uses pre-computed shape function weights for interpolation
+ * - Computes the position vector (x,y,z) for each quadrature point
+ */
 struct InterpolateQPPosition {
     Kokkos::View<size_t*>::const_type num_nodes_per_element;
     Kokkos::View<size_t*>::const_type num_qps_per_element;
-    Kokkos::View<double***>::const_type shape_interpolation;  // Num Nodes x Num Quadrature points
-    Kokkos::View<double** [7]>::const_type node_position_rotation;  // Node global position vector
-    Kokkos::View<double** [3]> qp_position;                         // quadrature point position
+    Kokkos::View<double***>::const_type shape_interpolation;  //< num_elem x num_nodes x num_qps
+    Kokkos::View<double** [7]>::const_type
+        node_position_rotation;  //< node position vector/generalized coordinates in global csys
+    Kokkos::View<double** [3]> qp_position;  //< quadrature point position - x, y, z (computed)
 
     KOKKOS_FUNCTION
     void operator()(const int i_elem) const {
         const auto num_nodes = num_nodes_per_element(i_elem);
         const auto num_qps = num_qps_per_element(i_elem);
-
         for (auto j = 0U; j < num_qps; ++j) {
             auto local_result = Kokkos::Array<double, 3>{};
             for (auto i = 0U; i < num_nodes; ++i) {

src/beams/interpolate_to_quadrature_points.hpp → src/elements/beams/interpolate_to_quadrature_points.hpp

@@ -8,26 +8,38 @@
 
 namespace openturbine {
 
+/**
+ * @brief Interpolates various quantities from nodes to quadrature points for beam elements
+ *
+ * This functor handles the interpolation of multiple quantities from nodal points to
+ * quadrature points (QPs) for beam elements, including:
+ * - Displacements (u) and their derivatives
+ * - Rotations (r) and their derivatives
+ * - Velocities (u_dot) and angular velocities (omega)
+ * - Accelerations (u_ddot) and angular accelerations (omega_dot)
+ * - Final positions (x) of quadrature points
+ */
 struct InterpolateToQuadraturePoints {
-    Kokkos::View<size_t*>::const_type num_nodes_per_element;
-    Kokkos::View<size_t*>::const_type num_qps_per_element;
-    Kokkos::View<double***>::const_type shape_interp;
-    Kokkos::View<double***>::const_type shape_deriv;
-    Kokkos::View<double**>::const_type qp_jacobian;
-    Kokkos::View<double** [7]>::const_type node_u;
-    Kokkos::View<double** [6]>::const_type node_u_dot;
-    Kokkos::View<double** [6]>::const_type node_u_ddot;
-    Kokkos::View<double** [3]>::const_type qp_x0;
-    Kokkos::View<double** [4]>::const_type qp_r0;
-    Kokkos::View<double** [3]> qp_u;
-    Kokkos::View<double** [3]> qp_uprime;
-    Kokkos::View<double** [4]> qp_r;
-    Kokkos::View<double** [4]> qp_rprime;
-    Kokkos::View<double** [3]> qp_u_dot;
-    Kokkos::View<double** [3]> qp_omega;
-    Kokkos::View<double** [3]> qp_u_ddot;
-    Kokkos::View<double** [3]> qp_omega_dot;
-    Kokkos::View<double** [7]> qp_x;
+    Kokkos::View<size_t*>::const_type num_nodes_per_element;  //< Number of nodes per element
+    Kokkos::View<size_t*>::const_type num_qps_per_element;    //< Number of QPs per element
+    Kokkos::View<double***>::const_type shape_interp;         //< shape functions at QPs
+    Kokkos::View<double***>::const_type shape_deriv;          //< shape function derivatives at QPs
+    Kokkos::View<double**>::const_type qp_jacobian;           //< Jacobian at QPs
+    Kokkos::View<double** [7]>::const_type node_u;            //< Nodal displacements
+    Kokkos::View<double** [6]>::const_type node_u_dot;        //< Nodal velocities
+    Kokkos::View<double** [6]>::const_type node_u_ddot;       //< Nodal accelerations
+    Kokkos::View<double** [3]>::const_type qp_x0;             //< Initial positions at QPs
+    Kokkos::View<double** [4]>::const_type qp_r0;             //< Initial rotations at QPs
+    // Output quantities at quadrature points
+    Kokkos::View<double** [3]> qp_u;          //< Interpolated displacements at QPs
+    Kokkos::View<double** [3]> qp_uprime;     //< Displacement derivatives at QPs
+    Kokkos::View<double** [4]> qp_r;          //< Interpolated rotations at QPs
+    Kokkos::View<double** [4]> qp_rprime;     //< Rotation derivatives at QPs
+    Kokkos::View<double** [3]> qp_u_dot;      //< Interpolated velocities at QPs
+    Kokkos::View<double** [3]> qp_omega;      //< Interpolated angular velocities at QPs
+    Kokkos::View<double** [3]> qp_u_ddot;     //< Interpolated accelerations at QPs
+    Kokkos::View<double** [3]> qp_omega_dot;  //< Interpolated angular accelerations at QPs
+    Kokkos::View<double** [7]> qp_x;          //< Final positions of quadrature points
 
     KOKKOS_FUNCTION
     void operator()(Kokkos::TeamPolicy<>::member_type member) const {