NonlinearElasticExample.hpp

Go to the documentation of this file.

/**
 * @file NonlinearElasticExample.hpp
 * @example mofem/tutorials/vec-2_nonlinear_elasticity/src/NonlinearElasticExample.hpp
 * @date 2025-10-5
 *
 * @copyright Anonymous authors (c) 2025 under the MIT license (see LICENSE for
 * details)
 */
 
#ifdef WITH_ADOL_C
  #include <MatOps.hpp>
  #include <MatElastic.hpp>
  #include <MatHuHu.hpp>
  #include <MatUmat.hpp>
  #include <MatGenericElastic.hpp>
  #include <MatOpsExampleGeneric.hpp>
  #ifdef WITH_ADOL_C_SIMPLE_DAMAGE
    #include <MatGenericElasticSimpleDamage.hpp>
  #endif
#endif
 
struct NonlinearElasticExample {
 
  NonlinearElasticExample(MoFEM::Interface &m_field) : mField(m_field) {}
 
  MoFEMErrorCode runProblem();
 
protected:
  MoFEM::Interface &mField;
  boost::shared_ptr<DomainEle> updateFe;
  boost::shared_ptr<PostProcEleDomain> postProcDomainFe;
  boost::shared_ptr<PostProcEleBdy> postProcBdyFe;
#ifdef AXISYMMETRIC
  static constexpr int MAT_DIM = 3;
  PetscBool useAdolcMaterial = PETSC_TRUE;
#else
  static constexpr int MAT_DIM = SPACE_DIM;
  PetscBool useAdolcMaterial = PETSC_FALSE;
#endif
  FieldApproximationBase approximationBase = AINSWORTH_LEGENDRE_BASE;
 
  MoFEMErrorCode readMesh();
  MoFEMErrorCode setupProblem();
  MoFEMErrorCode boundaryCondition();
  MoFEMErrorCode assembleSystemHencky();
  MoFEMErrorCode assembleSystemAdolc();
  MoFEMErrorCode opTest();
  MoFEMErrorCode TsSolve();
  MoFEMErrorCode gettingNorms();
  MoFEMErrorCode outputResults();
  MoFEMErrorCode checkResults();
 
  struct TSElasticPostStep {
 
    static MoFEMErrorCode
    postStepInitialise(NonlinearElasticExample *example_ptr) {
      elasticExamplePtr = example_ptr;
      return 0;
    }
    static MoFEMErrorCode postStepDestroy() {
      elasticExamplePtr = nullptr;
      return 0;
    }
    static MoFEMErrorCode preStepFun(TS ts) { return 0; }
    static MoFEMErrorCode postStepFun(TS ts) {
      MoFEMFunctionBegin;
      (void)ts;
      auto simple = elasticExamplePtr->mField.getInterface<Simple>();
      auto dm = simple->getDM();
      if (elasticExamplePtr->updateFe) {
        CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
                                        elasticExamplePtr->updateFe);
      }
      MoFEMFunctionReturn(0);
    }
 
  private:
    static inline NonlinearElasticExample *elasticExamplePtr;
  };
 
  std::vector<Tag>
      listTagsToTransfer; ///< list of tags to transfer to postprocessor
 
#ifdef WITH_ADOL_C
  #ifdef AXISYMMETRIC
  static inline constexpr int modelType = MatOps::MODEL_AXISYMMETRIC;
  #else
  static inline constexpr int modelType =
      (SPACE_DIM == 3) ? MatOps::MODEL_3D : MatOps::MODEL_2D_PLANE_STRAIN;
  #endif
  boost::shared_ptr<MatOps::PhysicalEquations> physicalEquationsPtr;
  boost::shared_ptr<MatOps::PhysicalEquations> physicalHuHuPtr;
#endif
 
  static auto get_skin(MoFEM::Interface &m_field, Range body_ents) {
    Skinner skin(&m_field.get_moab());
    Range skin_ents;
    CHK_MOAB_THROW(skin.find_skin(0, body_ents, false, skin_ents), "find_skin");
    return skin_ents;
  };
 
  static auto filter_true_skin(MoFEM::Interface &m_field, Range &&skin) {
    Range boundary_ents;
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
    CHK_MOAB_THROW(pcomm->filter_pstatus(skin,
                                         PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                         PSTATUS_NOT, -1, &boundary_ents),
                   "filter_pstatus");
    return boundary_ents;
  };
};
 
//! [Run problem]
MoFEMErrorCode NonlinearElasticExample::runProblem() {
  MoFEMFunctionBegin;
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR boundaryCondition();
 
  PetscBool test_op = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test_op", &test_op,
                             PETSC_NULLPTR);
 
  if (useAdolcMaterial == PETSC_TRUE) {
    CHKERR assembleSystemAdolc();
  } else {
    CHKERR assembleSystemHencky();
  }
 
  if (test_op == PETSC_TRUE) {
    CHKERR opTest();
  }
  CHKERR TsSolve();
  CHKERR gettingNorms();
  CHKERR outputResults();
  CHKERR checkResults();
  MoFEMFunctionReturn(0);
}
//! [Run problem]
 
//! [Read mesh]
MoFEMErrorCode NonlinearElasticExample::readMesh() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  CHKERR simple->getOptions();
  CHKERR simple->loadFile();
  CHKERR mField.getInterface<MeshsetsManager>()->setMeshsetFromFile();
  MoFEMFunctionReturn(0);
}
//! [Read mesh]
 
//! [Set up problem]
MoFEMErrorCode NonlinearElasticExample::setupProblem() {
  MoFEMFunctionBegin;
  Simple *simple = mField.getInterface<Simple>();
 
  enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
  const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
  PetscInt choice_base_value = AINSWORTH;
  CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-base", list_bases,
                              LASBASETOPT, &choice_base_value, PETSC_NULLPTR);
 
  switch (choice_base_value) {
  case AINSWORTH:
    approximationBase = AINSWORTH_LEGENDRE_BASE;
    MOFEM_LOG("WORLD", Sev::inform)
        << "Set AINSWORTH_LEGENDRE_BASE for displacements";
    break;
  case DEMKOWICZ:
    approximationBase = DEMKOWICZ_JACOBI_BASE;
    MOFEM_LOG("WORLD", Sev::inform)
        << "Set DEMKOWICZ_JACOBI_BASE for displacements";
    break;
  default:
    approximationBase = LASTBASE;
    break;
  }
 
  // Add field
  CHKERR simple->addDomainField("U", H1, approximationBase, SPACE_DIM);
  CHKERR simple->addBoundaryField("U", H1, approximationBase, SPACE_DIM);
  CHKERR simple->addDataField("GEOMETRY", H1, approximationBase, SPACE_DIM);
  int order = 2;
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
  CHKERR simple->setFieldOrder("U", order);
  CHKERR simple->setFieldOrder("GEOMETRY", 2);
  CHKERR simple->setUp();
 
  auto project_ho_geometry = [&]() {
    Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
    return mField.loop_dofs("GEOMETRY", ent_method);
  };
  CHKERR project_ho_geometry();
 
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-use_adolc_material",
                             &useAdolcMaterial, PETSC_NULLPTR);
