affine_constraints.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2023 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_affine_constraints_h
#define dealii_affine_constraints_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/exceptions.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/table.h>
#include <deal.II/base/template_constraints.h>
#include <deal.II/base/thread_local_storage.h>
 
#include <deal.II/lac/sparsity_pattern_base.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/vector_element_access.h>
 
#include <boost/range/iterator_range.hpp>
 
#include <algorithm>
#include <set>
#include <type_traits>
#include <utility>
#include <vector>
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declarations
#ifndef DOXYGEN
template <typename>
class FullMatrix;
class SparsityPattern;
class DynamicSparsityPattern;
class BlockSparsityPattern;
class BlockDynamicSparsityPattern;
template <typename number>
class SparseMatrix;
template <typename number>
class BlockSparseMatrix;
 
namespace internal
{
  namespace AffineConstraints
  {
    using size_type = types::global_dof_index;
 
    struct Distributing
    {
      Distributing(const size_type global_row = numbers::invalid_size_type,
                   const size_type local_row  = numbers::invalid_size_type);
 
      Distributing(const Distributing &in);
 
      Distributing &
      operator=(const Distributing &in);
 
      bool
      operator<(const Distributing &in) const
      {
        return global_row < in.global_row;
      }
 
      size_type         global_row;
      size_type         local_row;
      mutable size_type constraint_position;
    };
 
 
 
    template <typename number>
    struct DataCache
    {
      DataCache();
 
      void
      reinit();
 
      size_type
      insert_new_index(const std::pair<size_type, number> &pair);
 
      void
      append_index(const size_type                     index,
                   const std::pair<size_type, number> &pair);
 
      size_type
      get_size(const size_type index) const;
 
      const std::pair<size_type, number> *
      get_entry(const size_type index) const;
 
      size_type row_length;
 
      std::vector<std::pair<size_type, number>> data;
 
      std::vector<size_type> individual_size;
    };
 
 
 
    template <typename number>
    class GlobalRowsFromLocal
    {
    public:
      GlobalRowsFromLocal();
 
      void
      reinit(const size_type n_local_rows);
 
      void
      insert_index(const size_type global_row,
                   const size_type local_row,
                   const number    constraint_value);
      void
      sort();
 
      void
      print(std::ostream &os);
 
      size_type
      size() const;
 
      size_type
      size(const size_type counter_index) const;
 
      size_type
      global_row(const size_type counter_index) const;
 
      size_type &
      global_row(const size_type counter_index);
 
      size_type
      local_row(const size_type counter_index) const;
 
      size_type &
      local_row(const size_type counter_index);
 
      size_type
      local_row(const size_type counter_index,
                const size_type index_in_constraint) const;
 
      number
      constraint_value(const size_type counter_index,
                       const size_type index_in_constraint) const;
 
      bool
      have_indirect_rows() const;
 
      void
      insert_constraint(const size_type constrained_local_dof);
 
      size_type
      n_constraints() const;
 
      size_type
      n_inhomogeneities() const;
 
      void
      set_ith_constraint_inhomogeneous(const size_type i);
 
      size_type
      constraint_origin(size_type i) const;
 
      std::vector<Distributing> total_row_indices;
 
    private:
      DataCache<number> data_cache;
 
      size_type n_active_rows;
 
      size_type n_inhomogeneous_rows;
    };
 
 
 
    template <typename number>
    struct ScratchData
    {
      ScratchData()
        : in_use(false)
      {}
 
      ScratchData(const ScratchData &)
        : in_use(false)
      {}
 
      bool in_use;
 
      std::vector<std::pair<size_type, size_type>> new_entries;
 
      std::vector<size_type> rows;
 
      std::vector<size_type> columns;
 
      std::vector<number> values;
 
      std::vector<size_type> block_starts;
 
      std::vector<size_type> vector_indices;
 
      std::vector<number> vector_values;
 
      GlobalRowsFromLocal<number> global_rows;
 
      GlobalRowsFromLocal<number> global_columns;
    };
  } // namespace AffineConstraints
} // namespace internal
 
namespace internal
{
  namespace AffineConstraintsImplementation
  {
    template <class VectorType>
    void
    set_zero_all(const std::vector<types::global_dof_index> &cm,
                 VectorType &                                vec);
 
