EshelbianPlasticity.hpp

Go to the documentation of this file.

/**
 * \file EshelbianPlasticity.hpp
 * \brief Eshelbian plasticity interface
 *
 * \brief Problem implementation for mix element for large-strain elasticity
 *
 * For reference on mixed formulation see: \cite gopalakrishnan2012second and
 * \cite cockburn2010new
 *
 * \todo Implementation of plasticity
 */
 
// DO NOT DELETE MIGHT BE USEFULL
// #include <boost/multiprecision/mpfr.hpp>
// using mpfr50 = boost::multiprecision::number<
//     boost::multiprecision::mpfr_float_backend<50>>;
 
// inline mpfr50 exp_hi(const mpfr50 x) { return boost::multiprecision::exp(x);
// } inline mpfr50 log_hi(const mpfr50 x) { return
// boost::multiprecision::log(x); }
 
constexpr int SPACE_DIM = 3;
// constexpr auto A = AssemblyType::BLOCK_SCHUR;
constexpr auto A = AssemblyType::BLOCK_MAT;
 
#include "MatOps.hpp"
 
#ifndef __ESHELBIAN_PLASTICITY_HPP__
#define __ESHELBIAN_PLASTICITY_HPP__
 
#include "AnalyticalExprPython.hpp"
#include <bitset>
 
namespace EshelbianPlasticity {
 
using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
 
struct ContactTree;
 
enum SymmetrySelector { SYMMETRIC, NOT_SYMMETRIC, ANIT_SYMMETRIC };
enum RotSelector { SMALL_ROT, MODERATE_ROT, LARGE_ROT, NO_H1_CONFIGURATION };
enum StretchSelector {
  LINEAR,
  LOG /*linear extension*/,
  LOG_QUADRATIC /*quadratic extension*/,
  STRETCH_SELECTOR_LAST
};
enum EnergyReleaseSelector { GRIFFITH_FORCE, GRIFFITH_SKELETON };
 
using MatrixPtr = boost::shared_ptr<MatrixDouble>;
using VectorPtr = boost::shared_ptr<VectorDouble>;
 
using EntData = EntitiesFieldData::EntData;
using UserDataOperator = ForcesAndSourcesCore::UserDataOperator;
using VolUserDataOperator = VolumeElementForcesAndSourcesCore::UserDataOperator;
using FaceUserDataOperator = FaceElementForcesAndSourcesCore::UserDataOperator;
using EleOnSide = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
using SideEleOp = EleOnSide::UserDataOperator;
 
struct PhysicalEquations;
struct LogarithmicMaterialCurvature;
struct TopologicalData;
struct AuxiliaryLogarithmicStressData;
struct AuxiliaryLogarithmicStressMaterialData;
struct EshelbianCore;
 
struct AnalyticalExprPython;
 
struct MatOpsDataAtIntegrationPts
    : public boost::enable_shared_from_this<MatOpsDataAtIntegrationPts> {
 
  MatOpsDataAtIntegrationPts() {
    matOpsData = MatOps::createMatOpsDataPtr();
    matOpsData->insertCommonData("approxPAtPts");
    matOpsData->insertCommonData("approxP0AtPts");
    matOpsData->insertCommonData("divPAtPts");
    matOpsData->insertCommonData("wL2AtPts");
    matOpsData->insertCommonData("wL2DotAtPts");
    matOpsData->insertCommonData("wL2DotDotAtPts");
    matOpsData->insertCommonData("logStretchTensorAtPts");
    matOpsData->insertCommonData("logStretchTensor0AtPts");
    matOpsData->insertCommonData("logStretchDotTensorAtPts");
    matOpsData->insertCommonData("gradLogStretchDotTensorAtPts");
    matOpsData->insertCommonData("rotAxisAtPts");
    matOpsData->insertCommonData("rotAxisGradAtPts");
    matOpsData->insertCommonData("rotAxisDotAtPts");
    matOpsData->insertCommonData("rotAxisGradDotAtPts");
    matOpsData->insertCommonData("stretchH1AtPts");
    matOpsData->insertCommonData("wGradH1AtPts");
    matOpsData->insertCommonData("stretchTensorAtPts");
    matOpsData->insertCommonData("diffStretchH1AtPts");
    matOpsData->insertCommonData("diffStretchTensorAtPts");
    matOpsData->insertCommonData("hAtPts");
    matOpsData->insertCommonData("adjointPdstretchAtPts");
    matOpsData->insertCommonData("adjointPdUAtPts");
    matOpsData->insertCommonData("adjointPdUdOmegaAtPts");
    matOpsData->insertCommonData("adjointPdUdPAtPts");
    matOpsData->insertCommonData("eigenVals");
    matOpsData->insertCommonData("eigenVecs");
    matOpsData->insertCommonData("PAtPts");
    matOpsData->insertCommonData("PAtPts_du");
    matOpsData->insertCommonData("plasticH");
    matOpsData->insertCommonData("plasticFlow");
    matOpsData->insertCommonData("plasticF");
    matOpsData->insertCommonData("invPlasticF");
  }
 
  boost::shared_ptr<MatOps::MatOpsData> matOpsData;
};
 
struct DataAtIntegrationPts : public MatOpsDataAtIntegrationPts {
 
