AuxiliaryLogarithmicStress.cpp

Go to the documentation of this file.

/**
 * @file AuxiliaryLogarithmicStress.cpp
 * @brief Auxiliary logarithmic fields, material equations and postprocessing.
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <AuxiliaryLogarithmicStress.hpp>
#include <AuxiliaryLogarithmicStressOperators.hpp>
#include <AuxiliaryLogarithmicStressPostProc.hpp>
#include <EshelbianPlasticity.hpp>
#include <MatrixFunction.hpp>
 
#include <cmath>
#include <optional>
 
namespace EshelbianPlasticity {
namespace {
 
using Material = NeoHookeanLogarithmicMaterial;
using MaterialData = AuxiliaryLogarithmicStressMaterialData;
using OpBase = FormsIntegrators<VolUserDataOperator>::Assembly<A>::OpBase;
 
/** Nonassembling producer; the field evaluators own its input storage. */
struct OpReconstructAuxiliaryLogarithmicStretch : public VolUserDataOperator {
  OpReconstructAuxiliaryLogarithmicStretch(
      boost::shared_ptr<AuxiliaryLogarithmicStressData> values,
      MatrixPtr reconstructed_log_stretch)
      : VolUserDataOperator(NOSPACE, OPSPACE), valuesPtr(std::move(values)),
        logStretchPtr(std::move(reconstructed_log_stretch)) {}
 
  MoFEMErrorCode doWork(int, EntityType, EntData &) override {
    MoFEMFunctionBegin;
    const int nb_gauss = getGaussPts().size2();
    auto get_d =
        MatrixSizeHelper<GetFTensor2SymmetricFromMatType<3, -1, DL>, DL>::get(
            *valuesPtr->logDeviatorTensor, nb_gauss);
    auto get_theta =
        MatrixSizeHelper<GetFTensor1FromMatType<1, -1, DL>, DL>::get(
            *valuesPtr->logJacobian, nb_gauss);
    auto get_h =
        MatrixSizeHelper<GetFTensor2SymmetricFromMatType<3, -1, DL>, DL>::size(
            *logStretchPtr, nb_gauss);
    auto t_d = get_d();
    auto t_theta = get_theta();
    auto t_h = get_h();
    FTENSOR_INDEXES(3, i, j);
    constexpr auto t_identity = FTensor::Kronecker_Delta_symmetric<int>();
    for (int gg = 0; gg != nb_gauss; ++gg) {
      t_h(i, j) = t_d(i, j) + (t_theta(0) / 3.) * t_identity(i, j);
      ++t_d;
      ++t_theta;
      ++t_h;
    }
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<AuxiliaryLogarithmicStressData> valuesPtr;
  MatrixPtr logStretchPtr;
};
 
/** Raw symmetric coordinates are (00,01,02,11,12,22), without shear scaling.
 * Each column is a symmetric tensor's six distinct entries. Contracting the
 * full SymmLTensor with B here would incorrectly double the shear entries.
 */
const FTensor::Tensor2<double, 6, 5> &getPackedDeviatorBasis() {
  static const auto t_projection = [] {
    FTensor::Index<'a', 5> a;
    const auto &t_basis = Tensor2SymmetricDeviatorBasis::getBasis();
    const FTensor::Number<0> n0;
    const FTensor::Number<1> n1;
    const FTensor::Number<2> n2;
    const FTensor::Number<3> n3;
    const FTensor::Number<4> n4;
    const FTensor::Number<5> n5;
    FTensor::Tensor2<double, 6, 5> t_projection;
    t_projection(n0, a) = t_basis(n0, n0, a);
    t_projection(n1, a) = t_basis(n0, n1, a);
    t_projection(n2, a) = t_basis(n0, n2, a);
    t_projection(n3, a) = t_basis(n1, n1, a);
    t_projection(n4, a) = t_basis(n1, n2, a);
    t_projection(n5, a) = t_basis(n2, n2, a);
    return t_projection;
  }();
  return t_projection;
}
 
const FTensor::Tensor1<double, 6> &getPackedVolumeBasis() {
  static const FTensor::Tensor1<double, 6> t_projection(1. / 3., 0., 0.,
                                                        1. / 3., 0., 1. / 3.);
  return t_projection;
}
 
MoFEMErrorCode validateValues(const MatrixDouble &values) {
  MoFEMFunctionBegin;
  if (!std::all_of(values.data().begin(), values.data().end(),
                   [](const double value) { return std::isfinite(value); }))
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "Auxiliary logarithmic material or stress-work values are not "
            "representable");
  MoFEMFunctionReturn(0);
}
 
