affine_constraints.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 1998 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.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 <typename 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
 
 
 
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
  };
 
  AffineConstraints();
 
  DEAL_II_DEPRECATED_EARLY_WITH_COMMENT(
    "Use the constructor with two index set arguments.")
  explicit AffineConstraints(const IndexSet &locally_stored_constraints);
 
  AffineConstraints(const IndexSet &locally_owned_dofs,
                    const IndexSet &locally_stored_constraints);
 
  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();
 
  DEAL_II_DEPRECATED_EARLY_WITH_COMMENT(
    "Use the reinit() function with two index set arguments.")
  void
  reinit(const IndexSet &locally_stored_constraints);
 
  void
  reinit(const IndexSet &locally_owned_dofs,
         const IndexSet &locally_stored_constraints);
 
  bool
  can_store_line(const size_type line_n) const;
 
  const IndexSet &
  get_locally_owned_indices() const;
 
  const IndexSet &
  get_local_lines() const;
 
  DEAL_II_DEPRECATED_EARLY_WITH_COMMENT("Use get_view() and merge() instead.")
  void
  add_selected_constraints(const AffineConstraints &constraints_in,
                           const IndexSet          &filter);
 
  void
  add_constraint(
    const size_type                                      constrained_dof,
    const ArrayView<const std::pair<size_type, number>> &dependencies,
    const number                                         inhomogeneity = 0);
 
  void
  constrain_dof_to_zero(const size_type constrained_dof);
 
  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;
 
  template <typename other_number>
  void
  merge(
    const AffineConstraints<other_number> &other_constraints,
    const MergeConflictBehavior merge_conflict_behavior = no_conflicts_allowed,
    const bool                  allow_different_local_lines = false);
 
  void
  shift(const size_type offset);
 
  AffineConstraints
  get_view(const IndexSet &mask) const;
 
  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 <typename VectorType>
  void
  condense(VectorType &vec) const;
 
  template <typename VectorType>
  void
  condense(const VectorType &vec_ghosted, VectorType &output) const;
 
  template <typename VectorType>
  void
  condense(SparseMatrix<number> &matrix, VectorType &vector) const;
 
  template <typename BlockVectorType>
  void
  condense(BlockSparseMatrix<number> &matrix, BlockVectorType &vector) const;
 
  template <typename 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 <typename VectorType>
  void
  distribute_local_to_global(const size_type index,
                             const number    value,
                             VectorType     &global_vector) const;
 
  template <typename ForwardIteratorVec,
            typename ForwardIteratorInd,
            typename 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,
            typename VectorType>
  void
  get_dof_values(const VectorType  &global_vector,
                 ForwardIteratorInd local_indices_begin,
                 ForwardIteratorVec local_vector_begin,
                 ForwardIteratorVec local_vector_end) const;
 
  template <typename 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 typename AffineConstraints<
                     number>::ConstraintLine::Entries &entries       = {},
                   const number                        inhomogeneity = 0.0);
 
    ConstraintLine(const ConstraintLine &other) = default;
 
    ConstraintLine(ConstraintLine &&other) noexcept = default;
 
    ConstraintLine &
    operator=(const ConstraintLine &other) = default;
 
    ConstraintLine &
    operator=(ConstraintLine &&other) noexcept = default;
 
    std::size_t
    memory_consumption() const;
 
    template <class Archive>
    void
    serialize(Archive &ar, const unsigned int)
    {
      ar &index &entries &inhomogeneity;
    }
 
    friend void
    swap(ConstraintLine &l1, ConstraintLine &l2) noexcept
    {
      std::swap(l1.index, l2.index);
      std::swap(l1.entries, l2.entries);
      std::swap(l1.inhomogeneity, l2.inhomogeneity);
    }
  };
 
  using const_iterator = typename std::vector<ConstraintLine>::const_iterator;
 
  using LineRange = boost::iterator_range<const_iterator>;
 
  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 &constraints_to_make_consistent,
                              const MPI_Comm  mpi_communicator);
 
  DeclExceptionMsg(
    ExcMatrixIsClosed,
    "You are attempting an operation on an AffineConstraints object "
    "that requires the object to not be 'closed', i.e., for which you "
    "must not already have called the close() member function. But the "
    "object is already closed, and so the operation can not be "
    "performed.");
  DeclExceptionMsg(
    ExcMatrixNotClosed,
    "You are attempting an operation on an AffineConstraints object "
    "that requires the object to be 'closed', i.e., for which you "
    "needed to call the close() member function. But the object "
    "is not currently closed, and so the operation can not be "
    "performed.");
  DeclException1(ExcLineInexistent,
                 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 locally_owned_dofs;
 
  IndexSet local_lines;
 
  IndexSet needed_elements_for_distribute;
 
  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::bool_constant<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::bool_constant<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()
  : AffineConstraints<number>(IndexSet(), IndexSet())
{}
 
 
 
template <typename number>
inline AffineConstraints<number>::AffineConstraints(
  const IndexSet &locally_stored_constraints)
  : AffineConstraints<number>(locally_stored_constraints,
                              locally_stored_constraints)
{}
 
 
 
template <typename number>
inline AffineConstraints<number>::AffineConstraints(
  const IndexSet &locally_owned_dofs,
  const IndexSet &locally_stored_constraints)
  : lines()
  , locally_owned_dofs(locally_owned_dofs)
  , local_lines(locally_stored_constraints)
  , sorted(false)
{
  Assert(locally_owned_dofs.is_subset_of(locally_stored_constraints),
         ExcMessage("The set of locally stored constraints needs to be a "
                    "superset of the locally owned DoFs."));
 
  // 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)
  , locally_owned_dofs(affine_constraints.locally_owned_dofs)
  , local_lines(affine_constraints.local_lines)
  , needed_elements_for_distribute(
      affine_constraints.needed_elements_for_distribute)
  , 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 <typename 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_locally_owned_indices() const
{
  return locally_owned_dofs;
}
 
 
 
template <typename number>
inline const IndexSet &
AffineConstraints<number>::get_local_lines() const
{
  return local_lines;
}
 
 
 
template <typename number>
template <typename 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,
          typename 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(global_vector.has_ghost_elements() == false, ExcGhostsPresent());
  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,
          typename 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_v<decltype(check(std::declval<MatrixType *>())),
                       std::true_type>;
    };
 
    // 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.reserve(other.lines.size());
 
  for (const auto &l : other.lines)
    lines.emplace_back(l.index,
                       typename ConstraintLine::Entries(l.entries.begin(),
                                                        l.entries.end()),
                       l.inhomogeneity);
 
  lines_cache = other.lines_cache;
  local_lines = other.local_lines;
  sorted      = other.sorted;
 
  locally_owned_dofs             = other.locally_owned_dofs;
  needed_elements_for_distribute = other.needed_elements_for_distribute;
}
 
 
template <typename number>
template <typename other_number>
void
AffineConstraints<number>::merge(
  const AffineConstraints<other_number> &other_constraints,
  const MergeConflictBehavior            merge_conflict_behavior,
  const bool                             allow_different_local_lines)
{
  (void)allow_different_local_lines;
  Assert(allow_different_local_lines ||
           local_lines == other_constraints.local_lines,
         ExcMessage(
           "local_lines for this and the other objects are not the same "
           "although allow_different_local_lines is false."));
 
  // store the previous state with respect to sorting
  const bool object_was_sorted = sorted;
  sorted                       = false;
 
  // first action is to fold into the present object possible constraints
  // in the second object. we don't strictly need to do this any more since
  // the AffineConstraints container has learned to deal with chains of
  // constraints in the close() function, but we have traditionally done
  // this and it's not overly hard to do.
  //
  // for this, loop over all constraints and replace the constraint lines
  // with a new one where constraints are replaced if necessary.
  typename ConstraintLine::Entries tmp;
  for (ConstraintLine &line : lines)
    {
      tmp.clear();
      for (const std::pair<size_type, number> &entry : line.entries)
        {
          // if the present dof is not stored, or not constrained, or if we
          // won't take the constraint from the other object, then simply copy
          // it over
          if ((other_constraints.local_lines.size() != 0. &&
               other_constraints.local_lines.is_element(entry.first) ==
                 false) ||
              other_constraints.is_constrained(entry.first) == false ||
              ((merge_conflict_behavior != right_object_wins) &&
               other_constraints.is_constrained(entry.first) &&
               this->is_constrained(entry.first)))
            tmp.push_back(entry);
          else
            // otherwise resolve further constraints by replacing the old
            // entry by a sequence of new entries taken from the other
            // object, but with multiplied weights
            {
              const auto *other_entries =
                other_constraints.get_constraint_entries(entry.first);
              Assert(other_entries != nullptr, ExcInternalError());
 
              const number weight = entry.second;
 
              for (const auto &other_entry : *other_entries)
                tmp.emplace_back(other_entry.first,
                                 other_entry.second * weight);
 
              line.inhomogeneity +=
                other_constraints.get_inhomogeneity(entry.first) * weight;
            }
        }
      // finally exchange old and newly resolved line
      line.entries.swap(tmp);
    }
 
  if (local_lines.size() != 0)
    local_lines.add_indices(other_constraints.local_lines);
 
  {
    // do not bother to resize the lines cache exactly since it is pretty
    // cheap to adjust it along the way.
    std::fill(lines_cache.begin(),
              lines_cache.end(),
              numbers::invalid_size_type);
 
    // reset lines_cache for our own constraints
    size_type index = 0;
    for (const ConstraintLine &line : lines)
      {
        const size_type local_line_no = calculate_line_index(line.index);
        if (local_line_no >= lines_cache.size())
          lines_cache.resize(local_line_no + 1, numbers::invalid_size_type);
        lines_cache[local_line_no] = index++;
      }
 
    // Add other_constraints to lines cache and our list of constraints
    for (const auto &line : other_constraints.lines)
      {
        const size_type local_line_no = calculate_line_index(line.index);
        if (local_line_no >= lines_cache.size())
          {
            lines_cache.resize(local_line_no + 1, numbers::invalid_size_type);
            lines.emplace_back(line.index,
                               typename ConstraintLine::Entries(
                                 line.entries.begin(), line.entries.end()),
                               line.inhomogeneity);
            lines_cache[local_line_no] = index++;
          }
        else if (lines_cache[local_line_no] == numbers::invalid_size_type)
          {
            // there are no constraints for that line yet
            lines.emplace_back(line.index,
                               typename ConstraintLine::Entries(
                                 line.entries.begin(), line.entries.end()),
                               line.inhomogeneity);
            AssertIndexRange(local_line_no, lines_cache.size());
            lines_cache[local_line_no] = index++;
          }
        else
          {
            // we already store that line
            switch (merge_conflict_behavior)
              {
                case no_conflicts_allowed:
                  AssertThrow(false,
                              ExcDoFIsConstrainedFromBothObjects(line.index));
                  break;
 
                case left_object_wins:
                  // ignore this constraint
                  break;
 
                case right_object_wins:
                  AssertIndexRange(local_line_no, lines_cache.size());
                  lines[lines_cache[local_line_no]] = {
                    line.index,
                    typename ConstraintLine::Entries(line.entries.begin(),
                                                     line.entries.end()),
                    static_cast<number>(line.inhomogeneity)};
                  break;
 
                default:
                  DEAL_II_NOT_IMPLEMENTED();
              }
          }
      }
 
    // check that we set the pointers correctly
    for (size_type i = 0; i < lines_cache.size(); ++i)
      if (lines_cache[i] != numbers::invalid_size_type)
        Assert(i == calculate_line_index(lines[lines_cache[i]].index),
               ExcInternalError());
  }
 
  // if the object was sorted before, then make sure it is so afterward as
  // well. otherwise leave everything in the unsorted state
  if (object_was_sorted == true)
    close();
}
 
 
 
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)
{}
 
 
 
DEAL_II_NAMESPACE_CLOSE
 
#endif