    template <class T>
    void
    set_zero_all(const std::vector<types::global_dof_index> &cm,
                 ::Vector<T> &                         vec);
 
    template <class T>
    void
    set_zero_all(const std::vector<types::global_dof_index> &cm,
                 ::BlockVector<T> &                    vec);
  } // namespace AffineConstraintsImplementation
} // namespace internal
#endif
 
// TODO[WB]: We should have a function of the kind
//   AffineConstraints::add_constraint (const size_type constrained_dof,
//     const std::vector<std::pair<size_type, number> > &entries,
//     const number inhomogeneity = 0);
// rather than building up constraints piecemeal through add_line/add_entry
// etc. This would also eliminate the possibility of accidentally changing
// existing constraints into something pointless, see the discussion on the
// mailing list on "Tiny bug in interpolate_boundary_values" in Sept. 2010.
 
template <typename number = double>
class AffineConstraints : public Subscriptor
{
public:
  using size_type = types::global_dof_index;
 
  enum MergeConflictBehavior
  {
    no_conflicts_allowed,
 
    left_object_wins,
 
    right_object_wins
  };
 
  explicit AffineConstraints(const IndexSet &local_constraints = IndexSet());
 
  explicit AffineConstraints(const AffineConstraints &affine_constraints);
 
  AffineConstraints(AffineConstraints &&affine_constraints) noexcept =
    default; // NOLINT
 
  AffineConstraints &
  operator=(const AffineConstraints &) = delete;
 
  AffineConstraints &
  operator=(AffineConstraints &&affine_constraints) noexcept =
    default; // NOLINT
 
  template <typename other_number>
  void
  copy_from(const AffineConstraints<other_number> &other);
 
  void
  reinit(const IndexSet &local_constraints = IndexSet());
 
  bool
  can_store_line(const size_type line_n) const;
 
  const IndexSet &
  get_local_lines() const;
 
  void
  add_selected_constraints(const AffineConstraints &constraints_in,
                           const IndexSet &         filter);
 
  void
  add_line(const size_type line_n);
 
  void
  add_lines(const std::vector<bool> &lines);
 
  void
  add_lines(const std::set<size_type> &lines);
 
  void
  add_lines(const IndexSet &lines);
 
  void
  add_entry(const size_type constrained_dof_index,
            const size_type column,
            const number    weight);
 
  void
  add_entries(
    const size_type                                  constrained_dof_index,
    const std::vector<std::pair<size_type, number>> &col_weight_pairs);
 
  void
  set_inhomogeneity(const size_type constrained_dof_index, const number value);
 
  void
  close();
 
  bool
  is_closed() const;
 
  bool
  is_closed(const MPI_Comm comm) const;
 
  void
  merge(
    const AffineConstraints &   other_constraints,
    const MergeConflictBehavior merge_conflict_behavior = no_conflicts_allowed,
    const bool                  allow_different_local_lines = false);
 
  void
  shift(const size_type offset);
 
  void
  clear();
 
  size_type
  n_constraints() const;
 
  size_type
  n_identities() const;
 
  size_type
  n_inhomogeneities() const;
 
  bool
  is_constrained(const size_type line_n) const;
 
  bool
  is_identity_constrained(const size_type line_n) const;
 
  bool
  are_identity_constrained(const size_type line_n_1,
                           const size_type line_n_2) const;
 
  size_type
  max_constraint_indirections() const;
 
  bool
  is_inhomogeneously_constrained(const size_type index) const;
 
  bool
  has_inhomogeneities() const;
 
  const std::vector<std::pair<size_type, number>> *
  get_constraint_entries(const size_type line_n) const;
 
  number
  get_inhomogeneity(const size_type line_n) const;
 
  void
  print(std::ostream &out) const;
 
  void
  write_dot(std::ostream &) const;
 
  std::size_t
  memory_consumption() const;
 
  void
  resolve_indices(std::vector<types::global_dof_index> &indices) const;
 
  void
  condense(SparsityPattern &sparsity) const;
 
  void
  condense(BlockSparsityPattern &sparsity) const;
 
  void
  condense(DynamicSparsityPattern &sparsity) const;
 
