elementcoloring.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_FUFEM_PARALLEL_ELEMENTCOLORING_HH
#define DUNE_FUFEM_PARALLEL_ELEMENTCOLORING_HH
 
#include <cstddef>
#include <type_traits>
#include <utility>
#if __has_include(<execution>)
#include <execution>
#else
#include <version>
#endif
#include <vector>
#include <unordered_map>
#include <queue>
#include <algorithm>
#include <limits>
#include <iostream>
#include <initializer_list>
 
#include <dune/common/timer.hh>
#include <dune/common/rangeutilities.hh>
 
#include <dune/geometry/referenceelements.hh>
 
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/common/scsgmapper.hh>
 
 
 
namespace Dune::Fufem::Impl {
 
 
 
/* A class for building a container that looks like std::vector<std::vector<T>>
 * but uses flat storage. The class should be used as follows:
 *
 * // Create with given outer size
 * auto vec = FlatBlockedVector<T>(size);
 *
 * // Reserve size for all entries.
 * // It's allowed to do this in arbitrary order.
 * vec.capacity(2) = ...;
 * vec.capacity(1) = ...;
 * vec.capacity(2) += ...;
 *
 * // Initialize blocks
 * vec.initialize();
 *
 * // Fill blocks in arbitrary order
 * vec[i].push_back(value);
 * vec[i].push_back(value);
 *
 * // Swap or access entries using operator[] or iterators
 * // or resize within reserved capacity
 * vec[i][0] = ...;
 * auto&& vec_i = vec[i];
 * std::sort(vec_i.begin(), vec_i.end());
 * vec[i].resize(1);
 *
 * // Export a range of ranges. This will steal the internal data.
 * auto r = std::move(vec).toRange();
 *
 * // Use r like a nested range
 * for (auto&& ri : r)
 *   for (auto&& rij : ri)
 *     foo(rij);
 */
template<class T, class st = std::size_t>
class FlatBlockedVector
{
public:
  using size_type = st;
 
private:
 
  struct MutableBlock
  {
    void push_back(const T& t)
    {
      data_[end_] = t;
      ++end_;
    }
 
    T& operator[](size_type i)
    {
      return data_[begin_+i];
    }
 
    auto begin()
    {
      return data_.begin() + begin_;
    }
 
    auto end()
    {
      return data_.begin() + end_;
    }
 
    void resize(size_type size)
    {
      end_ = begin_+size;
    }
 
    std::vector<T>& data_;
    size_type& begin_;
    size_type& end_;
  };
 
public:
 
  FlatBlockedVector(size_type size):
    blockBegin_(size, 0),
    blockEnd_(size, 0)
  {
    blockEnd_.clear();
    blockEnd_.resize(size, 0);
  }
 
  size_type size() const
  {
    return blockBegin_.size();
  }
 
  size_type& capacity(size_type i)
  {
    return blockEnd_[i];
  }
 
  void initialize()
  {
    auto size = blockBegin_.size();
    blockBegin_[0] = 0;
    for(auto i : Dune::range(std::size_t(1), size))
    {
      blockBegin_[i] = blockBegin_[i-1] + capacity(i-1);
      blockEnd_[i-1] = blockBegin_[i-1];
    }
    data_.resize(blockBegin_[size-1] + capacity(size-1));
    blockEnd_[size-1] = blockBegin_[size-1];
  }
 
  auto operator[] (size_type i) const
  {
    return Dune::transformedRangeView(Dune::range(blockBegin_[i], blockEnd_[i]), [&](auto j) {
      return data_[j];
    });
  }
 
  auto operator[] (size_type i)
  {
    return MutableBlock{data_, blockBegin_[i], blockEnd_[i]};
  }
 
  auto toRange() &&
  {
    auto size = blockBegin_.size();
    return Dune::transformedRangeView(Dune::range(size),
      [setVector=std::move(*this)](auto i) {
        return setVector[i];
      });
  }
 
  std::vector<T> data_;
  std::vector<size_type> blockBegin_;
  std::vector<size_type> blockEnd_;
};
 
 
 
template<class V>
void sortAndUnifyVector(V&& v)
{
  std::sort(v.begin(), v.end());
  auto lastUnique = std::unique(v.begin(), v.end());
  v.resize(lastUnique-v.begin());
}
 
 
 
}
 
 
 
