Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update constraint assembly #364

Open
wants to merge 6 commits into
base: update_tangent_application
Choose a base branch
from
Open

Update constraint assembly #364

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
92b81b9
Removed the need to do a specific matrix transpose at constraint matr…
ddement Feb 24, 2025
0420859
Moved lambda residual contribution to be constraint-wise instead of u…
ddement Feb 26, 2025
35b3155
Changed constraint variable update to happen all in stack memory
ddement Feb 27, 2025
60e3c10
clang-format
ddement Feb 27, 2025
0bd8e69
Merge branch 'update_tangent_application' into update_constraint_asse…
ddement Feb 27, 2025
e749b12
clang-format and clang-tidy
ddement Feb 28, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions src/constraints/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,11 +5,14 @@ install(FILES
calculate_constraint_output.hpp
calculate_constraint_residual_gradient.hpp
calculate_fixed_bc_constraint.hpp
calculate_fixed_bc_3DOF_constraint.hpp
calculate_prescribed_bc_constraint.hpp
calculate_prescribed_bc_3DOF_constraint.hpp
calculate_revolute_joint_constraint.hpp
calculate_revolute_joint_force.hpp
calculate_revolute_joint_output.hpp
calculate_rigid_joint_constraint.hpp
calculate_rigid_joint_3DOF_constraint.hpp
calculate_rotation_control_constraint.hpp
constraint.hpp
constraints.hpp
Expand Down
592 changes: 518 additions & 74 deletions src/constraints/calculate_constraint_residual_gradient.hpp
View file
Open in desktop

Large diffs are not rendered by default.

65 changes: 65 additions & 0 deletions src/constraints/calculate_fixed_bc_3DOF_constraint.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,65 @@
#pragma once

#include <Kokkos_Core.hpp>

#include "constraints.hpp"
#include "math/matrix_operations.hpp"
#include "math/quaternion_operations.hpp"
#include "math/vector_operations.hpp"

namespace openturbine {

KOKKOS_FUNCTION
inline void CalculateFixedBC3DOFConstraint(
const Kokkos::View<double[3]>::const_type& X0, const Kokkos::View<double[7]>::const_type& node_u,
const Kokkos::View<double[6]>& residual_terms,
const Kokkos::View<double[6][6]>& target_gradient_terms
) {
constexpr auto u1_data = Kokkos::Array<double, 3>{0., 0., 0.};
constexpr auto R1_data = Kokkos::Array<double, 4>{1., 0., 0., 0.};
const auto u2_data = Kokkos::Array<double, 3>{node_u(0), node_u(1), node_u(2)};
const auto R2_data = Kokkos::Array<double, 4>{node_u(3), node_u(4), node_u(5), node_u(6)};
auto R1t_data = Kokkos::Array<double, 4>{};
auto R1_X0_data = Kokkos::Array<double, 4>{};
auto R2_R1t_data = Kokkos::Array<double, 4>{};
auto A_data = Kokkos::Array<double, 9>{};
auto C_data = Kokkos::Array<double, 9>{};
auto V3_data = Kokkos::Array<double, 3>{};

const auto u1 = Kokkos::View<double[3]>::const_type{u1_data.data()};
const auto R1 = Kokkos::View<double[4]>::const_type{R1_data.data()};
const auto u2 = Kokkos::View<double[3]>::const_type{u2_data.data()};
const auto R2 = Kokkos::View<double[4]>::const_type{R2_data.data()};
const auto R1t = Kokkos::View<double[4]>{R1t_data.data()};
const auto R1_X0 = Kokkos::View<double[4]>{R1_X0_data.data()};
const auto R2_R1t = Kokkos::View<double[4]>{R2_R1t_data.data()};
const auto A = Kokkos::View<double[3][3]>{A_data.data()};
const auto C = Kokkos::View<double[3][3]>{C_data.data()};
const auto V3 = Kokkos::View<double[3]>{V3_data.data()};

//----------------------------------------------------------------------
// Residual Vector
//----------------------------------------------------------------------

// Phi(0:3) = u2 + X0 - u1 - R1*X0
QuaternionInverse(R1, R1t);
RotateVectorByQuaternion(R1, X0, R1_X0);
for (int i = 0; i < 3; ++i) {
residual_terms(i) = u2(i) + X0(i) - u1(i) - R1_X0(i);
}

//----------------------------------------------------------------------
// Constraint Gradient Matrix
//----------------------------------------------------------------------

//---------------------------------
// Target Node
//---------------------------------

// B(0:3,0:3) = I
for (int i = 0; i < 3; ++i) {
target_gradient_terms(i, i) = 1.;
}
}

} // namespace openturbine
170 changes: 68 additions & 102 deletions src/constraints/calculate_fixed_bc_constraint.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,107 +9,73 @@

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;
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 {
const auto i_node2 = target_node_index(i_constraint);

// Initial difference between nodes
const auto X0_data = Kokkos::Array<double, 3>{
X0_(i_constraint, 0), X0_(i_constraint, 1), X0_(i_constraint, 2)
};
const auto X0 = View_3::const_type{X0_data.data()};

// Base node displacement
constexpr auto u1_data = Kokkos::Array<double, 3>{0., 0., 0.};
const auto u1 = View_3::const_type{u1_data.data()};
constexpr auto R1_data = Kokkos::Array<double, 4>{1., 0., 0., 0.};
const auto R1 = Kokkos::View<double[4]>::const_type{R1_data.data()};

// Target node displacement
const auto R2_data = Kokkos::Array<double, 4>{
node_u(i_node2, 3), node_u(i_node2, 4), node_u(i_node2, 5), node_u(i_node2, 6)
};
const auto R2 = Kokkos::View<double[4]>::const_type{R2_data.data()};
const auto u2_data =
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()};

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_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()};

