enhancedassumedstrains.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers mueller@ibb.uni-stuttgart.de
// SPDX-License-Identifier: LGPL-3.0-or-later
 
#pragma once
 
#if HAVE_DUNE_LOCALFEFUNCTIONS
  #include <dune/localfefunctions/derivativetransformators.hh>
  #include <dune/localfefunctions/linearAlgebraHelper.hh>
  #include <dune/localfefunctions/meta.hh>
 
  #include <ikarus/finiteelements/fehelper.hh>
  #include <ikarus/finiteelements/ferequirements.hh>
  #include <ikarus/finiteelements/mechanics/easvariants.hh>
  #include <ikarus/finiteelements/mixin.hh>
 
namespace Ikarus {
 
template <typename PreFE, typename FE>
class EnhancedAssumedStrains;
 
struct EnhancedAssumedStrainsPre
{
  int numberOfParameters{};
 
  template <typename PreFE, typename FE>
  using Skill = EnhancedAssumedStrains<PreFE, FE>;
};
 
template <typename PreFE, typename FE>
class EnhancedAssumedStrains
{
public:
  using Traits            = PreFE::Traits;
  using FERequirementType = typename Traits::FERequirementType;
  using LocalView         = typename Traits::LocalView;
  using Geometry          = typename Traits::Geometry;
  using GridView          = typename Traits::GridView;
  using Pre               = EnhancedAssumedStrainsPre;
 
  using SupportedResultTypes = std::tuple<decltype(makeRT<ResultTypes::linearStress>())>;
 
  explicit EnhancedAssumedStrains(const Pre& pre) { this->setEASType(pre.numberOfParameters); }
 
  bool isDisplacementBased() const { return easVariant_.isDisplacmentBased(); }
 
  const auto& easVariant() const { return easVariant_; }
 
  auto numberOfEASParameters() const { return easVariant_.numberOfEASParameters(); }
 
  template <template <typename, int, int> class RT>
  auto calculateAtImpl(const FERequirementType& req, const Dune::FieldVector<double, Traits::mydim>& local,
                       Dune::PriorityTag<2>) const {
    if (isDisplacementBased())
      return underlying().template calculateAtImpl<RT>(req, local, Dune::PriorityTag<1>());
 
    using namespace Dune::Indices;
    using namespace Dune::DerivativeDirections;
    auto resultVector = underlying().template calculateAtImpl<RT>(req, local, Dune::PriorityTag<1>());
    if constexpr (isSameResultType<RT, ResultTypes::linearStress>) {
      using RTWrapper = ResultWrapper<RT<typename Traits::ctype, Traits::mydim, Traits::worlddim>, ResultShape::Vector>;
      RTWrapper resultWrapper;
 
      auto calculateAtContribution = [&]<typename EAST>(const EAST& easFunction) {
        typename EAST::DType D;
        calculateDAndLMatrix(easFunction, req, D, L_);
        const auto ufunc = underlying().displacementFunction(req);
        const auto disp  = Dune::viewAsFlatEigenVector(ufunc.coefficientsRef());
        const auto C     = underlying().materialTangentFunction(req);
        const auto alpha = (-D.inverse() * L_ * disp).eval();
        const auto M     = easFunction.calcM(local);
        const auto CEval = C(local);
        auto easStress   = (CEval * M * alpha).eval();
        resultWrapper    = resultVector.asVec() + easStress;
      };
      easVariant_(calculateAtContribution);
      return resultWrapper;
    }
  }
 
  void setEASType(int numberOfEASParameters) {
    if (numberOfEASParameters != 0)
      easApplicabilityCheck();
    easVariant_.setEASType(numberOfEASParameters);
  }
 
protected:
  void bindImpl() {
    assert(underlying().localView().bound());
    easVariant_.bind(underlying().localView().element().geometry());
  }
 
public:
protected:
  inline void easApplicabilityCheck() const {
    const auto& numberOfNodes = underlying().numberOfNodes();
    assert(not(not((numberOfNodes == 4 and Traits::mydim == 2) or (numberOfNodes == 8 and Traits::mydim == 3)) and
               (not isDisplacementBased())) &&
           "EAS only supported for Q1 or H1 elements");
  }
 