/** Nonassembling producer. Geometry and total P (including bubbles) precede it.
 */
struct OpCalculateAuxiliaryLogarithmicMaterial : public VolUserDataOperator {
  OpCalculateAuxiliaryLogarithmicMaterial(
      boost::shared_ptr<DataAtIntegrationPts> data,
      boost::shared_ptr<MaterialData> material_data,
      AuxiliaryLogarithmicMaterialParameters parameters, const double alpha_u,
      const bool lhs)
      : VolUserDataOperator(NOSPACE, OPSPACE), dataAtPts(std::move(data)),
        materialData(std::move(material_data)),
        parameterFun(std::move(parameters)), alphaU(alpha_u),
        calculateLhs(lhs) {
    if (!std::isfinite(alphaU))
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Neo-Hookean deviator viscosity must be finite");
  }
 
  MoFEMErrorCode doWork(int, EntityType, EntData &) override {
    MoFEMFunctionBegin;
    const int nb_gauss = getGaussPts().size2();
    const auto parameters = parameterFun(getFEEntityHandle());
    const MoFEMErrorCode material_error =
        calculateMaterial(parameters, nb_gauss);
    if (material_error == MOFEM_OPERATION_UNSUCCESSFUL && !calculateLhs &&
        getSNESCtx() == SnesMethod::CTX_SNESSETFUNCTION &&
        getFEMethod()->snes) {
      // Let SNES reject an unrepresentable trial and release its work vectors.
      // Returning a raw callback error leaves the solver unusable on retry.
      CHKERR SNESSetFunctionDomainError(getFEMethod()->snes);
      materialData->residualDeviator->clear();
      materialData->residualJacobian->clear();
      materialData->residualStress->clear();
      MoFEMFunctionReturnHot(0);
    }
    CHKERR material_error;
    if (calculateLhs)
      CHKERR calculateTangent(nb_gauss);
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
  boost::shared_ptr<MaterialData> materialData;
  AuxiliaryLogarithmicMaterialParameters parameterFun;
  const double alphaU;
  bool calculateLhs;
 
  MoFEMErrorCode calculateMaterial(const Material::Parameters &parameters,
                                   const int nb_gauss) {
    MoFEMFunctionBegin;
    const auto fields = dataAtPts->auxiliaryData;
    auto get_d = MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::get(
        *fields->logDeviator, nb_gauss);
    auto get_theta =
        MatrixSizeHelper<GetFTensor1FromMatType<1, -1, DL>, DL>::get(
            *fields->logJacobian, nb_gauss);
    auto get_stress =
        MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::get(
            *fields->stress, nb_gauss);
    auto get_elastic_stress =
        MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::size(
            *materialData->elasticStress, nb_gauss);
    auto get_dm = MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::size(
        *materialData->materialDeviator, nb_gauss);
    auto get_compliance =
        MatrixSizeHelper<GetFTensor2FromMatType<5, 5, -1, DL>, DL>::size(
            *materialData->compliance, nb_gauss);
    auto get_energy =
        MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::size(
            *materialData->energy, nb_gauss);
    auto get_volume =
        MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::size(
            *materialData->volume, nb_gauss);
    auto get_rd = MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::size(
        *materialData->residualDeviator, nb_gauss);
    auto get_rtheta =
        MatrixSizeHelper<GetFTensor1FromMatType<1, -1, DL>, DL>::size(
            *materialData->residualJacobian, nb_gauss);
    auto get_rt = MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::size(
        *materialData->residualStress, nb_gauss);
    auto t_d = get_d();
    auto t_theta = get_theta();
    auto t_stress = get_stress();
    auto t_elastic_stress = get_elastic_stress();
    std::optional<decltype(t_d)> t_dot_d;
    if (alphaU != 0.)
      t_dot_d.emplace(
          MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::get(
              *fields->logDeviatorDot, nb_gauss)());
    auto t_dm = get_dm();
    auto t_compliance = get_compliance();
    auto t_energy = get_energy();
    auto t_volume = get_volume();
    auto t_rd = get_rd();
    auto t_rtheta = get_rtheta();
    auto t_rt = get_rt();
    auto t_work = dataAtPts->getFTensorAdjointPdU(nb_gauss);
    auto t_plastic_f = dataAtPts->getFTensorPlasticF(nb_gauss);
    const auto &t_projection = getPackedDeviatorBasis();
    const auto &t_volume_projection = getPackedVolumeBasis();
    FTensor::Index<'a', 5> a;
    FTensor::Index<'b', 5> b;
    FTensor::Index<'L', 6> L;
 
    for (int gg = 0; gg != nb_gauss; ++gg) {
      Material::Coordinates t_current_d, t_current_stress;
      t_current_d(a) = t_d(a);
      t_current_stress(a) = t_stress(a);
      if (t_dot_d)
        t_current_stress(a) -= alphaU * (*t_dot_d)(a);
      t_elastic_stress(a) = t_current_stress(a);
      Material::InverseState inverse;
      Material::VolumeState volume;
      CHKERR Material::evaluateInverse(parameters, t_current_stress, inverse);
      CHKERR Material::evaluateVolume(parameters, t_theta(0), volume);
      double deviatoric_energy;
      CHKERR Material::evaluateDeviator(parameters, t_current_d,
                                        deviatoric_energy);
      t_dm(a) = inverse.tMaterialDeviator(a);
      t_compliance(a, b) = inverse.tCompliance(a, b);
      t_energy(MaterialData::DEVIATORIC) = deviatoric_energy;
      t_energy(MaterialData::VOLUMETRIC) = volume.energy;
      t_energy(MaterialData::CONJUGATE) = inverse.conjugateEnergy;
      t_volume(MaterialData::JACOBIAN) = volume.jacobian;
      t_volume(MaterialData::FIRST_DERIVATIVE) = volume.firstDerivative;
      t_volume(MaterialData::SECOND_DERIVATIVE) = volume.secondDerivative;
 
      // P_bar=P Fp^T/Jp. Total Td balances stress work; the elastic part
      // Td-alphaU*dotD determines Dm through the pointwise material inverse.
      const double det_plastic_f = determinantTensor3by3(t_plastic_f);
      t_rd(a) = det_plastic_f * (t_stress(a) - t_projection(L, a) * t_work(L));
      t_rtheta(0) = det_plastic_f * (volume.firstDerivative -
                                     t_volume_projection(L) * t_work(L));
      t_rt(a) = det_plastic_f * (t_d(a) - inverse.tMaterialDeviator(a));
      ++t_d;
      ++t_theta;
      ++t_stress;
      ++t_elastic_stress;
      if (t_dot_d)
        ++*t_dot_d;
      ++t_dm;
      ++t_compliance;
      ++t_energy;
      ++t_volume;
      ++t_rd;
      ++t_rtheta;
      ++t_rt;
      ++t_work;
      ++t_plastic_f;
    }
    CHKERR validateValues(*materialData->residualDeviator);
    CHKERR validateValues(*materialData->residualJacobian);
    CHKERR validateValues(*materialData->residualStress);
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode calculateTangent(const int nb_gauss) {
    MoFEMFunctionBegin;
#ifndef NDEBUG
    if (getTSCtx() != TSMethod::CTX_TSSETIJACOBIAN)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Auxiliary tangent evaluation requires the current geometry "
              "Jacobian context");
    if (dataAtPts->nbUniq.size() != nb_gauss)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Auxiliary work curvature requires the current stretch spectrum");
#endif
    auto get_compliance =
        MatrixSizeHelper<GetFTensor2FromMatType<5, 5, -1, DL>, DL>::get(
            *materialData->compliance, nb_gauss);
    auto get_volume =
        MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::get(
            *materialData->volume, nb_gauss);
    auto get_dd =
        MatrixSizeHelper<GetFTensor2FromMatType<5, 5, -1, DL>, DL>::size(
            *materialData->tangentDeviator, nb_gauss);
    auto get_dtheta =
        MatrixSizeHelper<GetFTensor2FromMatType<5, 1, -1, DL>, DL>::size(
            *materialData->tangentDeviatorJacobian, nb_gauss);
    auto get_dt =
        MatrixSizeHelper<GetFTensor2FromMatType<5, 5, -1, DL>, DL>::size(
            *materialData->tangentDeviatorStress, nb_gauss);
    auto get_td =
        MatrixSizeHelper<GetFTensor2FromMatType<5, 5, -1, DL>, DL>::size(
            *materialData->tangentStressDeviator, nb_gauss);
    auto get_thetatheta =
        MatrixSizeHelper<GetFTensor2FromMatType<1, 1, -1, DL>, DL>::size(
            *materialData->tangentJacobian, nb_gauss);
    auto get_tt =
        MatrixSizeHelper<GetFTensor2FromMatType<5, 5, -1, DL>, DL>::size(
            *materialData->tangentStress, nb_gauss);
    auto get_domega =
        MatrixSizeHelper<GetFTensor2FromMatType<5, 3, -1, DL>, DL>::size(
            *materialData->deviatorRotation, nb_gauss);
    auto get_thetaomega =
        MatrixSizeHelper<GetFTensor2FromMatType<1, 3, -1, DL>, DL>::size(
            *materialData->jacobianRotation, nb_gauss);
    auto get_fd =
        MatrixSizeHelper<GetFTensor3FromMatType<3, 3, 5, -1, DL>, DL>::size(
            *materialData->deformationDeviator, nb_gauss);
    auto get_ftheta =
        MatrixSizeHelper<GetFTensor3FromMatType<3, 3, 1, -1, DL>, DL>::size(
            *materialData->deformationJacobian, nb_gauss);
    auto t_compliance = get_compliance();
    auto t_volume = get_volume();
    auto t_dd = get_dd();
    auto t_dtheta = get_dtheta();
    auto t_dt = get_dt();
    auto t_td = get_td();
    auto t_thetatheta = get_thetatheta();
    auto t_tt = get_tt();
    auto t_domega = get_domega();
    auto t_thetaomega = get_thetaomega();
    auto t_fd = get_fd();
    auto t_ftheta = get_ftheta();
    auto t_work = dataAtPts->getFTensorAdjointPdU(nb_gauss);
    auto t_work_omega = dataAtPts->getFTensorAdjointPdUdOmega(nb_gauss);
    auto t_fh = dataAtPts->getFTensorAdjointPdUdP(nb_gauss);
    auto t_adjoint_p = dataAtPts->getFTensorAdjointPdstretch(nb_gauss);
    auto t_eigen_values = dataAtPts->getFTensorEigenVals(nb_gauss);
    auto t_eigen_vectors = dataAtPts->getFTensorEigenVecs(nb_gauss);
    auto t_plastic_f = dataAtPts->getFTensorPlasticF(nb_gauss);
    const auto &t_basis = Tensor2SymmetricDeviatorBasis::getBasis();
    const auto &t_projection = getPackedDeviatorBasis();
    const auto &t_volume_projection = getPackedVolumeBasis();
    const EigenMatrix::Fun<double> exp_fun = [](const double value) {
      return std::exp(value);
    };
    FTENSOR_INDEXES(3, i, j, k, l);
    FTensor::Index<'a', 5> a;
    FTensor::Index<'b', 5> b;
    FTensor::Index<'L', 6> L;
    FTensor::Index<'z', 1> z;
    const FTensor::Tensor1<double, 1> t_scalar_component(1.);
    const FTensor::Number<0> n0;
    constexpr auto t_identity = FTensor::Kronecker_Delta<double>();
    const double alpha_u_shift = alphaU * getTSa();
 
    for (int gg = 0; gg != nb_gauss; ++gg) {
      FTensor::Tensor2_symmetric<double, 3> t_symmetric_adjoint;
      t_symmetric_adjoint(i, j) = (t_adjoint_p(i, j) || t_adjoint_p(j, i)) / 2.;
      const auto t_curvature = EigenMatrix::getDiffDiffMat(
          t_eigen_values, t_eigen_vectors, exp_fun, exp_fun, exp_fun,
          t_symmetric_adjoint, dataAtPts->nbUniq[gg]);
 
      const double det_plastic_f = determinantTensor3by3(t_plastic_f);
      t_dd(a, b) =
          -det_plastic_f *
          (t_basis(i, j, a) * (t_curvature(i, j, k, l) * t_basis(k, l, b)));
      // exp(D+theta*I/3) separates exactly, including noncommuting D modes.
      t_dtheta(a, n0) = -det_plastic_f * (t_projection(L, a) * t_work(L)) / 3.;
      t_thetatheta(0, 0) =
          det_plastic_f * (t_volume(MaterialData::SECOND_DERIVATIVE) -
                           (t_volume_projection(L) * t_work(L)) / 3.);
      t_dt(a, b) = det_plastic_f * t_identity(a, b);
      t_td(a, b) = det_plastic_f *
                   (t_identity(a, b) + alpha_u_shift * t_compliance(a, b));
      t_tt(a, b) = -det_plastic_f * t_compliance(a, b);
      t_domega(a, k) =
          -det_plastic_f * (t_projection(L, a) * t_work_omega(k, L));
      t_thetaomega(n0, k) =
          -det_plastic_f * (t_volume_projection(L) * t_work_omega(k, L));
      // P and bubble coefficients are reference stresses, so these blocks
      // differentiate Fe Fp, not Fe. Jp cancels the Piola pullback factor.
      FTensor::Tensor3<double, 3, 3, 6> t_reference_fh;
      t_reference_fh(i, j, L) = t_fh(i, k, L) * t_plastic_f(k, j);
      t_fd(i, j, a) = t_reference_fh(i, j, L) * t_projection(L, a);
      t_ftheta(i, j, z) = (t_reference_fh(i, j, L) * t_volume_projection(L)) *
                          t_scalar_component(z);
 
      ++t_compliance;
      ++t_volume;
      ++t_dd;
      ++t_dtheta;
      ++t_dt;
      ++t_td;
      ++t_thetatheta;
      ++t_tt;
      ++t_domega;
      ++t_thetaomega;
      ++t_fd;
      ++t_ftheta;
      ++t_work;
      ++t_work_omega;
      ++t_fh;
      ++t_adjoint_p;
      ++t_eigen_values;
      ++t_eigen_vectors;
      ++t_plastic_f;
    }
    for (const auto &values :
         {materialData->tangentDeviator, materialData->tangentDeviatorJacobian,
          materialData->tangentDeviatorStress,
          materialData->tangentStressDeviator, materialData->tangentJacobian,
          materialData->tangentStress, materialData->deviatorRotation,
          materialData->jacobianRotation, materialData->deformationDeviator,
          materialData->deformationJacobian})
      CHKERR validateValues(*values);
    MoFEMFunctionReturn(0);
  }
};
 
/** Scalar shape functions with a rectangular material coefficient block. */
template <int RowDim, int ColDim>
struct OpAuxiliaryCoefficientMass : public OpBase {
  OpAuxiliaryCoefficientMass(const std::string &row_field,
                             const std::string &col_field, MatrixPtr values,
                             const bool assemble_transpose = true)
      : OpBase(row_field, col_field, OPROWCOL), valuesPtr(std::move(values)) {
    sYmm = row_field == col_field;
    assembleTranspose = assemble_transpose && !sYmm;
  }
 
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data) override {
    MoFEMFunctionBegin;
    FTensor::Index<'a', RowDim> a;
    FTensor::Index<'b', ColDim> b;
    auto get_values =
        MatrixSizeHelper<GetFTensor2FromMatType<RowDim, ColDim, -1, DL>,
                         DL>::get(*valuesPtr, nbIntegrationPts);
    auto t_values = get_values();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_row = row_data.getFTensor0N();
    for (int gg = 0; gg != nbIntegrationPts; ++gg) {
      const double alpha = t_w * getMeasure();
      int rr = 0;
      for (; rr != nbRows / RowDim; ++rr) {
        auto t_col = col_data.getFTensor0N(gg, 0);
        auto t_m = getLocMat<RowDim, ColDim>(RowDim * rr);
        for (int cc = 0; cc != nbCols / ColDim; ++cc) {
          t_m(a, b) += (alpha * t_row * t_col) * t_values(a, b);
          ++t_m;
          ++t_col;
        }
        ++t_row;
      }
      for (; rr < nbRowBaseFunctions; ++rr)
        ++t_row;
      ++t_values;
      ++t_w;
    }
    MoFEMFunctionReturn(0);
  }
 
private:
  MatrixPtr valuesPtr;
};
 
