Reference documentation for deal.II version GIT 85919f3e70 2023-05-28 07:10:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
affine_constraints.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_affine_constraints_h
17 #define dealii_affine_constraints_h
18 
19 #include <deal.II/base/config.h>
20 
22 #include <deal.II/base/index_set.h>
24 #include <deal.II/base/table.h>
27 
29 #include <deal.II/lac/vector.h>
31 
32 #include <boost/range/iterator_range.hpp>
33 
34 #include <algorithm>
35 #include <set>
36 #include <type_traits>
37 #include <utility>
38 #include <vector>
39 
41 
42 // Forward declarations
43 #ifndef DOXYGEN
44 template <typename>
45 class FullMatrix;
46 class SparsityPattern;
50 template <typename number>
51 class SparseMatrix;
52 template <typename number>
53 class BlockSparseMatrix;
54 
55 namespace internal
56 {
57  namespace AffineConstraints
58  {
60 
78  struct Distributing
79  {
80  Distributing(const size_type global_row = numbers::invalid_size_type,
81  const size_type local_row = numbers::invalid_size_type);
82 
83  Distributing(const Distributing &in);
84 
85  Distributing &
86  operator=(const Distributing &in);
87 
88  bool
89  operator<(const Distributing &in) const
90  {
91  return global_row < in.global_row;
92  }
93 
94  size_type global_row;
95  size_type local_row;
96  mutable size_type constraint_position;
97  };
98 
99 
100 
112  template <typename number>
113  struct DataCache
114  {
115  DataCache();
116 
117  void
118  reinit();
119 
120  size_type
121  insert_new_index(const std::pair<size_type, number> &pair);
122 
123  void
124  append_index(const size_type index,
125  const std::pair<size_type, number> &pair);
126 
127  size_type
128  get_size(const size_type index) const;
129 
130  const std::pair<size_type, number> *
131  get_entry(const size_type index) const;
132 
133  size_type row_length;
134 
135  std::vector<std::pair<size_type, number>> data;
136 
137  std::vector<size_type> individual_size;
138  };
139 
140 
141 
170  template <typename number>
171  class GlobalRowsFromLocal
172  {
173  public:
177  GlobalRowsFromLocal();
178 
179  void
180  reinit(const size_type n_local_rows);
181 
182  void
183  insert_index(const size_type global_row,
184  const size_type local_row,
185  const number constraint_value);
186  void
187  sort();
188 
189  void
190  print(std::ostream &os);
191 
195  size_type
196  size() const;
197 
202  size_type
203  size(const size_type counter_index) const;
204 
208  size_type
209  global_row(const size_type counter_index) const;
210 
214  size_type &
215  global_row(const size_type counter_index);
216 
222  size_type
223  local_row(const size_type counter_index) const;
224 
228  size_type &
229  local_row(const size_type counter_index);
230 
236  size_type
237  local_row(const size_type counter_index,
238  const size_type index_in_constraint) const;
239 
244  number
245  constraint_value(const size_type counter_index,
246  const size_type index_in_constraint) const;
247 
253  bool
254  have_indirect_rows() const;
255 
260  void
261  insert_constraint(const size_type constrained_local_dof);
262 
268  size_type
269  n_constraints() const;
270 
275  size_type
276  n_inhomogeneities() const;
277 
284  void
285  set_ith_constraint_inhomogeneous(const size_type i);
286 
291  size_type
292  constraint_origin(size_type i) const;
293 
299  std::vector<Distributing> total_row_indices;
300 
301  private:
305  DataCache<number> data_cache;
306 
311  size_type n_active_rows;
312 
317  size_type n_inhomogeneous_rows;
318  };
319 
320 
321 
339  template <typename number>
340  struct ScratchData
341  {
345  ScratchData()
346  : in_use(false)
347  {}
348 
352  ScratchData(const ScratchData &)
353  : in_use(false)
354  {}
355 
359  bool in_use;
360 
364  std::vector<std::pair<size_type, size_type>> new_entries;
365 
369  std::vector<size_type> rows;
370 
374  std::vector<size_type> columns;
375 
379  std::vector<number> values;
380 
384  std::vector<size_type> block_starts;
385 
389  std::vector<size_type> vector_indices;
390 
394  std::vector<number> vector_values;
395 
399  GlobalRowsFromLocal<number> global_rows;
400 
404  GlobalRowsFromLocal<number> global_columns;
405  };
406  } // namespace AffineConstraints
407 } // namespace internal
408 
409 namespace internal
410 {
411  namespace AffineConstraintsImplementation
412  {
413  template <class VectorType>
414  void
415  set_zero_all(const std::vector<types::global_dof_index> &cm,
416  VectorType & vec);
417 
418  template <class T>
419  void
420  set_zero_all(const std::vector<types::global_dof_index> &cm,
421  ::Vector<T> & vec);
422 
423  template <class T>
424  void
425  set_zero_all(const std::vector<types::global_dof_index> &cm,
426  ::BlockVector<T> & vec);
427  } // namespace AffineConstraintsImplementation
428 } // namespace internal
429 #endif
430 
431 // TODO[WB]: We should have a function of the kind
432 // AffineConstraints::add_constraint (const size_type constrained_dof,
433 // const std::vector<std::pair<size_type, number> > &entries,
434 // const number inhomogeneity = 0);
435 // rather than building up constraints piecemeal through add_line/add_entry
436 // etc. This would also eliminate the possibility of accidentally changing
437 // existing constraints into something pointless, see the discussion on the
438 // mailing list on "Tiny bug in interpolate_boundary_values" in Sept. 2010.
439 
513 template <typename number = double>
515 {
516 public:
521 
528  {
534 
541 
548  };
549 
567  explicit AffineConstraints(const IndexSet &local_constraints = IndexSet());
568 
572  explicit AffineConstraints(const AffineConstraints &affine_constraints);
573 
577  AffineConstraints(AffineConstraints &&affine_constraints) noexcept =
578  default; // NOLINT
579 
590  operator=(const AffineConstraints &) = delete;
591 
596  operator=(AffineConstraints &&affine_constraints) noexcept =
597  default; // NOLINT
598 
605  template <typename other_number>
606  void
608 
615  void
616  reinit(const IndexSet &local_constraints = IndexSet());
617 
624  bool
625  can_store_line(const size_type line_n) const;
626 
633  const IndexSet &
635 
655  void
657  const IndexSet & filter);
658 
668  void
669  add_line(const size_type line_n);
670 
683  void
684  add_lines(const std::vector<bool> &lines);
685 
698  void
699  add_lines(const std::set<size_type> &lines);
700 
713  void
715 
731  void
732  add_entry(const size_type constrained_dof_index,
733  const size_type column,
734  const number weight);
735 
741  void
743  const size_type constrained_dof_index,
744  const std::vector<std::pair<size_type, number>> &col_weight_pairs);
745 
757  void
758  set_inhomogeneity(const size_type constrained_dof_index, const number value);
759 
781  void
782  close();
783 
789  bool
790  is_closed() const;
791 
797  bool
798  is_closed(const MPI_Comm comm) const;
799 
824  void
826  const AffineConstraints & other_constraints,
827  const MergeConflictBehavior merge_conflict_behavior = no_conflicts_allowed,
828  const bool allow_different_local_lines = false);
829 
842  void
843  shift(const size_type offset);
844 
852  void
853  clear();
854 
867  size_type
868  n_constraints() const;
869 
874  size_type
875  n_identities() const;
876 
882  size_type
884 
894  bool
895  is_constrained(const size_type line_n) const;
896 
908  bool
909  is_identity_constrained(const size_type line_n) const;
910 
917  bool
919  const size_type line_n_2) const;
920 
931  size_type
933 
938  bool
940 
946  bool
948 
953  const std::vector<std::pair<size_type, number>> *
954  get_constraint_entries(const size_type line_n) const;
955 
960  number
961  get_inhomogeneity(const size_type line_n) const;
962 
983  void
984  print(std::ostream &out) const;
985 
998  void
999  write_dot(std::ostream &) const;
1000 
1005  std::size_t
1007 
1014  void
1015  resolve_indices(std::vector<types::global_dof_index> &indices) const;
1016 
1042  void
1043  condense(SparsityPattern &sparsity) const;
1044 
1048  void
1049  condense(BlockSparsityPattern &sparsity) const;
1050 
1055  void
1057 
1062  void
1064 
1072  void
1074 
1078  void
1080 
1092  template <class VectorType>
1093  void
1094  condense(VectorType &vec) const;
1095 
1102  template <class VectorType>
1103  void
1104  condense(const VectorType &vec_ghosted, VectorType &output) const;
1105 
1118  template <class VectorType>
1119  void
1120  condense(SparseMatrix<number> &matrix, VectorType &vector) const;
1121 
1126  template <class BlockVectorType>
1127  void
1128  condense(BlockSparseMatrix<number> &matrix, BlockVectorType &vector) const;
1129 
1136  template <class VectorType>
1137  void
1138  set_zero(VectorType &vec) const;
1139 
1195  template <class InVector, class OutVector>
1196  void
1197  distribute_local_to_global(const InVector & local_vector,
1198  const std::vector<size_type> &local_dof_indices,
1199  OutVector & global_vector) const;
1200 
1248  template <typename VectorType>
1249  void
1251  const std::vector<size_type> &local_dof_indices,
1252  VectorType & global_vector,
1253  const FullMatrix<number> & local_matrix) const;
1254 
1274  template <typename VectorType>
1275  void
1277  const Vector<number> & local_vector,
1278  const std::vector<size_type> &local_dof_indices_row,
1279  const std::vector<size_type> &local_dof_indices_col,
1280  VectorType & global_vector,
1281  const FullMatrix<number> & local_matrix,
1282  bool diagonal = false) const;
1283 
1287  template <class VectorType>
1288  void
1290  const number value,
1291  VectorType & global_vector) const;
1292 
1325  template <typename ForwardIteratorVec,
1326  typename ForwardIteratorInd,
1327  class VectorType>
1328  void
1329  distribute_local_to_global(ForwardIteratorVec local_vector_begin,
1330  ForwardIteratorVec local_vector_end,
1331  ForwardIteratorInd local_indices_begin,
1332  VectorType & global_vector) const;
1333 
1385  template <typename MatrixType>
1386  void
1388  const std::vector<size_type> &local_dof_indices,
1389  MatrixType & global_matrix) const;
1390 
1418  template <typename MatrixType>
1419  void
1421  const std::vector<size_type> &row_indices,
1422  const std::vector<size_type> &col_indices,
1423  MatrixType & global_matrix) const;
1424 
1441  template <typename MatrixType>
1442  void
1444  const std::vector<size_type> &row_indices,
1445  const AffineConstraints &column_affine_constraints,
1446  const std::vector<size_type> &column_indices,
1447  MatrixType & global_matrix) const;
1448 
1469  template <typename MatrixType, typename VectorType>
1470  void
1472  const Vector<number> & local_vector,
1473  const std::vector<size_type> &local_dof_indices,
1474  MatrixType & global_matrix,
1475  VectorType & global_vector,
1476  bool use_inhomogeneities_for_rhs = false) const;
1477 
1531  void
1533  const std::vector<size_type> &local_dof_indices,
1534  SparsityPatternBase & sparsity_pattern,
1535  const bool keep_constrained_entries = true,
1536  const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1537 
1541  void
1543  const std::vector<size_type> &row_indices,
1544  const std::vector<size_type> &col_indices,
1545  SparsityPatternBase & sparsity_pattern,
1546  const bool keep_constrained_entries = true,
1547  const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1548 
1553  void
1555  const std::vector<size_type> & row_indices,
1556  const AffineConstraints<number> &col_constraints,
1557  const std::vector<size_type> & col_indices,
1558  SparsityPatternBase & sparsity_pattern,
1559  const bool keep_constrained_entries = true,
1560  const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1561 
1581  template <typename ForwardIteratorVec,
1582  typename ForwardIteratorInd,
1583  class VectorType>
1584  void
1585  get_dof_values(const VectorType & global_vector,
1586  ForwardIteratorInd local_indices_begin,
1587  ForwardIteratorVec local_vector_begin,
1588  ForwardIteratorVec local_vector_end) const;
1589 
1611  template <class VectorType>
1612  void
1613  distribute(VectorType &vec) const;
1614 
1623  {
1628  using Entries = std::vector<std::pair<size_type, number>>;
1629 
1636 
1645 
1650 
1655  const Entries & entries = {},
1656  const number & inhomogeneity = 0.0);
1657 
1661  template <typename ConstraintLineType>
1662  ConstraintLine(const ConstraintLineType &other);
1663 
1667  template <typename ConstraintLineType>
1668  ConstraintLine &
1669  operator=(const ConstraintLineType &other);
1670 
1678  bool
1679  operator<(const ConstraintLine &) const;
1680 
1686  bool
1687  operator==(const ConstraintLine &) const;
1688 
1693  std::size_t
1695 
1701  template <class Archive>
1702  void
1703  serialize(Archive &ar, const unsigned int)
1704  {
1705  ar &index &entries &inhomogeneity;
1706  }
1707  };
1708 
1712  using const_iterator = typename std::vector<ConstraintLine>::const_iterator;
1713 
1717  using LineRange = boost::iterator_range<const_iterator>;
1718 
1727  const LineRange
1728  get_lines() const;
1729 
1758  bool
1759  is_consistent_in_parallel(const std::vector<IndexSet> &locally_owned_dofs,
1760  const IndexSet & locally_active_dofs,
1761  const MPI_Comm mpi_communicator,
1762  const bool verbose = false) const;
1763 
1769  void
1770  make_consistent_in_parallel(const IndexSet &locally_owned_dofs,
1771  const IndexSet &locally_relevant_dofs,
1772  const MPI_Comm mpi_communicator);
1773 
1792  size_type,
1793  << "The specified line " << arg1 << " does not exist.");
1800  size_type,
1801  size_type,
1802  number,
1803  number,
1804  << "The entry for the indices " << arg1 << " and " << arg2
1805  << " already exists, but the values " << arg3 << " (old) and "
1806  << arg4 << " (new) differ "
1807  << "by " << (arg4 - arg3) << '.');
1814  int,
1815  int,
1816  << "You tried to constrain DoF " << arg1 << " to DoF " << arg2
1817  << ", but that one is also constrained. This is not allowed!");
1824  size_type,
1825  << "Degree of freedom " << arg1
1826  << " is constrained from both object in a merge operation.");
1833  size_type,
1834  << "In the given argument a degree of freedom is constrained "
1835  << "to another DoF with number " << arg1
1836  << ", which however is constrained by this object. This is not"
1837  << " allowed.");
1844  size_type,
1845  << "The index set given to this constraints object indicates "
1846  << "constraints for degree of freedom " << arg1
1847  << " should not be stored by this object, but a constraint "
1848  << "is being added.");
1849 
1856  size_type,
1857  size_type,
1858  << "The index set given to this constraints object indicates "
1859  << "constraints using degree of freedom " << arg2
1860  << " should not be stored by this object, but a constraint "
1861  << "for degree of freedom " << arg1 << " uses it.");
1862 
1869  int,
1870  int,
1871  << "While distributing the constraint for DoF " << arg1
1872  << ", it turns out that one of the processors "
1873  << "who own the " << arg2 << " degrees of freedom that x_"
1874  << arg1 << " is constrained against does not know about "
1875  << "the constraint on x_" << arg1
1876  << ". Did you not initialize the AffineConstraints container "
1877  << "with the appropriate locally_relevant set so "
1878  << "that every processor who owns a DoF that constrains "
1879  << "another DoF also knows about this constraint?");
1880 
1881  template <typename>
1882  friend class AffineConstraints;
1883 
1884 private:
1896  std::vector<ConstraintLine> lines;
1897 
1930  std::vector<size_type> lines_cache;
1931 
1938 
1942  bool sorted;
1943 
1945  internal::AffineConstraints::ScratchData<number>>
1947 
1952  size_type
1953  calculate_line_index(const size_type line_n) const;
1954 
1959  template <typename MatrixType, typename VectorType>
1960  void
1962  const Vector<number> & local_vector,
1963  const std::vector<size_type> &local_dof_indices,
1964  MatrixType & global_matrix,
1965  VectorType & global_vector,
1966  const bool use_inhomogeneities_for_rhs,
1967  const std::integral_constant<bool, false>) const;
1968 
1973  template <typename MatrixType, typename VectorType>
1974  void
1976  const Vector<number> & local_vector,
1977  const std::vector<size_type> &local_dof_indices,
1978  MatrixType & global_matrix,
1979  VectorType & global_vector,
1980  const bool use_inhomogeneities_for_rhs,
1981  const std::integral_constant<bool, true>) const;
1982 
1990  void
1991  make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
1992  internal::AffineConstraints::GlobalRowsFromLocal<number>
1993  &global_rows) const;
1994 
2002  void
2003  make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
2004  std::vector<size_type> & active_dofs) const;
2005 
2009  template <typename MatrixScalar, typename VectorScalar>
2012  const size_type i,
2013  const internal::AffineConstraints::GlobalRowsFromLocal<number> &global_rows,
2014  const Vector<VectorScalar> & local_vector,
2015  const std::vector<size_type> & local_dof_indices,
2016  const FullMatrix<MatrixScalar> &local_matrix) const;
2017 };
2018 
2019 /* ---------------- template and inline functions ----------------- */
2020 
2021 template <typename number>
2023  const IndexSet &local_constraints)
2024  : lines()
2025  , local_lines(local_constraints)
2026  , sorted(false)
2027 {
2028  // make sure the IndexSet is compressed. Otherwise this can lead to crashes
2029  // that are hard to find (only happen in release mode).
2030  // see tests/mpi/affine_constraints_crash_01
2032 }
2033 
2034 template <typename number>
2036  const AffineConstraints &affine_constraints)
2037  : Subscriptor()
2038  , lines(affine_constraints.lines)
2039  , lines_cache(affine_constraints.lines_cache)
2040  , local_lines(affine_constraints.local_lines)
2041  , sorted(affine_constraints.sorted)
2042 {}
2043 
2044 template <typename number>
2045 inline void
2047 {
2048  Assert(sorted == false, ExcMatrixIsClosed());
2049 
2050  // the following can happen when we compute with distributed meshes and dof
2051  // handlers and we constrain a degree of freedom whose number we don't have
2052  // locally. if we don't abort here the program will try to allocate several
2053  // terabytes of memory to resize the various arrays below :-)
2055  const size_type line_index = calculate_line_index(line_n);
2056 
2057  // check whether line already exists; it may, in which case we can just quit
2058  if (is_constrained(line_n))
2059  return;
2060 
2061  // if necessary enlarge vector of existing entries for cache
2062  if (line_index >= lines_cache.size())
2063  lines_cache.resize(std::max(2 * static_cast<size_type>(lines_cache.size()),
2064  line_index + 1),
2066 
2067  // push a new line to the end of the list
2068  lines.emplace_back();
2069  lines.back().index = line_n;
2070  lines.back().inhomogeneity = 0.;
2071  lines_cache[line_index] = lines.size() - 1;
2072 }
2073 
2074 
2075 
2076 template <typename number>
2077 inline void
2079  const size_type column,
2080  const number weight)
2081 {
2082  Assert(sorted == false, ExcMatrixIsClosed());
2083  Assert(constrained_dof_index != column,
2084  ExcMessage("Can't constrain a degree of freedom to itself"));
2085 
2086  // Ensure that the current line is present in the cache:
2087  const size_type line_index = calculate_line_index(constrained_dof_index);
2088  Assert(line_index < lines_cache.size(),
2089  ExcMessage("The current AffineConstraints does not contain the line "
2090  "for the current entry. Call AffineConstraints::add_line "
2091  "before calling this function."));
2092 
2093  // if in debug mode, check whether an entry for this column already exists
2094  // and if it's the same as the one entered at present
2095  //
2096  // in any case: exit the function if an entry for this column already
2097  // exists, since we don't want to enter it twice
2099  ExcInternalError());
2101  ExcColumnNotStoredHere(constrained_dof_index, column));
2102  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
2103  Assert(line_ptr->index == constrained_dof_index, ExcInternalError());
2104  for (const auto &p : line_ptr->entries)
2105  if (p.first == column)
2106  {
2107  Assert(std::abs(p.second - weight) < 1.e-14,
2109  constrained_dof_index, column, p.second, weight));
2110  return;
2111  }
2112 
2113  line_ptr->entries.emplace_back(column, weight);
2114 }
2115 
2116 
2117 
2118 template <typename number>
2119 inline void
2121  const size_type constrained_dof_index,
2122  const number value)
2123 {
2124  const size_type line_index = calculate_line_index(constrained_dof_index);
2125  Assert(line_index < lines_cache.size() &&
2126  lines_cache[line_index] != numbers::invalid_size_type,
2127  ExcMessage("call add_line() before calling set_inhomogeneity()"));
2128  Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
2129  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
2130  line_ptr->inhomogeneity = value;
2131 }
2132 
2133 
2134 
2135 template <typename number>
2136 template <class VectorType>
2137 inline void
2139 {
2140  // since lines is a private member, we cannot pass it to the functions
2141  // above. therefore, copy the content which is cheap
2142  std::vector<size_type> constrained_lines(lines.size());
2143  for (unsigned int i = 0; i < lines.size(); ++i)
2144  constrained_lines[i] = lines[i].index;
2145  internal::AffineConstraintsImplementation::set_zero_all(constrained_lines,
2146  vec);
2147 }
2148 
2149 template <typename number>
2152 {
2153  return lines.size();
2154 }
2155 
2156 template <typename number>
2159 {
2160  return std::count_if(lines.begin(),
2161  lines.end(),
2162  [](const ConstraintLine &line) {
2163  return (line.entries.size() == 1) &&
2164  (line.entries[0].second == number(1.));
2165  });
2166 }
2167 
2168 template <typename number>
2171 {
2172  return std::count_if(lines.begin(),
2173  lines.end(),
2174  [](const ConstraintLine &line) {
2175  return (line.inhomogeneity != number(0.));
2176  });
2177 }
2178 
2179 template <typename number>
2180 inline bool
2182 {
2183  if (lines.empty())
2184  return false;
2185 
2186  const size_type line_index = calculate_line_index(index);
2187  return ((line_index < lines_cache.size()) &&
2188  (lines_cache[line_index] != numbers::invalid_size_type));
2189 }
2190 
2191 template <typename number>
2192 inline bool
2194  const size_type line_n) const
2195 {
2196  // check whether the entry is constrained. could use is_constrained, but
2197  // that means computing the line index twice
2198  const size_type line_index = calculate_line_index(line_n);
2199  if (line_index >= lines_cache.size() ||
2200  lines_cache[line_index] == numbers::invalid_size_type)
2201  return false;
2202  else
2203  {
2204  Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
2205  return (lines[lines_cache[line_index]].inhomogeneity != number(0.));
2206  }
2207 }
2208 
2209 template <typename number>
2210 inline const std::vector<std::pair<types::global_dof_index, number>> *
2212 {
2213  if (lines.empty())
2214  return nullptr;
2215 
2216  // check whether the entry is constrained. could use is_constrained, but
2217  // that means computing the line index twice
2218  const size_type line_index = calculate_line_index(line_n);
2219  if (line_index >= lines_cache.size() ||
2220  lines_cache[line_index] == numbers::invalid_size_type)
2221  return nullptr;
2222  else
2223  return &lines[lines_cache[line_index]].entries;
2224 }
2225 
2226 template <typename number>
2227 inline number
2229 {
2230  // check whether the entry is constrained. could use is_constrained, but
2231  // that means computing the line index twice
2232  const size_type line_index = calculate_line_index(line_n);
2233  if (line_index >= lines_cache.size() ||
2234  lines_cache[line_index] == numbers::invalid_size_type)
2235  return 0;
2236  else
2237  return lines[lines_cache[line_index]].inhomogeneity;
2238 }
2239 
2240 template <typename number>
2243 {
2244  // IndexSet is unused (serial case)
2245  if (local_lines.size() == 0)
2246  return line_n;
2247 
2249 
2250  return local_lines.index_within_set(line_n);
2251 }
2252 
2253 template <typename number>
2254 inline bool
2256 {
2257  return local_lines.size() == 0 || local_lines.is_element(line_n);
2258 }
2259 
2260 template <typename number>
2261 inline const IndexSet &
2263 {
2264  return local_lines;
2265 }
2266 
2267 template <typename number>
2268 template <class VectorType>
2269 inline void
2271  const size_type index,
2272  const number value,
2273  VectorType & global_vector) const
2274 {
2275  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2276 
2277  if (is_constrained(index) == false)
2278  global_vector(index) += value;
2279  else
2280  {
2281  const ConstraintLine &position =
2283  for (size_type j = 0; j < position.entries.size(); ++j)
2284  global_vector(position.entries[j].first) +=
2285  value * position.entries[j].second;
2286  }
2287 }
2288 
2289 template <typename number>
2290 template <typename ForwardIteratorVec,
2291  typename ForwardIteratorInd,
2292  class VectorType>
2293 inline void
2295  ForwardIteratorVec local_vector_begin,
2296  ForwardIteratorVec local_vector_end,
2297  ForwardIteratorInd local_indices_begin,
2298  VectorType & global_vector) const
2299 {
2300  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2301  for (; local_vector_begin != local_vector_end;
2302  ++local_vector_begin, ++local_indices_begin)
2303  {
2304  if (is_constrained(*local_indices_begin) == false)
2305  internal::ElementAccess<VectorType>::add(*local_vector_begin,
2306  *local_indices_begin,
2307  global_vector);
2308  else
2309  {
2310  const ConstraintLine &position =
2311  lines[lines_cache[calculate_line_index(*local_indices_begin)]];
2312  for (size_type j = 0; j < position.entries.size(); ++j)
2314  (*local_vector_begin) * position.entries[j].second,
2315  position.entries[j].first,
2316  global_vector);
2317  }
2318  }
2319 }
2320 
2321 template <typename number>
2322 template <class InVector, class OutVector>
2323 inline void
2325  const InVector & local_vector,
2326  const std::vector<size_type> &local_dof_indices,
2327  OutVector & global_vector) const
2328 {
2329  Assert(local_vector.size() == local_dof_indices.size(),
2330  ExcDimensionMismatch(local_vector.size(), local_dof_indices.size()));
2331  distribute_local_to_global(local_vector.begin(),
2332  local_vector.end(),
2333  local_dof_indices.begin(),
2334  global_vector);
2335 }
2336 
2337 template <typename number>
2338 template <typename ForwardIteratorVec,
2339  typename ForwardIteratorInd,
2340  class VectorType>
2341 inline void
2343  const VectorType & global_vector,
2344  ForwardIteratorInd local_indices_begin,
2345  ForwardIteratorVec local_vector_begin,
2346  ForwardIteratorVec local_vector_end) const
2347 {
2348  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2349  for (; local_vector_begin != local_vector_end;
2350  ++local_vector_begin, ++local_indices_begin)
2351  {
2352  if (is_constrained(*local_indices_begin) == false)
2353  *local_vector_begin = global_vector(*local_indices_begin);
2354  else
2355  {
2356  const ConstraintLine &position =
2357  lines[lines_cache[calculate_line_index(*local_indices_begin)]];
2358  typename VectorType::value_type value = position.inhomogeneity;
2359  for (size_type j = 0; j < position.entries.size(); ++j)
2360  value += (global_vector(position.entries[j].first) *
2361  position.entries[j].second);
2362  *local_vector_begin = value;
2363  }
2364  }
2365 }
2366 
2367 // Forward declarations
2368 #ifndef DOXYGEN
2369 template <typename MatrixType>
2370 class BlockMatrixBase;
2371 template <typename SparsityPatternType>
2373 template <typename number>
2374 class BlockSparseMatrixEZ;
2375 
2376 namespace internal
2377 {
2378  namespace AffineConstraints
2379  {
2397  template <typename MatrixType>
2398  struct IsBlockMatrix
2399  {
2400  private:
2406  template <typename T>
2407  static std::true_type
2408  check(const BlockMatrixBase<T> *);
2409 
2414  template <typename T>
2415  static std::true_type
2416  check(const BlockSparseMatrixEZ<T> *);
2417 
2422  static std::false_type
2423  check(...);
2424 
2425  public:
2432  static const bool value =
2433  std::is_same<decltype(check(std::declval<MatrixType *>())),
2434  std::true_type>::value;
2435  };
2436 
2437  // instantiation of the static member
2438  template <typename MatrixType>
2439  const bool IsBlockMatrix<MatrixType>::value;
2440 
2441  } // namespace AffineConstraints
2442 } // namespace internal
2443 #endif
2444 
2445 
2446 
2447 template <typename number>
2448 template <typename other_number>
2449 inline void
2451  const AffineConstraints<other_number> &other)
2452 {
2453  lines.clear();
2454  lines.insert(lines.begin(), other.lines.begin(), other.lines.end());
2455  lines_cache = other.lines_cache;
2456  local_lines = other.local_lines;
2457  sorted = other.sorted;
2458 }
2459 
2460 
2461 
2462 template <typename number>
2463 template <typename MatrixType>
2464 inline void
2466  const FullMatrix<number> & local_matrix,
2467  const std::vector<size_type> &local_dof_indices,
2468  MatrixType & global_matrix) const
2469 {
2470  // create a dummy and hand on to the function actually implementing this
2471  // feature in the cm.templates.h file.
2473  distribute_local_to_global(
2474  local_matrix,
2475  dummy,
2476  local_dof_indices,
2477  global_matrix,
2478  dummy,
2479  false,
2480  std::integral_constant<
2481  bool,
2482  internal::AffineConstraints::IsBlockMatrix<MatrixType>::value>());
2483 }
2484 
2485 
2486 
2487 template <typename number>
2488 template <typename MatrixType, typename VectorType>
2489 inline void
2491  const FullMatrix<number> & local_matrix,
2492  const Vector<number> & local_vector,
2493  const std::vector<size_type> &local_dof_indices,
2494  MatrixType & global_matrix,
2495  VectorType & global_vector,
2496  bool use_inhomogeneities_for_rhs) const
2497 {
2498  // enter the internal function with the respective block information set,
2499  // the actual implementation follows in the cm.templates.h file.
2500  distribute_local_to_global(
2501  local_matrix,
2502  local_vector,
2503  local_dof_indices,
2504  global_matrix,
2505  global_vector,
2506  use_inhomogeneities_for_rhs,
2507  std::integral_constant<
2508  bool,
2509  internal::AffineConstraints::IsBlockMatrix<MatrixType>::value>());
2510 }
2511 
2512 
2513 
2514 template <typename number>
2516  const size_type & index,
2518  const number &inhomogeneity)
2519  : index(index)
2520  , entries(entries)
2521  , inhomogeneity(inhomogeneity)
2522 {}
2523 
2524 
2525 
2526 template <typename number>
2527 template <typename ConstraintLineType>
2529  const ConstraintLineType &other)
2530 {
2531  this->index = other.index;
2532 
2533  entries.clear();
2534  entries.insert(entries.begin(), other.entries.begin(), other.entries.end());
2535 
2536  this->inhomogeneity = other.inhomogeneity;
2537 }
2538 
2539 
2540 
2541 template <typename number>
2542 template <typename ConstraintLineType>
2545  const ConstraintLineType &other)
2546 {
2547  this->index = other.index;
2548 
2549  entries.clear();
2550  entries.insert(entries.begin(), other.entries.begin(), other.entries.end());
2551 
2552  this->inhomogeneity = other.inhomogeneity;
2553 
2554  return *this;
2555 }
2556 
2558 
2559 #endif
const LineRange get_lines() const
void condense(DynamicSparsityPattern &sparsity) const
void distribute_local_to_global(ForwardIteratorVec local_vector_begin, ForwardIteratorVec local_vector_end, ForwardIteratorInd local_indices_begin, VectorType &global_vector) const
bool is_closed() const
void merge(const AffineConstraints &other_constraints, const MergeConflictBehavior merge_conflict_behavior=no_conflicts_allowed, const bool allow_different_local_lines=false)
size_type n_inhomogeneities() const
types::global_dof_index size_type
bool has_inhomogeneities() const
void add_line(const size_type line_n)
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
void distribute_local_to_global(const FullMatrix< number > &local_matrix, const Vector< number > &local_vector, const std::vector< size_type > &local_dof_indices, MatrixType &global_matrix, VectorType &global_vector, const bool use_inhomogeneities_for_rhs, const std::integral_constant< bool, true >) const
void add_selected_constraints(const AffineConstraints &constraints_in, const IndexSet &filter)
void condense(VectorType &vec) const
bool can_store_line(const size_type line_n) const
AffineConstraints(AffineConstraints &&affine_constraints) noexcept=default
Threads::ThreadLocalStorage< internal::AffineConstraints::ScratchData< number > > scratch_data
friend class AffineConstraints
void add_entry(const size_type constrained_dof_index, const size_type column, const number weight)
void reinit(const IndexSet &local_constraints=IndexSet())
void make_consistent_in_parallel(const IndexSet &locally_owned_dofs, const IndexSet &locally_relevant_dofs, const MPI_Comm mpi_communicator)
number get_inhomogeneity(const size_type line_n) const
void make_sorted_row_list(const std::vector< size_type > &local_dof_indices, std::vector< size_type > &active_dofs) const
void get_dof_values(const VectorType &global_vector, ForwardIteratorInd local_indices_begin, ForwardIteratorVec local_vector_begin, ForwardIteratorVec local_vector_end) const
void distribute_local_to_global(const InVector &local_vector, const std::vector< size_type > &local_dof_indices, OutVector &global_vector) const
void make_sorted_row_list(const std::vector< size_type > &local_dof_indices, internal::AffineConstraints::GlobalRowsFromLocal< number > &global_rows) const
const IndexSet & get_local_lines() 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 shift(const size_type offset)
size_type n_identities() const
std::size_t memory_consumption() const
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
typename std::vector< ConstraintLine >::const_iterator const_iterator
void condense(SparseMatrix< number > &matrix, VectorType &vector) const
void set_inhomogeneity(const size_type constrained_dof_index, const number value)
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
bool is_closed(const MPI_Comm comm) const
void condense(BlockSparsityPattern &sparsity) const
void distribute_local_to_global(const FullMatrix< number > &local_matrix, const Vector< number > &local_vector, const std::vector< size_type > &local_dof_indices, MatrixType &global_matrix, VectorType &global_vector, const bool use_inhomogeneities_for_rhs, const std::integral_constant< bool, false >) const
void condense(SparsityPattern &sparsity) const
void distribute_local_to_global(const size_type index, const number value, VectorType &global_vector) const
size_type max_constraint_indirections() const
void add_lines(const std::set< size_type > &lines)
AffineConstraints(const AffineConstraints &affine_constraints)
void distribute(VectorType &vec) const
void add_lines(const IndexSet &lines)
void condense(BlockSparseMatrix< number > &matrix) 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 resolve_indices(std::vector< types::global_dof_index > &indices) const
std::vector< ConstraintLine > lines
bool is_constrained(const size_type line_n) const
void condense(BlockDynamicSparsityPattern &sparsity) const
void distribute_local_to_global(const FullMatrix< number > &local_matrix, const std::vector< size_type > &local_dof_indices, MatrixType &global_matrix) const
AffineConstraints & operator=(const AffineConstraints &)=delete
const std::vector< std::pair< size_type, number > > * get_constraint_entries(const size_type line_n) const
size_type calculate_line_index(const size_type line_n) const
void condense(SparseMatrix< number > &matrix) const
bool is_inhomogeneously_constrained(const size_type index) const
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 condense(const VectorType &vec_ghosted, VectorType &output) const
void copy_from(const AffineConstraints< other_number > &other)
std::vector< size_type > lines_cache
AffineConstraints(const IndexSet &local_constraints=IndexSet())
void condense(BlockSparseMatrix< number > &matrix, BlockVectorType &vector) 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
bool is_identity_constrained(const size_type line_n) const
size_type n_constraints() const
void print(std::ostream &out) const
void write_dot(std::ostream &) const
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
void set_zero(VectorType &vec) const
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
void add_lines(const std::vector< bool > &lines)
boost::iterator_range< const_iterator > LineRange
bool are_identity_constrained(const size_type line_n_1, const size_type line_n_2) const
AffineConstraints & operator=(AffineConstraints &&affine_constraints) noexcept=default
void add_entries(const size_type constrained_dof_index, const std::vector< std::pair< size_type, number >> &col_weight_pairs)
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
size_type size() const
Definition: index_set.h:1661
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1883
bool is_element(const size_type index) const
Definition: index_set.h:1776
void compress() const
Definition: index_set.h:1669
A class that provides a separate storage location on each thread that accesses the object.
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcLineInexistant(size_type arg1)
#define DeclException0(Exception0)
Definition: exceptions.h:465
static ::ExceptionBase & ExcIncorrectConstraint(int arg1, int arg2)
static ::ExceptionBase & ExcInternalError()
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:579
static ::ExceptionBase & ExcColumnNotStoredHere(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcDoFConstrainedToConstrainedDoF(int arg1, int arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcDoFIsConstrainedToConstrainedDoF(size_type arg1)
static ::ExceptionBase & ExcDoFIsConstrainedFromBothObjects(size_type arg1)
static ::ExceptionBase & ExcRowNotStoredHere(size_type arg1)
static ::ExceptionBase & ExcMatrixNotClosed()
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:533
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:510
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcMatrixIsClosed()
static ::ExceptionBase & ExcEntryAlreadyExists(size_type arg1, size_type arg2, number arg3, number arg4)
@ matrix
Contents is actually a matrix.
@ diagonal
Matrix is diagonal.
types::global_dof_index size_type
Definition: cuda_kernels.h:45
bool check(const ConstraintKinds kind_in, const unsigned int dim)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
const types::global_dof_index invalid_dof_index
Definition: types.h:233
const types::global_dof_index invalid_size_type
Definition: types.h:222
unsigned int global_dof_index
Definition: types.h:82
bool operator<(const ConstraintLine &) const
void serialize(Archive &ar, const unsigned int)
ConstraintLine(const size_type &index=numbers::invalid_dof_index, const Entries &entries={}, const number &inhomogeneity=0.0)
std::vector< std::pair< size_type, number > > Entries
std::size_t memory_consumption() const
bool operator==(const ConstraintLine &) const
ConstraintLine & operator=(const ConstraintLineType &other)
typename internal::ProductTypeImpl< typename std::decay< T >::type, typename std::decay< U >::type >::type type
static void add(const typename VectorType::value_type value, const types::global_dof_index i, VectorType &V)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
const MPI_Comm comm