affine_constraints.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2019 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/lac/vector.h>
#include <deal.II/lac/vector_element_access.h>

#include <boost/range/iterator_range.hpp>

#include <set>
#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 internals
{
  template <typename number>
  class GlobalRowsFromLocal;
}

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


template <typename number>
class AffineConstraints;
#endif

using ConstraintMatrix DEAL_II_DEPRECATED = AffineConstraints<double>;
// Note: Unfortunately, we cannot move this compatibility alias into
// constraint_matrix.h directly. This would break a lot of user projects
// that include constraint_matrix.h transitively due to various deal.II
// headers that include the file.


// 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) = default; // NOLINT

  AffineConstraints &
  operator=(const AffineConstraints &) = delete;

  AffineConstraints &
  operator=(AffineConstraints &&affine_constraints) = default; // NOLINT

  void
  copy_from(const AffineConstraints &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 line_n, const size_type column, const number value);

  void
  add_entries(const size_type                                  line_n,
              const std::vector<std::pair<size_type, number>> &col_val_pairs);

  void
  set_inhomogeneity(const size_type line_n, const number value);

  void
  close();

  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;

  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;

  template <typename SparsityPatternType>
  void
  add_entries_local_to_global(
    const std::vector<size_type> &local_dof_indices,
    SparsityPatternType &         sparsity_pattern,
    const bool                    keep_constrained_entries = true,
    const Table<2, bool> &        dof_mask = Table<2, bool>()) const;

  template <typename SparsityPatternType>
  void
  add_entries_local_to_global(
    const std::vector<size_type> &row_indices,
    const std::vector<size_type> &col_indices,
    SparsityPatternType &         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;

    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;

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

private:
  std::vector<ConstraintLine> lines;

  std::vector<size_type> lines_cache;

  IndexSet local_lines;

  bool sorted;

  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,
                             bool use_inhomogeneities_for_rhs,
                             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,
                             bool use_inhomogeneities_for_rhs,
                             std::integral_constant<bool, true>) const;

  template <typename SparsityPatternType>
  void
  add_entries_local_to_global(const std::vector<size_type> &local_dof_indices,
                              SparsityPatternType &         sparsity_pattern,
                              const bool            keep_constrained_entries,
                              const Table<2, bool> &dof_mask,
                              std::integral_constant<bool, false>) const;

  template <typename SparsityPatternType>
  void
  add_entries_local_to_global(const std::vector<size_type> &local_dof_indices,
                              SparsityPatternType &         sparsity_pattern,
                              const bool            keep_constrained_entries,
                              const Table<2, bool> &dof_mask,
                              std::integral_constant<bool, true>) const;

  void
  make_sorted_row_list(
    const std::vector<size_type> &          local_dof_indices,
    internals::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 internals::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 line_n,
                                     const size_type column,
                                     const number    value)
{
  Assert(sorted == false, ExcMatrixIsClosed());
  Assert(line_n != 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(line_n);
  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(line_n, column));
  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
  Assert(line_ptr->index == line_n, ExcInternalError());
  for (const auto &p : line_ptr->entries)
    if (p.first == column)
      {
        Assert(std::abs(p.second - value) < 1.e-14,
               ExcEntryAlreadyExists(line_n, column, p.second, value));
        return;
      }

  line_ptr->entries.emplace_back(column, value);
}

template <typename number>
inline void
AffineConstraints<number>::set_inhomogeneity(const size_type line_n,
                                             const number    value)
{
  const size_type line_index = calculate_line_index(line_n);
  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 bool
AffineConstraints<number>::is_constrained(const size_type index) const
{
  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
{
  // 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())
    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() || 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;
        }
    }
}

template <typename MatrixType>
class BlockMatrixBase;
template <typename SparsityPatternType>
class BlockSparsityPatternBase;
template <typename number>
class BlockSparseMatrixEZ;

template <typename MatrixType>
struct IsBlockMatrix
{
private:
  struct yes_type
  {
    char c[1];
  };
  struct no_type
  {
    char c[2];
  };

  template <typename T>
  static yes_type
  check_for_block_matrix(const BlockMatrixBase<T> *);

  template <typename T>
  static yes_type
  check_for_block_matrix(const BlockSparsityPatternBase<T> *);

  template <typename T>
  static yes_type
  check_for_block_matrix(const BlockSparseMatrixEZ<T> *);

  static no_type
  check_for_block_matrix(...);

public:
  static const bool value =
    (sizeof(check_for_block_matrix(static_cast<MatrixType *>(nullptr))) ==
     sizeof(yes_type));
};

// instantiation of the static member
template <typename MatrixType>
const bool IsBlockMatrix<MatrixType>::value;

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, 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, IsBlockMatrix<MatrixType>::value>());
}

template <typename number>
template <typename SparsityPatternType>
inline void
AffineConstraints<number>::add_entries_local_to_global(
  const std::vector<size_type> &local_dof_indices,
  SparsityPatternType &         sparsity_pattern,
  const bool                    keep_constrained_entries,
  const Table<2, bool> &        dof_mask) const
{
  // enter the internal function with the respective block information set,
  // the actual implementation follows in the cm.templates.h file.
  add_entries_local_to_global(
    local_dof_indices,
    sparsity_pattern,
    keep_constrained_entries,
    dof_mask,
    std::integral_constant<bool, IsBlockMatrix<SparsityPatternType>::value>());
}

DEAL_II_NAMESPACE_CLOSE

#endif