EshelbianMat.cpp

Go to the documentation of this file.

/**
 * \file EshelbianMat.cpp
 * \brief Implementation of automatic differentiation
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <MatElastic.hpp>
#include <MatrixFunction.hpp>
 
#include <AuxiliaryLogarithmicStressOperators.hpp>
#include <AuxiliaryLogarithmicStressPostProc.hpp>
#include <EshelbianAux.hpp>
#include <EshelbianPlasticity.hpp>
#include <PlasticIncrementalOptimization.hpp>
 
#include <EshelbianTopologicalDerivative.hpp>
 
#include <boost/math/constants/constants.hpp>
#include <regex>
 
namespace EshelbianPlasticity {
 
std::vector<PhysicalEquations::FieldDefinition>
PhysicalEquations::getMaterialFieldDefinitions(const EshelbianCore &ep) const {
  std::vector<FieldDefinition> fields;
  auto set_fields = [&]<FieldApproximationBase Base>() {
    MoFEMFunctionBegin;
    using Orders = EshelbianCore::FieldOrders<Base>;
    if (getFeatures().test(AUXILIARY_LOGARITHMIC_STRESS)) {
      fields = {{ep.logDeviator, Orders::logDeviator(ep.spaceOrder), 5},
                {ep.logJacobian, Orders::logJacobian(ep.spaceOrder), 1},
                {ep.auxiliaryLogStress,
                 Orders::auxiliaryLogStress(ep.spaceOrder), 5}};
    } else {
      fields = {{ep.stretchTensor,
                 (getFeatures() & noStretchMask).any()
                     ? -1
                     : Orders::stretch(ep.spaceOrder),
                 6}};
    }
    MoFEMFunctionReturn(0);
  };
  CHK_THROW_MESSAGE(ep.withFieldOrders(set_fields),
                    "Get material field orders failed");
  return fields;
}
 
std::vector<std::string>
PhysicalEquations::getMaterialFields(const EshelbianCore &ep) const {
  std::vector<std::string> fields;
  for (const auto &field : getMaterialFieldDefinitions(ep))
    if (field.order >= 0)
      fields.push_back(field.name);
  return fields;
}
 
std::vector<std::pair<std::string, std::string>>
PhysicalEquations::getMaterialCouplings(const EshelbianCore &ep) const {
  std::vector<std::pair<std::string, std::string>> blocks;
  // Eliminated stretch and the updated plastic flow direction couple P/bubble.
  if ((getFeatures() & noStretchMask).any() || ep.plasticVolume) {
    blocks.emplace_back(ep.bubbleField, ep.bubbleField);
    blocks.emplace_back(ep.bubbleField, ep.piolaStress);
  }
  for (const auto &field : getMaterialFields(ep))
    blocks.emplace_back(field, field);
  if (getFeatures().test(AUXILIARY_LOGARITHMIC_STRESS)) {
    blocks.emplace_back(ep.logDeviator, ep.logJacobian);
    blocks.emplace_back(ep.logDeviator, ep.auxiliaryLogStress);
    for (const auto &field : {ep.logDeviator, ep.logJacobian})
      for (const auto &mechanical :
           {ep.piolaStress, ep.bubbleField, ep.rotAxis})
        blocks.emplace_back(field, mechanical);
  } else if (getFeatures().test(DIRECT_STRETCH)) {
    for (const auto &field : {ep.piolaStress, ep.bubbleField, ep.rotAxis})
      blocks.emplace_back(ep.stretchTensor, field);
  }
  return blocks;
}
 
std::vector<std::pair<std::string, std::string>>
PhysicalEquations::getMaterialEmptyBlocks(const EshelbianCore &ep) const {
  std::vector<std::pair<std::string, std::string>> blocks;
  for (const auto &field : getMaterialFields(ep))
    blocks.emplace_back(ep.spatialL2Disp, field);
  if ((getFeatures() & noStretchMask).none()) {
    blocks.emplace_back(ep.bubbleField, ep.bubbleField);
    blocks.emplace_back(ep.piolaStress, ep.piolaStress);
    blocks.emplace_back(ep.bubbleField, ep.piolaStress);
  }
  if (getFeatures().test(AUXILIARY_LOGARITHMIC_STRESS))
    for (const auto &field :
         {ep.logJacobian, ep.piolaStress, ep.bubbleField, ep.rotAxis})
      blocks.emplace_back(ep.auxiliaryLogStress, field);
  return blocks;
}
 
MoFEMErrorCode
PhysicalEquations::checkOptions(const EshelbianCore &ep,
                                const MaterialModel material_model,
                                const Features features) {
  MoFEMFunctionBegin;
  if (!features.test(AUXILIARY_LOGARITHMIC_STRESS))
    MoFEMFunctionReturnHot(0);
 
  double alpha_grad_u = 0;
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, PETSC_NULLPTR,
                               "-neo_hookean_viscosity_alpha_grad_u",
                               &alpha_grad_u, PETSC_NULLPTR);
  PetscBool h1_update = ep.physicalH1Update;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-physical_h1_update", &h1_update, PETSC_NULLPTR);
  const bool plastic_solver =
      ep.plasticVolume &&
      (ep.solverType == EshelbianCore::IncrementalOptimization ||
       ep.solverType == EshelbianCore::TestEquilibratedMechanicalValue ||
       ep.solverType == EshelbianCore::TestIncrementalOptimizationLayout ||
       ep.solverType == EshelbianCore::TestIncrementalOptimizationTransaction ||
       ep.solverType ==
           EshelbianCore::TestIncrementalOptimizationObjectiveDerivative ||
       ep.solverType ==
           EshelbianCore::TestIncrementalOptimizationConstraintDerivative);
  const std::pair<bool, const char *> unsupported[] = {
      {material_model != Neohookean, "requires -material neo_hookean"},
      {ep.stretchSelector != LOG, "requires -stretches log"},
      {ep.rotSelector != LARGE_ROT, "requires -rotations large"},
      {ep.symmetrySelector != SYMMETRIC,
       "requires the symmetric row convention"},
      {ep.gradApproximator != NO_H1_CONFIGURATION, "requires -grad no_h1"},
      {alpha_grad_u != 0,
       "requires zero -neo_hookean_viscosity_alpha_grad_u"},
      {!plastic_solver &&
           (ep.solverType != EshelbianCore::TimeSolver || ep.physicalTimeFlg),
       "requires time_solver without dynamic relaxation or plastic "
       "incremental optimization"},
      {h1_update == PETSC_TRUE, "does not yet support -physical_h1_update"},
      {ep.crackingOn || ep.interfaceCrack,
       "does not yet support -cracking_on or -cohesive_interface_on"},
      {!ep.internalStressTagName.empty(),
       "does not yet support internal stress"},
      {ep.meshTransferSourceMeshFileSpecified == PETSC_TRUE,
       "does not yet support material mesh transfer"}};
  for (const auto &[invalid, reason] : unsupported) {
    if (invalid)
      SETERRQ(ep.mField.get_comm(), MOFEM_NOT_IMPLEMENTED,
              "-auxiliary_logarithmic_stress %s", reason);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhysicalEquations::checkSetup(const EshelbianCore &ep) const {
  MoFEMFunctionBegin;
  if (!getFeatures().test(AUXILIARY_LOGARITHMIC_STRESS))
    MoFEMFunctionReturnHot(0);
 
  Range sdf_contact_faces;
  for (const auto &block :
       ep.mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
           std::regex("CONTACT_SDF(.*)"))) {
    CHKERR block->getMeshsetIdEntitiesByDimension(
        ep.mField.get_moab(), SPACE_DIM - 1, sdf_contact_faces, true);
  }
  const auto nonempty = [](const auto &ptr) { return ptr && !ptr->empty(); };
  const std::pair<bool, const char *> unsupported_data[] = {
      {nonempty(ep.externalStrainVecPtr), "external strain prescriptions"},
      {nonempty(ep.bcSpatialPressureVecPtr), "follower pressure loads"},
      {nonempty(ep.contactFaces) || !sdf_contact_faces.empty(),
       "contact meshsets"},
      {nonempty(ep.frontEdges) || nonempty(ep.crackFaces) ||
           nonempty(ep.interfaceFaces),
       "front, crack or interface meshsets"}};
  for (const auto &[local_active, capability] : unsupported_data) {
    int active = local_active;
    CHKERR MPI_Allreduce(MPI_IN_PLACE, &active, 1, MPI_INT, MPI_MAX,
                         ep.mField.get_comm());
    if (active)
      SETERRQ(ep.mField.get_comm(), MOFEM_NOT_IMPLEMENTED,
              "-auxiliary_logarithmic_stress does not yet support %s",
              capability);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhysicalEquations::pushMaterialFields(
    const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &pipeline,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr, const FieldState state) {
  MoFEMFunctionBegin;
  data_ptr->physicsPtr = ep.physicalEquations;
  const bool previous = state == PREVIOUS;
  auto strain_ptr = previous ? data_ptr->getLogStretchTensor0AtPts()
                             : data_ptr->getLogStretchTensorAtPts();
  auto values = previous ? ep.solTSStep : SmartPetscObj<Vec>();
  if (getFeatures().test(AUXILIARY_LOGARITHMIC_STRESS)) {
    auto &fields =
        previous ? data_ptr->auxiliaryData0 : data_ptr->auxiliaryData;
    if (!fields)
      fields = boost::make_shared<AuxiliaryLogarithmicStressData>();
    CHKERR pushAuxiliaryLogarithmicFields(ep, pipeline, fields, strain_ptr,
                                          values);
  } else if ((getFeatures() & noStretchMask).any()) {
    CHKERR pushStretchFromStress(ep, pipeline, data_ptr, strain_ptr,
                                 previous ? data_ptr->getApproxP0AtPts()
                                          : data_ptr->getApproxPAtPts());
  } else {
    pipeline.push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
        ep.stretchTensor, strain_ptr, values, MBTET));
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhysicalEquations::pushMaterialRates(
    const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &pipeline,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr) {
  MoFEMFunctionBegin;
  if (getFeatures().test(DIRECT_STRETCH)) {
    pipeline.push_back(new OpCalculateTensor2SymmetricFieldValuesDot<3>(
        ep.stretchTensor, data_ptr->getLogStretchDotTensorAtPts(), MBTET));
    pipeline.push_back(new OpCalculateVectorFieldGradientDot<6, 3>(
        ep.stretchTensor, data_ptr->getGradLogStretchDotTensorAtPts(), MBTET));
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhysicalEquations::pushStretchFromStress(
    const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &pipeline,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
    boost::shared_ptr<MatrixDouble> strain_ptr,
    boost::shared_ptr<MatrixDouble> stress_ptr) {
  MoFEMFunctionBegin;
  data_ptr->physicsPtr = ep.physicalEquations;
  auto external_pressure_ptr = boost::make_shared<VectorDouble>();
  if (auto op = returnOpCalculateExternalPressure(
          external_pressure_ptr, ep.externalStrainVecPtr, ep.timeScaleMap))
    pipeline.push_back(op);
  pipeline.push_back(returnOpCalculateStretchFromStress(
      std::move(data_ptr), ep.physicalEquations, std::move(strain_ptr),
      std::move(stress_ptr), std::move(external_pressure_ptr)));
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhysicalEquations::pushMaterialEvaluation(
    const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &pipeline,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr, const bool rhs,
    const bool lhs) {
  MoFEMFunctionBegin;
  if (getFeatures().test(AUXILIARY_LOGARITHMIC_STRESS)) {
    if (rhs || lhs)
      CHKERR pushAuxiliaryLogarithmicMaterialEvaluation(ep, pipeline, data_ptr,
                                                        lhs);
  } else if (getFeatures().test(DIRECT_STRETCH)) {
    pipeline.push_back(
        returnOpJacobian(rhs, lhs, data_ptr, ep.physicalEquations));
  }
  MoFEMFunctionReturn(0);
}
 
/** Material Helmholtz derivative with respect to the explicit elastic state. */
struct OpCalculateHelmholtzStateGradient : public VolUserDataOperator {
  OpCalculateHelmholtzStateGradient(
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<MatrixDouble> state_gradient_ptr)
      : VolUserDataOperator(NOSPACE, OPSPACE), dataAtPts(std::move(data_ptr)),
        stateGradientPtr(std::move(state_gradient_ptr)) {
    if (!dataAtPts || !stateGradientPtr)
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "Helmholtz state-gradient operator has a null data pointer");
  }
 
  MoFEMErrorCode doWork(int, EntityType, EntData &) override {
    MoFEMFunctionBegin;
    constexpr int state_size = 6;
    const int nb_integration_pts = getGaussPts().size2();
    using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
    auto get_state_gradient =
        MatrixSizeHelper<GetFTensor1FromMatType<state_size, state_size, DL>,
                         DL>::size(*stateGradientPtr, nb_integration_pts);
    auto t_state_gradient = get_state_gradient();
    auto t_work = dataAtPts->getFTensorAdjointPdU(nb_integration_pts);
    auto t_plastic_f = dataAtPts->getFTensorPlasticF(nb_integration_pts);
    FTENSOR_INDEX(state_size, L);
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      // At equilibrium the direct stretch row identifies this assembled
      // stress-work derivative with the material gradient. Jp=1 for the
      // trace-free plastic approximation.
      const double det_plastic_f = determinantTensor3by3(t_plastic_f);
      t_state_gradient(L) = det_plastic_f * t_work(L);
      ++t_state_gradient;
      ++t_work;
      ++t_plastic_f;
    }
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
  boost::shared_ptr<MatrixDouble> stateGradientPtr;
};
 