auto C_data = Kokkos::Array<double, 9>{};
const auto C = View_3x3{C_data.data()};

auto V3_data = Kokkos::Array<double, 3>{};
const auto V3 = View_3{V3_data.data()};

//----------------------------------------------------------------------
// 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);
for (int i = 0; i < 3; ++i) {
residual_terms(i_constraint, i) = u2(i) + X0(i) - u1(i) - R1_X0(i);
}

// Angular residual
// Phi(3:6) = axial(R2*inv(RC)*inv(R1))
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);
}
KOKKOS_FUNCTION
inline void CalculateFixedBCConstraint(
const Kokkos::View<double[3]>::const_type& X0, const Kokkos::View<double[7]>::const_type& node_u,
const Kokkos::View<double[6]>& residual_terms,
const Kokkos::View<double[6][6]>& target_gradient_terms
) {
constexpr auto u1_data = Kokkos::Array<double, 3>{0., 0., 0.};
constexpr auto R1_data = Kokkos::Array<double, 4>{1., 0., 0., 0.};
const auto u2_data = Kokkos::Array<double, 3>{node_u(0), node_u(1), node_u(2)};
const auto R2_data = Kokkos::Array<double, 4>{node_u(3), node_u(4), node_u(5), node_u(6)};
auto R1t_data = Kokkos::Array<double, 4>{};
auto R1_X0_data = Kokkos::Array<double, 4>{};
auto R2_R1t_data = Kokkos::Array<double, 4>{};
auto A_data = Kokkos::Array<double, 9>{};
auto C_data = Kokkos::Array<double, 9>{};
auto V3_data = Kokkos::Array<double, 3>{};

const auto u1 = Kokkos::View<double[3]>::const_type{u1_data.data()};
const auto R1 = Kokkos::View<double[4]>::const_type{R1_data.data()};
const auto u2 = Kokkos::View<double[3]>::const_type{u2_data.data()};
const auto R2 = Kokkos::View<double[4]>::const_type{R2_data.data()};
const auto R1t = Kokkos::View<double[4]>{R1t_data.data()};
const auto R1_X0 = Kokkos::View<double[4]>{R1_X0_data.data()};
const auto R2_R1t = Kokkos::View<double[4]>{R2_R1t_data.data()};
const auto A = Kokkos::View<double[3][3]>{A_data.data()};
const auto C = Kokkos::View<double[3][3]>{C_data.data()};
const auto V3 = Kokkos::View<double[3]>{V3_data.data()};

//----------------------------------------------------------------------
// Residual Vector
//----------------------------------------------------------------------

// Phi(0:3) = u2 + X0 - u1 - R1*X0
QuaternionInverse(R1, R1t);
RotateVectorByQuaternion(R1, X0, R1_X0);
for (int i = 0; i < 3; ++i) {
residual_terms(i) = u2(i) + X0(i) - u1(i) - R1_X0(i);
}

// Angular residual
// Phi(3:6) = axial(R2*inv(RC)*inv(R1))
QuaternionCompose(R2, R1t, R2_R1t);
QuaternionToRotationMatrix(R2_R1t, C);
AxialVectorOfMatrix(C, V3);
for (int i = 0; i < 3; ++i) {
residual_terms(i + 3) = V3(i);
}

//----------------------------------------------------------------------
// Constraint Gradient Matrix
//----------------------------------------------------------------------

//---------------------------------
// Target Node
//---------------------------------

// B(0:3,0:3) = I
for (int i = 0; i < 3; ++i) {
target_gradient_terms(i, i) = 1.;
}

// 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 + 3, j + 3) = A(j, i);
}