  void
  condense(BlockDynamicSparsityPattern &sparsity) const;
 
  void
  condense(SparseMatrix<number> &matrix) const;
 
  void
  condense(BlockSparseMatrix<number> &matrix) const;
 
  template <class VectorType>
  void
  condense(VectorType &vec) const;
 
  template <class VectorType>
  void
  condense(const VectorType &vec_ghosted, VectorType &output) const;
 
  template <class VectorType>
  void
  condense(SparseMatrix<number> &matrix, VectorType &vector) const;
 
  template <class BlockVectorType>
  void
  condense(BlockSparseMatrix<number> &matrix, BlockVectorType &vector) const;
 
  template <class VectorType>
  void
  set_zero(VectorType &vec) const;
 
  template <class InVector, class OutVector>
  void
  distribute_local_to_global(const InVector &              local_vector,
                             const std::vector<size_type> &local_dof_indices,
                             OutVector &                   global_vector) const;
 
  template <typename VectorType>
  void
  distribute_local_to_global(const Vector<number> &        local_vector,
                             const std::vector<size_type> &local_dof_indices,
                             VectorType &                  global_vector,
                             const FullMatrix<number> &    local_matrix) const;
 
  template <typename VectorType>
  void
  distribute_local_to_global(
    const Vector<number> &        local_vector,
    const std::vector<size_type> &local_dof_indices_row,
    const std::vector<size_type> &local_dof_indices_col,
    VectorType &                  global_vector,
    const FullMatrix<number> &    local_matrix,
    bool                          diagonal = false) const;
 
  template <class VectorType>
  void
  distribute_local_to_global(const size_type index,
                             const number    value,
                             VectorType &    global_vector) const;
 
  template <typename ForwardIteratorVec,
            typename ForwardIteratorInd,
            class VectorType>
  void
  distribute_local_to_global(ForwardIteratorVec local_vector_begin,
                             ForwardIteratorVec local_vector_end,
                             ForwardIteratorInd local_indices_begin,
                             VectorType &       global_vector) const;
 
  template <typename MatrixType>
  void
  distribute_local_to_global(const FullMatrix<number> &    local_matrix,
                             const std::vector<size_type> &local_dof_indices,
                             MatrixType &                  global_matrix) const;
 
  template <typename MatrixType>
  void
  distribute_local_to_global(const FullMatrix<number> &    local_matrix,
                             const std::vector<size_type> &row_indices,
                             const std::vector<size_type> &col_indices,
                             MatrixType &                  global_matrix) const;
 
  template <typename MatrixType>
  void
  distribute_local_to_global(const FullMatrix<number> &    local_matrix,
                             const std::vector<size_type> &row_indices,
                             const AffineConstraints &column_affine_constraints,
                             const std::vector<size_type> &column_indices,
                             MatrixType &                  global_matrix) const;
 
  template <typename MatrixType, typename VectorType>
  void
  distribute_local_to_global(const FullMatrix<number> &    local_matrix,
                             const Vector<number> &        local_vector,
                             const std::vector<size_type> &local_dof_indices,
                             MatrixType &                  global_matrix,
                             VectorType &                  global_vector,
                             bool use_inhomogeneities_for_rhs = false) const;
 
  void
  add_entries_local_to_global(
    const std::vector<size_type> &local_dof_indices,
    SparsityPatternBase &         sparsity_pattern,
    const bool                    keep_constrained_entries = true,
    const Table<2, bool> &        dof_mask = Table<2, bool>()) const;
 
  void
  add_entries_local_to_global(
    const std::vector<size_type> &row_indices,
    const std::vector<size_type> &col_indices,
    SparsityPatternBase &         sparsity_pattern,
    const bool                    keep_constrained_entries = true,
    const Table<2, bool> &        dof_mask = Table<2, bool>()) const;
 
  void
  add_entries_local_to_global(
    const std::vector<size_type> &   row_indices,
    const AffineConstraints<number> &col_constraints,
    const std::vector<size_type> &   col_indices,
    SparsityPatternBase &            sparsity_pattern,
    const bool                       keep_constrained_entries = true,
    const Table<2, bool> &           dof_mask = Table<2, bool>()) const;
 
