materialhelpers.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers mueller@ibb.uni-stuttgart.de
// SPDX-License-Identifier: LGPL-3.0-or-later
 
#pragma once
 
#include <ranges>
 
#include <dune/common/float_cmp.hh>
 
#include <Eigen/Dense>
 
#include <ikarus/utils/concepts.hh>
#include <ikarus/utils/tensorutils.hh>
 
namespace Ikarus::Materials {
struct MatrixIndexPair
{
  Eigen::Index row; 
  Eigen::Index col; 
};
} // namespace Ikarus::Materials
namespace Ikarus::Materials::Impl {
 
template <size_t size>
consteval auto createfreeVoigtIndices(const std::array<MatrixIndexPair, size>& fixed) {
  std::array<size_t, 6 - size> res{};
  std::array<size_t, size> voigtFixedIndices;
  std::ranges::transform(fixed, voigtFixedIndices.begin(), [](auto pair) { return toVoigt(pair.row, pair.col); });
  std::ranges::sort(voigtFixedIndices);
  std::ranges::set_difference(std::ranges::iota_view(size_t(0), size_t(6)), voigtFixedIndices, res.begin());
  std::ranges::sort(res);
  return res;
}
 
template <size_t size>
consteval auto createFixedVoigtIndices(const std::array<MatrixIndexPair, size>& fixed) {
  std::array<size_t, size> fixedIndices;
  std::ranges::transform(fixed, fixedIndices.begin(), [](auto pair) { return toVoigt(pair.row, pair.col); });
  std::ranges::sort(fixedIndices);
  return fixedIndices;
}
 
template <size_t size>
constexpr size_t countDiagonalIndices(const std::array<MatrixIndexPair, size>& fixed) {
  size_t count = 0;
  for (auto v : fixed) {
    if (v.col == v.row)
      ++count;
  }
  return count;
}
 
template <typename Derived>
decltype(auto) maybeFromVoigt(const Eigen::MatrixBase<Derived>& E) {
  if constexpr (Concepts::EigenVector<Derived>) { // receiving vector means Voigt notation
    return fromVoigt(E.derived(), true);
  } else
    return E.derived();
}
 
template <typename ScalarType>
void checkPositiveOrAbort(ScalarType det) {
  if (Dune::FloatCmp::le(static_cast<double>(det), 0.0, 1e-10)) {
    std::cerr << "Determinant of right Cauchy Green tensor C must be greater than zero. detC = " +
                     std::to_string(static_cast<double>(det));
    abort();
  }
}
 
template <typename ScalarType, typename Derived>
auto principalStretches(const Eigen::MatrixBase<Derived>& C, int options = Eigen::ComputeEigenvectors) {
  Eigen::SelfAdjointEigenSolver<Derived> eigensolver{};
 
  eigensolver.compute(C, options);
 
  if (eigensolver.info() != Eigen::Success)
    DUNE_THROW(Dune::MathError, "Failed to compute eigenvalues and eigenvectors of C.");
 
  auto& eigenvalues  = eigensolver.eigenvalues();
  auto& eigenvectors = options == Eigen::ComputeEigenvectors ? eigensolver.eigenvectors() : Derived::Zero();
 
  auto principalStretches = eigenvalues.cwiseSqrt().eval();
  return std::make_pair(principalStretches, eigenvectors);
}
 
template <Concepts::EigenVector3 Vector>
inline Vector::Scalar determinantFromPrincipalValues(const Vector& principalValues) {
  return principalValues.prod();
}
 
template <Concepts::EigenVector3 Vector>
inline Vector deviatoricStretches(const Vector& lambda) {
  using ScalarType = typename Vector::Scalar;
  ScalarType J     = determinantFromPrincipalValues(lambda);
  ScalarType Jmod  = pow(J, -1.0 / 3.0);
  return Jmod * lambda;
}
 
template <Concepts::EigenVector3 Vector>
inline Vector invariants(const Vector& lambda) {
  using ScalarType   = typename Vector::Scalar;
  auto lambdaSquared = lambda.array().square();
  auto invariants    = Vector::Zero().eval();
 
  invariants[0] = lambdaSquared.sum();
  invariants[1] =
      lambdaSquared[0] * lambdaSquared[1] + lambdaSquared[1] * lambdaSquared[2] + lambdaSquared[0] * lambdaSquared[2];
  invariants[2] = determinantFromPrincipalValues(lambdaSquared);
 
  return invariants;
}
} // namespace Ikarus::Materials::Impl