HMHNeohookean.cpp

Go to the documentation of this file.

/**
 * @file HMHNeohookean.cpp
 * @brief Direct logarithmic Neo-Hookean material adapter.
 */
 
#include <AuxiliaryLogarithmicStressOperators.hpp>
#include <NeoHookeanLogarithmicMaterial.hpp>
#include <cmath>
#include <optional>
#include <regex>
#include <sstream>
 
namespace EshelbianPlasticity {
 
struct HMHNeohookean : public PhysicalEquations {
  using Material = NeoHookeanLogarithmicMaterial;
  using PackedStress = FTensor::Tensor1<double, size_symm>;
  using PackedTangent = FTensor::Tensor2<double, size_symm, size_symm>;
  using FormBase =
      typename FormsIntegrators<VolUserDataOperator>::Assembly<A>::OpBase;
 
  MoFEMErrorCode pushAuxiliaryLogarithmicMaterialEvaluation(
      const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &pipeline,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr, bool lhs) override {
    MoFEMFunctionBegin;
    auto parameters = prepareAuxiliaryEvaluation(data_ptr);
    CHKERR EshelbianPlasticity::pushAuxiliaryLogarithmicMaterialEvaluation(
        ep, pipeline, data_ptr, data_ptr->auxiliaryMaterialData,
        std::move(parameters), lhs);
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode pushAuxiliaryLogarithmicStressOps(
      const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &pipeline,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr, bool lhs) override {
    MoFEMFunctionBegin;
    CHKERR pushAuxiliaryLogarithmicMaterialOps(
        ep, pipeline, data_ptr->auxiliaryMaterialData, lhs);
    MoFEMFunctionReturn(0);
  }
 
  HMHNeohookean(MoFEM::Interface &m_field, const double c10,
                const double bulk_modulus, const Features features)
      : PhysicalEquations(Neohookean, features), mField(m_field),
        defaultParameters{c10, bulk_modulus} {
    if ((getFeatures() & noStretchMask).any())
      CHK_THROW_MESSAGE(
          MOFEM_NOT_IMPLEMENTED,
          "Neo-Hookean -no_stretch physical-stress recovery was "
          "retired; omit -no_stretch and solve the material fields");
    if (EshelbianCore::stretchSelector != LOG)
      CHK_THROW_MESSAGE(
          MOFEM_NOT_IMPLEMENTED,
          "Neo-Hookean requires natural logarithms: -stretches log");
    CHK_THROW_MESSAGE(getOptions(), "Neo-Hookean options are invalid");
    CHK_THROW_MESSAGE(extractBlockData(),
                      "Neo-Hookean material blocks are invalid");
  }
 
  Material::Parameters getParameters(const EntityHandle entity) const {
    for (const auto &block : blockData)
      if (block.entities.find(entity) != block.entities.end())
        return block.parameters;
    if (!blockData.empty())
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "MAT_NEOHOOKEAN blocks must cover every material element");
    return defaultParameters;
  }
 