#ifndef WITH_ADOL_C
  if (useAdolcMaterial == PETSC_TRUE) {
    SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
            "ADOL-C support is not enabled. Please reconfigure with "
            "-DWITH_ADOL_C=ON and recompile to use AdolC material model.");
  }
#endif
  if (useAdolcMaterial == PETSC_TRUE) {
    MOFEM_LOG("WORLD", Sev::inform) << "Use ADOL-C material model";
  } else {
    MOFEM_LOG("WORLD", Sev::inform) << "Use Hencky material model";
  }
 
  auto tag_meshest = [&]() {
    MoFEMFunctionBegin;
 
    auto set_block = [&](auto name, int dim) {
      std::map<int, Range> map;
      auto set_tag_impl = [&](auto name) {
        MoFEMFunctionBegin;
        auto mesh_mng = mField.getInterface<MeshsetsManager>();
        auto bcs = mesh_mng->getCubitMeshsetPtr(
 
            std::regex((boost::format("%s(.*)") % name).str())
 
        );
        std::map<int, int> ids_map;
        int bit_id = 1;
        for (auto bc : bcs) {
          int id = bc->getMeshsetId();
          ids_map[id] = bit_id;
          bit_id <<= 1;
        }
        for (auto bc : bcs) {
          Range r;
          CHKERR bc->getMeshsetIdEntitiesByDimension(mField.get_moab(), dim, r,
                                                     true);
          map[ids_map[bc->getMeshsetId()]] = r;
          MOFEM_LOG("EXAMPLE", Sev::inform)
              << "Block " << name << " id " << bc->getMeshsetId() << " : "
              << ids_map[bc->getMeshsetId()] << " has " << r.size()
              << " entities";
        }
        MoFEMFunctionReturn(0);
      };
 
      CHKERR set_tag_impl(name);
 
      return std::make_pair(name, map);
    };
 
    auto set_skin = [&](auto &&map) {
      for (auto &m : map.second) {
        auto s = filter_true_skin(mField, get_skin(mField, m.second));
        m.second.swap(s);
        MOFEM_LOG("EXAMPLE", Sev::inform)
            << "Skin for block " << map.first << " id " << m.first << " has "
            << m.second.size() << " entities";
      }
      return map;
    };
 
    auto set_tag = [&](auto &&map) {
      Tag th;
      auto name = map.first;
      int def_val[] = {-1};
      CHK_MOAB_THROW(mField.get_moab().tag_get_handle(
                         name, 1, MB_TYPE_INTEGER, th,
                         MB_TAG_SPARSE | MB_TAG_CREAT, def_val),
                     "create tag");
      for (auto &m : map.second) {
        int id = m.first;
        for (auto ent : m.second) {
          int current_id;
          CHK_MOAB_THROW(
              mField.get_moab().tag_get_data(th, &ent, 1, &current_id),
              "get tag data");
          if (current_id != -1) {
            id |= current_id;
          }
          CHK_MOAB_THROW(mField.get_moab().tag_set_data(th, &ent, 1, &id),
                         "set tag data");
        }
      }
      return th;
    };
 
    if (SPACE_DIM == 3) {
      listTagsToTransfer.push_back(set_tag(set_block("MAT_", 3)));
      listTagsToTransfer.push_back(set_tag(set_skin(set_block("MAT_", 3))));
    } else {
      listTagsToTransfer.push_back(set_tag(set_block("MAT_", 2)));
    }
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR tag_meshest();
 
  MoFEMFunctionReturn(0);
}
//! [Set up problem]
 
//! [Boundary condition]
 
MoFEMErrorCode NonlinearElasticExample::boundaryCondition() {
  MoFEMFunctionBegin;
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
  auto simple = mField.getInterface<Simple>();
  auto bc_mng = mField.getInterface<BcManager>();
  auto time_scale = boost::make_shared<ExampleTimeScale>();
 
  CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpForce>::add(
      pipeline_mng->getOpBoundaryRhsPipeline(), mField, "U", {time_scale},
      "FORCE", "PRESSURE", Sev::inform);
 
  //! [Define gravity vector]
  CHKERR DomainNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
      pipeline_mng->getOpDomainRhsPipeline(), mField, "U", {time_scale},
      "BODY_FORCE", Sev::inform);
 
  // Essential BC
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
                                           "U", 0, 0);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
                                           "U", 1, 1);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
                                           "U", 2, 2);
  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
      simple->getProblemName(), "U");
 
  MoFEMFunctionReturn(0);
}
//! [Boundary condition]
 
