gridcoloring.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright © DUNE 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_ASSEMBLER_PARALLEL_GRIDCOLORING_HH
#define DUNE_ASSEMBLER_PARALLEL_GRIDCOLORING_HH
 
#include <algorithm>
#include <cstddef>
#include <iostream>
#include <unordered_map>
#include <utility>
#include <vector>
 
#include <dune/common/rangeutilities.hh>
#include <dune/common/timer.hh>
 
#include <dune/geometry/referenceelements.hh>
 
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/common/scsgmapper.hh>
 
#include <dune/assembler/parallel/graphcoloring.hh>
#include <dune/assembler/parallel/graphutilities.hh>
 
namespace Dune::Assembler::Experimental {
 
 
 
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 coloredElementRange (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> coloredElementSeedRange(colorCount);
  auto elementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>(gridView, mcmgElementLayout());
  for (const auto& element : elements(gridView))
    coloredElementSeedRange[coloring[elementMapper.index(element)]].push_back(element.seed());
 
  // For convenience transform from range of seeds to range of entities
  // This is safe 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(coloredElementSeedRange), [seedToElement](const auto& seedRange) {
    return Dune::transformedRangeView(seedRange, seedToElement);
  });
}
 
 
 
template <class GridView, class size_type=uint32_t>
auto coloredElementRange (const GridView& gridView, bool verbose=false)
{
  Dune::Timer timer;
 
  timer.reset();
  auto adjacency = Dune::Assembler::Experimental::vertexBasedElementAdjacency(gridView, verbose);
  if (verbose)
    std::cout << "Building adjacency map took " << timer.elapsed() << "s." << std::endl;
 
  timer.reset();
  auto coloring = Dune::Assembler::Experimental::computeColoring(adjacency);
  if (verbose)
    std::cout << "Building coloring took " << timer.elapsed() << "s." << std::endl;
 
  timer.reset();
  auto coloredElementRange = Dune::Assembler::Experimental::coloredElementRange(gridView, coloring);
  if (verbose)
  {
    std::cout << "Building colored element range took " << timer.elapsed() << "s." << std::endl;
    std::cout << "Number of used colors is " << coloredElementRange.size() << std::endl;
    std::size_t elementCount = 0;
    for (const auto& color : coloredElementRange)
      elementCount += color.size();
    std::cout << "Total number of elements in the range is " << elementCount << std::endl;
  }
 
  return coloredElementRange;
}
 
 
 
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;
}
 
 
} // end namespace Dune::Assembler::Experimental
 
#endif // DUNE_ASSEMBLER_PARALLEL_GRIDCOLORING_HH