xref: /libCEED/examples/solids/qfunctions/finite-strain-mooney-rivlin.h (revision d4cc18453651bd0f94c1a2e078b2646a92dafdcc)

// Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors.
// All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
//
// SPDX-License-Identifier: BSD-2-Clause
//
// This file is part of CEED:  http://github.com/ceed

/// @file
/// Hyperelasticity, finite strain for solid mechanics example using PETSc

#include <ceed/types.h>
#ifndef CEED_RUNNING_JIT_PASS
#include <math.h>
#endif

// -----------------------------------------------------------------------------
// Mooney-Rivlin context
#ifndef PHYSICS_STRUCT_MR
#define PHYSICS_STRUCT_MR
typedef struct Physics_private_MR *Physics_MR;

struct Physics_private_MR {
  // material properties for MR
  CeedScalar mu_1;
  CeedScalar mu_2;
  CeedScalar lambda;
};
#endif

// -----------------------------------------------------------------------------
// Series approximation of log1p()
//  log1p() is not vectorized in libc
//
//  The series expansion is accurate to 1e-7 in the range sqrt(2)/2 < J < sqrt(2), with machine precision accuracy near J=1.
//  The initialization extends this range to 0.35 ~= sqrt(2)/4 < J < sqrt(2)*2 ~= 2.83, which should be sufficient for applications of the Neo-Hookean
//  model.
// -----------------------------------------------------------------------------
#ifndef LOG1P_SERIES_SHIFTED
#define LOG1P_SERIES_SHIFTED
CEED_QFUNCTION_HELPER CeedScalar log1p_series_shifted(CeedScalar x) {
  const CeedScalar left = sqrt(2.) / 2 - 1, right = sqrt(2.) - 1;
  CeedScalar       sum = 0;
  if (1) {           // Disable if the smaller range sqrt(2)/2 < J < sqrt(2) is sufficient
    if (x < left) {  // Replace if with while for arbitrary range (may hurt vectorization)
      sum -= log(2.) / 2;
      x = 1 + 2 * x;
    } else if (right < x) {
      sum += log(2.) / 2;
      x = (x - 1) / 2;
    }
  }
  CeedScalar       y  = x / (2. + x);
  const CeedScalar y2 = y * y;
  sum += y;
  y *= y2;
  sum += y / 3;
  y *= y2;
  sum += y / 5;
  y *= y2;
  sum += y / 7;
  return 2 * sum;
};
#endif

// -----------------------------------------------------------------------------
// Compute det F - 1
// -----------------------------------------------------------------------------
#ifndef DETJM1
#define DETJM1
CEED_QFUNCTION_HELPER CeedScalar computeJM1(const CeedScalar grad_u[3][3]) {
  return grad_u[0][0] * (grad_u[1][1] * grad_u[2][2] - grad_u[1][2] * grad_u[2][1]) +
         grad_u[0][1] * (grad_u[1][2] * grad_u[2][0] - grad_u[1][0] * grad_u[2][2]) +
         grad_u[0][2] * (grad_u[1][0] * grad_u[2][1] - grad_u[2][0] * grad_u[1][1]) + grad_u[0][0] + grad_u[1][1] + grad_u[2][2] +
         grad_u[0][0] * grad_u[1][1] + grad_u[0][0] * grad_u[2][2] + grad_u[1][1] * grad_u[2][2] - grad_u[0][1] * grad_u[1][0] -
         grad_u[0][2] * grad_u[2][0] - grad_u[1][2] * grad_u[2][1];
};
#endif

