resultevaluators.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 <dune/common/math.hh>
 
#include <Eigen/Geometry>
 
#include <ikarus/finiteelements/mechanics/materials.hh>
#include <ikarus/utils/tensorutils.hh>
 
namespace Ikarus::ResultEvaluators {
 
struct VonMises
{
  template <typename R>
  double operator()(const R& resultArray, [[maybe_unused]] const auto& pos, [[maybe_unused]] const auto& fe,
                    [[maybe_unused]] const int comp) const {
    auto sigma = fromVoigt(resultArray, false);
 
    auto I2 = 1.0 / 2.0 * (sigma.squaredNorm() - 1.0 / 3.0 * Dune::power(sigma.trace(), 2));
    return std::sqrt(3.0 * I2);
  }
 
  constexpr static std::string name() { return "VonMises"; }
 
  constexpr static int ncomps() { return 1; }
};
 
struct HydrostaticStress
{
  template <typename R>
  double operator()(const R& resultArray, [[maybe_unused]] const auto& pos, [[maybe_unused]] const auto& fe,
                    [[maybe_unused]] const int comp) const {
    const auto sigma = fromVoigt(resultArray, false);
    return 1.0 / sigma.rows() * sigma.trace();
  }
 
  constexpr static std::string name() { return "HydrostaticStress"; }
 
  constexpr static int ncomps() { return 1; }
};
 
template <int dim>
requires(dim == 2 or dim == 3)
struct PrincipalStress
{
  double operator()(const auto& resultArray, [[maybe_unused]] const auto& pos, [[maybe_unused]] const auto& fe,
                    const int comp) const {
    auto mat = fromVoigt(resultArray, false);
    Eigen::SelfAdjointEigenSolver<decltype(mat)> eigensolver(mat, Eigen::EigenvaluesOnly);
    return eigensolver.eigenvalues().reverse()[comp];
  }
 
  constexpr static std::string name() { return "PrincipalStress"; }
 
  constexpr static int ncomps() { return dim; }
};
 
struct Triaxiality
{
  template <typename R>
  double operator()(const R& resultArray, [[maybe_unused]] const auto& pos, [[maybe_unused]] const auto& fe,
                    [[maybe_unused]] const int comp) const {
    auto sigeq = VonMises{}(resultArray, pos, fe, 0);
    auto sigm  = HydrostaticStress{}(resultArray, pos, fe, 0);
    return sigm / sigeq;
  }
  constexpr static std::string name() { return "Triaxiality"; }
 
  constexpr static int ncomps() { return 1; }
};
 
struct PolarStress
{
  PolarStress(const Dune::FieldVector<double, 2>& origin)
      : origin_(origin) {}
 
  template <typename R>
  double operator()(const R& resultArray, const auto& pos, const auto& fe, const int comp) const {
    static_assert(R::CompileTimeTraits::RowsAtCompileTime == 3, "PolarStress is only valid for 2D.");
    if (comp > 2)
      DUNE_THROW(Dune::RangeError, "PolarStress: Comp out of range.");
 
    // Offset to center the coordinate system in the reference geometry
    Dune::FieldVector<double, 2> posGlobal = fe.geometry().global(pos) - origin_;
    auto theta                             = std::atan2(posGlobal[1], posGlobal[0]);
 
    const auto sigma = fromVoigt(resultArray, false);
    Eigen::Rotation2D<double> r(theta);
 
    // deliberately not evaluating this that it stays an expression for below
    auto polarStress = r.inverse() * sigma * r;
    switch (comp) {
      case 0:
        return polarStress(0, 0);
      case 1:
        return polarStress(1, 1);
      case 2:
        return polarStress(1, 0);
      default:
        __builtin_unreachable();
    }
  }
 
  constexpr static std::string name() { return "PolarStress"; }
 
  constexpr static int ncomps() { return 3; }
 
private:
  Dune::FieldVector<double, 2> origin_;
};
 
} // namespace Ikarus::ResultEvaluators