repartition.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <deal.II/base/exceptions.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/numbers.h>
#include <deal.II/base/utilities.h>
 
#include <deal.II/distributed/tria.h>
 
#include <deal.II/dofs/dof_handler.h>
 
#include <deal.II/grid/tria.h>
 
#include <deal.II/non_matching/mesh_classifier.h>
 
#include <deal.II/numerics/solution_transfer.h>
 
#include <meltpooldg/cut/cut_data.hpp>
#include <meltpooldg/utilities/attach_vectors.hpp>
 
#include <boost/range/iterator_range.hpp>
#include <boost/signals2/connection.hpp>
 
#include <algorithm>
#include <functional>
#include <memory>
#include <vector>
 
namespace MeltPoolDG::CutUtil
{
  template <int dim>
  class PartitionWeights
  {
  public:
    using CellIterator = typename dealii::Triangulation<dim>::cell_iterator;
 
    explicit PartitionWeights(const dealii::Triangulation<dim> &tria_in,
                              const CutRepartitionData         &repartition_data_in)
      : tria(&tria_in)
      , repartition_data(repartition_data_in)
    {
      AssertThrow(
        (dynamic_cast<const dealii::parallel::distributed::Triangulation<dim> *>(&tria_in) !=
         nullptr),
        dealii::ExcMessage(
          "The weighted repartitioning requires a parallel::distributed::Triangulation."));
    }
 
    void
    update(const dealii::NonMatching::MeshClassifier<dim> &mesh_classifier)
    {
      // resize weight vector and initialize to 1
      weights.assign(tria->n_active_cells(), 1);
 
      // loop over cells and set the partitioning weight
      for (const auto &cell : tria->active_cell_iterators())
        if (cell->is_locally_owned())
          weights[cell->active_cell_index()] =
            get_partition_weight(mesh_classifier.location_to_level_set(cell));
 
      if (connections.empty())
        connect_to_triangulation();
    }
 
    unsigned int
    get_weight(const CellIterator &cell) const
    {
      if (cell->has_children())
        {
          unsigned int weight = 0;
          for (unsigned int c = 0; c < cell->n_children(); ++c)
            weight = std::max(weight, get_weight(cell->child(c)));
          return weight;
        }
 
      const unsigned int index = cell->active_cell_index();
 
      // assert that the cached weights are up-to-date with the current triangulation
      Assert(index < weights.size(),
             dealii::ExcMessage("The cached partitioning weights are outdated. Call "
                                "PartitionWeights::update() after each change of the "
                                "triangulation."));
 
      return weights[index];
    }
 
    dealii::Utilities::MPI::MinMaxAvg
    compute_load_statistics(const MPI_Comm mpi_comm) const
    {
      double local_load = 0.;
      for (const auto &cell : tria->active_cell_iterators())
        if (cell->is_locally_owned())
          local_load += get_weight(cell);
 
      return dealii::Utilities::MPI::min_max_avg(local_load, mpi_comm);
    }
 
  private:
    const dealii::Triangulation<dim> *tria;
 
    const CutRepartitionData repartition_data;
 
    std::vector<unsigned int> weights;
 
    unsigned int weights_handle = dealii::numbers::invalid_unsigned_int;
 
    std::vector<boost::signals2::scoped_connection> connections;
 
    void
    connect_to_triangulation()
    {
      auto &signals = tria->signals;
 
      connections.emplace_back(signals.weight.connect(
        [this](const CellIterator &cell, const dealii::CellStatus) -> unsigned int {
          return this->get_weight(cell);
        }));
 
      for (auto *signal :
           {&signals.pre_distributed_refinement, &signals.pre_distributed_repartition})
        connections.emplace_back(signal->connect([this]() { this->prepare_transfer(); }));
 
      for (auto *signal :
           {&signals.post_distributed_refinement, &signals.post_distributed_repartition})
        connections.emplace_back(signal->connect([this]() { this->finish_transfer(); }));
    }
 
    void
    prepare_transfer()
    {
      Assert(weights_handle == dealii::numbers::invalid_unsigned_int, dealii::ExcInternalError());
 
      // get_weight() returns the maximum weight of the children for cells whose children will be
      // coarsened
      weights_handle = const_cast<dealii::Triangulation<dim> &>(*tria).register_data_attach(
        [this](const CellIterator &cell, const dealii::CellStatus) {
          return dealii::Utilities::pack(this->get_weight(cell), /*allow_compression=*/false);
        },
        /*returns_variable_size_data=*/false);
    }
 
    void
    finish_transfer()
    {
      Assert(weights_handle != dealii::numbers::invalid_unsigned_int, dealii::ExcInternalError());
 
      std::vector<unsigned int> new_weights(tria->n_active_cells(), 1);
 
      const_cast<dealii::Triangulation<dim> &>(*tria).notify_ready_to_unpack(
        weights_handle,
        [&new_weights](const CellIterator                                             &cell,
                       const dealii::CellStatus                                        status,
                       const boost::iterator_range<std::vector<char>::const_iterator> &data) {
          const auto weight = dealii::Utilities::unpack<unsigned int>(data.begin(),
                                                                      data.end(),
                                                                      /*allow_compression=*/false);
 
          // the children of a refined cell inherit the weight of their parent
          if (status == dealii::CellStatus::cell_will_be_refined)
            for (const auto &child : cell->child_iterators())
              new_weights[child->active_cell_index()] = weight;
          else
            new_weights[cell->active_cell_index()] = weight;
        });
 
      weights.swap(new_weights);
      weights_handle = dealii::numbers::invalid_unsigned_int;
    }
 