//! [Push operators to pipeline]
MoFEMErrorCode NonlinearElasticExample::assembleSystemHencky() {
  MoFEMFunctionBegin;
 
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
 
  auto integration_rule = [](int, int, int approx_order) {
    return 2 * (approx_order - 1) + 1;
  };
 
  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
  CHKERR pipeline_mng->setEvaluationIntegrationRule(integration_rule);
 
  auto add_domain_ops_lhs = [&](auto &pip) {
    MoFEMFunctionBegin;
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
                                                          "GEOMETRY");
    CHKERR
    HenckyOps::opFactoryDomainLhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
        mField, pip, "U", "MAT_ELASTIC", Sev::inform);
    MoFEMFunctionReturn(0);
  };
 
  auto add_domain_ops_rhs = [&](auto &pip) {
    MoFEMFunctionBegin;
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
                                                          "GEOMETRY");
    CHKERR
    HenckyOps::opFactoryDomainRhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
        mField, pip, "U", "MAT_ELASTIC", Sev::inform);
    MoFEMFunctionReturn(0);
  };
 
  CHKERR add_domain_ops_lhs(pipeline_mng->getOpDomainLhsPipeline());
  CHKERR add_domain_ops_rhs(pipeline_mng->getOpDomainRhsPipeline());
 
  // push operators to evaluation pipeline
  CHKERR add_domain_ops_rhs(pipeline_mng->getOpEvaluationPipeline());
  auto &domain_reaction_fe = pipeline_mng->getEvaluationFE();
  domain_reaction_fe->postProcessHook =
      EssentialPreProcReaction<DisplacementCubitBcData>(mField,
                                                        domain_reaction_fe);
 
  PetscBool post_proc_vol;
  PetscBool post_proc_skin;
 
  if constexpr (SPACE_DIM == 2) {
    post_proc_vol = PETSC_TRUE;
    post_proc_skin = PETSC_FALSE;
  } else {
    post_proc_vol = PETSC_FALSE;
    post_proc_skin = PETSC_TRUE;
  }
 
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_vol",
                             &post_proc_vol, PETSC_NULLPTR);
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_skin",
                             &post_proc_skin, PETSC_NULLPTR);
 
  // Setup postprocessing
  auto create_post_proc_fe = [&]() {
    auto post_proc_ele_domain = [this](auto &pip_domain) {
      CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip_domain, {H1},
                                                            "GEOMETRY");
      auto common_ptr = commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
          mField, pip_domain, "U", "MAT_ELASTIC", Sev::inform);
      return common_ptr;
    };
 
    auto post_proc_map = [&](auto &pip, auto u_ptr, auto common_ptr) {
      MoFEMFunctionBegin;
 
      using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
      pip->getOpPtrVector().push_back(
 
          new OpPPMap(pip->getPostProcMesh(), pip->getMapGaussPts(), {},
                      {{"U", u_ptr}},
                      {{"GRAD", common_ptr->matGradPtr},
                       {"FIRST_PIOLA", common_ptr->getMatFirstPiolaStress()}},
                      {}));
      MoFEMFunctionReturn(0);
    };
 
    auto push_post_proc_bdy = [&](auto &pip_bdy) {
      if (post_proc_skin == PETSC_FALSE)
        return boost::shared_ptr<PostProcEleBdy>();
      auto simple = mField.getInterface<Simple>();
      auto domain_fe_name = simple->getDomainFEName();
      auto u_ptr = boost::make_shared<MatrixDouble>();
      pip_bdy->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
      auto op_loop_side =
          new OpLoopSide<SideEle>(mField, domain_fe_name, SPACE_DIM);
      auto common_ptr = post_proc_ele_domain(op_loop_side->getOpPtrVector());
      pip_bdy->getOpPtrVector().push_back(op_loop_side);
      CHKERR post_proc_map(pip_bdy, u_ptr, common_ptr);
      return pip_bdy;
    };
 
    auto push_post_proc_domain = [&](auto &pip_domain) {
      if (post_proc_vol == PETSC_FALSE)
        return boost::shared_ptr<PostProcEleDomain>();
 
      auto u_ptr = boost::make_shared<MatrixDouble>();
      pip_domain->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
      auto common_ptr = post_proc_ele_domain(pip_domain->getOpPtrVector());
 
      CHKERR post_proc_map(pip_domain, u_ptr, common_ptr);
 
      return pip_domain;
    };
 
    auto post_proc_fe_domain = boost::make_shared<PostProcEleDomain>(mField);
    auto post_proc_fe_bdy = boost::make_shared<PostProcEleBdy>(mField);
 
    return std::make_pair(push_post_proc_domain(post_proc_fe_domain),
                          push_post_proc_bdy(post_proc_fe_bdy));
  };
 
  auto post_proc_pair = create_post_proc_fe();
  postProcDomainFe = post_proc_pair.first;
  postProcBdyFe = post_proc_pair.second;
 
  MoFEMFunctionReturn(0);
}
//! [Push operators to pipeline]
 