  /** Map derivatives from five orthonormal coordinates and theta to the six
   * raw symmetric coefficients. SymmLTensor supplies both shear occurrences.
   */
  static MoFEMErrorCode
  evaluateDirect(const Material::Parameters &parameters,
                 const Material::SymmetricTensor &t_log_stretch, double &energy,
                 PackedStress *stress_ptr = nullptr,
                 PackedTangent *tangent_ptr = nullptr) {
    MoFEMFunctionBegin;
    FTENSOR_INDEXES(3, i, j);
    FTENSOR_INDEXES(size_symm, L, J);
    FTensor::Index<'a', 5> a;
    FTensor::Index<'b', 5> b;
    const auto t_coordinates = Material::getCoordinates(t_log_stretch);
    Material::Coordinates t_stress;
    Material::Tangent t_hessian;
    Material::VolumeState volume;
    CHKERR Material::evaluateDeviator(parameters, t_coordinates, energy,
                                      stress_ptr ? &t_stress : nullptr,
                                      tangent_ptr ? &t_hessian : nullptr);
    CHKERR Material::evaluateVolume(parameters, t_log_stretch(i, i), volume);
    energy += volume.energy;
    if (!std::isfinite(energy))
      SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
              "Neo-Hookean total strain energy is not representable");
    const auto t_packed_basis = FTensor::SymmLTensor<double, 3>();
    constexpr auto t_identity = FTensor::Kronecker_Delta_symmetric<int>();
    if (stress_ptr) {
      const auto t_deviatoric_stress = Material::getTensor(t_stress);
      Material::SymmetricTensor t_total_stress;
      t_total_stress(i, j) =
          t_deviatoric_stress(i, j) + volume.firstDerivative * t_identity(i, j);
      (*stress_ptr)(L) = t_packed_basis(i, j, L) * t_total_stress(i, j);
      for (int component = 0; component != size_symm; ++component)
        if (!std::isfinite((*stress_ptr)(component)))
          SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                  "Neo-Hookean packed material stress is not representable");
    }
    if (tangent_ptr) {
      const auto &t_basis = Tensor2SymmetricDeviatorBasis::getBasis();
      FTensor::Tensor2<double, 5, size_symm> t_projection;
      t_projection(a, L) = t_basis(i, j, a) * t_packed_basis(i, j, L);
      PackedStress t_trace;
      t_trace(L) = t_packed_basis(i, j, L) * t_identity(i, j);
      FTensor::Tensor2<double, 5, size_symm> t_hessian_projection;
      t_hessian_projection(a, J) = t_hessian(a, b) * t_projection(b, J);
      (*tangent_ptr)(L, J) = t_projection(a, L) * t_hessian_projection(a, J) +
                             volume.secondDerivative * t_trace(L) * t_trace(J);
      for (int row = 0; row != size_symm; ++row)
        for (int col = 0; col != size_symm; ++col)
          if (!std::isfinite((*tangent_ptr)(row, col)))
            SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                    "Neo-Hookean packed material tangent is not representable");
    }
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode evaluateElasticLogarithmicMaterialHessian(
      EntityHandle entity, int, DataAtIntegrationPts &,
      const Material::SymmetricTensor &elastic_log_stretch,
      Material::Tangent &hessian) override {
    MoFEMFunctionBegin;
    double energy = 0.;
    CHKERR Material::evaluateDeviator(
        getParameters(entity), Material::getCoordinates(elastic_log_stretch),
        energy, nullptr, &hessian);
    MoFEMFunctionReturn(0);
  }
 
  bool providesHelmholtzCurvature() const override { return true; }
 
  MoFEMErrorCode
  evaluateHelmholtzCurvature(EntityHandle entity, int point,
                             DataAtIntegrationPts &data,
                             LogarithmicMaterialCurvature &c) override {
    MoFEMFunctionBegin;
    using namespace LogarithmicDirectionalDerivatives;
    FTENSOR_INDEXES(3, i, j);
    FTensor::Index<'a', 5> a;
    FTensor::Index<'b', 5> b;
    FTensor::Index<'A', 6> A;
    const auto parameters = getParameters(entity);
    const bool auxiliary = getFeatures().test(AUXILIARY_LOGARITHMIC_STRESS);
    const double jac = std::exp(c.stretch(i, i));
    const double g2 = parameters.bulkModulus * jac * (2 * jac - 1);
    const double g3 = parameters.bulkModulus * jac * (4 * jac - 1);
    const auto basis = FTensor::SymmLTensor<double, 3>();
    c.stressResult(a) = 0.;
    const auto curvature =
        deviatorCurvature(parameters.c10, c.stretch, c.adjoint, c.direction);
    for (int n = 0; n != 6; ++n) {
      FTensor::Tensor1<double, 6> unit;
      unit(A) = 0.;
      unit(n) = 1.;
      Material::SymmetricTensor e;
      e(i, j) = basis(i, j, A) * unit(A);
      c.stretchResult(n) =
          curvature.second(i, j) * e(i, j) -
          (auxiliary ? 0. : curvature.third(i, j) * e(i, j)) +
          (g2 - c.adjoint(i, i) * g3) * c.direction(i, i) * e(i, i);
    }
    if (auxiliary) {
      using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
      const int count = data.auxiliaryMaterialData->compliance->size1();
      auto compliance =
          MatrixSizeHelper<GetFTensor2FromMatType<5, 5, -1, DL>, DL>::get(
              *data.auxiliaryMaterialData->compliance, count)();
      auto dm = MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::get(
          *data.auxiliaryMaterialData->materialDeviator, count)();
      for (int gg = 0; gg != point; ++gg) {
        ++compliance;
        ++dm;
      }
      Material::Coordinates p, s, e;
      p(a) = compliance(a, b) * c.stressAdjoint(b);
      s(a) = compliance(a, b) * c.stressDirection(b);
      const auto h = Tensor2SymmetricDeviatorBasis::getTensor(dm);
      const auto pt = Material::getTensor(p);
      const auto st = Material::getTensor(s);
      const auto stress_curvature =
          deviatorCurvature(parameters.c10, h, pt, st);
      for (int n = 0; n != 5; ++n) {
        Material::Coordinates unit;
        unit(a) = 0.;
        unit(n) = 1.;
        e(a) = compliance(a, b) * unit(b);
        c.stressResult(n) =
            -stress_curvature.third(i, j) * Material::getTensor(e)(i, j);
      }
    }
    MoFEMFunctionReturn(0);
  }
 
  struct OpJacobian : public EshelbianPlasticity::OpJacobian {
    using EshelbianPlasticity::OpJacobian::OpJacobian;
    MoFEMErrorCode evaluateRhs(EntData &) override { return 0; }
    MoFEMErrorCode evaluateLhs(EntData &) override { return 0; }
  };
 
  UserDataOperator *
  returnOpJacobian(const bool eval_rhs, const bool eval_lhs,
                   boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                   boost::shared_ptr<PhysicalEquations> physics_ptr) override {
    return new OpJacobian(eval_rhs, eval_lhs, std::move(data_ptr),
                          std::move(physics_ptr));
  }
 
  struct OpEnergy : public VolUserDataOperator {
    OpEnergy(boost::shared_ptr<DataAtIntegrationPts> data_ptr,
             boost::shared_ptr<double> total_energy_ptr)
        : VolUserDataOperator(NOSPACE, OPSPACE), dataAtPts(std::move(data_ptr)),
          totalEnergyPtr(std::move(total_energy_ptr)) {}
 
    MoFEMErrorCode doWork(int, EntityType, EntData &) override {
      MoFEMFunctionBegin;
      const auto material_ptr = getMaterial(dataAtPts);
      const auto parameters = material_ptr->getParameters(getFEEntityHandle());
      const int nb_integration_pts = getGaussPts().size2();
      auto t_eigenvalues = dataAtPts->getFTensorEigenVals(nb_integration_pts);
      dataAtPts->energyAtPts.resize(nb_integration_pts, false);
      auto t_energy = getFTensor0FromVec(dataAtPts->energyAtPts);
      auto t_weight = getFTensor0IntegrationWeight();
      auto t_plasticF = dataAtPts->getFTensorPlasticF(nb_integration_pts);
      double element_energy = 0.;
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        // Isotropic energy uses the principal logarithms supplied by either
        // configuration pipeline, independently of its displacement gradient.
        const Material::SymmetricTensor t_log_stretch(
            t_eigenvalues(0), 0., 0., t_eigenvalues(1), 0., t_eigenvalues(2));
        double energy;
        CHKERR evaluateDirect(parameters, t_log_stretch, energy);
        t_energy = energy;
        // Pull the intermediate-volume density back to the reference cell.
        const double det_plasticF = determinantTensor3by3(t_plasticF);
        element_energy += t_weight * energy * det_plasticF;
        ++t_eigenvalues;
        ++t_energy;
        ++t_weight;
        ++t_plasticF;
      }
      if (totalEnergyPtr)
        *totalEnergyPtr += getMeasure() * element_energy;
      MoFEMFunctionReturn(0);
    }
 
  private:
    boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
    boost::shared_ptr<double> totalEnergyPtr;
  };
 
  VolUserDataOperator *returnOpCalculateHelmholtzFreeEnergy(
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<double> total_energy_ptr) override {
    return new OpEnergy(std::move(data_ptr), std::move(total_energy_ptr));
  }
 
  bool providesHelmholtzFreeEnergy() const override { return true; }
 
  struct OpAuxiliaryHelmholtzGradient : public VolUserDataOperator {
    OpAuxiliaryHelmholtzGradient(
        boost::shared_ptr<DataAtIntegrationPts> data_ptr,
        boost::shared_ptr<MatrixDouble> deviator_gradient_ptr,
        boost::shared_ptr<MatrixDouble> volume_gradient_ptr)
        : VolUserDataOperator(NOSPACE, OPSPACE), dataAtPts(std::move(data_ptr)),
          deviatorGradientPtr(std::move(deviator_gradient_ptr)),
          volumeGradientPtr(std::move(volume_gradient_ptr)) {}
 
    MoFEMErrorCode doWork(int, EntityType, EntData &) override {
      MoFEMFunctionBegin;
      const auto material_ptr = getMaterial(dataAtPts);
      const auto parameters = material_ptr->getParameters(getFEEntityHandle());
      const int nb_gauss = getGaussPts().size2();
      const auto fields = dataAtPts->auxiliaryData;
      using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
      auto t_d = MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::get(
          *fields->logDeviator, nb_gauss)();
      auto t_theta =
          MatrixSizeHelper<GetFTensor1FromMatType<1, -1, DL>, DL>::get(
              *fields->logJacobian, nb_gauss)();
      auto t_deviator_gradient =
          MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::size(
              *deviatorGradientPtr, nb_gauss)();
      auto t_volume_gradient =
          MatrixSizeHelper<GetFTensor1FromMatType<1, -1, DL>, DL>::size(
              *volumeGradientPtr, nb_gauss)();
      FTENSOR_INDEX(5, L);
      for (int gg = 0; gg != nb_gauss; ++gg) {
        Material::Coordinates t_coordinates, t_stress;
        t_coordinates(L) = t_d(L);
        double energy;
        Material::VolumeState volume;
        CHKERR Material::evaluateDeviator(parameters, t_coordinates, energy,
                                          &t_stress);
        CHKERR Material::evaluateVolume(parameters, t_theta(0), volume);
        // Differentiate the physical energy at D_k. Weak copy closure does
        // not make the independent T_d its pointwise derivative.
        t_deviator_gradient(L) = t_stress(L);
        t_volume_gradient(0) = volume.firstDerivative;
        ++t_d;
        ++t_theta;
        ++t_deviator_gradient;
        ++t_volume_gradient;
      }
      MoFEMFunctionReturn(0);
    }
 
  private:
    boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
    boost::shared_ptr<MatrixDouble> deviatorGradientPtr;
    boost::shared_ptr<MatrixDouble> volumeGradientPtr;
  };
 
  MoFEMErrorCode pushHelmholtzStateGradient(
      const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &pipeline,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr) override {
    MoFEMFunctionBegin;
    if (!getFeatures().test(AUXILIARY_LOGARITHMIC_STRESS)) {
      CHKERR PhysicalEquations::pushHelmholtzStateGradient(ep, pipeline,
                                                           data_ptr);
      MoFEMFunctionReturnHot(0);
    }
    auto deviator_gradient = boost::make_shared<MatrixDouble>();
    auto volume_gradient = boost::make_shared<MatrixDouble>();
    pipeline.push_back(new OpAuxiliaryHelmholtzGradient(
        data_ptr, deviator_gradient, volume_gradient));
    using OpDeviatorGradient = FormsIntegrators<VolUserDataOperator>::Assembly<
        A>::LinearForm<GAUSS>::OpBaseTimesVector<1, 5, 5>;
    using OpVolumeGradient = FormsIntegrators<VolUserDataOperator>::Assembly<
        A>::LinearForm<GAUSS>::OpBaseTimesVector<1, 1, 1>;
    pipeline.push_back(
        new OpDeviatorGradient(ep.logDeviator, deviator_gradient));
    pipeline.push_back(new OpVolumeGradient(ep.logJacobian, volume_gradient));
    MoFEMFunctionReturn(0);
  }
 
  struct OpResidual : public FormBase {
    OpResidual(const std::string &field,
               boost::shared_ptr<DataAtIntegrationPts> data_ptr,
               const double alpha_u)
        : FormBase(field, field, FormBase::OPROW),
          dataAtPts(std::move(data_ptr)), alphaU(alpha_u) {
      CHK_THROW_MESSAGE(checkDirectConfiguration(),
                        "Unsupported direct Neo-Hookean configuration");
      if (!std::isfinite(alphaU))
        CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                          "Neo-Hookean stretch viscosity must be finite");
    }
 
    MoFEMErrorCode iNtegrate(EntData &row_data) override {
      MoFEMFunctionBegin;
      FTENSOR_INDEXES(3, i, j);
      FTENSOR_INDEX(size_symm, L);
      const auto t_packed_basis = FTensor::SymmLTensor<double, 3>();
      const auto material_ptr = getMaterial(dataAtPts);
      const auto parameters = material_ptr->getParameters(getFEEntityHandle());
      const double alpha_grad_u = material_ptr->alphaGradU;
      auto t_h = dataAtPts->getFTensorLogStretch(nbIntegrationPts);
      auto t_physical = dataAtPts->getFTensorAdjointPdU(nbIntegrationPts);
      // Inactive viscosity terms do not introduce rate-data dependencies.
      std::optional<decltype(dataAtPts->getFTensorLogStretchDot(
          nbIntegrationPts))>
          t_dot_h;
      std::optional<decltype(dataAtPts->getFTensorGradLogStretchDot(
          nbIntegrationPts))>
          t_grad_dot_h;
      if (alphaU != 0.)
        t_dot_h.emplace(dataAtPts->getFTensorLogStretchDot(nbIntegrationPts));
      if (alpha_grad_u != 0.)
        t_grad_dot_h.emplace(
            dataAtPts->getFTensorGradLogStretchDot(nbIntegrationPts));
      auto t_row = row_data.getFTensor0N();
      auto t_grad_row = row_data.getFTensor1DiffN<3>();
      auto t_weight = getFTensor0IntegrationWeight();
      for (int gg = 0; gg != nbIntegrationPts; ++gg) {
        Material::SymmetricTensor t_log_stretch;
        t_log_stretch(i, j) = t_h(i, j);
        PackedStress t_stress;
        double energy;
        CHKERR evaluateDirect(parameters, t_log_stretch, energy, &t_stress);
        const double alpha = getMeasure() * t_weight;
        PackedStress t_residual;
        t_residual(L) = t_stress(L) - t_physical(L);
        if (t_dot_h)
          t_residual(L) +=
              alphaU * (t_packed_basis(i, j, L) * (*t_dot_h)(i, j));
        FTensor::Tensor2<double, size_symm, 3> t_gradient;
        t_gradient(L, i) = 0.;
        if (t_grad_dot_h)
          t_gradient(L, i) = alpha_grad_u * (*t_grad_dot_h)(L, i);
        auto t_nf = getNf<size_symm>();
        int rr = 0;
        for (; rr != nbRows / size_symm; ++rr) {
          t_nf(L) += alpha *
                     (t_row * t_residual(L) + t_grad_row(i) * t_gradient(L, i));
          ++t_nf;
          ++t_row;
          ++t_grad_row;
        }
        for (; rr != nbRowBaseFunctions; ++rr) {
          ++t_row;
          ++t_grad_row;
        }
        ++t_weight;
        ++t_h;
        ++t_physical;
        if (t_dot_h)
          ++*t_dot_h;
        if (t_grad_dot_h)
          ++*t_grad_dot_h;
      }
      MoFEMFunctionReturn(0);
    }
 
  private:
    boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
    const double alphaU;
  };
 
  VolUserDataOperator *
  returnOpSpatialPhysical(const std::string &field,
                          boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                          const double alpha_u) override {
    return new OpResidual(field, std::move(data_ptr), alpha_u);
  }
 
  struct OpTangent : public FormBase {
    OpTangent(const std::string &row_field, const std::string &col_field,
              boost::shared_ptr<DataAtIntegrationPts> data_ptr,
              const double alpha_u)
        : FormBase(row_field, col_field, FormBase::OPROWCOL),
          dataAtPts(std::move(data_ptr)), alphaU(alpha_u) {
      sYmm = false;
      CHK_THROW_MESSAGE(checkDirectConfiguration(),
                        "Unsupported direct Neo-Hookean configuration");
      if (!std::isfinite(alphaU))
        CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                          "Neo-Hookean stretch viscosity must be finite");
    }
 
    MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data) override {
      MoFEMFunctionBegin;
      FTENSOR_INDEXES(3, i, j, k, l);
      FTENSOR_INDEXES(size_symm, L, J);
      const auto material_ptr = getMaterial(dataAtPts);
      const auto parameters = material_ptr->getParameters(getFEEntityHandle());
      const auto t_packed_basis = FTensor::SymmLTensor<double, 3>();
      constexpr auto t_identity = FTensor::Kronecker_Delta<int>();
      PackedTangent t_viscous_metric;
      t_viscous_metric(L, J) =
          t_packed_basis(i, j, L) * t_packed_basis(i, j, J);
      const double ts_a = getTSa();
      const double alpha_grad_u = material_ptr->alphaGradU;
      auto t_h = dataAtPts->getFTensorLogStretch(nbIntegrationPts);
      auto t_pullback = dataAtPts->getFTensorAdjointPdstretch(nbIntegrationPts);
      auto t_eigenvalues = dataAtPts->getFTensorEigenVals(nbIntegrationPts);
      auto t_eigenvectors = dataAtPts->getFTensorEigenVecs(nbIntegrationPts);
      auto t_row = row_data.getFTensor0N();
      auto t_grad_row = row_data.getFTensor1DiffN<3>();
      auto t_weight = getFTensor0IntegrationWeight();
#ifndef NDEBUG
      if (dataAtPts->nbUniq.size() != nbIntegrationPts)
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Neo-Hookean tangent requires current geometric eigendata");
#endif
      for (int gg = 0; gg != nbIntegrationPts; ++gg) {
        Material::SymmetricTensor t_log_stretch;
        t_log_stretch(i, j) = t_h(i, j);
        PackedTangent t_tangent;
        double energy;
        CHKERR evaluateDirect(parameters, t_log_stretch, energy, nullptr,
                              &t_tangent);
        Material::SymmetricTensor t_symmetric_pullback;
        t_symmetric_pullback(i, j) =
            (t_pullback(i, j) || t_pullback(j, i)) / 2.;
        const auto t_work_curvature = EigenMatrix::getDiffDiffMat(
            t_eigenvalues, t_eigenvectors, EshelbianCore::f, EshelbianCore::d_f,
            EshelbianCore::dd_f, t_symmetric_pullback, dataAtPts->nbUniq[gg]);
        t_tangent(L, J) -=
            t_packed_basis(i, j, L) *
            (t_work_curvature(i, j, k, l) * t_packed_basis(k, l, J));
        t_tangent(L, J) += (alphaU * ts_a) * t_viscous_metric(L, J);
        const double alpha = getMeasure() * t_weight;
        int rr = 0;
        for (; rr != nbRows / size_symm; ++rr) {
          auto t_col = col_data.getFTensor0N(gg, 0);
          auto t_grad_col = col_data.getFTensor1DiffN<3>(gg, 0);
          auto t_matrix = getLocMat<size_symm>(size_symm * rr);
          for (int cc = 0; cc != nbCols / size_symm; ++cc) {
            t_matrix(L, J) += alpha * t_row * t_col * t_tangent(L, J);
            const double gradient_weight =
                alpha * alpha_grad_u * ts_a * (t_grad_row(i) * t_grad_col(i));
            t_matrix(L, J) += gradient_weight * t_identity(L, J);
            ++t_matrix;
            ++t_col;
            ++t_grad_col;
          }
          ++t_row;
          ++t_grad_row;
        }
        for (; rr != nbRowBaseFunctions; ++rr) {
          ++t_row;
          ++t_grad_row;
        }
        ++t_weight;
        ++t_h;
        ++t_pullback;
        ++t_eigenvalues;
        ++t_eigenvectors;
      }
      MoFEMFunctionReturn(0);
    }
 
  private:
    boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
    const double alphaU;
  };
 
  VolUserDataOperator *returnOpSpatialPhysical_du_du(
      std::string row_field, std::string col_field,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      const double alpha_u) override {
    return new OpTangent(row_field, col_field, std::move(data_ptr), alpha_u);
  }
 
  struct OpExternalStrain : public FormBase {
    OpExternalStrain(
        const std::string &field,
        boost::shared_ptr<ExternalStrainVec> external_strain_ptr,
        std::map<std::string, boost::shared_ptr<ScalingMethod>> scaling_methods)
        : FormBase(field, field, FormBase::OPROW),
          externalStrainPtr(std::move(external_strain_ptr)),
          scalingMethods(std::move(scaling_methods)) {}
 
    MoFEMErrorCode iNtegrate(EntData &row_data) override {
      MoFEMFunctionBegin;
      if (!externalStrainPtr)
        MoFEMFunctionReturnHot(0);
      const double time = EshelbianCore::physicalTimeFlg
                              ? EshelbianCore::currentPhysicalTime
                              : getFEMethod()->ts_t;
      double pressure = 0.;
      for (const auto &block : *externalStrainPtr) {
        if (block.blockName.find("ANALYTICAL_EXTERNALSTRAIN") !=
            std::string::npos)
          SETERRQ(
              PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
              "Analytical external strain is not implemented for Neo-Hookean");
        if (block.ents.find(getFEEntityHandle()) == block.ents.end())
          continue;
        const auto scaling = scalingMethods.find(block.blockName);
        double scale = 1.;
        if (scaling != scalingMethods.end())
          scale = scaling->second->getScale(time);
        else
          MOFEM_LOG("SELF", Sev::warning)
              << "No scaling method found for " << block.blockName;
        pressure += 3. * block.bulkModulusK * block.val * scale;
      }
      if (!std::isfinite(pressure))
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Neo-Hookean prescribed external pressure must be finite");
      FTENSOR_INDEXES(3, i, j);
      FTENSOR_INDEX(size_symm, L);
      const auto t_packed_basis = FTensor::SymmLTensor<double, 3>();
      constexpr auto t_identity = FTensor::Kronecker_Delta_symmetric<int>();
      PackedStress t_load;
      t_load(L) = pressure * (t_packed_basis(i, j, L) * t_identity(i, j));
      auto t_row = row_data.getFTensor0N();
      auto t_weight = getFTensor0IntegrationWeight();
      for (int gg = 0; gg != nbIntegrationPts; ++gg) {
        const double alpha = getMeasure() * t_weight;
        auto t_nf = getNf<size_symm>();
        int rr = 0;
        for (; rr != nbRows / size_symm; ++rr) {
          t_nf(L) -= alpha * t_row * t_load(L);
          ++t_nf;
          ++t_row;
        }
        for (; rr != nbRowBaseFunctions; ++rr)
          ++t_row;
        ++t_weight;
      }
      MoFEMFunctionReturn(0);
    }
 
  private:
    boost::shared_ptr<ExternalStrainVec> externalStrainPtr;
    std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethods;
  };
 
  VolUserDataOperator *returnOpSpatialPhysicalExternalStrain(
      const std::string &field, boost::shared_ptr<DataAtIntegrationPts>,
      boost::shared_ptr<ExternalStrainVec> external_strain_ptr,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> scaling_methods)
      override {
    return new OpExternalStrain(field, std::move(external_strain_ptr),
                                std::move(scaling_methods));
  }
 
  VolUserDataOperator *returnOpCalculateExternalPressure(
      VectorPtr pressure_ptr,
      boost::shared_ptr<ExternalStrainVec> external_strain_ptr,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> scaling_methods)
      override {
    return new OpCalculateExternalPressure(std::move(pressure_ptr),
                                           std::move(external_strain_ptr),
                                           std::move(scaling_methods));
  }
 
  VolUserDataOperator *
  returnOpCalculateStretchFromStress(boost::shared_ptr<DataAtIntegrationPts>,
                                     boost::shared_ptr<PhysicalEquations>,
                                     boost::shared_ptr<MatrixDouble>) override {
    return retiredInverse();
  }
  VolUserDataOperator *
  returnOpCalculateStretchFromStress(boost::shared_ptr<DataAtIntegrationPts>,
                                     boost::shared_ptr<PhysicalEquations>,
                                     boost::shared_ptr<MatrixDouble>,
                                     VectorPtr) override {
    return retiredInverse();
  }
  VolUserDataOperator *returnOpCalculateStretchFromStress(
      boost::shared_ptr<DataAtIntegrationPts>,
      boost::shared_ptr<PhysicalEquations>, boost::shared_ptr<MatrixDouble>,
      boost::shared_ptr<MatrixDouble>, VectorPtr) override {
    return retiredInverse();
  }
  VolUserDataOperator *returnOpCalculateVarStretchFromStress(
      boost::shared_ptr<DataAtIntegrationPts>,
      boost::shared_ptr<PhysicalEquations>) override {
    return retiredInverse();
  }
 