template <int BASE_DIM, int FIELD_DIM>
using OpBaseTensorBase = FormsIntegrators<VolUserDataOperator>::Assembly<
    A>::BiLinearForm<GAUSS>::OpBaseTensorBase<BASE_DIM, FIELD_DIM, SPACE_DIM>;
 
template <int Dim>
using OpAuxiliaryResidual = FormsIntegrators<VolUserDataOperator>::Assembly<
    A>::LinearForm<GAUSS>::OpBaseTimesVector<1, Dim, -1>;
 
/** D_m is computed by the material map rather than interpolated as a field. */
struct OpFormatDeviator : public VolUserDataOperator {
  OpFormatDeviator(const std::string &name, MatrixPtr coordinates,
                   MaterialPostProcData &output)
      : VolUserDataOperator(NOSPACE, OPSPACE),
        coordinatesPtr(std::move(coordinates)),
        tensorPtr(boost::make_shared<MatrixDouble>()) {
    output.symmetricFields[name] = tensorPtr;
  }
 
  MoFEMErrorCode doWork(int, EntityType, EntData &) override {
    MoFEMFunctionBegin;
    const int nb_gauss = getGaussPts().size2();
    auto get_coordinates =
        MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::get(
            *coordinatesPtr, nb_gauss);
    auto get_tensor =
        MatrixSizeHelper<GetFTensor2SymmetricFromMatType<3, -1, DL>, DL>::size(
            *tensorPtr, nb_gauss);
    auto t_coordinates = get_coordinates();
    auto t_tensor = get_tensor();
    FTENSOR_INDEXES(3, i, j);
    for (int gg = 0; gg != nb_gauss; ++gg) {
      const auto t_deviator =
          Tensor2SymmetricDeviatorBasis::getTensor(t_coordinates);
      t_tensor(i, j) = t_deviator(i, j);
      ++t_coordinates;
      ++t_tensor;
    }
    MoFEMFunctionReturn(0);
  }
 
private:
  MatrixPtr coordinatesPtr;
  MatrixPtr tensorPtr;
};
 
