rotation_split_atom.cpp

Go to the documentation of this file.

/** Regression checks for product L2 rotations, derivatives and TS branching. */
#include <array>
#include <numbers>
#include <MoFEM.hpp>
using namespace MoFEM;
#include "AuxiliaryLogarithmicStressTestSetup.hpp"
#include <L2Rotation.hpp>
#include <LogarithmicDirectionalDerivatives.hpp>
#include <TSElasticPreAndPostStep.hpp>
using namespace EshelbianPlasticity;
 
namespace {
 
struct OpCheckRotationLayout : VolUserDataOperator {
  OpCheckRotationLayout(const std::string &field,
                        boost::shared_ptr<PetscInt> visits)
      : VolUserDataOperator(field, OPROW), visitsPtr(std::move(visits)) {
    doEntities.fill(false);
    doEntities[MBTET] = true;
  }
 
  MoFEMErrorCode doWork(int, EntityType, EntData &data) override {
    MoFEMFunctionBegin;
    const int dofs = data.getFieldData().size();
    if (!dofs)
      MoFEMFunctionReturnHot(0);
    const int points = getGaussPts().size2();
    if (dofs < 3 || dofs % 3 || points <= 1 || data.getN().size2() < dofs / 3)
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Incomplete tetrahedral L2 rotation block");
    const auto &field_dofs = data.getFieldDofs();
    for (int c = 0; c != 3; ++c)
      if (field_dofs[c]->getDofOrder() != 0 ||
          field_dofs[c]->getDofCoeffIdx() != c)
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "First rotation block is not the constant vector mode");
    const double phi0 = data.getN()(0, 0);
    if (!std::isfinite(phi0) || std::abs(phi0) < 1.e-12)
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Constant rotation basis has invalid normalisation");
    FTENSOR_INDEX(3, i);
    for (int gg = 0; gg != points; ++gg) {
      auto gradient = data.getFTensor1DiffN<3>(gg, 0);
      if (std::abs(data.getN()(gg, 0) - phi0) > 1.e-13 ||
          gradient(i) * gradient(i) > 1.e-26)
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "First L2 rotation basis is not spatially constant");
    }
    ++*visitsPtr;
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<PetscInt> visitsPtr;
};
 
/** Compare a split vector with ordinary interpolation and known c0. */
struct OpCheckRotationValues : VolUserDataOperator {
  OpCheckRotationValues(MatrixPtr total, MatrixPtr constant, MatrixPtr higher,
                        MatrixPtr reference, double scale)
      : VolUserDataOperator(NOSPACE, OPSPACE), totalPtr(total),
        constantPtr(constant), higherPtr(higher), referencePtr(reference),
        stateScale(scale) {}
  MoFEMErrorCode doWork(int, EntityType, EntData &) override {
    MoFEMFunctionBegin;
    using View = GetFTensor1FromMatType<3, -1, DL>;
    const int points = getGaussPts().size2();
    auto t = MatrixSizeHelper<View, DL>::get(*totalPtr, points)();
    auto a = MatrixSizeHelper<View, DL>::get(*constantPtr, points)();
    auto q = MatrixSizeHelper<View, DL>::get(*higherPtr, points)();
    auto r = MatrixSizeHelper<View, DL>::get(*referencePtr, points)();
    FTENSOR_INDEX(3, i);
    const FTensor::Tensor1<double, 3> expected(
        .4 * stateScale, -.3 * stateScale, .2 * stateScale);
    for (int gg = 0; gg != points; ++gg) {
      FTensor::Tensor1<double, 3> error;
      error(i) = a(i) - expected(i);
      if (error(i) * error(i) > 1.e-24)
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "Incorrect constant rotation vector");
      error(i) = t(i) - r(i);
      if (error(i) * error(i) > 1.e-24)
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "Split total differs from core field evaluation");
      error(i) = t(i) - a(i) - q(i);
      if (error(i) * error(i) > 1.e-24)
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "Rotation decomposition does not reconstruct the total");
      ++t;
      ++a;
      ++q;
      ++r;
    }
    MoFEMFunctionReturn(0);
  }
  MatrixPtr totalPtr, constantPtr, higherPtr, referencePtr;
  double stateScale;
};
 