private:
  AuxiliaryLogarithmicMaterialParameters prepareAuxiliaryEvaluation(
      const boost::shared_ptr<DataAtIntegrationPts> &data_ptr) {
    if (!data_ptr)
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Auxiliary material evaluation requires field data");
    if (!data_ptr->auxiliaryMaterialData)
      data_ptr->auxiliaryMaterialData =
          boost::make_shared<AuxiliaryLogarithmicStressMaterialData>();
    return [this](EntityHandle entity) { return getParameters(entity); };
  }
 
  static MoFEMErrorCode checkDirectConfiguration() {
    MoFEMFunctionBegin;
    if (EshelbianCore::gradApproximator != NO_H1_CONFIGURATION)
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
              "Neo-Hookean direct material equations require -grad no_h1");
    if (EshelbianCore::rotSelector != LARGE_ROT &&
        EshelbianCore::rotSelector != SMALL_ROT)
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
              "Neo-Hookean direct material equations require "
              "-rotations large or -rotations small");
    MoFEMFunctionReturn(0);
  }
 
  static boost::shared_ptr<HMHNeohookean>
  getMaterial(const boost::shared_ptr<DataAtIntegrationPts> &data_ptr) {
    if (!data_ptr)
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Neo-Hookean integration-point data is missing");
    const auto material_ptr =
        boost::dynamic_pointer_cast<HMHNeohookean>(data_ptr->physicsPtr);
    if (!material_ptr)
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Neo-Hookean material pointer is missing");
    return material_ptr;
  }
 
  static VolUserDataOperator *retiredInverse() {
    CHK_THROW_MESSAGE(MOFEM_NOT_IMPLEMENTED,
                      "Neo-Hookean physical-stress-to-stretch recovery was "
                      "retired; solve the logarithmic material fields");
    return nullptr;
  }
 
  MoFEMErrorCode getOptions() {
    MoFEMFunctionBegin;
    PetscOptionsBegin(mField.get_comm(), "neo_hookean_", "Neo-Hookean material",
                      "none");
    CHKERR PetscOptionsScalar("-c10", "C10", "", defaultParameters.c10,
                              &defaultParameters.c10, nullptr);
    CHKERR PetscOptionsScalar("-K", "Bulk modulus", "",
                              defaultParameters.bulkModulus,
                              &defaultParameters.bulkModulus, nullptr);
    CHKERR PetscOptionsScalar("-viscosity_alpha_grad_u",
                              "Logarithmic-stretch-gradient rate viscosity", "",
                              alphaGradU, &alphaGradU, nullptr);
    PetscOptionsEnd();
    CHKERR Material::validateParameters(defaultParameters);
    if (!std::isfinite(alphaGradU))
      SETERRQ(
          mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
          "Neo-Hookean logarithmic-stretch-gradient viscosity must be finite");
    char *options_ptr = nullptr;
    CHKERR PetscOptionsGetAll(nullptr, &options_ptr);
    std::istringstream options(options_ptr ? options_ptr : "");
    CHKERR PetscFree(options_ptr);
    std::string option;
    while (options >> option)
      if (option.rfind("-nh_stretch_", 0) == 0 ||
          option == "-neo_hookean_min_eigen_value")
        SETERRQ(mField.get_comm(), MOFEM_NOT_IMPLEMENTED,
                "Option %s belongs to the retired Neo-Hookean physical inverse "
                "or tangent projection; remove it",
                option.c_str());
    MOFEM_LOG("EP", Sev::inform) << "Neo-Hookean C10=" << defaultParameters.c10
                                 << " K=" << defaultParameters.bulkModulus
                                 << " alphaGradU=" << alphaGradU;
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode extractBlockData() {
    MoFEMFunctionBegin;
    auto meshsets = mField.getInterface<MeshsetsManager>();
    auto json_config = mField.getInterface<JsonConfigManager>();
    for (const auto block :
         meshsets->getCubitMeshsetPtr(std::regex("MAT_NEOHOOKEAN(.*)"))) {
      Material::Parameters parameters;
      const auto json_parameters = json_config->getParamsFromBlockset(
          "MAT_NEOHOOKEAN", block->getMeshsetId());
      if (!json_parameters.empty()) {
        if (json_parameters.size() != 2 || !json_parameters.count("c10") ||
            !json_parameters.count("k"))
          SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
                  "MAT_NEOHOOKEAN JSON block must contain exactly c10,k");
        parameters = {json_parameters.at("c10"), json_parameters.at("k")};
      } else {
        std::vector<double> attributes;
        CHKERR block->getAttributes(attributes);
        if (attributes.size() < 2)
          SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
                  "MAT_NEOHOOKEAN block requires C10,K attributes");
        parameters = {attributes[0], attributes[1]};
      }
      CHKERR Material::validateParameters(parameters);
      Range entities;
      CHKERR mField.get_moab().get_entities_by_handle(block->getMeshset(),
                                                      entities, true);
      blockData.push_back({parameters, std::move(entities)});
      MOFEM_LOG("EP", Sev::inform)
          << "MAT_NEOHOOKEAN " << block->getMeshsetId()
          << " C10=" << parameters.c10 << " K=" << parameters.bulkModulus;
    }
    MoFEMFunctionReturn(0);
  }
 
  struct BlockData {
    Material::Parameters parameters;
    Range entities;
  };
  MoFEM::Interface &mField;
  Material::Parameters defaultParameters;
  double alphaGradU = 0.;
  std::vector<BlockData> blockData;
};
 
} // namespace EshelbianPlasticity