matrix_operations.hpp

Go to the documentation of this file.

#pragma once
 
#include <array>
#include <ranges>
 
#include <Kokkos_Core.hpp>
 
#include "math/quaternion_operations.hpp"
 
namespace kynema::math {
 
template <typename Matrix>
KOKKOS_INLINE_FUNCTION void AX_Matrix(const Matrix& A, const Matrix& AX_A) {
    double trace{0.};
    for (auto i = 0; i < A.extent_int(0); ++i) {
        trace += A(i, i);
    }
 
    trace /= 2.;
 
    for (auto i = 0; i < A.extent_int(0); ++i) {
        for (auto j = 0; j < A.extent_int(1); ++j) {
            AX_A(i, j) = -A(i, j) / 2.;
        }
        AX_A(i, i) += trace;
    }
}
 
template <typename Matrix, typename Vector>
KOKKOS_INLINE_FUNCTION void AxialVectorOfMatrix(const Matrix& m, const Vector& v) {
    v(0) = (m(2, 1) - m(1, 2)) / 2.;
    v(1) = (m(0, 2) - m(2, 0)) / 2.;
    v(2) = (m(1, 0) - m(0, 1)) / 2.;
}
 
inline std::array<std::array<double, 6>, 6> RotateMatrix6(
    const std::array<std::array<double, 6>, 6>& m, const std::array<double, 4>& q
) {
    const auto rm = QuaternionToRotationMatrix(q);
    std::array<std::array<double, 6>, 6> r{{
        {rm[0][0], rm[0][1], rm[0][2], 0., 0., 0.},
        {rm[1][0], rm[1][1], rm[1][2], 0., 0., 0.},
        {rm[2][0], rm[2][1], rm[2][2], 0., 0., 0.},
        {0., 0., 0., rm[0][0], rm[0][1], rm[0][2]},
        {0., 0., 0., rm[1][0], rm[1][1], rm[1][2]},
        {0., 0., 0., rm[2][0], rm[2][1], rm[2][2]},
    }};
 
    // matmul(r,m)
    std::array<std::array<double, 6>, 6> mt{};
    for (auto i : std::views::iota(0U, 6U)) {
        for (auto j : std::views::iota(0U, 6U)) {
            mt[i][j] = 0.;
            for (auto k : std::views::iota(0U, 6U)) {
                mt[i][j] += r[i][k] * m[k][j];
            }
        }
    }
 
    // matmul(matmul(r,m),r^T)
    std::array<std::array<double, 6>, 6> mo{};
    for (auto i : std::views::iota(0U, 6U)) {
        for (auto j : std::views::iota(0U, 6U)) {
            mo[i][j] = 0.;
            for (auto k : std::views::iota(0U, 6U)) {
                mo[i][j] += mt[i][k] * r[j][k];
            }
        }
    }
 
    return mo;
}
 
}  // namespace kynema::math