// -----------------------------------------------------------------------------
// Compute matrix^(-1), where matrix is symetric, returns array of 6
// -----------------------------------------------------------------------------
#ifndef MatinvSym
#define MatinvSym
CEED_QFUNCTION_HELPER int computeMatinvSym(const CeedScalar A[3][3], const CeedScalar detA, CeedScalar Ainv[6]) {
  // Compute A^(-1) : A-Inverse
  CeedScalar B[6] = {
      A[1][1] * A[2][2] - A[1][2] * A[2][1], /* *NOPAD* */
      A[0][0] * A[2][2] - A[0][2] * A[2][0], /* *NOPAD* */
      A[0][0] * A[1][1] - A[0][1] * A[1][0], /* *NOPAD* */
      A[0][2] * A[1][0] - A[0][0] * A[1][2], /* *NOPAD* */
      A[0][1] * A[1][2] - A[0][2] * A[1][1], /* *NOPAD* */
      A[0][2] * A[2][1] - A[0][1] * A[2][2]  /* *NOPAD* */
  };
  for (CeedInt m = 0; m < 6; m++) Ainv[m] = B[m] / (detA);

  return CEED_ERROR_SUCCESS;
}
#endif
// -----------------------------------------------------------------------------
// Common computations between FS and dFS
// -----------------------------------------------------------------------------
CEED_QFUNCTION_HELPER int commonFSMR(const CeedScalar mu_1, const CeedScalar mu_2, const CeedScalar lambda, const CeedScalar grad_u[3][3],
                                     CeedScalar Swork[6], CeedScalar Cwork[6], CeedScalar Cinvwork[6], CeedScalar *logJ) {
  // E - Green-Lagrange strain tensor
  //     E = 1/2 (grad_u + grad_u^T + grad_u^T*grad_u)
  const CeedInt indj[6] = {0, 1, 2, 1, 0, 0}, indk[6] = {0, 1, 2, 2, 2, 1};
  CeedScalar    E2work[6];
  for (CeedInt m = 0; m < 6; m++) {
    E2work[m] = grad_u[indj[m]][indk[m]] + grad_u[indk[m]][indj[m]];
    for (CeedInt n = 0; n < 3; n++) E2work[m] += grad_u[n][indj[m]] * grad_u[n][indk[m]];
  }
  CeedScalar E2[3][3] = {
      {E2work[0], E2work[5], E2work[4]},
      {E2work[5], E2work[1], E2work[3]},
      {E2work[4], E2work[3], E2work[2]}
  };

  // C : right Cauchy-Green tensor
  // C = I + 2E
  const CeedScalar C[3][3] = {
      {1 + E2[0][0], E2[0][1],     E2[0][2]    },
      {E2[0][1],     1 + E2[1][1], E2[1][2]    },
      {E2[0][2],     E2[1][2],     1 + E2[2][2]}
  };

  Cwork[0] = C[0][0];
  Cwork[1] = C[1][1];
  Cwork[2] = C[2][2];
  Cwork[3] = C[1][2];
  Cwork[4] = C[0][2];
  Cwork[5] = C[0][1];
  // Compute invariants
  // I_1 = trace(C)
  const CeedScalar I_1 = C[0][0] + C[1][1] + C[2][2];
  // J-1
  const CeedScalar Jm1 = computeJM1(grad_u);
  // J = Jm1 + 1
  // Compute C^(-1) : C-Inverse
  const CeedScalar detC = (Jm1 + 1.) * (Jm1 + 1.);
  computeMatinvSym(C, detC, Cinvwork);

  // Compute the Second Piola-Kirchhoff (S)
  // S = (lambda*logJ - mu_1 -2*mu_2)*Cinvwork +(mu_1+mu_2*I_1)*I3-mu_2*Cwork
  // *1 for indices 0-2 for I_3

  *logJ = log1p_series_shifted(Jm1);
  for (CeedInt i = 0; i < 6; i++) {
    Swork[i] = (lambda * *logJ - mu_1 - 2 * mu_2) * Cinvwork[i] + (mu_1 + mu_2 * I_1) * (i < 3)  // identity I_3
               - mu_2 * Cwork[i];
  }

  return CEED_ERROR_SUCCESS;
}

