dof_tools.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_dof_tools_h
#define dealii_dof_tools_h


#include <deal.II/base/config.h>

#include <deal.II/base/exceptions.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/point.h>

#include <deal.II/dofs/dof_handler.h>

#include <deal.II/fe/component_mask.h>

#include <deal.II/hp/dof_handler.h>

#include <deal.II/lac/affine_constraints.h>

#include <map>
#include <ostream>
#include <set>
#include <vector>


DEAL_II_NAMESPACE_OPEN

class BlockMask;
template <int dim, typename RangeNumberType>
class Function;
template <int dim, int spacedim>
class FiniteElement;
namespace hp
{
  template <int dim, int spacedim>
  class MappingCollection;
  template <int dim, int spacedim>
  class FECollection;
} // namespace hp
template <class MeshType>
class InterGridMap;
template <int dim, int spacedim>
class Mapping;
class SparsityPattern;
template <int dim, class T>
class Table;
template <typename Number>
class Vector;

namespace GridTools
{
  template <typename CellIterator>
  struct PeriodicFacePair;
}


namespace DoFTools
{
  enum Coupling
  {
    none,
    always,
    nonzero
  };

  template <int dim, int spacedim>
  void
  convert_couplings_to_blocks(const hp::DoFHandler<dim, spacedim> &dof_handler,
                              const Table<2, Coupling> &table_by_component,
                              std::vector<Table<2, Coupling>> &tables_by_block);

  template <int dim, int spacedim>
  void
  convert_couplings_to_blocks(const DoFHandler<dim, spacedim> &dof_handler,
                              const Table<2, Coupling> &table_by_component,
                              std::vector<Table<2, Coupling>> &tables_by_block);

  template <int dim, int spacedim>
  Table<2, Coupling>
  dof_couplings_from_component_couplings(
    const FiniteElement<dim, spacedim> &fe,
    const Table<2, Coupling> &          component_couplings);

  template <int dim, int spacedim>
  std::vector<Table<2, Coupling>>
  dof_couplings_from_component_couplings(
    const hp::FECollection<dim, spacedim> &fe,
    const Table<2, Coupling> &             component_couplings);
  template <int dim, int spacedim>
  DEAL_II_DEPRECATED unsigned int
  max_dofs_per_cell(const DoFHandler<dim, spacedim> &dh);

  template <int dim, int spacedim>
  DEAL_II_DEPRECATED unsigned int
  max_dofs_per_cell(const hp::DoFHandler<dim, spacedim> &dh);


  template <int dim, int spacedim>
  DEAL_II_DEPRECATED unsigned int
  max_dofs_per_face(const DoFHandler<dim, spacedim> &dh);

  template <int dim, int spacedim>
  DEAL_II_DEPRECATED unsigned int
  max_dofs_per_face(const hp::DoFHandler<dim, spacedim> &dh);

  template <int dim, int spacedim>
  DEAL_II_DEPRECATED unsigned int
  max_dofs_per_vertex(const DoFHandler<dim, spacedim> &dh);

  template <int dim, int spacedim>
  DEAL_II_DEPRECATED unsigned int
  max_dofs_per_vertex(const hp::DoFHandler<dim, spacedim> &dh);

  template <int dim, int spacedim>
  DEAL_II_DEPRECATED unsigned int
  n_components(const DoFHandler<dim, spacedim> &dh);

  template <int dim, int spacedim>
  DEAL_II_DEPRECATED unsigned int
  n_components(const hp::DoFHandler<dim, spacedim> &dh);

  template <int dim, int spacedim>
  DEAL_II_DEPRECATED bool
  fe_is_primitive(const DoFHandler<dim, spacedim> &dh);

  template <int dim, int spacedim>
  DEAL_II_DEPRECATED bool
  fe_is_primitive(const hp::DoFHandler<dim, spacedim> &dh);

  template <typename DoFHandlerType,
            typename SparsityPatternType,
            typename number = double>
  void
  make_sparsity_pattern(
    const DoFHandlerType &           dof_handler,
    SparsityPatternType &            sparsity_pattern,
    const AffineConstraints<number> &constraints = AffineConstraints<number>(),
    const bool                       keep_constrained_dofs = true,
    const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id);

