constraint_matrix.h

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


#ifndef dealii__constraint_matrix_h
#define dealii__constraint_matrix_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/template_constraints.h>
#include <deal.II/base/thread_local_storage.h>

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

#include <vector>
#include <map>
#include <set>
#include <utility>
#include <complex>


DEAL_II_NAMESPACE_OPEN

template<int dim, class T> class Table;
template <typename> class FullMatrix;
class SparsityPattern;
class DynamicSparsityPattern;
class BlockSparsityPattern;
class BlockDynamicSparsityPattern;
template <typename number> class SparseMatrix;
template <typename number> class BlockSparseMatrix;
class BlockIndices;


namespace internals
{
  class GlobalRowsFromLocal;
}


//TODO[WB]: We should have a function of the kind
//   ConstraintMatrix::add_constraint (const size_type constrained_dof,
//                const std::vector<std::pair<size_type, double> > &entries,
//                const double 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.

class ConstraintMatrix : public Subscriptor
{
public:
  typedef types::global_dof_index size_type;

  enum MergeConflictBehavior
  {
    no_conflicts_allowed,

    left_object_wins,

    right_object_wins
  };

  explicit ConstraintMatrix (const IndexSet &local_constraints = IndexSet());

  explicit ConstraintMatrix (const ConstraintMatrix &constraint_matrix);

  void reinit (const IndexSet &local_constraints = IndexSet());

  bool can_store_line (const size_type line_index) const;

  const IndexSet &get_local_lines() const;

  void add_selected_constraints (const ConstraintMatrix &constraints_in,
                                 const IndexSet         &filter);

  void add_line (const size_type line);

  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,
                  const size_type column,
                  const double value);

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

  void set_inhomogeneity (const size_type line,
                          const double    value);

  void close ();

  void merge (const ConstraintMatrix &other_constraints,
              const MergeConflictBehavior merge_conflict_behavior = no_conflicts_allowed);

  void shift (const size_type offset);

  void clear ();

  size_type n_constraints () const;

  bool is_constrained (const size_type index) const;

  bool is_identity_constrained (const size_type index) const;

  bool are_identity_constrained (const size_type index1,
                                 const size_type index2) 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,double> > *
  get_constraint_entries (const size_type line) const;

  double get_inhomogeneity (const size_type line) const;

  void print (std::ostream &) 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;

  template<typename number>
  void condense (SparseMatrix<number> &matrix) const;

  template <typename number>
  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<typename number, class VectorType>
  void condense (SparseMatrix<number> &matrix,
                 VectorType           &vector) const;

  template <typename number, 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, typename LocalType>
  void
  distribute_local_to_global (const Vector<LocalType>      &local_vector,
                              const std::vector<size_type> &local_dof_indices,
                              VectorType                   &global_vector,
                              const FullMatrix<LocalType>  &local_matrix) const;

  template <class VectorType>
  void
  distribute_local_to_global (const size_type index,
                              const double    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<typename MatrixType::value_type> &local_matrix,
                              const std::vector<size_type> &local_dof_indices,
                              MatrixType                   &global_matrix) const;

  template <typename MatrixType>
  void
  distribute_local_to_global (const FullMatrix<typename MatrixType::value_type> &local_matrix,
                              const std::vector<size_type> &row_indices,
                              const std::vector<size_type> &col_indices,
                              MatrixType                   &global_matrix) const;

  template <typename MatrixType, typename VectorType>
  void
  distribute_local_to_global (const FullMatrix<typename MatrixType::value_type> &local_matrix,
                              const Vector<typename VectorType::value_type>     &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                 = default_empty_table) 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                 = default_empty_table) 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;

  DeclException0 (ExcMatrixIsClosed);
  DeclException0 (ExcMatrixNotClosed);
  DeclException1 (ExcLineInexistant,
                  size_type,
                  << "The specified line " << arg1
                  << " does not exist.");
  DeclException4 (ExcEntryAlreadyExists,
                  size_type, size_type, double, double,
                  << "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 constraint matrix 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 constraint matrix 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 ConstraintMatrix "
                  << "with the appropriate locally_relevant set so "
                  << "that every processor who owns a DoF that constrains "
                  << "another DoF also knows about this constraint?");