// -----------------------------------------------------------------------------
// Residual evaluation for hyperelasticity, finite strain
// -----------------------------------------------------------------------------
CEED_QFUNCTION(ElasFSResidual_MR)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
  // Inputs
  const CeedScalar(*ug)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0], (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1];

  // Outputs
  CeedScalar(*dvdX)[3][CEED_Q_VLA] = (CeedScalar(*)[3][CEED_Q_VLA])out[0];
  // Store grad_u for HyperFSdF (Jacobian of HyperFSF)
  CeedScalar(*grad_u)[3][CEED_Q_VLA] = (CeedScalar(*)[3][CEED_Q_VLA])out[1];

  // Context
  const Physics_MR context = (Physics_MR)ctx;
  const CeedScalar mu_1    = context->mu_1;
  const CeedScalar mu_2    = context->mu_2;
  const CeedScalar lambda  = context->lambda;

  // Formulation Terminology:
  //  I3    : 3x3 Identity matrix
  //  C     : right Cauchy-Green tensor
  //  C_inv : inverse of C
  //  F     : deformation gradient
  //  S     : 2nd Piola-Kirchhoff
  //  P     : 1st Piola-Kirchhoff

  // Quadrature Point Loop
  CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
    // Read spatial derivatives of u
    const CeedScalar du[3][3] = {
        {ug[0][0][i], ug[1][0][i], ug[2][0][i]},
        {ug[0][1][i], ug[1][1][i], ug[2][1][i]},
        {ug[0][2][i], ug[1][2][i], ug[2][2][i]}
    };
    // -- Qdata
    const CeedScalar wdetJ      = q_data[0][i];
    const CeedScalar dXdx[3][3] = {
        {q_data[1][i], q_data[2][i], q_data[3][i]},
        {q_data[4][i], q_data[5][i], q_data[6][i]},
        {q_data[7][i], q_data[8][i], q_data[9][i]}
    };

    // Compute grad_u
    //   dXdx = (dx/dX)^(-1)
    // Apply dXdx to du = grad_u
    for (CeedInt j = 0; j < 3; j++) {    // Component
      for (CeedInt k = 0; k < 3; k++) {  // Derivative
        grad_u[j][k][i] = 0;
        for (CeedInt m = 0; m < 3; m++) grad_u[j][k][i] += dXdx[m][k] * du[j][m];
      }
    }

    // I3 : 3x3 Identity matrix
    // Compute The Deformation Gradient : F = I3 + grad_u
    const CeedScalar F[3][3] = {
        {grad_u[0][0][i] + 1, grad_u[0][1][i],     grad_u[0][2][i]    },
        {grad_u[1][0][i],     grad_u[1][1][i] + 1, grad_u[1][2][i]    },
        {grad_u[2][0][i],     grad_u[2][1][i],     grad_u[2][2][i] + 1}
    };

    const CeedScalar tempgradu[3][3] = {
        {grad_u[0][0][i], grad_u[0][1][i], grad_u[0][2][i]},
        {grad_u[1][0][i], grad_u[1][1][i], grad_u[1][2][i]},
        {grad_u[2][0][i], grad_u[2][1][i], grad_u[2][2][i]}
    };

    // Common components of finite strain calculations
    CeedScalar Swork[6], Cwork[6], Cinvwork[6], logJ;
    commonFSMR(mu_1, mu_2, lambda, tempgradu, Swork, Cwork, Cinvwork, &logJ);

    // Second Piola-Kirchhoff (S)
    const CeedScalar S[3][3] = {
        {Swork[0], Swork[5], Swork[4]},
        {Swork[5], Swork[1], Swork[3]},
        {Swork[4], Swork[3], Swork[2]}
    };

    // Compute the First Piola-Kirchhoff : P = F*S
    CeedScalar P[3][3];
    for (CeedInt j = 0; j < 3; j++) {
      for (CeedInt k = 0; k < 3; k++) {
        P[j][k] = 0;
        for (CeedInt m = 0; m < 3; m++) P[j][k] += F[j][m] * S[m][k];
      }
    }

    // Apply dXdx^T and weight to P (First Piola-Kirchhoff)
    for (CeedInt j = 0; j < 3; j++) {    // Component
      for (CeedInt k = 0; k < 3; k++) {  // Derivative
        dvdX[k][j][i] = 0;
        for (CeedInt m = 0; m < 3; m++) dvdX[k][j][i] += dXdx[k][m] * P[j][m] * wdetJ;
      }
    }
  }  // End of Quadrature Point Loop

  return CEED_ERROR_SUCCESS;
}