MoFEMErrorCode checkRotationVectors(EshelbianCore &ep) {
  MoFEMFunctionBegin;
  const auto make_vector = [&](double scale, SmartPetscObj<Vec> &vector) {
    MoFEMFunctionBegin;
    vector = createDMVector(ep.dmElastic);
    CHKERR VecZeroEntries(vector);
    PetscInt first, last;
    CHKERR VecGetOwnershipRange(vector, &first, &last);
    const Problem *problem;
    CHKERR DMMoFEMGetProblemPtr(ep.dmElastic, &problem);
    const double c0[] = {.4, -.3, .2};
    for (const auto &dof : *problem->getNumeredRowDofsPtr()) {
      const auto index = dof->getPetscGlobalDofIdx();
      if (dof->getName() != ep.rotAxis || index < first || index >= last)
        continue;
      const auto component = dof->getDofCoeffIdx();
      const double value =
          dof->getDofOrder() == 0 ? c0[component] : .001 * (index % 13 - 6);
      CHKERR VecSetValue(vector, index, scale * value, INSERT_VALUES);
    }
    CHKERR VecAssemblyBegin(vector);
    CHKERR VecAssemblyEnd(vector);
    CHKERR VecGhostUpdateBegin(vector, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(vector, INSERT_VALUES, SCATTER_FORWARD);
    MoFEMFunctionReturn(0);
  };
  SmartPetscObj<Vec> current, previous, direction;
  CHKERR make_vector(1., current);
  CHKERR make_vector(-.7, previous);
  CHKERR make_vector(.13, direction);
  CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, current, INSERT_VALUES,
                                 SCATTER_REVERSE);
  auto fe = boost::make_shared<VolumeElementForcesAndSourcesCore>(ep.mField);
  fe->getUserPolynomialBase() =
      boost::make_shared<CGGUserPolynomialBase>(nullptr, true);
  fe->getRuleHook = [](int, int, int order) { return 2 * order + 2; };
  CHKERR MoFEM::AddHOOps<3, 3, 3>::add(fe->getOpPtrVector(), {L2},
                                       ep.materialH1Positions);
  for (const auto &state :
       {std::make_pair(SmartPetscObj<Vec>(), 1.), std::make_pair(previous, -.7),
        std::make_pair(direction, .13)}) {
    auto total = boost::make_shared<MatrixDouble>();
    auto constant = boost::make_shared<MatrixDouble>();
    auto higher = boost::make_shared<MatrixDouble>();
    auto reference = boost::make_shared<MatrixDouble>();
    fe->getOpPtrVector().push_back(new OpCalculateL2RotationValues(
        ep.rotAxis, total, constant, higher, state.first));
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
        ep.rotAxis, reference, state.first, MBTET));
    fe->getOpPtrVector().push_back(new OpCheckRotationValues(
        total, constant, higher, reference, state.second));
  }
  CHKERR DMoFEMLoopFiniteElements(ep.dmElastic, ep.elementVolumeName, fe);
  MoFEMFunctionReturn(0);
}
 
/** Contract the same constant/higher Hessian blocks used by mechanics. */
LogarithmicDirectionalDerivatives::Tensor contractRotationHessian(
    const LogarithmicDirectionalDerivatives::SplitRotation &rotation,
    const LogarithmicDirectionalDerivatives::Vector &da,
    const LogarithmicDirectionalDerivatives::Vector &dq,
    const LogarithmicDirectionalDerivatives::Vector &ea,
    const LogarithmicDirectionalDerivatives::Vector &eq) {
  FTENSOR_INDEXES(3, i, j, k, l);
  const auto aa = rotation.secondTensor(0, 0);
  const auto aq = rotation.secondTensor(0, 1);
  const auto qa = rotation.secondTensor(1, 0);
  const auto qq = rotation.secondTensor(1, 1);
  LogarithmicDirectionalDerivatives::Tensor result;
  result(i, j) =
      aa(i, j, k, l) * da(k) * ea(l) + aq(i, j, k, l) * da(k) * eq(l) +
      qa(i, j, k, l) * dq(k) * ea(l) + qq(i, j, k, l) * dq(k) * eq(l);
  return result;
}
 
