Problem Statement

This example comes from a mesh preprocessing step in a finite element library. For some background: consider the simple triangle mesh

There are many different kinds of connectivity

• element to face

• element to edge

• element to vertex

• face to edge

• face to vertex

• edge to vertex

as well as the corresponding inverse connectivities (e.g. vertex to face). All of this connectivity information can be derived from the element-to-vertex connectivity table, so most mesh file formats choose to only store the element-to-vertex information to keep the file size small.

In this document, we'll walk through a few algorithms (on CPU and GPU) to derive the other kinds of connectivity information from the element-to-vertex data:

In practice we're mostly interested in 3D meshes, but to keep this document's explanations and diagrams simple, we'll be considering a 2D triangle mesh. In this example, we'll assume we have a mesh given in terms of triangle-to-vertex connectivity and we would like to determine the edge-to-vertex and triangle-to-edge information.

The overall algorithmic ideas generalize immediately to 3D meshes, and if the reader wishes to see the actual implementation for tetrahedron meshes, please see the attached source code. The performance of these different approaches will be measured on a mesh with 1572864 tetrahedra and 274625 vertices.

CPU Reference Implementation

To me, the simplest approach is to first derive the edge-to-vertex information and then register which of those edges belong to each element. We can do this by iterating over each triangle in the mesh and using std::unordered_map to insert and deduplicate edges, by using the edge's pair of vertex ids as the key type:

​x
// std::unordered_map requires a hash function for the key types
struct array_hasher {
  template<size_t n>
  std::size_t operator()(const std::array< uint32_t, n > & arr) const {
    uint32_t seed = 0;
    for(const auto elem : arr) {
      seed ^= std::hash<uint32_t>()(elem) + 0x9e3779b9 + (seed<<6) + (seed>>2);
    }
    return seed;
  }
};

using edge_map = std::unordered_map< 
  std::array< uint32_t, 2 >, 
  uint32_t, 
  array_hasher 
>;

The uint32_t value for the edge_map holds the edge's index (which we increment every time a new edge is inserted).

As we iterate over the elements, we try to insert each of that element's edges into the edge_map (note: sort the vertex ids for an edge before inserting it, so that we can detected duplicates with different orientations). Each time a new edge is inserted into the edge_map, we increment a counter tracking the total number of unique edges. The element then records that edge's index in its own element-to-edge connectivity table.

xxxxxxxxxx
// our edge ordering convention for the triangle
constexpr uint32_t local_edge_ids[3][2] = {{0, 1}, {1, 2}, {2, 0}};

edge_map sorted_edges;          // sorted vertex ids -> edge index
std::vector< edge > edges;      // edge     -> vertices
std::vector< tri_edges > conn;  // triangle -> edges

for (i : triangles) {
  for (j : {0, 1, 2}) {

    // the vertices of this triangle's local edge j
    edge e = {tris[i][local_edge_ids[j][0]], tris[i][local_edge_ids[j][1]]};

    // sorting the keys before insertion accounts edges with different orientations
    auto [it, is_new] = sorted_edges.insert({sorted(e), edges.size()});

    // the first time an edge is visited, it gets added to the edge map
    if (is_new) edges.push_back(e);

    // record the triangle->edge connectivity
    conn[i][j] = it->second;
    
  }
}

Running that loop on the mesh from before, one iteration at a time:

The runtime of this reference implementation on a mesh is 2400 ms.

CPU Parallel Implementation

The previous implementation is slow, but part of that is expected since it only runs on a single thread. However, it's not entirely obvious how to leverage multiple cores for the algorithm as written: std::unordered_map is not thread safe under insertion.

Luckily, this problem is relatively straightforward to solve: there are several thread safe map-like containers available in the open source community (abseil, gtl, tbb, etc). We'll use one from the library gtl in this example, gtl::parallel_flat_hash_map.

There are still more details involved in updating the algorithm to work in parallel (e.g. how to assign edge ids in a way that is both deterministic and in agreement with the serial CPU implementation), but this is more of a CUDA-focused document, so we'll gloss over those details in this document. If the reader is so inclined, please see the attached code for the complete implementation details on how to adapt the serial algorithm to use the gtl container.

The runtime of our parallel variant of the CPU code is: 314 ms (~7.6x speedup with 32 threads)

GPU Implementation

Getting this code running on the GPU isn't as simple as switching out a container. Although there are map-like containers in the cuCollections library, they likely wouldn't perform well in this situation since every single thread is trying to insert and access the same container. Instead, let's consider a different algorithm altogether.

Instead of identifying and deduplicating edges incrementally, we can also do it in parallel by breaking it up into separate steps. The first step in this process is to emit metadata from each triangle/edge pair

By construction, this step is embarrassingly parallel: each thread is responsible for processing a particular triangle/edge pair and writing out the corresponding data to the appropriate location in the output array.

The next step is to sort these rows by the keys (and break ties based on the element index). This has the effect of grouping duplicate keys, so we can easily assign a deterministic edge ordering and determine whether an edge is on the interior or boundary of the mesh (by checking if that key appears once or twice).

From here, we can choose to keep the first key (i.e. the element with the lowest index "wins"), and then prefix sum those row markings to determine the number of unique edges and determine the edge ids.

On its face, this approach may sound wasteful since it involves 4 steps (emit, sort, mark, prefix sum) instead of 1, but the runtime tells a different story: 5.5ms (~57x faster than the 32-thread version)!

Why is it so much faster? There may be some work inefficiency, but each of these operations is embarassingly parallelizable: no mutexes, atomics, or data races. As a result, each operation can take full advantage of the GPU hardware.

Summary

A common mistake I see developers make when migrating code to the GPU is trying to directly map every CPU part of the algorithm over to the GPU, one-to-one. In practice, this rarely works well. Many CPU programming idioms are unnatural or impractical to implement on a GPU, and trying to do so is like trying to jam a square peg in a round hole. It's important to distinguish between what an algorithm does versus how an algorithm does it. A good GPU algorithm may look fundamentally different than the CPU version, and that's okay!

The sourcefiles for this code are available here