parallel neighbor mapping

[starkekr]

I am looking to develop a parallel neighbor mapping using the existing functionality in mfem. I have a phasefield code that I am trying to implement a grain tracker for. To do this I have to identify features (grains) in the mesh. Right now I have a code which correctly identifies features in the nonconformal mesh. These feature contain a list of local element ids. I have the code working in serial but have some problems in parallel. In parallel I can correctly identify the features (create a list of element ids), but I am struggling to sync these feature together into a global feature because I need to create a map to connect elements across processor boundaries. I have attached an image to show an example of the problem

On the left is the grain structure (showing the boundaries) and the right shows the grain ids which are incorrectly segmented across parallel boundaries.

[v-dobrev]

Here is one potential option: if your mapping (in serial) is feature_idx -> list of element_idx, you can try switch to the transpose mapping element_idx -> feature_idx which can be stored as an order 0 (p=0) L2 GridFunction. In parallel this becomes a ParGridFunction that can then be "updated" when the mesh is modified similar to other fields to maintain the element_idx -> feature_idx mapping.

[starkekr]

@paul-hilscher Is this algorithm implemented in mfem somewhere? I am more interested in, like you said, phase 3, particularly interested in forming the 1=2 && 1=4 connection in figure 3. I was going about this by using shared faces to form this connection. Can you expand on this part please?

I guess I dont really understand how using L2 order 0 as suggested above helps here.

Maybe shared faces are the wrong data object here. I will also be using mesh refinement so maybe the ghost elements will be needed here as well.

[paul-hilscher]

@starkekr ,

I will also be using mesh refinement so maybe the ghost elements will be needed here as well.

yes, ghost elements are required. L2/Order 0 provides you functionality to exchange those ghost values across processor boundaries.

Is this algorithm implemented in mfem somewhere?

No it is part of my private project. You can find the code below. Note:

• it uses some internal functions but it should be easy to replace this
• I have not tested below code for NC Meshes (it should still work tough assuming the loop over shared faces remains valid)
• Consider below code same license as mfem

main function

  mfem::L2_FECollection       fec(0, mesh.getParMesh().Dimension());
  mfem::ParFiniteElementSpace fespace(&mesh.getParMesh(), &fec);

  mfem::ParGridFunction gf(&fespace);

  // gf of type L2 order 0 has values only in the "cell center"
  gf = ... // Data

  assign_connected_components(gf);

ConnectedComponents.hpp

#include "mfem.hpp"

//! @brief A class implementing the Union-Find (Disjoint Set Union) data structure.
//!
//! @details This class provides the necessary functions to perform union and find operations on disjoint sets.
//!
//! if grid_function values with 0.0 are considered as background, and values > 0.0 are considered as foreground.
//!
//! Based on "Harrison, C., 2011. Data-parallel mesh connected components labeling and analysis."
//! Note that we support exchanges of halos which simplifies it significantly.
//!
void assign_connected_components(mfem::ParGridFunction &grid_function);

ConnectedComponents.cpp

#include <cstddef>
#include <functional>
#include <initializer_list>
#include <iterator>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>

#include "Common/Parallel/Parallel.hpp"
#include "DisjointSet.hpp"
#include "mfem.hpp"
#include "mpi.h"

namespace
{

  // Custom hash function for std::pair<int, int>
  struct PairHash {
    template<typename T1, typename T2>
    std::size_t operator()(std::pair<T1, T2> const &pair) const
    {
      // NOLINTBEGIN(hicpp-signed-bitwise)
      return std::hash<T1>()(pair.first) ^ (std::hash<T2>()(pair.second) << 1);
      // NOLINTEND(hicpp-signed-bitwise)
    }
  };

  using LabelPairSets = std::unordered_set<std::pair<int, int>, PairHash>;

  // Marks the local connected components
  int phase_1_local_components(mfem::ParMesh &par_mesh, mfem::ParGridFunction &grid_function,
                               std::function<bool(double)> const &is_label)
  {
    int current_component = -1;

    int const num_elements = par_mesh.GetNE();

    // Visited array to track processed elements
    std::vector<bool> visited(num_elements, false);

    // Queue for BFS
    std::vector<int> element_queue;

    mfem::Array<int> faces, orientations;

    for(int i = 0; i < num_elements; ++i) {
      // Skip already processed elements
      if(visited[i]) {
        continue;
      }

      // Skip elements not part of a sphere (sphere_field <= 0)
      if(!is_label(grid_function[i])) {
        grid_function[i] = -1;
        continue;
      }

      // Start a new component
      current_component++;

      element_queue.push_back(i);

      while(!element_queue.empty()) {
        int const j = element_queue.back();
        element_queue.pop_back();

        if(visited[j]) {
          continue;
        }

        visited[j] = true;

        grid_function[j] = current_component;

        // Get the faces of the current element
        par_mesh.GetElementFaces(j, faces, orientations);

        // Process neighbors through shared faces
        for(int const f: faces) {
          // Get elements sharing this face
          int elem1 = -1, elem2 = -1;
          par_mesh.GetFaceElements(f, &elem1, &elem2);

          // Check both elements adjacent to the face
          for(int const k: {elem1, elem2}) {
            if(k < 0 || k == j || visited[k]) {
              continue;
            }

            if(is_label(grid_function[k])) {
              element_queue.push_back(k);
            }
          }
        }
      }
    }

    return current_component + 1;
  }

