mapping_dof_vector.h

/*  ______________________________________________________________________
 *
 *  ExaDG - High-Order Discontinuous Galerkin for the Exa-Scale
 *
 *  Copyright (C) 2021 by the ExaDG authors
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program. If not, see <https://www.gnu.org/licenses/>.
 *  ______________________________________________________________________
 */
 
#ifndef EXADG_GRID_MAPPING_DOF_VECTOR_H_
#define EXADG_GRID_MAPPING_DOF_VECTOR_H_
 
// deal.II
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/mg_level_object.h>
#include <deal.II/base/multithread_info.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q_cache.h>
#include <deal.II/lac/la_parallel_block_vector.h>
#include <deal.II/multigrid/mg_transfer_matrix_free.h>
 
// ExaDG
#include <exadg/grid/grid_data.h>
#include <exadg/solvers_and_preconditioners/multigrid/transfer.h>
 
namespace ExaDG
{
template<int dim, typename Number>
class MappingDoFVector
{
public:
  typedef dealii::LinearAlgebra::distributed::Vector<Number> VectorType;

  MappingDoFVector(unsigned int const mapping_degree)
  {
    mapping_q_cache = std::make_shared<dealii::MappingQCache<dim>>(mapping_degree);
 
    hierarchic_to_lexicographic_numbering =
      dealii::FETools::hierarchic_to_lexicographic_numbering<dim>(mapping_degree);
    lexicographic_to_hierarchic_numbering =
      dealii::Utilities::invert_permutation(hierarchic_to_lexicographic_numbering);
  }

  virtual ~MappingDoFVector()
  {
  }

  std::shared_ptr<dealii::Mapping<dim> const>
  get_mapping() const
  {
    AssertThrow(mapping_q_cache.get(),
                dealii::ExcMessage("Mapping object mapping_q_cache is not initialized."));
 
    return mapping_q_cache;
  }
 
  std::shared_ptr<dealii::MappingQCache<dim>>
  get_mapping_q_cache() const
  {
    AssertThrow(mapping_q_cache.get(),
                dealii::ExcMessage("Mapping object mapping_q_cache is not initialized."));
 
    return mapping_q_cache;
  }

  void
  fill_grid_coordinates_vector(VectorType &                    grid_coordinates,
                               dealii::DoFHandler<dim> const & dof_handler) const
  {
    AssertThrow(mapping_q_cache.get(),
                dealii::ExcMessage("Mapping object mapping_q_cache is not initialized."));
 
    // use the deformed state described by the dealii::MappingQCache object
    fill_grid_coordinates_vector(*mapping_q_cache, grid_coordinates, dof_handler);
  }

  void
  fill_grid_coordinates_vector(dealii::Mapping<dim> const &    mapping,
                               VectorType &                    grid_coordinates,
                               dealii::DoFHandler<dim> const & dof_handler) const
  {
    if(grid_coordinates.size() != dof_handler.n_dofs())
    {
      dealii::IndexSet const relevant_dofs_grid =
        dealii::DoFTools::extract_locally_relevant_dofs(dof_handler);
      grid_coordinates.reinit(dof_handler.locally_owned_dofs(),
                              relevant_dofs_grid,
                              dof_handler.get_mpi_communicator());
    }
    else
    {
      grid_coordinates = 0;
    }
 
    AssertThrow(get_element_type(dof_handler.get_triangulation()) == ElementType::Hypercube,
                dealii::ExcMessage("Only implemented for hypercube elements."));
 
    dealii::FiniteElement<dim> const & fe = dof_handler.get_fe();
 
    std::vector<std::array<unsigned int, dim>> component_to_system_index(
      fe.base_element(0).dofs_per_cell);
 
    if(fe.dofs_per_vertex > 0) // dealii::FE_Q
    {
      for(unsigned int i = 0; i < fe.dofs_per_cell; ++i)
      {
        component_to_system_index
          [hierarchic_to_lexicographic_numbering[fe.system_to_component_index(i).second]]
          [fe.system_to_component_index(i).first] = i;
      }
    }
    else // dealii::FE_DGQ
    {
      for(unsigned int i = 0; i < fe.dofs_per_cell; ++i)
      {
        component_to_system_index[fe.system_to_component_index(i).second]
                                 [fe.system_to_component_index(i).first] = i;
      }
    }
 
    // Set up dealii::FEValues with FE_Nothing and the Gauss-Lobatto quadrature to
    // reduce setup cost, as we only use the geometry information (this means
    // we need to call fe_values.reinit(cell) with Triangulation::cell_iterator
    // rather than dealii::DoFHandler::cell_iterator).
    dealii::FE_Nothing<dim> fe_nothing;
    dealii::FEValues<dim>   fe_values(mapping,
                                    fe_nothing,
                                    dealii::QGaussLobatto<dim>(fe.degree + 1),
                                    dealii::update_quadrature_points);
 
    std::vector<dealii::types::global_dof_index> dof_indices(fe.dofs_per_cell);
 
    for(auto const & cell : dof_handler.active_cell_iterators())
    {
      if(not cell->is_artificial())
      {
        fe_values.reinit(typename dealii::Triangulation<dim>::cell_iterator(cell));
        cell->get_dof_indices(dof_indices);
        for(unsigned int i = 0; i < fe_values.n_quadrature_points; ++i)
        {
          dealii::Point<dim> const point = fe_values.quadrature_point(i);
          for(unsigned int d = 0; d < dim; ++d)
          {
            if(grid_coordinates.get_partitioner()->in_local_range(
                 dof_indices[component_to_system_index[i][d]]))
            {
              grid_coordinates(dof_indices[component_to_system_index[i][d]]) = point[d];
            }
          }
        }
      }
    }
 
    grid_coordinates.update_ghost_values();
  }

