ElasticAdaptiveExample.hpp

Go to the documentation of this file.

/**
 * @file ElasticAdaptiveExample.hpp
 * @brief elastic example with adaptive order refinement
 *
 * @copyright Copyright (c) 2026
 *
 */
 
#include <ElasticExample.hpp>
 
struct ElasticAdaptiveExample : ElasticExample {
 
  ElasticAdaptiveExample(MoFEM::Interface &m_field) : ElasticExample(m_field) {}
 
  MoFEMErrorCode runProblem();
 
private:
  MoFEMErrorCode readMesh();
  MoFEMErrorCode setupAdaptivity();
  MoFEMErrorCode resetPipelines();
  MoFEMErrorCode computeErrorNorms();
  MoFEMErrorCode outputResults(int ref_iter = -1);
  MoFEMErrorCode refineOrder(int ref_level);
  MoFEMErrorCode computeErrorIndicators();
  MoFEMErrorCode refineGeometry();
  MoFEMErrorCode calculateGeometryError();
  MoFEMErrorCode kspSetUpAndSolve(SmartPetscObj<KSP> solver) override;
 
  int oRder = 2;
  PetscBool refGeom = PETSC_FALSE;
  int refAdaptNum = 0;
  double refAdaptThreshold = 1;
  double refGeomThreshold = 1;
  std::array<double, 2> meanError = {0.0, 0.0};
  double errorIndicMean = 0.0;
 
  int atomTest = 0;
};
 
enum IndicInds { ERROR_TOT = 0, COUNT_ENTS, ERROR_NORM, LAST_IND };
 
#include <AdaptiveOrderRef.hpp>
#include <HighOrderGeometryRef.hpp>
 
//! [Run problem]
MoFEMErrorCode ElasticAdaptiveExample::runProblem() {
  MoFEMFunctionBegin;
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR setupAdaptivity();
  CHKERR boundaryCondition();
  CHKERR assembleSystem();
  CHKERR solveSystem();
  if (atomTest == 3)
    CHKERR computeErrorNorms();
  CHKERR checkResults();
 
  if (refAdaptNum > 0) {
    CHKERR computeErrorIndicators();
    CHKERR outputResults(0);
  } else {
    CHKERR outputResults();
  }
 
  for (int it = 0; it < refAdaptNum; ++it) {
    CHKERR refineOrder(it + 1);
    CHKERR resetPipelines();
    CHKERR assembleSystem();
    CHKERR solveSystem();
    CHKERR checkResults();
    if (it == refAdaptNum - 1 && atomTest == 4)
      CHKERR computeErrorNorms();
    CHKERR computeErrorIndicators();
    CHKERR outputResults(it + 1);
  }
  MoFEMFunctionReturn(0);
}
//! [Run problem]
 
//! [Read mesh]
MoFEMErrorCode ElasticAdaptiveExample::readMesh() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  CHKERR simple->getOptions();
 
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-ref_geom", &refGeom,
                             PETSC_NULLPTR);
 
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-ref_adapt_num", &refAdaptNum,
                            PETSC_NULLPTR);
 
  if ((refAdaptNum > 0) && refGeom) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Adaptive order and geometry refinement cannot be used together");
  }
 
  simple->getAddSkeletonFE() = true;
  CHKERR simple->loadFile();
 
  // Add meshsets if config file provided
  CHKERR mField.getInterface<MeshsetsManager>()->setMeshsetFromFile();
 
  auto update_ghost_ents = [&]() {
    MoFEMFunctionBegin;
 
    for (auto m : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
             std::regex((boost::format("%s(.*)") % "MAT_ELASTIC").str()))
 
    ) {
      Range ents;
      CHKERR mField.get_moab().get_entities_by_handle(m->getMeshset(), ents,
                                                      true);
      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(ents);
      CHKERR mField.get_moab().add_entities(m->getMeshset(), ents);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR update_ghost_ents();
 
  MoFEMFunctionReturn(0);
}
//! [Read mesh]
 
//! [Set up problem]
MoFEMErrorCode ElasticAdaptiveExample::setupAdaptivity() {
  MoFEMFunctionBegin;
 
  Range domainEntities;
  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
                                                     domainEntities);
  Tag th_order;
  CHKERR getTagHandle(mField, "ORDER", MB_TYPE_INTEGER, th_order);
 
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &oRder, PETSC_NULLPTR);
  for (auto ent : domainEntities) {
    CHKERR mField.get_moab().tag_set_data(th_order, &ent, 1, &oRder);
  }
 
  CHKERR copyTagOnSkin(mField, "ORDER", MB_TYPE_INTEGER);
 
  if (refGeom) {
    CHKERR calculateGeometryError();
    CHKERR refineGeometry(); // optional ho refinement
  }
 
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-test", &atomTest,
                            PETSC_NULLPTR);
 
  MoFEMFunctionReturn(0);
}
//! [Set up problem]
 
//! [Reset pipelines]
MoFEMErrorCode ElasticAdaptiveExample::resetPipelines() {
  MoFEMFunctionBegin;
  auto pip = mField.getInterface<PipelineManager>();
 
  pip->getDomainLhsFE().reset();
  pip->getDomainRhsFE().reset();
  pip->getBoundaryRhsFE().reset();
  pip->getBoundaryLhsFE().reset();
  pip->getSkeletonLhsFE().reset();
  pip->getSkeletonRhsFE().reset();
 
  MoFEMFunctionReturn(0);
}
//! [Reset pipelines]
 