MoFEMErrorCode NonlinearElasticExample::opTest() {
  MoFEMFunctionBegin;
 
  // get operators tester
  auto simple = mField.getInterface<Simple>();
  auto opt = mField.getInterface<OperatorsTester>(); // get interface to
                                                     // OperatorsTester
  auto pip = mField.getInterface<PipelineManager>(); // get interface to
                                                     // pipeline manager
 
  constexpr double eps = 1e-9;
 
  auto x = opt->setRandomFields(simple->getDM(), {
 
                                                     {"U", {-1e-6, 1e-6}}
 
                                                 });
 
  auto dot_x = opt->setRandomFields(simple->getDM(), {
 
                                                         {"U", {-1, 1}}
 
                                                     });
 
  auto diff_x = opt->setRandomFields(simple->getDM(), {
 
                                                          {"U", {-1, 1}}
 
                                                      });
 
  auto test_domain_ops = [&](auto fe_name, auto lhs_pipeline,
                             auto rhs_pipeline) {
    MoFEMFunctionBegin;
 
    auto diff_res = opt->checkCentralFiniteDifference(
        simple->getDM(), fe_name, rhs_pipeline, lhs_pipeline, x, dot_x,
        SmartPetscObj<Vec>(), diff_x, 0, 0.5, eps);
 
    // Calculate norm of difference between directional derivative calculated
    // from finite difference, and tangent matrix.
    double fnorm;
    CHKERR VecNorm(diff_res, NORM_2, &fnorm);
    MOFEM_LOG_C("EXAMPLE", Sev::inform,
                "Test consistency of tangent matrix %3.4e", fnorm);
 
    CHKERR mField.getInterface<TestingInterface>()
        ->checkError<TestType::DERIVATIVE>("domain tangent", fnorm, 0.,
                                           {1e-5, 0.});
 
    MoFEMFunctionReturn(0);
  };
 
  MOFEM_LOG("EXAMPLE", Sev::inform)
      << "Test operators with finite difference of directional "
         "derivative";
  CHKERR test_domain_ops(simple->getDomainFEName(), pip->getDomainLhsFE(),
                         pip->getDomainRhsFE());
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode NonlinearElasticExample::assembleSystemAdolc() {
  MoFEMFunctionBegin;
 
  auto *simple = mField.getInterface<Simple>();
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
 
  auto integration_rule = [](int, int, int approx_order) {
    return 2 * (approx_order - 1) + 1;
  };
 
  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
  CHKERR pipeline_mng->setEvaluationIntegrationRule(integration_rule);
 
#ifdef WITH_ADOL_C
  auto mat_ops_data_ptr = MatOps::createMatOpsDataPtr();
 
  #ifdef WITH_ADOL_C_UMAT
  {
    constexpr int umat_model_type =
        (SPACE_DIM == 3) ? MatOps::MODEL_3D : MatOps::MODEL_2D_PLANE_STRAIN;
    auto umat_physical_equations_ptr =
        MatOps::createMatOpsPhysicalEquationsPtr<MatOps::UMAT, umat_model_type>(
            mat_ops_data_ptr, MatOps::MatOpsTagsRegistry::setTagName("Umat"));
 
    CHKERR umat_physical_equations_ptr->getOptions(&mField);
    CHKERR umat_physical_equations_ptr->recordTape();
 
    physicalEquationsPtr = umat_physical_equations_ptr;
  }
  #else
  if (!physicalEquationsPtr) {
    physicalEquationsPtr =
        MatOps::createMatOpsPhysicalEquationsPtr<MatOps::PIOLA_RESPONSE::META,
                                                 modelType>(
            mat_ops_data_ptr,
            MatOps::MatOpsTagsRegistry::setTagName("MetaPiolaResponse", 0));
  }
 
    #ifdef WITH_ADOL_C_GENERIC_ELASTIC
  {
    auto generic_physical_equations_ptr =
        MatOps::createMatOpsPhysicalEquationsPtr<MatOps::GenericElastic,
                                                 modelType>(
            mat_ops_data_ptr,
            MatOps::MatOpsTagsRegistry::setTagName("GenericElastic"));
 
    MatOpsExampleGeneric();
    // In example, lets do a simple linear elastic material, so we can verify the tangent matrix
    generic_physical_equations_ptr->hookEvaluateVariable =
        MatOpsExampleGeneric::hookEvaluateVariable;
    generic_physical_equations_ptr->hookEvaluateDerivatives =
        MatOpsExampleGeneric::hookEvaluateDerivatives;
    generic_physical_equations_ptr->hookUpdateState =
        MatOpsExampleGeneric::hookUpdateState;
 
    CHKERR generic_physical_equations_ptr->getOptions(&mField);
    CHKERR generic_physical_equations_ptr->recordTape();
 
    auto &material_map = getMetaPiolaResponseMap(
        boost::static_pointer_cast<MatOps::MatPiolaResponse>(
            physicalEquationsPtr));
    material_map[generic_physical_equations_ptr->tAg] =
        boost::static_pointer_cast<MatOps::MatPiolaResponse>(
            generic_physical_equations_ptr);
    physicalEquationsPtr->tAg = generic_physical_equations_ptr->tAg;
  }
    #endif // WITH_ADOL_C_GENERIC_ELASTIC
 
    #ifdef WITH_ADOL_C_SIMPLE_DAMAGE
  {
    auto generic_physical_equations_ptr =
        MatOps::createMatOpsPhysicalEquationsPtr<MatOps::GenericElastic,
                                                 modelType>(
            mat_ops_data_ptr, MatOps::MatOpsTagsRegistry::setTagName(
                                  "GenericElasticSimpleDamage"));
 
    CHKERR GenericElasticExamples::configureSimpleDamage(
        mField, boost::dynamic_pointer_cast<MatOps::GenericElastic>(
                    generic_physical_equations_ptr));
 
    CHKERR generic_physical_equations_ptr->getOptions(&mField);
    CHKERR generic_physical_equations_ptr->recordTape();
 
    auto &material_map = getMetaPiolaResponseMap(
        boost::static_pointer_cast<MatOps::MatPiolaResponse>(
            physicalEquationsPtr));
    material_map[generic_physical_equations_ptr->tAg] =
        boost::static_pointer_cast<MatOps::MatPiolaResponse>(
            generic_physical_equations_ptr);
    physicalEquationsPtr->tAg = generic_physical_equations_ptr->tAg;
  }
    #endif // WITH_ADOL_C_SIMPLE_DAMAGE
 
  if (physicalEquationsPtr) {
    CHKERR physicalEquationsPtr->getOptions(&mField);
    CHKERR physicalEquationsPtr->recordTape();
  }
 
  if (!physicalHuHuPtr && !mField.getInterface<MeshsetsManager>()
                               ->getCubitMeshsetPtr(std::regex("MAT_HUHU(.*)"))
                               .empty()) {
    const int tag_hu_hu = MatOps::MatOpsTagsRegistry::setTagName("HuHU");
    physicalHuHuPtr =
        MatOps::createMatOpsPhysicalEquationsPtr<MatOps::HUHU, modelType>(
            MatOps::createMatOpsDataPtr(), tag_hu_hu);
    CHKERR physicalHuHuPtr->getOptions(&mField);
    CHKERR physicalHuHuPtr->recordTape();
  }
  #endif   //  WITH_ADOL_C_UMAT
 
  auto add_domain_ops_lhs = [&](auto &pip) {
    MoFEMFunctionBegin;
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
                                                          "GEOMETRY");
    if (physicalEquationsPtr) {
      CHKERR MatOps::OpMaterialFactory<MatOps::PIOLA_RESPONSE, modelType>::
          template opLhsFactory<PETSC, GAUSS, DomainEleOp>(
              mField, pip, "U", physicalEquationsPtr);
    }
    if (physicalHuHuPtr) {
      CHKERR MatOps::OpMaterialFactory<MatOps::HUHU, modelType>::
          template opLhsFactory<PETSC, GAUSS, DomainEleOp>(
              mField, pip, simple->getDomainFEName(), "U", physicalHuHuPtr);
    }
    MoFEMFunctionReturn(0);
  };
 
  auto add_domain_ops_rhs = [&](auto &pip) {
    MoFEMFunctionBegin;
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
                                                          "GEOMETRY");
    if (physicalEquationsPtr) {
      CHKERR MatOps::OpMaterialFactory<MatOps::PIOLA_RESPONSE, modelType>::
          template opRhsFactory<PETSC, GAUSS, DomainEleOp>(
              mField, pip, "U", physicalEquationsPtr);
    }
    if (physicalHuHuPtr) {
      CHKERR MatOps::OpMaterialFactory<MatOps::HUHU, modelType>::
          template opRhsFactory<PETSC, GAUSS, DomainEleOp>(
              mField, pip, simple->getDomainFEName(), "U", physicalHuHuPtr);
    }
    MoFEMFunctionReturn(0);
  };
 
  auto add_domain_ops_update = [&](auto &pip) {
    MoFEMFunctionBegin;
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
                                                          "GEOMETRY");
    if (physicalEquationsPtr) {
      CHKERR MatOps::OpMaterialFactory<MatOps::PIOLA_RESPONSE, modelType>::
          opUpdateFactory(mField, pip, "U", physicalEquationsPtr);
    }
    MoFEMFunctionReturn(0);
  };
 
  CHKERR add_domain_ops_lhs(pipeline_mng->getOpDomainLhsPipeline());
  CHKERR add_domain_ops_rhs(pipeline_mng->getOpDomainRhsPipeline());
 
  updateFe = boost::make_shared<DomainEle>(mField);
  updateFe->getRuleHook = integration_rule;
  CHKERR add_domain_ops_update(updateFe->getOpPtrVector());
 
  CHKERR add_domain_ops_rhs(pipeline_mng->getOpEvaluationPipeline());
  auto &domain_reaction_fe = pipeline_mng->getEvaluationFE();
  domain_reaction_fe->postProcessHook =
      EssentialPreProcReaction<DisplacementCubitBcData>(mField,
                                                        domain_reaction_fe);
 
  PetscBool post_proc_vol;
  PetscBool post_proc_skin;
 
  if constexpr (SPACE_DIM == 2) {
    post_proc_vol = PETSC_TRUE;
    post_proc_skin = PETSC_FALSE;
  } else {
    post_proc_vol = PETSC_FALSE;
    post_proc_skin = PETSC_TRUE;
  }
 
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_vol",
                             &post_proc_vol, PETSC_NULLPTR);
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_skin",
                             &post_proc_skin, PETSC_NULLPTR);
 
  // Setup postprocessing
  auto create_post_proc_fe = [&]() {
    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
    auto post_proc_ele_domain = [&](auto &pip_domain) {
      MoFEMFunctionBegin;
      CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip_domain, {H1},
                                                            "GEOMETRY");
      if (physicalEquationsPtr) {
        CHKERR MatOps::OpMaterialFactory<MatOps::PIOLA_RESPONSE, modelType>::
            opPostProcFactory(mField, pip_domain, "U", physicalEquationsPtr);
      }
      MoFEMFunctionReturn(0);
    };
 
    auto post_proc_map = [&](auto &pip, auto u_ptr) {
      MoFEMFunctionBegin;
      if (physicalEquationsPtr) {
        auto grad_ptr =
            physicalEquationsPtr->matOpsDataPtr->getCommonDataPtr("grad");
        auto first_piola_ptr =
            physicalEquationsPtr->matOpsDataPtr->getCommonDataPtr("P");
        pip->getOpPtrVector().push_back(new OpPPMap(pip->getPostProcMesh(),
                                                    pip->getMapGaussPts(), {},
                                                    {{"U", u_ptr}}, {}, {}));
        pip->getOpPtrVector().push_back(
            new OpPostProcMapInMoab<MAT_DIM, MAT_DIM>(
                pip->getPostProcMesh(), pip->getMapGaussPts(), {}, {},
                {{"GRAD", grad_ptr}, {"FIRST_PIOLA", first_piola_ptr}}, {}));
      }
      MoFEMFunctionReturn(0);
    };
 
    auto push_post_proc_bdy = [&](auto &pip_bdy) {
      if (post_proc_skin == PETSC_FALSE)
        return boost::shared_ptr<PostProcEleBdy>();
      auto simple = mField.getInterface<Simple>();
      auto domain_fe_name = simple->getDomainFEName();
      auto u_ptr = boost::make_shared<MatrixDouble>();
      pip_bdy->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
      auto op_loop_side =
          new OpLoopSide<SideEle>(mField, domain_fe_name, SPACE_DIM);
      CHKERR post_proc_ele_domain(op_loop_side->getOpPtrVector());
      pip_bdy->getOpPtrVector().push_back(op_loop_side);
      CHKERR post_proc_map(pip_bdy, u_ptr);
      return pip_bdy;
    };
 
    auto push_post_proc_domain = [&](auto &pip_domain) {
      if (post_proc_vol == PETSC_FALSE)
        return boost::shared_ptr<PostProcEleDomain>();
      auto u_ptr = boost::make_shared<MatrixDouble>();
      pip_domain->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
      CHKERR post_proc_ele_domain(pip_domain->getOpPtrVector());
 
      if (physicalEquationsPtr) {
        auto state_tags = physicalEquationsPtr->matOpsDataPtr->getStateTags();
        if (!state_tags.empty()) {
          MatOps::StateTags selected_state_tags;
          MatOps::StateTags statev_tags;
          selected_state_tags.reserve(state_tags.size());
          statev_tags.reserve(state_tags.size());
 
          constexpr auto statev_prefix = "UMAT_STATEV_";
          constexpr auto statev_prefix_len = sizeof("UMAT_STATEV_") - 1;
          auto get_statev_index = [&](const std::string &name) {
            return std::stoi(name.substr(statev_prefix_len));
          };
 
          for (const auto &state_tag : state_tags) {
            const auto &name = state_tag.name;
 
            if (name == "UMAT_DFGRD0" || name == "UMAT_STRESS") {
              selected_state_tags.push_back(state_tag);
              continue;
            }
 
            if (name.rfind(statev_prefix, 0) == 0) {
              statev_tags.push_back(state_tag);
            }
          }
 
          std::sort(statev_tags.begin(), statev_tags.end(),
                    [&](const auto &lhs, const auto &rhs) {
                      return get_statev_index(lhs.name) <
                             get_statev_index(rhs.name);
                    });
          if (statev_tags.size() > 10)
            statev_tags.resize(10);
          selected_state_tags.insert(selected_state_tags.end(),
                                     statev_tags.begin(), statev_tags.end());
 
          if (selected_state_tags.empty())
            return pip_domain;
 
          auto simple = mField.getInterface<Simple>();
          int state_order = 2;
          CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &state_order,
                                    PETSC_NULLPTR);
          auto op_loop_this = new OpLoopThis<DomainEle>(
              mField, simple->getDomainFEName(), Sev::noisy);
          op_loop_this->getThisFEPtr()->getRuleHook = integration_rule;
          pip_domain->getOpPtrVector().push_back(op_loop_this);
          CHKERR addTagDGProjectionOps<SPACE_DIM, SPACE_DIM>(
              mField.get_moab(), pip_domain->getOpPtrVector(),
              op_loop_this->getThisFEPtr()->getOpPtrVector(),
              pip_domain->getPostProcMesh(), pip_domain->getMapGaussPts(),
              selected_state_tags, state_order, approximationBase);
        }
      }
 
      CHKERR post_proc_map(pip_domain, u_ptr);
 
      return pip_domain;
    };
 
    auto post_proc_fe_domain = boost::make_shared<PostProcEleDomain>(mField);
    auto post_proc_fe_bdy = boost::make_shared<PostProcEleBdy>(mField);
 
    return std::make_pair(push_post_proc_domain(post_proc_fe_domain),
                          push_post_proc_bdy(post_proc_fe_bdy));
  };
 
  auto post_proc_pair = create_post_proc_fe();
  postProcDomainFe = post_proc_pair.first;
  postProcBdyFe = post_proc_pair.second;
 