  void
  initialize_mapping_from_dof_vector(std::shared_ptr<dealii::Mapping<dim> const> mapping,
                                     VectorType const &              displacement_vector,
                                     dealii::DoFHandler<dim> const & dof_handler)
  {
    AssertThrow(dealii::MultithreadInfo::n_threads() == 1, dealii::ExcNotImplemented());
 
    AssertThrow(mapping_q_cache.get(),
                dealii::ExcMessage("Mapping object mapping_q_cache is not initialized."));
 
    VectorType displacement_vector_ghosted;
    if(dof_handler.n_dofs() > 0 and displacement_vector.size() == dof_handler.n_dofs())
    {
      dealii::IndexSet const locally_relevant_dofs =
        dealii::DoFTools::extract_locally_relevant_dofs(dof_handler);
      displacement_vector_ghosted.reinit(dof_handler.locally_owned_dofs(),
                                         locally_relevant_dofs,
                                         dof_handler.get_mpi_communicator());
      displacement_vector_ghosted.copy_locally_owned_data_from(displacement_vector);
      displacement_vector_ghosted.update_ghost_values();
    }
 
    AssertThrow(get_element_type(dof_handler.get_triangulation()) == ElementType::Hypercube,
                dealii::ExcMessage("Only implemented for hypercube elements."));
 
    std::shared_ptr<dealii::FEValues<dim>> fe_values;
 
    // Set up dealii::FEValues with FE_Nothing and the Gauss-Lobatto quadrature to
    // reduce setup cost, as we only use the geometry information (this means
    // we need to call fe_values.reinit(cell) with Triangulation::cell_iterator
    // rather than dealii::DoFHandler::cell_iterator).
    dealii::FE_Nothing<dim> fe_nothing;
 
    if(mapping.get() != 0)
    {
      fe_values = std::make_shared<dealii::FEValues<dim>>(*mapping,
                                                          fe_nothing,
                                                          dealii::QGaussLobatto<dim>(
                                                            mapping_q_cache->get_degree() + 1),
                                                          dealii::update_quadrature_points);
    }
 
    // take the grid coordinates described by mapping and add deformation described by displacement
    // vector
    mapping_q_cache->initialize(
      dof_handler.get_triangulation(),
      [&](const typename dealii::Triangulation<dim>::cell_iterator & cell_tria)
        -> std::vector<dealii::Point<dim>> {
        unsigned int const scalar_dofs_per_cell =
          dealii::Utilities::pow(mapping_q_cache->get_degree() + 1, dim);
 
        std::vector<dealii::Point<dim>> grid_coordinates(scalar_dofs_per_cell);
 
        if(mapping.get() != 0)
        {
          fe_values->reinit(cell_tria);
          // extract displacement and add to original position
          for(unsigned int i = 0; i < scalar_dofs_per_cell; ++i)
          {
            grid_coordinates[i] =
              fe_values->quadrature_point(this->hierarchic_to_lexicographic_numbering[i]);
          }
        }
 
        // if this function is called with an empty dof-vector, this indicates that the
        // displacements are zero and the points do not have to be moved
        if(dof_handler.n_dofs() > 0 and displacement_vector.size() > 0 and
           cell_tria->is_active() and not(cell_tria->is_artificial()))
        {
          typename dealii::DoFHandler<dim>::cell_iterator cell(&cell_tria->get_triangulation(),
                                                               cell_tria->level(),
                                                               cell_tria->index(),
                                                               &dof_handler);
 
          dealii::FiniteElement<dim> const & fe = dof_handler.get_fe();
          AssertThrow(fe.element_multiplicity(0) == dim,
                      dealii::ExcMessage("Expected finite element with dim components."));
 
          std::vector<dealii::types::global_dof_index> dof_indices(fe.dofs_per_cell);
          cell->get_dof_indices(dof_indices);
 
          for(unsigned int i = 0; i < dof_indices.size(); ++i)
          {
            std::pair<unsigned int, unsigned int> const id = fe.system_to_component_index(i);
 
            if(fe.dofs_per_vertex > 0) // dealii::FE_Q
            {
              grid_coordinates[id.second][id.first] += displacement_vector_ghosted(dof_indices[i]);
            }
            else // dealii::FE_DGQ
            {
              grid_coordinates[this->lexicographic_to_hierarchic_numbering[id.second]][id.first] +=
                displacement_vector_ghosted(dof_indices[i]);
            }
          }
        }
 
        return grid_coordinates;
      });
  }
 