  template <typename ForwardIteratorVec,
            typename ForwardIteratorInd,
            class VectorType>
  void
  get_dof_values(const VectorType & global_vector,
                 ForwardIteratorInd local_indices_begin,
                 ForwardIteratorVec local_vector_begin,
                 ForwardIteratorVec local_vector_end) const;
 
  template <class VectorType>
  void
  distribute(VectorType &vec) const;
 
  struct ConstraintLine
  {
    using Entries = std::vector<std::pair<size_type, number>>;
 
    size_type index;
 
    Entries entries;
 
    number inhomogeneity;
 
    ConstraintLine(const size_type &index         = numbers::invalid_dof_index,
                   const Entries &  entries       = {},
                   const number &   inhomogeneity = 0.0);
 
    template <typename ConstraintLineType>
    ConstraintLine(const ConstraintLineType &other);
 
    template <typename ConstraintLineType>
    ConstraintLine &
    operator=(const ConstraintLineType &other);
 
    bool
    operator<(const ConstraintLine &) const;
 
    bool
    operator==(const ConstraintLine &) const;
 
    std::size_t
    memory_consumption() const;
 
    template <class Archive>
    void
    serialize(Archive &ar, const unsigned int)
    {
      ar &index &entries &inhomogeneity;
    }
  };
 
  using const_iterator = typename std::vector<ConstraintLine>::const_iterator;
 
  using LineRange = boost::iterator_range<const_iterator>;
 
  const LineRange
  get_lines() const;
 
  bool
  is_consistent_in_parallel(const std::vector<IndexSet> &locally_owned_dofs,
                            const IndexSet &             locally_active_dofs,
                            const MPI_Comm               mpi_communicator,
                            const bool                   verbose = false) const;
 
  void
  make_consistent_in_parallel(const IndexSet &locally_owned_dofs,
                              const IndexSet &locally_relevant_dofs,
                              const MPI_Comm  mpi_communicator);
 
  DeclException0(ExcMatrixIsClosed);
  DeclException0(ExcMatrixNotClosed);
  DeclException1(ExcLineInexistant,
                 size_type,
                 << "The specified line " << arg1 << " does not exist.");
  DeclException4(ExcEntryAlreadyExists,
                 size_type,
                 size_type,
                 number,
                 number,
                 << "The entry for the indices " << arg1 << " and " << arg2
                 << " already exists, but the values " << arg3 << " (old) and "
                 << arg4 << " (new) differ "
                 << "by " << (arg4 - arg3) << '.');
  DeclException2(ExcDoFConstrainedToConstrainedDoF,
                 int,
                 int,
                 << "You tried to constrain DoF " << arg1 << " to DoF " << arg2
                 << ", but that one is also constrained. This is not allowed!");
  DeclException1(ExcDoFIsConstrainedFromBothObjects,
                 size_type,
                 << "Degree of freedom " << arg1
                 << " is constrained from both object in a merge operation.");
  DeclException1(ExcDoFIsConstrainedToConstrainedDoF,
                 size_type,
                 << "In the given argument a degree of freedom is constrained "
                 << "to another DoF with number " << arg1
                 << ", which however is constrained by this object. This is not"
                 << " allowed.");
  DeclException1(ExcRowNotStoredHere,
                 size_type,
                 << "The index set given to this constraints object indicates "
                 << "constraints for degree of freedom " << arg1
                 << " should not be stored by this object, but a constraint "
                 << "is being added.");
 
  DeclException2(ExcColumnNotStoredHere,
                 size_type,
                 size_type,
                 << "The index set given to this constraints object indicates "
                 << "constraints using degree of freedom " << arg2
                 << " should not be stored by this object, but a constraint "
                 << "for degree of freedom " << arg1 << " uses it.");
 
  DeclException2(ExcIncorrectConstraint,
                 int,
                 int,
                 << "While distributing the constraint for DoF " << arg1
                 << ", it turns out that one of the processors "
                 << "who own the " << arg2 << " degrees of freedom that x_"
                 << arg1 << " is constrained against does not know about "
                 << "the constraint on x_" << arg1
                 << ". Did you not initialize the AffineConstraints container "
                 << "with the appropriate locally_relevant set so "
                 << "that every processor who owns a DoF that constrains "
                 << "another DoF also knows about this constraint?");
 