struct OpFormatMaterial : public VolUserDataOperator {
  OpFormatMaterial(boost::shared_ptr<DataAtIntegrationPts> data,
                   MaterialPostProcData &output)
      : VolUserDataOperator(NOSPACE, OPSPACE), dataAtPts(std::move(data)) {
    physicalEnergyPtr = boost::make_shared<VectorDouble>();
    mixedEnergyPtr = boost::make_shared<VectorDouble>();
    jacobianPtr = boost::make_shared<VectorDouble>();
    output.scalarFields["PhysicalEnergy"] = physicalEnergyPtr;
    output.scalarFields["MixedStoredEnergy"] = mixedEnergyPtr;
    output.scalarFields["J"] = jacobianPtr;
  }
 
  MoFEMErrorCode doWork(int, EntityType, EntData &) override {
    MoFEMFunctionBegin;
    const int nb_gauss = getGaussPts().size2();
    const auto fields = dataAtPts->auxiliaryData;
    const auto material = dataAtPts->auxiliaryMaterialData;
    auto get_energy =
        MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::get(
            *material->energy, nb_gauss);
    auto get_volume =
        MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::get(
            *material->volume, nb_gauss);
    for (const auto &values_ptr :
         {physicalEnergyPtr, mixedEnergyPtr, jacobianPtr})
      values_ptr->resize(nb_gauss, false);
    CHKERR evaluateAuxiliaryMixedStoredEnergy(*fields, *material,
                                              *mixedEnergyPtr);
    dataAtPts->energyAtPts.resize(nb_gauss, false);
    auto t_material_energy = get_energy();
    auto t_volume = get_volume();
    auto t_physical = getFTensor0FromVec(*physicalEnergyPtr);
    auto t_jacobian = getFTensor0FromVec(*jacobianPtr);
    auto t_energy = getFTensor0FromVec(dataAtPts->energyAtPts);
    for (int gg = 0; gg != nb_gauss; ++gg) {
      const double physical_energy =
          t_material_energy(MaterialData::DEVIATORIC) +
          t_material_energy(MaterialData::VOLUMETRIC);
      if (!std::isfinite(physical_energy))
        SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                "Logarithmic material diagnostic output is not representable");
      t_physical = physical_energy;
      t_energy = physical_energy;
      t_jacobian = t_volume(MaterialData::JACOBIAN);
      ++t_material_energy;
      ++t_volume;
      ++t_physical;
      ++t_energy;
      ++t_jacobian;
    }
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
  VectorPtr physicalEnergyPtr;
  VectorPtr mixedEnergyPtr;
  VectorPtr jacobianPtr;
};
 
} // namespace
 
