HighOrderGeometryRef.hpp

Go to the documentation of this file.

/**
 * @file HighOrderGeometryRef.hpp
 * @brief Higher Order Geometry refinement
 * @date 2025-05-1
 * 
 * @copyright Copyright (c) 2025
 * 
 */
 
/**
 * @brief Calculate geometry error and save on a tag
 */
 
MoFEMErrorCode ElasticAdaptiveExample::calculateGeometryError() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
 
  double def_val_tag = 0.0;
  Tag tag_err;
  CHKERR mField.get_moab().tag_get_handle("GEOMETRY_ERROR", 1, MB_TYPE_DOUBLE,
                                          tag_err, MB_TAG_CREAT | MB_TAG_SPARSE,
                                          &def_val_tag);
 
  auto LoCoordsPtr = boost::make_shared<MatrixDouble>();
  auto *op_ptr = new BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE);
  op_ptr->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
                              EntityType type,
                              EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
    auto op_ptr = static_cast<BoundaryEleOp *>(base_op_ptr);
    *LoCoordsPtr = op_ptr->getCoordsAtGaussPts();
    MoFEMFunctionReturn(0);
  };
 
  auto *op_ptr_ho = new BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE);
  op_ptr_ho->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
                                 EntityType type,
                                 EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
    auto op_ptr_ho = static_cast<BoundaryEleOp *>(base_op_ptr);
    auto fe_ent = op_ptr_ho->getFEEntityHandle();
    FTensor::Index<'i', SPACE_DIM> i;
    MatrixDouble &lo_mat = *LoCoordsPtr;
    MatrixDouble &ho_mat = op_ptr_ho->getCoordsAtGaussPts();
    auto t_lo_coords = getFTensor1FromMat<SPACE_DIM>(lo_mat);
    auto t_ho_coords = getFTensor1FromMat<SPACE_DIM>(ho_mat);
    auto t_w = op_ptr_ho->getFTensor0IntegrationWeight();
    double err_int = 0;
    const auto nb_gauss_pts = op_ptr_ho->getGaussPts().size2();
    // E = sqrt( int (x-xh)*(x-xh) da )  / face_area
    for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
      err_int += t_w * ((t_lo_coords(i) - t_ho_coords(i)) *
                        (t_lo_coords(i) - t_ho_coords(i)));
      ++t_lo_coords;
      ++t_ho_coords;
      ++t_w;
    }
    const double err_measure = sqrt(err_int) / op_ptr_ho->getMeasure();
    if (err_measure > std::numeric_limits<float>::epsilon()) {
      meanError[IndicInds::ERROR_TOT] += err_measure;
      meanError[IndicInds::COUNT_ENTS] += 1.0;
    }
 
    CHKERR mField.get_moab().tag_set_data(tag_err, &fe_ent, 1, &err_measure);
    Range adj_tet;
    CHKERR mField.get_moab().get_adjacencies(&fe_ent, 1, SPACE_DIM, false,
                                             adj_tet);
    if (adj_tet.size() != 1)
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "There should be only one ent adjacent to the boundary ent");
    auto tet_ent = *adj_tet.begin();
    double tet_err_measure = 0;
    CHKERR mField.get_moab().tag_get_data(tag_err, &tet_ent, 1,
                                          &tet_err_measure);
    if (err_measure > tet_err_measure) {
      CHKERR mField.get_moab().tag_set_data(tag_err, &tet_ent, 1, &err_measure);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto fe = boost::make_shared<BoundaryEle>(mField);
  fe->getOpPtrVector().push_back(op_ptr);
 
  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      fe->getOpPtrVector(), {NOSPACE}, "GEOMETRY");
  fe->getOpPtrVector().push_back(op_ptr_ho);
  auto dm = simple->getDM();
  CHK_THROW_MESSAGE(
      DMoFEMLoopFiniteElements(dm, simple->getBoundaryFEName(), fe),
      "calc metric");
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ElasticAdaptiveExample::refineGeometry() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
 
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, "", "-ref_geom_threshold",
                             &refGeomThreshold, PETSC_NULLPTR);
 
  if (Tag tag_err = 0; mField.get_moab().tag_get_handle(
                           "GEOMETRY_ERROR", tag_err) == MB_SUCCESS) {
 
    std::array<double, 2> mean_error_glob = {0.0, 0.0};
    MPI_Allreduce(&meanError[0], &mean_error_glob[0], 2, MPI_DOUBLE, MPI_SUM,
                  mField.get_comm());
    if (mean_error_glob[IndicInds::COUNT_ENTS] > 0) {
      refGeomThreshold *= (mean_error_glob[IndicInds::ERROR_TOT] /
                           mean_error_glob[IndicInds::COUNT_ENTS]);
    }
 
    MOFEM_LOG("WORLD", Sev::inform)
        << "Geometry refinement threshold = " << refGeomThreshold;
 
    auto get_ents_to_refine = [&]() {
      Range outer_ref_ents;
 
      auto *op_ptr_ref = new BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE);
      op_ptr_ref->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
                                      EntityType type,
                                      EntitiesFieldData::EntData &data) {
        MoFEMFunctionBegin;
        auto op_ptr_ref = static_cast<BoundaryEleOp *>(base_op_ptr);
        auto fe_ent = op_ptr_ref->getFEEntityHandle();
        double err_measure = 0;
        CHKERR mField.get_moab().tag_get_data(tag_err, &fe_ent, 1,
                                              &err_measure);
        if (err_measure > refGeomThreshold)
          outer_ref_ents.insert(fe_ent);
 
        MoFEMFunctionReturn(0);
      };
 
      auto fe = boost::make_shared<BoundaryEle>(mField);
      fe->getOpPtrVector().push_back(op_ptr_ref);
      auto dm = simple->getDM();
      CHK_THROW_MESSAGE(
          DMoFEMLoopFiniteElements(dm, simple->getBoundaryFEName(), fe),
          "ref ents");
      return outer_ref_ents;
    };
 
    auto ents_to_ref = get_ents_to_refine();
 
    Range all_domain_ents;
    CHKERR mField.get_moab().get_entities_by_handle(simple->getMeshset(),
                                                    all_domain_ents);
    Range all_boundary_ents;
    CHKERR mField.get_moab().get_entities_by_handle(
        simple->getBoundaryMeshSet(), all_boundary_ents);
 
    Range tets_to_refine[2];
    CHKERR mField.get_moab().get_adjacencies(ents_to_ref, SPACE_DIM, false,
                                             tets_to_refine[1],
                                             moab::Interface::UNION);
    Range ent_conn;
    CHKERR mField.get_moab().get_adjacencies(tets_to_refine[1], SPACE_DIM - 1,
                                             true, ent_conn,
                                             moab::Interface::UNION);
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(ent_conn);
 
    CHKERR mField.get_moab().get_adjacencies(
        ent_conn, SPACE_DIM, true, tets_to_refine[0], moab::Interface::UNION);
    tets_to_refine[0] = subtract(tets_to_refine[0], tets_to_refine[1]);
 