// -----------------------------------------------------------------------------
// Jacobian evaluation for hyperelasticity, finite strain
// -----------------------------------------------------------------------------
CEED_QFUNCTION(ElasFSJacobian_MR)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
  // Inputs
  const CeedScalar(*deltaug)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0],
        (*q_data)[CEED_Q_VLA]               = (const CeedScalar(*)[CEED_Q_VLA])in[1];
  // grad_u is used for hyperelasticity (non-linear)
  const CeedScalar(*grad_u)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[2];

  // Outputs
  CeedScalar(*deltadvdX)[3][CEED_Q_VLA] = (CeedScalar(*)[3][CEED_Q_VLA])out[0];

  // Context
  const Physics_MR context = (Physics_MR)ctx;
  const CeedScalar mu_1    = context->mu_1;
  const CeedScalar mu_2    = context->mu_2;
  const CeedScalar lambda  = context->lambda;

  // Quadrature Point Loop
  CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
    // Read spatial derivatives of delta_u
    const CeedScalar deltadu[3][3] = {
        {deltaug[0][0][i], deltaug[1][0][i], deltaug[2][0][i]},
        {deltaug[0][1][i], deltaug[1][1][i], deltaug[2][1][i]},
        {deltaug[0][2][i], deltaug[1][2][i], deltaug[2][2][i]}
    };
    // -- Qdata
    const CeedScalar wdetJ      = q_data[0][i];
    const CeedScalar dXdx[3][3] = {
        {q_data[1][i], q_data[2][i], q_data[3][i]},
        {q_data[4][i], q_data[5][i], q_data[6][i]},
        {q_data[7][i], q_data[8][i], q_data[9][i]}
    };

    // Compute graddeltau
    //   dXdx = (dx/dX)^(-1)
    // Apply dXdx to deltadu = graddelta
    // this is dF
    CeedScalar graddeltau[3][3];
    for (CeedInt j = 0; j < 3; j++) {    // Component
      for (CeedInt k = 0; k < 3; k++) {  // Derivative
        graddeltau[j][k] = 0;
        for (CeedInt m = 0; m < 3; m++) graddeltau[j][k] += dXdx[m][k] * deltadu[j][m];
      }
    }

    // I3 : 3x3 Identity matrix
    // Deformation Gradient : F = I3 + grad_u
    const CeedScalar F[3][3] = {
        {grad_u[0][0][i] + 1, grad_u[0][1][i],     grad_u[0][2][i]    },
        {grad_u[1][0][i],     grad_u[1][1][i] + 1, grad_u[1][2][i]    },
        {grad_u[2][0][i],     grad_u[2][1][i],     grad_u[2][2][i] + 1}
    };

    const CeedScalar tempgradu[3][3] = {
        {grad_u[0][0][i], grad_u[0][1][i], grad_u[0][2][i]},
        {grad_u[1][0][i], grad_u[1][1][i], grad_u[1][2][i]},
        {grad_u[2][0][i], grad_u[2][1][i], grad_u[2][2][i]}
    };

    // Common components of finite strain calculations
    CeedScalar Swork[6], Cwork[6], Cinvwork[6], logJ;
    commonFSMR(mu_1, mu_2, lambda, tempgradu, Swork, Cwork, Cinvwork, &logJ);

    // dE - Green-Lagrange strain tensor
    const CeedInt indj[6] = {0, 1, 2, 1, 0, 0}, indk[6] = {0, 1, 2, 2, 2, 1};
    CeedScalar    dEwork[6];
    for (CeedInt m = 0; m < 6; m++) {
      dEwork[m] = 0;
      for (CeedInt n = 0; n < 3; n++) dEwork[m] += (graddeltau[n][indj[m]] * F[n][indk[m]] + F[n][indj[m]] * graddeltau[n][indk[m]]) / 2.;
    }
    CeedScalar dE[3][3] = {
        {dEwork[0], dEwork[5], dEwork[4]},
        {dEwork[5], dEwork[1], dEwork[3]},
        {dEwork[4], dEwork[3], dEwork[2]}
    };
    // C : right Cauchy-Green tensor
    // C^(-1) : C-Inverse
    const CeedScalar C[3][3] = {
        {Cwork[0], Cwork[5], Cwork[4]},
        {Cwork[5], Cwork[1], Cwork[3]},
        {Cwork[4], Cwork[3], Cwork[2]}
    };
    const CeedScalar C_inv[3][3] = {
        {Cinvwork[0], Cinvwork[5], Cinvwork[4]},
        {Cinvwork[5], Cinvwork[1], Cinvwork[3]},
        {Cinvwork[4], Cinvwork[3], Cinvwork[2]}
    };
    // -- C_inv:dE
    CeedScalar Cinv_contract_dE = 0;
    for (CeedInt j = 0; j < 3; j++) {
      for (CeedInt k = 0; k < 3; k++) Cinv_contract_dE += C_inv[j][k] * dE[j][k];
    }

    // -- C:dE
    CeedScalar C_contract_dE = 0;
    for (CeedInt j = 0; j < 3; j++) {
      for (CeedInt k = 0; k < 3; k++) C_contract_dE += C[j][k] * dE[j][k];
    }

    // -- dE*C_inv
    CeedScalar dE_Cinv[3][3];
    for (CeedInt j = 0; j < 3; j++) {
      for (CeedInt k = 0; k < 3; k++) {
        dE_Cinv[j][k] = 0;
        for (CeedInt m = 0; m < 3; m++) dE_Cinv[j][k] += dE[j][m] * C_inv[m][k];
      }
    }

    // -- C_inv*dE*C_inv
    CeedScalar Cinv_dE_Cinv[3][3];
    // This product is symmetric and we only use the upper-triangular part below, but naively compute the whole thing here
    for (CeedInt j = 0; j < 3; j++) {
      for (CeedInt k = 0; k < 3; k++) {
        Cinv_dE_Cinv[j][k] = 0;
        for (CeedInt m = 0; m < 3; m++) Cinv_dE_Cinv[j][k] += C_inv[j][m] * dE_Cinv[m][k];
      }
    }

    // Compute dS = (mu_2)*((2*I_3:dE)*I_3 - dE) + 2*d*Cinv_dE_Cinv + lambda*Cinv_contract_dE*Cinvwork - 2*lambda*logJ*Cinv_dE_Cinv;
    // (2*I_3:dE)*I_3 - dE = 2*trace(dE)*I_3 - dE = 2trace(dE) - dE on the diagonal
    // (2*I_3:dE)*I_3 - dE = -dE elsewhere
    // CeedScalar J = Jm1 + 1;
    CeedScalar tr_dE = dE[0][0] + dE[1][1] + dE[2][2];
    CeedScalar dSwork[6];
    for (CeedInt i = 0; i < 6; i++) {
      dSwork[i] = lambda * Cinv_contract_dE * Cinvwork[i] + 2 * (mu_1 + 2 * mu_2 - lambda * logJ) * Cinv_dE_Cinv[indj[i]][indk[i]] +
                  2 * mu_2 * (tr_dE * (i < 3) - dEwork[i]);
    }

    CeedScalar dS[3][3] = {
        {dSwork[0], dSwork[5], dSwork[4]},
        {dSwork[5], dSwork[1], dSwork[3]},
        {dSwork[4], dSwork[3], dSwork[2]}
    };
    // Second Piola-Kirchhoff (S)
    const CeedScalar S[3][3] = {
        {Swork[0], Swork[5], Swork[4]},
        {Swork[5], Swork[1], Swork[3]},
        {Swork[4], Swork[3], Swork[2]}
    };
    // dP = dPdF:dF = dF*S + F*dS
    CeedScalar dP[3][3];
    for (CeedInt j = 0; j < 3; j++) {
      for (CeedInt k = 0; k < 3; k++) {
        dP[j][k] = 0;
        for (CeedInt m = 0; m < 3; m++) dP[j][k] += graddeltau[j][m] * S[m][k] + F[j][m] * dS[m][k];
      }
    }

    // Apply dXdx^T and weight
    for (CeedInt j = 0; j < 3; j++) {    // Component
      for (CeedInt k = 0; k < 3; k++) {  // Derivative
        deltadvdX[k][j][i] = 0;
        for (CeedInt m = 0; m < 3; m++) deltadvdX[k][j][i] += dXdx[k][m] * dP[j][m] * wdetJ;
      }
    }
  }  // End of Quadrature Point Loop

  return CEED_ERROR_SUCCESS;
}