  template <typename DoFHandlerType,
            typename SparsityPatternType,
            typename number = double>
  void
  make_sparsity_pattern(
    const DoFHandlerType &           dof_handler,
    const Table<2, Coupling> &       coupling,
    SparsityPatternType &            sparsity_pattern,
    const AffineConstraints<number> &constraints = AffineConstraints<number>(),
    const bool                       keep_constrained_dofs = true,
    const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id);

  template <typename DoFHandlerType, typename SparsityPatternType>
  void
  make_sparsity_pattern(const DoFHandlerType &dof_row,
                        const DoFHandlerType &dof_col,
                        SparsityPatternType & sparsity);

  template <typename DoFHandlerType, typename SparsityPatternType>
  void
  make_flux_sparsity_pattern(const DoFHandlerType &dof_handler,
                             SparsityPatternType & sparsity_pattern);

  template <typename DoFHandlerType,
            typename SparsityPatternType,
            typename number>
  void
  make_flux_sparsity_pattern(
    const DoFHandlerType &           dof_handler,
    SparsityPatternType &            sparsity_pattern,
    const AffineConstraints<number> &constraints,
    const bool                       keep_constrained_dofs = true,
    const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id);


  template <typename DoFHandlerType, typename SparsityPatternType>
  void
  make_flux_sparsity_pattern(
    const DoFHandlerType &    dof,
    SparsityPatternType &     sparsity,
    const Table<2, Coupling> &cell_integrals_mask,
    const Table<2, Coupling> &face_integrals_mask,
    const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id);


  template <typename DoFHandlerType,
            typename SparsityPatternType,
            typename number>
  void
  make_flux_sparsity_pattern(const DoFHandlerType &           dof,
                             SparsityPatternType &            sparsity,
                             const AffineConstraints<number> &constraints,
                             const bool                keep_constrained_dofs,
                             const Table<2, Coupling> &couplings,
                             const Table<2, Coupling> &face_couplings,
                             const types::subdomain_id subdomain_id);

  template <typename DoFHandlerType, typename SparsityPatternType>
  void
  make_boundary_sparsity_pattern(
    const DoFHandlerType &                      dof,
    const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
    SparsityPatternType &                       sparsity_pattern);

  template <typename DoFHandlerType,
            typename SparsityPatternType,
            typename number>
  void
  make_boundary_sparsity_pattern(
    const DoFHandlerType &dof,
    const std::map<types::boundary_id,
                   const Function<DoFHandlerType::space_dimension, number> *>
      &                                         boundary_ids,
    const std::vector<types::global_dof_index> &dof_to_boundary_mapping,
    SparsityPatternType &                       sparsity);

  template <typename DoFHandlerType, typename number>
  void
  make_hanging_node_constraints(const DoFHandlerType &     dof_handler,
                                AffineConstraints<number> &constraints);

  template <int dim, int spacedim>
  void
  compute_intergrid_constraints(
    const DoFHandler<dim, spacedim> &              coarse_grid,
    const unsigned int                             coarse_component,
    const DoFHandler<dim, spacedim> &              fine_grid,
    const unsigned int                             fine_component,
    const InterGridMap<DoFHandler<dim, spacedim>> &coarse_to_fine_grid_map,
    AffineConstraints<double> &                    constraints);


  template <int dim, int spacedim>
  void
  compute_intergrid_transfer_representation(
    const DoFHandler<dim, spacedim> &              coarse_grid,
    const unsigned int                             coarse_component,
    const DoFHandler<dim, spacedim> &              fine_grid,
    const unsigned int                             fine_component,
    const InterGridMap<DoFHandler<dim, spacedim>> &coarse_to_fine_grid_map,
    std::vector<std::map<types::global_dof_index, float>>
      &transfer_representation);