  boost::shared_ptr<AuxiliaryLogarithmicStressData> auxiliaryData;
  boost::shared_ptr<AuxiliaryLogarithmicStressData> auxiliaryData0;
  boost::shared_ptr<AuxiliaryLogarithmicStressMaterialData>
      auxiliaryMaterialData;
 
  MatrixDouble tractionAtPts;
 
  MatrixDouble wH1AtPts;
  MatrixDouble XH1AtPts;
  MatrixDouble contactL2AtPts;
  MatrixDouble logStretch2H1AtPts;
  MatrixDouble logStretchTotalTensorAtPts;
 
  // Constant-mode and higher-mode covectors share the existing vector DOFs.
  struct RotationLinearization {
    MatrixDouble h, work, previousWork, stretch, piola, stretchRotation;
  };
  RotationLinearization constantRotationLinearization;
  MatrixDouble constantRotationSecond, constantHigherRotationSecond,
      higherConstantRotationSecond;
  MatrixPtr constantRotation = boost::make_shared<MatrixDouble>();
  MatrixPtr higherRotation = boost::make_shared<MatrixDouble>();
  MatrixPtr previousConstantRotation = boost::make_shared<MatrixDouble>();
  MatrixPtr previousHigherRotation = boost::make_shared<MatrixDouble>();
 
  MatrixDouble hdOmegaAtPts;
  MatrixDouble hdLogStretchAtPts;
  MatrixDouble leviKirchhoffAtPts;
  MatrixDouble leviKirchhoff0AtPts;
  MatrixDouble leviKirchhoffdOmegaAtPts;
  MatrixDouble leviKirchhoffdLogStreatchAtPts;
  MatrixDouble leviKirchhoffPAtPts;
  MatrixDouble rotMatAtPts;
  VectorDouble energyAtPts; //< this is density of energy at integration points
 
  MatrixPtr varConstantRotation = boost::make_shared<MatrixDouble>();
  MatrixPtr varHigherRotation = boost::make_shared<MatrixDouble>();
  MatrixDouble varRotAxis;
  MatrixDouble varGradRotAxis;
  MatrixDouble varLogStreach;
  MatrixDouble varPiola;
  MatrixDouble varDivPiola;
  MatrixDouble varWL2;
  MatrixDouble varHybridDispAtPts;
 
  VectorInt nbUniq;
  MatrixDouble eigenValsC;
  MatrixDouble eigenVecsC;
  VectorInt nbUniqC;
 
  MatrixDouble matD;
  MatrixDouble matAxiatorD;
  MatrixDouble matDeviatorD;
  MatrixDouble matInvD;
 
  MatrixDouble internalStressAtPts;
 
  MatrixDouble hybridDispAtPts;
  MatrixDouble gradHybridDispAtPts;
  MatrixDouble faceMaterialForceAtPts;
  VectorDouble normalPressureAtPts;
 
  MatrixDouble gradientAtPts;
 
  VectorDouble muAtPts;
  VectorDouble lambdaAtPts;
  VectorDouble youngModulusAtPts;
  double piolaScale = 1.;
 
  static constexpr int SizeSymm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
 
  template <class Tensor>
  static inline auto getFTensor(MatrixDouble &m, const int nb_gauss_pts) {
    return MatrixSizeHelper<Tensor, DL>::get(m, nb_gauss_pts)();
  }
 
  template <class Tensor>
  static inline auto getFTensor(MatrixPtr m, const int nb_gauss_pts) {
    return getFTensor<Tensor>(*m, nb_gauss_pts);
  }
 