MoFEMErrorCode PhysicalEquations::pushHelmholtzStateGradient(
    const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &pipeline,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr) {
  MoFEMFunctionBegin;
  auto gradient_at_points = boost::make_shared<MatrixDouble>();
  pipeline.push_back(
      new OpCalculateHelmholtzStateGradient(data_ptr, gradient_at_points));
  using OpAssembleGradient = FormsIntegrators<VolUserDataOperator>::Assembly<
      A>::LinearForm<GAUSS>::OpBaseTimesVector<1, 6, 6>;
  pipeline.push_back(
      new OpAssembleGradient(ep.stretchTensor, gradient_at_points));
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhysicalEquations::pushAuxiliaryLogarithmicStressOps(
    const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &,
    boost::shared_ptr<DataAtIntegrationPts>, bool) {
  MoFEMFunctionBegin;
  SETERRQ(ep.mField.get_comm(), MOFEM_NOT_IMPLEMENTED,
          "Selected material does not implement auxiliary logarithmic stress");
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhysicalEquations::pushAuxiliaryLogarithmicMaterialEvaluation(
    const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &,
    boost::shared_ptr<DataAtIntegrationPts>, bool) {
  MoFEMFunctionBegin;
  SETERRQ(ep.mField.get_comm(), MOFEM_NOT_IMPLEMENTED,
          "Selected material does not evaluate auxiliary logarithmic stress");
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhysicalEquations::pushMaterialResidual(
    const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &pipeline,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr) {
  MoFEMFunctionBegin;
  if ((getFeatures() & noStretchMask).any())
    MoFEMFunctionReturnHot(0);
  if (getFeatures().test(AUXILIARY_LOGARITHMIC_STRESS)) {
    CHKERR pushAuxiliaryLogarithmicStressOps(ep, pipeline, data_ptr, false);
    MoFEMFunctionReturnHot(0);
  }
 
  if (!ep.internalStressTagName.empty()) {
    switch (ep.meshTransferInterpOrder) {
    case 0:
      pipeline.push_back(
          new OpGetInternalStress<0>(data_ptr, ep.internalStressTagName));
      break;
    case 1:
      pipeline.push_back(
          new OpGetInternalStress<1>(data_ptr, ep.internalStressTagName));
      break;
    default:
      SETERRQ(ep.mField.get_comm(), MOFEM_NOT_IMPLEMENTED,
              "Unsupported mesh transfer interpolation order %d, for "
              "internal stress",
              ep.meshTransferInterpOrder);
    }
    TimeScale::ScalingFun default_scaling = [](double) { return 1; };
    auto time_scale =
        boost::make_shared<EshelbianCore::DynamicRelaxationTimeScale>(
            "internal_stress_history.txt", false, default_scaling);
    if (ep.internalStressVoigt)
      pipeline.push_back(new OpSpatialPhysicalInternalStress<true>(
          ep.stretchTensor, data_ptr, time_scale));
    else
      pipeline.push_back(new OpSpatialPhysicalInternalStress<false>(
          ep.stretchTensor, data_ptr, time_scale));
  }
  if (auto op = returnOpSpatialPhysicalExternalStrain(
          ep.stretchTensor, data_ptr, ep.externalStrainVecPtr, ep.timeScaleMap))
    pipeline.push_back(op);
  else if (ep.externalStrainVecPtr && !ep.externalStrainVecPtr->empty())
    SETERRQ(
        ep.mField.get_comm(), MOFEM_NOT_IMPLEMENTED,
        "OpSpatialPhysicalExternalStrain not implemented for this material");
 
  pipeline.push_back(
      returnOpSpatialPhysical(ep.stretchTensor, data_ptr, ep.alphaU));
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhysicalEquations::pushMaterialTangent(
    const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &pipeline,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr) {
  MoFEMFunctionBegin;
  if ((getFeatures() & noStretchMask).any()) {
    const bool heterogeneous = getFeatures().test(NON_HOMOGENEOUS_MATERIAL);
    pipeline.push_back(new OpSpatialConsistency_dP_dP(
        ep.piolaStress, ep.piolaStress, data_ptr, heterogeneous));
    pipeline.push_back(new OpSpatialConsistency_dBubble_dP(
        ep.bubbleField, ep.piolaStress, data_ptr, heterogeneous));
    pipeline.push_back(new OpSpatialConsistency_dBubble_dBubble(
        ep.bubbleField, ep.bubbleField, data_ptr, heterogeneous));
  } else if (getFeatures().test(AUXILIARY_LOGARITHMIC_STRESS)) {
    CHKERR pushAuxiliaryLogarithmicStressOps(ep, pipeline, data_ptr, true);
  } else {
    pipeline.push_back(returnOpSpatialPhysical_du_du(
        ep.stretchTensor, ep.stretchTensor, data_ptr, ep.alphaU));
    pipeline.push_back(new OpSpatialPhysical_du_dP(
        ep.stretchTensor, ep.piolaStress, data_ptr, true));
    pipeline.push_back(new OpSpatialPhysical_du_dBubble(
        ep.stretchTensor, ep.bubbleField, data_ptr, true));
    pipeline.push_back(
        new OpSpatialPhysical_du_domega(ep.stretchTensor, ep.rotAxis, data_ptr,
                                        ep.symmetrySelector == SYMMETRIC));
    if (ep.symmetrySelector > SYMMETRIC)
      pipeline.push_back(new OpSpatialRotation_domega_du(
          ep.rotAxis, ep.stretchTensor, data_ptr, false));
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhysicalEquations::pushEnergyEvaluation(
    const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &pipeline,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
    boost::shared_ptr<double> total_energy_ptr) {
  MoFEMFunctionBegin;
  if (auto op =
          returnOpCalculateHelmholtzFreeEnergy(data_ptr, total_energy_ptr)) {
    // Energy uses the physical strain recovered from stress for these models.
    if (usesStressBasedStrain())
      CHKERR pushStretchFromStress(ep, pipeline, data_ptr,
                                   data_ptr->getLogStretchTotalTensorAtPts());
    pipeline.push_back(op);
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhysicalEquations::pushMaterialForceFields(
    const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &pipeline,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr) {
  MoFEMFunctionBegin;
  if (usesStressBasedStrain()) {
    CHKERR pushStretchFromStress(ep, pipeline, data_ptr);
  } else {
    // These fields represent mechanical stretch; this force evaluation does
    // not reconstruct physical strain for imposed internal stress or strain.
    CHKERR pushMaterialFields(ep, pipeline, data_ptr);
  }
  MoFEMFunctionReturn(0);
}
 
UserDataOperator *PhysicalEquations::returnOpJacobian(
    const bool eval_rhs, const bool eval_lhs,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
    boost::shared_ptr<PhysicalEquations> physics_ptr) {
  return nullptr;
}
 
VolUserDataOperator *PhysicalEquations::returnOpSpatialPhysical(
    const std::string &field_name,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha_u) {
  return nullptr;
}
 
VolUserDataOperator *PhysicalEquations::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) {
  return nullptr;
}
 
VolUserDataOperator *PhysicalEquations::returnOpSpatialPhysical_du_du(
    std::string row_field, std::string col_field,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha) {
  return nullptr;
}
 
VolUserDataOperator *PhysicalEquations::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) {
  return nullptr;
}
 
VolUserDataOperator *PhysicalEquations::returnOpCalculateHelmholtzFreeEnergy(
    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
    boost::shared_ptr<double> total_helmholtz_free_energy_ptr) {
  return nullptr;
}
 
bool PhysicalEquations::providesHelmholtzFreeEnergy() const { return false; }
 
VolUserDataOperator *PhysicalEquations::returnOpCalculateStretchFromStress(
    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
    boost::shared_ptr<PhysicalEquations> physics_ptr,
    boost::shared_ptr<MatrixDouble> strain_ptr) {
  return nullptr;
}
 
VolUserDataOperator *PhysicalEquations::returnOpCalculateStretchFromStress(
    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
    boost::shared_ptr<PhysicalEquations> physics_ptr,
    boost::shared_ptr<MatrixDouble> strain_ptr,
    VectorPtr external_pressure_ptr) {
  return returnOpCalculateStretchFromStress(
      std::move(data_ptr), std::move(physics_ptr), std::move(strain_ptr));
}
 
VolUserDataOperator *PhysicalEquations::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) {
  return returnOpCalculateStretchFromStress(
      std::move(data_ptr), std::move(physics_ptr), std::move(strain_ptr),
      std::move(external_pressure_ptr));
}
 
VolUserDataOperator *PhysicalEquations::returnOpCalculateExternalPressure(
    VectorPtr external_pressure_ptr,
    boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
    std::map<std::string, boost::shared_ptr<ScalingMethod>> smv) {
  return nullptr;
}
 
VolUserDataOperator *PhysicalEquations::returnOpCalculateVarStretchFromStress(
    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
    boost::shared_ptr<PhysicalEquations> physics_ptr) {
  return nullptr;
};
 
MoFEMErrorCode PhysicalEquations::pushPostProc(
    const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &pipeline,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
    MaterialPostProcData &output) {
  MoFEMFunctionBegin;
  if ((getFeatures() & noStretchMask).any()) {
    pipeline.push_back(
        returnOpJacobian(true, false, data_ptr, ep.physicalEquations));
  } else {
    CHKERR pushMaterialEvaluation(ep, pipeline, data_ptr, true, false);
  }
  if (getFeatures().test(AUXILIARY_LOGARITHMIC_STRESS)) {
    CHKERR pushAuxiliaryLogarithmicPostProc(pipeline, data_ptr, output);
  } else if (auto op =
                 returnOpCalculateHelmholtzFreeEnergy(data_ptr, nullptr)) {
    pipeline.push_back(op);
  }
  output.symmetricFields["LogSpatialStretch"] =
      data_ptr->getLogStretchTensorAtPts();
  output.symmetricFields["SpatialStretch"] = data_ptr->getStretchTensorAtPts();
  if ((getFeatures() & noStretchMask).none())
    output.vectorFields["EiegnLogStreach"] = data_ptr->getEigenVals();
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhysicalEquations::pushMaterialVariation(
    const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &pipeline,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
    SmartPetscObj<Vec> variation, MaterialPostProcData *output) {
  MoFEMFunctionBegin;
  if (getFeatures().test(AUXILIARY_LOGARITHMIC_STRESS)) {
    auto values = boost::make_shared<AuxiliaryLogarithmicStressData>();
    CHKERR pushAuxiliaryLogarithmicFields(
        ep, pipeline, values, data_ptr->getVarLogStreachPts(), variation);
    if (output)
      CHKERR pushAuxiliaryLogarithmicStatePostProc(pipeline, values, "Var",
                                                   *output);
  } else if ((getFeatures() & noStretchMask).any()) {
    pipeline.push_back(
        returnOpCalculateVarStretchFromStress(data_ptr, ep.physicalEquations));
  } else {
    pipeline.push_back(new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
        ep.stretchTensor, data_ptr->getVarLogStreachPts(), variation, MBTET));
  }
  if (output)
    output->symmetricFields["VarLogSpatialStretch"] =
        data_ptr->getVarLogStreachPts();
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhysicalEquations::pushPostProcResidual(
    const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &pipeline,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
    SmartPetscObj<Vec> residual, MaterialPostProcData &output) {
  MoFEMFunctionBegin;
  if (getFeatures().test(AUXILIARY_LOGARITHMIC_STRESS)) {
    auto values = boost::make_shared<AuxiliaryLogarithmicStressData>();
    auto residual_stretch = boost::make_shared<MatrixDouble>();
    CHKERR pushAuxiliaryLogarithmicFields(ep, pipeline, values,
                                          residual_stretch, residual);
    CHKERR pushAuxiliaryLogarithmicStatePostProc(pipeline, values, "Res",
                                                 output);
  } else if (getFeatures().test(DIRECT_STRETCH)) {
    auto residual_stretch = boost::make_shared<MatrixDouble>();
    pipeline.push_back(new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
        ep.stretchTensor, residual_stretch, residual, MBTET));
    output.symmetricFields["ResLogSpatialStretch"] = residual_stretch;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhysicalEquations::pushSkeletonEvaluation(
    const EshelbianCore &ep, VolumeElementForcesAndSourcesCoreOnSide &fe,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
    SmartPetscObj<Vec> residual, MaterialPostProcData &output) {
  MoFEMFunctionBegin;
  auto &pipeline = fe.getOpPtrVector();
  fe.getUserPolynomialBase() = boost::make_shared<CGGUserPolynomialBase>(
      nullptr, !getFeatures().test(AUXILIARY_LOGARITHMIC_STRESS));
  CHKERR EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      pipeline, {HDIV, H1, L2}, ep.materialH1Positions, ep.frontAdjEdges);
  pipeline.push_back(new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
      ep.piolaStress, data_ptr->getApproxPAtPts()));
  pipeline.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
      ep.rotAxis, data_ptr->getRotAxisAtPts(), MBTET));
  pipeline.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
      ep.spatialL2Disp, data_ptr->getSmallWL2AtPts(), MBTET));
 
  if (getFeatures().test(AUXILIARY_LOGARITHMIC_STRESS)) {
    const bool add_bubble = ep.mField.check_field(ep.bubbleField);
    if (add_bubble)
      pipeline.push_back(new OpCalculateHTensorTensorField<3, 3>(
          ep.bubbleField, data_ptr->getApproxPAtPts(), MBMAXTYPE));
    CHKERR VecSetDM(ep.solTSStep.get(), PETSC_NULLPTR);
    pipeline.push_back(new OpCalculateHVecTensorField<3, 3>(
        ep.piolaStress, data_ptr->getApproxP0AtPts(), nullptr, ep.solTSStep));
    if (add_bubble)
      pipeline.push_back(new OpCalculateHTensorTensorField<3, 3>(
          ep.bubbleField, data_ptr->getApproxP0AtPts(), nullptr, ep.solTSStep,
          MBMAXTYPE));
    CHKERR pushMaterialFields(ep, pipeline, data_ptr, PREVIOUS);
    CHKERR pushMaterialFields(ep, pipeline, data_ptr, CURRENT);
    pipeline.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
        ep.rotAxis, data_ptr->getRotAxis0AtPts(), ep.solTSStep, MBTET));
    pipeline.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
        ep.spatialH1Disp, data_ptr->getSmallWH1AtPts()));
    pipeline.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        ep.spatialH1Disp, data_ptr->getSmallWGradH1AtPts()));
    CHKERR addCalculatePlasticLogarithmicStretchFieldValues(
        pipeline, ep.plasticHField, data_ptr->getPlasticH(), MBTET);
    pipeline.push_back(new OpCalculateRotationAndSpatialGradient(data_ptr));
    CHKERR pushMaterialEvaluation(ep, pipeline, data_ptr, true, false);
    CHKERR pushAuxiliaryLogarithmicPostProc(pipeline, data_ptr, output);
    output.symmetricFields["SpatialStretch"] =
        data_ptr->getStretchTensorAtPts();
    if (residual)
      CHKERR pushPostProcResidual(ep, pipeline, data_ptr, residual, output);
  } else {
    CHKERR pushMaterialFields(ep, pipeline, data_ptr, CURRENT);
  }
  output.symmetricFields["LogSpatialStretch"] =
      data_ptr->getLogStretchTensorAtPts();
  MoFEMFunctionReturn(0);
}
 
} // namespace EshelbianPlasticity
 