/** Independent central differences of values/derivatives, through order 3. */
MoFEMErrorCode checkSplitRotationDerivatives() {
  MoFEMFunctionBegin;
  using namespace LogarithmicDirectionalDerivatives;
  FTENSOR_INDEXES(3, i, j);
  struct Direction {
    Vector a, q;
  };
  const Vector zero(0., 0., 0.);
  const Direction directions[] = {{Vector(.3, -.2, .1), zero},
                                  {zero, Vector(-.1, .4, .2)},
                                  {Vector(.2, .1, -.3), Vector(-.3, .2, .1)}};
  const Direction bases[] = {
      {zero, zero},
      {Vector(.7, -.2, .4), Vector(-.3, .5, .1)},
      {Vector(0., 0., 6.1), Vector(.2, -.1, .08)},
      {Vector(1.e-8, -2.e-8, 3.e-8), Vector(.1, .2, -.1)}};
  const double step = 4.e-6;
  const auto check = [&](const auto &exact, const auto &fd,
                         const char *description) {
    MoFEMFunctionBegin;
    Tensor error;
    error(i, j) = exact(i, j) - fd(i, j);
    if (error(i, j) * error(i, j) > 1.e-14)
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Split rotation %s finite difference failed", description);
    MoFEMFunctionReturn(0);
  };
  for (const auto &base : bases) {
    const SplitRotation rotation(base.a, base.q);
    for (const auto &d : directions) {
      Vector plus_a, minus_a, plus_q, minus_q;
      plus_a(i) = base.a(i) + step * d.a(i);
      minus_a(i) = base.a(i) - step * d.a(i);
      plus_q(i) = base.q(i) + step * d.q(i);
      minus_q(i) = base.q(i) - step * d.q(i);
      const SplitRotation plus(plus_a, plus_q), minus(minus_a, minus_q);
      Tensor fd;
      fd(i, j) = (plus.value()(i, j) - minus.value()(i, j)) / (2 * step);
      CHKERR check(rotation.first(d.a, d.q), fd, "first");
      for (const auto &e : directions) {
        fd(i, j) = (plus.first(e.a, e.q)(i, j) - minus.first(e.a, e.q)(i, j)) /
                   (2 * step);
        CHKERR check(contractRotationHessian(rotation, d.a, d.q, e.a, e.q), fd,
                     "second");
        for (const auto &f : directions) {
          LogarithmicDirectionalDerivatives::Matrix<3> a, q, result;
          a[0](i, j) = skew(base.a)(i, j);
          q[0](i, j) = skew(base.q)(i, j);
          a[1](i, j) = skew(d.a)(i, j);
          q[1](i, j) = skew(d.q)(i, j);
          a[2](i, j) = skew(e.a)(i, j);
          q[2](i, j) = skew(e.q)(i, j);
          a[4](i, j) = skew(f.a)(i, j);
          q[4](i, j) = skew(f.q)(i, j);
          rotation.fill(a, q, result);
          fd(i, j) =
              (contractRotationHessian(plus, e.a, e.q, f.a, f.q)(i, j) -
               contractRotationHessian(minus, e.a, e.q, f.a, f.q)(i, j)) /
              (2 * step);
          CHKERR check(result[7], fd, "third");
          CHKERR check(result[3],
                       contractRotationHessian(rotation, d.a, d.q, e.a, e.q),
                       "jet second");
        }
      }
    }
  }
  const Vector a(0., 0., 6.1), q(.2, -.1, .08);
  const Vector shortened(0., 0., 6.1 - 2 * std::numbers::pi);
  CHKERR check(SplitRotation(a, q).value(), SplitRotation(shortened, q).value(),
               "branch value");
  MoFEMFunctionReturn(0);
}
 
struct OpCheckShortenedRotation : VolUserDataOperator {
  OpCheckShortenedRotation(std::array<MatrixPtr, 2> constant,
                           std::array<MatrixPtr, 2> higher)
      : VolUserDataOperator(NOSPACE, OPSPACE), constantPtr(constant),
        higherPtr(higher) {}
  MoFEMErrorCode doWork(int, EntityType, EntData &) override {
    MoFEMFunctionBegin;
    using namespace LogarithmicDirectionalDerivatives;
    FTENSOR_INDEXES(3, i, j);
    const int points = getGaussPts().size2();
    using View = GetFTensor1FromMatType<3, -1, DL>;
    auto a =
        std::array{MatrixSizeHelper<View, DL>::get(*constantPtr[0], points)(),
                   MatrixSizeHelper<View, DL>::get(*constantPtr[1], points)()};
    auto q =
        std::array{MatrixSizeHelper<View, DL>::get(*higherPtr[0], points)(),
                   MatrixSizeHelper<View, DL>::get(*higherPtr[1], points)()};
    for (int gg = 0; gg != points; ++gg) {
      Vector old_a, new_a, old_q, new_q;
      old_a(i) = a[0](i);
      new_a(i) = a[1](i);
      old_q(i) = q[0](i);
      new_q(i) = q[1](i);
      Tensor error;
      error(i, j) = SplitRotation(old_a, old_q).value()(i, j) -
                    SplitRotation(new_a, new_q).value()(i, j);
      if (error(i, j) * error(i, j) > 1.e-25 || new_a.l2() >= std::numbers::pi)
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "Constant branch update changed the product rotation");
      ++a[0];
      ++a[1];
      ++q[0];
      ++q[1];
    }
    MoFEMFunctionReturn(0);
  }
  std::array<MatrixPtr, 2> constantPtr, higherPtr;
};
 
