DeformationGradient.h

Go to the documentation of this file.

 
#ifndef PBAT_FEM_DEFORMATIONGRADIENT_H
#define PBAT_FEM_DEFORMATIONGRADIENT_H
 
#include "Concepts.h"
#include "pbat/Aliases.h"
#include "pbat/math/linalg/mini/BinaryOperations.h"
#include "pbat/math/linalg/mini/Concepts.h"
#include "pbat/math/linalg/mini/Matrix.h"
 
namespace pbat::fem {

template <CElement TElement, class TDerivedx, class TDerivedX>
auto DeformationGradient(
    Eigen::MatrixBase<TDerivedx> const& x,
    Eigen::MatrixBase<TDerivedX> const& GP) -> Matrix<TDerivedx::RowsAtCompileTime, TElement::kDims>
{
    return x * GP;
}

template <
    CElement TElement,
    int Dims,
    math::linalg::mini::CMatrix TMatrixGF,
    math::linalg::mini::CMatrix TMatrixGP,
    class ScalarType = typename TMatrixGF::ScalarType>
auto GradientSegmentWrtDofs(TMatrixGF const& GF, TMatrixGP const& GP, auto i)
    -> math::linalg::mini::SVector<ScalarType, Dims>
{
    using namespace math::linalg::mini;
    SVector<ScalarType, Dims> dPsidx = Zeros<ScalarType, Dims, 1>{};
    for (auto k = 0; k < Dims; ++k)
    {
        dPsidx += GP(i, k) * GF.template Slice<Dims, 1>(k * Dims, 0);
    }
    return dPsidx;
}

template <
    CElement TElement,
    int Dims,
    math::linalg::mini::CMatrix TMatrixGF,
    math::linalg::mini::CMatrix TMatrixGP,
    class ScalarType = typename TMatrixGF::ScalarType>
auto GradientWrtDofs(TMatrixGF const& GF, TMatrixGP const& GP)
    -> math::linalg::mini::SVector<ScalarType, TElement::kNodes * Dims>
{
    auto constexpr kRows = TElement::kNodes * Dims;
    using namespace math::linalg::mini;
    SVector<ScalarType, kRows> dPsidx = Zeros<ScalarType, kRows, 1>{};
    for (auto k = 0; k < Dims; ++k)
    {
        for (auto i = 0; i < TElement::kNodes; ++i)
        {
            dPsidx.template Slice<Dims, 1>(i * Dims, 0) +=
                GP(i, k) * GF.template Slice<Dims, 1>(k * Dims, 0);
        }
    }
    return dPsidx;
}

template <
    CElement TElement,
    int Dims,
    math::linalg::mini::CMatrix TMatrixHF,
    math::linalg::mini::CMatrix TMatrixGP,
    class ScalarType = typename TMatrixHF::ScalarType>
auto HessianBlockWrtDofs(TMatrixHF const& HF, TMatrixGP const& GP, auto i, auto j)
    -> math::linalg::mini::SMatrix<ScalarType, Dims, Dims>
{
    using namespace math::linalg::mini;
    SMatrix<ScalarType, Dims, Dims> d2Psidx2 = Zeros<ScalarType, Dims, Dims>{};
    for (auto kj = 0; kj < Dims; ++kj)
    {
        for (auto ki = 0; ki < Dims; ++ki)
        {
            d2Psidx2 += GP(i, ki) * GP(j, kj) * HF.template Slice<Dims, Dims>(ki * Dims, kj * Dims);
        }
    }
    return d2Psidx2;
}

template <
    CElement TElement,
    int Dims,
    math::linalg::mini::CMatrix TMatrixHF,
    math::linalg::mini::CMatrix TMatrixGP,
    class ScalarType = typename TMatrixHF::ScalarType>
auto HessianWrtDofs(TMatrixHF const& HF, TMatrixGP const& GP)
    -> math::linalg::mini::SMatrix<ScalarType, TElement::kNodes * Dims, TElement::kNodes * Dims>
{
    auto constexpr kRows = TElement::kNodes * Dims;
    auto constexpr kCols = TElement::kNodes * Dims;
    using namespace math::linalg::mini;
    SMatrix<ScalarType, kRows, kCols> d2Psidx2 = Zeros<ScalarType, kRows, kCols>{};
    for (auto kj = 0; kj < Dims; ++kj)
    {
        for (auto ki = 0; ki < Dims; ++ki)
        {
            for (auto j = 0; j < TElement::kNodes; ++j)
            {
                for (auto i = 0; i < TElement::kNodes; ++i)
                {
                    d2Psidx2.template Slice<Dims, Dims>(i * Dims, j * Dims) +=
                        GP(i, ki) * GP(j, kj) * HF.template Slice<Dims, Dims>(ki * Dims, kj * Dims);
                }
            }
        }
    }
    return d2Psidx2;
}
 
} // namespace pbat::fem
 
#endif // PBAT_FEM_DEFORMATIONGRADIENT_H