//----------------------------------------------------------------------
// Constraint Gradient Matrix
//----------------------------------------------------------------------

//---------------------------------
// Target Node
//---------------------------------

// B(0:3,0:3) = I
for (int i = 0; i < 3; ++i) {
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)))
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
65 changes: 65 additions & 0 deletions src/constraints/calculate_prescribed_bc_3DOF_constraint.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,65 @@
#pragma once

#include <Kokkos_Core.hpp>

#include "constraints.hpp"
#include "math/matrix_operations.hpp"
#include "math/quaternion_operations.hpp"
#include "math/vector_operations.hpp"

namespace openturbine {

KOKKOS_FUNCTION
inline void CalculatePrescribedBC3DOFConstraint(
const Kokkos::View<double[3]>::const_type& X0, const Kokkos::View<double[7]>::const_type& inputs,
const Kokkos::View<double[7]>::const_type& node_u, const Kokkos::View<double[6]>& residual_terms,
const Kokkos::View<double[6][6]>& target_gradient_terms
) {
const auto u1_data = Kokkos::Array<double, 3>{inputs(0), inputs(1), inputs(2)};
const auto R1_data = Kokkos::Array<double, 4>{inputs(3), inputs(4), inputs(5), inputs(6)};
const auto u2_data = Kokkos::Array<double, 3>{node_u(0), node_u(1), node_u(2)};
const auto R2_data = Kokkos::Array<double, 4>{node_u(3), node_u(4), node_u(5), node_u(6)};
auto R1t_data = Kokkos::Array<double, 4>{};
auto R1_X0_data = Kokkos::Array<double, 4>{};
auto R2_R1t_data = Kokkos::Array<double, 4>{};
auto A_data = Kokkos::Array<double, 9>{};
auto C_data = Kokkos::Array<double, 9>{};
auto V3_data = Kokkos::Array<double, 3>{};

const auto u1 = Kokkos::View<double[3]>::const_type{u1_data.data()};
const auto R1 = Kokkos::View<double[4]>::const_type{R1_data.data()};
const auto u2 = Kokkos::View<double[3]>::const_type{u2_data.data()};
const auto R2 = Kokkos::View<double[4]>::const_type{R2_data.data()};
const auto R1t = Kokkos::View<double[4]>{R1t_data.data()};
const auto R1_X0 = Kokkos::View<double[4]>{R1_X0_data.data()};
const auto R2_R1t = Kokkos::View<double[4]>{R2_R1t_data.data()};
const auto A = Kokkos::View<double[3][3]>{A_data.data()};
const auto C = Kokkos::View<double[3][3]>{C_data.data()};
const auto V3 = Kokkos::View<double[3]>{V3_data.data()};

//----------------------------------------------------------------------
// Residual Vector
//----------------------------------------------------------------------

// Phi(0:3) = u2 + X0 - u1 - R1*X0
QuaternionInverse(R1, R1t);
RotateVectorByQuaternion(R1, X0, R1_X0);
for (int i = 0; i < 3; ++i) {
residual_terms(i) = u2(i) + X0(i) - u1(i) - R1_X0(i);
}

//----------------------------------------------------------------------
// Constraint Gradient Matrix
//----------------------------------------------------------------------

//---------------------------------
// Target Node
//---------------------------------

// B(0:3,0:3) = I
for (int i = 0; i < 3; ++i) {
target_gradient_terms(i, i) = 1.;
}
}

} // namespace openturbine
Loading