Enable definition of constraints from a 6-DOF (base) node to a 3-DOF (target) node

src/constraints/calculate_constraint_force.hpp

@@ -12,23 +12,18 @@ namespace openturbine {
 struct CalculateConstraintForce {
     Kokkos::View<ConstraintType*>::const_type type;
     Kokkos::View<size_t*>::const_type target_node_index;
-    Kokkos::View<double* [3][3]>::const_type axes;
-    Kokkos::View<double* [7]>::const_type inputs;
-    Kokkos::View<double* [7]>::const_type node_u;
-    Kokkos::View<double* [6]> system_residual_terms;
+    View_Nx3x3::const_type axes;
+    View_Nx7::const_type inputs;
+    View_Nx7::const_type node_u;
+    View_Nx6 system_residual_terms;
 
     KOKKOS_FUNCTION
     void operator()(const int i_constraint) const {
-        switch (type(i_constraint)) {
-            case ConstraintType::kRevoluteJoint: {
-                // Applies the torque from a revolute joint constraint to the system residual
-                CalculateRevoluteJointForce{
-                    target_node_index, axes, inputs, node_u, system_residual_terms
-                }(i_constraint);
-            } break;
-            default: {
-                // Do nothing
-            } break;
+        if (type(i_constraint) == ConstraintType::kRevoluteJoint) {
+            // Applies the torque from a revolute joint constraint to the system residual
+            CalculateRevoluteJointForce{i_constraint, target_node_index,    axes, inputs,
+                                        node_u,       system_residual_terms}();
+            return;
         }
     }
 };

src/constraints/calculate_constraint_output.hpp

@@ -10,23 +10,19 @@ namespace openturbine {
 struct CalculateConstraintOutput {
     Kokkos::View<ConstraintType*>::const_type type;
     Kokkos::View<size_t*>::const_type target_node_index;
-    Kokkos::View<double* [3][3]>::const_type axes;
-    Kokkos::View<double* [7]>::const_type node_x0;     // Initial position
-    Kokkos::View<double* [7]>::const_type node_u;      // Displacement
-    Kokkos::View<double* [6]>::const_type node_udot;   // Velocity
-    Kokkos::View<double* [6]>::const_type node_uddot;  // Acceleration
-    Kokkos::View<double* [3]> outputs;
+    View_Nx3x3::const_type axes;
+    View_Nx7::const_type node_x0;     // Initial position
+    View_Nx7::const_type node_u;      // Displacement
+    View_Nx6::const_type node_udot;   // Velocity
+    View_Nx6::const_type node_uddot;  // Acceleration
+    View_Nx3 outputs;
 
     KOKKOS_FUNCTION
     void operator()(const int i_constraint) const {
-        switch (type(i_constraint)) {
-            case ConstraintType::kRevoluteJoint: {
-                CalculateRevoluteJointOutput{target_node_index, axes,       node_x0, node_u,
-                                             node_udot,         node_uddot, outputs}(i_constraint);
-            } break;
-            default: {
-                // Do nothing
-            } break;
+        if (type(i_constraint) == ConstraintType::kRevoluteJoint) {
+            CalculateRevoluteJointOutput{i_constraint, target_node_index, axes,       node_x0,
+                                         node_u,       node_udot,         node_uddot, outputs}();
+            return;
         }
     }
 };

src/constraints/calculate_constraint_residual_gradient.hpp

@@ -16,57 +16,83 @@ namespace openturbine {
 
 struct CalculateConstraintResidualGradient {
     Kokkos::View<ConstraintType*>::const_type type;
+    Kokkos::View<Kokkos::pair<size_t, size_t>*>::const_type target_node_col_range;
     Kokkos::View<size_t*>::const_type base_node_index;
     Kokkos::View<size_t*>::const_type target_node_index;
-    Kokkos::View<double* [3]>::const_type X0_;
-    Kokkos::View<double* [3][3]>::const_type axes;
-    Kokkos::View<double* [7]>::const_type constraint_inputs;
-    Kokkos::View<double* [7]>::const_type node_u;
-    Kokkos::View<double* [6]> residual_terms;
-    Kokkos::View<double* [6][6]> base_gradient_terms;
-    Kokkos::View<double* [6][6]> target_gradient_terms;
+    View_Nx3::const_type X0_;
+    View_Nx3x3::const_type axes;
+    View_Nx7::const_type constraint_inputs;
+    View_Nx7::const_type node_u;
+    View_Nx6 residual_terms;
+    View_Nx6x6 base_gradient_terms;
+    View_Nx6x6 target_gradient_terms;
 
     KOKKOS_FUNCTION
     void operator()(const int i_constraint) const {
-        if (type(i_constraint) == ConstraintType::kFixedBC) {
-            CalculateFixedBCConstraint{target_node_index, X0_,
-                                       constraint_inputs, node_u,
-                                       residual_terms,    target_gradient_terms}(i_constraint);
-        } else if (type(i_constraint) == ConstraintType::kPrescribedBC) {
-            CalculatePrescribedBCConstraint{target_node_index, X0_,
-                                            constraint_inputs, node_u,
-                                            residual_terms,    target_gradient_terms}(i_constraint);
-        } else if (type(i_constraint) == ConstraintType::kRigidJoint) {
-            CalculateRigidJointConstraint{base_node_index,
-                                          target_node_index,
-                                          X0_,
-                                          constraint_inputs,
-                                          node_u,
-                                          residual_terms,
-                                          base_gradient_terms,
-                                          target_gradient_terms}(i_constraint);
-        } else if (type(i_constraint) == ConstraintType::kRevoluteJoint) {
-            CalculateRevoluteJointConstraint{base_node_index,
-                                             target_node_index,
-                                             X0_,
-                                             axes,
-                                             constraint_inputs,
-                                             node_u,
-                                             residual_terms,
-                                             base_gradient_terms,
-                                             target_gradient_terms}(i_constraint);
-        } else if (type(i_constraint) == ConstraintType::kRotationControl) {
-            CalculateRotationControlConstraint{base_node_index,
-                                               target_node_index,
-                                               X0_,
-                                               axes,
-                                               constraint_inputs,
-                                               node_u,
-                                               residual_terms,
-                                               base_gradient_terms,
-                                               target_gradient_terms}(i_constraint);
+        auto constraint_type = type(i_constraint);
+        if (constraint_type == ConstraintType::kFixedBC ||
+            constraint_type == ConstraintType::kFixedBC3DOFs) {
+            CalculateFixedBCConstraint{
+                i_constraint, target_node_col_range, target_node_index,     X0_, constraint_inputs,
+                node_u,       residual_terms,        target_gradient_terms,
+            }();
+            return;
+        };
+        if (constraint_type == ConstraintType::kPrescribedBC ||
+            constraint_type == ConstraintType::kPrescribedBC3DOFs) {
+            CalculatePrescribedBCConstraint{
+                i_constraint, target_node_col_range, target_node_index,     X0_, constraint_inputs,
+                node_u,       residual_terms,        target_gradient_terms,
+            }();
+            return;
+        };
+        if (constraint_type == ConstraintType::kRigidJoint ||
+            constraint_type == ConstraintType::kRigidJoint6DOFsTo3DOFs) {
+            CalculateRigidJointConstraint{
+                i_constraint,
+                target_node_col_range,
+                base_node_index,
+                target_node_index,
+                X0_,
+                constraint_inputs,
+                node_u,
+                residual_terms,
+                base_gradient_terms,
+                target_gradient_terms,
+            }();
+            return;
+        };
+        if (constraint_type == ConstraintType::kRevoluteJoint) {
+            CalculateRevoluteJointConstraint{
+                i_constraint,
+                base_node_index,
+                target_node_index,
+                X0_,
+                axes,
+                constraint_inputs,
+                node_u,
+                residual_terms,
+                base_gradient_terms,
+                target_gradient_terms,
+            }();
+            return;
+        };
+        if (constraint_type == ConstraintType::kRotationControl) {
+            CalculateRotationControlConstraint{
+                i_constraint,
+                base_node_index,
+                target_node_index,
+                X0_,
+                axes,
+                constraint_inputs,
+                node_u,
+                residual_terms,
+                base_gradient_terms,
+                target_gradient_terms,
+            }();
+            return;
         }
-    }
+    };
 };
 
 }  // namespace openturbine

src/constraints/calculate_fixed_bc_constraint.hpp

@@ -11,15 +11,17 @@
 namespace openturbine {
 
 struct CalculateFixedBCConstraint {
+    int i_constraint;
+    Kokkos::View<Kokkos::pair<size_t, size_t>*>::const_type target_node_col_range;
     Kokkos::View<size_t*>::const_type target_node_index;
-    Kokkos::View<double* [3]>::const_type X0_;
-    Kokkos::View<double* [7]>::const_type constraint_inputs;
-    Kokkos::View<double* [7]>::const_type node_u;
-    Kokkos::View<double* [6]> residual_terms;
-    Kokkos::View<double* [6][6]> target_gradient_terms;
+    View_Nx3::const_type X0_;
+    View_Nx7::const_type constraint_inputs;
+    View_Nx7::const_type node_u;
+    View_Nx6 residual_terms;
+    View_Nx6x6 target_gradient_terms;
 
     KOKKOS_FUNCTION
-    void operator()(const int i_constraint) const {
+    void operator()() const {
         const auto i_node2 = target_node_index(i_constraint);
 
         // Initial difference between nodes
@@ -43,21 +45,14 @@ struct CalculateFixedBCConstraint {
             Kokkos::Array<double, 3>{node_u(i_node2, 0), node_u(i_node2, 1), node_u(i_node2, 2)};
         const auto u2 = View_3::const_type{u2_data.data()};
 
-        // Rotation control
-        constexpr auto RCt_data = Kokkos::Array<double, 4>{1., 0., 0., 0.};
-        const auto RCt = Kokkos::View<double[4]>::const_type{RCt_data.data()};
-
         auto R1t_data = Kokkos::Array<double, 4>{};
         const auto R1t = Kokkos::View<double[4]>{R1t_data.data()};
 
         auto R1_X0_data = Kokkos::Array<double, 4>{};
         const auto R1_X0 = Kokkos::View<double[4]>{R1_X0_data.data()};
 
-        auto R2_RCt_data = Kokkos::Array<double, 4>{};
-        const auto R2_RCt = Kokkos::View<double[4]>{R2_RCt_data.data()};
-
-        auto R2_RCt_R1t_data = Kokkos::Array<double, 4>{};
-        const auto R2_RCt_R1t = Kokkos::View<double[4]>{R2_RCt_R1t_data.data()};
+        auto R2_R1t_data = Kokkos::Array<double, 4>{};
+        const auto R2_R1t = Kokkos::View<double[4]>{R2_R1t_data.data()};
 
         auto A_data = Kokkos::Array<double, 9>{};
         const auto A = View_3x3{A_data.data()};
@@ -72,6 +67,9 @@ struct CalculateFixedBCConstraint {
         // Residual Vector
         //----------------------------------------------------------------------
 
+        const auto target_node_cols =
+            target_node_col_range(i_constraint).second - target_node_col_range(i_constraint).first;
+
         // Phi(0:3) = u2 + X0 - u1 - R1*X0
         QuaternionInverse(R1, R1t);
         RotateVectorByQuaternion(R1, X0, R1_X0);
@@ -81,12 +79,13 @@ struct CalculateFixedBCConstraint {
 
         // Angular residual
         // Phi(3:6) = axial(R2*inv(RC)*inv(R1))
-        QuaternionCompose(R2, RCt, R2_RCt);
-        QuaternionCompose(R2_RCt, R1t, R2_RCt_R1t);
-        QuaternionToRotationMatrix(R2_RCt_R1t, C);
-        AxialVectorOfMatrix(C, V3);
-        for (int i = 0; i < 3; ++i) {
-            residual_terms(i_constraint, i + 3) = V3(i);
+        if (target_node_cols == 6) {
+            QuaternionCompose(R2, R1t, R2_R1t);
+            QuaternionToRotationMatrix(R2_R1t, C);
+            AxialVectorOfMatrix(C, V3);
+            for (int i = 0; i < 3; ++i) {
+                residual_terms(i_constraint, i + 3) = V3(i);
+            }
         }
 
         //----------------------------------------------------------------------
@@ -103,13 +102,15 @@ struct CalculateFixedBCConstraint {
         }
 
         // B(3:6,3:6) = AX(R1*RC*inv(R2)) = transpose(AX(R2*inv(RC)*inv(R1)))
-        AX_Matrix(C, A);
-        for (int i = 0; i < 3; ++i) {
-            for (int j = 0; j < 3; ++j) {
-                target_gradient_terms(i_constraint, i + 3, j + 3) = A(j, i);
+        if (target_node_cols == 6) {
+            AX_Matrix(C, A);
+            for (int i = 0; i < 3; ++i) {
+                for (int j = 0; j < 3; ++j) {
+                    target_gradient_terms(i_constraint, i + 3, j + 3) = A(j, i);
+                }
             }
         }
-    }
+    };
 };
 
 }  // namespace openturbine

src/constraints/calculate_prescribed_bc_constraint.hpp

@@ -11,15 +11,17 @@
 namespace openturbine {
 
 struct CalculatePrescribedBCConstraint {
+    int i_constraint;
+    Kokkos::View<Kokkos::pair<size_t, size_t>*>::const_type target_node_col_range;
     Kokkos::View<size_t*>::const_type target_node_index;
-    Kokkos::View<double* [3]>::const_type X0_;
-    Kokkos::View<double* [7]>::const_type constraint_inputs;
-    Kokkos::View<double* [7]>::const_type node_u;
-    Kokkos::View<double* [6]> residual_terms;
-    Kokkos::View<double* [6][6]> target_gradient_terms;
+    View_Nx3::const_type X0_;
+    View_Nx7::const_type constraint_inputs;
+    View_Nx7::const_type node_u;
+    View_Nx6 residual_terms;
+    View_Nx6x6 target_gradient_terms;
 
     KOKKOS_FUNCTION
-    void operator()(const int i_constraint) const {
+    void operator()() const {
         const auto i_node2 = target_node_index(i_constraint);
 
         // Initial difference between nodes
@@ -50,21 +52,14 @@ struct CalculatePrescribedBCConstraint {
             Kokkos::Array<double, 3>{node_u(i_node2, 0), node_u(i_node2, 1), node_u(i_node2, 2)};
         const auto u2 = View_3::const_type{u2_data.data()};
 
-        // Rotation control
-        constexpr auto RCt_data = Kokkos::Array<double, 4>{1., 0., 0., 0.};
-        const auto RCt = Kokkos::View<double[4]>::const_type{RCt_data.data()};
-
         auto R1t_data = Kokkos::Array<double, 4>{};
         const auto R1t = Kokkos::View<double[4]>{R1t_data.data()};
 
         auto R1_X0_data = Kokkos::Array<double, 4>{};
         const auto R1_X0 = Kokkos::View<double[4]>{R1_X0_data.data()};
 
-        auto R2_RCt_data = Kokkos::Array<double, 4>{};
-        const auto R2_RCt = Kokkos::View<double[4]>{R2_RCt_data.data()};
-
-        auto R2_RCt_R1t_data = Kokkos::Array<double, 4>{};
-        const auto R2_RCt_R1t = Kokkos::View<double[4]>{R2_RCt_R1t_data.data()};
+        auto R2_R1t_data = Kokkos::Array<double, 4>{};
+        const auto R2_R1t = Kokkos::View<double[4]>{R2_R1t_data.data()};
 
         auto A_data = Kokkos::Array<double, 9>{};
         const auto A = View_3x3{A_data.data()};
@@ -79,6 +74,9 @@ struct CalculatePrescribedBCConstraint {
         // Residual Vector
         //----------------------------------------------------------------------
 
+        const auto n_target_node_cols =
+            target_node_col_range(i_constraint).second - target_node_col_range(i_constraint).first;
+
         // Phi(0:3) = u2 + X0 - u1 - R1*X0
         QuaternionInverse(R1, R1t);
         RotateVectorByQuaternion(R1, X0, R1_X0);
@@ -87,13 +85,17 @@ struct CalculatePrescribedBCConstraint {
         }
 
         // Angular residual
-        // Phi(3:6) = axial(R2*inv(RC)*inv(R1))
-        QuaternionCompose(R2, RCt, R2_RCt);
-        QuaternionCompose(R2_RCt, R1t, R2_RCt_R1t);
-        QuaternionToRotationMatrix(R2_RCt_R1t, C);
-        AxialVectorOfMatrix(C, V3);
         for (int i = 0; i < 3; ++i) {
-            residual_terms(i_constraint, i + 3) = V3(i);
+            residual_terms(i_constraint, i + 3) = 0.;
+        }
+        // Phi(3:6) = axial(R2*inv(RC)*inv(R1))
+        if (n_target_node_cols == 6) {
+            QuaternionCompose(R2, R1t, R2_R1t);
+            QuaternionToRotationMatrix(R2_R1t, C);
+            AxialVectorOfMatrix(C, V3);
+            for (int i = 0; i < 3; ++i) {
+                residual_terms(i_constraint, i + 3) = V3(i);
+            }
         }
 
         //----------------------------------------------------------------------
@@ -109,11 +111,19 @@ struct CalculatePrescribedBCConstraint {
             target_gradient_terms(i_constraint, i, i) = 1.;
         }
 
-        // B(3:6,3:6) = AX(R1*RC*inv(R2)) = transpose(AX(R2*inv(RC)*inv(R1)))
-        AX_Matrix(C, A);
+        // B(3:6,3:6) -> initialize to 0
         for (int i = 0; i < 3; ++i) {
             for (int j = 0; j < 3; ++j) {
-                target_gradient_terms(i_constraint, i + 3, j + 3) = A(j, i);
+                target_gradient_terms(i_constraint, i + 3, j + 3) = 0.;
+            }
+        }
+        // B(3:6,3:6) = AX(R1*RC*inv(R2)) = transpose(AX(R2*inv(RC)*inv(R1)))
+        if (n_target_node_cols == 6) {
+            AX_Matrix(C, A);
+            for (int i = 0; i < 3; ++i) {
+                for (int j = 0; j < 3; ++j) {
+                    target_gradient_terms(i_constraint, i + 3, j + 3) = A(j, i);
+                }
             }
         }
     }