LogarithmicDirectionalDerivatives.hpp

Go to the documentation of this file.

/**
 * \file LogarithmicDirectionalDerivatives.hpp
 * \brief Mixed directional derivatives of matrix exponentials for the
 * reduced Hessian of the logarithmic plastic parametrisation.
 *
 * A jet, Matrix<N>, holds exp(A) together with its derivatives in up to N = 3
 * directions. Block `mask` is the mixed derivative in the directions whose
 * bits are set, block 0 is the value; for N = 3 with directions (p, s, e),
 * block 3 is D^2 exp(A)[p, s] and block 7 is D^3 exp(A)[p, s, e]. Jets form
 * an algebra: product() of the jets of two factors is the jet of their
 * product, which gives the mixed derivatives of F = R U F^p from the jets of
 * the three exponentials without differentiating the product by hand.
 *
 * The exponentials themselves are closed-form:
 * - SymmetricExponential, for symmetric arguments (stretches), uses the
 *   eigen-decomposition. The value is EigenMatrix::getMat
 *   (MatrixFunction.hpp), the spectral form of isotropic tensor functions
 *   used in computational plasticity \cite miehe1998comparison
 *   \cite miehe2001algorithms. The derivatives are the divided-difference
 *   (Daleckii--Krein) representation of the Fr\'echet derivatives in the
 *   eigenbasis \cite daleckii1965integration, \cite bhatia1997matrix
 *   (Chapter V), \cite higham2008functions (Chapter 3), extended to second
 *   and third order as in \cite higham2014higher. The divided differences of
 *   exp at the eigenvalues are evaluated from the series in the complete
 *   homogeneous symmetric polynomials of the shifts from their mean when the
 *   spread is small, and recursively otherwise, which is the stable
 *   treatment of nearly coincident points \cite mccurdy1984accurate,
 *   \cite deboor2005divided, \cite opitz1964steigungsmatrizen.
 * - RotationExponential, for exp(skew(w)), uses Rodrigues' formula and its
 *   derivatives from LieGroups::SO3 (Lie.hpp) \cite sola2018micro, with the
 *   small-angle series for the coefficient functions
 *   \cite grassia1998practical.
 * No series is used for the matrix exponential itself; for the background
 * on why series and scaling-and-squaring are a poor fit for derivatives see
 * \cite moler2003nineteen and \cite higham2008functions (Chapter 10).
 *
 * Usage. Directions are written into the single-bit blocks of an argument
 * jet, the base into block 0, and fill() writes the complete jet:
 * \code
 * using namespace LogarithmicDirectionalDerivatives;
 * FTENSOR_INDEXES(3, i, j);
 *
 * Matrix<3> u, f, r;                       // arguments
 * u[0](i, j) = t_H(i, j);                   // base: total log stretch
 * u[1](i, j) = t_ph(i, j);                  // direction p (adjoint)
 * u[2](i, j) = t_sh(i, j);                  // direction s
 * u[4](i, j) = t_e(i, j);                   // direction e (basis function)
 * f[0](i, j) = t_Hp(i, j);                  // plastic log stretch
 * f[2](i, j) = t_vh(i, j);                  // only a direction s
 * r[0](i, j) = skew(t_omega)(i, j);         // rotation vector as skew block
 * r[1](i, j) = skew(t_pr)(i, j);
 *
 * const SymmetricExponential stretch(u[0]), plastic(f[0]);
 * const RotationExponential rotation(axial(r[0]));
 * Matrix<3> U, Fp, R;
 * stretch.fill(u, U);      // U[7] = D^3 exp(H)[ph, sh, e]
 * plastic.fill(f, Fp);     // blocks with bit 1 or 4 are zero
 * rotation.fill(r, R);     // derivatives in the axial vectors
 *
 * const auto F = product(product(R, U), Fp); // jet of F = R U Fp
 * double d2 = t_P(i, j) * F[3](i, j);        // P : D^2 F[p, s]
 * double d3 = t_P(i, j) * F[7](i, j);        // P : D^3 F[p, s, e]
 *
 * // exp(-Hp) from the same eigenvectors, and a stand-alone derivative
 * const auto inverse = plastic.negated();
 * const auto dexp = stretch.first(t_e);     // D exp(H)[e]
 *
 * // Neo-Hookean deviatoric curvature, f = C10 tr exp(2 dev H):
 * const auto c = deviatorCurvature(c10, t_H, t_ph, t_sh);
 * double D2f = c.second(i, j) * t_e(i, j);  // D^2 f[s, e]
 * double D3f = c.third(i, j) * t_e(i, j);   // D^3 f[p, s, e]
 * \endcode
 * A zero direction block gives zero derivatives, so the same argument jet
 * can be reused with one direction changed. A base object is constructed
 * once per quadrature point and its fill() is called for every direction
 * set; the cost of fill() is a few hundred flops.
 */
 