private:

  struct ConstraintLine
  {
    typedef std::vector<std::pair<size_type,double> > Entries;

    size_type line;

    Entries entries;

    double inhomogeneity;

    bool operator < (const ConstraintLine &) const;

    bool operator == (const ConstraintLine &) const;

    std::size_t memory_consumption () const;
  };

  std::vector<ConstraintLine> lines;

  std::vector<size_type> lines_cache;

  IndexSet local_lines;

  bool sorted;

  size_type calculate_line_index (const size_type line) const;

  static bool check_zero_weight (const std::pair<size_type, double> &p);

  static const Table<2,bool> default_empty_table;

  template <typename MatrixType, typename VectorType>
  void
  distribute_local_to_global (const FullMatrix<typename MatrixType::value_type>  &local_matrix,
                              const Vector<typename VectorType::value_type>      &local_vector,
                              const std::vector<size_type> &local_dof_indices,
                              MatrixType                   &global_matrix,
                              VectorType                   &global_vector,
                              bool                          use_inhomogeneities_for_rhs,
                              internal::bool2type<false>) const;

  template <typename MatrixType, typename VectorType>
  void
  distribute_local_to_global (const FullMatrix<typename MatrixType::value_type>  &local_matrix,
                              const Vector<typename VectorType::value_type>      &local_vector,
                              const std::vector<size_type> &local_dof_indices,
                              MatrixType                   &global_matrix,
                              VectorType                   &global_vector,
                              bool                          use_inhomogeneities_for_rhs,
                              internal::bool2type<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,
                               internal::bool2type<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,
                               internal::bool2type<true>) const;

  void
  make_sorted_row_list (const std::vector<size_type>   &local_dof_indices,
                        internals::GlobalRowsFromLocal &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 LocalType>
  LocalType
  resolve_vector_entry (const size_type                       i,
                        const internals::GlobalRowsFromLocal &global_rows,
                        const Vector<LocalType>              &local_vector,
                        const std::vector<size_type>         &local_dof_indices,
                        const FullMatrix<LocalType>          &local_matrix) const;

  ConstraintMatrix &operator= (const ConstraintMatrix &other);
};



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

inline
ConstraintMatrix::ConstraintMatrix (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/constraint_matrix_crash_01
  local_lines.compress();
}



inline
ConstraintMatrix::ConstraintMatrix (const ConstraintMatrix &constraint_matrix)
  :
  Subscriptor (),
  lines (constraint_matrix.lines),
  lines_cache (constraint_matrix.lines_cache),
  local_lines (constraint_matrix.local_lines),
  sorted (constraint_matrix.sorted)
{}


inline
void
ConstraintMatrix::add_line (const size_type line)
{
  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 != numbers::invalid_size_type,
          ExcInternalError());
  const size_type line_index = calculate_line_index (line);

  // check whether line already exists; it may, in which case we can just quit
  if (is_constrained(line))
    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.push_back (ConstraintLine());
  lines.back().line = line;
  lines.back().inhomogeneity = 0.;
  lines_cache[line_index] = lines.size()-1;
}



inline
void
ConstraintMatrix::add_entry (const size_type line,
                             const size_type column,
                             const double    value)
{
  Assert (sorted==false, ExcMatrixIsClosed());
  Assert (line != column,
          ExcMessage ("Can't constrain a degree of freedom to itself"));

  // 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[calculate_line_index(line)] != numbers::invalid_size_type,
          ExcInternalError());
  Assert (!local_lines.size() || local_lines.is_element(column),
          ExcColumnNotStoredHere(line, column));
  ConstraintLine *line_ptr = &lines[lines_cache[calculate_line_index(line)]];
  Assert (line_ptr->line == line, ExcInternalError());
  for (ConstraintLine::Entries::const_iterator
       p=line_ptr->entries.begin();
       p != line_ptr->entries.end(); ++p)
    if (p->first == column)
      {
        Assert (std::fabs(p->second - value) < 1.e-14,
                ExcEntryAlreadyExists(line, column, p->second, value));
        return;
      }

  line_ptr->entries.push_back (std::make_pair(column,value));
}