  template <typename ScalarType>
  void calculateMatrixImpl(
      const FERequirementType& par, typename Traits::template MatrixType<> K,
      const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
    using namespace Dune::DerivativeDirections;
    using namespace Dune;
    easApplicabilityCheck();
    if (isDisplacementBased())
      return;
 
    decltype(auto) LMat = [this]() -> decltype(auto) {
      if constexpr (std::is_same_v<ScalarType, double>)
        return [this]() -> Eigen::MatrixXd& { return L_; }();
      else
        return Eigen::MatrixX<ScalarType>{};
    }();
 
    auto calculateMatrixContribution = [&]<typename EAST>(const EAST& easFunction) {
      typename EAST::DType D;
      calculateDAndLMatrix(easFunction, par, D, LMat, dx);
 
      K.template triangularView<Eigen::Upper>() -= LMat.transpose() * D.inverse() * LMat;
      K.template triangularView<Eigen::StrictlyLower>() = K.transpose();
    };
    easVariant_(calculateMatrixContribution);
  }
 
  template <typename ScalarType>
  inline ScalarType calculateScalarImpl(
      const FERequirementType& par,
      const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
    easApplicabilityCheck();
    if (isDisplacementBased())
      return 0.0;
    DUNE_THROW(Dune::NotImplemented,
               "EAS element do not support any scalar calculations, i.e. they are not derivable from a potential");
  }
 
  template <typename ScalarType>
  void calculateVectorImpl(
      const FERequirementType& par, typename Traits::template VectorType<ScalarType> force,
      const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
    easApplicabilityCheck();
    using namespace Dune;
    using namespace Dune::DerivativeDirections;
    const auto uFunction      = underlying().displacementFunction(par, dx);
    auto strainFunction       = underlying().strainFunction(par, dx);
    const auto& numberOfNodes = underlying().numberOfNodes();
 
    auto calculateForceContribution = [&]<typename EAST>(const EAST& easFunction) {
      typename EAST::DType D;
      calculateDAndLMatrix(easFunction, par, D, L_);
 
      decltype(auto) LMat = [this]() -> decltype(auto) {
        if constexpr (std::is_same_v<ScalarType, double>)
          return [this]() -> Eigen::MatrixXd& { return L_; }();
        else
          return Eigen::MatrixX<ScalarType>{};
      }();
      const auto disp  = Dune::viewAsFlatEigenVector(uFunction.coefficientsRef());
      const auto alpha = (-D.inverse() * L_ * disp).eval();
      const auto geo   = underlying().localView().element().geometry();
      auto C           = underlying().materialTangentFunction(par);
 
      for (const auto& [gpIndex, gp] : strainFunction.viewOverIntegrationPoints()) {
        const auto M            = easFunction.calcM(gp.position());
        const double intElement = geo.integrationElement(gp.position()) * gp.weight();
        const auto CEval        = C(gpIndex);
        auto stresses           = (CEval * M * alpha).eval();
        for (size_t i = 0; i < numberOfNodes; ++i) {
          const auto bopI = strainFunction.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
          force.template segment<Traits::worlddim>(Traits::worlddim * i) += bopI.transpose() * stresses * intElement;
        }
      }
    };
    easVariant_(calculateForceContribution);
  }
 
private:
  EAS::Impl::EASVariant<Geometry> easVariant_;
  mutable Eigen::MatrixXd L_;
 
  //> CRTP
  const auto& underlying() const { return static_cast<const FE&>(*this); }
  auto& underlying() { return static_cast<FE&>(*this); }
 
  template <typename ScalarType, int enhancedStrainSize>
  void calculateDAndLMatrix(
      const auto& easFunction, const auto& par, Eigen::Matrix<double, enhancedStrainSize, enhancedStrainSize>& DMat,
      Eigen::MatrixX<ScalarType>& LMat,
      const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
    using namespace Dune;
    using namespace Dune::DerivativeDirections;
 
    auto strainFunction      = underlying().strainFunction(par, dx);
    const auto C             = underlying().materialTangentFunction(par);
    const auto geo           = underlying().localView().element().geometry();
    const auto numberOfNodes = underlying().numberOfNodes();
    DMat.setZero();
    LMat.setZero(enhancedStrainSize, underlying().localView().size());
    for (const auto& [gpIndex, gp] : strainFunction.viewOverIntegrationPoints()) {
      const auto M            = easFunction.calcM(gp.position());
      const auto CEval        = C(gpIndex);
      const double detJTimesW = geo.integrationElement(gp.position()) * gp.weight();
      DMat += M.transpose() * CEval * M * detJTimesW;
      for (size_t i = 0U; i < numberOfNodes; ++i) {
        const size_t I = Traits::worlddim * i;
        const auto Bi  = strainFunction.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
        LMat.template block<enhancedStrainSize, Traits::worlddim>(0, I) += M.transpose() * CEval * Bi * detJTimesW;
      }
    }
  }
};
 
auto eas(int numberOfEASParameters = 0) {
  EnhancedAssumedStrainsPre pre(numberOfEASParameters);
 
  return pre;
}
 
} // namespace Ikarus
 
#else
  #error EnhancedAssumedStrains depends on dune-localfefunctions, which is not included
#endif