  inline auto getFTensorApproxP(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("approxPAtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorApproxP0(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("approxP0AtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorDivP(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("divPAtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorSmallWL2(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("wL2AtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorSmallWL2Dot(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("wL2DotAtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorSmallWL2DotDot(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("wL2DotDotAtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorLogStretch(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2SymmetricFromMatType<SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("logStretchTensorAtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorLogStretch0(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2SymmetricFromMatType<SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("logStretchTensor0AtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorLogStretchDot(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2SymmetricFromMatType<SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("logStretchDotTensorAtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorGradLogStretchDot(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2FromMatType<SizeSymm, SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("gradLogStretchDotTensorAtPts"),
        nb_gauss_pts);
  }
 
  inline auto getFTensorRotAxis(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("rotAxisAtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorRotAxisGrad(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("rotAxisGradAtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorRotAxisDot(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("rotAxisDotAtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorRotAxisGradDot(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("rotAxisGradDotAtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorStretchH1(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("stretchH1AtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorSmallWGradH1(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("wGradH1AtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorStretch(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2SymmetricFromMatType<SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("stretchTensorAtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorDiffStretchH1(const int nb_gauss_pts) {
    return getFTensor<GetFTensor4FromMatType<SPACE_DIM, SPACE_DIM, SPACE_DIM,
                                             SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("diffStretchH1AtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorDiffStretch(const int nb_gauss_pts) {
    return getFTensor<GetFTensor4DdgFromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("diffStretchTensorAtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorSmallH(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("hAtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorPlasticH(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2SymmetricFromMatType<SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("plasticH"), nb_gauss_pts);
  }
 
  inline auto getFTensorPlasticF(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2SymmetricFromMatType<SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("plasticF"), nb_gauss_pts);
  }
 
  inline auto getFTensorInvPlasticF(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("invPlasticF"), nb_gauss_pts);
  }
 
  inline auto getFTensorAdjointPdstretch(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("adjointPdstretchAtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorAdjointPdU(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SizeSymm, -1, DL>>(
        matOpsData->getCommonDataPtr("adjointPdUAtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorAdjointPdUdOmega(const int nb_gauss_pts,
                                         const int mode = 1,
                                         const int point = 0) {
    return GetFTensor2FromMatImpl<3, SizeSymm, -1, DL, MatrixDouble>::get(
        mode ? *getAdjointPdUdOmegaAtPts()
             : constantRotationLinearization.stretchRotation,
        point, 0);
  }
 
  inline auto getFTensorAdjointPdUdP(const int nb_gauss_pts) {
    return getFTensor<
        GetFTensor3FromMatType<SPACE_DIM, SPACE_DIM, SizeSymm, -1, DL>>(
        matOpsData->getCommonDataPtr("adjointPdUdPAtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorEigenVals(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("eigenVals"), nb_gauss_pts);
  }
 
  inline auto getFTensorEigenVecs(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("eigenVecs"), nb_gauss_pts);
  }
 
  inline MatrixPtr getApproxPAtPts() {
    return matOpsData->getCommonDataPtr("approxPAtPts");
  }
 
  inline MatrixPtr getApproxP0AtPts() {
    return matOpsData->getCommonDataPtr("approxP0AtPts");
  }
 
  inline MatrixPtr getPAtPts() {
    return matOpsData->getCommonDataPtr("PAtPts");
  }
 
  inline MatrixPtr getDiffStretchH1AtPts() {
    return matOpsData->getCommonDataPtr("diffStretchH1AtPts");
  }
 
  inline MatrixPtr getStretchH1AtPts() {
    return matOpsData->getCommonDataPtr("stretchH1AtPts");
  }
 
  inline MatrixPtr getAdjointPdstretchAtPts() {
    return matOpsData->getCommonDataPtr("adjointPdstretchAtPts");
  }
 
  inline MatrixPtr getAdjointPdUAtPts() {
    return matOpsData->getCommonDataPtr("adjointPdUAtPts");
  }
 
  inline MatrixPtr getAdjointPdUdOmegaAtPts() {
    return matOpsData->getCommonDataPtr("adjointPdUdOmegaAtPts");
  }
 
  inline MatrixPtr getAdjointPdUdPAtPts() {
    return matOpsData->getCommonDataPtr("adjointPdUdPAtPts");
  }
 
  inline MatrixPtr getEigenVals() {
    return matOpsData->getCommonDataPtr("eigenVals");
  }
 
  inline MatrixPtr getEigenVecs() {
    return matOpsData->getCommonDataPtr("eigenVecs");
  }
 
  inline MatrixPtr getDivPAtPts() {
    return matOpsData->getCommonDataPtr("divPAtPts");
  }
 
  inline MatrixPtr getSmallWL2AtPts() {
    return matOpsData->getCommonDataPtr("wL2AtPts");
  }
 
  inline MatrixPtr getSmallWL2DotAtPts() {
    return matOpsData->getCommonDataPtr("wL2DotAtPts");
  }
 
  inline MatrixPtr getSmallWL2DotDotAtPts() {
    return matOpsData->getCommonDataPtr("wL2DotDotAtPts");
  }
 
  inline MatrixPtr getLogStretchTensorAtPts() {
    return matOpsData->getCommonDataPtr("logStretchTensorAtPts");
  }
 
  inline MatrixPtr getLogStretchTotalTensorAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &logStretchTotalTensorAtPts);
  }
 
  inline MatrixPtr getLogStretchTensor0AtPts() {
    return matOpsData->getCommonDataPtr("logStretchTensor0AtPts");
  }
 
  inline MatrixPtr getStretchTensorAtPts() {
    return matOpsData->getCommonDataPtr("stretchTensorAtPts");
  }
 
  inline MatrixPtr getDiffStretchTensorAtPts() {
    return matOpsData->getCommonDataPtr("diffStretchTensorAtPts");
  }
 
  inline MatrixPtr getDeformationGradient() {
    return matOpsData->getCommonDataPtr("hAtPts");
  }
 
  inline MatrixPtr getPlasticH() {
    return matOpsData->getCommonDataPtr("plasticH");
  }
 
  inline MatrixPtr getPlasticFlow() {
    return matOpsData->getCommonDataPtr("plasticFlow");
  }
 
  inline MatrixPtr getPlasticF() {
    return matOpsData->getCommonDataPtr("plasticF");
  }
 
  inline MatrixPtr getInvPlasticF() {
    return matOpsData->getCommonDataPtr("invPlasticF");
  }
 
  inline MatrixPtr getLogStretchDotTensorAtPts() {
    return matOpsData->getCommonDataPtr("logStretchDotTensorAtPts");
  }
 
  inline MatrixPtr getGradLogStretchDotTensorAtPts() {
    return matOpsData->getCommonDataPtr("gradLogStretchDotTensorAtPts");
  }
 
  inline MatrixPtr getRotAxisAtPts() {
    return matOpsData->getCommonDataPtr("rotAxisAtPts");
  }
 
  inline MatrixPtr getRotAxisGradAtPts() {
    return matOpsData->getCommonDataPtr("rotAxisGradAtPts");
  }
 
  inline MatrixPtr getRotAxisDotAtPts() {
    return matOpsData->getCommonDataPtr("rotAxisDotAtPts");
  }
 
  inline MatrixPtr getRotAxisGradDotAtPts() {
    return matOpsData->getCommonDataPtr("rotAxisGradDotAtPts");
  }
 
  inline MatrixPtr getMatDPtr() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matD);
  }
 
  inline MatrixPtr getMatInvDPtr() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matInvD);
  }
 
  inline MatrixPtr getSmallWH1AtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wH1AtPts);
  }
 
  inline MatrixPtr getLargeXH1AtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &XH1AtPts);
  }
 
  inline MatrixPtr getContactL2AtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &contactL2AtPts);
  }
 
  inline MatrixPtr getSmallWGradH1AtPts() {
    return matOpsData->getCommonDataPtr("wGradH1AtPts");
  }
 
  inline MatrixPtr getLeviKirchhoffAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &leviKirchhoffAtPts);
  };
 
  inline MatrixPtr getVarRotAxisPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varRotAxis);
  }
 
  inline MatrixPtr getVarGradRotAxisPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varGradRotAxis);
  }
 