  template <typename>
  friend class AffineConstraints;
 
private:
  std::vector<ConstraintLine> lines;
 
  std::vector<size_type> lines_cache;
 
  IndexSet local_lines;
 
  bool sorted;
 
  mutable Threads::ThreadLocalStorage<
    internal::AffineConstraints::ScratchData<number>>
    scratch_data;
 
  size_type
  calculate_line_index(const size_type line_n) const;
 
  template <typename MatrixType, typename VectorType>
  void
  distribute_local_to_global(const FullMatrix<number> &    local_matrix,
                             const Vector<number> &        local_vector,
                             const std::vector<size_type> &local_dof_indices,
                             MatrixType &                  global_matrix,
                             VectorType &                  global_vector,
                             const bool use_inhomogeneities_for_rhs,
                             const std::integral_constant<bool, false>) const;
 
  template <typename MatrixType, typename VectorType>
  void
  distribute_local_to_global(const FullMatrix<number> &    local_matrix,
                             const Vector<number> &        local_vector,
                             const std::vector<size_type> &local_dof_indices,
                             MatrixType &                  global_matrix,
                             VectorType &                  global_vector,
                             const bool use_inhomogeneities_for_rhs,
                             const std::integral_constant<bool, true>) const;
 
  void
  make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
                       internal::AffineConstraints::GlobalRowsFromLocal<number>
                         &global_rows) const;
 
  void
  make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
                       std::vector<size_type> &      active_dofs) const;
 
  template <typename MatrixScalar, typename VectorScalar>
  typename ProductType<VectorScalar, MatrixScalar>::type
  resolve_vector_entry(
    const size_type                                                 i,
    const internal::AffineConstraints::GlobalRowsFromLocal<number> &global_rows,
    const Vector<VectorScalar> &    local_vector,
    const std::vector<size_type> &  local_dof_indices,
    const FullMatrix<MatrixScalar> &local_matrix) const;
};
 
/* ---------------- template and inline functions ----------------- */
 
template <typename number>
inline AffineConstraints<number>::AffineConstraints(
  const IndexSet &local_constraints)
  : lines()
  , local_lines(local_constraints)
  , sorted(false)
{
  // make sure the IndexSet is compressed. Otherwise this can lead to crashes
  // that are hard to find (only happen in release mode).
  // see tests/mpi/affine_constraints_crash_01
  local_lines.compress();
}
 
template <typename number>
inline AffineConstraints<number>::AffineConstraints(
  const AffineConstraints &affine_constraints)
  : Subscriptor()
  , lines(affine_constraints.lines)
  , lines_cache(affine_constraints.lines_cache)
  , local_lines(affine_constraints.local_lines)
  , sorted(affine_constraints.sorted)
{}
 
template <typename number>
inline void
AffineConstraints<number>::add_line(const size_type line_n)
{
  Assert(sorted == false, ExcMatrixIsClosed());
 
  // the following can happen when we compute with distributed meshes and dof
  // handlers and we constrain a degree of freedom whose number we don't have
  // locally. if we don't abort here the program will try to allocate several
  // terabytes of memory to resize the various arrays below :-)
  Assert(line_n != numbers::invalid_size_type, ExcInternalError());
  const size_type line_index = calculate_line_index(line_n);
 
  // check whether line already exists; it may, in which case we can just quit
  if (is_constrained(line_n))
    return;
 
  // if necessary enlarge vector of existing entries for cache
  if (line_index >= lines_cache.size())
    lines_cache.resize(std::max(2 * static_cast<size_type>(lines_cache.size()),
                                line_index + 1),
                       numbers::invalid_size_type);
 
  // push a new line to the end of the list
  lines.emplace_back();
  lines.back().index         = line_n;
  lines.back().inhomogeneity = 0.;
  lines_cache[line_index]    = lines.size() - 1;
}
 
 
 