#else
  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
          "ADOL-C support is not enabled. Please reconfigure with "
          "-DWITH_ADOL_C=ON and recompile to use AdolC material model.");
 
#endif // WITH_ADOL
  MoFEMFunctionReturn(0);
}
 
//! [TS Solve]
MoFEMErrorCode NonlinearElasticExample::TsSolve() {
  MoFEMFunctionBegin;
  auto *simple = mField.getInterface<Simple>();
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
 
  auto dm = simple->getDM();
  auto ts = pipeline_mng->createTSIM();
 
  auto add_extra_finite_elements_to_solver_pipelines = [&]() {
    MoFEMFunctionBegin;
 
    auto pre_proc_ptr = boost::make_shared<FEMethod>();
    auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
    auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
 
    auto time_scale = boost::make_shared<ExampleTimeScale>();
 
    auto get_bc_hook_rhs = [this, pre_proc_ptr, time_scale]() {
      MoFEMFunctionBeginHot;
      CHKERR EssentialPreProc<DisplacementCubitBcData>(mField, pre_proc_ptr,
                                                       {time_scale}, false)();
      MoFEMFunctionReturnHot(0);
    };
 
    pre_proc_ptr->preProcessHook = get_bc_hook_rhs;
 
    auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {
      MoFEMFunctionBegin;
      CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
          mField, post_proc_rhs_ptr, nullptr, Sev::verbose)();
      CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(
          mField, post_proc_rhs_ptr, 1.)();
      MoFEMFunctionReturn(0);
    };
    auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
      MoFEMFunctionBegin;
      CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
          mField, post_proc_lhs_ptr, 1.)();
      MoFEMFunctionReturn(0);
    };
    post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
    post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs;
 
    // This is low level pushing finite elements (pipelines) to solver
    auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
    ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
    ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
    ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
    ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
    MoFEMFunctionReturn(0);
  };
 
  // Add extra finite elements to SNES solver pipelines to resolve essential
  // boundary conditions
  CHKERR add_extra_finite_elements_to_solver_pipelines();
 
  auto create_monitor_fe = [dm](auto &&post_proc_fe) {
    return boost::make_shared<Monitor>(dm.get(), post_proc_fe);
  };
 
  // Set monitor which postprocessing results and saves them to the hard drive
  boost::shared_ptr<FEMethod> null_fe;
  if (postProcDomainFe)
    CHKERR postProcDomainFe->setTagsToTransfer(
        std::vector<Tag>(listTagsToTransfer));
  if (postProcBdyFe)
    CHKERR postProcBdyFe->setTagsToTransfer(
        std::vector<Tag>(listTagsToTransfer));
  auto monitor_ptr =
      create_monitor_fe(std::make_pair(postProcDomainFe, postProcBdyFe));
  CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
                             null_fe, monitor_ptr);
 
  // Set time solver
  double ftime = 1;
  CHKERR TSSetMaxTime(ts, ftime);
  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
 
  auto B = createDMMatrix(dm);
  CHKERR TSSetI2Jacobian(ts, B, B, PETSC_NULLPTR, PETSC_NULLPTR);
  auto D = createDMVector(simple->getDM());
  CHKERR TSSetSolution(ts, D);
  CHKERR TSSetFromOptions(ts);
 
  PetscBool test_steps = PETSC_FALSE;
  PetscBool test_piola_stress = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(nullptr, "", "-test_piola_stress",
                             &test_piola_stress, nullptr);
  CHKERR PetscOptionsGetBool(nullptr, "", "-test_step_monitor", &test_steps,
                             nullptr);
 
  const bool check_steps = test_steps || test_piola_stress;
 
  auto check_piola_stress = [this](PetscReal time) {
    MoFEMFunctionBegin;
    auto simple = mField.getInterface<Simple>();
 
    FTensor::Tensor2<double, MAT_DIM, MAT_DIM> t_P_expected;
    FTensor::Tensor2<double, MAT_DIM, MAT_DIM> t_P_diff;
    FTENSOR_INDEX(MAT_DIM, i);
    FTENSOR_INDEX(MAT_DIM, J);
 
    auto displacement_ptr = boost::make_shared<MatrixDouble>();
 
    // These homogeneous cube fixtures prescribe Fzz = 1 - 0.1 * load time.
    const double applied_F_zz = 1. - 0.1 * time;
 
    t_P_expected(i, J) = 0.;
    auto set_linear_elastic_piola = [&](const double strain_zz,
                                        const double d_strain_d_F_zz) {
      constexpr double E = 1.;
      constexpr double nu = 0.3;
      const double lambda = E * nu / ((1. + nu) * (1. - 2. * nu));
      const double mu = E / (2. * (1. + nu));
      for (int d = 0; d != SPACE_DIM - 1; ++d)
        t_P_expected(d, d) = lambda * strain_zz;
      t_P_expected(SPACE_DIM - 1, SPACE_DIM - 1) =
          (lambda + 2. * mu) * strain_zz * d_strain_d_F_zz;
    };
 
#ifndef AXISYMMETRIC
    auto set_neohookean_piola = [&](const double F_zz, const double J) {
      constexpr double c10 = 1.;
      constexpr double K = 1.;
      double I1 = (SPACE_DIM - 1) + F_zz * F_zz;
      if constexpr (SPACE_DIM == 2)
        I1 += 1.; // Plane strain has F33 = 1.
      const double J_to_minus_two_thirds = std::pow(J, -2. / 3.);
      for (int d = 0; d != SPACE_DIM; ++d) {
        const double F_dd = d == SPACE_DIM - 1 ? F_zz : 1.;
        const double inv_F_dd = 1. / F_dd;
        t_P_expected(d, d) =
            2. * c10 * J_to_minus_two_thirds * (F_dd - (I1 / 3.) * inv_F_dd) +
            K * J * (J - 1.) * inv_F_dd;
      }
    };
 
#endif
 
    auto eval_piola_at_point = [&](boost::shared_ptr<MatrixDouble> &mat_p,
                                   std::array<double, 3> field_eval_coords) {
      MoFEMFunctionBegin;
      CHKERR mField.getInterface<TestingInterface>()
          ->addFieldEvaluator<SPACE_DIM>(
              simple->getProblemName(), simple->getDomainFEName(),
              field_eval_coords, [&](auto &pipeline) {
                MoFEMFunctionBegin;
                CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
                    pipeline, {H1}, "GEOMETRY");
 
                pipeline.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
                    "U", displacement_ptr));
 
                if (useAdolcMaterial == PETSC_FALSE) {
                  auto common_ptr =
                      commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
                          mField, pipeline, "U", "MAT_ELASTIC", Sev::inform);
                  mat_p = common_ptr->getMatFirstPiolaStress();
                } else {
#ifdef WITH_ADOL_C
                  CHKERR MatOps::OpMaterialFactory<
                      MatOps::PIOLA_RESPONSE,
                      modelType>::opPostProcFactory(mField, pipeline, "U",
                                                    physicalEquationsPtr);
                  mat_p = physicalEquationsPtr->matOpsDataPtr->getCommonDataPtr(
                      "P");
#else
              SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                      "ADOL-C support is not enabled. Please reconfigure with "
                      "-DWITH_ADOL_C=ON and recompile to use MatOps material model.");