  inline MatrixPtr getVarLogStreachPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varLogStreach);
  }
 
  inline MatrixPtr getVarPiolaPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varPiola);
  }
 
  inline MatrixPtr getDivVarPiolaPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varDivPiola);
  }
 
  inline MatrixPtr getVarWL2Pts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varWL2);
  }
 
  inline MatrixPtr getVarHybridDispAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &varHybridDispAtPts);
  }
 
  inline MatrixPtr getInternalStressAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &internalStressAtPts);
  }
 
  inline MatrixPtr getHybridDispAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &hybridDispAtPts);
  }
 
  inline MatrixPtr getGradHybridDispAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &gradHybridDispAtPts);
  }
 
  inline MatrixPtr getGradientAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &gradientAtPts);
  }
 
  inline MatrixPtr getRotMatAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &rotMatAtPts);
  }
 
  inline auto getFTensorEigenValsC(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(eigenValsC,
                                                                 nb_gauss_pts);
  }
 
  inline auto getFTensorEigenVecsC(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        eigenVecsC, nb_gauss_pts);
  }
 
  inline auto getFTensorTraction(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(tractionAtPts,
                                                                 nb_gauss_pts);
  }
 
  inline auto getFTensorLogStretch2H1(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2SymmetricFromMatType<SPACE_DIM, -1, DL>>(
        logStretch2H1AtPts, nb_gauss_pts);
  }
 
  inline auto getFTensorLogStretchTotal(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2SymmetricFromMatType<SPACE_DIM, -1, DL>>(
        logStretchTotalTensorAtPts, nb_gauss_pts);
  }
 
  inline auto getFTensorSmallHdOmega(const int nb_gauss_pts, const int mode = 1,
                                     const int point = 0) {
    return GetFTensor3FromMatImpl<3, 3, 3, -1, DL, MatrixDouble>::get(
        mode ? hdOmegaAtPts : constantRotationLinearization.h, point, 0);
  }
 
  inline auto getFTensorSmallHdLogStretch(const int nb_gauss_pts) {
    return getFTensor<
        GetFTensor3FromMatType<SPACE_DIM, SPACE_DIM, SizeSymm, -1, DL>>(
        hdLogStretchAtPts, nb_gauss_pts);
  }
 
  inline auto getFTensorLeviKirchhoff(const int nb_gauss_pts,
                                      const int mode = 1, const int point = 0) {
    return GetFTensor1FromMatImpl<3, -1, DL, MatrixDouble>::get(
        mode ? leviKirchhoffAtPts : constantRotationLinearization.work, point,
        0);
  }
 
  inline auto getFTensorLeviKirchhoff0(const int nb_gauss_pts,
                                       const int mode = 1,
                                       const int point = 0) {
    return GetFTensor1FromMatImpl<3, -1, DL, MatrixDouble>::get(
        mode ? leviKirchhoff0AtPts : constantRotationLinearization.previousWork,
        point, 0);
  }
 
  inline auto getFTensorLeviKirchhoffdOmega(const int nb_gauss_pts,
                                            const int row_mode = 1,
                                            const int col_mode = 1,
                                            const int point = 0) {
    auto &matrix = row_mode ? (col_mode ? leviKirchhoffdOmegaAtPts
                                        : higherConstantRotationSecond)
                            : (col_mode ? constantHigherRotationSecond
                                        : constantRotationSecond);
    return GetFTensor2FromMatImpl<3, 3, -1, DL, MatrixDouble>::get(matrix,
                                                                   point, 0);
  }
 
  inline auto getFTensorLeviKirchhoffdLogStretch(const int nb_gauss_pts,
                                                 const int mode = 1,
                                                 const int point = 0) {
    return GetFTensor2FromMatImpl<3, SizeSymm, -1, DL, MatrixDouble>::get(
        mode ? leviKirchhoffdLogStreatchAtPts
             : constantRotationLinearization.stretch,
        point, 0);
  }
 
  inline auto getFTensorLeviKirchhoffP(const int nb_gauss_pts,
                                       const int mode = 1,
                                       const int point = 0) {
    return GetFTensor3FromMatImpl<3, 3, 3, -1, DL, MatrixDouble>::get(
        mode ? leviKirchhoffPAtPts : constantRotationLinearization.piola, point,
        0);
  }
 
  inline auto getFTensorRotMat(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        rotMatAtPts, nb_gauss_pts);
  }
 
  inline auto getFTensorVarRotAxis(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(varRotAxis,
                                                                 nb_gauss_pts);
  }
 
  inline auto getFTensorVarGradRotAxis(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        varGradRotAxis, nb_gauss_pts);
  }
 
  inline auto getFTensorVarLogStreach(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SizeSymm, -1, DL>>(varLogStreach,
                                                                nb_gauss_pts);
  }
 
  inline auto getFTensorVarPiola(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        varPiola, nb_gauss_pts);
  }
 
  inline auto getFTensorDivVarPiola(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(varDivPiola,
                                                                 nb_gauss_pts);
  }
 
  inline auto getFTensorVarWL2(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(varWL2,
                                                                 nb_gauss_pts);
  }
 
  inline auto getFTensorInternalStress(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        internalStressAtPts, nb_gauss_pts);
  }
 
  inline auto getFTensorInternalStressVec(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM * SPACE_DIM, -1, DL>>(
        internalStressAtPts, nb_gauss_pts);
  }
 
  inline auto getFTensorSmallWH1(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(wH1AtPts,
                                                                 nb_gauss_pts);
  }
 
  inline auto getFTensorSmallHybridDisp(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(
        hybridDispAtPts, nb_gauss_pts);
  }
 
  inline auto getFTensorGradHybridDisp(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        gradHybridDispAtPts, nb_gauss_pts);
  }
 
  inline auto getFTensorFaceMaterialForce(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(
        faceMaterialForceAtPts, nb_gauss_pts);
  }
 
  boost::shared_ptr<PhysicalEquations> physicsPtr;
};
 
