entities_petsc_vector.cpp

Go to the documentation of this file.

/** \file entities_petsc_vector.cpp
 * \brief Exchange vectors with imported partition labels in serial and parallel.
 */
#include <MoFEM.hpp>
using namespace MoFEM;
 
namespace {
 
MoFEMErrorCode checkExchange(moab::Interface &moab, ParallelComm &pcomm,
                             bool include_vertices) {
  MoFEMFunctionBegin;
  constexpr int nb_coeffs = 2;
  auto exchange = CommInterface::createEntitiesPetscVector(
      PETSC_COMM_WORLD, moab,
      [include_vertices](Range ents) {
        auto selected = ents.subset_by_type(MBTRI);
        if (include_vertices)
          selected.merge(ents.subset_by_type(MBVERTEX));
        return selected;
      },
      nb_coeffs, Sev::inform, 0, include_vertices);
  const auto &owned = exchange.first.first;
  const auto &ghosts = exchange.first.second;
  const int entities_per_part = include_vertices ? 4 : 1;
  const int expected_owned = entities_per_part * (pcomm.size() == 1 ? 2 : 1);
  const int expected_ghosts = pcomm.size() == 1 ? 0 : entities_per_part;
  if (owned.size() != expected_owned || ghosts.size() != expected_ghosts)
    SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
            "Unexpected owned/ghost entity counts: %d/%d, expected %d/%d",
            static_cast<int>(owned.size()), static_cast<int>(ghosts.size()),
            expected_owned, expected_ghosts);
  PetscInt local_size, global_size;
  CHKERR VecGetLocalSize(exchange.second, &local_size);
  CHKERR VecGetSize(exchange.second, &global_size);
  if (local_size != nb_coeffs * expected_owned ||
      global_size != nb_coeffs * entities_per_part * 2)
    SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
            "Exchange vector omits mesh entities");
 
  const auto gid_tag = moab.globalId_tag();
  std::vector<int> gids(owned.size());
  CHKERR moab.tag_get_data(gid_tag, owned, gids.data());
  PetscScalar *values;
  CHKERR VecGetArray(exchange.second, &values);
  for (size_t e = 0; e != gids.size(); ++e)
    for (int c = 0; c != nb_coeffs; ++c)
      values[nb_coeffs * e + c] = nb_coeffs * gids[e] + c;
  CHKERR VecRestoreArray(exchange.second, &values);
  CHKERR VecGhostUpdateBegin(exchange.second, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(exchange.second, INSERT_VALUES, SCATTER_FORWARD);
 
  // Check entity identity as well as vector sizes: a valid scatter can still
  // exchange values for the wrong entity if the GLOBAL_ID mapping is incorrect.
  gids.resize(ghosts.size());
  CHKERR moab.tag_get_data(gid_tag, ghosts, gids.data());
  Vec local_form;
  const PetscScalar *local_values;
  CHKERR VecGhostGetLocalForm(exchange.second, &local_form);
  CHKERR VecGetArrayRead(local_form, &local_values);
  bool valid = true;
  for (size_t e = 0; e != gids.size(); ++e)
    for (int c = 0; c != nb_coeffs; ++c)
      valid = valid && local_values[local_size + nb_coeffs * e + c] ==
                           nb_coeffs * gids[e] + c;
  CHKERR VecRestoreArrayRead(local_form, &local_values);
  CHKERR VecGhostRestoreLocalForm(exchange.second, &local_form);
  if (!valid)
    SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
            "Ghost values do not match the owning entities");
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkSerialFiniteElements(moab::Interface &moab) {
  MoFEMFunctionBegin;
  Range faces, edges;
  CHKERR moab.get_entities_by_type(0, MBTRI, faces);
  CHKERR moab.get_adjacencies(faces, 1, true, edges, moab::Interface::UNION);
  int gid = 9;
  for (auto edge : edges) {
    CHKERR moab.tag_set_data(moab.globalId_tag(), &edge, 1, &gid);
    ++gid;
  }
  MoFEM::Core core(moab);
  MoFEM::Interface &m_field = core;
  const auto level = BitRefLevel().set(0);
  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(0, 2,
                                                                    level);
  CHKERR m_field.add_field("TEST_FIELD", L2, AINSWORTH_LEGENDRE_BASE, 1);
  CHKERR m_field.add_ents_to_field_by_type(0, MBTRI, "TEST_FIELD");
  CHKERR m_field.set_field_order(0, MBTRI, "TEST_FIELD", 0);
  CHKERR m_field.add_finite_element("TEST_FE");
  CHKERR m_field.modify_finite_element_add_field_row("TEST_FE", "TEST_FIELD");
  CHKERR m_field.modify_finite_element_add_field_col("TEST_FE", "TEST_FIELD");
  CHKERR m_field.modify_finite_element_add_field_data("TEST_FE", "TEST_FIELD");
  CHKERR m_field.add_ents_to_finite_element_by_type(0, MBTRI, "TEST_FE");
  CHKERR m_field.add_problem("TEST_PROBLEM");
  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "TEST_FE");
  CHKERR m_field.modify_problem_ref_level_add_bit("TEST_PROBLEM", level);
  CHKERR m_field.build_fields();
  CHKERR m_field.build_finite_elements();
  CHKERR m_field.build_adjacencies(level);
  auto *problems = m_field.getInterface<ProblemsManager>();
  CHKERR problems->buildProblem("TEST_PROBLEM", true);
  CHKERR problems->partitionSimpleProblem("TEST_PROBLEM");
  CHKERR problems->partitionFiniteElements("TEST_PROBLEM", true);
  const Problem *problem;
  CHKERR m_field.get_problem("TEST_PROBLEM", &problem);
  if (problem->numeredFiniteElementsPtr->size() != 2)
    SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
            "Serial problem lost a finite element from the imported mesh");
  for (const auto &fe : *problem->numeredFiniteElementsPtr)
    if (fe->getPart() != 0)
      SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
              "Serial finite element still belongs to imported partition %d",
              static_cast<int>(fe->getPart()));
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode runTest(moab::Interface &moab, ParallelComm &pcomm) {
  MoFEMFunctionBegin;
  if (pcomm.size() > 2)
    SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
            "This fixture requires one or two MPI ranks");
  const auto gid_tag = moab.globalId_tag();
  const auto part_tag = pcomm.part_tag();
  // Two disconnected triangles replicate a mesh imported with two partitions.
  // Keep both labels even in serial to exercise the reported failure.
  for (int part = 0; part != 2; ++part) {
    const double x = 2. * part;
    const double coords[] = {x, 0., 0., x + 1., 0., 0., x, 1., 0.};
    Range vertices;
    CHKERR moab.create_vertices(coords, 3, vertices);
    const EntityHandle conn[] = {vertices[0], vertices[1], vertices[2]};
    EntityHandle face, meshset;
    CHKERR moab.create_element(MBTRI, conn, 3, face);
    CHKERR moab.create_meshset(MESHSET_SET, meshset);
    CHKERR moab.add_entities(meshset, &face, 1);
    CHKERR moab.tag_set_data(part_tag, &meshset, 1, &part);
    CHKERR moab.tag_set_data(pcomm.partition_tag(), &face, 1, &part);
    const int face_gid = part + 1;
    const int vertex_gids[] = {3 + 3 * part, 4 + 3 * part, 5 + 3 * part};
    CHKERR moab.tag_set_data(gid_tag, &face, 1, &face_gid);
    CHKERR moab.tag_set_data(gid_tag, vertices, vertex_gids);
  }
  CHKERR checkExchange(moab, pcomm, false);
  CHKERR checkExchange(moab, pcomm, true);
  if (pcomm.size() == 1)
    CHKERR checkSerialFiniteElements(moab);
  MoFEMFunctionReturn(0);
}
 
} // namespace
 
int main(int argc, char *argv[]) {
  MoFEM::Core::Initialize(&argc, &argv, nullptr,
                          "Entity exchange with imported partition labels\n");
  try {
    moab::Core moab_core;
    ParallelComm pcomm(&moab_core, PETSC_COMM_WORLD);
    CHKERR runTest(moab_core, pcomm);
  }
  CATCH_ERRORS;
  CHKERR MoFEM::Core::Finalize();
  return 0;
}