// ! [KSP set up]
MoFEMErrorCode
ElasticAdaptiveExample::kspSetUpAndSolve(SmartPetscObj<KSP> solver) {
  MoFEMFunctionBegin;
 
  MOFEM_LOG_CHANNEL("TIMER");
  MOFEM_LOG_TAG("TIMER", "timer");
 
  DM dm;
  CHKERR KSPGetDM(solver, &dm);
  auto D = createDMVector(dm);
  auto F = vectorDuplicate(D);
 
  CHKERR VecZeroEntries(D);
  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
 
  BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
  MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp";
  CHKERR KSPSetUp(solver);
  MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp <= Done";
  MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve";
  CHKERR KSPSolve(solver, F, D);
  MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve <= Done";
 
  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
 
  MoFEMFunctionReturn(0);
}
// ! [KSP set up]
 
inline auto sqr = [](double v) { return v * v; };
 
inline MoFEM::VectorFunc analyticalDisplacement = [](double x, double y,
                                                     double z) {
  double E = 1e3;
  double nu = 0.3;
 
  double r = std::sqrt(x * x + y * y);
 
  double a = 0.5;
  double b = 1.0;
 
  double pa = 2;
  double pb = 1;
 
  double A = (pa * sqr(a) - pb * sqr(b)) / (sqr(b) - sqr(a));
  double B = sqr(a) * sqr(b) * (pa - pb) / (sqr(b) - sqr(a));
 
  VectorDouble u;
  u.resize(SPACE_DIM);
  std::vector<double> v_u;
 
  double u_r = (1. + nu) / E * ((1. - 2. * nu) * A * r + B / r);
 
  if (SPACE_DIM == 2)
    v_u = {u_r / r * x, u_r / r * y};
  else
    v_u = {u_r / r * x, u_r / r * y, 0.0};
 
  for (int i = 0; i < SPACE_DIM; ++i)
    u(i) = v_u[i];
 
  return u;
};
 
//! [Postprocess results]
MoFEMErrorCode ElasticAdaptiveExample::outputResults(int ref_iter) {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
 
  CHKERR resetPipelines();
 
  std::ostringstream strm;
  strm << "out_elastic";
  if (ref_iter >= 0) {
    strm << "_" << ref_iter;
  }
  strm << ".h5m";
 
  CHKERR ElasticPostProc::postProcessElasticResults<SPACE_DIM, DomainEle,
                                                    BoundaryEle, SideEle>(
      mField, simple->getDM(), simple->getDomainFEName(), strm.str(), {},
      {"GEOMETRY_ERROR", "SOLUTION_ERROR", "ORDER"}, Sev::verbose);
  MoFEMFunctionReturn(0);
}
//! [Postprocess results]
 
//! [Error Norms]
MoFEMErrorCode ElasticAdaptiveExample::computeErrorNorms() {
  MOFEM_LOG_CHANNEL("WORLD");
  auto pip = mField.getInterface<PipelineManager>();
  MoFEMFunctionBegin;
 
  CHKERR resetPipelines();
 
  auto integration_rule = [](int, int, int p_data) { return p_data + 1; };
  CHKERR pip->setDomainRhsIntegrationRule(integration_rule);
 
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      pip->getOpDomainRhsPipeline(), {H1}, "GEOMETRY");
 
  auto u_ptr = boost::make_shared<MatrixDouble>();
  pip->getOpDomainRhsPipeline().push_back(
      new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
 
  auto analytical_disp_ptr = boost::make_shared<MatrixDouble>();
 
  pip->getOpDomainRhsPipeline().push_back(
      new OpGetTensor1fromFunc<SPACE_DIM, SPACE_DIM>(analytical_disp_ptr,
                                                     analyticalDisplacement));
 
  auto normsVec = createVectorMPI(mField.get_comm(),
                                  (mField.get_comm_rank() == 0) ? 1 : 0, 1);
 
  pip->getOpDomainRhsPipeline().push_back(new OpCalcNormL2Tensor1<SPACE_DIM>(
      u_ptr, normsVec, 0, analytical_disp_ptr));
 
  CHKERR VecZeroEntries(normsVec);
  CHKERR pip->loopFiniteElements();
  CHKERR VecAssemblyBegin(normsVec);
  CHKERR VecAssemblyEnd(normsVec);
 
  if (mField.get_comm_rank() == 0) {
    const double *norms;
    CHKERR VecGetArrayRead(normsVec, &norms);
    MOFEM_TAG_AND_LOG("SELF", Sev::inform, "Errors")
        << "Displacement error L2 norm: " << std::scientific
        << std::sqrt(norms[0]);
 
    if (atomTest == 3 && std::sqrt(norms[0]) > 2.195e-06) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "atom test failed: displacement error L2 norm is too high: "
              "%3.6e",
              std::sqrt(norms[0]));
    }
 
    if (atomTest == 4 && std::sqrt(norms[0]) > 1.995e-07) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "atom test failed: displacement error L2 norm is too high: "
              "%3.6e",
              std::sqrt(norms[0]));
    }
 
    CHKERR VecRestoreArrayRead(normsVec, &norms);
  }
 
  MoFEMFunctionReturn(0);
}
//! [Error Norms]