#ifndef ESHELBIAN_LOGARITHMIC_DIRECTIONAL_DERIVATIVES_HPP
#define ESHELBIAN_LOGARITHMIC_DIRECTIONAL_DERIVATIVES_HPP
 
#include <Lie.hpp>
#include <MatrixFunction.hpp>
 
namespace EshelbianPlasticity {
 
/** Local directions for D(J_z - R_z^T p)[s,v], with p held fixed. */
struct LogarithmicMaterialCurvature {
  FTensor::Tensor2_symmetric<double, 3> stretch, adjoint, direction;
  FTensor::Tensor1<double, 5> stressAdjoint, stressDirection, stressResult;
  FTensor::Tensor1<double, 6> stretchResult;
};
 
namespace LogarithmicDirectionalDerivatives {
using Tensor = FTensor::Tensor2<double, 3, 3>;
using Symmetric = FTensor::Tensor2_symmetric<double, 3>;
using Vector = FTensor::Tensor1<double, 3>;
 
/** Cursor onto one 3x3 block of a jet; ++ moves to the next block. */
using View = FTensor::Tensor2<FTensor::CursorPtr<double *, 9>, 3, 3>;
 
/** Mixed derivatives in independent nilpotent directions. Bit masks denote
 * subsets of directions; block `mask` is the derivative in the directions
 * whose bits are set and block 0 is the value. All 2^N blocks live in one
 * contiguous array, block `mask` at doubles [9 mask, 9 mask + 9), accessed
 * through FTensor cursor views, so no block is copied and a loop over masks
 * is a cursor increment. */
template <int N> struct Matrix {
  static constexpr int count = 1 << N;
  static constexpr int blockSize = 9;
  static constexpr int size = count * blockSize;
  std::array<double, size> data;
 
  Matrix() { data.fill(0.); }
 
  /** View of block `mask`. Views are handles, so a const jet returns the
   * same view type, as FTensor pointer tensors do. The view is built directly
   * from the pointer: getFTensor2FromPtr goes through a ublas shallow array
   * adaptor whose shared_array allocates, which is prohibitive here. */
  View operator[](int mask) const {
    return View(const_cast<double *>(data.data()) + blockSize * mask);
  }
  /** Cursor at block 0; ++ advances by one block. */
  View begin() const { return (*this)[0]; }
};
 
/** result = a b over the jet algebra; result must be distinct from a and b,
 * which may alias each other. */
template <int N>
inline void product(const Matrix<N> &a, const Matrix<N> &b, Matrix<N> &result) {
  FTENSOR_INDEXES(3, i, j, k);
  auto t_r = result.begin();
  for (int mask = 0; mask != Matrix<N>::count; ++mask, ++t_r) {
    // Accumulate in a local block so the compiler need not assume that the
    // result aliases a or b. Subsets run from the mask itself down to the
    // empty set.
    Tensor block;
    block(i, j) = a[mask](i, k) * b[0](k, j);
    for (int s = (mask - 1) & mask; mask != 0; s = (s - 1) & mask) {
      block(i, j) += a[s](i, k) * b[mask ^ s](k, j);
      if (!s)
        break;
    }
    t_r(i, j) = block(i, j);
  }
}
 
template <int N> Matrix<N> product(const Matrix<N> &a, const Matrix<N> &b) {
  Matrix<N> result;
  product(a, b, result);
  return result;
}
 
template <class V> Tensor skew(const V &v) {
  FTENSOR_INDEXES(3, i, j, k);
  FTensor::Levi_Civita<double> epsilon;
  Tensor result;
  result(i, j) = epsilon(i, j, k) * v(k);
  return result;
}
 
/** Axial vector of a skew block, skew(v)(i,j) = epsilon(i,j,k) v(k). */
inline Vector axial(View block) {
  FTENSOR_INDEXES(3, i, j, k);
  FTensor::Levi_Civita<double> epsilon;
  Vector result;
  result(k) = 0.5 * (epsilon(i, j, k) * block(i, j));
  return result;
}
 
namespace detail {
constexpr int factorialTableSize = 32;
constexpr std::array<double, factorialTableSize> inverseFactorials() {
  std::array<double, factorialTableSize> table{};
  double f = 1.;
  for (int n = 0; n != factorialTableSize; ++n) {
    table[n] = 1. / f;
    f *= n + 1;
  }
  return table;
}
inline constexpr auto inverseFactorial = inverseFactorials();
 
/** Spread below which divided differences of exp use the series. */
constexpr double closeSpread = 0.1;
 
/** Divided difference exp[p_0, ..., p_{K-1}] of close points as the series
 * e^mean sum_n h_n(shifts) / (n + K - 1)! in the complete homogeneous
 * symmetric polynomials h_n of the shifts from the mean; with spread 0.1 the
 * term n is below 1e-17 for n >= 11. Exact for coincident points. */
template <std::size_t K>
double expDividedDifferenceSeries(const std::array<double, K> &points) {
  double mean = 0.;
  for (auto p : points)
    mean += p;
  mean /= double(K);
  constexpr int terms = 12;
  std::array<double, terms> h{};
  h[0] = 1.;
  for (auto p : points) {
    const double d = p - mean;
    for (int n = 1; n != terms; ++n)
      h[n] += d * h[n - 1];
  }
  double sum = 0.;
  for (int n = terms - 1; n >= 0; --n)
    sum += h[n] * inverseFactorial[n + K - 1];
  return std::exp(mean) * sum;
}
} // namespace detail
 
/** Exponential of a symmetric matrix with its directional derivatives up to
 * third order from the eigen-decomposition. The value is
 * EigenMatrix::getMat (MatrixFunction.hpp). The derivatives use the
 * divided-difference expansion in the eigenbasis (components X_ab of a
 * direction, exp[...] divided differences of exp at the eigenvalues):
 *   (D exp[E])_ab      = exp[l_a,l_b] E_ab,
 *   (D^2 exp[E,F])_ab  = sum_c exp[l_a,l_c,l_b] (E_ac F_cb + F_ac E_cb),
 *   (D^3 exp[E,F,G])_ab = sum_{c,d} exp[l_a,l_c,l_d,l_b]
 *                          sum_{perm (X,Y,Z) of (E,F,G)} X_ac Y_cd Z_db.
 * The divided differences are evaluated as a short series when eigenvalues
 * are close, so coincident and nearly coincident eigenvalues are exact; the
 * confluent formulas of EigenMatrix::getDiffMat and getDiffDiffMat apply only
 * to eigenvalues equal to 1e-14 and lose about eps/gap^2 in between. The
 * decomposition and the difference tables are formed once per base and
 * reused for every direction set. */
struct SymmetricExponential {
  template <class A> explicit SymmetricExponential(const A &a) {
    FTENSOR_INDEXES(3, i, j);
    vectors(i, j) = a(i, j);
    CHK_THROW_MESSAGE(MoFEM::computeEigenValuesSymmetric(vectors, values),
                      "Cannot diagonalise symmetric exponential argument");
    value = EigenMatrix::getMat(values, vectors, expFun());
    prepareDifferences();
  }
 
