AdaptiveOrderRef.hpp

Go to the documentation of this file.

/**
 * @file AdaptiveOrderRef.hpp
 * @brief Calculate error indicators
 * @date 2024-04-14
 *
 * @copyright Copyright (c) 2025
 *
 */
 
static MoFEMErrorCode getTagHandle(MoFEM::Interface &m_field, const char *name,
                                   DataType type, Tag &tag_handle) {
  MoFEMFunctionBegin;
  int int_val = 0;
  double double_val = 0;
  switch (type) {
  case MB_TYPE_INTEGER:
    CHKERR m_field.get_moab().tag_get_handle(
        name, 1, type, tag_handle, MB_TAG_CREAT | MB_TAG_SPARSE, &int_val);
    break;
  case MB_TYPE_DOUBLE:
    CHKERR m_field.get_moab().tag_get_handle(
        name, 1, type, tag_handle, MB_TAG_CREAT | MB_TAG_SPARSE, &double_val);
    break;
  default:
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
             "Wrong data type %d for tag", type);
  }
  MoFEMFunctionReturn(0);
}
 
static MoFEMErrorCode getDomainEntsPart(MoFEM::Interface &m_field,
                                        Range &domain_ents_part) {
  MoFEMFunctionBegin;
  auto simple = m_field.getInterface<Simple>();
  Range domain_ents;
  CHKERR m_field.get_moab().get_entities_by_dimension(simple->getMeshset(),
                                                      SPACE_DIM, domain_ents);
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
  CHKERR pcomm->filter_pstatus(domain_ents, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1,
                               &domain_ents_part);
  MoFEMFunctionReturn(0);
}
 
static MoFEMErrorCode copyTagOnSkin(MoFEM::Interface &m_field, const char *name,
                                    DataType type) {
  MoFEMFunctionBegin;
  auto simple = m_field.getInterface<Simple>();
 
  Range all_boundary_ents;
  CHKERR m_field.get_moab().get_entities_by_handle(simple->getBoundaryMeshSet(),
                                                   all_boundary_ents);
 
  Range domain_ents_part;
  CHKERR getDomainEntsPart(m_field, domain_ents_part);
  Range domain_faces, boundary_faces;
  CHKERR m_field.get_moab().get_adjacencies(domain_ents_part, SPACE_DIM - 1,
                                            true, domain_faces,
                                            moab::Interface::UNION);
  boundary_faces = intersect(domain_faces, all_boundary_ents);
 
  Tag th;
  CHKERR getTagHandle(m_field, name, type, th);
 
  for (auto &face_ent : boundary_faces) {
    Range adj_tet;
    CHKERR m_field.get_moab().get_adjacencies(&face_ent, 1, SPACE_DIM, true,
                                              adj_tet, moab::Interface::UNION);
    if (adj_tet.size() != 1)
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "There should be only one ent adjacent to the boundary ent");
 
    std::variant<int, double> data;
 
    switch (type) {
    case MB_TYPE_INTEGER:
      data = int{};
      break;
    case MB_TYPE_DOUBLE:
      data = double{};
      break;
    default:
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
               "Wrong data type %d for tag", type);
    }
 
    std::visit(
        [&](auto &val) {
          CHKERR m_field.get_moab().tag_get_data(th, adj_tet, &val);
          CHKERR m_field.get_moab().tag_set_data(th, &face_ent, 1, &val);
        },
        data);
  }
 
  MoFEMFunctionReturn(0);
}
 
struct OpCalculateTractionJump : public SideEleOp {
  OpCalculateTractionJump(const std::string field_name,
                          boost::shared_ptr<MatrixDouble> stress_ptr,
                          boost::shared_ptr<MatrixDouble> traction_ptr)
      : SideEleOp(field_name, SideEleOp::OPROW, false), stressPtr(stress_ptr),
        tractionJumpPtr(traction_ptr) {
    std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
    doEntities[MBVERTEX] = true;
  }
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
    MoFEMFunctionBegin;
    FTensor::Index<'i', SPACE_DIM> i;
    FTensor::Index<'j', SPACE_DIM> j;
    auto nb_gauss_pts = getGaussPts().size2();
 
    // auto n_total_side_eles = getLoopSize();
    auto n_InLoop = getNinTheLoop();
    auto e_sense = getSkeletonSense();
 
    if (n_InLoop == 0) {
      tractionJumpPtr->resize(nb_gauss_pts, SPACE_DIM, false);
      tractionJumpPtr->clear();
    }
 
    auto t_normal = getFTensor1NormalsAtGaussPts();
    auto t_stress = getFTensor2SymmetricFromMat<SPACE_DIM>(*stressPtr);
    auto t_traction = getFTensor1FromMat<SPACE_DIM>(*tractionJumpPtr);
 