namespace Dune {
namespace Fufem {
 
 
 
template <class Adjacency, class SeedNodesContainer>
auto computeColoring(const Adjacency& adjacency, SeedNodesContainer&& seedNodes)
{
 
  // Hack to access i-th entry without adjacency[i].
  // This is necessary to support TransformedRangeView,
  // because the latter does not provide operator[]
  // in dune-2.9 although the range is random-access.
  auto adjacentNodes = [&](auto i) {
    return adjacency.begin()[i];
  };
 
  using size_type = std::decay_t<decltype(*adjacentNodes(0).begin())>;
 
  // Use magic values to tag fully unprocessed and queued (but non-colored) entries.
  size_type nonProcessedStatus = std::numeric_limits<size_type>::max();
  size_type queuedStatus = std::numeric_limits<size_type>::max()-1;
 
  // Initialize result vector with unprocessed state for all entries.
  std::vector<size_type> color(adjacency.size(), nonProcessedStatus);
 
  // Queue for the advancing front method
  auto queue = std::queue<size_type>();
 
  // Some helper methods to make code more expressive
  auto isColored = [&](const auto& i) {
    return color[i]<queuedStatus;
  };
  auto isNonProcessed = [&](const auto& i) {
    return color[i]==nonProcessedStatus;
  };
  auto queueNode = [&](const auto& i) {
    color[i] = queuedStatus;
    queue.push(i);
  };
 
  // Fill queue with seed nodes (or {0} if seed node container is empty)
  for(size_type i : seedNodes)
    queueNode(i);
  if (queue.empty())
    queueNode(0);
 
  // A temporary data structure to store the colors of already
  // colored neighbors.
  //
  // There's various possible data structures here. The most
  // obvious would be a set-like container storing the used
  // color indices. This could be (with increasing performance)
  // std::set, std::unordered_set, std::vector (keep sorted and unique),
  // std::vector (just append, sort and remove duplicates afterwards).
  // Since we expect that the total number of colors is small,
  // it is significantly faster to simply mark used colors in
  // a flag vector. Because memory usage is not an issue here,
  // we use std::vector<char> which measurably outperformes
  // std::vector<bool> due to faster bit access.
  std::vector<char> colorIsUsed;
  colorIsUsed.resize(32);
 
  // If the graph has several connected components that are not
  // reachable by the seed nodes, then we need to do several sweeps
  // of the advancing front method.
  while (not queue.empty())
  {
    // Use advancing front method to traverse graph
    while (not queue.empty())
    {
      // Process next entry from queue
      auto i = queue.front();
      queue.pop();
 
      // Collect colors of already colored neighbors,
      // put nonprocessed neighbors into queue,
      // and do nothing with already queued neighbors.
      colorIsUsed.assign(colorIsUsed.size(), false);
      for (const auto& neighbor: adjacentNodes(i))
      {
        if (isColored(neighbor))
          colorIsUsed[color[neighbor]] = true;
        else if (isNonProcessed(neighbor))
          queueNode(neighbor);
      }
 
      // Assign first color that is not used by any processed neighbor to the current node
      color[i] = std::find(colorIsUsed.begin(), colorIsUsed.end(), false) - colorIsUsed.begin();
 
      // If necessary: Resize flag vector
      if (color[i]==colorIsUsed.size())
        colorIsUsed.resize(color[i]+1);
    }
 
    // If there are still some uncolored nodes, then the graph has connected
    // components that are not reachable by the provided seed nodes. Hence we
    // do another advancing front sweep starting from the first uncolored node.
    auto firstNonProcessedIt = std::find(color.begin(), color.end(), nonProcessedStatus);
    if (firstNonProcessedIt!=color.end())
      queueNode(firstNonProcessedIt - color.begin());
  }
 
  return color;
}
 
 
 
template <class Adjacency, class T>
auto computeColoring(const Adjacency& adjacency, std::initializer_list<T> seedIndices)
{
  return computeColoring(adjacency, std::vector<T>{seedIndices});
}
 
 
 
template <class Adjacency>
auto computeColoring(const Adjacency& adjacency)
{
  return computeColoring(adjacency, {0});
}
 
 
 
template <class CliqueList>
auto cliqueListToAdjacency(const CliqueList& cliqueList, std::size_t size, bool verbose=false)
{
  using size_type = std::decay_t<decltype(*(cliqueList.begin()->begin()))>;
 
  // In verbode mode: Compute size of the largest provided clique.
  // This is a lower bound for the clique number (largest existing clique)
  // which itself is a lower bound for the chromatic number.
  // Notice that even for planar 2d grids the element graph is not planar
  // and thus we may need more than 4 colors.
  if (verbose)
  {
    size_type maxCliqueSize = 0;
    for(const auto& clique : cliqueList)
      maxCliqueSize = std::max(size_type(maxCliqueSize), size_type(clique.size()));
    std::cout << "Largest found element clique (lower bound of required colors) has size " << maxCliqueSize << std::endl;
  }
 
  // Theoretically we could use a std::set for each adjacency[k]
  // and simply insert all entries from the cliques into this set.
  // However, simply pushing back entries and removing duplicates
  // afterwards is faster.
  // Furthermore using windows into a flat std::vector<size_type>
  // is faster than using a std::vector<std::vector<size_type>>.
  // The latter is provided by the FlatBlockedVector helper class.
  auto adjacency = Impl::FlatBlockedVector<size_type, size_type>(size);
 
  // We first must compute an upper bound for the
  // degree of each vertex to reserve the appropriate
  // block in the flat vector.
  for (const auto& clique : cliqueList)
    for (auto k : clique)
      adjacency.capacity(k) += clique.size()-1;
  adjacency.initialize();
 
  // Add edges to adjacency list
  for (const auto& clique : cliqueList)
    for (auto k : clique)
      for (auto l : clique)
        if (k!=l)
          adjacency[k].push_back(l);
 
  // Remove duplicate entries in inner container
  {
    auto indices = Dune::range(adjacency.size());
#if __cpp_lib_execution >= 201603L
    std::for_each(std::execution::par, indices.begin(), indices.end(), [&](auto i) {
#else
    std::for_each(indices.begin(), indices.end(), [&](auto i) {
#endif
      Impl::sortAndUnifyVector(adjacency[i]);
    });
  }
 
  return std::move(adjacency).toRange();
}
 
 
 
template <class GridView, class ElementMapper, class size_type=uint32_t>
auto vertexBasedElementAdjacency(const GridView& gridView, const ElementMapper& elementMapper, bool verbose=false)
{
  static const auto dim = GridView::dimension;
 
  auto subEntityMapper = Dune::SingleCodimSingleGeomTypeMapper<GridView, dim>(gridView);
 
  // Compute a map of vertex index to vector of adjacent elements indices.
  auto vertexCliques = Impl::FlatBlockedVector<size_type>(subEntityMapper.size());
 
  // Precompute and reserve block sizes
  for (const auto& element : elements(gridView))
    for (std::size_t i = 0; i<element.subEntities(dim); ++i)
      ++vertexCliques.capacity(subEntityMapper.subIndex(element, i, dim));
  vertexCliques.initialize();
 
  for (const auto& element : elements(gridView))
  {
    auto elementIndex = elementMapper.index(element);
    for (std::size_t i = 0; i<element.subEntities(dim); ++i)
      vertexCliques[subEntityMapper.subIndex(element, i, dim)].push_back(elementIndex);
  }
 
  return cliqueListToAdjacency(std::move(vertexCliques).toRange(), elementMapper.size(), verbose);
}
 
 
 
template <class GridView, class size_type=uint32_t>
auto vertexBasedElementAdjacency(const GridView& gridView, bool verbose=false)
{
  auto elementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>(gridView, mcmgElementLayout());
  return vertexBasedElementAdjacency(gridView, elementMapper, verbose);
}
 
 
 
template <class GridView, class SubEntityLayout, class ElementMapper, class size_type=uint32_t>
auto subEntityBasedElementAdjacency(const GridView& gridView, const SubEntityLayout& subEntityLayout, const ElementMapper& elementMapper)
{
  static const auto dim = GridView::dimension;
 
  auto subEntityMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>(gridView, subEntityLayout);
 
  std::vector<std::vector<size_type>> subEntityCliques(subEntityMapper.size());
  for (const auto& element : elements(gridView))
  {
    auto elementIndex = elementMapper.index(element);
    auto re = Dune::referenceElement<double,dim>(element.type());
    for (std::size_t codim = 1; codim <= dim; ++codim)
      for (std::size_t i = 0; i<re.size(codim); ++i)
        if (subEntityLayout(re.type(i, codim), GridView::dimension))
          subEntityCliques[subEntityMapper.subIndex(element, i, codim)].push_back(elementIndex);
  }
 
  return cliqueListToAdjacency(subEntityCliques, elementMapper.size());
}
 
 
 
template <class GridView, class SubEntityLayout, class size_type=uint32_t>
auto subEntityBasedElementAdjacency(const GridView& gridView, const SubEntityLayout& subEntityLayout)
{
  auto elementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>(gridView, mcmgElementLayout());
  return subEntityBasedElementAdjacency(gridView, elementMapper, subEntityLayout);
}
 
 
 
template<class GridView, class size_type>
auto coloredGridViewPartition(const GridView& gridView, const std::vector<size_type>& coloring)
{
  size_type colorCount = *std::max_element(coloring.begin(), coloring.end()) + 1;
 
  using ElementSeed = typename GridView::Grid::template Codim<0>::EntitySeed;
  using ElementSeedRange = std::vector<ElementSeed>;
 
  std::vector<ElementSeedRange> elementSeedPartition(colorCount);
  auto elementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>(gridView, mcmgElementLayout());
  for (const auto& element : elements(gridView))
    elementSeedPartition[coloring[elementMapper.index(element)]].push_back(element.seed());
 
  // For convenience transform from range of seeds to range of entities
  // This is save because transformedRangeView() stores copies of passed
  // r-values. Furthermore it transforms random access ranges to random
  // access ranges, such that the resulting entity range can be split
  // cheaply.
  auto seedToElement = [gridPtr = &gridView.grid()] (const auto& seed) {
    return gridPtr->entity(seed);
  };
  return Dune::transformedRangeView(std::move(elementSeedPartition), [seedToElement](const auto& seedRange) {
    return Dune::transformedRangeView(seedRange, seedToElement);
  });
}
 
 
 
template<class GridView, class size_type=uint32_t>
auto coloredGridViewPartition(const GridView& gridView, bool verbose=false)
{
  Dune::Timer timer;
 
  timer.reset();
  auto adjacency = Dune::Fufem::vertexBasedElementAdjacency(gridView, verbose);
  if (verbose)
    std::cout << "Building adjacency map took " << timer.elapsed() << "s." << std::endl;
 
  timer.reset();
  auto coloring = Dune::Fufem::computeColoring(adjacency);
  if (verbose)
    std::cout << "Building coloring took " << timer.elapsed() << "s." << std::endl;
 
  timer.reset();
  auto elementPartition = Dune::Fufem::coloredGridViewPartition(gridView, coloring);
  if (verbose)
  {
    std::cout << "Building colored element partition took " << timer.elapsed() << "s." << std::endl;
    std::cout << "Number of used colors is " << elementPartition.size() << std::endl;
    std::size_t elementCount = 0;
    for(const auto& color : elementPartition)
      elementCount += color.size();
    std::cout << "Total number of elements in the partition is " << elementCount << std::endl;
  }
 
  return elementPartition;
}
 
 
 
template<class ElementRange, class Basis>
bool checkOverlapFree(const ElementRange& elementRange, const Basis& basis) {
 
  bool overlapFree = true;
 
  auto elementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<typename Basis::GridView>(basis.gridView(), mcmgElementLayout());
 
  using MI = typename Basis::MultiIndex;
  std::unordered_map<MI, std::size_t> processedDOFs;
 
  auto localView = basis.localView();
  for (const auto& element : elementRange)
  {
    localView.bind(element);
 
    // Check if any of the DOFs of the present element showed up earlier.
    for (size_t i=0; i<localView.tree().size(); ++i)
    {
      auto localIndex = localView.tree().localIndex(i);
      auto globalIndex = localView.index(localIndex);
      auto it = processedDOFs.find(globalIndex);
      if (it != processedDOFs.end())
      {
        std::cout << "Overlap detected! Index [" << globalIndex << "] shows up in element "
          << it->second << " and " << elementMapper.index(element) << std::endl;
        overlapFree =false;
      }
    }
 
    for (size_t i=0; i<localView.tree().size(); ++i)
    {
      auto localIndex = localView.tree().localIndex(i);
      auto globalIndex = localView.index(localIndex);
      processedDOFs[globalIndex] = elementMapper.index(element);
    }
  }
  return overlapFree;
}
 
 
 
} // namespace Fufem
} // namespace Dune
 
 
#endif // DUNE_FUFEM_PARALLEL_ELEMENTCOLORING_HH