MoFEMErrorCode pushAuxiliaryLogarithmicFields(
    const EshelbianCore &ep, boost::ptr_deque<UserDataOperator> &pipeline,
    boost::shared_ptr<AuxiliaryLogarithmicStressData> values,
    MatrixPtr reconstructed_log_stretch, SmartPetscObj<Vec> state) {
  MoFEMFunctionBegin;
  if (!ep.physicalEquations->getFeatures().test(
          PhysicalEquations::AUXILIARY_LOGARITHMIC_STRESS) ||
      !values || !reconstructed_log_stretch)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Auxiliary field evaluation requires the auxiliary layout and "
            "valid output storage");
  pipeline.push_back(new OpCalculateTensor2SymmetricDeviatorFieldValues<3>(
      ep.logDeviator, values->logDeviatorTensor, state, MBTET, 0,
      values->logDeviator));
  pipeline.push_back(new OpCalculateVectorFieldValues<1>(
      ep.logJacobian, values->logJacobian, state, MBTET));
  pipeline.push_back(new OpCalculateTensor2SymmetricDeviatorFieldValues<3>(
      ep.auxiliaryLogStress, values->stressTensor, state, MBTET, 0,
      values->stress));
  pipeline.push_back(new OpReconstructAuxiliaryLogarithmicStretch(
      std::move(values), std::move(reconstructed_log_stretch)));
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode evaluateAuxiliaryMixedStoredEnergy(
    const AuxiliaryLogarithmicStressData &fields,
    const AuxiliaryLogarithmicStressMaterialData &material,
    VectorDouble &mixed_energy) {
  MoFEMFunctionBegin;
  const int nb_gauss = fields.logDeviator->size1();
  auto get_d = MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::get(
      *fields.logDeviator, nb_gauss);
  auto get_stress =
      MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::get(
          *material.elasticStress, nb_gauss);
  auto get_material_energy =
      MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::get(
          *material.energy, nb_gauss);
  mixed_energy.resize(nb_gauss, false);
  auto t_d = get_d();
  auto t_stress = get_stress();
  auto t_material_energy = get_material_energy();
  auto t_mixed = getFTensor0FromVec(mixed_energy);
  FTensor::Index<'a', 5> a;
  for (int gg = 0; gg != nb_gauss; ++gg) {
    const double energy = t_stress(a) * t_d(a) -
                          t_material_energy(MaterialData::CONJUGATE) +
                          t_material_energy(MaterialData::VOLUMETRIC);
    if (!std::isfinite(energy))
      SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
              "Mixed stored-energy contribution is not representable");
    t_mixed = energy;
    ++t_d;
    ++t_stress;
    ++t_material_energy;
    ++t_mixed;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode pushAuxiliaryLogarithmicMaterialEvaluation(
    const EshelbianCore &ep,
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
    boost::shared_ptr<DataAtIntegrationPts> data,
    boost::shared_ptr<AuxiliaryLogarithmicStressMaterialData> material_data,
    AuxiliaryLogarithmicMaterialParameters parameters, const bool lhs) {
  MoFEMFunctionBegin;
  if (!data || !data->auxiliaryData || !material_data || !parameters)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Auxiliary material operators require field data, material data "
            "and an element parameter callback");
  if (!ep.physicalEquations->getFeatures().test(
          PhysicalEquations::AUXILIARY_LOGARITHMIC_STRESS) ||
      EshelbianCore::gradApproximator != NO_H1_CONFIGURATION ||
      EshelbianCore::symmetrySelector != SYMMETRIC)
    SETERRQ(ep.mField.get_comm(), MOFEM_NOT_IMPLEMENTED,
            "Auxiliary material operators require the symmetric "
            "no_h1 formulation");
  // Mesh-only diagnostics have no TS rate vector. The existing evaluator
  // clears their rate data, rather than reusing scratch from a previous solve.
  if (ep.alphaU != 0.)
    pipeline.push_back(new OpCalculateVectorFieldValuesDot<5>(
        ep.logDeviator, data->auxiliaryData->logDeviatorDot, MBTET, false));
  pipeline.push_back(new OpCalculateAuxiliaryLogarithmicMaterial(
      std::move(data), std::move(material_data), std::move(parameters),
      ep.alphaU, lhs));
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode pushAuxiliaryLogarithmicMaterialOps(
    const EshelbianCore &ep,
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
    boost::shared_ptr<AuxiliaryLogarithmicStressMaterialData> material_data,
    const bool lhs) {
  MoFEMFunctionBegin;
  if (!lhs) {
    pipeline.push_back(new OpAuxiliaryResidual<5>(
        ep.logDeviator, material_data->residualDeviator));
    pipeline.push_back(new OpAuxiliaryResidual<1>(
        ep.logJacobian, material_data->residualJacobian));
    pipeline.push_back(new OpAuxiliaryResidual<5>(
        ep.auxiliaryLogStress, material_data->residualStress));
    MoFEMFunctionReturnHot(0);
  }
 
  pipeline.push_back(new OpAuxiliaryCoefficientMass<5, 5>(
      ep.logDeviator, ep.logDeviator, material_data->tangentDeviator));
  pipeline.push_back(new OpAuxiliaryCoefficientMass<5, 1>(
      ep.logDeviator, ep.logJacobian, material_data->tangentDeviatorJacobian));
  pipeline.push_back(new OpAuxiliaryCoefficientMass<1, 1>(
      ep.logJacobian, ep.logJacobian, material_data->tangentJacobian));
  pipeline.push_back(new OpAuxiliaryCoefficientMass<5, 5>(
      ep.logDeviator, ep.auxiliaryLogStress,
      material_data->tangentDeviatorStress, false));
  pipeline.push_back(new OpAuxiliaryCoefficientMass<5, 5>(
      ep.auxiliaryLogStress, ep.logDeviator,
      material_data->tangentStressDeviator, false));
  pipeline.push_back(new OpAuxiliaryCoefficientMass<5, 5>(
      ep.auxiliaryLogStress, ep.auxiliaryLogStress,
      material_data->tangentStress));
  pipeline.push_back(new OpAuxiliaryCoefficientMass<5, 3>(
      ep.logDeviator, ep.rotAxis, material_data->deviatorRotation));
  pipeline.push_back(new OpAuxiliaryCoefficientMass<1, 3>(
      ep.logJacobian, ep.rotAxis, material_data->jacobianRotation));
  // The mixed potential contributes -P:F, including the tensor bubble field.
  constexpr auto work_sign = []() constexpr { return -1.; };
  pipeline.push_back(new OpBaseTensorBase<3, 5>(
      ep.logDeviator, ep.piolaStress, material_data->deformationDeviator,
      work_sign, true));
  pipeline.push_back(new OpBaseTensorBase<3, 1>(
      ep.logJacobian, ep.piolaStress, material_data->deformationJacobian,
      work_sign, true));
  pipeline.push_back(new OpBaseTensorBase<9, 5>(
      ep.logDeviator, ep.bubbleField, material_data->deformationDeviator,
      work_sign, true));
  pipeline.push_back(new OpBaseTensorBase<9, 1>(
      ep.logJacobian, ep.bubbleField, material_data->deformationJacobian,
      work_sign, true));
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode pushAuxiliaryLogarithmicStatePostProc(
    boost::ptr_deque<UserDataOperator> &,
    boost::shared_ptr<AuxiliaryLogarithmicStressData> values,
    const std::string &prefix, MaterialPostProcData &output) {
  MoFEMFunctionBegin;
  if (!values)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Auxiliary postprocessing requires evaluated material fields");
  output.symmetricFields[prefix + "D"] = values->logDeviatorTensor;
  output.symmetricFields[prefix + "Td"] = values->stressTensor;
  output.scalarFields[prefix + "theta"] = values->logJacobian;
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
pushAuxiliaryLogarithmicPostProc(boost::ptr_deque<UserDataOperator> &pipeline,
                                 boost::shared_ptr<DataAtIntegrationPts> data,
                                 MaterialPostProcData &output) {
  MoFEMFunctionBegin;
  if (!data || !data->auxiliaryData || !data->auxiliaryMaterialData)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Auxiliary postprocessing requires field and material producers");
  CHKERR pushAuxiliaryLogarithmicStatePostProc(pipeline, data->auxiliaryData,
                                               "", output);
  pipeline.push_back(new OpFormatDeviator(
      "D_m", data->auxiliaryMaterialData->materialDeviator, output));
  pipeline.push_back(new OpFormatMaterial(std::move(data), output));
  MoFEMFunctionReturn(0);
}
 
} // namespace EshelbianPlasticity