  template <typename FaceIterator>
  void
  make_periodicity_constraints(
    const FaceIterator &                         face_1,
    const typename identity<FaceIterator>::type &face_2,
    AffineConstraints<double> &                  constraints,
    const ComponentMask &            component_mask   = ComponentMask(),
    const bool                       face_orientation = true,
    const bool                       face_flip        = false,
    const bool                       face_rotation    = false,
    const FullMatrix<double> &       matrix           = FullMatrix<double>(),
    const std::vector<unsigned int> &first_vector_components =
      std::vector<unsigned int>());



  template <typename DoFHandlerType>
  void
  make_periodicity_constraints(
    const std::vector<
      GridTools::PeriodicFacePair<typename DoFHandlerType::cell_iterator>>
      &                              periodic_faces,
    AffineConstraints<double> &      constraints,
    const ComponentMask &            component_mask = ComponentMask(),
    const std::vector<unsigned int> &first_vector_components =
      std::vector<unsigned int>());



  template <typename DoFHandlerType>
  void
  make_periodicity_constraints(
    const DoFHandlerType &     dof_handler,
    const types::boundary_id   b_id1,
    const types::boundary_id   b_id2,
    const int                  direction,
    AffineConstraints<double> &constraints,
    const ComponentMask &      component_mask = ComponentMask());



  template <typename DoFHandlerType>
  void
  make_periodicity_constraints(
    const DoFHandlerType &     dof_handler,
    const types::boundary_id   b_id,
    const int                  direction,
    AffineConstraints<double> &constraints,
    const ComponentMask &      component_mask = ComponentMask());

  template <int dim, int spacedim>
  void
  extract_hanging_node_dofs(const DoFHandler<dim, spacedim> &dof_handler,
                            std::vector<bool> &              selected_dofs);

  template <int dim, int spacedim>
  IndexSet
  extract_hanging_node_dofs(const DoFHandler<dim, spacedim> &dof_handler);

  template <int dim, int spacedim>
  void
  extract_dofs(const DoFHandler<dim, spacedim> &dof_handler,
               const ComponentMask &            component_mask,
               std::vector<bool> &              selected_dofs);

  template <int dim, int spacedim>
  void
  extract_dofs(const hp::DoFHandler<dim, spacedim> &dof_handler,
               const ComponentMask &                component_mask,
               std::vector<bool> &                  selected_dofs);

  template <int dim, int spacedim>
  void
  extract_dofs(const DoFHandler<dim, spacedim> &dof_handler,
               const BlockMask &                block_mask,
               std::vector<bool> &              selected_dofs);

  template <int dim, int spacedim>
  void
  extract_dofs(const hp::DoFHandler<dim, spacedim> &dof_handler,
               const BlockMask &                    block_mask,
               std::vector<bool> &                  selected_dofs);

  template <typename DoFHandlerType>
  void
  extract_level_dofs(const unsigned int    level,
                     const DoFHandlerType &dof,
                     const ComponentMask & component_mask,
                     std::vector<bool> &   selected_dofs);

  template <typename DoFHandlerType>
  void
  extract_level_dofs(const unsigned int    level,
                     const DoFHandlerType &dof,
                     const BlockMask &     component_mask,
                     std::vector<bool> &   selected_dofs);

  template <typename DoFHandlerType>
  void
  extract_boundary_dofs(const DoFHandlerType &              dof_handler,
                        const ComponentMask &               component_mask,
                        std::vector<bool> &                 selected_dofs,
                        const std::set<types::boundary_id> &boundary_ids =
                          std::set<types::boundary_id>());

  template <typename DoFHandlerType>
  void
  extract_boundary_dofs(const DoFHandlerType &              dof_handler,
                        const ComponentMask &               component_mask,
                        IndexSet &                          selected_dofs,
                        const std::set<types::boundary_id> &boundary_ids =
                          std::set<types::boundary_id>());

  template <typename DoFHandlerType>
  void
  extract_dofs_with_support_on_boundary(
    const DoFHandlerType &              dof_handler,
    const ComponentMask &               component_mask,
    std::vector<bool> &                 selected_dofs,
    const std::set<types::boundary_id> &boundary_ids =
      std::set<types::boundary_id>());