#endif
                }
                MoFEMFunctionReturn(0);
              });
      MoFEMFunctionReturn(0);
    };
 
    boost::shared_ptr<MatrixDouble> mat_p;
    using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
    // Infer the homogeneous fixture's response from the configured material
    // dispatcher. The default MetaPiolaResponse tag represents Neo-Hookean.
    std::string material_name = "Hencky";
#ifdef WITH_ADOL_C
    if (useAdolcMaterial) {
      int material_tag = physicalEquationsPtr->tAg;
      bool found_range = false;
      if (physicalEquationsPtr->tagVsRangePtr) {
        for (const auto &[tag, range] : *physicalEquationsPtr->tagVsRangePtr) {
          if (range.empty())
            continue;
          if (found_range && tag != material_tag)
            SETERRQ(mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
                    "The analytical fixture requires a homogeneous material");
          material_tag = tag;
          found_range = true;
        }
      }
      material_name = MatOps::MatOpsTagsRegistry::getTagName(material_tag);
    }
#endif
 
#ifdef AXISYMMETRIC
    if (material_name != "MetaPiolaResponse")
      SETERRQ(
          mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
          "The axisymmetric analytical fixture requires Neo-Hookean material");
#endif
 
    if (material_name == "Hencky") {
      // Homogeneous cube: E = 1, nu = 0.3.
      set_linear_elastic_piola(std::log(applied_F_zz), 1. / applied_F_zz);
    } else if (material_name == "MetaPiolaResponse") {
#ifdef AXISYMMETRIC
      // Uniform expansion, c10 = K = 1, including hoop stretch. The
      // deviatoric stress vanishes; P = K J (J - 1) F^{-T}.
      const double stretch = 1. + 0.01 * time;
      const double jacobian = stretch * stretch * stretch;
      const double piola = jacobian * (jacobian - 1.) / stretch;
      FTensor::Kronecker_Delta<double> delta;
      t_P_expected(i, J) = piola * delta(i, J);
      CHKERR eval_piola_at_point(mat_p, {1.37, 0.41, 0.});
      auto t_u = MatrixSizeHelper<GetFTensor1FromMatType<2, -1, DL>, DL>::get(
          *displacement_ptr, 1)();
      const double u_err = std::hypot(t_u(0) - (stretch - 1.) * 1.37,
                                      t_u(1) - (stretch - 1.) * 0.41);
      CHKERR mField.getInterface<TestingInterface>()
          ->checkError<TestType::POINTWISE>("axisymmetric displacement", u_err,
                                            0., {1e-8, 0.});
#else
      set_neohookean_piola(applied_F_zz, applied_F_zz);
#endif
    } else if (material_name == "Umat") {
      // Reference for this tutorial's linear-elastic UMAT fixture only.
      char strain_type[32] = "small";
      CHKERR PetscOptionsGetString(nullptr, "", "-umat_strain_type",
                                   strain_type, sizeof(strain_type), nullptr);
      if (std::string(strain_type) != "small")
        SETERRQ(mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
                "The analytical UMAT fixture requires small strain");
      set_linear_elastic_piola(applied_F_zz - 1., 1.);
    } else {
      SETERRQ(mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "No analytical fixture for material %s", material_name.c_str());
    }
 
    if (!mat_p)
      CHKERR eval_piola_at_point(mat_p, {0.123, -0.087, 0.193});
    auto t_P =
        MatrixSizeHelper<GetFTensor2FromMatType<MAT_DIM, MAT_DIM, -1, DL>,
                         DL>::get(*mat_p, 1)();
    t_P_diff(i, J) = t_P(i, J) - t_P_expected(i, J);
    const double p_err = std::sqrt(t_P_diff(i, J) * t_P_diff(i, J));
    for (int r = 0; r != MAT_DIM; ++r)
      for (int c = 0; c != MAT_DIM; ++c)
        MOFEM_LOG("EXAMPLE", Sev::inform)
            << "Piola stress (" << material_name << ") P(" << r << ", " << c
            << ") actual " << std::scientific << t_P(r, c) << " expected "
            << t_P_expected(r, c) << " diff " << t_P_diff(r, c);
    MOFEM_LOG("EXAMPLE", Sev::inform) << "Piola stress (" << material_name
                                      << ") error " << std::scientific << p_err;
    CHKERR mField.getInterface<TestingInterface>()
        ->checkError<TestType::POINTWISE>("Piola stress", p_err, 0.,
                                          {1e-8, 0.});
 
    MoFEMFunctionReturn(0);
  };
 
  if (check_steps) {
    CHKERR mField.getInterface<TestingInterface>()->addTSMonitor(
        ts, [this, check_piola_stress](TS, PetscInt step, PetscReal time,
                                       Vec solution) {
          MoFEMFunctionBegin;
          auto simple = mField.getInterface<Simple>();
          CHKERR VecGhostUpdateBegin(solution, INSERT_VALUES, SCATTER_FORWARD);
          CHKERR VecGhostUpdateEnd(solution, INSERT_VALUES, SCATTER_FORWARD);
          CHKERR DMoFEMMeshToLocalVector(simple->getDM(), solution,
                                         INSERT_VALUES, SCATTER_REVERSE);
          CHKERR check_piola_stress(time);
          MOFEM_LOG("EXAMPLE", Sev::inform)
              << "Numerical checks passed at TS step " << step << ", time "
              << time;
          MoFEMFunctionReturn(0);
        });
  }
 
  CHKERR TSElasticPostStep::postStepInitialise(this);
  CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
  CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
  const auto solve_error = TSSolve(ts, nullptr);
  CHKERR TSElasticPostStep::postStepDestroy();
  CHKERR solve_error;
 
  // Solver convergence is cheap to check and carries no magic numbers, so it
  // stands alone right after the solve rather than being deferred to
  // checkResults() through stored state.
  {
    PetscReal expected_time;
    CHKERR TSGetMaxTime(ts, &expected_time);
    auto testing = mField.getInterface<TestingInterface>();
    CHKERR testing->checkSolver(ts, expected_time, {1e-12, 1e-12},
                                "nonlinear elasticity");
  }
 
  CHKERR TSGetTime(ts, &ftime);
 
  PetscInt steps, snesfails, rejects, nonlinits, linits;
  CHKERR TSGetStepNumber(ts, &steps);
  CHKERR TSGetSNESFailures(ts, &snesfails);
  CHKERR TSGetStepRejections(ts, &rejects);
  CHKERR TSGetSNESIterations(ts, &nonlinits);
  CHKERR TSGetKSPIterations(ts, &linits);
  MOFEM_LOG_C("EXAMPLE", Sev::inform,
              "steps %d (%d rejected, %d SNES fails), ftime %g, nonlinits "
              "%d, linits %d",
              steps, rejects, snesfails, ftime, nonlinits, linits);
 
  MoFEMFunctionReturn(0);
}
//! [TS Solve]
 
