KelvinVector.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#pragma once
 
#include <Eigen/Core>
 
#include "BaseLib/Error.h"
 
namespace MathLib
{
namespace KelvinVector
{
constexpr int kelvin_vector_dimensions(int const displacement_dim)
{
    if (displacement_dim == 2)
    {
        return 4;
    }
    else if (displacement_dim == 3)
    {
        return 6;
    }
    OGS_FATAL(
        "Cannot convert displacement dimension {} to kelvin vector dimension.",
        displacement_dim);
}
 
//
// Kelvin vector and matrix templates for given displacement dimension.
//

template <int DisplacementDim>
using KelvinVectorType =
    Eigen::Matrix<double, kelvin_vector_dimensions(DisplacementDim), 1,
                  Eigen::ColMajor>;

template <int DisplacementDim>
using KelvinMatrixType =
    Eigen::Matrix<double, kelvin_vector_dimensions(DisplacementDim),
                  kelvin_vector_dimensions(DisplacementDim), Eigen::RowMajor>;

template <int DisplacementDim>
constexpr auto KVzero()
{
    return KelvinVectorType<DisplacementDim>::Zero();
}

template <int DisplacementDim>
constexpr auto KMzero()
{
    return KelvinMatrixType<DisplacementDim>::Zero();
}

template <int DisplacementDim>
constexpr auto KVnan()
{
    return KelvinVectorType<DisplacementDim>::Constant(
        std::numeric_limits<double>::quiet_NaN());
}

template <int DisplacementDim>
constexpr auto KMnan()
{
    return KelvinMatrixType<DisplacementDim>::Constant(
        std::numeric_limits<double>::quiet_NaN());
}

template <int KelvinVectorSize>
struct Invariants final
{
    static_assert(KelvinVectorSize == 4 || KelvinVectorSize == 6,
                  "KelvinVector invariants for vectors of size different than "
                  "4 or 6 is not allowed.");
    static Eigen::Matrix<double, KelvinVectorSize, KelvinVectorSize> const
        deviatoric_projection;
    static Eigen::Matrix<double, KelvinVectorSize, KelvinVectorSize> const
        spherical_projection;
    static Eigen::Matrix<double, KelvinVectorSize, 1> const identity2;

    static Eigen::Matrix<double, KelvinVectorSize, 1> const ones2;

    static double determinant(
        Eigen::Matrix<double, KelvinVectorSize, 1> const& v);

    static double equivalentStress(
        Eigen::Matrix<double, KelvinVectorSize, 1> const& deviatoric_v);

    static double FrobeniusNorm(
        Eigen::Matrix<double, KelvinVectorSize, 1> const& deviatoric_v);

    static double J2(
        Eigen::Matrix<double, KelvinVectorSize, 1> const& deviatoric_v);

    static double J3(
        Eigen::Matrix<double, KelvinVectorSize, 1> const& deviatoric_v);

    static double trace(Eigen::Matrix<double, KelvinVectorSize, 1> const& v);

    static Eigen::Vector3d diagonal(
        Eigen::Matrix<double, KelvinVectorSize, 1> const& v);
};
 
//
// Inverses of a Kelvin vector.
//

template <int KelvinVectorSize>
Eigen::Matrix<double, KelvinVectorSize, 1, Eigen::ColMajor, KelvinVectorSize, 1>
inverse(Eigen::Matrix<double,
                      KelvinVectorSize,
                      1,
                      Eigen::ColMajor,
                      KelvinVectorSize,
                      1> const& v);

template <int KelvinVectorSize>
Eigen::Matrix<double, 3, 3> kelvinVectorToTensor(Eigen::Matrix<double,
                                                               KelvinVectorSize,
                                                               1,
                                                               Eigen::ColMajor,
                                                               KelvinVectorSize,
                                                               1> const& v);

template <int DisplacementDim>
KelvinVectorType<DisplacementDim> tensorToKelvin(
    Eigen::Matrix<double, 3, 3> const& m);

template <int KelvinVectorSize>
Eigen::Matrix<double, KelvinVectorSize, 1, Eigen::ColMajor, KelvinVectorSize, 1>
kelvinVectorToSymmetricTensor(Eigen::Matrix<double,
                                            KelvinVectorSize,
                                            1,
                                            Eigen::ColMajor,
                                            KelvinVectorSize,
                                            1> const& v);

template <typename Derived>
Eigen::Matrix<double, Eigen::MatrixBase<Derived>::RowsAtCompileTime, 1>
symmetricTensorToKelvinVector(Eigen::MatrixBase<Derived> const& v)
{
    static_assert(
        (Eigen::MatrixBase<Derived>::ColsAtCompileTime == 1) ||
            (Eigen::MatrixBase<Derived>::ColsAtCompileTime == Eigen::Dynamic),
        "KelvinVector must be a column vector");
    if (v.cols() != 1)
    {
        OGS_FATAL(
            "KelvinVector must be a column vector, but input has {:d} columns.",
            v.cols());
    }
 
    Eigen::Matrix<double, Eigen::MatrixBase<Derived>::RowsAtCompileTime, 1>
        result;
    if (v.rows() == 4)
    {
        result.resize(4, 1);
        result << v[0], v[1], v[2], v[3] * std::sqrt(2.);
    }
    else if (v.rows() == 6)
    {
        result.resize(6, 1);
        result << v[0], v[1], v[2], v[3] * std::sqrt(2.), v[4] * std::sqrt(2.),
            v[5] * std::sqrt(2.);
    }
    else
    {
        OGS_FATAL(
            "Symmetric tensor to Kelvin vector conversion expected an input "
            "vector of size 4 or 6, but a vector of size {:d} was given.",
            v.size());
    }
    return result;
}

template <int DisplacementDim>
KelvinVectorType<DisplacementDim> symmetricTensorToKelvinVector(
    std::vector<double> const& values)
{
    constexpr int kelvin_vector_size =
        kelvin_vector_dimensions(DisplacementDim);
 
    if (values.size() != kelvin_vector_size)
    {
        OGS_FATAL(
            "Symmetric tensor to Kelvin vector conversion expected an input "
            "vector of size {:d}, but a vector of size {:d} was given.",
            kelvin_vector_size, values.size());
    }
 
    return symmetricTensorToKelvinVector(
        Eigen::Map<typename MathLib::KelvinVector::KelvinVectorType<
            DisplacementDim> const>(
            values.data(), kelvin_vector_dimensions(DisplacementDim), 1));
}

template <int DisplacementDim>
Eigen::Matrix<double, DisplacementDim,
              kelvin_vector_dimensions(DisplacementDim)>
liftVectorToKelvin(Eigen::Matrix<double, DisplacementDim, 1> const& v);

template <int DisplacementDim>
Eigen::Matrix<double, DisplacementDim, 1> reduceKelvinToVector(
    Eigen::Matrix<double, DisplacementDim,
                  kelvin_vector_dimensions(DisplacementDim)> const& m);

template <int DisplacementDim>
KelvinMatrixType<DisplacementDim> fourthOrderRotationMatrix(
    Eigen::Matrix<double, DisplacementDim, DisplacementDim, Eigen::ColMajor,
                  DisplacementDim, DisplacementDim> const& transformation);
 
}  // namespace KelvinVector
}  // namespace MathLib
 
#include "KelvinVector-impl.h"