struct OpJacobian;
 
// Forward declarations
struct ExternalStrain;
using ExternalStrainVec = std::vector<ExternalStrain>;
 
struct MaterialPostProcData {
  using Map = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
  Map::DataMapVec scalarFields;
  Map::DataMapMat vectorFields;
  Map::DataMapMat symmetricFields;
};
 
struct PhysicalEquations {
 
  enum MaterialModel {
    StVenantKirchhoff,
    MooneyRivlin,
    Hencky,
    Neohookean,
    Meta, //< Metamaterial model, that is using off the shelf MatOps
          // functionality
    LastMaterial
  };
 
  enum Feature {
    DIRECT_STRETCH = 0,           ///< Direct stretch formulation
    NO_STRETCH_LINEAR,            ///< No-stretch linear formulation
    NO_STRETCH_NONLINEAR,         ///< No-stretch nonlinear formulation
    AUXILIARY_LOGARITHMIC_STRESS, ///< Auxiliary logarithmic stress formulation
    NON_HOMOGENEOUS_MATERIAL,
    LAST_FEATURE
  };
 
  using Features = std::bitset<LAST_FEATURE>;
  inline static const Features noStretchMask =
      Features{}.set(NO_STRETCH_LINEAR).set(NO_STRETCH_NONLINEAR);
 