template <typename number>
inline void
AffineConstraints<number>::add_entry(const size_type constrained_dof_index,
                                     const size_type column,
                                     const number    weight)
{
  Assert(sorted == false, ExcMatrixIsClosed());
  Assert(constrained_dof_index != column,
         ExcMessage("Can't constrain a degree of freedom to itself"));
 
  // Ensure that the current line is present in the cache:
  const size_type line_index = calculate_line_index(constrained_dof_index);
  Assert(line_index < lines_cache.size(),
         ExcMessage("The current AffineConstraints does not contain the line "
                    "for the current entry. Call AffineConstraints::add_line "
                    "before calling this function."));
 
  // if in debug mode, check whether an entry for this column already exists
  // and if it's the same as the one entered at present
  //
  // in any case: exit the function if an entry for this column already
  // exists, since we don't want to enter it twice
  Assert(lines_cache[line_index] != numbers::invalid_size_type,
         ExcInternalError());
  Assert(!local_lines.size() || local_lines.is_element(column),
         ExcColumnNotStoredHere(constrained_dof_index, column));
  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
  Assert(line_ptr->index == constrained_dof_index, ExcInternalError());
  for (const auto &p : line_ptr->entries)
    if (p.first == column)
      {
        Assert(std::abs(p.second - weight) < 1.e-14,
               ExcEntryAlreadyExists(
                 constrained_dof_index, column, p.second, weight));
        return;
      }
 
  line_ptr->entries.emplace_back(column, weight);
}
 
 
 
template <typename number>
inline void
AffineConstraints<number>::set_inhomogeneity(
  const size_type constrained_dof_index,
  const number    value)
{
  const size_type line_index = calculate_line_index(constrained_dof_index);
  Assert(line_index < lines_cache.size() &&
           lines_cache[line_index] != numbers::invalid_size_type,
         ExcMessage("call add_line() before calling set_inhomogeneity()"));
  Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
  line_ptr->inhomogeneity  = value;
}
 
 
 
template <typename number>
template <class VectorType>
inline void
AffineConstraints<number>::set_zero(VectorType &vec) const
{
  // since lines is a private member, we cannot pass it to the functions
  // above. therefore, copy the content which is cheap
  std::vector<size_type> constrained_lines(lines.size());
  for (unsigned int i = 0; i < lines.size(); ++i)
    constrained_lines[i] = lines[i].index;
  internal::AffineConstraintsImplementation::set_zero_all(constrained_lines,
                                                          vec);
}
 
template <typename number>
inline types::global_dof_index
AffineConstraints<number>::n_constraints() const
{
  return lines.size();
}
 
template <typename number>
inline types::global_dof_index
AffineConstraints<number>::n_identities() const
{
  return std::count_if(lines.begin(),
                       lines.end(),
                       [](const ConstraintLine &line) {
                         return (line.entries.size() == 1) &&
                                (line.entries[0].second == number(1.));
                       });
}
 
template <typename number>
inline types::global_dof_index
AffineConstraints<number>::n_inhomogeneities() const
{
  return std::count_if(lines.begin(),
                       lines.end(),
                       [](const ConstraintLine &line) {
                         return (line.inhomogeneity != number(0.));
                       });
}
 
template <typename number>
inline bool
AffineConstraints<number>::is_constrained(const size_type index) const
{
  if (lines.empty())
    return false;
 
  const size_type line_index = calculate_line_index(index);
  return ((line_index < lines_cache.size()) &&
          (lines_cache[line_index] != numbers::invalid_size_type));
}
 
template <typename number>
inline bool
AffineConstraints<number>::is_inhomogeneously_constrained(
  const size_type line_n) const
{
  // check whether the entry is constrained. could use is_constrained, but
  // that means computing the line index twice
  const size_type line_index = calculate_line_index(line_n);
  if (line_index >= lines_cache.size() ||
      lines_cache[line_index] == numbers::invalid_size_type)
    return false;
  else
    {
      Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
      return (lines[lines_cache[line_index]].inhomogeneity != number(0.));
    }
}
 
template <typename number>
inline const std::vector<std::pair<types::global_dof_index, number>> *
AffineConstraints<number>::get_constraint_entries(const size_type line_n) const
{
  if (lines.empty())
    return nullptr;
 
  // check whether the entry is constrained. could use is_constrained, but
  // that means computing the line index twice
  const size_type line_index = calculate_line_index(line_n);
  if (line_index >= lines_cache.size() ||
      lines_cache[line_index] == numbers::invalid_size_type)
    return nullptr;
  else
    return &lines[lines_cache[line_index]].entries;
}
 
