What is this code doing?

```
template<class LO, class GO, class Node>
void generate_fem2d_q1_graph(size_t numCells1D, RCP<const Comm<int> > comm , GraphPack<4,LO,GO,Node> & pack) {

  // We assume a NxN Q1 fem grid. Cells are partitioned among ranks.
  // Cell/Nodes ids are as follows (assuming N=2):
  //
  //  0---1---2
  //  | 0 | 1 |
  //  3---4---5
  //  | 2 | 3 |
  //  6---7---8
  //
  // The order in which we store the ids within a cell is irrelevant (so long as consistent),
  // since the local pattern is full anyways.
  //
  // NOTE: this routine sets overlapMap NOT to be the "col" map. Instead, since we partition
  //       cells, it contains all nodes belonging to local cells. If all cells are "fully" owned,
  //       then locally overlapMap=uniqueMap, while the "col" map of the graph will ALWAYS
  //       have something more (the halo).

  const size_t rank    = comm->getRank();
  const size_t numProc = comm->getSize();

  const size_t numGlobalCells = numCells1D*numCells1D;
  const size_t numNodes1D     = numCells1D+1;

  size_t numMyCells  = numGlobalCells / numProc;
  size_t remainder   = numGlobalCells % numProc;
  size_t myCellStart = numMyCells*rank;
  if (rank < remainder) {
    ++numMyCells;
    myCellStart += rank;
  } else {
    myCellStart += remainder;
  }

  // Add also repeated GIDs, we'll remove duplicates later
  Teuchos::Array<GO> ovNodes;
  pack.element2node.resize(numMyCells);
  for (size_t cell=0; cell<numMyCells; ++cell) {
    //  0---1---2
    //  | 0 | 1 |
    //  3---4---5
    //  | 2 | 3 |
    //  6---7---8

    auto cellId = myCellStart+cell;

    auto icell = cellId / numCells1D;
    auto jcell = cellId % numCells1D;

    auto offset = icell*numNodes1D;

    // Store for the map
    ovNodes.append(offset+jcell);
    ovNodes.append(offset+jcell+1);
    ovNodes.append(offset+jcell+numNodes1D);
    ovNodes.append(offset+jcell+numNodes1D+1);

    // Store for the assembly
    pack.element2node[cell].resize(4);
    pack.element2node[cell][0] = offset + jcell;
    pack.element2node[cell][1] = offset + jcell + 1;
    pack.element2node[cell][2] = offset + jcell + numNodes1D;
    pack.element2node[cell][3] = offset + jcell + numNodes1D + 1;
  }

  // Remove duplicates. Note: std::unique needs consecutive duplicates, so sort first
  auto start = ovNodes.data();
  auto end   = ovNodes.data()+ovNodes.size();
  std::sort(start,end);
  auto new_end = std::unique(start,end);
  auto new_size = std::distance(start, new_end);
  ovNodes.resize(new_size);

  // Create overlap map, and use Tpetra utility to create the unique one
  pack.overlapMap = createNonContigMapWithNode<LO,GO,Node>(ovNodes(),comm);
  pack.uniqueMap  = createOneToOne(pack.overlapMap);

  // Kokkos version of the element2node array
  Kokkos::resize(pack.k_element2node,numMyCells);
  auto k_e2n = pack.k_element2node;
  Kokkos::parallel_for(Kokkos::RangePolicy<typename Node::execution_space>(0,numMyCells),
                       KOKKOS_LAMBDA(const size_t cell) {
      auto cellId = myCellStart+cell;

      auto icell = cellId / numCells1D;
      auto jcell = cellId % numCells1D;

      auto offset = icell*numNodes1D;

      k_e2n(icell,0) = offset + jcell;
      k_e2n(icell,1) = offset + jcell + 1;
      k_e2n(icell,2) = offset + jcell + numNodes1D;
      k_e2n(icell,3) = offset + jcell + numNodes1D + 1;
    });
  Kokkos::fence();
}
```