  enum FieldState {
    CURRENT, ///< Current state of the field
    PREVIOUS ///< Previous state of the field
  };
 
  PhysicalEquations(MaterialModel model, Features features)
      : materialModel(model), materialFeatures(features) {}
  virtual ~PhysicalEquations() = default;
 
  /// Select and configure the material before registering fields and DMs.
  static MoFEMErrorCode create(EshelbianCore &ep);
 
  /// Check material capabilities after geometry and boundary setup.
  MoFEMErrorCode checkSetup(const EshelbianCore &ep) const;
 
  MaterialModel getMaterialModel() const { return materialModel; }
  /// Test individual bits or compare (getFeatures() & mask) with a feature
  /// mask.
  const Features &getFeatures() const { return materialFeatures; }
 
  struct FieldDefinition {
    std::string name;
    int order;
    int coefficients;
  };
 
  /// USER_BASE L2 fields; negative order retains an inactive field
  /// registration.
  std::vector<FieldDefinition>
  getMaterialFieldDefinitions(const EshelbianCore &ep) const;
 
  /// Active material unknowns used by the element, DMs and field split.
  std::vector<std::string> getMaterialFields(const EshelbianCore &ep) const;
 
  /// Material-dependent blocks; each pair also represents its transpose.
  std::vector<std::pair<std::string, std::string>>
  getMaterialCouplings(const EshelbianCore &ep) const;
  std::vector<std::pair<std::string, std::string>>
  getMaterialEmptyBlocks(const EshelbianCore &ep) const;
 
  /// Evaluate the selected fields, or recover stretch from the matching stress.
  MoFEMErrorCode
  pushMaterialFields(const EshelbianCore &ep,
                     boost::ptr_deque<UserDataOperator> &pipeline,
                     boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                     FieldState state = CURRENT);
 
  MoFEMErrorCode
  pushMaterialRates(const EshelbianCore &ep,
                    boost::ptr_deque<UserDataOperator> &pipeline,
                    boost::shared_ptr<DataAtIntegrationPts> data_ptr);
 
  /// Evaluate constitutive data after kinematics and before material assembly.
  MoFEMErrorCode pushMaterialEvaluation(
      const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &pipeline,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr, bool rhs, bool lhs);
 
  /// Contract material curvature for the incremental material-energy adjoint.
  virtual bool providesHelmholtzCurvature() const { return false; }
 
  virtual MoFEMErrorCode
  evaluateHelmholtzCurvature(EntityHandle entity, int gauss_point,
                             DataAtIntegrationPts &data,
                             LogarithmicMaterialCurvature &curvature);
 
  /// Pure material Hessian in five orthonormal deviatoric log-strain
  /// coordinates. Evaluated at the current elastic base, without any adjoint
  /// contraction.
  virtual MoFEMErrorCode evaluateElasticLogarithmicMaterialHessian(
      EntityHandle entity, int gauss_point, DataAtIntegrationPts &data,
      const FTensor::Tensor2_symmetric<double, 3> &elastic_log_stretch,
      FTensor::Tensor2<double, 5, 5> &hessian);
 
  /// Assemble the physical Helmholtz state derivative at equilibrium.
  virtual MoFEMErrorCode
  pushHelmholtzStateGradient(const EshelbianCore &ep,
                             boost::ptr_deque<UserDataOperator> &pipeline,
                             boost::shared_ptr<DataAtIntegrationPts> data_ptr);
 