template <typename number>
inline number
AffineConstraints<number>::get_inhomogeneity(const size_type line_n) const
{
  // check whether the entry is constrained. could use is_constrained, but
  // that means computing the line index twice
  const size_type line_index = calculate_line_index(line_n);
  if (line_index >= lines_cache.size() ||
      lines_cache[line_index] == numbers::invalid_size_type)
    return 0;
  else
    return lines[lines_cache[line_index]].inhomogeneity;
}
 
template <typename number>
inline types::global_dof_index
AffineConstraints<number>::calculate_line_index(const size_type line_n) const
{
  // IndexSet is unused (serial case)
  if (local_lines.size() == 0)
    return line_n;
 
  Assert(local_lines.is_element(line_n), ExcRowNotStoredHere(line_n));
 
  return local_lines.index_within_set(line_n);
}
 
template <typename number>
inline bool
AffineConstraints<number>::can_store_line(size_type line_n) const
{
  return local_lines.size() == 0 || local_lines.is_element(line_n);
}
 
template <typename number>
inline const IndexSet &
AffineConstraints<number>::get_local_lines() const
{
  return local_lines;
}
 
template <typename number>
template <class VectorType>
inline void
AffineConstraints<number>::distribute_local_to_global(
  const size_type index,
  const number    value,
  VectorType &    global_vector) const
{
  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
 
  if (is_constrained(index) == false)
    global_vector(index) += value;
  else
    {
      const ConstraintLine &position =
        lines[lines_cache[calculate_line_index(index)]];
      for (size_type j = 0; j < position.entries.size(); ++j)
        global_vector(position.entries[j].first) +=
          value * position.entries[j].second;
    }
}
 
template <typename number>
template <typename ForwardIteratorVec,
          typename ForwardIteratorInd,
          class VectorType>
inline void
AffineConstraints<number>::distribute_local_to_global(
  ForwardIteratorVec local_vector_begin,
  ForwardIteratorVec local_vector_end,
  ForwardIteratorInd local_indices_begin,
  VectorType &       global_vector) const
{
  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
  for (; local_vector_begin != local_vector_end;
       ++local_vector_begin, ++local_indices_begin)
    {
      if (is_constrained(*local_indices_begin) == false)
        internal::ElementAccess<VectorType>::add(*local_vector_begin,
                                                 *local_indices_begin,
                                                 global_vector);
      else
        {
          const ConstraintLine &position =
            lines[lines_cache[calculate_line_index(*local_indices_begin)]];
          for (size_type j = 0; j < position.entries.size(); ++j)
            internal::ElementAccess<VectorType>::add(
              (*local_vector_begin) * position.entries[j].second,
              position.entries[j].first,
              global_vector);
        }
    }
}
 
template <typename number>
template <class InVector, class OutVector>
inline void
AffineConstraints<number>::distribute_local_to_global(
  const InVector &              local_vector,
  const std::vector<size_type> &local_dof_indices,
  OutVector &                   global_vector) const
{
  Assert(local_vector.size() == local_dof_indices.size(),
         ExcDimensionMismatch(local_vector.size(), local_dof_indices.size()));
  distribute_local_to_global(local_vector.begin(),
                             local_vector.end(),
                             local_dof_indices.begin(),
                             global_vector);
}
 
template <typename number>
template <typename ForwardIteratorVec,
          typename ForwardIteratorInd,
          class VectorType>
inline void
AffineConstraints<number>::get_dof_values(
  const VectorType & global_vector,
  ForwardIteratorInd local_indices_begin,
  ForwardIteratorVec local_vector_begin,
  ForwardIteratorVec local_vector_end) const
{
  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
  for (; local_vector_begin != local_vector_end;
       ++local_vector_begin, ++local_indices_begin)
    {
      if (is_constrained(*local_indices_begin) == false)
        *local_vector_begin = global_vector(*local_indices_begin);
      else
        {
          const ConstraintLine &position =
            lines[lines_cache[calculate_line_index(*local_indices_begin)]];
          typename VectorType::value_type value = position.inhomogeneity;
          for (size_type j = 0; j < position.entries.size(); ++j)
            value += (global_vector(position.entries[j].first) *
                      position.entries[j].second);
          *local_vector_begin = value;
        }
    }
}
 