  /** exp(-A) from the same eigenvectors, without a second diagonalisation. */
  SymmetricExponential negated() const {
    SymmetricExponential result(*this);
    FTENSOR_INDEXES(3, i);
    result.values(i) = -values(i);
    result.value = EigenMatrix::getMat(result.values, result.vectors, expFun());
    result.prepareDifferences();
    return result;
  }
 
  /** D exp(A)[e] */
  Symmetric first(const Symmetric &e) const {
    return fromEigenbasis(firstInBasis(toEigenbasis(e)));
  }
 
  /** D^2 exp(A)[e1, e2] */
  Symmetric second(const Symmetric &e1, const Symmetric &e2) const {
    return fromEigenbasis(secondInBasis(toEigenbasis(e1), toEigenbasis(e2)));
  }
 
  /** D^3 exp(A)[e1, e2, e3] */
  Symmetric third(const Symmetric &e1, const Symmetric &e2,
                  const Symmetric &e3) const {
    return fromEigenbasis(
        thirdInBasis(toEigenbasis(e1), toEigenbasis(e2), toEigenbasis(e3)));
  }
 
  /** Symmetric part of a jet block. */
  static Symmetric symmetric(View block) {
    FTENSOR_INDEXES(3, i, j);
    Symmetric result;
    result(i, j) = (block(i, j) || block(j, i)) / 2.;
    return result;
  }
 
