Skip to content

Commit

Permalink
Merge pull request #237 from Exawind/qp_to_elem_shape
Browse files Browse the repository at this point in the history 
Changes shape, shape derivative, weights, and jacobians to use per-element Views
  • Loading branch information
@ddement
ddement authored Aug 13, 2024
2 parents 5840474 + 7e3adb0 commit afea395
Show file tree
Hide file tree
Showing 55 changed files with 1,325 additions and 1,508 deletions.
157 changes: 79 additions & 78 deletions src/restruct_poc/beams/beams.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,46 +47,47 @@ struct Beams {
View_Nx6 node_FX; // External forces

// Quadrature point data
View_N qp_weight; // Integration weights
View_N qp_jacobian; // Jacobian vector
View_Nx6x6 qp_Mstar; // Mass matrix in material frame
View_Nx6x6 qp_Cstar; // Stiffness matrix in material frame
View_Nx3 qp_x0; // Initial position
View_Nx3 qp_x0_prime; // Initial position derivative
View_Nx4 qp_r0; // Initial rotation
View_Nx3 qp_u; // State: translation displacement
View_Nx3 qp_u_prime; // State: translation displacement derivative
View_Nx3 qp_u_dot; // State: translation velocity
View_Nx3 qp_u_ddot; // State: translation acceleration
View_Nx4 qp_r; // State: rotation
View_Nx4 qp_r_prime; // State: rotation derivative
View_Nx3 qp_omega; // State: angular velocity
View_Nx3 qp_omega_dot; // State: position/rotation
View_Nx3x4 qp_E; // Quaternion derivative
View_Nx3x3 qp_eta_tilde; //
View_Nx3x3 qp_omega_tilde; //
View_Nx3x3 qp_omega_dot_tilde; //
View_Nx3x3 qp_x0pupss; //
View_Nx3x3 qp_M_tilde; //
View_Nx3x3 qp_N_tilde; //
View_Nx3 qp_eta; //
View_Nx3x3 qp_rho; //
View_Nx6 qp_strain; // Strain
View_Nx6 qp_Fc; // Elastic force
View_Nx6 qp_Fd; // Elastic force
View_Nx6 qp_Fi; // Inertial force
View_Nx6 qp_Fg; // Gravity force
View_Nx6x6 qp_RR0; // Global rotation
View_Nx6x6 qp_Muu; // Mass in global frame
View_Nx6x6 qp_Cuu; // Stiffness in global frame
View_Nx6x6 qp_Ouu; // Linearization matrices
View_Nx6x6 qp_Puu; // Linearization matrices
View_Nx6x6 qp_Quu; // Linearization matrices
View_Nx6x6 qp_Guu; // Linearization matrices
View_Nx6x6 qp_Kuu; // Linearization matrices
View_NxN qp_weight; // Integration weights
View_NxN qp_jacobian; // Jacobian vector
Kokkos::View<double** [6][6]> qp_Mstar; // Mass matrix in material frame
Kokkos::View<double** [6][6]> qp_Cstar; // Stiffness matrix in material frame
Kokkos::View<double** [3]> qp_x0; // Initial position
Kokkos::View<double** [3]> qp_x0_prime; // Initial position derivative
Kokkos::View<double** [4]> qp_r0; // Initial rotation
Kokkos::View<double** [3]> qp_u; // State: translation displacement
Kokkos::View<double** [3]> qp_u_prime; // State: translation displacement derivative
Kokkos::View<double** [3]> qp_u_dot; // State: translation velocity
Kokkos::View<double** [3]> qp_u_ddot; // State: translation acceleration
Kokkos::View<double** [4]> qp_r; // State: rotation
Kokkos::View<double** [4]> qp_r_prime; // State: rotation derivative
Kokkos::View<double** [3]> qp_omega; // State: angular velocity
Kokkos::View<double** [3]> qp_omega_dot; // State: position/rotation
Kokkos::View<double** [3][4]> qp_E; // Quaternion derivative
Kokkos::View<double** [3][3]> qp_eta_tilde; //
Kokkos::View<double** [3][3]> qp_omega_tilde; //
Kokkos::View<double** [3][3]> qp_omega_dot_tilde; //
Kokkos::View<double** [3][3]> qp_x0pupss; //
Kokkos::View<double** [3][3]> qp_M_tilde; //
Kokkos::View<double** [3][3]> qp_N_tilde; //
Kokkos::View<double** [3]> qp_eta; //
Kokkos::View<double** [3][3]> qp_rho; //
Kokkos::View<double** [6]> qp_strain; // Strain
Kokkos::View<double** [6]> qp_Fc; // Elastic force
Kokkos::View<double** [6]> qp_Fd; // Elastic force
Kokkos::View<double** [6]> qp_Fi; // Inertial force
Kokkos::View<double** [6]> qp_Fg; // Gravity force
Kokkos::View<double** [6][6]> qp_RR0; // Global rotation
Kokkos::View<double** [6][6]> qp_Muu; // Mass in global frame
Kokkos::View<double** [6][6]> qp_Cuu; // Stiffness in global frame
Kokkos::View<double** [6][6]> qp_Ouu; // Linearization matrices
Kokkos::View<double** [6][6]> qp_Puu; // Linearization matrices
Kokkos::View<double** [6][6]> qp_Quu; // Linearization matrices
Kokkos::View<double** [6][6]> qp_Guu; // Linearization matrices
Kokkos::View<double** [6][6]> qp_Kuu; // Linearization matrices

View_NxN shape_interp; // shape function matrix for interpolation [Nodes x QPs]
View_NxN shape_deriv; // shape function matrix for derivative interp [Nodes x QPs]
Kokkos::View<double***> shape_interp; // shape function matrix for interpolation [Nodes x QPs]
Kokkos::View<double***>
shape_deriv; // shape function matrix for derivative interp [Nodes x QPs]

// Constructor which initializes views based on given sizes
Beams(
Expand All @@ -112,45 +113,45 @@ struct Beams {
node_FI("node_force_inertial", n_nodes),
node_FX("node_force_external", n_nodes),
// Quadrature Point data
qp_weight("qp_weight", n_qps),
qp_jacobian("qp_jacobian", n_qps),
qp_Mstar("qp_Mstar", n_qps),
qp_Cstar("qp_Cstar", n_qps),
qp_x0("qp_x0", n_qps),
qp_x0_prime("qp_x0_prime", n_qps),
qp_r0("qp_r0", n_qps),
qp_u("qp_u", n_qps),
qp_u_prime("qp_u_prime", n_qps),
qp_u_dot("qp_u_dot", n_qps),
qp_u_ddot("qp_u_ddot", n_qps),
qp_r("qp_r", n_qps),
qp_r_prime("qp_r_prime", n_qps),
qp_omega("qp_omega", n_qps),
qp_omega_dot("qp_omega_dot", n_qps),
qp_E("qp_E", n_qps),
qp_eta_tilde("R1_3x3", n_qps),
qp_omega_tilde("R1_3x3", n_qps),
qp_omega_dot_tilde("R1_3x3", n_qps),
qp_x0pupss("R1_3x3", n_qps),
qp_M_tilde("R1_3x3", n_qps),
qp_N_tilde("R1_3x3", n_qps),
qp_eta("V_3", n_qps),
qp_rho("R1_3x3", n_qps),
qp_strain("qp_strain", n_qps),
qp_Fc("qp_Fc", n_qps),
qp_Fd("qp_Fd", n_qps),
qp_Fi("qp_Fi", n_qps),
qp_Fg("qp_Fg", n_qps),
qp_RR0("qp_RR0", n_qps),
qp_Muu("qp_Muu", n_qps),
qp_Cuu("qp_Cuu", n_qps),
qp_Ouu("qp_Ouu", n_qps),
qp_Puu("qp_Puu", n_qps),
qp_Quu("qp_Quu", n_qps),
qp_Guu("qp_Guu", n_qps),
qp_Kuu("qp_Kuu", n_qps),
shape_interp("shape_interp", n_nodes, max_e_qps),
shape_deriv("deriv_interp", n_nodes, max_e_qps) {}
qp_weight("qp_weight", n_beams, max_elem_qps),
qp_jacobian("qp_jacobian", n_beams, max_elem_qps),
qp_Mstar("qp_Mstar", n_beams, max_elem_qps),
qp_Cstar("qp_Cstar", n_beams, max_elem_qps),
qp_x0("qp_x0", n_beams, max_elem_qps),
qp_x0_prime("qp_x0_prime", n_beams, max_elem_qps),
qp_r0("qp_r0", n_beams, max_elem_qps),
qp_u("qp_u", n_beams, max_elem_qps),
qp_u_prime("qp_u_prime", n_beams, max_elem_qps),
qp_u_dot("qp_u_dot", n_beams, max_elem_qps),
qp_u_ddot("qp_u_ddot", n_beams, max_elem_qps),
qp_r("qp_r", n_beams, max_elem_qps),
qp_r_prime("qp_r_prime", n_beams, max_elem_qps),
qp_omega("qp_omega", n_beams, max_elem_qps),
qp_omega_dot("qp_omega_dot", n_beams, max_elem_qps),
qp_E("qp_E", n_beams, max_elem_qps),
qp_eta_tilde("R1_3x3", n_beams, max_elem_qps),
qp_omega_tilde("R1_3x3", n_beams, max_elem_qps),
qp_omega_dot_tilde("R1_3x3", n_beams, max_elem_qps),
qp_x0pupss("R1_3x3", n_beams, max_elem_qps),
qp_M_tilde("R1_3x3", n_beams, max_elem_qps),
qp_N_tilde("R1_3x3", n_beams, max_elem_qps),
qp_eta("V_3", n_beams, max_elem_qps),
qp_rho("R1_3x3", n_beams, max_elem_qps),
qp_strain("qp_strain", n_beams, max_elem_qps),
qp_Fc("qp_Fc", n_beams, max_elem_qps),
qp_Fd("qp_Fd", n_beams, max_elem_qps),
qp_Fi("qp_Fi", n_beams, max_elem_qps),
qp_Fg("qp_Fg", n_beams, max_elem_qps),
qp_RR0("qp_RR0", n_beams, max_elem_qps),
qp_Muu("qp_Muu", n_beams, max_elem_qps),
qp_Cuu("qp_Cuu", n_beams, max_elem_qps),
qp_Ouu("qp_Ouu", n_beams, max_elem_qps),
qp_Puu("qp_Puu", n_beams, max_elem_qps),
qp_Quu("qp_Quu", n_beams, max_elem_qps),
qp_Guu("qp_Guu", n_beams, max_elem_qps),
qp_Kuu("qp_Kuu", n_beams, max_elem_qps),
shape_interp("shape_interp", n_beams, max_elem_nodes, max_elem_qps),
shape_deriv("deriv_interp", n_beams, max_elem_nodes, max_elem_qps) {}
};

} // namespace openturbine
21 changes: 13 additions & 8 deletions src/restruct_poc/beams/calculate_jacobian.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,20 +11,25 @@ namespace openturbine {

struct CalculateJacobian {
Kokkos::View<Beams::ElemIndices*>::const_type elem_indices; // Element indices
View_NxN::const_type shape_derivative; // Num Nodes x Num Quadrature points
View_Nx7::const_type node_position_rotation; // Node global position/rotation vector
View_Nx3 qp_position_derivative; // quadrature point position derivative
View_N qp_jacobian; // Jacobians
Kokkos::View<double***>::const_type shape_derivative; // Num Nodes x Num Quadrature points
View_Nx7::const_type node_position_rotation; // Node global position/rotation vector
Kokkos::View<double** [3]> qp_position_derivative; // quadrature point position derivative
View_NxN qp_jacobian; // Jacobians

KOKKOS_FUNCTION
void operator()(const int i_elem) const {
const auto& idx = elem_indices[i_elem];
const auto shape_deriv =
Kokkos::subview(shape_derivative, idx.node_range, idx.qp_shape_range);
const auto qp_pos_deriv = Kokkos::subview(qp_position_derivative, idx.qp_range, Kokkos::ALL);
const auto shape_deriv = Kokkos::subview(
shape_derivative, i_elem, Kokkos::make_pair(size_t{0U}, idx.num_nodes),
idx.qp_shape_range
);
const auto qp_pos_deriv = Kokkos::subview(
qp_position_derivative, i_elem, Kokkos::pair(size_t{0U}, idx.num_qps), Kokkos::ALL
);
const auto node_pos =
Kokkos::subview(node_position_rotation, idx.node_range, Kokkos::make_pair(0, 3));
const auto qp_jacob = Kokkos::subview(qp_jacobian, idx.qp_range);
const auto qp_jacob =
Kokkos::subview(qp_jacobian, i_elem, Kokkos::make_pair(size_t{0U}, idx.num_qps));

InterpVector3(shape_deriv, node_pos, qp_pos_deriv);

Expand Down
18 changes: 13 additions & 5 deletions src/restruct_poc/beams/create_beams.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,11 +62,19 @@ inline Beams CreateBeams(const BeamsInput& beams_input) {
PopulateElementViews(
beams_input.elements[i], // Element inputs
Kokkos::subview(host_node_x0, idx.node_range, Kokkos::ALL),
Kokkos::subview(host_qp_weight, idx.qp_range),
Kokkos::subview(host_qp_Mstar, idx.qp_range, Kokkos::ALL, Kokkos::ALL),
Kokkos::subview(host_qp_Cstar, idx.qp_range, Kokkos::ALL, Kokkos::ALL),
Kokkos::subview(host_shape_interp, idx.node_range, idx.qp_shape_range),
Kokkos::subview(host_shape_deriv, idx.node_range, idx.qp_shape_range)
Kokkos::subview(host_qp_weight, i, Kokkos::pair(size_t{0U}, idx.num_qps)),
Kokkos::subview(
host_qp_Mstar, i, Kokkos::pair(size_t{0U}, idx.num_qps), Kokkos::ALL, Kokkos::ALL
),
Kokkos::subview(
host_qp_Cstar, i, Kokkos::pair(size_t{0U}, idx.num_qps), Kokkos::ALL, Kokkos::ALL
),
Kokkos::subview(
host_shape_interp, i, Kokkos::pair(size_t{0U}, idx.num_nodes), idx.qp_shape_range
),
Kokkos::subview(
host_shape_deriv, i, Kokkos::pair(size_t{0U}, idx.num_nodes), idx.qp_shape_range
)
);
}

Expand Down
14 changes: 9 additions & 5 deletions src/restruct_poc/beams/interpolate_QP_position.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,18 +10,22 @@ namespace openturbine {

struct InterpolateQPPosition {
Kokkos::View<Beams::ElemIndices*>::const_type elem_indices; // Element indices
View_NxN::const_type shape_interpolation; // Num Nodes x Num Quadrature points
Kokkos::View<double***>::const_type shape_interpolation; // Num Nodes x Num Quadrature points
View_Nx7::const_type node_position_rotation; // Node global position vector
View_Nx3 qp_position; // quadrature point position
Kokkos::View<double** [3]> qp_position; // quadrature point position

KOKKOS_FUNCTION
void operator()(const int i_elem) const {
const auto& idx = elem_indices[i_elem];
const auto shape_interp =
Kokkos::subview(shape_interpolation, idx.node_range, idx.qp_shape_range);
const auto shape_interp = Kokkos::subview(
shape_interpolation, i_elem, Kokkos::make_pair(size_t{0U}, idx.num_nodes),
idx.qp_shape_range
);
const auto node_pos =
Kokkos::subview(node_position_rotation, idx.node_range, Kokkos::make_pair(0, 3));
const auto qp_pos = Kokkos::subview(qp_position, idx.qp_range, Kokkos::ALL);
const auto qp_pos = Kokkos::subview(
qp_position, i_elem, Kokkos::make_pair(size_t{0U}, idx.num_qps), Kokkos::ALL
);

for (auto j = 0U; j < idx.num_qps; ++j) {
auto local_result = Kokkos::Array<double, 3>{};
Expand Down
14 changes: 9 additions & 5 deletions src/restruct_poc/beams/interpolate_QP_rotation.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,18 +11,22 @@ namespace openturbine {

struct InterpolateQPRotation {
Kokkos::View<Beams::ElemIndices*>::const_type elem_indices; // Element indices
View_NxN::const_type shape_interpolation; // Num Nodes x Num Quadrature points
Kokkos::View<double***>::const_type shape_interpolation; // Num Nodes x Num Quadrature points
View_Nx7::const_type node_position_rotation; // Node global position vector
View_Nx4 qp_rotation; // quadrature point rotation
Kokkos::View<double** [4]> qp_rotation; // quadrature point rotation

KOKKOS_FUNCTION
void operator()(const int i_elem) const {
const auto& idx = elem_indices[i_elem];
const auto shape_interp =
Kokkos::subview(shape_interpolation, idx.node_range, idx.qp_shape_range);
const auto shape_interp = Kokkos::subview(
shape_interpolation, i_elem, Kokkos::make_pair(size_t{0U}, idx.num_nodes),
idx.qp_shape_range
);
const auto node_rot =
Kokkos::subview(node_position_rotation, idx.node_range, Kokkos::make_pair(3, 7));
const auto qp_rot = Kokkos::subview(qp_rotation, idx.qp_range, Kokkos::ALL);
const auto qp_rot = Kokkos::subview(
qp_rotation, i_elem, Kokkos::make_pair(size_t{0U}, idx.num_qps), Kokkos::ALL
);

InterpQuaternion(shape_interp, node_rot, qp_rot);
}
Expand Down
11 changes: 7 additions & 4 deletions src/restruct_poc/beams/interpolate_QP_rotation_derivative.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,15 +10,18 @@ namespace openturbine {

struct InterpolateQPRotationDerivative {
Kokkos::View<Beams::ElemIndices*>::const_type elem_indices; // Element indices
View_NxN::const_type shape_derivative; // Num Nodes x Num Quadrature points
Kokkos::View<double***>::const_type shape_derivative; // Num Nodes x Num Quadrature points
View_N::const_type qp_jacobian; // Jacobians
View_Nx7::const_type node_pos_rot; // Node global position/rotation vector
View_Nx4 qp_rot_deriv; // quadrature point rotation derivative
View_Nx7::const_type node_pos_rot; // Node global position/rotation vector
Kokkos::View<double** [4]> qp_rot_deriv; // quadrature point rotation derivative

KOKKOS_FUNCTION
void operator()(const int i_elem) const {
auto& idx = elem_indices[i_elem];
auto shape_deriv = Kokkos::subview(shape_derivative, idx.node_range, idx.qp_shape_range);
auto shape_deriv = Kokkos::subview(
shape_derivative, i_elem, Kokkos::make_pair(size_t{0U}, idx.num_nodes),
idx.qp_shape_range
);
auto qp_rot_deriv = Kokkos::subview(qp_rot_deriv, idx.qp_range, Kokkos::ALL);
auto node_rot = Kokkos::subview(node_pos_rot, idx.node_range, Kokkos::make_pair(3, 7));
auto qp_jac = Kokkos::subview(qp_jacobian, idx.qp_range);
Expand Down
Loading

0 comments on commit afea395

Please sign in to comment.