// -----------------------------------------------------------------------------
// Strain energy computation for hyperelasticity, finite strain
// -----------------------------------------------------------------------------
CEED_QFUNCTION(ElasFSEnergy_MR)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
  // Inputs
  const CeedScalar(*ug)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0], (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1];

  // Outputs
  CeedScalar(*energy) = (CeedScalar(*))out[0];

  // Context
  const Physics_MR context = (Physics_MR)ctx;
  const CeedScalar mu_1    = context->mu_1;
  const CeedScalar mu_2    = context->mu_2;
  const CeedScalar lambda  = context->lambda;

  // Quadrature Point Loop
  CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
    // Read spatial derivatives of u
    const CeedScalar du[3][3] = {
        {ug[0][0][i], ug[1][0][i], ug[2][0][i]},
        {ug[0][1][i], ug[1][1][i], ug[2][1][i]},
        {ug[0][2][i], ug[1][2][i], ug[2][2][i]}
    };
    // -- Qdata
    const CeedScalar wdetJ      = q_data[0][i];
    const CeedScalar dXdx[3][3] = {
        {q_data[1][i], q_data[2][i], q_data[3][i]},
        {q_data[4][i], q_data[5][i], q_data[6][i]},
        {q_data[7][i], q_data[8][i], q_data[9][i]}
    };
    // Compute grad_u
    //   dXdx = (dx/dX)^(-1)
    // Apply dXdx to du = grad_u
    CeedScalar grad_u[3][3];
    for (int j = 0; j < 3; j++) {    // Component
      for (int k = 0; k < 3; k++) {  // Derivative
        grad_u[j][k] = 0;
        for (int m = 0; m < 3; m++) grad_u[j][k] += dXdx[m][k] * du[j][m];
      }
    }

    // E - Green-Lagrange strain tensor
    // E = 1/2 (grad_u + grad_u^T + grad_u^T*grad_u)
    const CeedInt indj[6] = {0, 1, 2, 1, 0, 0}, indk[6] = {0, 1, 2, 2, 2, 1};
    CeedScalar    E2work[6];
    for (CeedInt m = 0; m < 6; m++) {
      E2work[m] = grad_u[indj[m]][indk[m]] + grad_u[indk[m]][indj[m]];
      for (CeedInt n = 0; n < 3; n++) E2work[m] += grad_u[n][indj[m]] * grad_u[n][indk[m]];
    }
    CeedScalar E2[3][3] = {
        {E2work[0], E2work[5], E2work[4]},
        {E2work[5], E2work[1], E2work[3]},
        {E2work[4], E2work[3], E2work[2]}
    };

    // C : right Cauchy-Green tensor
    // C = I + 2E
    const CeedScalar C[3][3] = {
        {1 + E2[0][0], E2[0][1],     E2[0][2]    },
        {E2[0][1],     1 + E2[1][1], E2[1][2]    },
        {E2[0][2],     E2[1][2],     1 + E2[2][2]}
    };
    // Compute CC = C*C = C^2
    CeedScalar CC[3][3];
    for (CeedInt j = 0; j < 3; j++) {
      for (CeedInt k = 0; k < 3; k++) {
        CC[j][k] = 0;
        for (CeedInt m = 0; m < 3; m++) CC[j][k] += C[j][m] * C[m][k];
      }
    }

    const CeedScalar Jm1 = computeJM1(grad_u);
    // CeedScalar J = Jm1 + 1;
    // Compute invariants
    // I_1 = trace(C)
    const CeedScalar I_1 = C[0][0] + C[1][1] + C[2][2];
    // Trace(C^2)
    const CeedScalar tr_CC = CC[0][0] + CC[1][1] + CC[2][2];
    // I_2 = 0.5(I_1^2 - trace(C^2))
    const CeedScalar I_2 = 0.5 * (I_1 * I_1 - tr_CC);

    const CeedScalar logJ = log1p_series_shifted(Jm1);
    // Strain energy Phi(E) for Mooney-Rivlin
    energy[i] = (0.5 * lambda * (logJ) * (logJ) - (mu_1 + 2 * mu_2) * logJ + (mu_1 / 2.) * (I_1 - 3) + (mu_2 / 2.) * (I_2 - 3)) * wdetJ;

  }  // End of Quadrature Point Loop

  return CEED_ERROR_SUCCESS;
}