  /** Fill the jet of exp at this base: block 0 the value, single-bit blocks
   * the first derivatives in the (symmetric) direction blocks of `argument`,
   * two-bit blocks the second and the three-bit block the third derivative.
   * A zero direction block gives zero derivatives. */
  template <int N>
  void fill(const Matrix<N> &argument, Matrix<N> &result) const {
    FTENSOR_INDEXES(3, i, j);
    result.data.fill(0.);
    result[0](i, j) = value(i, j);
    std::array<Tensor, N> e;
    std::array<bool, N> active{};
    for (int d = 0; d != N; ++d) {
      const auto block = symmetric(argument[1 << d]);
      active[d] = block(i, j) * block(i, j) > 0.;
      if (active[d]) {
        e[d] = toEigenbasis(block);
        result[1 << d](i, j) = fromEigenbasis(firstInBasis(e[d]))(i, j);
      }
    }
    for (int d1 = 0; d1 != N; ++d1)
      for (int d2 = d1 + 1; d2 != N; ++d2)
        if (active[d1] && active[d2])
          result[(1 << d1) | (1 << d2)](i, j) =
              fromEigenbasis(secondInBasis(e[d1], e[d2]))(i, j);
    if constexpr (N == 3)
      if (active[0] && active[1] && active[2])
        result[7](i, j) = fromEigenbasis(thirdInBasis(e[0], e[1], e[2]))(i, j);
  }
 
  FTensor::Tensor1<double, 3> values;
  /** Eigenvectors in rows (MoFEM convention). */
  FTensor::Tensor2<double, 3, 3> vectors;
  Symmetric value;
 
private:
  static const EigenMatrix::Fun<double> &expFun() {
    static const EigenMatrix::Fun<double> f = [](double x) {
      return std::exp(x);
    };
    return f;
  }
 
  Tensor toEigenbasis(const Symmetric &e) const {
    FTENSOR_INDEXES(3, i, j);
    FTensor::Index<'a', 3> a;
    FTensor::Index<'b', 3> b;
    Tensor c;
    c(a, b) = vectors(a, i) * (e(i, j) * vectors(b, j));
    return c;
  }
 
  Symmetric fromEigenbasis(const Tensor &t) const {
    FTENSOR_INDEXES(3, i, j);
    FTensor::Index<'a', 3> a;
    FTensor::Index<'b', 3> b;
    Tensor full;
    full(i, j) = vectors(a, i) * (t(a, b) * vectors(b, j));
    Symmetric result;
    result(i, j) = (full(i, j) || full(j, i)) / 2.;
    return result;
  }
 