#include <impl/EshelbianMatCore.cpp>
#include <impl/HMHHencky.cpp>
#include <impl/HMHNeohookean.cpp>
namespace EshelbianPlasticity {
 
MoFEMErrorCode PhysicalEquations::create(EshelbianCore &ep) {
  MoFEMFunctionBegin;
  const char *list_materials[] = {"stvenant_kirchhoff", "mooney_rivlin",
                                  "hencky", "neo_hookean", "meta"};
  PetscInt choice_material = MooneyRivlin;
  PetscBool no_stretch = PETSC_FALSE;
  PetscBool auxiliary_logarithmic_stress = PETSC_FALSE;
  PetscOptionsBegin(ep.mField.get_comm(), "", "Material physical equations",
                    "none");
  CHKERR PetscOptionsEList("-material", "material model", "", list_materials,
                           LastMaterial, list_materials[choice_material],
                           &choice_material, PETSC_NULLPTR);
  CHKERR PetscOptionsBool("-no_stretch", "do not solve for stretch", "",
                          no_stretch, &no_stretch, PETSC_NULLPTR);
  CHKERR PetscOptionsBool(
      "-auxiliary_logarithmic_stress",
      "independent logarithmic deviator, volume and auxiliary stress", "",
      PETSC_FALSE, &auxiliary_logarithmic_stress, PETSC_NULLPTR);
  PetscOptionsEnd();
 
  const auto material_model = static_cast<MaterialModel>(choice_material);
  if (auxiliary_logarithmic_stress && no_stretch &&
      material_model == Neohookean)
    SETERRQ(
        ep.mField.get_comm(), MOFEM_NOT_IMPLEMENTED,
        "-auxiliary_logarithmic_stress cannot be combined with -no_stretch");
  Features features;
  features.set(auxiliary_logarithmic_stress
                   ? AUXILIARY_LOGARITHMIC_STRESS
                   : (no_stretch ? NO_STRETCH_LINEAR : DIRECT_STRETCH));
  CHKERR checkOptions(ep, material_model, features);
  if (!ep.dataAtPts)
    ep.dataAtPts = boost::make_shared<DataAtIntegrationPts>();
  const std::map<std::string, std::string> material_model_names_dictionary = {
      {"grad", "stretchH1AtPts"}, {"P", "PAtPts"}, {"P_dF", "PAtPts_du"}};
  MatOps::MatPiolaResponse::useDeformationGradient = true;
  boost::shared_ptr<MatOps::PhysicalEquations> mat_physical_equations_ptr;
 
  switch (material_model) {
  case StVenantKirchhoff:
    MOFEM_LOG("EP", Sev::inform) << "St. Venant-Kirchhoff material model";
    mat_physical_equations_ptr = MatOps::createMatOpsPhysicalEquationsPtr<
        MatOps::ELASTICITY::STVENANTKIRCHHOFF, MatOps::MODEL_3D>(
        MatOps::createMatOpsDataPtr(ep.dataAtPts->matOpsData,
                                    material_model_names_dictionary),
        MatOps::MatOpsTagsRegistry::setTagName("StVenantKirchhoff", 0));
    break;
  case MooneyRivlin:
    MOFEM_LOG("EP", Sev::inform) << "Mooney-Rivlin material model";
    mat_physical_equations_ptr = MatOps::createMatOpsPhysicalEquationsPtr<
        MatOps::ELASTICITY::MOONEYRIVLINWRIGGERSEQ63, MatOps::MODEL_3D>(
        MatOps::createMatOpsDataPtr(ep.dataAtPts->matOpsData,
                                    material_model_names_dictionary),
        MatOps::MatOpsTagsRegistry::setTagName("MooneyRivlinWriggersEq63", 0));
    break;
  case Hencky:
    MOFEM_LOG("EP", Sev::inform) << "Hencky material model";
    ep.physicalEquations =
        boost::make_shared<HMHHencky>(ep.mField, 5., 0.25, features);
    break;
  case Neohookean:
    MOFEM_LOG("EP", Sev::inform) << "Neo-Hookean material model";
    ep.physicalEquations =
        boost::make_shared<HMHNeohookean>(ep.mField, 1., 2., features);
    break;
  case Meta:
    MOFEM_LOG("EP", Sev::inform) << "Meta Material model";
    mat_physical_equations_ptr =
        MatOps::createMatOpsPhysicalEquationsPtr<MatOps::PIOLA_RESPONSE::META,
                                                 MatOps::MODEL_3D>(
            MatOps::createMatOpsDataPtr(ep.dataAtPts->matOpsData,
                                        material_model_names_dictionary),
            MatOps::MatOpsTagsRegistry::setTagName("MetaPiolaResponse", 0));
    break;
  default:
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown material");
  }
 
  if (mat_physical_equations_ptr) {
    CHKERR mat_physical_equations_ptr->getOptions(&ep.mField);
    CHKERR mat_physical_equations_ptr->recordTape();
    if (material_model == Meta)
      ep.physicalEquations = boost::make_shared<MatPhysicalEquations>(
          boost::static_pointer_cast<MatOps::MatPiolaResponse>(
              mat_physical_equations_ptr),
          material_model, features);
    else
      ep.physicalEquations = boost::make_shared<MatElasticPhysicalEquations>(
          boost::static_pointer_cast<MatOps::MatElastic>(
              mat_physical_equations_ptr),
          material_model, features);
  }
  ep.dataAtPts->physicsPtr = ep.physicalEquations;
  MOFEM_LOG("EP", Sev::inform) << "No stretch: -no_stretch "
                               << ((ep.physicalEquations->getFeatures() &
                                    PhysicalEquations::noStretchMask)
                                           .any()
                                       ? "yes"
                                       : "no");
  MoFEMFunctionReturn(0);
}
 
}; // namespace EshelbianPlasticity