  // Step 2: Re-number the components globally
  // We collect all labels and ensure that each label is globally unique
  //
  int phase_2_relabel_for_cross_processor_comparison(mfem::ParMesh &par_mesh, mfem::ParGridFunction &grid_function,
                                                     int num_local_components)
  {
    int rank_offset = 0;

    // MPI_Exscan computes the exclusive scan (also known as a prefix sum).
    // For a process i, the result (rank_offset) is the sum of all
    // num_local_components values from processes with ranks 0 to  i-1.
    MPI_Exscan(&num_local_components, &rank_offset, 1, MPI_INT, MPI_SUM, par_mesh.GetComm());

    // Apply the global offset to local components
    for(int i = 0; i < par_mesh.GetNE(); ++i) {
      if(grid_function[i] >= 0) {
        grid_function[i] += rank_offset;
      }
    }

    return rank_offset;
  }

  // We use ghost cells here to exchange the local union.
  std::unordered_set<std::pair<int, int>, PairHash> phase_3_prime_merge_labels_across_processors(
    mfem::ParMesh &par_mesh, mfem::ParGridFunction &components,

    int local_number_components, int rank_offset)

  {
    std::unordered_set<std::pair<int, int>, PairHash> local_unions_set;

    for(int i = 0; i < local_number_components; i++) {
      local_unions_set.insert({rank_offset + i, rank_offset + i});
    }

    par_mesh.ExchangeFaceNbrData();
    components.ExchangeFaceNbrData();

    auto const &components_ghost = components.FaceNbrData();

    // We do not need to look over local faces as those have been addressed
    // in the first step but now we iterative over faces and merge the labels between the primary and
    // ghost cells
    int const num_shared_faces = par_mesh.GetNSharedFaces();

    for(int sf = 0; sf < num_shared_faces; sf++) {
      int const face_id = par_mesh.GetSharedFace(sf);

      mfem::FaceElementTransformations const *ftr = par_mesh.GetFaceElementTransformations(face_id);

      if(ftr != nullptr) {
        // If you have a negative element index el_no, then the "face neighbor
        // index" is given by nbr_el_no = -1 - el_no. The "shifted face neighbor
        // index" is given by pmesh.GetNE() + nbr_el_no. This should work in
        // GetValue.
        // @todo: what is the "shifted face neighbor index"?

        int const local_cell_index = ftr->Elem1No;  // Primary cell index
        int const ghost_cell_index = ftr->Elem2No;  // Ghost cell index (might be -1 if
                                                    // no ghost element exists)

        int const local_cell_label = static_cast<int>(components[local_cell_index]);
        int const ghost_cell_label = static_cast<int>(components_ghost[-1 - ghost_cell_index]);

        if((local_cell_label != -1) && (ghost_cell_label != -1)) {
          local_unions_set.insert({local_cell_label, ghost_cell_label});
        }
      }
    }

    return local_unions_set;
  }

  void phase_4_final_assignment_of_labels(mfem::ParGridFunction                                   &grid_function,
                                          std::unordered_set<std::pair<int, int>, PairHash> const &local_unions)
  {
    // Step 1: Gather all unions from all processes
    int const size = Parallel::GetNumProcesses();

    // Serialize local unions into a vector of pairs
    std::vector<std::pair<int, int>> local_union_list(local_unions.begin(), local_unions.end());

    int local_size = local_union_list.size();

    std::vector<int> union_sizes(size);
    MPI_Allgather(&local_size, 1, MPI_INT, union_sizes.data(), 1, MPI_INT, MPI_COMM_WORLD);

    // Calculate displacements for global union list
    std::vector<int> displs(size, 0);
    int              total_union_size = 0;
    for(int i = 0; i < size; ++i) {
      displs[i] = total_union_size;
      total_union_size += union_sizes[i];
    }

    // Gather all unions to construct the global union list
    std::vector<std::pair<int, int>> global_union_list(total_union_size);
    MPI_Allgatherv(local_union_list.data(), local_size, MPI_2INT, global_union_list.data(), union_sizes.data(),
                   displs.data(), MPI_2INT, MPI_COMM_WORLD);

    auto new_labels = calculate_disjoint_sets(global_union_list);

    // Assign the compacted labels back to the grid_function grid function
    for(int i = 0; i < grid_function.Size(); ++i) {
      int const old_label = grid_function[i];

      if(old_label >= 0) {
        grid_function[i] = new_labels[old_label];
      }
    }
  }

}  // namespace