  MoFEMErrorCode
  pushStretchFromStress(const EshelbianCore &ep,
                        boost::ptr_deque<UserDataOperator> &pipeline,
                        boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                        boost::shared_ptr<MatrixDouble> strain_ptr = nullptr,
                        boost::shared_ptr<MatrixDouble> stress_ptr = nullptr);
 
  MoFEMErrorCode
  pushMaterialResidual(const EshelbianCore &ep,
                       boost::ptr_deque<UserDataOperator> &pipeline,
                       boost::shared_ptr<DataAtIntegrationPts> data_ptr);
 
  MoFEMErrorCode
  pushMaterialTangent(const EshelbianCore &ep,
                      boost::ptr_deque<UserDataOperator> &pipeline,
                      boost::shared_ptr<DataAtIntegrationPts> data_ptr);
 
  MoFEMErrorCode
  pushEnergyEvaluation(const EshelbianCore &ep,
                       boost::ptr_deque<UserDataOperator> &pipeline,
                       boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                       boost::shared_ptr<double> total_energy_ptr);
 
  MoFEMErrorCode
  pushMaterialForceFields(const EshelbianCore &ep,
                          boost::ptr_deque<UserDataOperator> &pipeline,
                          boost::shared_ptr<DataAtIntegrationPts> data_ptr);
 
  MoFEMErrorCode pushPostProc(const EshelbianCore &ep,
                              boost::ptr_deque<UserDataOperator> &pipeline,
                              boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                              MaterialPostProcData &output);
 
  MoFEMErrorCode pushMaterialVariation(
      const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &pipeline,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      SmartPetscObj<Vec> variation, MaterialPostProcData *output = nullptr);
 
  MoFEMErrorCode pushPostProcResidual(
      const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &pipeline,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      SmartPetscObj<Vec> residual, MaterialPostProcData &output);
 
  MoFEMErrorCode pushSkeletonEvaluation(
      const EshelbianCore &ep, VolumeElementForcesAndSourcesCoreOnSide &fe,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      SmartPetscObj<Vec> residual, MaterialPostProcData &output);
 
  /// Assemble the selected material field group after kinematic reconstruction.
  virtual MoFEMErrorCode pushAuxiliaryLogarithmicStressOps(
      const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &pipeline,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr, bool lhs);
 
  /// Evaluate the auxiliary material copy and energies without assembly.
  virtual MoFEMErrorCode pushAuxiliaryLogarithmicMaterialEvaluation(
      const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &pipeline,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr, bool lhs = false);
 
  virtual UserDataOperator *
  returnOpJacobian(const bool eval_rhs, const bool eval_lhs,
                   boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                   boost::shared_ptr<PhysicalEquations> physics_ptr);
 
  virtual VolUserDataOperator *
  returnOpSpatialPhysical(const std::string &field_name,
                          boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                          const double alpha_u);
 
  virtual VolUserDataOperator *returnOpSpatialPhysicalExternalStrain(
      const std::string &field_name,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv);
 
  virtual VolUserDataOperator *returnOpSpatialPhysical_du_du(
      std::string row_field, std::string col_field,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha);
 
  virtual VolUserDataOperator *
  returnOpTopoSpatialPhysical(const std::string &field_name,
                              boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                              SmartPetscObj<Vec> assemble_vec,
                              boost::shared_ptr<TopologicalData> topo_ptr,
                              const double alpha_u,
                              boost::shared_ptr<double> J_ptr);
 
  virtual VolUserDataOperator *returnOpCalculateHelmholtzFreeEnergy(
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<double> total_helmholtz_free_energy_ptr);
 
  /** True when the material supplies its Helmholtz free energy. */
  virtual bool providesHelmholtzFreeEnergy() const;
 
  virtual VolUserDataOperator *returnOpCalculateStretchFromStress(
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<PhysicalEquations> physics_ptr,
      boost::shared_ptr<MatrixDouble> strain_ptr = nullptr);
 
  virtual VolUserDataOperator *returnOpCalculateStretchFromStress(
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<PhysicalEquations> physics_ptr,
      boost::shared_ptr<MatrixDouble> strain_ptr,
      VectorPtr external_pressure_ptr);
 
  virtual VolUserDataOperator *returnOpCalculateStretchFromStress(
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<PhysicalEquations> physics_ptr,
      boost::shared_ptr<MatrixDouble> strain_ptr,
      boost::shared_ptr<MatrixDouble> stress_ptr,
      VectorPtr external_pressure_ptr);
 
  virtual VolUserDataOperator *returnOpCalculateExternalPressure(
      VectorPtr external_pressure_ptr,
      boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv);
 