  template <typename DoFHandlerType, typename number = double>
  IndexSet
  extract_dofs_with_support_contained_within(
    const DoFHandlerType &dof_handler,
    const std::function<
      bool(const typename DoFHandlerType::active_cell_iterator &)> &predicate,
    const AffineConstraints<number> &constraints = AffineConstraints<number>());

  template <typename DoFHandlerType>
  void
  extract_constant_modes(const DoFHandlerType &          dof_handler,
                         const ComponentMask &           component_mask,
                         std::vector<std::vector<bool>> &constant_modes);

  template <typename DoFHandlerType>
  void
  extract_subdomain_dofs(const DoFHandlerType &    dof_handler,
                         const types::subdomain_id subdomain_id,
                         std::vector<bool> &       selected_dofs);


  template <typename DoFHandlerType>
  DEAL_II_DEPRECATED void
  extract_locally_owned_dofs(const DoFHandlerType &dof_handler,
                             IndexSet &            dof_set);


  template <typename DoFHandlerType>
  void
  extract_locally_active_dofs(const DoFHandlerType &dof_handler,
                              IndexSet &            dof_set);

  template <typename DoFHandlerType>
  void
  extract_locally_relevant_dofs(const DoFHandlerType &dof_handler,
                                IndexSet &            dof_set);

  template <typename DoFHandlerType>
  std::vector<IndexSet>
  locally_owned_dofs_per_subdomain(const DoFHandlerType &dof_handler);

  template <typename DoFHandlerType>
  std::vector<IndexSet>
  locally_relevant_dofs_per_subdomain(const DoFHandlerType &dof_handler);


  template <typename DoFHandlerType>
  void
  extract_locally_relevant_level_dofs(const DoFHandlerType &dof_handler,
                                      const unsigned int    level,
                                      IndexSet &            dof_set);


  template <typename DoFHandlerType>
  void
  get_subdomain_association(const DoFHandlerType &            dof_handler,
                            std::vector<types::subdomain_id> &subdomain);

  template <typename DoFHandlerType>
  unsigned int
  count_dofs_with_subdomain_association(const DoFHandlerType &    dof_handler,
                                        const types::subdomain_id subdomain);

  template <typename DoFHandlerType>
  void
  count_dofs_with_subdomain_association(
    const DoFHandlerType &     dof_handler,
    const types::subdomain_id  subdomain,
    std::vector<unsigned int> &n_dofs_on_subdomain);

  template <typename DoFHandlerType>
  IndexSet
  dof_indices_with_subdomain_association(const DoFHandlerType &    dof_handler,
                                         const types::subdomain_id subdomain);
  // @}

  template <typename DoFHandlerType>
  std::vector<types::global_dof_index>
  get_dofs_on_patch(
    const std::vector<typename DoFHandlerType::active_cell_iterator> &patch);

  template <int dim, int spacedim>
  void
  make_cell_patches(
    SparsityPattern &                block_list,
    const DoFHandler<dim, spacedim> &dof_handler,
    const unsigned int               level,
    const std::vector<bool> &        selected_dofs = std::vector<bool>(),
    const types::global_dof_index    offset        = 0);

  template <typename DoFHandlerType>
  std::vector<unsigned int>
  make_vertex_patches(SparsityPattern &     block_list,
                      const DoFHandlerType &dof_handler,
                      const unsigned int    level,
                      const bool            interior_dofs_only,
                      const bool            boundary_patches       = false,
                      const bool            level_boundary_patches = false,
                      const bool            single_cell_patches    = false,
                      const bool            invert_vertex_mapping  = false);

  template <typename DoFHandlerType>
  std::vector<unsigned int>
  make_vertex_patches(SparsityPattern &     block_list,
                      const DoFHandlerType &dof_handler,
                      const unsigned int    level,
                      const BlockMask &     exclude_boundary_dofs = BlockMask(),
                      const bool            boundary_patches      = false,
                      const bool            level_boundary_patches = false,
                      const bool            single_cell_patches    = false,
                      const bool            invert_vertex_mapping  = false);

  template <typename DoFHandlerType>
  void
  make_child_patches(SparsityPattern &     block_list,
                     const DoFHandlerType &dof_handler,
                     const unsigned int    level,
                     const bool            interior_dofs_only,
                     const bool            boundary_dofs = false);

