dof_tools.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 2016 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#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/table.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/point.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/sparsity_pattern.h>
#include <deal.II/dofs/function_map.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/fe.h>
#include <deal.II/fe/component_mask.h>
#include <deal.II/hp/mapping_collection.h>

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

DEAL_II_NAMESPACE_OPEN

template<int dim, class T> class Table;
class SparsityPattern;
template <typename number> class Vector;
template <int dim, typename Number> class Function;
template <int dim, int spacedim> class FiniteElement;
template <int dim, int spacedim> class DoFHandler;
namespace hp
{
  template <int dim, int spacedim> class DoFHandler;
  template <int dim, int spacedim> class MappingCollection;
}
class ConstraintMatrix;
template <class MeshType> class InterGridMap;
template <int dim, int spacedim> class Mapping;

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

//TODO: map_support_points_to_dofs should generate a multimap, rather than just a map, since several dofs may be located at the same support point

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

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

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


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

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

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

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

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

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

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

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

  template <typename DoFHandlerType, typename SparsityPatternType>
  void
  make_sparsity_pattern (const DoFHandlerType      &dof_handler,
                         SparsityPatternType       &sparsity_pattern,
                         const ConstraintMatrix    &constraints           = ConstraintMatrix(),
                         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_handler,
                         const Table<2, Coupling> &coupling,
                         SparsityPatternType      &sparsity_pattern,
                         const ConstraintMatrix   &constraints           = ConstraintMatrix(),
                         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>
  void
  make_flux_sparsity_pattern (const DoFHandlerType      &dof_handler,
                              SparsityPatternType       &sparsity_pattern,
                              const ConstraintMatrix    &constraints,
                              const bool                 keep_constrained_dofs = true,
                              const types::subdomain_id  subdomain_id          = numbers::invalid_unsigned_int);

  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);

  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>
  void
  make_boundary_sparsity_pattern
  (const DoFHandlerType                                              &dof,
   const typename FunctionMap<DoFHandlerType::space_dimension>::type &boundary_ids,
   const std::vector<types::global_dof_index>                        &dof_to_boundary_mapping,
   SparsityPatternType                                               &sparsity);

  template <typename DoFHandlerType>
  void
  make_hanging_node_constraints (const DoFHandlerType &dof_handler,
                                 ConstraintMatrix     &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,
                                 ConstraintMatrix                              &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,
   ::ConstraintMatrix                    &constraint_matrix,
   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,
   ::ConstraintMatrix        &constraint_matrix,
   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,
   ::ConstraintMatrix &constraint_matrix,
   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,
   ::ConstraintMatrix &constraint_matrix,
   const ComponentMask      &component_mask = ComponentMask());

  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 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>
  void
  extract_constant_modes (const DoFHandlerType            &dof_handler,
                          const ComponentMask             &component_mask,
                          std::vector<std::vector<bool> > &constant_modes);

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

  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>
  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, class Sparsity>
  void make_cell_patches(Sparsity                &block_list,
                         const DoFHandlerType    &dof_handler,
                         const unsigned int       level,
                         const std::vector<bool> &selected_dofs = std::vector<bool>(),
                         types::global_dof_index  offset        = 0);

  template <typename DoFHandlerType>
  void 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);

  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>
  std::vector<types::global_dof_index>
  get_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 <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, template <int, int> class DoFHandlerType>
  void
  make_zero_boundary_constraints (const DoFHandlerType<dim,spacedim> &dof,
                                  const types::boundary_id            boundary_id,
                                  ConstraintMatrix                   &zero_boundary_constraints,
                                  const ComponentMask                &component_mask = ComponentMask());

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


  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);
  DeclException0 (ExcFiniteElementsDontMatch);
  DeclException0 (ExcGridNotCoarser);
  DeclException0 (ExcGridsDontMatch);
  DeclException0 (ExcNoFESelected);
  DeclException0 (ExcInvalidBoundaryIndicator);
}



/* ------------------------- 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().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().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().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;
  }
}

#endif

DEAL_II_NAMESPACE_CLOSE

#endif