  Tensor firstInBasis(const Tensor &x) const {
    Tensor t;
    for (int a = 0; a != 3; ++a)
      for (int b = 0; b != 3; ++b)
        t(a, b) = second_[a][b] * x(a, b);
    return t;
  }
 
  Tensor secondInBasis(const Tensor &x, const Tensor &y) const {
    Tensor t;
    for (int a = 0; a != 3; ++a)
      for (int b = 0; b != 3; ++b) {
        double sum = 0.;
        for (int c = 0; c != 3; ++c)
          sum += third_[a][c][b] * (x(a, c) * y(c, b) + y(a, c) * x(c, b));
        t(a, b) = sum;
      }
    return t;
  }
 
  Tensor thirdInBasis(const Tensor &x, const Tensor &y, const Tensor &z) const {
    Tensor t;
    for (int a = 0; a != 3; ++a)
      for (int b = 0; b != 3; ++b) {
        double sum = 0.;
        for (int c = 0; c != 3; ++c)
          for (int d = 0; d != 3; ++d)
            sum += fourth_[a][c][d][b] *
                   (x(a, c) * y(c, d) * z(d, b) + x(a, c) * z(c, d) * y(d, b) +
                    y(a, c) * x(c, d) * z(d, b) + y(a, c) * z(c, d) * x(d, b) +
                    z(a, c) * x(c, d) * y(d, b) + z(a, c) * y(c, d) * x(d, b));
        t(a, b) = sum;
      }
    return t;
  }
 
  /** Divided differences of exp at the eigenvalues for every index tuple.
   * The eigenvalues are ascending (LAPACK dsyev), so for a sorted tuple the
   * extreme points are its first and last entry and the recursion
   * exp[p_0..p_k] = (exp[p_1..p_k] - exp[p_0..p_{k-1}]) / (p_k - p_0) divides
   * by the largest spread; tuples with a small spread use the series. */
  void prepareDifferences() {
#ifndef NDEBUG
    if (values(0) > values(1) || values(1) > values(2))
      THROW_MESSAGE("Eigenvalues are expected in ascending order");
#endif
    const auto l = [&](int a) { return values(a); };
    const auto close = [](double spread) {
      return spread <= detail::closeSpread;
    };
    const double e[3] = {std::exp(l(0)), std::exp(l(1)), std::exp(l(2))};
    for (int a = 0; a != 3; ++a)
      for (int b = a; b != 3; ++b)
        second_[a][b] = second_[b][a] =
            close(l(b) - l(a)) ? detail::expDividedDifferenceSeries(
                                     std::array<double, 2>{l(a), l(b)})
                               : (e[b] - e[a]) / (l(b) - l(a));
    for (int a = 0; a != 3; ++a)
      for (int b = a; b != 3; ++b)
        for (int c = b; c != 3; ++c)
          third_[a][b][c] =
              close(l(c) - l(a))
                  ? detail::expDividedDifferenceSeries(
                        std::array<double, 3>{l(a), l(b), l(c)})
                  : (second_[b][c] - second_[a][b]) / (l(c) - l(a));
    for (int a = 0; a != 3; ++a)
      for (int b = 0; b != 3; ++b)
        for (int c = 0; c != 3; ++c) {
          std::array<int, 3> s{a, b, c};
          std::sort(s.begin(), s.end());
          third_[a][b][c] = third_[s[0]][s[1]][s[2]];
        }
    for (int a = 0; a != 3; ++a)
      for (int b = a; b != 3; ++b)
        for (int c = b; c != 3; ++c)
          for (int d = c; d != 3; ++d)
            fourth_[a][b][c][d] =
                close(l(d) - l(a))
                    ? detail::expDividedDifferenceSeries(
                          std::array<double, 4>{l(a), l(b), l(c), l(d)})
                    : (third_[b][c][d] - third_[a][b][c]) / (l(d) - l(a));
    for (int a = 0; a != 3; ++a)
      for (int b = 0; b != 3; ++b)
        for (int c = 0; c != 3; ++c)
          for (int d = 0; d != 3; ++d) {
            std::array<int, 4> s{a, b, c, d};
            std::sort(s.begin(), s.end());
            fourth_[a][b][c][d] = fourth_[s[0]][s[1]][s[2]][s[3]];
          }
  }
 