    for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
      FTensor::Tensor1<double, 3> t_unit_normal;
      t_unit_normal(i) = t_normal(i);
      t_unit_normal.normalize();
 
      t_traction(i) -= e_sense * (t_stress(i, j) * t_unit_normal(j));
 
      ++t_normal;
      ++t_stress;
      ++t_traction;
    }
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<MatrixDouble> stressPtr;
  boost::shared_ptr<MatrixDouble> tractionJumpPtr;
};
 
MoFEMErrorCode ElasticAdaptiveExample::refineOrder(int ref_level) {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto comm = mField.getInterface<CommInterface>();
 
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, "", "-ref_adapt_threshold",
                             &refAdaptThreshold, PETSC_NULLPTR);
 
  Tag th_error_ind, th_order;
  CHKERR getTagHandle(mField, "SOLUTION_ERROR", MB_TYPE_DOUBLE, th_error_ind);
  CHKERR getTagHandle(mField, "ORDER", MB_TYPE_INTEGER, th_order);
 
  std::vector<Range> refinement_levels(ref_level + 2);
 
  Range domain_ents_part;
  CHKERR getDomainEntsPart(mField, domain_ents_part);
 
  for (auto ent : domain_ents_part) {
    double err_indic = 0;
    CHKERR mField.get_moab().tag_get_data(th_error_ind, &ent, 1, &err_indic);
    int order;
    CHKERR mField.get_moab().tag_get_data(th_order, &ent, 1, &order);
 
    if (err_indic > refAdaptThreshold * errorIndicMean) {
      order += 1;
      refinement_levels[order - oRder].insert(ent);
    }
  }
 
  for (int ll = refinement_levels.size() - 1; ll > 1; ll--) {
    CHKERR comm->synchroniseEntities(refinement_levels[ll]);
 
    Range face_ents;
    CHKERR mField.get_moab().get_adjacencies(refinement_levels[ll],
                                             SPACE_DIM - 1, true, face_ents,
                                             moab::Interface::UNION);
    Range tets;
    CHKERR mField.get_moab().get_adjacencies(face_ents, SPACE_DIM, true, tets,
                                             moab::Interface::UNION);
    tets = subtract(tets, refinement_levels[ll]);
    refinement_levels[ll - 1].merge(tets);
  }
 
  for (int ll = 1; ll < refinement_levels.size(); ll++) {
 
    if (oRder + ll > 8) {
      MOFEM_LOG("WORLD", Sev::warning)
          << "setting approximation order higher than 8 is not currently "
             "supported"
          << endl;
    } else {
 
      for (auto ent : refinement_levels[ll]) {
        int order = oRder + ll;
        CHKERR mField.get_moab().tag_set_data(th_order, &ent, 1, &order);
      }
 
      Range ents_to_refine;
      ents_to_refine.merge(refinement_levels[ll]);
      for (int d = 1; d < SPACE_DIM; ++d) {
        Range lower_dim_ents;
        CHKERR mField.get_moab().get_adjacencies(refinement_levels[ll], d, true,
                                                 lower_dim_ents,
                                                 moab::Interface::UNION);
        ents_to_refine.merge(lower_dim_ents);
      }
 
      CHKERR comm->synchroniseEntities(ents_to_refine);
 
      CHKERR simple->setFieldOrder("U", oRder + ll, &ents_to_refine);
    }
  }
 
  CHKERR copyTagOnSkin(mField, "ORDER", MB_TYPE_INTEGER);
 
  CHKERR comm->synchroniseFieldEntities("U");
  CHKERR simple->reSetUp();
  
  CHKERR boundaryCondition();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ElasticAdaptiveExample::computeErrorIndicators() {
  MoFEMFunctionBegin;
  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
  auto simple = mField.getInterface<Simple>();
 
  Tag th_error_ind;
  CHKERR getTagHandle(mField, "SOLUTION_ERROR", MB_TYPE_DOUBLE, th_error_ind);
 
  auto integration_rule = [](int, int, int p_data) { return 2 * p_data + 1; };
  CHKERR pipeline_mng->setSkeletonRhsIntegrationRule(integration_rule);
  auto &skeleton_pipeline = pipeline_mng->getOpSkeletonRhsPipeline();
 
  if (skeleton_pipeline.empty()) {
 
    auto op_loop_side =
        new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
 
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        op_loop_side->getOpPtrVector(), {H1}, "GEOMETRY");
 
    // Use HookeOps commonDataFactory to get matStrainPtr and matCauchyStressPtr
    auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEle>(
        mField, op_loop_side->getOpPtrVector(), "U", "MAT_ELASTIC",
        Sev::verbose);
 
    auto mat_jump_ptr = boost::make_shared<MatrixDouble>();
    op_loop_side->getOpPtrVector().push_back(new OpCalculateTractionJump(
        "U", common_ptr->getMatCauchyStress(), mat_jump_ptr));
 
    auto *op_ptr = new BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE);
    op_ptr->doWorkRhsHook = [this, th_error_ind,
                             mat_jump_ptr](DataOperator *base_op_ptr, int side,
                                           EntityType type,
                                           EntitiesFieldData::EntData &data) {
      MoFEMFunctionBegin;
      FTensor::Index<'i', SPACE_DIM> i;
      auto op_ptr = static_cast<BoundaryEleOp *>(base_op_ptr);
      auto fe_ent = op_ptr->getFEEntityHandle();
 
      auto t_jump = getFTensor1FromMat<SPACE_DIM>(*mat_jump_ptr);
      auto t_w = op_ptr->getFTensor0IntegrationWeight();
      double err_int = 0;
      const auto nb_gauss_pts = op_ptr->getGaussPts().size2();
      for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
        err_int += t_w * (t_jump(i) * t_jump(i));
        ++t_jump;
        ++t_w;
      }
      double err_measure = sqrt(err_int);
      CHKERR mField.get_moab().tag_set_data(th_error_ind, &fe_ent, 1,
                                            &err_measure);
 
      MoFEMFunctionReturn(0);
    };
 
    skeleton_pipeline.push_back(op_loop_side);
    skeleton_pipeline.push_back(op_ptr);
  }
 
  Range fe_ents;
  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM - 1, fe_ents);
  double zero = 0;
  CHKERR mField.get_moab().tag_clear_data(th_error_ind, fe_ents, &zero);
 
  CHKERR DMoFEMLoopFiniteElements(simple->getDM(), simple->getSkeletonFEName(),
                                  pipeline_mng->getSkeletonRhsFE());
 
  Range domain_ents_part;
  CHKERR getDomainEntsPart(mField, domain_ents_part);
 
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
  Range empty_range;
  CHKERR pcomm->reduce_tags(th_error_ind, MPI_SUM, empty_range);
 
  std::array<double, 3> error_indic_loc = {0.0, 0.0, 0.0};
 
  for (auto domain_ent : domain_ents_part) {
    Range face_ents;
    CHKERR mField.get_moab().get_adjacencies(
        &domain_ent, 1, SPACE_DIM - 1, true, face_ents, moab::Interface::UNION);
 
    double err_indic_sum = 0;
    for (auto face_ent : face_ents) {
      double err_indic = 0;
      CHKERR mField.get_moab().tag_get_data(th_error_ind, &face_ent, 1,
                                            &err_indic);
      err_indic_sum += err_indic;
    }
    CHKERR mField.get_moab().tag_set_data(th_error_ind, &domain_ent, 1,
                                          &err_indic_sum);
 
    error_indic_loc[IndicInds::ERROR_TOT] += err_indic_sum;
 
    const EntityHandle *vert_conn;
    int vert_num;
    CHKERR mField.get_moab().get_connectivity(domain_ent, vert_conn, vert_num,
                                              true);
    std::vector<double> vpos(3 * vert_num);
    CHKERR mField.get_moab().get_coords(vert_conn, vert_num, &vpos[0]);
    double vol = Tools::tetVolume(vpos.data());
 
    error_indic_loc[IndicInds::ERROR_NORM] += err_indic_sum * vol;
  }
 
  CHKERR copyTagOnSkin(mField, "SOLUTION_ERROR", MB_TYPE_DOUBLE);
 
  error_indic_loc[IndicInds::COUNT_ENTS] =
      static_cast<double>(domain_ents_part.size());
 
  std::array<double, 3> error_indic_glob = {0.0, 0.0, 0.0};
 
  MPI_Allreduce(&error_indic_loc[0], &error_indic_glob[0], 3, MPI_DOUBLE,
                MPI_SUM, mField.get_comm());
 
  errorIndicMean = error_indic_glob[IndicInds::ERROR_TOT] /
                   error_indic_glob[IndicInds::COUNT_ENTS];
 
  MOFEM_LOG("WORLD", Sev::inform) << "Global error indicator (norm): "
                                  << error_indic_glob[IndicInds::ERROR_NORM];
  MOFEM_LOG("WORLD", Sev::inform)
      << "Global error indicator (mean): " << errorIndicMean;
 
  MoFEMFunctionReturn(0);
}