NeoHookeanLogarithmicMaterial.cpp

Go to the documentation of this file.

/**
 * @file NeoHookeanLogarithmicMaterial.cpp
 * @brief Unique auxiliary logarithmic-stress inverse and its derivatives.
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <MatrixFunction.hpp>
#include <NeoHookeanLogarithmicMaterial.hpp>
#include <lapack_wrap.h>
 
#include <algorithm>
#include <array>
#include <cmath>
#include <limits>
 
namespace EshelbianPlasticity {
namespace {
 
using Material = NeoHookeanLogarithmicMaterial;
using PrincipalValues = FTensor::Tensor1<double, 3>;
using Eigenvectors = FTensor::Tensor2<double, 3, 3>;
 
MoFEMErrorCode validateCoordinates(const Material::Coordinates &t_values) {
  MoFEMFunctionBegin;
  for (int aa = 0; aa != 5; ++aa) {
    if (!std::isfinite(t_values(aa)))
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Neo-Hookean logarithmic coordinates must be finite");
  }
  MoFEMFunctionReturn(0);
}
 
// Scale before LAPACK so its absolute tolerances also resolve very small data.
MoFEMErrorCode evaluateEigenvectors(const Material::Coordinates &t_values,
                                    PrincipalValues &t_eigenvalues,
                                    Eigenvectors &t_eigenvectors,
                                    double &scale) {
  MoFEMFunctionBegin;
  FTENSOR_INDEXES(3, i, j);
  FTensor::Index<'A', 5> A;
  CHKERR validateCoordinates(t_values);
  scale = 0.;
  for (int aa = 0; aa != 5; ++aa)
    scale = std::max(scale, std::abs(t_values(aa)));
  if (scale == 0.) {
    t_eigenvalues(i) = 0.;
    constexpr auto t_identity = FTensor::Kronecker_Delta<int>();
    t_eigenvectors(i, j) = t_identity(i, j);
  } else {
    Material::Coordinates t_scaled;
    t_scaled(A) = t_values(A) / scale;
    t_eigenvectors(i, j) = Material::getTensor(t_scaled)(i, j);
    CHKERR computeEigenValuesSymmetric(t_eigenvectors, t_eigenvalues);
    t_eigenvalues(i) -=
        t_eigenvalues(0) / 3. + t_eigenvalues(1) / 3. + t_eigenvalues(2) / 3.;
  }
  MoFEMFunctionReturn(0);
}
 
double getLogSum(const double a, const double b) {
  const double maximum = std::max(a, b);
  return maximum + std::log1p(std::exp(std::min(a, b) - maximum));
}
 
MoFEMErrorCode evaluateEnergy(const double mu,
                              const PrincipalValues &t_log_values,
                              double &energy) {
  MoFEMFunctionBegin;
  energy = 0.;
  for (int aa = 0; aa != 3; ++aa) {
    const double x = t_log_values(aa);
    if (std::abs(x) <= 1.e-3) {
      // Scale before squaring so a small strain retains its quadratic energy.
      const double remainder =
          0.5 + x * (1. / 6. + x * (1. / 24. + x * (1. / 120. + x / 720.)));
      energy += ((mu / 2.) * x) * x * remainder;
    } else if (x <= 1.) {
      energy += (mu / 2.) * (std::expm1(x) - x);
    } else {
      // Combine exp(x)-1-x before exponentiating a potentially large term.
      energy += std::exp(std::log(mu) + x - std::log(2.) +
                         std::log1p(-(1. + x) * std::exp(-x)));
    }
  }
  if (!std::isfinite(energy))
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "Neo-Hookean deviatoric energy is not representable");
  MoFEMFunctionReturn(0);
}
 
Material::SymmetricTensor getPrincipalMode(const Eigenvectors &t_vectors,
                                           const int aa, const int bb) {
  FTENSOR_INDEXES(3, i, j);
  Material::SymmetricTensor t_mode;
  t_mode(i, j) = (t_vectors(aa, i) * t_vectors(bb, j) ||
                  t_vectors(bb, i) * t_vectors(aa, j)) /
                 (aa == bb ? 2. : std::sqrt(2.));
  return t_mode;
}
 
MoFEMErrorCode addTangentMode(Material::Tangent &t_tangent,
                              const Material::SymmetricTensor &t_mode,
                              const double coefficient) {
  MoFEMFunctionBegin;
  FTensor::Index<'A', 5> A;
  FTensor::Index<'B', 5> B;
  if (!std::isfinite(coefficient) || coefficient <= 0.)
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "Neo-Hookean material derivative is not representable");
  const auto t_coordinates = Material::getCoordinates(t_mode);
  t_tangent(A, B) += coefficient * t_coordinates(A) * t_coordinates(B);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode validateTangent(const Material::Tangent &t_tangent) {
  MoFEMFunctionBegin;
  for (int aa = 0; aa != 5; ++aa)
    for (int bb = 0; bb != 5; ++bb)
      if (!std::isfinite(t_tangent(aa, bb)))
        SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                "Neo-Hookean material derivative is not representable");
  MoFEMFunctionReturn(0);
}
 
} // namespace
 
MoFEMErrorCode Material::validateParameters(const Parameters &parameters) {
  MoFEMFunctionBegin;
  if (!std::isfinite(parameters.c10) || parameters.c10 <= 0. ||
      !std::isfinite(2. * parameters.c10) ||
      !std::isfinite(parameters.bulkModulus) || parameters.bulkModulus <= 0.)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Neo-Hookean C10 and K must be finite and positive, with finite "
            "mu = 2*C10");
  MoFEMFunctionReturn(0);
}
 
Material::SymmetricTensor Material::getTensor(const Coordinates &t_coordinates,
                                              const double theta) {
  return Tensor2SymmetricDeviatorBasis::getTensor(t_coordinates, theta);
}
 
Material::Coordinates
Material::getCoordinates(const SymmetricTensor &t_tensor) {
  return Tensor2SymmetricDeviatorBasis::getCoordinates(t_tensor);
}
 
MoFEMErrorCode Material::evaluateDeviator(const Parameters &parameters,
                                          const Coordinates &t_deviator,
                                          double &energy,
                                          Coordinates *stress_ptr,
                                          Tangent *hessian_ptr) {
  MoFEMFunctionBegin;
  FTENSOR_INDEXES(3, i, j);
  FTensor::Index<'A', 5> A;
  FTensor::Index<'B', 5> B;
#ifndef NDEBUG
  CHKERR validateParameters(parameters);
#endif
  const double mu = 2. * parameters.c10;
  const double log_mu = std::log(mu);
  PrincipalValues t_log_values;
  Eigenvectors t_vectors;
  double scale;
  CHKERR evaluateEigenvectors(t_deviator, t_log_values, t_vectors, scale);
  t_log_values(i) *= 2. * scale;
  for (int aa = 0; aa != 3; ++aa)
    if (!std::isfinite(t_log_values(aa)))
      SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
              "Neo-Hookean principal logarithm is not representable");
  CHKERR evaluateEnergy(mu, t_log_values, energy);
 
  if (stress_ptr) {
    // Remove the spherical part before restoring the stress scale. expm1
    // retains small differences even when all principal stretches are near one.
    const double maximum =
        std::max({t_log_values(0), t_log_values(1), t_log_values(2)});
    PrincipalValues t_deviatoric_values;
    for (int aa = 0; aa != 3; ++aa)
      t_deviatoric_values(aa) = std::expm1(t_log_values(aa) - maximum);
    const double mean = t_deviatoric_values(0) / 3. +
                        t_deviatoric_values(1) / 3. +
                        t_deviatoric_values(2) / 3.;
    t_deviatoric_values(i) -= mean;
    auto principal_stress = [log_mu, maximum](const double value) {
      return value == 0. ? 0.
                         : std::copysign(std::exp(log_mu + maximum +
                                                  std::log(std::abs(value))),
                                         value);
    };
    const auto t_stress =
        EigenMatrix::getMat(t_deviatoric_values, t_vectors, principal_stress);
    *stress_ptr = getCoordinates(t_stress);
    CHKERR validateCoordinates(*stress_ptr);
  }
 
  if (hessian_ptr) {
    auto &t_hessian = *hessian_ptr;
    t_hessian(A, B) = 0.;
    for (int aa = 0; aa != 3; ++aa) {
      CHKERR addTangentMode(t_hessian, getPrincipalMode(t_vectors, aa, aa),
                            std::exp(std::log(2.) + log_mu + t_log_values(aa)));
      for (int bb = 0; bb != aa; ++bb) {
        const double difference = std::abs(t_log_values(aa) - t_log_values(bb));
        const double factor =
            difference == 0. ? 1. : -std::expm1(-difference) / difference;
        const double coefficient = std::exp(
            std::log(2.) + log_mu +
            std::max(t_log_values(aa), t_log_values(bb)) + std::log(factor));
        CHKERR addTangentMode(t_hessian, getPrincipalMode(t_vectors, aa, bb),
                              coefficient);
      }
    }
    CHKERR validateTangent(t_hessian);
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode Material::evaluateVolume(const Parameters &parameters,
                                        const double theta,
                                        VolumeState &state) {
  MoFEMFunctionBegin;
#ifndef NDEBUG
  CHKERR validateParameters(parameters);
#endif
  if (!std::isfinite(theta))
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Neo-Hookean logarithmic volume must be finite");
  state.jacobian = std::exp(theta);
  const double jacobian_minus_one = std::expm1(theta);
  state.energy =
      0.5 * (parameters.bulkModulus * jacobian_minus_one) * jacobian_minus_one;
  const double k_j = parameters.bulkModulus * state.jacobian;
  state.firstDerivative = k_j * jacobian_minus_one;
  state.secondDerivative = state.firstDerivative + k_j * state.jacobian;
  if (!std::isfinite(state.jacobian) || state.jacobian <= 0. ||
      !std::isfinite(state.energy) || !std::isfinite(state.firstDerivative) ||
      !std::isfinite(state.secondDerivative))
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "Neo-Hookean volumetric state is not representable");
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode Material::evaluateInverse(const Parameters &parameters,
                                         const Coordinates &t_stress,
                                         InverseState &state) {
  MoFEMFunctionBegin;
  FTENSOR_INDEXES(3, i, j);
  FTensor::Index<'A', 5> A;
  FTensor::Index<'B', 5> B;
#ifndef NDEBUG
  CHKERR validateParameters(parameters);
#endif
  const double mu = 2. * parameters.c10;
  const double log_mu = std::log(mu);
  PrincipalValues t_values;
  Eigenvectors t_vectors;
  double scale;
  CHKERR evaluateEigenvectors(t_stress, t_values, t_vectors, scale);
 
  // z = log(zeta/mu) retains the positive bracket even when zeta underflows.
  const double minimum = std::min({t_values(0), t_values(1), t_values(2)});
  PrincipalValues t_log_gaps;
  double residual_at_zero = 0.;
  for (int aa = 0; aa != 3; ++aa) {
    const double difference = t_values(aa) - minimum;
    t_log_gaps(aa) = difference == 0. || scale == 0.
                         ? -std::numeric_limits<double>::infinity()
                         : std::log(difference) + std::log(scale) - log_mu;
    residual_at_zero += getLogSum(t_log_gaps(aa), 0.);
  }
  double lower = -residual_at_zero;
  double upper = 0.;
  double z = lower / 3.;
  PrincipalValues t_log_values;
  constexpr int maximum_iterations = 100;
  bool converged = false;
  for (int iteration = 0; iteration != maximum_iterations; ++iteration) {
    double residual = 0.;
    double derivative = 0.;
    double residual_scale = 1.;
    for (int aa = 0; aa != 3; ++aa) {
      t_log_values(aa) = getLogSum(t_log_gaps(aa), z);
      residual += t_log_values(aa);
      derivative += std::exp(z - t_log_values(aa));
      residual_scale = std::max(residual_scale, std::abs(t_log_values(aa)));
    }
    state.inverseResidual = residual;
    state.inverseIterations = iteration + 1;
    if (std::abs(residual) <=
        32. * std::numeric_limits<double>::epsilon() * residual_scale) {
      converged = true;
      break;
    }
    if (residual > 0.)
      upper = z;
    else
      lower = z;
    const double newton = z - residual / derivative;
    z = newton > lower && newton < upper ? newton : (lower + upper) / 2.;
  }
  if (!converged)
    SETERRQ(
        PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
        "Neo-Hookean auxiliary logarithmic-stress inverse did not converge");
 
  const double mean_log = state.inverseResidual / 3.;
  t_log_values(i) -= mean_log;
  const auto t_material_deviator = EigenMatrix::getMat(
      t_log_values, t_vectors, [](const double value) { return value / 2.; });
  state.tMaterialDeviator = getCoordinates(t_material_deviator);
  CHKERR evaluateEnergy(mu, t_log_values, state.deviatoricEnergy);
  state.conjugateEnergy =
      t_stress(A) * state.tMaterialDeviator(A) - state.deviatoricEnergy;
  if (!std::isfinite(state.conjugateEnergy))
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "Neo-Hookean conjugate energy is not representable");
 
  // Differentiate Td = mu*dev(exp(2*Dm)), including the scalar trace
  // constraint. Its diagonal compliance is half the weighted variance of
  // principal tensor increments. Merge the two smallest weights first: the
  // within-pair variance and the variance of that pair's mean against the third
  // principal direction give two independent modes without subtracting large
  // matrices.
  PrincipalValues t_log_weights;
  t_log_weights(i) = -log_mu - t_log_values(i);
  const double log_weight_sum = getLogSum(
      getLogSum(t_log_weights(0), t_log_weights(1)), t_log_weights(2));
  auto &t_compliance = state.tCompliance;
  t_compliance(A, B) = 0.;
  std::array<int, 3> weight_order{0, 1, 2};
  std::sort(weight_order.begin(), weight_order.end(),
            [&t_log_weights](const int aa, const int bb) {
              return t_log_weights(aa) < t_log_weights(bb);
            });
  const int first = weight_order[0];
  const int second = weight_order[1];
  const int third = weight_order[2];
  const double log_pair_sum =
      getLogSum(t_log_weights(first), t_log_weights(second));
  const auto t_first_mode = getPrincipalMode(t_vectors, first, first);
  const auto t_second_mode = getPrincipalMode(t_vectors, second, second);
  const auto t_third_mode = getPrincipalMode(t_vectors, third, third);
  SymmetricTensor t_difference;
  t_difference(i, j) = t_first_mode(i, j) - t_second_mode(i, j);
  CHKERR addTangentMode(t_compliance, t_difference,
                        std::exp(-std::log(2.) + t_log_weights(first) +
                                 t_log_weights(second) - log_pair_sum));
  t_difference(i, j) =
      std::exp(t_log_weights(first) - log_pair_sum) * t_first_mode(i, j) +
      std::exp(t_log_weights(second) - log_pair_sum) * t_second_mode(i, j) -
      t_third_mode(i, j);
  CHKERR addTangentMode(t_compliance, t_difference,
                        std::exp(-std::log(2.) + log_pair_sum +
                                 t_log_weights(third) - log_weight_sum));
  for (int aa = 0; aa != 3; ++aa) {
    for (int bb = 0; bb != aa; ++bb) {
      const double difference = std::abs(t_log_values(aa) - t_log_values(bb));
      const double factor =
          difference == 0. ? 1. : difference / -std::expm1(-difference);
      const double shear_coefficient = std::exp(
          -std::log(2.) - log_mu -
          std::max(t_log_values(aa), t_log_values(bb)) + std::log(factor));
      CHKERR addTangentMode(t_compliance, getPrincipalMode(t_vectors, aa, bb),
                            shear_coefficient);
    }
  }
  CHKERR validateCoordinates(state.tMaterialDeviator);
  CHKERR validateTangent(t_compliance);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode Material::evaluateFenchelGap(const Parameters &parameters,
                                            const Coordinates &t_deviator,
                                            const Coordinates &t_stress,
                                            const InverseState &inverse,
                                            double &gap,
                                            double *energy_ptr) {
  MoFEMFunctionBegin;
  FTensor::Index<'A', 5> A;
  CHKERR validateCoordinates(t_stress);
  double energy;
  CHKERR evaluateDeviator(parameters, t_deviator, energy);
  gap = energy + inverse.conjugateEnergy - t_stress(A) * t_deviator(A);
  if (!std::isfinite(gap))
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "Neo-Hookean Fenchel gap is not representable");
  if (energy_ptr)
    *energy_ptr = energy;
  MoFEMFunctionReturn(0);
}
 
} // namespace EshelbianPlasticity