void assign_connected_components(mfem::ParGridFunction &grid_function)
{
  // Simple function where we only consider values > 0 to be foreground
  auto is_label = [](double value) -> bool { return value > 0.0; };

  mfem::ParMesh &par_mesh = *grid_function.ParFESpace()->GetParMesh();

  // Phase 1: Identify the components within a processor
  int const local_number_components = phase_1_local_components(par_mesh, grid_function, is_label);

  // Phase 2: Re-number the components globally
  // We collect all labels and ensure that each label is globally unique
  int const rank_offset =
    phase_2_relabel_for_cross_processor_comparison(par_mesh, grid_function, local_number_components);

  // Phase 3: Merging of labels across processors add identity
  auto const local_unions =
    phase_3_prime_merge_labels_across_processors(par_mesh, grid_function, local_number_components, rank_offset);

  // Phase 4: Final assignment of labels across processors
  phase_4_final_assignment_of_labels(grid_function, local_unions);
}

DisjointSet.hpp

#pragma once

#include <algorithm>
#include <cstdint>
#include <unordered_map>
#include <vector>

//! @brief A class implementing the Union-Find (Disjoint Set Union) data structure.
//!
//!
class UnionFind
{
 public:
  //! @brief Constructs a Union-Find structure with a given size.
  //!
  //! @param size The number of elements in the Union-Find structure.
  //!
  explicit UnionFind(int size) : parent(size), rank(size, 0)
  {
    for(int i = 0; i < size; ++i) {
      parent[i] = i;
    }
  }

  //! @brief Finds the root of an element with path compression.
  //!
  //! @param x The element whose root is to be found.
  //!
  //! @return The root of the element.
  //!
  int find(int x)
  {
    if(parent[x] != x) {
      parent[x] = find(parent[x]);  // Path compression
    }
    return parent[x];
  }

  //! @brief Unites two elements into the same set.
  //!
  //! @param x The first element.
  //! @param y The second element.
  //!
  void unite(int x, int y)
  {
    int root_x = find(x);
    int root_y = find(y);
    if(root_x != root_y) {
      if(rank[root_x] > rank[root_y]) {
        parent[root_y] = root_x;
      } else if(rank[root_x] < rank[root_y]) {
        parent[root_x] = root_y;
      } else {
        parent[root_y] = root_x;
        rank[root_x]++;
      }
    }
  }

 private:
  // The parent vector for storing the set representatives.
  std::vector<int> parent;

  // The rank vector for optimizing union operations.
  std::vector<int> rank;
};

//! @brief Processes a global union list and assigns new labels to connected components.
//!
//! @param global_union_list A list of pairs representing union operations.
//! @return A mapping from original indices to new labels.
//!
inline std::unordered_map<int, int> calculate_disjoint_sets(std::vector<std::pair<int, int>> &global_union_list)
{
  // Find the maximum element to initialize Union-Find
  int max_element = 0;
  for(auto const &pair: global_union_list) {
    max_element = std::max({max_element, pair.first, pair.second});
  }

  UnionFind uf(max_element + 1);

  // Perform union operations
  for(auto const &pair: global_union_list) {
    uf.unite(pair.first, pair.second);
  }

  // Assign labels
  std::unordered_map<int, int> new_labels;
  int                          current_label = 0;
  std::unordered_map<int, int> label_map;

  for(int i = 0; i <= max_element; ++i) {
    int root = uf.find(i);
    if(label_map.find(root) == label_map.end()) {
      label_map[root] = current_label++;
    }
    new_labels[i] = label_map[root];
  }

  return new_labels;
}

Hope that helps.

Paul

[paul-hilscher]

@v-dobrev ,

would there be some interest on your end if I bring in the algorithm from https://www.osti.gov/servlets/purl/1048292 as a miniapp ? Those ghost cell exchange questions pop up occasionally so imho a specific example inside mfem would be helpful to the community.

[starkekr]

Thank you! It will take me some time to digest this but I really appreciate it. Yes I think so!

[v-dobrev]

Hi @paul-hilscher,

Yes, I think adding such a miniapp will be useful. Let me know if there's anything I can do to help.

[starkekr]

I get an error:

Verification failed: (el_to_face != NULL) is false:
--> el_to_face not generated
... in function: void mfem::Mesh::GetElementFaces(int, mfem::Array&, mfem::Array&) const
... in file: mesh/mesh.cpp:7203

in the call _pmesh->GetElementFaces(j, faces, orientations);
Do i need to do anything to generate this Table?

I am setting up the mesh like this:

// 1. Initialize MPI and HYPRE.
Mpi::Init(argc, argv);
Hypre::Init();

// 2. Parse command line options.
int order = 1;

OptionsParser args(argc, argv);
args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
args.ParseCheck();

Mesh *cart_mesh;
int nx = 5;
cart_mesh = new Mesh(Mesh::MakeCartesian2D(nx, nx, Element::QUADRILATERAL,false,1,1,true));
cart_mesh->EnsureNCMesh(true);
cart_mesh->UniformRefinement(1);
// 4. Define a parallel mesh by a partitioning of the serial mesh. Refine
//    this mesh once in parallel to increase the resolution.
ParMesh mesh(MPI_COMM_WORLD, *cart_mesh);
mesh.EnsureNCMesh(true);
mesh.UniformRefinement(1);
mesh.Finalize();
Array<int> faces, orientations;
mesh.GetElementFaces(0, faces, orientations);

[v-dobrev]

Try using true instead of false in the 4th argument of Mesh::MakeCartesian2D -- I think that will fix the error.