  /** exp[l_a,l_b], exp[l_a,l_b,l_c], exp[l_a,l_b,l_c,l_d]. */
  double second_[3][3];
  double third_[3][3][3];
  double fourth_[3][3][3][3];
};
 
/** Rotation exp(skew(w)) and its directional derivatives up to third order
 * from LieGroups::SO3 (Lie.hpp): exp, diffExp, diffDiffExp and
 * diffDiffDiffExp contracted with the axial direction vectors. The value and
 * the first two derivative tensors are formed once per base and reused for
 * every direction set. */
struct RotationExponential {
  explicit RotationExponential(const Vector &w) : omega(w) {
    FTENSOR_INDEXES(3, i, j, k, l);
    const double theta = std::sqrt(omega(i) * omega(i));
    R(i, j) = LieGroups::SO3::exp(omega, theta)(i, j);
    dR(i, j, k) = LieGroups::SO3::diffExp(omega, theta)(i, j, k);
    ddR(i, j, k, l) = LieGroups::SO3::diffDiffExp(omega, theta)(i, j, k, l);
  }
 
  Tensor value() const { return R; }
 
  Tensor first(const Vector &a) const {
    FTENSOR_INDEXES(3, i, j, k);
    Tensor r;
    r(i, j) = dR(i, j, k) * a(k);
    return r;
  }
 
  Tensor second(const Vector &a, const Vector &b) const {
    FTENSOR_INDEXES(3, i, j, k, l);
    Tensor r;
    r(i, j) = ddR(i, j, k, l) * a(k) * b(l);
    return r;
  }
 
  Tensor third(const Vector &a, const Vector &b, const Vector &c) const {
    FTENSOR_INDEXES(3, i, j, k, l);
    const double theta = std::sqrt(omega(i) * omega(i));
    const auto dddR = LieGroups::SO3::diffDiffDiffExp(omega, theta, c);
    Tensor r;
    r(i, j) = dddR(i, j, k, l) * a(k) * b(l);
    return r;
  }
 
  /** Fill the jet of exp(skew(w)): block 0 the value and the other blocks
   * the derivatives in the axial vectors of the skew direction blocks of
   * `argument`. A zero direction block gives zero derivatives. */
  template <int N>
  void fill(const Matrix<N> &argument, Matrix<N> &result) const {
    FTENSOR_INDEXES(3, i, j);
    result.data.fill(0.);
    result[0](i, j) = R(i, j);
    std::array<Vector, N> v;
    std::array<bool, N> active{};
    for (int d = 0; d != N; ++d) {
      v[d] = axial(argument[1 << d]);
      active[d] = v[d](i) * v[d](i) > 0.;
      if (active[d])
        result[1 << d](i, j) = first(v[d])(i, j);
    }
    for (int d1 = 0; d1 != N; ++d1)
      for (int d2 = d1 + 1; d2 != N; ++d2)
        if (active[d1] && active[d2])
          result[(1 << d1) | (1 << d2)](i, j) = second(v[d1], v[d2])(i, j);
    if constexpr (N == 3)
      if (active[0] && active[1] && active[2])
        result[7](i, j) = third(v[0], v[1], v[2])(i, j);
  }
 
  Vector omega;
  Tensor R;
  FTensor::Tensor3<double, 3, 3, 3> dR;
  FTensor::Tensor4<double, 3, 3, 3, 3> ddR;
};
 
/** Product interpolation exp(skew(a)) exp(skew(q)).
 * Used by L2 rotation kinematics and reduced curvature. Constant and higher
 * directions stay distinct; the jet product supplies every mixed derivative.
 */
struct SplitRotation {
  SplitRotation(const Vector &a, const Vector &q) : constant(a), higher(q) {}
 