  template <typename DoFHandlerType>
  void
  make_single_patch(SparsityPattern &     block_list,
                    const DoFHandlerType &dof_handler,
                    const unsigned int    level,
                    const bool            interior_dofs_only = false);

  template <typename DoFHandlerType>
  void
  count_dofs_per_component(
    const DoFHandlerType &                dof_handler,
    std::vector<types::global_dof_index> &dofs_per_component,
    const bool                            vector_valued_once = false,
    std::vector<unsigned int> target_component = std::vector<unsigned int>());

  template <typename DoFHandlerType>
  void
  count_dofs_per_block(const DoFHandlerType &                dof,
                       std::vector<types::global_dof_index> &dofs_per_block,
                       const std::vector<unsigned int> &     target_block =
                         std::vector<unsigned int>());

  template <typename DoFHandlerType>
  void
  get_active_fe_indices(const DoFHandlerType &     dof_handler,
                        std::vector<unsigned int> &active_fe_indices);

  template <typename DoFHandlerType>
  unsigned int
  count_dofs_on_patch(
    const std::vector<typename DoFHandlerType::active_cell_iterator> &patch);

  template <typename DoFHandlerType>
  void
  map_dof_to_boundary_indices(const DoFHandlerType &                dof_handler,
                              std::vector<types::global_dof_index> &mapping);

  template <typename DoFHandlerType>
  void
  map_dof_to_boundary_indices(const DoFHandlerType &              dof_handler,
                              const std::set<types::boundary_id> &boundary_ids,
                              std::vector<types::global_dof_index> &mapping);

  template <int dim, int spacedim>
  void
  map_dofs_to_support_points(const Mapping<dim, spacedim> &   mapping,
                             const DoFHandler<dim, spacedim> &dof_handler,
                             std::vector<Point<spacedim>> &   support_points);

  template <int dim, int spacedim>
  void
  map_dofs_to_support_points(
    const ::hp::MappingCollection<dim, spacedim> &mapping,
    const hp::DoFHandler<dim, spacedim> &               dof_handler,
    std::vector<Point<spacedim>> &                      support_points);

  template <int dim, int spacedim>
  void
  map_dofs_to_support_points(
    const Mapping<dim, spacedim> &                      mapping,
    const DoFHandler<dim, spacedim> &                   dof_handler,
    std::map<types::global_dof_index, Point<spacedim>> &support_points);

  template <int dim, int spacedim>
  void
  map_dofs_to_support_points(
    const ::hp::MappingCollection<dim, spacedim> &mapping,
    const hp::DoFHandler<dim, spacedim> &               dof_handler,
    std::map<types::global_dof_index, Point<spacedim>> &support_points);


  template <typename DoFHandlerType, class Comp>
  void
  map_support_points_to_dofs(
    const Mapping<DoFHandlerType::dimension, DoFHandlerType::space_dimension>
      &                   mapping,
    const DoFHandlerType &dof_handler,
    std::map<Point<DoFHandlerType::space_dimension>,
             types::global_dof_index,
             Comp> &      point_to_index_map);
  template <typename DoFHandlerType, typename Number>
  void
  distribute_cell_to_dof_vector(const DoFHandlerType &dof_handler,
                                const Vector<Number> &cell_data,
                                Vector<double> &      dof_data,
                                const unsigned int    component = 0);


  template <int spacedim>
  void
  write_gnuplot_dof_support_point_info(
    std::ostream &                                            out,
    const std::map<types::global_dof_index, Point<spacedim>> &support_points);


  template <int dim,
            int spacedim,
            template <int, int> class DoFHandlerType,
            typename number>
  void
  make_zero_boundary_constraints(
    const DoFHandlerType<dim, spacedim> &dof,
    const types::boundary_id             boundary_id,
    AffineConstraints<number> &          zero_boundary_constraints,
    const ComponentMask &                component_mask = ComponentMask());

  template <int dim,
            int spacedim,
            template <int, int> class DoFHandlerType,
            typename number>
  void
  make_zero_boundary_constraints(
    const DoFHandlerType<dim, spacedim> &dof,
    AffineConstraints<number> &          zero_boundary_constraints,
    const ComponentMask &                component_mask = ComponentMask());