// Verify that the Schur shift changes only preconditioning, including rebuilds.
MoFEMErrorCode checkRotationPreconditioner(EshelbianCore &ep, TS ts, Vec state,
                                           Vec rates, Mat jacobian,
                                           Mat pc_mat) {
  MoFEMFunctionBegin;
  if (ep.alphaOmegaPrec == 0.)
    MoFEMFunctionReturnHot(0);
  auto exact = vectorDuplicate(state);
  auto shifted = vectorDuplicate(state);
  auto residual = vectorDuplicate(state);
  CHKERR MatMult(jacobian, state, exact);
  CHKERR MatMult(pc_mat, state, shifted);
  CHKERR TSComputeIFunction(ts, 0., state, rates, residual, PETSC_FALSE);
  auto schur_shifted = matDuplicate(ep.S, MAT_COPY_VALUES);
 
  auto &pipeline = ep.elasticFeLhs->getOpPtrVector();
  if (!dynamic_cast<OpSpatialRotationPrec *>(&pipeline.back()))
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Expected rotation preconditioner at the end of the LHS pipeline");
  pipeline.pop_back();
  CHKERR TSComputeIJacobian(ts, 0., state, rates, 1., jacobian, pc_mat,
                            PETSC_FALSE);
  auto base = vectorDuplicate(state);
  CHKERR MatMult(jacobian, state, base);
  CHKERR VecAXPY(base, -1., exact);
  double error;
  CHKERR VecNorm(base, NORM_2, &error);
  if (error > 1e-10)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Rotation preconditioning changed the exact Jacobian: %g", error);
  CHKERR TSComputeIFunction(ts, 0., state, rates, base, PETSC_FALSE);
  CHKERR VecAXPY(base, -1., residual);
  CHKERR VecNorm(base, NORM_2, &error);
  if (error > 1e-12)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Rotation preconditioning changed the residual: %g", error);
 
  CHKERR MatMult(pc_mat, state, base);
  CHKERR VecAXPY(base, -1., shifted);
  double shift_norm, work;
  CHKERR VecNorm(base, NORM_2, &shift_norm);
  CHKERR VecDot(state, base, &work);
  if (!std::isfinite(shift_norm) || shift_norm < 1e-10 || work >= 0.)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Expected a nonzero positive rotation-gradient preconditioner");
  IS rotation_raw;
  CHKERR DMMoFEMGetFieldIS(ep.dmElastic, RowColData::ROW, ep.rotAxis.c_str(),
                           &rotation_raw);
  SmartPetscObj<IS> rotation_is(rotation_raw);
  Vec rotation_sub;
  CHKERR VecGetSubVector(base, rotation_is, &rotation_sub);
  CHKERR VecZeroEntries(rotation_sub);
  CHKERR VecRestoreSubVector(base, rotation_is, &rotation_sub);
  CHKERR VecNorm(base, NORM_2, &error);
  if (error > 1e-10)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Rotation preconditioner modified another field: %g", error);
  CHKERR MatAXPY(schur_shifted, -1., ep.S, SAME_NONZERO_PATTERN);
  CHKERR MatNorm(schur_shifted, NORM_FROBENIUS, &error);
  if (!std::isfinite(error) || error < 1e-10)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Schur preconditioner did not include the rotation shift");
 
  pipeline.push_back(
      new OpSpatialRotationPrec(ep.rotAxis, ep.dataAtPts, ep.alphaOmegaPrec));
  CHKERR TSComputeIJacobian(ts, 0., state, rates, 1., jacobian, pc_mat,
                            PETSC_FALSE);
  CHKERR MatMult(pc_mat, state, base);
  CHKERR VecAXPY(base, -1., shifted);
  CHKERR VecNorm(base, NORM_2, &error);
  if (error > 1e-10)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Rotation preconditioner changed after identical reassembly: %g",
            error);
 
  auto constant = vectorDuplicate(state);
  CHKERR VecZeroEntries(constant);
  const Problem *problem;
  CHKERR DMMoFEMGetProblemPtr(ep.dmElastic, &problem);
  PetscInt first, last;
  CHKERR VecGetOwnershipRange(constant, &first, &last);
  for (const auto &dof : *problem->getNumeredRowDofsPtr()) {
    const auto index = dof->getPetscGlobalDofIdx();
    if (dof->getName() == ep.rotAxis && dof->getDofOrder() == 0 &&
        index >= first && index < last)
      CHKERR VecSetValue(constant, index, 1., INSERT_VALUES);
  }
  CHKERR VecAssemblyBegin(constant);
  CHKERR VecAssemblyEnd(constant);
  CHKERR MatMult(pc_mat, constant, base);
  CHKERR MatMult(jacobian, constant, exact);
  CHKERR VecAXPY(base, -1., exact);
  CHKERR VecNorm(base, NORM_2, &error);
  if (error > 1e-10)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Gradient preconditioner penalised constant L2 rotation: %g",
            error);
 
  Mat a00;
  CHKERR MatNestGetSubMat(pc_mat, 1, 1, &a00);
  Vec x_raw, b_raw;
  CHKERR MatCreateVecs(a00, &x_raw, &b_raw);
  SmartPetscObj<Vec> x(x_raw), b(b_raw);
  auto solution = vectorDuplicate(x);
  CHKERR VecSet(x, 1.);
  CHKERR MatMult(a00, x, b);
  CHKERR MatSolve(a00, b, solution);
  double rhs_norm;
  CHKERR VecNorm(b, NORM_2, &rhs_norm);
  CHKERR MatMult(a00, solution, x);
  CHKERR VecAXPY(x, -1., b);
  CHKERR VecNorm(x, NORM_2, &error);
  error /= std::max(rhs_norm, std::numeric_limits<double>::min());
  if (!std::isfinite(error) || error > 1e-10)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Shifted A00 solve does not match its operator: %g", error);
  MOFEM_LOG("EP", Sev::inform)
      << "Rotation preconditioner: exact residual/Jacobian preserved, shift "
      << shift_norm << ", relative A00 residual " << error;
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkConstantBranchUpdate(EshelbianCore &ep) {
  MoFEMFunctionBegin;
  auto state = createDMVector(ep.dmElastic);
  CHKERR VecZeroEntries(state);
  const Problem *problem;
  CHKERR DMMoFEMGetProblemPtr(ep.dmElastic, &problem);
  PetscInt first, last;
  CHKERR VecGetOwnershipRange(state, &first, &last);
  for (const auto &dof : *problem->getNumeredRowDofsPtr()) {
    const auto index = dof->getPetscGlobalDofIdx();
    if (dof->getName() != ep.rotAxis || index < first || index >= last)
      continue;
    const double value = dof->getDofOrder() == 0
                             ? (dof->getDofCoeffIdx() == 2 ? 6.1 : 0.)
                             : .003 * (index % 11 - 5);
    CHKERR VecSetValue(state, index, value, INSERT_VALUES);
  }
  CHKERR VecAssemblyBegin(state);
  CHKERR VecAssemblyEnd(state);
  CHKERR VecGhostUpdateBegin(state, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(state, INSERT_VALUES, SCATTER_FORWARD);
  auto old = vectorDuplicate(state);
  CHKERR VecCopy(state, old);
  CHKERR VecCopy(state, ep.solTSStep);
  CHKERR VecScale(ep.solTSStep, .7);
  auto history = vectorDuplicate(ep.solTSStep);
  CHKERR VecCopy(ep.solTSStep, history);
  ep.dataAtPts = boost::make_shared<DataAtIntegrationPts>();
  auto ts = createTS(ep.mField.get_comm());
  CHKERR TSSetDM(ts, ep.dmElastic);
  CHKERR TSSetType(ts, TSBEULER);
  CHKERR TSSetSolution(ts, state);
  CHKERR TSElasticPreAndPostStep::setPreStep(&ep, ts);
  if (ep.gradApproximator != NO_H1_CONFIGURATION) {
    CHKERR TSPreStep(ts);
    PetscBool equal;
    CHKERR VecEqual(state, old, &equal);
    if (!equal)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Pre-step changed H1 incremental kinematics before transfer");
    CHKERR VecEqual(ep.solTSStep, history, &equal);
    if (!equal)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Pre-step changed H1 physical history");
    MoFEMFunctionReturnHot(0);
  }
  PetscBool check_schur = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(nullptr, nullptr, "-test_rotation_schur",
                             &check_schur, nullptr);
  boost::shared_ptr<EshelbianCore::SetUpSchur> schur;
  if (check_schur) {
    CHKERR ep.setElasticElementOps(1);
    CHKERR ep.setElasticElementToTs(ep.dmElastic);
    CHKERR TSSetOptionsPrefix(ts, "elastic_");
    CHKERR TSSetFromOptions(ts);
    schur = EshelbianCore::SetUpSchur::createSetUpSchur(ep.mField, &ep);
    CHKERR schur->setUp(ts);
  }
  const auto check_schur_hierarchy = [&]() {
    MoFEMFunctionBegin;
    if (!check_schur)
      MoFEMFunctionReturnHot(0);
    SNES snes;
    KSP ksp, *subksp;
    PC pc, a00_pc, schur_pc;
    CHKERR TSGetSNES(ts, &snes);
    CHKERR SNESGetKSP(snes, &ksp);
    CHKERR KSPGetPC(ksp, &pc);
    // Reassemble through the same TS Jacobian callback used by Newton solves.
    auto rates = vectorDuplicate(state);
    CHKERR VecZeroEntries(rates);
    Mat jacobian, preconditioner;
    CHKERR TSGetIJacobian(ts, &jacobian, &preconditioner, nullptr, nullptr);
    CHKERR TSComputeIJacobian(ts, 0., state, rates, 1., jacobian,
                              preconditioner, PETSC_FALSE);
    CHKERR checkRotationPreconditioner(ep, ts, state, rates, jacobian,
                                       preconditioner);
    CHKERR KSPSetOperators(ksp, jacobian, preconditioner);
    CHKERR KSPSetUp(ksp);
    CHKERR PCFieldSplitSchurGetSubKSP(pc, nullptr, &subksp);
    CHKERR KSPGetPC(subksp[0], &a00_pc);
    CHKERR KSPGetPC(subksp[1], &schur_pc);
    CHKERR PetscFree(subksp);
    PetscBool shell;
    CHKERR PetscObjectTypeCompare(getPetscObject(a00_pc), PCSHELL, &shell);
    if (!shell)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Rotation shortening destroyed the Schur A00 shell PC");
    DM dm;
    CHKERR PCGetDM(schur_pc, &dm);
    const Problem *problem;
    CHKERR DMMoFEMGetProblemPtr(dm, &problem);
    if (problem->getName() != "SUB_SCHUR")
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Rotation shortening lost the Schur DM");
    MoFEMFunctionReturn(0);
  };
  CHKERR check_schur_hierarchy();
  SNES snes;
  CHKERR TSGetSNES(ts, &snes);
  CHKERR SNESSetLagJacobian(snes, 3);
  CHKERR SNESSetLagPreconditioner(snes, 4);
  CHKERR VecCopy(old, state);
  CHKERR TSPreStep(ts);
  CHKERR check_schur_hierarchy();
  PetscBool equal;
  CHKERR VecEqual(state, old, &equal);
  if (equal)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Expected a constant rotation branch change");
  auto difference = vectorDuplicate(state);
  CHKERR VecCopy(state, difference);
  CHKERR VecAXPY(difference, -1., old);
  const double *array;
  CHKERR VecGetArrayRead(difference, &array);
  for (const auto &dof : *problem->getNumeredRowDofsPtr()) {
    const auto index = dof->getPetscGlobalDofIdx();
    if (index < first || index >= last)
      continue;
    if ((dof->getName() != ep.rotAxis || dof->getDofOrder() != 0) &&
        array[index - first] != 0.)
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Branch update changed a higher mode or another field");
  }
  CHKERR VecRestoreArrayRead(difference, &array);
  auto shortened = vectorDuplicate(state);
  CHKERR VecCopy(state, shortened);
  CHKERR TSPreStep(ts);
  CHKERR VecEqual(state, shortened, &equal);
  if (!equal)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Short constant rotation should not be changed twice");
  CHKERR check_schur_hierarchy();
  PetscInt lag;
  CHKERR SNESGetLagJacobian(snes, &lag);
  if (lag != 3)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Rotation shortening changed the configured Jacobian lag");
  CHKERR SNESGetLagPreconditioner(snes, &lag);
  if (lag != 4)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Rotation shortening changed the configured preconditioner lag");
  CHKERR VecEqual(ep.solTSStep, history, &equal);
  if (!equal)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "TS pre-step replaced the distinct previous physical state");
  for (double *penalty : {&ep.alphaR, &ep.alphaViscousR}) {
    const double saved = *penalty;
    *penalty = .3;
    CHKERR VecCopy(old, state);
    CHKERR TSPreStep(ts);
    CHKERR VecEqual(state, old, &equal);
    *penalty = saved;
    if (!equal)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Raw coordinate penalty must prevent branch shortening");
  }
  CHKERR TSPreStep(ts);
  auto fe = boost::make_shared<VolumeElementForcesAndSourcesCore>(ep.mField);
  fe->getUserPolynomialBase() =
      boost::make_shared<CGGUserPolynomialBase>(nullptr, true);
  fe->getRuleHook = [](int, int, int order) { return 2 * order + 2; };
  CHKERR MoFEM::AddHOOps<3, 3, 3>::add(fe->getOpPtrVector(), {L2},
                                       ep.materialH1Positions);
  std::array<MatrixPtr, 2> a, q;
  for (int mode = 0; mode != 2; ++mode) {
    a[mode] = boost::make_shared<MatrixDouble>();
    q[mode] = boost::make_shared<MatrixDouble>();
    fe->getOpPtrVector().push_back(new OpCalculateL2RotationValues(
        ep.rotAxis, boost::make_shared<MatrixDouble>(), a[mode], q[mode],
        mode ? state : old));
  }
  fe->getOpPtrVector().push_back(new OpCheckShortenedRotation(a, q));
  CHKERR DMoFEMLoopFiniteElements(ep.dmElastic, ep.elementVolumeName, fe);
  MoFEMFunctionReturn(0);
}
 