    unsigned int
    get_partition_weight(const dealii::NonMatching::LocationToLevelSet location) const
    {
      switch (location)
        {
          case dealii::NonMatching::LocationToLevelSet::inside:
            return repartition_data.weights.inside;
 
          case dealii::NonMatching::LocationToLevelSet::outside:
            return repartition_data.weights.outside;
 
          case dealii::NonMatching::LocationToLevelSet::intersected:
            return repartition_data.weights.intersected;
 
          default:
            AssertThrow(false, dealii::ExcMessage("Unknown location of the cell."));
        }
    }
  };
 
  template <int dim, typename VectorType>
  bool
  repartition_triangulation(
    [[maybe_unused]] dealii::Triangulation<dim>                            &tria,
    [[maybe_unused]] const AttachDoFHandlerAndVectorsType<dim, VectorType> &attach_vectors,
    [[maybe_unused]] const std::function<void()>   &distribute_level_set_dofs,
    [[maybe_unused]] const dealii::DoFHandler<dim> &level_set_dof_handler,
    [[maybe_unused]] VectorType                    &level_set_dof_vector,
    [[maybe_unused]] const std::function<void()>   &setup_dof_system,
    [[maybe_unused]] const std::function<void()>   &post = {})
  {
    if constexpr (dim == 1)
      {
        // parallel::distributed::Triangulation is not available for dim = 1
        return false;
      }
    else
      {
        // check that the triangulation is of type parallel::distributed::Triangulation
        auto *distributed_tria =
          dynamic_cast<dealii::parallel::distributed::Triangulation<dim> *>(&tria);
 
        if (distributed_tria == nullptr)
          return false;
 
        // 1. prepare the solution transfer
        DoFHandlerAndVectorDataType<dim, VectorType> data;
        if (attach_vectors)
          {
            attach_vectors(data);
            data.shrink_to_fit();
          }
 
        const unsigned int n_dof_handlers = data.size();
 
        // create a solution transfer for all DoFHandlers
        std::vector<std::unique_ptr<dealii::SolutionTransfer<dim, VectorType>>> solution_transfers(
          n_dof_handlers);
        std::vector<std::vector<const VectorType *>> old_vectors(n_dof_handlers);
        std::vector<std::vector<bool>>               had_ghost_values(n_dof_handlers);
 
        for (unsigned int j = 0; j < n_dof_handlers; ++j)
          {
            std::vector<VectorType *> vectors;
            data[j].second(vectors);
 
            for (const auto &v : vectors)
              {
                had_ghost_values[j].push_back(v->has_ghost_elements());
                // ghost values are required to read the DoF values of the locally owned cells
                v->update_ghost_values();
                old_vectors[j].push_back(v);
              }
 
            solution_transfers[j] =
              std::make_unique<dealii::SolutionTransfer<dim, VectorType>>(*data[j].first);
            solution_transfers[j]->prepare_for_coarsening_and_refinement(old_vectors[j]);
          }
 
        level_set_dof_vector.update_ghost_values();
        dealii::SolutionTransfer<dim, VectorType> level_set_transfer(level_set_dof_handler);
        level_set_transfer.prepare_for_coarsening_and_refinement(level_set_dof_vector);
 
        // 2. repartition the triangulation according to the connected cell weights; the active FE
        // indices of hp-DoFHandlers are transferred automatically
        distributed_tria->repartition();
 
        // 3. transfer the level set first, so that the cut operations can classify the cells
        distribute_level_set_dofs();
        level_set_dof_vector.zero_out_ghost_values();
        level_set_transfer.interpolate(level_set_dof_vector);
        level_set_dof_vector.update_ghost_values();
 
        // 4. update dof-related data to match the new partitioning
        setup_dof_system();
 
        // 5. transfer the attached vectors; collect the target vectors again, since they might have
        // been reinitialized/reallocated in setup_dof_system()
        if (n_dof_handlers > 0)
          {
            DoFHandlerAndVectorDataType<dim, VectorType> new_data;
            attach_vectors(new_data);
 
            for (unsigned int j = 0; j < n_dof_handlers; ++j)
              {
                std::vector<VectorType *> new_vectors;
                new_data[j].second(new_vectors);
 
                for (const auto &v : new_vectors)
                  v->zero_out_ghost_values();
 
                solution_transfers[j]->interpolate(new_vectors);
 
                for (unsigned int i = 0; i < new_vectors.size(); ++i)
                  if (had_ghost_values[j][i])
                    new_vectors[i]->update_ghost_values();
              }
          }
 
        // 6. post-processing
        if (post)
          post();
 
        return true;
      }
  }
 
  inline bool
  do_repartition(const unsigned int n_time_step, const CutRepartitionData &repartition_data)
  {
    return repartition_data.enable and (n_time_step % repartition_data.every_n_step == 0);
  }
} // namespace MeltPoolDG::CutUtil