  Tensor value() const {
    FTENSOR_INDEXES(3, i, j, k);
    Tensor result;
    result(i, j) = constant.R(i, k) * higher.R(k, j);
    return result;
  }
 
  Tensor first(const Vector &da, const Vector &dq) const {
    FTENSOR_INDEXES(3, i, j, k);
    const auto dA = constant.first(da), dQ = higher.first(dq);
    Tensor result;
    result(i, j) = dA(i, k) * higher.R(k, j) + constant.R(i, k) * dQ(k, j);
    return result;
  }
 
  /** Coordinate derivative tensor for one of the two three-vector factors. */
  auto firstTensor(const int mode) const {
    FTENSOR_INDEXES(3, i, j, k, m);
    FTensor::Tensor3<double, 3, 3, 3> result;
    if (!mode)
      result(i, j, k) = constant.dR(i, m, k) * higher.R(m, j);
    else
      result(i, j, k) = constant.R(i, m) * higher.dR(m, j, k);
    return result;
  }
 
  /** Constant/constant, constant/higher, higher/constant or higher/higher. */
  auto secondTensor(const int row_mode, const int col_mode) const {
    FTENSOR_INDEXES(3, i, j, k, l, m);
    FTensor::Tensor4<double, 3, 3, 3, 3> result;
    if (!row_mode && !col_mode)
      result(i, j, k, l) = constant.ddR(i, m, k, l) * higher.R(m, j);
    else if (!row_mode)
      result(i, j, k, l) = constant.dR(i, m, k) * higher.dR(m, j, l);
    else if (!col_mode)
      result(i, j, k, l) = constant.dR(i, m, l) * higher.dR(m, j, k);
    else
      result(i, j, k, l) = constant.R(i, m) * higher.ddR(m, j, k, l);
    return result;
  }
 
  template <int N>
  void fill(const Matrix<N> &a, const Matrix<N> &q, Matrix<N> &result) const {
    Matrix<N> A, Q;
    constant.fill(a, A);
    higher.fill(q, Q);
    result = product(A, Q);
  }
 
  RotationExponential constant, higher;
};
 
/** Deviatoric curvature of f = C10 tr exp(2 dev H) as symmetric tensors to be
 * contracted with a stretch direction e: second(i,j) e(i,j) = D^2 f[s, e] and
 * third(i,j) e(i,j) = D^3 f[p, s, e]. */
struct DeviatorCurvature {
  Symmetric second, third;
};
 
inline DeviatorCurvature deviatorCurvature(double c10, const Symmetric &h,
                                           const Symmetric &p,
                                           const Symmetric &s) {
  FTENSOR_INDEXES(3, i, j, k, l);
  constexpr auto t_diff_deviator = FTensor::diff_deviator<int, 3>();
  const auto deviator = [&](const Symmetric &t) {
    Symmetric result;
    result(i, j) = 2. * (t_diff_deviator(i, j, k, l) * t(k, l));
    return result;
  };
  const auto a = deviator(h), dp = deviator(p), ds = deviator(s);
  const SymmetricExponential exponential(a);
  // D^k tr exp(A)[..., E] = D^(k-1) exp(A)[...] : E with E = 2 dev e, so the
  // tensors to contract with e are twice the deviators of the derivatives.
  DeviatorCurvature result;
  result.second(i, j) = c10 * deviator(exponential.first(ds))(i, j);
  result.third(i, j) = c10 * deviator(exponential.second(dp, ds))(i, j);
  return result;
}
} // namespace LogarithmicDirectionalDerivatives
} // namespace EshelbianPlasticity
 
#endif // ESHELBIAN_LOGARITHMIC_DIRECTIONAL_DERIVATIVES_HPP