#ifndef NDEBUG
    CHKERR CutMeshInterface::SaveData(mField.get_moab())(
        "ents_to_refine_ho_2.vtk", tets_to_refine[0]);
    CHKERR CutMeshInterface::SaveData(mField.get_moab())(
        "ents_to_refine_ho_3.vtk", tets_to_refine[1]);
#endif
 
    for (int i = 0; i < 2; ++i) {
      Range ents_to_refine;
      ents_to_refine.merge(tets_to_refine[i]);
      for (int d = 1; d < SPACE_DIM; ++d) {
        Range lower_dim_ents;
        CHKERR mField.get_moab().get_adjacencies(
            tets_to_refine[i], d, true, lower_dim_ents, moab::Interface::UNION);
        ents_to_refine.merge(lower_dim_ents);
      }
      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
          ents_to_refine);
 
      CHKERR simple->setFieldOrder("U", oRder + i + 1, &ents_to_refine);
    }
 
    CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities("U");
    CHKERR simple->reSetUp();
 
    // save order on tags (volumes and boundary)
    int orders[] = {oRder, oRder + 1, oRder + 2};
    Tag tag_order;
 
    CHKERR getTagHandle(mField, "ORDER", MB_TYPE_INTEGER, tag_order);
 
    for (int i = 0; i < 2; ++i) {
      Range adj_faces;
      CHKERR mField.get_moab().get_adjacencies(tets_to_refine[i], SPACE_DIM - 1,
                                               true, adj_faces,
                                               moab::Interface::UNION);
      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
          adj_faces);
 
      adj_faces = intersect(adj_faces, all_boundary_ents);
      CHKERR mField.get_moab().tag_clear_data(tag_order, tets_to_refine[i],
                                              &orders[i + 1]);
      CHKERR mField.get_moab().tag_clear_data(tag_order, adj_faces,
                                              &orders[i + 1]);
    }
  }
  MoFEMFunctionReturn(0);
}