// -----------------------------------------------------------------------------
// Nodal diagnostic quantities for hyperelasticity, finite strain
// -----------------------------------------------------------------------------
CEED_QFUNCTION(ElasFSDiagnostic_MR)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
  // Inputs
  const CeedScalar(*u)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], (*ug)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[1],
        (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2];

  // Outputs
  CeedScalar(*diagnostic)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];

  // Context
  const Physics_MR context = (Physics_MR)ctx;
  const CeedScalar mu_1    = context->mu_1;
  const CeedScalar mu_2    = context->mu_2;
  const CeedScalar lambda  = context->lambda;

  // Quadrature Point Loop
  CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
    // Read spatial derivatives of u
    const CeedScalar du[3][3] = {
        {ug[0][0][i], ug[1][0][i], ug[2][0][i]},
        {ug[0][1][i], ug[1][1][i], ug[2][1][i]},
        {ug[0][2][i], ug[1][2][i], ug[2][2][i]}
    };
    // -- Qdata
    const CeedScalar dXdx[3][3] = {
        {q_data[1][i], q_data[2][i], q_data[3][i]},
        {q_data[4][i], q_data[5][i], q_data[6][i]},
        {q_data[7][i], q_data[8][i], q_data[9][i]}
    };

    // Compute grad_u
    //   dXdx = (dx/dX)^(-1)
    // Apply dXdx to du = grad_u
    CeedScalar grad_u[3][3];
    for (int j = 0; j < 3; j++) {    // Component
      for (int k = 0; k < 3; k++) {  // Derivative
        grad_u[j][k] = 0;
        for (int m = 0; m < 3; m++) grad_u[j][k] += dXdx[m][k] * du[j][m];
      }
    }

    // E - Green-Lagrange strain tensor
    //     E = 1/2 (grad_u + grad_u^T + grad_u^T*grad_u)
    const CeedInt indj[6] = {0, 1, 2, 1, 0, 0}, indk[6] = {0, 1, 2, 2, 2, 1};
    CeedScalar    E2work[6];
    for (CeedInt m = 0; m < 6; m++) {
      E2work[m] = grad_u[indj[m]][indk[m]] + grad_u[indk[m]][indj[m]];
      for (CeedInt n = 0; n < 3; n++) E2work[m] += grad_u[n][indj[m]] * grad_u[n][indk[m]];
    }
    CeedScalar E2[3][3] = {
        {E2work[0], E2work[5], E2work[4]},
        {E2work[5], E2work[1], E2work[3]},
        {E2work[4], E2work[3], E2work[2]}
    };

    // Displacement
    diagnostic[0][i] = u[0][i];
    diagnostic[1][i] = u[1][i];
    diagnostic[2][i] = u[2][i];

    // Pressure
    const CeedScalar Jm1  = computeJM1(grad_u);
    const CeedScalar logJ = log1p_series_shifted(Jm1);
    diagnostic[3][i]      = -lambda * logJ;

    // Stress tensor invariants
    diagnostic[4][i] = (E2[0][0] + E2[1][1] + E2[2][2]) / 2.;
    diagnostic[5][i] = 0.;
    for (CeedInt j = 0; j < 3; j++) {
      for (CeedInt m = 0; m < 3; m++) diagnostic[5][i] += E2[j][m] * E2[m][j] / 4.;
    }
    diagnostic[6][i] = Jm1 + 1;

    // C : right Cauchy-Green tensor
    // C = I + 2E
    const CeedScalar C[3][3] = {
        {1 + E2[0][0], E2[0][1],     E2[0][2]    },
        {E2[0][1],     1 + E2[1][1], E2[1][2]    },
        {E2[0][2],     E2[1][2],     1 + E2[2][2]}
    };
    // Compute CC = C*C = C^2
    CeedScalar CC[3][3];
    for (CeedInt j = 0; j < 3; j++) {
      for (CeedInt k = 0; k < 3; k++) {
        CC[j][k] = 0;
        for (CeedInt m = 0; m < 3; m++) CC[j][k] += C[j][m] * C[m][k];
      }
    }

    // CeedScalar J = Jm1 + 1;
    // Compute invariants
    // I_1 = trace(C)
    const CeedScalar I_1 = C[0][0] + C[1][1] + C[2][2];
    // Trace(C^2)
    const CeedScalar tr_CC = CC[0][0] + CC[1][1] + CC[2][2];
    // I_2 = 0.5(I_1^2 - trace(C^2))
    const CeedScalar I_2 = 0.5 * (pow(I_1, 2) - tr_CC);

    // Strain energy
    diagnostic[7][i] = (0.5 * lambda * logJ * logJ - (mu_1 + 2 * mu_2) * logJ + (mu_1 / 2.) * (I_1 - 3) + (mu_2 / 2.) * (I_2 - 3));
  }  // End of Quadrature Point Loop

  return CEED_ERROR_SUCCESS;
}
// -----------------------------------------------------------------------------