  virtual VolUserDataOperator *returnOpCalculateVarStretchFromStress(
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<PhysicalEquations> physics_ptr);
 
private:
  bool usesStressBasedStrain() const {
    return (materialFeatures & noStretchMask).any() || materialModel == Hencky;
  }
 
  static MoFEMErrorCode checkOptions(const EshelbianCore &ep,
                                     MaterialModel material_model,
                                     Features features);
  MaterialModel materialModel;
 
protected:
  Features materialFeatures;
};
 
struct BcDisp {
  BcDisp(std::string name, std::vector<double> attr, Range faces,
         std::string load_history_file = "");
  std::string blockName;
  std::string loadHistoryFile;
  Range faces;
  VectorDouble3 vals;
  VectorInt3 flags;
};
using BcDispVec = std::vector<BcDisp>;
 
struct BcRot {
  BcRot(std::string name, std::vector<double> attr, Range faces,
        std::string load_history_file = "");
  std::string blockName;
  std::string loadHistoryFile;
  Range faces;
  VectorDouble vals;
  double theta;
};
using BcRotVec = std::vector<BcRot>;
 
typedef std::vector<Range> TractionFreeBc;
 
struct TractionBc {
  TractionBc(std::string name, std::vector<double> attr, Range faces,
             std::string load_history_file = "");
  std::string blockName;
  std::string loadHistoryFile;
  Range faces;
  VectorDouble3 vals;
  VectorInt3 flags;
};
using TractionBcVec = std::vector<TractionBc>;
 
struct NormalDisplacementBc {
  NormalDisplacementBc(std::string name, std::vector<double> attr, Range faces,
                       std::string load_history_file = "");
  std::string blockName;
  std::string loadHistoryFile;
  Range faces;
  double val;
};
using NormalDisplacementBcVec = std::vector<NormalDisplacementBc>;
 
struct SpringBc {
  SpringBc(std::string name, std::vector<double> attr, Range faces);
  std::string blockName;
  Range faces;
  double normalStiffness;
  double tangentialStiffness;
};
using SpringBcVec = std::vector<SpringBc>;
 
struct AnalyticalDisplacementBc {
  AnalyticalDisplacementBc(std::string name, std::vector<double> attr,
                           Range faces, std::string load_history_file = "");
  std::string blockName;
  Range faces;
  VectorInt3 flags;
};
using AnalyticalDisplacementBcVec = std::vector<AnalyticalDisplacementBc>;
 
struct AnalyticalTractionBc {
  AnalyticalTractionBc(std::string name, std::vector<double> attr, Range faces,
                       std::string load_history_file = "");
  std::string blockName;
  Range faces;
  VectorInt3 flags;
};
using AnalyticalTractionBcVec = std::vector<AnalyticalTractionBc>;
 
struct PressureBc {
  PressureBc(std::string name, std::vector<double> attr, Range faces,
             std::string load_history_file = "");
  std::string blockName;
  std::string loadHistoryFile;
  Range faces;
  double val;
};
using PressureBcVec = std::vector<PressureBc>;
 
struct ExternalStrain {
  ExternalStrain(std::string name, std::vector<double> attr, Range ents,
                 std::string load_history_file = "");
  std::string blockName;
  std::string loadHistoryFile;
  Range ents;
  double val;
  double bulkModulusK;
};
 
template <typename OP_PTR>
std::tuple<std::string, MatrixDouble>
getAnalyticalExpr(OP_PTR op_ptr, MatrixDouble &analytical_expr,
                  const std::string block_name);
 
template <typename OP_PTR>
std::tuple<std::string, VectorDouble>
getAnalyticalExternalStrain(OP_PTR op_ptr, VectorDouble &analytical_expr,
                            const std::string block_name);
 
template <typename OP_PTR>
VectorDouble getAnalyticalElastic(OP_PTR op_ptr, const std::string block_name);
 
MatrixDouble analytical_expr_function(double delta_t, double t,
                                      int nb_gauss_pts,
                                      MatrixDouble &m_ref_coords,
                                      MatrixDouble &m_ref_normals,
                                      const std::string block_name);
 
VectorDouble analytical_externalstrain_function(double delta_t, double t,
                                                int nb_gauss_pts,
                                                MatrixDouble &m_ref_coords,
                                                const std::string block_name);
 
VectorDouble analytical_elastic_function(double delta_t, double t,
                                         int nb_gauss_pts,
                                         MatrixDouble &m_ref_coords,
                                         const std::string block_name);
 
#include "CGGUserPolynomialBase.hpp"
#include "EshelbianCore.hpp"
#include "EshelbianOperators.hpp"
 
} // namespace EshelbianPlasticity
 
#endif //__ESHELBIAN_PLASTICITY_HPP__