  std::vector<unsigned int> hierarchic_to_lexicographic_numbering;
  std::vector<unsigned int> lexicographic_to_hierarchic_numbering;
 
protected:
  std::shared_ptr<dealii::MappingQCache<dim>> mapping_q_cache;
};
 
 
namespace MappingTools
{
template<int dim, typename Number>
void
initialize_coarse_mappings_from_mapping_dof_vector(
  std::vector<std::shared_ptr<MappingDoFVector<dim, Number>>> & coarse_mappings,
  unsigned int const                                            degree_coarse_mappings,
  std::shared_ptr<MappingDoFVector<dim, Number> const> const &  fine_mapping,
  dealii::Triangulation<dim> const &                            triangulation)
{
  AssertThrow(dealii::MultithreadInfo::n_threads() == 1, dealii::ExcNotImplemented());
 
  std::shared_ptr<dealii::MappingQCache<dim>> mapping_q_cache = fine_mapping->get_mapping_q_cache();
  AssertThrow(mapping_q_cache.get(),
              dealii::ExcMessage("Shared pointer mapping_q_cache is invalid."));
 
  dealii::FESystem<dim>   fe(dealii::FE_Q<dim>(degree_coarse_mappings), dim);
  dealii::DoFHandler<dim> dof_handler(triangulation);
  dof_handler.distribute_dofs(fe);
  dof_handler.distribute_mg_dofs();
 
  std::shared_ptr<MappingDoFVector<dim, Number>> mapping_dof_vector_all_levels =
    std::make_shared<MappingDoFVector<dim, Number>>(degree_coarse_mappings);
 
  // fill a dof vector with grid coordinates of the fine level using degree_coarse_mappings
  typedef dealii::LinearAlgebra::distributed::Vector<Number> VectorType;
  VectorType                                                 grid_coordinates_fine_level;
 
  {
    mapping_dof_vector_all_levels->fill_grid_coordinates_vector(*fine_mapping->get_mapping(),
                                                                grid_coordinates_fine_level,
                                                                dof_handler);
  }
 
  // project the solution onto all coarse levels of the triangulation using degree_coarse_mappings
  dealii::MGLevelObject<VectorType> grid_coordinates_all_levels,
    grid_coordinates_all_levels_ghosted;
  unsigned int const n_levels = triangulation.n_global_levels();
  grid_coordinates_all_levels.resize(0, n_levels - 1);
  grid_coordinates_all_levels_ghosted.resize(0, n_levels - 1);
 
  dealii::MGTransferMatrixFree<dim, Number> transfer;
  transfer.build(dof_handler);
  transfer.interpolate_to_mg(dof_handler, grid_coordinates_all_levels, grid_coordinates_fine_level);
 
  // ghosting
  for(unsigned int level = 0; level < n_levels; level++)
  {
    dealii::IndexSet const relevant_dofs =
      dealii::DoFTools::extract_locally_relevant_level_dofs(dof_handler, level);
 
    grid_coordinates_all_levels_ghosted[level].reinit(
      dof_handler.locally_owned_mg_dofs(level),
      relevant_dofs,
      grid_coordinates_all_levels[level].get_mpi_communicator());
 
    grid_coordinates_all_levels_ghosted[level].copy_locally_owned_data_from(
      grid_coordinates_all_levels[level]);
 
    grid_coordinates_all_levels_ghosted[level].update_ghost_values();
  }
 
  // Call the initialize() function of dealii::MappingQCache, which initializes the mapping for all
  // levels according to grid_coordinates_all_levels_ghosted.
  mapping_dof_vector_all_levels->get_mapping_q_cache()->initialize(
    dof_handler.get_triangulation(),
    [&](const typename dealii::Triangulation<dim>::cell_iterator & cell_tria)
      -> std::vector<dealii::Point<dim>> {
      unsigned int const level = cell_tria->level();
 
      typename dealii::DoFHandler<dim>::cell_iterator cell(&cell_tria->get_triangulation(),
                                                           level,
                                                           cell_tria->index(),
                                                           &dof_handler);
 
      AssertThrow(fe.element_multiplicity(0) == dim,
                  dealii::ExcMessage("Expected finite element with dim components."));
 
      unsigned int const scalar_dofs_per_cell = dealii::Utilities::pow(fe.degree + 1, dim);
 
      std::vector<dealii::Point<dim>> grid_coordinates(scalar_dofs_per_cell);
 
      if(cell->level_subdomain_id() != dealii::numbers::artificial_subdomain_id)
      {
        std::vector<dealii::types::global_dof_index> dof_indices(fe.dofs_per_cell);
        cell->get_mg_dof_indices(dof_indices);
 
        for(unsigned int i = 0; i < dof_indices.size(); ++i)
        {
          std::pair<unsigned int, unsigned int> const id = fe.system_to_component_index(i);
 
          if(fe.dofs_per_vertex > 0) // dealii::FE_Q
          {
            grid_coordinates[id.second][id.first] =
              grid_coordinates_all_levels_ghosted[level](dof_indices[i]);
          }
          else // dealii::FE_DGQ
          {
            grid_coordinates[mapping_dof_vector_all_levels
                               ->lexicographic_to_hierarchic_numbering[id.second]][id.first] =
              grid_coordinates_all_levels_ghosted[level](dof_indices[i]);
          }
        }
      }
 
      return grid_coordinates;
    });
 
 
  AssertThrow(
    coarse_mappings.size() == n_levels - 1,
    dealii::ExcMessage(
      "coarse_mappings does not have correct size relative to the number of levels of the triangulation."));
 
  // Finally, let all coarse grid mappings point to the same MappingDoFVector object. Using the
  // same Mapping object for all multigrid h-levels is some form of legacy code. The class
  // dealii::MatrixFree can internally extract the coarse-level mapping information (provided
  // through the fine-level Mapping object).
  for(unsigned int h_level = 0; h_level < coarse_mappings.size(); ++h_level)
  {
    coarse_mappings[h_level] = mapping_dof_vector_all_levels;
  }
}

template<int dim, typename Number>
void
initialize_coarse_mappings_from_mapping_dof_vector(
  std::vector<std::shared_ptr<MappingDoFVector<dim, Number>>> &          coarse_mappings,
  unsigned int const                                                     degree_coarse_mappings,
  std::shared_ptr<MappingDoFVector<dim, Number> const> const &           fine_mapping,
  std::shared_ptr<dealii::Triangulation<dim> const> const &              fine_triangulation,
  std::vector<std::shared_ptr<dealii::Triangulation<dim> const>> const & coarse_triangulations)
{
  std::shared_ptr<dealii::MappingQCache<dim>> mapping_q_cache = fine_mapping->get_mapping_q_cache();
  AssertThrow(mapping_q_cache.get(),
              dealii::ExcMessage("Shared pointer mapping_q_cache is invalid."));
 
  // setup dof-handlers and constraints for all levels using degree_coarse_mappings
  dealii::FESystem<dim>                          fe(dealii::FE_Q<dim>(degree_coarse_mappings), dim);
  unsigned int const                             n_h_levels = coarse_triangulations.size() + 1;
  std::vector<dealii::DoFHandler<dim>>           dof_handlers_all_levels(n_h_levels);
  std::vector<dealii::AffineConstraints<Number>> constraints_all_levels(n_h_levels);
  for(unsigned int h_level = 0; h_level < n_h_levels; ++h_level)
  {
    if(h_level == n_h_levels - 1)
      dof_handlers_all_levels[h_level].reinit(*fine_triangulation);
    else
      dof_handlers_all_levels[h_level].reinit(*coarse_triangulations[h_level]);
    dof_handlers_all_levels[h_level].distribute_dofs(fe);
    // constraints are irrelevant for interpolation
    constraints_all_levels[h_level].close();
  }
 
  // fill a dof vector with grid coordinates of the fine level using degree_coarse_mappings
  typedef dealii::LinearAlgebra::distributed::Vector<Number> VectorType;
  VectorType                                                 grid_coordinates_fine_level;
 
  {
    std::shared_ptr<MappingDoFVector<dim, Number>> mapping_dof_vector_fine_level =
      std::make_shared<MappingDoFVector<dim, Number>>(degree_coarse_mappings);
 
    auto const & dof_handler_fine_level = dof_handlers_all_levels[n_h_levels - 1];
    mapping_dof_vector_fine_level->fill_grid_coordinates_vector(*fine_mapping->get_mapping(),
                                                                grid_coordinates_fine_level,
                                                                dof_handler_fine_level);
  }
 
  // create transfer objects
  dealii::MGLevelObject<dealii::MGTwoLevelTransfer<dim, VectorType>> transfers(0, n_h_levels - 1);
  for(unsigned int h_level = 1; h_level < n_h_levels; ++h_level)
  {
    transfers[h_level].reinit_geometric_transfer(dof_handlers_all_levels[h_level],
                                                 dof_handlers_all_levels[h_level - 1],
                                                 constraints_all_levels[h_level],
                                                 constraints_all_levels[h_level - 1]);
  }
 
  // a function that initializes the dof-vector for a given level and dof_handler
  const std::function<void(unsigned int const, VectorType &)> initialize_dof_vector =
    [&](unsigned int const h_level, VectorType & vector) {
      dealii::IndexSet const locally_relevant_dofs =
        dealii::DoFTools::extract_locally_relevant_dofs(dof_handlers_all_levels[h_level]);
      vector.reinit(dof_handlers_all_levels[h_level].locally_owned_dofs(),
                    locally_relevant_dofs,
                    dof_handlers_all_levels[h_level].get_mpi_communicator());
    };
 
  dealii::MGTransferGlobalCoarsening<dim, VectorType> mg_transfer_global_coarsening(
    transfers, initialize_dof_vector);
 
  // Transfer grid coordinates to coarser h-levels.
  // The dealii::DoFHandler object will not be used for global coarsening.
  dealii::DoFHandler<dim>           dof_handler_dummy;
  dealii::MGLevelObject<VectorType> grid_coordinates_all_levels(0, n_h_levels - 1);
  mg_transfer_global_coarsening.interpolate_to_mg(dof_handler_dummy,
                                                  grid_coordinates_all_levels,
                                                  grid_coordinates_fine_level);
 
  // initialize mapping for all coarse h-levels using the dof-vectors with grid coordinates
  AssertThrow(coarse_mappings.size() == n_h_levels - 1,
              dealii::ExcMessage(
                "coarse_mappings does not have correct size relative to coarse_triangulations."));
 
  for(unsigned int h_level = 0; h_level < coarse_mappings.size(); ++h_level)
  {
    coarse_mappings[h_level] =
      std::make_shared<MappingDoFVector<dim, Number>>(degree_coarse_mappings);
 
    // grid_coordinates_all_levels describes absolute coordinates -> use an uninitialized mapping
    // in order to interpret the grid coordinates vector as absolute coordinates and not as
    // displacements.
    std::shared_ptr<dealii::Mapping<dim> const> mapping_dummy;
    coarse_mappings[h_level]->initialize_mapping_from_dof_vector(
      mapping_dummy, grid_coordinates_all_levels[h_level], dof_handlers_all_levels[h_level]);
  }
}

template<int dim, typename Number>
void
initialize_coarse_mappings(
  std::vector<std::shared_ptr<MappingDoFVector<dim, Number>>> &          coarse_mappings,
  unsigned int const                                                     degree_coarse_mappings,
  std::shared_ptr<MappingDoFVector<dim, Number> const> const &           fine_mapping,
  std::shared_ptr<dealii::Triangulation<dim> const> const &              fine_triangulation,
  std::vector<std::shared_ptr<dealii::Triangulation<dim> const>> const & coarse_triangulations)
{
  if(fine_mapping.get())
  {
    if(coarse_triangulations.size() > 0)
    {
      MappingTools::initialize_coarse_mappings_from_mapping_dof_vector<dim, Number>(
        coarse_mappings,
        degree_coarse_mappings,
        fine_mapping,
        fine_triangulation,
        coarse_triangulations);
    }
    else
    {
      MappingTools::initialize_coarse_mappings_from_mapping_dof_vector<dim, Number>(
        coarse_mappings, degree_coarse_mappings, fine_mapping, *fine_triangulation);
    }
  }
}
 
} // namespace MappingTools
} // namespace ExaDG
 
#endif /* EXADG_GRID_MAPPING_DOF_VECTOR_H_ */