// Forward declarations
#ifndef DOXYGEN
template <typename MatrixType>
class BlockMatrixBase;
template <typename SparsityPatternType>
class BlockSparsityPatternBase;
template <typename number>
class BlockSparseMatrixEZ;
 
namespace internal
{
  namespace AffineConstraints
  {
    template <typename MatrixType>
    struct IsBlockMatrix
    {
    private:
      template <typename T>
      static std::true_type
      check(const BlockMatrixBase<T> *);
 
      template <typename T>
      static std::true_type
      check(const BlockSparseMatrixEZ<T> *);
 
      static std::false_type
      check(...);
 
    public:
      static const bool value =
        std::is_same<decltype(check(std::declval<MatrixType *>())),
                     std::true_type>::value;
    };
 
    // instantiation of the static member
    template <typename MatrixType>
    const bool IsBlockMatrix<MatrixType>::value;
 
  } // namespace AffineConstraints
} // namespace internal
#endif
 
 
 
template <typename number>
template <typename other_number>
inline void
AffineConstraints<number>::copy_from(
  const AffineConstraints<other_number> &other)
{
  lines.clear();
  lines.insert(lines.begin(), other.lines.begin(), other.lines.end());
  lines_cache = other.lines_cache;
  local_lines = other.local_lines;
  sorted      = other.sorted;
}
 
 
 
template <typename number>
template <typename MatrixType>
inline void
AffineConstraints<number>::distribute_local_to_global(
  const FullMatrix<number> &    local_matrix,
  const std::vector<size_type> &local_dof_indices,
  MatrixType &                  global_matrix) const
{
  // create a dummy and hand on to the function actually implementing this
  // feature in the cm.templates.h file.
  Vector<typename MatrixType::value_type> dummy(0);
  distribute_local_to_global(
    local_matrix,
    dummy,
    local_dof_indices,
    global_matrix,
    dummy,
    false,
    std::integral_constant<
      bool,
      internal::AffineConstraints::IsBlockMatrix<MatrixType>::value>());
}
 
 
 
template <typename number>
template <typename MatrixType, typename VectorType>
inline void
AffineConstraints<number>::distribute_local_to_global(
  const FullMatrix<number> &    local_matrix,
  const Vector<number> &        local_vector,
  const std::vector<size_type> &local_dof_indices,
  MatrixType &                  global_matrix,
  VectorType &                  global_vector,
  bool                          use_inhomogeneities_for_rhs) const
{
  // enter the internal function with the respective block information set,
  // the actual implementation follows in the cm.templates.h file.
  distribute_local_to_global(
    local_matrix,
    local_vector,
    local_dof_indices,
    global_matrix,
    global_vector,
    use_inhomogeneities_for_rhs,
    std::integral_constant<
      bool,
      internal::AffineConstraints::IsBlockMatrix<MatrixType>::value>());
}
 
 
 
template <typename number>
inline AffineConstraints<number>::ConstraintLine::ConstraintLine(
  const size_type &                                                  index,
  const typename AffineConstraints<number>::ConstraintLine::Entries &entries,
  const number &inhomogeneity)
  : index(index)
  , entries(entries)
  , inhomogeneity(inhomogeneity)
{}
 
 
 
template <typename number>
template <typename ConstraintLineType>
inline AffineConstraints<number>::ConstraintLine::ConstraintLine(
  const ConstraintLineType &other)
{
  this->index = other.index;
 
  entries.clear();
  entries.insert(entries.begin(), other.entries.begin(), other.entries.end());
 
  this->inhomogeneity = other.inhomogeneity;
}
 
 
 
template <typename number>
template <typename ConstraintLineType>
inline typename AffineConstraints<number>::ConstraintLine &
AffineConstraints<number>::ConstraintLine::operator=(
  const ConstraintLineType &other)
{
  this->index = other.index;
 
  entries.clear();
  entries.insert(entries.begin(), other.entries.begin(), other.entries.end());
 
  this->inhomogeneity = other.inhomogeneity;
 
  return *this;
}
 
DEAL_II_NAMESPACE_CLOSE
 
#endif