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"
 
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 TopologicalData;
 
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("rotAxis0AtPts");
    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");
  }
 
  boost::shared_ptr<MatOps::MatOpsData> matOpsData;
};
 
struct DataAtIntegrationPts : public MatOpsDataAtIntegrationPts {
 
  MatrixDouble tractionAtPts;
 
  MatrixDouble wH1AtPts;
  MatrixDouble XH1AtPts;
  MatrixDouble contactL2AtPts;
  MatrixDouble logStretch2H1AtPts;
  MatrixDouble logStretchTotalTensorAtPts;
 
  MatrixDouble hdOmegaAtPts;
  MatrixDouble hdLogStretchAtPts;
  MatrixDouble leviKirchhoffAtPts;
  MatrixDouble leviKirchhoff0AtPts;
  MatrixDouble leviKirchhoffdOmegaAtPts;
  MatrixDouble leviKirchhoffdLogStreatchAtPts;
  MatrixDouble leviKirchhoffPAtPts;
  MatrixDouble rotMatAtPts;
  MatrixDouble SigmaAtPts;
  VectorDouble energyAtPts; //< this is density of energy at integration points
 
  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;
 
  double mu;
  double lambda;
  VectorDouble muAtPts;
  VectorDouble lambdaAtPts;
  VectorDouble youngModulusAtPts;
  double piolaScale = 1.;
 
  inline auto getPiolaScalePtr() {
    return boost::shared_ptr<double>(shared_from_this(), &piolaScale);
  }
 
  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 getFTensorP(const int nb_gauss_pts) {
    return getFTensor<
        GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("PAtPts"), nb_gauss_pts);
  }
 
  inline auto getFTensorPDu(const int nb_gauss_pts) {
    return getFTensor<GetFTensor4FromMatType<
        SPACE_DIM, SPACE_DIM, SPACE_DIM, SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("PAtPts_du"), 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 getFTensorRotAxis0(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(
        matOpsData->getCommonDataPtr("rotAxis0AtPts"), 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 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) {
    return getFTensor<GetFTensor2FromMatType<SPACE_DIM, SizeSymm, -1, DL>>(
        matOpsData->getCommonDataPtr("adjointPdUdOmegaAtPts"),
        nb_gauss_pts);
  }
 
  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 getPAtPtsDu() {
    return matOpsData->getCommonDataPtr("PAtPts_du");
  }
 
  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 getLogStretchDotTensorAtPts() {
    return matOpsData->getCommonDataPtr("logStretchDotTensorAtPts");
  }
 
  inline MatrixPtr getGradLogStretchDotTensorAtPts() {
    return matOpsData->getCommonDataPtr("gradLogStretchDotTensorAtPts");
  }
 
  inline MatrixPtr getRotAxisAtPts() {
    return matOpsData->getCommonDataPtr("rotAxisAtPts");
  }
 
  inline MatrixPtr getRotAxis0AtPts() {
    return matOpsData->getCommonDataPtr("rotAxis0AtPts");
  }
 
  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 getMatAxiatorDPtr() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matAxiatorD);
  }
 
  inline MatrixPtr getMatDeviatorDPtr() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matDeviatorD);
  }
 
  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) {
    return getFTensor<GetFTensor3FromMatType<SPACE_DIM, SPACE_DIM, SPACE_DIM,
                                            -1, DL>>(hdOmegaAtPts,
                                                     nb_gauss_pts);
  }
 
  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) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(
        leviKirchhoffAtPts, nb_gauss_pts);
  }
 
  inline auto getFTensorLeviKirchhoff0(const int nb_gauss_pts) {
    return getFTensor<GetFTensor1FromMatType<SPACE_DIM, -1, DL>>(
        leviKirchhoff0AtPts, nb_gauss_pts);
  }
 
  inline auto getFTensorLeviKirchhoffdOmega(const int nb_gauss_pts) {
    return getFTensor<
        GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        leviKirchhoffdOmegaAtPts, nb_gauss_pts);
  }
 
  inline auto getFTensorLeviKirchhoffdLogStretch(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2FromMatType<SPACE_DIM, SizeSymm, -1, DL>>(
        leviKirchhoffdLogStreatchAtPts, nb_gauss_pts);
  }
 
  inline auto getFTensorLeviKirchhoffP(const int nb_gauss_pts) {
    return getFTensor<GetFTensor3FromMatType<SPACE_DIM, SPACE_DIM, SPACE_DIM,
                                            -1, DL>>(leviKirchhoffPAtPts,
                                                     nb_gauss_pts);
  }
 
  inline auto getFTensorRotMat(const int nb_gauss_pts) {
    return getFTensor<
        GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(rotMatAtPts,
                                                              nb_gauss_pts);
  }
 
  inline auto getFTensorSigma(const int nb_gauss_pts) {
    return getFTensor<
        GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(SigmaAtPts,
                                                              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);
  }
 
  inline auto getFTensorGradient(const int nb_gauss_pts) {
    return getFTensor<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>>(
        gradientAtPts, nb_gauss_pts);
  }
 
  boost::shared_ptr<PhysicalEquations> physicsPtr;
};
 
struct OpJacobian;
 
// Forward declarations
struct ExternalStrain;
using ExternalStrainVec = std::vector<ExternalStrain>;
 
struct PhysicalEquations {
 
  PhysicalEquations() = default;
  virtual ~PhysicalEquations() = default;
 
  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 *
  returnOpCalculateEnergy(boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                          boost::shared_ptr<double> total_energy_ptr);
 
  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 *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);
};
 
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;
};
 
/**
 * @brief Push pointwise external-pressure evaluation before stress recovery.
 *
 * Materials which do not use an external pressure in their local recovery
 * return a null pressure operator and retain their existing behaviour.
 */
template <typename OpVector>
inline void pushOpCalculateStretchFromStress(
    OpVector &op_vector, boost::shared_ptr<PhysicalEquations> physics_ptr,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
    boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
    const std::map<std::string, boost::shared_ptr<ScalingMethod>> &smv,
    boost::shared_ptr<MatrixDouble> strain_ptr = nullptr) {
  auto external_pressure_ptr = boost::make_shared<VectorDouble>();
  if (auto op = physics_ptr->returnOpCalculateExternalPressure(
          external_pressure_ptr, std::move(external_strain_vec_ptr), smv)) {
    op_vector.push_back(op);
  }
  op_vector.push_back(physics_ptr->returnOpCalculateStretchFromStress(
      std::move(data_ptr), physics_ptr, std::move(strain_ptr),
      std::move(external_pressure_ptr)));
}
 
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__