inline
void
ConstraintMatrix::set_inhomogeneity (const size_type line,
                                     const double    value)
{
  const size_type line_index = calculate_line_index(line);
  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;
}



inline
types::global_dof_index
ConstraintMatrix::n_constraints () const
{
  return lines.size();
}



inline
bool
ConstraintMatrix::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));
}



inline
bool
ConstraintMatrix::is_inhomogeneously_constrained (const size_type index) 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(index);
  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 == 0);
    }
}



inline
const std::vector<std::pair<types::global_dof_index,double> > *
ConstraintMatrix::get_constraint_entries (const size_type line) 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);
  if (line_index >= lines_cache.size() ||
      lines_cache[line_index] == numbers::invalid_size_type)
    return 0;
  else
    return &lines[lines_cache[line_index]].entries;
}



inline
double
ConstraintMatrix::get_inhomogeneity (const size_type line) 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);
  if (line_index >= lines_cache.size() ||
      lines_cache[line_index] == numbers::invalid_size_type)
    return 0;
  else
    return lines[lines_cache[line_index]].inhomogeneity;
}



inline types::global_dof_index
ConstraintMatrix::calculate_line_index (const size_type line) const
{
  //IndexSet is unused (serial case)
  if (!local_lines.size())
    return line;

  Assert(local_lines.is_element(line),
         ExcRowNotStoredHere(line));

  return local_lines.index_within_set(line);
}



inline bool
ConstraintMatrix::can_store_line (size_type line_index) const
{
  return !local_lines.size() || local_lines.is_element(line_index);
}



inline
const IndexSet &
ConstraintMatrix::get_local_lines () const
{
  return local_lines;
}



template <class VectorType>
inline
void ConstraintMatrix::distribute_local_to_global (
  const size_type index,
  const double    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 ForwardIteratorVec, typename ForwardIteratorInd,
          class VectorType>
inline
void ConstraintMatrix::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)
        global_vector(*local_indices_begin) += *local_vector_begin;
      else
        {
          const ConstraintLine &position =
            lines[lines_cache[calculate_line_index(*local_indices_begin)]];
          for (size_type j=0; j<position.entries.size(); ++j)
            global_vector(position.entries[j].first)
            += *local_vector_begin * position.entries[j].second;
        }
    }
}


template <class InVector, class OutVector>
inline
void
ConstraintMatrix::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 ForwardIteratorVec, typename ForwardIteratorInd,
          class VectorType>
inline
void ConstraintMatrix::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
                                    ((MatrixType *)0))
                             ==
                             sizeof(yes_type));
};


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


template <typename MatrixType>
inline
void
ConstraintMatrix::
distribute_local_to_global (const FullMatrix<typename MatrixType::value_type>     &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,
                              ::internal::bool2type<IsBlockMatrix<MatrixType>::value>());
}




template <typename MatrixType, typename VectorType>
inline
void
ConstraintMatrix::
distribute_local_to_global (const FullMatrix<typename MatrixType::value_type>     &local_matrix,
                            const Vector<typename VectorType::value_type>         &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,
                              ::internal::bool2type<IsBlockMatrix<MatrixType>::value>());
}




template <typename SparsityPatternType>
inline
void
ConstraintMatrix::
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,
                               internal::bool2type<IsBlockMatrix<SparsityPatternType>::value>());
}


DEAL_II_NAMESPACE_CLOSE

#endif