  DeclException0(ExcFiniteElementsDontMatch);
  DeclException0(ExcGridNotCoarser);
  DeclException0(ExcGridsDontMatch);
  DeclException0(ExcNoFESelected);
  DeclException0(ExcInvalidBoundaryIndicator);
} // namespace DoFTools



/* ------------------------- inline functions -------------- */

#ifndef DOXYGEN

namespace DoFTools
{
  inline Coupling
  operator|=(Coupling &c1, const Coupling c2)
  {
    if (c2 == always)
      c1 = always;
    else if (c1 != always && c2 == nonzero)
      return c1 = nonzero;
    return c1;
  }


  inline Coupling
  operator|(const Coupling c1, const Coupling c2)
  {
    if (c1 == always || c2 == always)
      return always;
    if (c1 == nonzero || c2 == nonzero)
      return nonzero;
    return none;
  }


  // ---------------------- inline and template functions --------------------

  template <int dim, int spacedim>
  inline unsigned int
  max_dofs_per_cell(const DoFHandler<dim, spacedim> &dh)
  {
    return dh.get_fe().dofs_per_cell;
  }


  template <int dim, int spacedim>
  inline unsigned int
  max_dofs_per_face(const DoFHandler<dim, spacedim> &dh)
  {
    return dh.get_fe().dofs_per_face;
  }


  template <int dim, int spacedim>
  inline unsigned int
  max_dofs_per_vertex(const DoFHandler<dim, spacedim> &dh)
  {
    return dh.get_fe().dofs_per_vertex;
  }


  template <int dim, int spacedim>
  inline unsigned int
  n_components(const DoFHandler<dim, spacedim> &dh)
  {
    return dh.get_fe().n_components();
  }



  template <int dim, int spacedim>
  inline bool
  fe_is_primitive(const DoFHandler<dim, spacedim> &dh)
  {
    return dh.get_fe().is_primitive();
  }


  template <int dim, int spacedim>
  inline unsigned int
  max_dofs_per_cell(const hp::DoFHandler<dim, spacedim> &dh)
  {
    return dh.get_fe_collection().max_dofs_per_cell();
  }


  template <int dim, int spacedim>
  inline unsigned int
  max_dofs_per_face(const hp::DoFHandler<dim, spacedim> &dh)
  {
    return dh.get_fe_collection().max_dofs_per_face();
  }


  template <int dim, int spacedim>
  inline unsigned int
  max_dofs_per_vertex(const hp::DoFHandler<dim, spacedim> &dh)
  {
    return dh.get_fe_collection().max_dofs_per_vertex();
  }


  template <int dim, int spacedim>
  inline unsigned int
  n_components(const hp::DoFHandler<dim, spacedim> &dh)
  {
    return dh.get_fe(0).n_components();
  }


  template <int dim, int spacedim>
  inline bool
  fe_is_primitive(const hp::DoFHandler<dim, spacedim> &dh)
  {
    return dh.get_fe(0).is_primitive();
  }


  template <typename DoFHandlerType, class Comp>
  void
  map_support_points_to_dofs(
    const Mapping<DoFHandlerType::dimension, DoFHandlerType::space_dimension>
      &                   mapping,
    const DoFHandlerType &dof_handler,
    std::map<Point<DoFHandlerType::space_dimension>,
             types::global_dof_index,
             Comp> &      point_to_index_map)
  {
    // let the checking of arguments be
    // done by the function first
    // called
    std::vector<Point<DoFHandlerType::space_dimension>> support_points(
      dof_handler.n_dofs());
    map_dofs_to_support_points(mapping, dof_handler, support_points);
    // now copy over the results of the
    // previous function into the
    // output arg
    point_to_index_map.clear();
    for (types::global_dof_index i = 0; i < dof_handler.n_dofs(); ++i)
      point_to_index_map[support_points[i]] = i;
  }
} // namespace DoFTools

#endif

DEAL_II_NAMESPACE_CLOSE

#endif