Reference documentation for deal.II version GIT 5dcc62ab57 2022-07-04 21:05:02+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 
28 #include <deal.II/lac/vector.h>
30 
31 #include <boost/range/iterator_range.hpp>
32 
33 #include <algorithm>
34 #include <set>
35 #include <type_traits>
36 #include <utility>
37 #include <vector>
38 
40 
41 // Forward declarations
42 #ifndef DOXYGEN
43 template <typename>
44 class FullMatrix;
45 class SparsityPattern;
49 template <typename number>
50 class SparseMatrix;
51 template <typename number>
52 class BlockSparseMatrix;
53 
54 namespace internal
55 {
56  namespace AffineConstraints
57  {
59 
77  struct Distributing
78  {
79  Distributing(const size_type global_row = numbers::invalid_size_type,
80  const size_type local_row = numbers::invalid_size_type);
81 
82  Distributing(const Distributing &in);
83 
84  Distributing &
85  operator=(const Distributing &in);
86 
87  bool
88  operator<(const Distributing &in) const
89  {
90  return global_row < in.global_row;
91  }
92 
93  size_type global_row;
94  size_type local_row;
95  mutable size_type constraint_position;
96  };
97 
98 
99 
111  template <typename number>
112  struct DataCache
113  {
114  DataCache();
115 
116  void
117  reinit();
118 
119  size_type
120  insert_new_index(const std::pair<size_type, number> &pair);
121 
122  void
123  append_index(const size_type index,
124  const std::pair<size_type, number> &pair);
125 
126  size_type
127  get_size(const size_type index) const;
128 
129  const std::pair<size_type, number> *
130  get_entry(const size_type index) const;
131 
132  size_type row_length;
133 
134  std::vector<std::pair<size_type, number>> data;
135 
136  std::vector<size_type> individual_size;
137  };
138 
139 
140 
169  template <typename number>
170  class GlobalRowsFromLocal
171  {
172  public:
176  GlobalRowsFromLocal();
177 
178  void
179  reinit(const size_type n_local_rows);
180 
181  void
182  insert_index(const size_type global_row,
183  const size_type local_row,
184  const number constraint_value);
185  void
186  sort();
187 
188  void
189  print(std::ostream &os);
190 
194  size_type
195  size() const;
196 
201  size_type
202  size(const size_type counter_index) const;
203 
207  size_type
208  global_row(const size_type counter_index) const;
209 
213  size_type &
214  global_row(const size_type counter_index);
215 
221  size_type
222  local_row(const size_type counter_index) const;
223 
227  size_type &
228  local_row(const size_type counter_index);
229 
235  size_type
236  local_row(const size_type counter_index,
237  const size_type index_in_constraint) const;
238 
243  number
244  constraint_value(const size_type counter_index,
245  const size_type index_in_constraint) const;
246 
252  bool
253  have_indirect_rows() const;
254 
259  void
260  insert_constraint(const size_type constrained_local_dof);
261 
267  size_type
268  n_constraints() const;
269 
274  size_type
275  n_inhomogeneities() const;
276 
283  void
284  set_ith_constraint_inhomogeneous(const size_type i);
285 
290  size_type
291  constraint_origin(size_type i) const;
292 
298  std::vector<Distributing> total_row_indices;
299 
300  private:
304  DataCache<number> data_cache;
305 
310  size_type n_active_rows;
311 
316  size_type n_inhomogeneous_rows;
317  };
318 
319 
320 
338  template <typename number>
339  struct ScratchData
340  {
344  ScratchData()
345  : in_use(false)
346  {}
347 
351  ScratchData(const ScratchData &)
352  : in_use(false)
353  {}
354 
358  bool in_use;
359 
363  std::vector<size_type> columns;
364 
368  std::vector<number> values;
369 
373  std::vector<size_type> block_starts;
374 
378  std::vector<size_type> vector_indices;
379 
383  std::vector<number> vector_values;
384 
388  GlobalRowsFromLocal<number> global_rows;
389 
393  GlobalRowsFromLocal<number> global_columns;
394  };
395  } // namespace AffineConstraints
396 } // namespace internal
397 
398 namespace internal
399 {
400  namespace AffineConstraintsImplementation
401  {
402  template <class VectorType>
403  void
404  set_zero_all(const std::vector<types::global_dof_index> &cm,
405  VectorType & vec);
406 
407  template <class T>
408  void
409  set_zero_all(const std::vector<types::global_dof_index> &cm,
410  ::Vector<T> & vec);
411 
412  template <class T>
413  void
414  set_zero_all(const std::vector<types::global_dof_index> &cm,
415  ::BlockVector<T> & vec);
416  } // namespace AffineConstraintsImplementation
417 } // namespace internal
418 #endif
419 
420 // TODO[WB]: We should have a function of the kind
421 // AffineConstraints::add_constraint (const size_type constrained_dof,
422 // const std::vector<std::pair<size_type, number> > &entries,
423 // const number inhomogeneity = 0);
424 // rather than building up constraints piecemeal through add_line/add_entry
425 // etc. This would also eliminate the possibility of accidentally changing
426 // existing constraints into something pointless, see the discussion on the
427 // mailing list on "Tiny bug in interpolate_boundary_values" in Sept. 2010.
428 
502 template <typename number = double>
504 {
505 public:
510 
517  {
523 
530 
537  };
538 
556  explicit AffineConstraints(const IndexSet &local_constraints = IndexSet());
557 
561  explicit AffineConstraints(const AffineConstraints &affine_constraints);
562 
566  AffineConstraints(AffineConstraints &&affine_constraints) noexcept =
567  default; // NOLINT
568 
579  operator=(const AffineConstraints &) = delete;
580 
585  operator=(AffineConstraints &&affine_constraints) noexcept =
586  default; // NOLINT
587 
594  template <typename other_number>
595  void
597 
604  void
605  reinit(const IndexSet &local_constraints = IndexSet());
606 
613  bool
614  can_store_line(const size_type line_n) const;
615 
622  const IndexSet &
624 
644  void
646  const IndexSet & filter);
647 
657  void
658  add_line(const size_type line_n);
659 
672  void
673  add_lines(const std::vector<bool> &lines);
674 
687  void
688  add_lines(const std::set<size_type> &lines);
689 
702  void
704 
720  void
721  add_entry(const size_type constrained_dof_index,
722  const size_type column,
723  const number weight);
724 
730  void
732  const size_type constrained_dof_index,
733  const std::vector<std::pair<size_type, number>> &col_weight_pairs);
734 
746  void
747  set_inhomogeneity(const size_type constrained_dof_index, const number value);
748 
770  void
771  close();
772 
778  bool
779  is_closed() const;
780 
786  bool
787  is_closed(const MPI_Comm &comm) const;
788 
813  void
815  const AffineConstraints & other_constraints,
816  const MergeConflictBehavior merge_conflict_behavior = no_conflicts_allowed,
817  const bool allow_different_local_lines = false);
818 
831  void
832  shift(const size_type offset);
833 
841  void
842  clear();
843 
856  size_type
857  n_constraints() const;
858 
863  size_type
864  n_identities() const;
865 
871  size_type
873 
883  bool
884  is_constrained(const size_type line_n) const;
885 
897  bool
898  is_identity_constrained(const size_type line_n) const;
899 
906  bool
908  const size_type line_n_2) const;
909 
920  size_type
922 
927  bool
929 
935  bool
937 
942  const std::vector<std::pair<size_type, number>> *
943  get_constraint_entries(const size_type line_n) const;
944 
949  number
950  get_inhomogeneity(const size_type line_n) const;
951 
972  void
973  print(std::ostream &out) const;
974 
987  void
988  write_dot(std::ostream &) const;
989 
994  std::size_t
996 
1003  void
1004  resolve_indices(std::vector<types::global_dof_index> &indices) const;
1005 
1031  void
1032  condense(SparsityPattern &sparsity) const;
1033 
1037  void
1038  condense(BlockSparsityPattern &sparsity) const;
1039 
1044  void
1046 
1051  void
1053 
1061  void
1063 
1067  void
1069 
1081  template <class VectorType>
1082  void
1083  condense(VectorType &vec) const;
1084 
1091  template <class VectorType>
1092  void
1093  condense(const VectorType &vec_ghosted, VectorType &output) const;
1094 
1107  template <class VectorType>
1108  void
1109  condense(SparseMatrix<number> &matrix, VectorType &vector) const;
1110 
1115  template <class BlockVectorType>
1116  void
1117  condense(BlockSparseMatrix<number> &matrix, BlockVectorType &vector) const;
1118 
1125  template <class VectorType>
1126  void
1127  set_zero(VectorType &vec) const;
1128 
1184  template <class InVector, class OutVector>
1185  void
1186  distribute_local_to_global(const InVector & local_vector,
1187  const std::vector<size_type> &local_dof_indices,
1188  OutVector & global_vector) const;
1189 
1237  template <typename VectorType>
1238  void
1240  const std::vector<size_type> &local_dof_indices,
1241  VectorType & global_vector,
1242  const FullMatrix<number> & local_matrix) const;
1243 
1263  template <typename VectorType>
1264  void
1266  const Vector<number> & local_vector,
1267  const std::vector<size_type> &local_dof_indices_row,
1268  const std::vector<size_type> &local_dof_indices_col,
1269  VectorType & global_vector,
1270  const FullMatrix<number> & local_matrix,
1271  bool diagonal = false) const;
1272 
1276  template <class VectorType>
1277  void
1279  const number value,
1280  VectorType & global_vector) const;
1281 
1314  template <typename ForwardIteratorVec,
1315  typename ForwardIteratorInd,
1316  class VectorType>
1317  void
1318  distribute_local_to_global(ForwardIteratorVec local_vector_begin,
1319  ForwardIteratorVec local_vector_end,
1320  ForwardIteratorInd local_indices_begin,
1321  VectorType & global_vector) const;
1322 
1374  template <typename MatrixType>
1375  void
1377  const std::vector<size_type> &local_dof_indices,
1378  MatrixType & global_matrix) const;
1379 
1407  template <typename MatrixType>
1408  void
1410  const std::vector<size_type> &row_indices,
1411  const std::vector<size_type> &col_indices,
1412  MatrixType & global_matrix) const;
1413 
1430  template <typename MatrixType>
1431  void
1433  const std::vector<size_type> &row_indices,
1434  const AffineConstraints &column_affine_constraints,
1435  const std::vector<size_type> &column_indices,
1436  MatrixType & global_matrix) const;
1437 
1458  template <typename MatrixType, typename VectorType>
1459  void
1461  const Vector<number> & local_vector,
1462  const std::vector<size_type> &local_dof_indices,
1463  MatrixType & global_matrix,
1464  VectorType & global_vector,
1465  bool use_inhomogeneities_for_rhs = false) const;
1466 
1520  template <typename SparsityPatternType>
1521  void
1523  const std::vector<size_type> &local_dof_indices,
1524  SparsityPatternType & sparsity_pattern,
1525  const bool keep_constrained_entries = true,
1526  const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1527 
1531  template <typename SparsityPatternType>
1532  void
1534  const std::vector<size_type> &row_indices,
1535  const std::vector<size_type> &col_indices,
1536  SparsityPatternType & sparsity_pattern,
1537  const bool keep_constrained_entries = true,
1538  const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1539 
1544  template <typename SparsityPatternType>
1545  void
1547  const std::vector<size_type> & row_indices,
1548  const AffineConstraints<number> &col_constraints,
1549  const std::vector<size_type> & col_indices,
1550  SparsityPatternType & sparsity_pattern,
1551  const bool keep_constrained_entries = true,
1552  const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1553 
1573  template <typename ForwardIteratorVec,
1574  typename ForwardIteratorInd,
1575  class VectorType>
1576  void
1577  get_dof_values(const VectorType & global_vector,
1578  ForwardIteratorInd local_indices_begin,
1579  ForwardIteratorVec local_vector_begin,
1580  ForwardIteratorVec local_vector_end) const;
1581 
1603  template <class VectorType>
1604  void
1605  distribute(VectorType &vec) const;
1606 
1615  {
1620  using Entries = std::vector<std::pair<size_type, number>>;
1621 
1628 
1637 
1642 
1647  const Entries & entries = {},
1648  const number & inhomogeneity = 0.0);
1649 
1653  template <typename ConstraintLineType>
1654  ConstraintLine(const ConstraintLineType &other);
1655 
1659  template <typename ConstraintLineType>
1660  ConstraintLine &
1661  operator=(const ConstraintLineType &other);
1662 
1670  bool
1671  operator<(const ConstraintLine &) const;
1672 
1678  bool
1679  operator==(const ConstraintLine &) const;
1680 
1685  std::size_t
1687 
1693  template <class Archive>
1694  void
1695  serialize(Archive &ar, const unsigned int)
1696  {
1697  ar &index &entries &inhomogeneity;
1698  }
1699  };
1700 
1704  using const_iterator = typename std::vector<ConstraintLine>::const_iterator;
1705 
1709  using LineRange = boost::iterator_range<const_iterator>;
1710 
1719  const LineRange
1720  get_lines() const;
1721 
1750  bool
1751  is_consistent_in_parallel(const std::vector<IndexSet> &locally_owned_dofs,
1752  const IndexSet & locally_active_dofs,
1753  const MPI_Comm & mpi_communicator,
1754  const bool verbose = false) const;
1755 
1761  void
1762  make_consistent_in_parallel(const IndexSet &locally_owned_dofs,
1763  const IndexSet &locally_relevant_dofs,
1764  const MPI_Comm mpi_communicator);
1765 
1784  size_type,
1785  << "The specified line " << arg1 << " does not exist.");
1792  size_type,
1793  size_type,
1794  number,
1795  number,
1796  << "The entry for the indices " << arg1 << " and " << arg2
1797  << " already exists, but the values " << arg3 << " (old) and "
1798  << arg4 << " (new) differ "
1799  << "by " << (arg4 - arg3) << '.');
1806  int,
1807  int,
1808  << "You tried to constrain DoF " << arg1 << " to DoF " << arg2
1809  << ", but that one is also constrained. This is not allowed!");
1816  size_type,
1817  << "Degree of freedom " << arg1
1818  << " is constrained from both object in a merge operation.");
1825  size_type,
1826  << "In the given argument a degree of freedom is constrained "
1827  << "to another DoF with number " << arg1
1828  << ", which however is constrained by this object. This is not"
1829  << " allowed.");
1836  size_type,
1837  << "The index set given to this constraints object indicates "
1838  << "constraints for degree of freedom " << arg1
1839  << " should not be stored by this object, but a constraint "
1840  << "is being added.");
1841 
1848  size_type,
1849  size_type,
1850  << "The index set given to this constraints object indicates "
1851  << "constraints using degree of freedom " << arg2
1852  << " should not be stored by this object, but a constraint "
1853  << "for degree of freedom " << arg1 << " uses it.");
1854 
1861  int,
1862  int,
1863  << "While distributing the constraint for DoF " << arg1
1864  << ", it turns out that one of the processors "
1865  << "who own the " << arg2 << " degrees of freedom that x_"
1866  << arg1 << " is constrained against does not know about "
1867  << "the constraint on x_" << arg1
1868  << ". Did you not initialize the AffineConstraints container "
1869  << "with the appropriate locally_relevant set so "
1870  << "that every processor who owns a DoF that constrains "
1871  << "another DoF also knows about this constraint?");
1872 
1873  template <typename>
1874  friend class AffineConstraints;
1875 
1876 private:
1888  std::vector<ConstraintLine> lines;
1889 
1922  std::vector<size_type> lines_cache;
1923 
1930 
1934  bool sorted;
1935 
1937  internal::AffineConstraints::ScratchData<number>>
1939 
1944  size_type
1945  calculate_line_index(const size_type line_n) const;
1946 
1951  template <typename MatrixType, typename VectorType>
1952  void
1954  const Vector<number> & local_vector,
1955  const std::vector<size_type> &local_dof_indices,
1956  MatrixType & global_matrix,
1957  VectorType & global_vector,
1958  const bool use_inhomogeneities_for_rhs,
1959  const std::integral_constant<bool, false>) const;
1960 
1965  template <typename MatrixType, typename VectorType>
1966  void
1968  const Vector<number> & local_vector,
1969  const std::vector<size_type> &local_dof_indices,
1970  MatrixType & global_matrix,
1971  VectorType & global_vector,
1972  const bool use_inhomogeneities_for_rhs,
1973  const std::integral_constant<bool, true>) const;
1974 
1979  template <typename SparsityPatternType>
1980  void
1981  add_entries_local_to_global(const std::vector<size_type> &local_dof_indices,
1982  SparsityPatternType & sparsity_pattern,
1983  const bool keep_constrained_entries,
1984  const Table<2, bool> &dof_mask,
1985  const std::integral_constant<bool, false>) const;
1986 
1991  template <typename SparsityPatternType>
1992  void
1993  add_entries_local_to_global(const std::vector<size_type> &local_dof_indices,
1994  SparsityPatternType & sparsity_pattern,
1995  const bool keep_constrained_entries,
1996  const Table<2, bool> &dof_mask,
1997  const std::integral_constant<bool, true>) const;
1998 
2006  void
2007  make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
2008  internal::AffineConstraints::GlobalRowsFromLocal<number>
2009  &global_rows) const;
2010 
2018  void
2019  make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
2020  std::vector<size_type> & active_dofs) const;
2021 
2025  template <typename MatrixScalar, typename VectorScalar>
2028  const size_type i,
2029  const internal::AffineConstraints::GlobalRowsFromLocal<number> &global_rows,
2030  const Vector<VectorScalar> & local_vector,
2031  const std::vector<size_type> & local_dof_indices,
2032  const FullMatrix<MatrixScalar> &local_matrix) const;
2033 };
2034 
2035 /* ---------------- template and inline functions ----------------- */
2036 
2037 template <typename number>
2039  const IndexSet &local_constraints)
2040  : lines()
2041  , local_lines(local_constraints)
2042  , sorted(false)
2043 {
2044  // make sure the IndexSet is compressed. Otherwise this can lead to crashes
2045  // that are hard to find (only happen in release mode).
2046  // see tests/mpi/affine_constraints_crash_01
2048 }
2049 
2050 template <typename number>
2052  const AffineConstraints &affine_constraints)
2053  : Subscriptor()
2054  , lines(affine_constraints.lines)
2055  , lines_cache(affine_constraints.lines_cache)
2056  , local_lines(affine_constraints.local_lines)
2057  , sorted(affine_constraints.sorted)
2058 {}
2059 
2060 template <typename number>
2061 inline void
2063 {
2064  Assert(sorted == false, ExcMatrixIsClosed());
2065 
2066  // the following can happen when we compute with distributed meshes and dof
2067  // handlers and we constrain a degree of freedom whose number we don't have
2068  // locally. if we don't abort here the program will try to allocate several
2069  // terabytes of memory to resize the various arrays below :-)
2071  const size_type line_index = calculate_line_index(line_n);
2072 
2073  // check whether line already exists; it may, in which case we can just quit
2074  if (is_constrained(line_n))
2075  return;
2076 
2077  // if necessary enlarge vector of existing entries for cache
2078  if (line_index >= lines_cache.size())
2079  lines_cache.resize(std::max(2 * static_cast<size_type>(lines_cache.size()),
2080  line_index + 1),
2082 
2083  // push a new line to the end of the list
2084  lines.emplace_back();
2085  lines.back().index = line_n;
2086  lines.back().inhomogeneity = 0.;
2087  lines_cache[line_index] = lines.size() - 1;
2088 }
2089 
2090 
2091 
2092 template <typename number>
2093 inline void
2095  const size_type column,
2096  const number weight)
2097 {
2098  Assert(sorted == false, ExcMatrixIsClosed());
2099  Assert(constrained_dof_index != column,
2100  ExcMessage("Can't constrain a degree of freedom to itself"));
2101 
2102  // Ensure that the current line is present in the cache:
2103  const size_type line_index = calculate_line_index(constrained_dof_index);
2104  Assert(line_index < lines_cache.size(),
2105  ExcMessage("The current AffineConstraints does not contain the line "
2106  "for the current entry. Call AffineConstraints::add_line "
2107  "before calling this function."));
2108 
2109  // if in debug mode, check whether an entry for this column already exists
2110  // and if it's the same as the one entered at present
2111  //
2112  // in any case: exit the function if an entry for this column already
2113  // exists, since we don't want to enter it twice
2115  ExcInternalError());
2117  ExcColumnNotStoredHere(constrained_dof_index, column));
2118  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
2119  Assert(line_ptr->index == constrained_dof_index, ExcInternalError());
2120  for (const auto &p : line_ptr->entries)
2121  if (p.first == column)
2122  {
2123  Assert(std::abs(p.second - weight) < 1.e-14,
2125  constrained_dof_index, column, p.second, weight));
2126  return;
2127  }
2128 
2129  line_ptr->entries.emplace_back(column, weight);
2130 }
2131 
2132 
2133 
2134 template <typename number>
2135 inline void
2137  const size_type constrained_dof_index,
2138  const number value)
2139 {
2140  const size_type line_index = calculate_line_index(constrained_dof_index);
2141  Assert(line_index < lines_cache.size() &&
2142  lines_cache[line_index] != numbers::invalid_size_type,
2143  ExcMessage("call add_line() before calling set_inhomogeneity()"));
2144  Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
2145  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
2146  line_ptr->inhomogeneity = value;
2147 }
2148 
2149 
2150 
2151 template <typename number>
2152 template <class VectorType>
2153 inline void
2155 {
2156  // since lines is a private member, we cannot pass it to the functions
2157  // above. therefore, copy the content which is cheap
2158  std::vector<size_type> constrained_lines(lines.size());
2159  for (unsigned int i = 0; i < lines.size(); ++i)
2160  constrained_lines[i] = lines[i].index;
2161  internal::AffineConstraintsImplementation::set_zero_all(constrained_lines,
2162  vec);
2163 }
2164 
2165 template <typename number>
2168 {
2169  return lines.size();
2170 }
2171 
2172 template <typename number>
2175 {
2176  return std::count_if(lines.begin(),
2177  lines.end(),
2178  [](const ConstraintLine &line) {
2179  return (line.entries.size() == 1) &&
2180  (line.entries[0].second == number(1.));
2181  });
2182 }
2183 
2184 template <typename number>
2187 {
2188  return std::count_if(lines.begin(),
2189  lines.end(),
2190  [](const ConstraintLine &line) {
2191  return (line.inhomogeneity != number(0.));
2192  });
2193 }
2194 
2195 template <typename number>
2196 inline bool
2198 {
2199  const size_type line_index = calculate_line_index(index);
2200  return ((line_index < lines_cache.size()) &&
2201  (lines_cache[line_index] != numbers::invalid_size_type));
2202 }
2203 
2204 template <typename number>
2205 inline bool
2207  const size_type line_n) const
2208 {
2209  // check whether the entry is constrained. could use is_constrained, but
2210  // that means computing the line index twice
2211  const size_type line_index = calculate_line_index(line_n);
2212  if (line_index >= lines_cache.size() ||
2213  lines_cache[line_index] == numbers::invalid_size_type)
2214  return false;
2215  else
2216  {
2217  Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
2218  return (lines[lines_cache[line_index]].inhomogeneity != number(0.));
2219  }
2220 }
2221 
2222 template <typename number>
2223 inline const std::vector<std::pair<types::global_dof_index, number>> *
2225 {
2226  // check whether the entry is constrained. could use is_constrained, but
2227  // that means computing the line index twice
2228  const size_type line_index = calculate_line_index(line_n);
2229  if (line_index >= lines_cache.size() ||
2230  lines_cache[line_index] == numbers::invalid_size_type)
2231  return nullptr;
2232  else
2233  return &lines[lines_cache[line_index]].entries;
2234 }
2235 
2236 template <typename number>
2237 inline number
2239 {
2240  // check whether the entry is constrained. could use is_constrained, but
2241  // that means computing the line index twice
2242  const size_type line_index = calculate_line_index(line_n);
2243  if (line_index >= lines_cache.size() ||
2244  lines_cache[line_index] == numbers::invalid_size_type)
2245  return 0;
2246  else
2247  return lines[lines_cache[line_index]].inhomogeneity;
2248 }
2249 
2250 template <typename number>
2253 {
2254  // IndexSet is unused (serial case)
2255  if (local_lines.size() == 0)
2256  return line_n;
2257 
2259 
2260  return local_lines.index_within_set(line_n);
2261 }
2262 
2263 template <typename number>
2264 inline bool
2266 {
2267  return local_lines.size() == 0 || local_lines.is_element(line_n);
2268 }
2269 
2270 template <typename number>
2271 inline const IndexSet &
2273 {
2274  return local_lines;
2275 }
2276 
2277 template <typename number>
2278 template <class VectorType>
2279 inline void
2281  const size_type index,
2282  const number value,
2283  VectorType & global_vector) const
2284 {
2285  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2286 
2287  if (is_constrained(index) == false)
2288  global_vector(index) += value;
2289  else
2290  {
2291  const ConstraintLine &position =
2293  for (size_type j = 0; j < position.entries.size(); ++j)
2294  global_vector(position.entries[j].first) +=
2295  value * position.entries[j].second;
2296  }
2297 }
2298 
2299 template <typename number>
2300 template <typename ForwardIteratorVec,
2301  typename ForwardIteratorInd,
2302  class VectorType>
2303 inline void
2305  ForwardIteratorVec local_vector_begin,
2306  ForwardIteratorVec local_vector_end,
2307  ForwardIteratorInd local_indices_begin,
2308  VectorType & global_vector) const
2309 {
2310  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2311  for (; local_vector_begin != local_vector_end;
2312  ++local_vector_begin, ++local_indices_begin)
2313  {
2314  if (is_constrained(*local_indices_begin) == false)
2315  internal::ElementAccess<VectorType>::add(*local_vector_begin,
2316  *local_indices_begin,
2317  global_vector);
2318  else
2319  {
2320  const ConstraintLine &position =
2321  lines[lines_cache[calculate_line_index(*local_indices_begin)]];
2322  for (size_type j = 0; j < position.entries.size(); ++j)
2324  (*local_vector_begin) * position.entries[j].second,
2325  position.entries[j].first,
2326  global_vector);
2327  }
2328  }
2329 }
2330 
2331 template <typename number>
2332 template <class InVector, class OutVector>
2333 inline void
2335  const InVector & local_vector,
2336  const std::vector<size_type> &local_dof_indices,
2337  OutVector & global_vector) const
2338 {
2339  Assert(local_vector.size() == local_dof_indices.size(),
2340  ExcDimensionMismatch(local_vector.size(), local_dof_indices.size()));
2341  distribute_local_to_global(local_vector.begin(),
2342  local_vector.end(),
2343  local_dof_indices.begin(),
2344  global_vector);
2345 }
2346 
2347 template <typename number>
2348 template <typename ForwardIteratorVec,
2349  typename ForwardIteratorInd,
2350  class VectorType>
2351 inline void
2353  const VectorType & global_vector,
2354  ForwardIteratorInd local_indices_begin,
2355  ForwardIteratorVec local_vector_begin,
2356  ForwardIteratorVec local_vector_end) const
2357 {
2358  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2359  for (; local_vector_begin != local_vector_end;
2360  ++local_vector_begin, ++local_indices_begin)
2361  {
2362  if (is_constrained(*local_indices_begin) == false)
2363  *local_vector_begin = global_vector(*local_indices_begin);
2364  else
2365  {
2366  const ConstraintLine &position =
2367  lines[lines_cache[calculate_line_index(*local_indices_begin)]];
2368  typename VectorType::value_type value = position.inhomogeneity;
2369  for (size_type j = 0; j < position.entries.size(); ++j)
2370  value += (global_vector(position.entries[j].first) *
2371  position.entries[j].second);
2372  *local_vector_begin = value;
2373  }
2374  }
2375 }
2376 
2377 // Forward declarations
2378 #ifndef DOXYGEN
2379 template <typename MatrixType>
2380 class BlockMatrixBase;
2381 template <typename SparsityPatternType>
2383 template <typename number>
2384 class BlockSparseMatrixEZ;
2385 
2386 namespace internal
2387 {
2388  namespace AffineConstraints
2389  {
2407  template <typename MatrixType>
2408  struct IsBlockMatrix
2409  {
2410  private:
2416  template <typename T>
2417  static std::true_type
2418  check(const BlockMatrixBase<T> *);
2419 
2424  template <typename T>
2425  static std::true_type
2426  check(const BlockSparseMatrixEZ<T> *);
2427 
2432  static std::false_type
2433  check(...);
2434 
2435  public:
2442  static const bool value =
2443  std::is_same<decltype(check(std::declval<MatrixType *>())),
2444  std::true_type>::value;
2445  };
2446 
2447  // instantiation of the static member
2448  template <typename MatrixType>
2449  const bool IsBlockMatrix<MatrixType>::value;
2450 
2451 
2460  template <typename MatrixType>
2461  struct IsBlockSparsityPattern
2462  {
2463  private:
2468  template <typename T>
2469  static std::true_type
2471 
2476  static std::false_type
2477  check(...);
2478 
2479  public:
2485  static const bool value =
2486  std::is_same<decltype(check(std::declval<MatrixType *>())),
2487  std::true_type>::value;
2488  };
2489 
2490  // instantiation of the static member
2491  template <typename MatrixType>
2492  const bool IsBlockSparsityPattern<MatrixType>::value;
2493 
2494  } // namespace AffineConstraints
2495 } // namespace internal
2496 #endif
2497 
2498 
2499 
2500 template <typename number>
2501 template <typename other_number>
2502 inline void
2504  const AffineConstraints<other_number> &other)
2505 {
2506  lines.clear();
2507  lines.insert(lines.begin(), other.lines.begin(), other.lines.end());
2508  lines_cache = other.lines_cache;
2509  local_lines = other.local_lines;
2510  sorted = other.sorted;
2511 }
2512 
2513 
2514 
2515 template <typename number>
2516 template <typename MatrixType>
2517 inline void
2519  const FullMatrix<number> & local_matrix,
2520  const std::vector<size_type> &local_dof_indices,
2521  MatrixType & global_matrix) const
2522 {
2523  // create a dummy and hand on to the function actually implementing this
2524  // feature in the cm.templates.h file.
2526  distribute_local_to_global(
2527  local_matrix,
2528  dummy,
2529  local_dof_indices,
2530  global_matrix,
2531  dummy,
2532  false,
2533  std::integral_constant<
2534  bool,
2535  internal::AffineConstraints::IsBlockMatrix<MatrixType>::value>());
2536 }
2537 
2538 
2539 
2540 template <typename number>
2541 template <typename MatrixType, typename VectorType>
2542 inline void
2544  const FullMatrix<number> & local_matrix,
2545  const Vector<number> & local_vector,
2546  const std::vector<size_type> &local_dof_indices,
2547  MatrixType & global_matrix,
2548  VectorType & global_vector,
2549  bool use_inhomogeneities_for_rhs) const
2550 {
2551  // enter the internal function with the respective block information set,
2552  // the actual implementation follows in the cm.templates.h file.
2553  distribute_local_to_global(
2554  local_matrix,
2555  local_vector,
2556  local_dof_indices,
2557  global_matrix,
2558  global_vector,
2559  use_inhomogeneities_for_rhs,
2560  std::integral_constant<
2561  bool,
2562  internal::AffineConstraints::IsBlockMatrix<MatrixType>::value>());
2563 }
2564 
2565 
2566 
2567 template <typename number>
2568 template <typename SparsityPatternType>
2569 inline void
2571  const std::vector<size_type> &local_dof_indices,
2572  SparsityPatternType & sparsity_pattern,
2573  const bool keep_constrained_entries,
2574  const Table<2, bool> & dof_mask) const
2575 {
2576  // enter the internal function with the respective block information set,
2577  // the actual implementation follows in the cm.templates.h file.
2578  add_entries_local_to_global(
2579  local_dof_indices,
2580  sparsity_pattern,
2581  keep_constrained_entries,
2582  dof_mask,
2583  std::integral_constant<bool,
2584  internal::AffineConstraints::IsBlockSparsityPattern<
2585  SparsityPatternType>::value>());
2586 }
2587 
2588 
2589 
2590 template <typename number>
2592  const size_type & index,
2594  const number &inhomogeneity)
2595  : index(index)
2596  , entries(entries)
2597  , inhomogeneity(inhomogeneity)
2598 {}
2599 
2600 
2601 
2602 template <typename number>
2603 template <typename ConstraintLineType>
2605  const ConstraintLineType &other)
2606 {
2607  this->index = other.index;
2608 
2609  entries.clear();
2610  entries.insert(entries.begin(), other.entries.begin(), other.entries.end());
2611 
2612  this->inhomogeneity = other.inhomogeneity;
2613 }
2614 
2615 
2616 
2617 template <typename number>
2618 template <typename ConstraintLineType>
2621  const ConstraintLineType &other)
2622 {
2623  this->index = other.index;
2624 
2625  entries.clear();
2626  entries.insert(entries.begin(), other.entries.begin(), other.entries.end());
2627 
2628  this->inhomogeneity = other.inhomogeneity;
2629 
2630  return *this;
2631 }
2632 
2634 
2635 #endif
bool is_closed(const MPI_Comm &comm) const
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 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 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
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
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 condense(SparsityPattern &sparsity) const
void add_entries_local_to_global(const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, SparsityPatternType &sparsity_pattern, const bool keep_constrained_entries=true, const Table< 2, bool > &dof_mask=Table< 2, bool >()) const
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, SparsityPatternType &sparsity_pattern, const bool keep_constrained_entries=true, const Table< 2, bool > &dof_mask=Table< 2, bool >()) 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_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 std::integral_constant< bool, false >) 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 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
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, const std::integral_constant< bool, true >) 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
void add_entries_local_to_global(const std::vector< size_type > &local_dof_indices, SparsityPatternType &sparsity_pattern, const bool keep_constrained_entries=true, const Table< 2, bool > &dof_mask=Table< 2, bool >()) const
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:1635
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1850
bool is_element(const size_type index) const
Definition: index_set.h:1743
void compress() const
Definition: index_set.h:1643
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:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcLineInexistant(size_type arg1)
#define DeclException0(Exception0)
Definition: exceptions.h:464
static ::ExceptionBase & ExcIncorrectConstraint(int arg1, int arg2)
static ::ExceptionBase & ExcInternalError()
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:578
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:1473
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:532
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
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:216
const types::global_dof_index invalid_size_type
Definition: types.h:210
unsigned int global_dof_index
Definition: types.h:76
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