//! [Getting norms]
MoFEMErrorCode NonlinearElasticExample::gettingNorms() {
  MoFEMFunctionBegin;
 
  auto simple = mField.getInterface<Simple>();
  auto dm = simple->getDM();
 
  auto T = createDMVector(simple->getDM());
  CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
                                 SCATTER_FORWARD);
  double nrm2;
  CHKERR VecNorm(T, NORM_2, &nrm2);
  MOFEM_LOG("EXAMPLE", Sev::inform) << "Solution norm " << nrm2;
 
  auto post_proc_norm_fe = boost::make_shared<DomainEle>(mField);
 
  auto post_proc_norm_rule_hook = [](int, int, int p) -> int { return 2 * p; };
  post_proc_norm_fe->getRuleHook = post_proc_norm_rule_hook;
 
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      post_proc_norm_fe->getOpPtrVector(), {H1}, "GEOMETRY");
 
  enum NORMS { U_NORM_L2 = 0, PIOLA_NORM, LAST_NORM };
  auto norms_vec =
      createVectorMPI(mField.get_comm(),
                      (mField.get_comm_rank() == 0) ? LAST_NORM : 0, LAST_NORM);
  CHKERR VecZeroEntries(norms_vec);
 
  auto u_ptr = boost::make_shared<MatrixDouble>();
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
 
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalcNormL2Tensor1<SPACE_DIM>(u_ptr, norms_vec, U_NORM_L2));
 
  if (useAdolcMaterial == PETSC_TRUE) {
#ifdef WITH_ADOL_C
    if (physicalEquationsPtr) {
      CHKERR MatOps::OpMaterialFactory<MatOps::PIOLA_RESPONSE, modelType>::
          opPostProcFactory(mField, post_proc_norm_fe->getOpPtrVector(), "U",
                            physicalEquationsPtr);
      auto m_P = physicalEquationsPtr->matOpsDataPtr->getCommonDataPtr("P");
      post_proc_norm_fe->getOpPtrVector().push_back(
          new OpCalcNormL2Tensor2<MAT_DIM, MAT_DIM>(m_P, norms_vec,
                                                    PIOLA_NORM));
    }
#else
    SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
            "ADOL-C support is not enabled. Please reconfigure with "
            "-DWITH_ADOL_C=ON and recompile to use AdolC material model.");