/** Verify product-rotation tangents with a nontrivial known H1 deformation. */
MoFEMErrorCode checkH1RotationJacobian(EshelbianCore &ep) {
  MoFEMFunctionBegin;
  CHKERR ep.mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(
      [&](boost::shared_ptr<FieldEntity> entity) {
        MoFEMFunctionBegin;
        auto values = entity->getEntFieldData();
        for (int c = 0; c != values.size(); ++c)
          values[c] = 0.;
        if (entity->getEntType() == MBVERTEX) {
          double coordinates[3];
          const auto handle = entity->getEnt();
          CHKERR ep.mField.get_moab().get_coords(&handle, 1, coordinates);
          auto x = getFTensor1FromPtr<3>(coordinates);
          auto displacement = getFTensor1FromPtr<3>(&values[0]);
          displacement(0) = .04 * x(0) + .02 * x(1);
          displacement(1) = -.03 * x(1) + .01 * x(2);
          displacement(2) = .02 * x(2);
        }
        MoFEMFunctionReturn(0);
      },
      ep.spatialH1Disp);
  auto state = createDMVector(ep.dmElastic);
  CHKERR VecZeroEntries(state);
  const Problem *problem;
  CHKERR DMMoFEMGetProblemPtr(ep.dmElastic, &problem);
  PetscInt first, last;
  CHKERR VecGetOwnershipRange(state, &first, &last);
  for (const auto &dof : *problem->getNumeredRowDofsPtr()) {
    const auto index = dof->getPetscGlobalDofIdx();
    if (index < first || index >= last)
      continue;
    const auto &name = dof->getName();
    if (name == ep.rotAxis || name == ep.stretchTensor ||
        name == ep.piolaStress)
      CHKERR VecSetValue(state, index, .01 * (index % 13 - 6), INSERT_VALUES);
  }
  CHKERR VecAssemblyBegin(state);
  CHKERR VecAssemblyEnd(state);
  CHKERR VecGhostUpdateBegin(state, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(state, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecCopy(state, ep.solTSStep);
  auto rates = vectorDuplicate(state);
  CHKERR VecZeroEntries(rates);
  boost::shared_ptr<VolumeElementForcesAndSourcesCore> rhs, lhs;
  CHKERR ep.setVolumeElementOps(1, true, false, rhs, lhs);
  auto tester = ep.mField.getInterface<OperatorsTester>();
  auto jacobian = tester->assembleMat(ep.dmElastic, ep.elementVolumeName, lhs,
                                      state, rates, rates, 1., 1., {});
  auto residual = tester->assembleVec(ep.dmElastic, ep.elementVolumeName, rhs,
                                      state, rates, rates, 1., 1., {});
  auto direction = vectorDuplicate(state);
  auto error = vectorDuplicate(state);
  auto action = vectorDuplicate(state);
  auto trial = vectorDuplicate(state);
  const OperatorsTester::VectorFunction evaluate = [&](Vec input, Vec output) {
    MoFEMFunctionBegin;
    CHKERR VecCopy(input, trial);
    CHKERR VecGhostUpdateBegin(trial, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(trial, INSERT_VALUES, SCATTER_FORWARD);
    auto value = tester->assembleVec(ep.dmElastic, ep.elementVolumeName, rhs,
                                     trial, rates, rates, 1., 1., {});
    CHKERR VecCopy(value, output);
    MoFEMFunctionReturn(0);
  };
  for (const auto &field : {ep.rotAxis, ep.stretchTensor, ep.piolaStress}) {
    CHKERR VecCopy(state, direction);
    // Keep only the selected unknown block to expose each cross coupling.
    double *values;
    CHKERR VecGetArray(direction, &values);
    for (const auto &dof : *problem->getNumeredRowDofsPtr()) {
      const auto index = dof->getPetscGlobalDofIdx();
      if (index >= first && index < last && dof->getName() != field)
        values[index - first] = 0.;
    }
    CHKERR VecRestoreArray(direction, &values);
    CHKERR tester->checkVectorCentralFiniteDifference(
        state, direction, residual, jacobian, 2.e-6, evaluate, error);
    CHKERR MatMult(jacobian, direction, action);
    PetscReal error_norm, action_norm;
    CHKERR VecNorm(error, NORM_INFINITY, &error_norm);
    CHKERR VecNorm(action, NORM_INFINITY, &action_norm);
    if (!std::isfinite(error_norm) || error_norm > 2.e-8 + 2.e-6 * action_norm)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "H1 split Jacobian for %s: error %g, action %g", field.c_str(),
              error_norm, action_norm);
    CHKERR PetscPrintf(ep.mField.get_comm(), "H1 split Jacobian %s: %.3e\n",
                       field.c_str(), error_norm);
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode runRotationSplitAtom(EshelbianCore &ep) {
  MoFEMFunctionBegin;
  const auto *field = ep.mField.get_field_structure(ep.rotAxis);
  if (field->getSpace() != L2 || field->getApproxBase() != USER_BASE ||
      field->getNbOfCoeffs() != 3)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Rotation split requires three-component USER_BASE L2");
  auto fe = boost::make_shared<VolumeElementForcesAndSourcesCore>(ep.mField);
  fe->getUserPolynomialBase() =
      boost::make_shared<CGGUserPolynomialBase>(nullptr, true);
  fe->getRuleHook = [](int, int, int order) { return 2 * order + 2; };
  CHKERR MoFEM::AddHOOps<3, 3, 3>::add(fe->getOpPtrVector(), {L2},
                                       ep.materialH1Positions);
  auto visits = boost::make_shared<PetscInt>(0);
  fe->getOpPtrVector().push_back(new OpCheckRotationLayout(ep.rotAxis, visits));
  CHKERR DMoFEMLoopFiniteElements(ep.dmElastic, ep.elementVolumeName, fe);
  PetscInt global_visits;
  CHKERR MPI_Allreduce(visits.get(), &global_visits, 1, MPIU_INT, MPI_SUM,
                       ep.mField.get_comm());
  if (!global_visits)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "No rotation layout evaluations executed");
  CHKERR checkRotationVectors(ep);
  CHKERR checkSplitRotationDerivatives();
  CHKERR checkConstantBranchUpdate(ep);
  PetscBool check_h1 = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(nullptr, nullptr, "-test_rotation_h1_jacobian",
                             &check_h1, nullptr);
  if (check_h1)
    CHKERR checkH1RotationJacobian(ep);
  CHKERR PetscPrintf(ep.mField.get_comm(),
                     "Rotation constant-mode layout verified on %d elements\n",
                     static_cast<int>(global_visits));
  MoFEMFunctionReturn(0);
}
} // namespace
 
int main(int argc, char *argv[]) {
  MoFEM::Core::Initialize(&argc, &argv, nullptr,
                          "Verify L2 constant-higher rotation split.\n");
  auto log = logging::core::get();
  for (const auto *name : {"EP", "EPSELF", "EPSYNC"}) {
    log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), name));
    LogManager::setLog(name);
  }
  try {
    CHKERR DMRegister_MGViaApproxOrders("DMMOFEM_MG");
    CHKERR runAuxiliaryLogarithmicStressTest(runRotationSplitAtom);
  }
  CATCH_ERRORS;
  MoFEM::Core::Finalize();
  return 0;
}