#endif
  } else {
    auto common_ptr = commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
        mField, post_proc_norm_fe->getOpPtrVector(), "U", "MAT_ELASTIC",
        Sev::inform);
    post_proc_norm_fe->getOpPtrVector().push_back(
        new OpCalcNormL2Tensor2<SPACE_DIM, SPACE_DIM>(
            common_ptr->getMatFirstPiolaStress(), norms_vec, PIOLA_NORM));
  }
 
  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
                                  post_proc_norm_fe);
 
  CHKERR VecAssemblyBegin(norms_vec);
  CHKERR VecAssemblyEnd(norms_vec);
 
  MOFEM_LOG_CHANNEL("SELF"); // Clear channel from old tags
  if (mField.get_comm_rank() == 0) {
    const double *norms;
    CHKERR VecGetArrayRead(norms_vec, &norms);
    MOFEM_TAG_AND_LOG("SELF", Sev::inform, "example")
        << "norm_u: " << std::scientific << std::sqrt(norms[U_NORM_L2]);
    MOFEM_TAG_AND_LOG("SELF", Sev::inform, "example")
        << "norm_piola: " << std::scientific << std::sqrt(norms[PIOLA_NORM]);
    CHKERR VecRestoreArrayRead(norms_vec, &norms);
  }
 
  MoFEMFunctionReturn(0);
}
//! [Getting norms]
 
//! [Postprocessing results]
MoFEMErrorCode NonlinearElasticExample::outputResults() {
  MoFEMFunctionBegin;
  MoFEMFunctionReturn(0);
}
//! [Postprocessing results]
 
//! [Check]
MoFEMErrorCode NonlinearElasticExample::checkResults() {
  MoFEMFunctionBegin;
  auto testing = mField.getInterface<TestingInterface>();
 
  // -test_solution_norm <value> both supplies the golden value and is the
  // trigger for this check, so it is read right here, where it is used,
  // rather than through addTestCommandLine's own "-test_nb"/"-test_name"
  // selectors, and rather than a class member set earlier in runProblem().
  PetscReal expected_solution_norm = 0;
  PetscBool test_solution_norm = PETSC_FALSE;
  CHKERR PetscOptionsGetReal(nullptr, "", "-test_solution_norm",
                             &expected_solution_norm, &test_solution_norm);
  if (test_solution_norm) {
    CHKERR testing->addTest(1, "solution coefficient norm", [&]() {
      MoFEMFunctionBegin;
      auto simple = mField.getInterface<Simple>();
      auto T = createDMVector(simple->getDM());
      CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
                                     SCATTER_FORWARD);
      double nrm2;
      CHKERR VecNorm(T, NORM_2, &nrm2);
      MOFEM_LOG("EXAMPLE", Sev::verbose) << "Regression norm " << nrm2;
      CHKERR testing->checkValue<TestType::SCALAR>("solution coefficient norm",
                                                   nrm2, expected_solution_norm,
                                                   {1e-2, 0.});
      MoFEMFunctionReturn(0);
    });
  }
 
  CHKERR testing->runTests();
  MoFEMFunctionReturn(0);
}
//! [Check]