Reference documentation for deal.II version 9.3.0
\(\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 - 2021 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 <set>
34 #include <type_traits>
35 #include <utility>
36 #include <vector>
37 
39 
40 // Forward declarations
41 #ifndef DOXYGEN
42 template <typename>
43 class FullMatrix;
44 class SparsityPattern;
48 template <typename number>
49 class SparseMatrix;
50 template <typename number>
51 class BlockSparseMatrix;
52 
53 namespace internal
54 {
55  namespace AffineConstraints
56  {
58 
76  struct Distributing
77  {
78  Distributing(const size_type global_row = numbers::invalid_size_type,
79  const size_type local_row = numbers::invalid_size_type);
80 
81  Distributing(const Distributing &in);
82 
83  Distributing &
84  operator=(const Distributing &in);
85 
86  bool
87  operator<(const Distributing &in) const
88  {
89  return global_row < in.global_row;
90  }
91 
92  size_type global_row;
93  size_type local_row;
94  mutable size_type constraint_position;
95  };
96 
97 
98 
110  template <typename number>
111  struct DataCache
112  {
113  DataCache();
114 
115  void
116  reinit();
117 
118  size_type
119  insert_new_index(const std::pair<size_type, number> &pair);
120 
121  void
122  append_index(const size_type index,
123  const std::pair<size_type, number> &pair);
124 
125  size_type
126  get_size(const size_type index) const;
127 
128  const std::pair<size_type, number> *
129  get_entry(const size_type index) const;
130 
131  size_type row_length;
132 
133  std::vector<std::pair<size_type, number>> data;
134 
135  std::vector<size_type> individual_size;
136  };
137 
138 
139 
168  template <typename number>
169  class GlobalRowsFromLocal
170  {
171  public:
175  GlobalRowsFromLocal();
176 
177  void
178  reinit(const size_type n_local_rows);
179 
180  void
181  insert_index(const size_type global_row,
182  const size_type local_row,
183  const number constraint_value);
184  void
185  sort();
186 
187  void
188  print(std::ostream &os);
189 
193  size_type
194  size() const;
195 
200  size_type
201  size(const size_type counter_index) const;
202 
206  size_type
207  global_row(const size_type counter_index) const;
208 
212  size_type &
213  global_row(const size_type counter_index);
214 
220  size_type
221  local_row(const size_type counter_index) const;
222 
226  size_type &
227  local_row(const size_type counter_index);
228 
234  size_type
235  local_row(const size_type counter_index,
236  const size_type index_in_constraint) const;
237 
242  number
243  constraint_value(const size_type counter_index,
244  const size_type index_in_constraint) const;
245 
251  bool
252  have_indirect_rows() const;
253 
258  void
259  insert_constraint(const size_type constrained_local_dof);
260 
266  size_type
267  n_constraints() const;
268 
273  size_type
274  n_inhomogeneities() const;
275 
282  void
283  set_ith_constraint_inhomogeneous(const size_type i);
284 
289  size_type
290  constraint_origin(size_type i) const;
291 
297  std::vector<Distributing> total_row_indices;
298 
299  private:
303  DataCache<number> data_cache;
304 
309  size_type n_active_rows;
310 
315  size_type n_inhomogeneous_rows;
316  };
317 
318 
319 
337  template <typename number>
338  struct ScratchData
339  {
343  ScratchData()
344  : in_use(false)
345  {}
346 
350  ScratchData(const ScratchData &)
351  : in_use(false)
352  {}
353 
357  bool in_use;
358 
362  std::vector<size_type> columns;
363 
367  std::vector<number> values;
368 
372  std::vector<size_type> block_starts;
373 
377  std::vector<size_type> vector_indices;
378 
382  std::vector<number> vector_values;
383 
387  GlobalRowsFromLocal<number> global_rows;
388 
392  GlobalRowsFromLocal<number> global_columns;
393  };
394  } // namespace AffineConstraints
395 } // namespace internal
396 
397 namespace internal
398 {
399  namespace AffineConstraintsImplementation
400  {
401  template <class VectorType>
402  void
403  set_zero_all(const std::vector<types::global_dof_index> &cm,
404  VectorType & vec);
405 
406  template <class T>
407  void
408  set_zero_all(const std::vector<types::global_dof_index> &cm,
409  ::Vector<T> & vec);
410 
411  template <class T>
412  void
413  set_zero_all(const std::vector<types::global_dof_index> &cm,
414  ::BlockVector<T> & vec);
415  } // namespace AffineConstraintsImplementation
416 } // namespace internal
417 
418 
419 template <typename number>
420 class AffineConstraints;
421 #endif
422 
423 // TODO[WB]: We should have a function of the kind
424 // AffineConstraints::add_constraint (const size_type constrained_dof,
425 // const std::vector<std::pair<size_type, number> > &entries,
426 // const number inhomogeneity = 0);
427 // rather than building up constraints piecemeal through add_line/add_entry
428 // etc. This would also eliminate the possibility of accidentally changing
429 // existing constraints into something pointless, see the discussion on the
430 // mailing list on "Tiny bug in interpolate_boundary_values" in Sept. 2010.
431 
505 template <typename number = double>
507 {
508 public:
513 
520  {
526 
533 
539  right_object_wins
540  };
541 
559  explicit AffineConstraints(const IndexSet &local_constraints = IndexSet());
560 
564  explicit AffineConstraints(const AffineConstraints &affine_constraints);
565 
569  AffineConstraints(AffineConstraints &&affine_constraints) noexcept =
570  default; // NOLINT
571 
582  operator=(const AffineConstraints &) = delete;
583 
588  operator=(AffineConstraints &&affine_constraints) noexcept =
589  default; // NOLINT
590 
597  template <typename other_number>
598  void
599  copy_from(const AffineConstraints<other_number> &other);
600 
607  void
608  reinit(const IndexSet &local_constraints = IndexSet());
609 
616  bool
617  can_store_line(const size_type line_n) const;
618 
625  const IndexSet &
626  get_local_lines() const;
627 
647  void
648  add_selected_constraints(const AffineConstraints &constraints_in,
649  const IndexSet & filter);
650 
660  void
661  add_line(const size_type line_n);
662 
675  void
676  add_lines(const std::vector<bool> &lines);
677 
690  void
691  add_lines(const std::set<size_type> &lines);
692 
705  void
706  add_lines(const IndexSet &lines);
707 
723  void
724  add_entry(const size_type constrained_dof_index,
725  const size_type column,
726  const number weight);
727 
733  void
734  add_entries(
735  const size_type constrained_dof_index,
736  const std::vector<std::pair<size_type, number>> &col_weight_pairs);
737 
749  void
750  set_inhomogeneity(const size_type constrained_dof_index, const number value);
751 
773  void
774  close();
775 
800  void
801  merge(
802  const AffineConstraints & other_constraints,
803  const MergeConflictBehavior merge_conflict_behavior = no_conflicts_allowed,
804  const bool allow_different_local_lines = false);
805 
818  void
819  shift(const size_type offset);
820 
828  void
829  clear();
830 
843  size_type
844  n_constraints() const;
845 
855  bool
856  is_constrained(const size_type line_n) const;
857 
869  bool
870  is_identity_constrained(const size_type line_n) const;
871 
878  bool
879  are_identity_constrained(const size_type line_n_1,
880  const size_type line_n_2) const;
881 
892  size_type
893  max_constraint_indirections() const;
894 
899  bool
900  is_inhomogeneously_constrained(const size_type index) const;
901 
907  bool
908  has_inhomogeneities() const;
909 
914  const std::vector<std::pair<size_type, number>> *
915  get_constraint_entries(const size_type line_n) const;
916 
921  number
922  get_inhomogeneity(const size_type line_n) const;
923 
944  void
945  print(std::ostream &out) const;
946 
959  void
960  write_dot(std::ostream &) const;
961 
966  std::size_t
967  memory_consumption() const;
968 
975  void
976  resolve_indices(std::vector<types::global_dof_index> &indices) const;
977 
1003  void
1004  condense(SparsityPattern &sparsity) const;
1005 
1009  void
1010  condense(BlockSparsityPattern &sparsity) const;
1011 
1016  void
1017  condense(DynamicSparsityPattern &sparsity) const;
1018 
1023  void
1024  condense(BlockDynamicSparsityPattern &sparsity) const;
1025 
1033  void
1034  condense(SparseMatrix<number> &matrix) const;
1035 
1039  void
1040  condense(BlockSparseMatrix<number> &matrix) const;
1041 
1053  template <class VectorType>
1054  void
1055  condense(VectorType &vec) const;
1056 
1063  template <class VectorType>
1064  void
1065  condense(const VectorType &vec_ghosted, VectorType &output) const;
1066 
1079  template <class VectorType>
1080  void
1081  condense(SparseMatrix<number> &matrix, VectorType &vector) const;
1082 
1087  template <class BlockVectorType>
1088  void
1089  condense(BlockSparseMatrix<number> &matrix, BlockVectorType &vector) const;
1090 
1097  template <class VectorType>
1098  void
1099  set_zero(VectorType &vec) const;
1100 
1156  template <class InVector, class OutVector>
1157  void
1158  distribute_local_to_global(const InVector & local_vector,
1159  const std::vector<size_type> &local_dof_indices,
1160  OutVector & global_vector) const;
1161 
1209  template <typename VectorType>
1210  void
1211  distribute_local_to_global(const Vector<number> & local_vector,
1212  const std::vector<size_type> &local_dof_indices,
1213  VectorType & global_vector,
1214  const FullMatrix<number> & local_matrix) const;
1215 
1235  template <typename VectorType>
1236  void
1237  distribute_local_to_global(
1238  const Vector<number> & local_vector,
1239  const std::vector<size_type> &local_dof_indices_row,
1240  const std::vector<size_type> &local_dof_indices_col,
1241  VectorType & global_vector,
1242  const FullMatrix<number> & local_matrix,
1243  bool diagonal = false) const;
1244 
1248  template <class VectorType>
1249  void
1250  distribute_local_to_global(const size_type index,
1251  const number value,
1252  VectorType & global_vector) const;
1253 
1286  template <typename ForwardIteratorVec,
1287  typename ForwardIteratorInd,
1288  class VectorType>
1289  void
1290  distribute_local_to_global(ForwardIteratorVec local_vector_begin,
1291  ForwardIteratorVec local_vector_end,
1292  ForwardIteratorInd local_indices_begin,
1293  VectorType & global_vector) const;
1294 
1346  template <typename MatrixType>
1347  void
1348  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1349  const std::vector<size_type> &local_dof_indices,
1350  MatrixType & global_matrix) const;
1351 
1379  template <typename MatrixType>
1380  void
1381  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1382  const std::vector<size_type> &row_indices,
1383  const std::vector<size_type> &col_indices,
1384  MatrixType & global_matrix) const;
1385 
1402  template <typename MatrixType>
1403  void
1404  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1405  const std::vector<size_type> &row_indices,
1406  const AffineConstraints &column_affine_constraints,
1407  const std::vector<size_type> &column_indices,
1408  MatrixType & global_matrix) const;
1409 
1430  template <typename MatrixType, typename VectorType>
1431  void
1432  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1433  const Vector<number> & local_vector,
1434  const std::vector<size_type> &local_dof_indices,
1435  MatrixType & global_matrix,
1436  VectorType & global_vector,
1437  bool use_inhomogeneities_for_rhs = false) const;
1438 
1492  template <typename SparsityPatternType>
1493  void
1494  add_entries_local_to_global(
1495  const std::vector<size_type> &local_dof_indices,
1496  SparsityPatternType & sparsity_pattern,
1497  const bool keep_constrained_entries = true,
1498  const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1499 
1503  template <typename SparsityPatternType>
1504  void
1505  add_entries_local_to_global(
1506  const std::vector<size_type> &row_indices,
1507  const std::vector<size_type> &col_indices,
1508  SparsityPatternType & sparsity_pattern,
1509  const bool keep_constrained_entries = true,
1510  const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1511 
1531  template <typename ForwardIteratorVec,
1532  typename ForwardIteratorInd,
1533  class VectorType>
1534  void
1535  get_dof_values(const VectorType & global_vector,
1536  ForwardIteratorInd local_indices_begin,
1537  ForwardIteratorVec local_vector_begin,
1538  ForwardIteratorVec local_vector_end) const;
1539 
1561  template <class VectorType>
1562  void
1563  distribute(VectorType &vec) const;
1564 
1573  {
1578  using Entries = std::vector<std::pair<size_type, number>>;
1579 
1586 
1595 
1600 
1605  const Entries & entries = {},
1606  const number & inhomogeneity = 0.0);
1607 
1611  template <typename ConstraintLineType>
1612  ConstraintLine(const ConstraintLineType &other);
1613 
1617  template <typename ConstraintLineType>
1618  ConstraintLine &
1619  operator=(const ConstraintLineType &other);
1620 
1628  bool
1629  operator<(const ConstraintLine &) const;
1630 
1636  bool
1637  operator==(const ConstraintLine &) const;
1638 
1643  std::size_t
1644  memory_consumption() const;
1645 
1651  template <class Archive>
1652  void
1653  serialize(Archive &ar, const unsigned int)
1654  {
1655  ar &index &entries &inhomogeneity;
1656  }
1657  };
1658 
1662  using const_iterator = typename std::vector<ConstraintLine>::const_iterator;
1663 
1667  using LineRange = boost::iterator_range<const_iterator>;
1668 
1677  const LineRange
1678  get_lines() const;
1679 
1708  bool
1709  is_consistent_in_parallel(const std::vector<IndexSet> &locally_owned_dofs,
1710  const IndexSet & locally_active_dofs,
1711  const MPI_Comm & mpi_communicator,
1712  const bool verbose = false) const;
1713 
1719  DeclException0(ExcMatrixIsClosed);
1725  DeclException0(ExcMatrixNotClosed);
1732  size_type,
1733  << "The specified line " << arg1 << " does not exist.");
1739  DeclException4(ExcEntryAlreadyExists,
1740  size_type,
1741  size_type,
1742  number,
1743  number,
1744  << "The entry for the indices " << arg1 << " and " << arg2
1745  << " already exists, but the values " << arg3 << " (old) and "
1746  << arg4 << " (new) differ "
1747  << "by " << (arg4 - arg3) << ".");
1753  DeclException2(ExcDoFConstrainedToConstrainedDoF,
1754  int,
1755  int,
1756  << "You tried to constrain DoF " << arg1 << " to DoF " << arg2
1757  << ", but that one is also constrained. This is not allowed!");
1763  DeclException1(ExcDoFIsConstrainedFromBothObjects,
1764  size_type,
1765  << "Degree of freedom " << arg1
1766  << " is constrained from both object in a merge operation.");
1772  DeclException1(ExcDoFIsConstrainedToConstrainedDoF,
1773  size_type,
1774  << "In the given argument a degree of freedom is constrained "
1775  << "to another DoF with number " << arg1
1776  << ", which however is constrained by this object. This is not"
1777  << " allowed.");
1783  DeclException1(ExcRowNotStoredHere,
1784  size_type,
1785  << "The index set given to this constraints object indicates "
1786  << "constraints for degree of freedom " << arg1
1787  << " should not be stored by this object, but a constraint "
1788  << "is being added.");
1789 
1795  DeclException2(ExcColumnNotStoredHere,
1796  size_type,
1797  size_type,
1798  << "The index set given to this constraints object indicates "
1799  << "constraints using degree of freedom " << arg2
1800  << " should not be stored by this object, but a constraint "
1801  << "for degree of freedom " << arg1 << " uses it.");
1802 
1808  DeclException2(ExcIncorrectConstraint,
1809  int,
1810  int,
1811  << "While distributing the constraint for DoF " << arg1
1812  << ", it turns out that one of the processors "
1813  << "who own the " << arg2 << " degrees of freedom that x_"
1814  << arg1 << " is constrained against does not know about "
1815  << "the constraint on x_" << arg1
1816  << ". Did you not initialize the AffineConstraints container "
1817  << "with the appropriate locally_relevant set so "
1818  << "that every processor who owns a DoF that constrains "
1819  << "another DoF also knows about this constraint?");
1820 
1821  template <typename>
1822  friend class AffineConstraints;
1823 
1824 private:
1836  std::vector<ConstraintLine> lines;
1837 
1870  std::vector<size_type> lines_cache;
1871 
1878 
1882  bool sorted;
1883 
1885  internal::AffineConstraints::ScratchData<number>>
1887 
1892  size_type
1893  calculate_line_index(const size_type line_n) const;
1894 
1899  template <typename MatrixType, typename VectorType>
1900  void
1901  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1902  const Vector<number> & local_vector,
1903  const std::vector<size_type> &local_dof_indices,
1904  MatrixType & global_matrix,
1905  VectorType & global_vector,
1906  const bool use_inhomogeneities_for_rhs,
1907  const std::integral_constant<bool, false>) const;
1908 
1913  template <typename MatrixType, typename VectorType>
1914  void
1915  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1916  const Vector<number> & local_vector,
1917  const std::vector<size_type> &local_dof_indices,
1918  MatrixType & global_matrix,
1919  VectorType & global_vector,
1920  const bool use_inhomogeneities_for_rhs,
1921  const std::integral_constant<bool, true>) const;
1922 
1927  template <typename SparsityPatternType>
1928  void
1929  add_entries_local_to_global(const std::vector<size_type> &local_dof_indices,
1930  SparsityPatternType & sparsity_pattern,
1931  const bool keep_constrained_entries,
1932  const Table<2, bool> &dof_mask,
1933  const std::integral_constant<bool, false>) const;
1934 
1939  template <typename SparsityPatternType>
1940  void
1941  add_entries_local_to_global(const std::vector<size_type> &local_dof_indices,
1942  SparsityPatternType & sparsity_pattern,
1943  const bool keep_constrained_entries,
1944  const Table<2, bool> &dof_mask,
1945  const std::integral_constant<bool, true>) const;
1946 
1954  void
1955  make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
1956  internal::AffineConstraints::GlobalRowsFromLocal<number>
1957  &global_rows) const;
1958 
1966  void
1967  make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
1968  std::vector<size_type> & active_dofs) const;
1969 
1973  template <typename MatrixScalar, typename VectorScalar>
1976  const size_type i,
1977  const internal::AffineConstraints::GlobalRowsFromLocal<number> &global_rows,
1978  const Vector<VectorScalar> & local_vector,
1979  const std::vector<size_type> & local_dof_indices,
1980  const FullMatrix<MatrixScalar> &local_matrix) const;
1981 };
1982 
1983 /* ---------------- template and inline functions ----------------- */
1984 
1985 template <typename number>
1987  const IndexSet &local_constraints)
1988  : lines()
1989  , local_lines(local_constraints)
1990  , sorted(false)
1991 {
1992  // make sure the IndexSet is compressed. Otherwise this can lead to crashes
1993  // that are hard to find (only happen in release mode).
1994  // see tests/mpi/affine_constraints_crash_01
1996 }
1997 
1998 template <typename number>
2000  const AffineConstraints &affine_constraints)
2001  : Subscriptor()
2002  , lines(affine_constraints.lines)
2003  , lines_cache(affine_constraints.lines_cache)
2004  , local_lines(affine_constraints.local_lines)
2005  , sorted(affine_constraints.sorted)
2006 {}
2007 
2008 template <typename number>
2009 inline void
2011 {
2012  Assert(sorted == false, ExcMatrixIsClosed());
2013 
2014  // the following can happen when we compute with distributed meshes and dof
2015  // handlers and we constrain a degree of freedom whose number we don't have
2016  // locally. if we don't abort here the program will try to allocate several
2017  // terabytes of memory to resize the various arrays below :-)
2019  const size_type line_index = calculate_line_index(line_n);
2020 
2021  // check whether line already exists; it may, in which case we can just quit
2022  if (is_constrained(line_n))
2023  return;
2024 
2025  // if necessary enlarge vector of existing entries for cache
2026  if (line_index >= lines_cache.size())
2027  lines_cache.resize(std::max(2 * static_cast<size_type>(lines_cache.size()),
2028  line_index + 1),
2030 
2031  // push a new line to the end of the list
2032  lines.emplace_back();
2033  lines.back().index = line_n;
2034  lines.back().inhomogeneity = 0.;
2035  lines_cache[line_index] = lines.size() - 1;
2036 }
2037 
2038 
2039 
2040 template <typename number>
2041 inline void
2043  const size_type column,
2044  const number weight)
2045 {
2046  Assert(sorted == false, ExcMatrixIsClosed());
2047  Assert(constrained_dof_index != column,
2048  ExcMessage("Can't constrain a degree of freedom to itself"));
2049 
2050  // Ensure that the current line is present in the cache:
2051  const size_type line_index = calculate_line_index(constrained_dof_index);
2052  Assert(line_index < lines_cache.size(),
2053  ExcMessage("The current AffineConstraints does not contain the line "
2054  "for the current entry. Call AffineConstraints::add_line "
2055  "before calling this function."));
2056 
2057  // if in debug mode, check whether an entry for this column already exists
2058  // and if it's the same as the one entered at present
2059  //
2060  // in any case: exit the function if an entry for this column already
2061  // exists, since we don't want to enter it twice
2063  ExcInternalError());
2065  ExcColumnNotStoredHere(constrained_dof_index, column));
2066  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
2067  Assert(line_ptr->index == constrained_dof_index, ExcInternalError());
2068  for (const auto &p : line_ptr->entries)
2069  if (p.first == column)
2070  {
2071  Assert(std::abs(p.second - weight) < 1.e-14,
2073  constrained_dof_index, column, p.second, weight));
2074  return;
2075  }
2076 
2077  line_ptr->entries.emplace_back(column, weight);
2078 }
2079 
2080 
2081 
2082 template <typename number>
2083 inline void
2085  const size_type constrained_dof_index,
2086  const number value)
2087 {
2088  const size_type line_index = calculate_line_index(constrained_dof_index);
2089  Assert(line_index < lines_cache.size() &&
2090  lines_cache[line_index] != numbers::invalid_size_type,
2091  ExcMessage("call add_line() before calling set_inhomogeneity()"));
2092  Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
2093  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
2094  line_ptr->inhomogeneity = value;
2095 }
2096 
2097 
2098 
2099 template <typename number>
2100 template <class VectorType>
2101 inline void
2103 {
2104  // since lines is a private member, we cannot pass it to the functions
2105  // above. therefore, copy the content which is cheap
2106  std::vector<size_type> constrained_lines(lines.size());
2107  for (unsigned int i = 0; i < lines.size(); ++i)
2108  constrained_lines[i] = lines[i].index;
2109  internal::AffineConstraintsImplementation::set_zero_all(constrained_lines,
2110  vec);
2111 }
2112 
2113 template <typename number>
2116 {
2117  return lines.size();
2118 }
2119 
2120 template <typename number>
2121 inline bool
2123 {
2124  const size_type line_index = calculate_line_index(index);
2125  return ((line_index < lines_cache.size()) &&
2126  (lines_cache[line_index] != numbers::invalid_size_type));
2127 }
2128 
2129 template <typename number>
2130 inline bool
2132  const size_type line_n) const
2133 {
2134  // check whether the entry is constrained. could use is_constrained, but
2135  // that means computing the line index twice
2136  const size_type line_index = calculate_line_index(line_n);
2137  if (line_index >= lines_cache.size() ||
2138  lines_cache[line_index] == numbers::invalid_size_type)
2139  return false;
2140  else
2141  {
2142  Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
2143  return !(lines[lines_cache[line_index]].inhomogeneity == number(0.));
2144  }
2145 }
2146 
2147 template <typename number>
2148 inline const std::vector<std::pair<types::global_dof_index, number>> *
2150 {
2151  // check whether the entry is constrained. could use is_constrained, but
2152  // that means computing the line index twice
2153  const size_type line_index = calculate_line_index(line_n);
2154  if (line_index >= lines_cache.size() ||
2155  lines_cache[line_index] == numbers::invalid_size_type)
2156  return nullptr;
2157  else
2158  return &lines[lines_cache[line_index]].entries;
2159 }
2160 
2161 template <typename number>
2162 inline number
2164 {
2165  // check whether the entry is constrained. could use is_constrained, but
2166  // that means computing the line index twice
2167  const size_type line_index = calculate_line_index(line_n);
2168  if (line_index >= lines_cache.size() ||
2169  lines_cache[line_index] == numbers::invalid_size_type)
2170  return 0;
2171  else
2172  return lines[lines_cache[line_index]].inhomogeneity;
2173 }
2174 
2175 template <typename number>
2178 {
2179  // IndexSet is unused (serial case)
2180  if (!local_lines.size())
2181  return line_n;
2182 
2184 
2185  return local_lines.index_within_set(line_n);
2186 }
2187 
2188 template <typename number>
2189 inline bool
2191 {
2192  return !local_lines.size() || local_lines.is_element(line_n);
2193 }
2194 
2195 template <typename number>
2196 inline const IndexSet &
2198 {
2199  return local_lines;
2200 }
2201 
2202 template <typename number>
2203 template <class VectorType>
2204 inline void
2206  const size_type index,
2207  const number value,
2208  VectorType & global_vector) const
2209 {
2210  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2211 
2212  if (is_constrained(index) == false)
2213  global_vector(index) += value;
2214  else
2215  {
2216  const ConstraintLine &position =
2218  for (size_type j = 0; j < position.entries.size(); ++j)
2219  global_vector(position.entries[j].first) +=
2220  value * position.entries[j].second;
2221  }
2222 }
2223 
2224 template <typename number>
2225 template <typename ForwardIteratorVec,
2226  typename ForwardIteratorInd,
2227  class VectorType>
2228 inline void
2230  ForwardIteratorVec local_vector_begin,
2231  ForwardIteratorVec local_vector_end,
2232  ForwardIteratorInd local_indices_begin,
2233  VectorType & global_vector) const
2234 {
2235  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2236  for (; local_vector_begin != local_vector_end;
2237  ++local_vector_begin, ++local_indices_begin)
2238  {
2239  if (is_constrained(*local_indices_begin) == false)
2240  internal::ElementAccess<VectorType>::add(*local_vector_begin,
2241  *local_indices_begin,
2242  global_vector);
2243  else
2244  {
2245  const ConstraintLine &position =
2246  lines[lines_cache[calculate_line_index(*local_indices_begin)]];
2247  for (size_type j = 0; j < position.entries.size(); ++j)
2249  (*local_vector_begin) * position.entries[j].second,
2250  position.entries[j].first,
2251  global_vector);
2252  }
2253  }
2254 }
2255 
2256 template <typename number>
2257 template <class InVector, class OutVector>
2258 inline void
2260  const InVector & local_vector,
2261  const std::vector<size_type> &local_dof_indices,
2262  OutVector & global_vector) const
2263 {
2264  Assert(local_vector.size() == local_dof_indices.size(),
2265  ExcDimensionMismatch(local_vector.size(), local_dof_indices.size()));
2266  distribute_local_to_global(local_vector.begin(),
2267  local_vector.end(),
2268  local_dof_indices.begin(),
2269  global_vector);
2270 }
2271 
2272 template <typename number>
2273 template <typename ForwardIteratorVec,
2274  typename ForwardIteratorInd,
2275  class VectorType>
2276 inline void
2278  const VectorType & global_vector,
2279  ForwardIteratorInd local_indices_begin,
2280  ForwardIteratorVec local_vector_begin,
2281  ForwardIteratorVec local_vector_end) const
2282 {
2283  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2284  for (; local_vector_begin != local_vector_end;
2285  ++local_vector_begin, ++local_indices_begin)
2286  {
2287  if (is_constrained(*local_indices_begin) == false)
2288  *local_vector_begin = global_vector(*local_indices_begin);
2289  else
2290  {
2291  const ConstraintLine &position =
2292  lines[lines_cache[calculate_line_index(*local_indices_begin)]];
2293  typename VectorType::value_type value = position.inhomogeneity;
2294  for (size_type j = 0; j < position.entries.size(); ++j)
2295  value += (global_vector(position.entries[j].first) *
2296  position.entries[j].second);
2297  *local_vector_begin = value;
2298  }
2299  }
2300 }
2301 
2302 template <typename MatrixType>
2304 template <typename SparsityPatternType>
2306 template <typename number>
2308 
2309 namespace internal
2310 {
2312  {
2330  template <typename MatrixType>
2332  {
2333  private:
2339  template <typename T>
2340  static std::true_type
2341  check(const BlockMatrixBase<T> *);
2342 
2347  template <typename T>
2348  static std::true_type
2349  check(const BlockSparseMatrixEZ<T> *);
2350 
2355  static std::false_type
2356  check(...);
2357 
2358  public:
2365  static const bool value =
2366  std::is_same<decltype(check(std::declval<MatrixType *>())),
2367  std::true_type>::value;
2368  };
2369 
2370  // instantiation of the static member
2371  template <typename MatrixType>
2373 
2374 
2383  template <typename MatrixType>
2385  {
2386  private:
2391  template <typename T>
2392  static std::true_type
2393  check(const BlockSparsityPatternBase<T> *);
2394 
2399  static std::false_type
2400  check(...);
2401 
2402  public:
2408  static const bool value =
2409  std::is_same<decltype(check(std::declval<MatrixType *>())),
2410  std::true_type>::value;
2411  };
2412 
2413  // instantiation of the static member
2414  template <typename MatrixType>
2416 
2417  } // namespace AffineConstraints
2418 } // namespace internal
2419 
2420 
2421 
2422 template <typename number>
2423 template <typename other_number>
2424 inline void
2426  const AffineConstraints<other_number> &other)
2427 {
2428  lines.clear();
2429  lines.insert(lines.begin(), other.lines.begin(), other.lines.end());
2430  lines_cache = other.lines_cache;
2431  local_lines = other.local_lines;
2432  sorted = other.sorted;
2433 }
2434 
2435 
2436 
2437 template <typename number>
2438 template <typename MatrixType>
2439 inline void
2441  const FullMatrix<number> & local_matrix,
2442  const std::vector<size_type> &local_dof_indices,
2443  MatrixType & global_matrix) const
2444 {
2445  // create a dummy and hand on to the function actually implementing this
2446  // feature in the cm.templates.h file.
2447  Vector<typename MatrixType::value_type> dummy(0);
2449  local_matrix,
2450  dummy,
2451  local_dof_indices,
2452  global_matrix,
2453  dummy,
2454  false,
2455  std::integral_constant<
2456  bool,
2458 }
2459 
2460 
2461 
2462 template <typename number>
2463 template <typename MatrixType, typename VectorType>
2464 inline void
2466  const FullMatrix<number> & local_matrix,
2467  const Vector<number> & local_vector,
2468  const std::vector<size_type> &local_dof_indices,
2469  MatrixType & global_matrix,
2470  VectorType & global_vector,
2471  bool use_inhomogeneities_for_rhs) const
2472 {
2473  // enter the internal function with the respective block information set,
2474  // the actual implementation follows in the cm.templates.h file.
2476  local_matrix,
2477  local_vector,
2478  local_dof_indices,
2479  global_matrix,
2480  global_vector,
2481  use_inhomogeneities_for_rhs,
2482  std::integral_constant<
2483  bool,
2485 }
2486 
2487 
2488 
2489 template <typename number>
2490 template <typename SparsityPatternType>
2491 inline void
2493  const std::vector<size_type> &local_dof_indices,
2494  SparsityPatternType & sparsity_pattern,
2495  const bool keep_constrained_entries,
2496  const Table<2, bool> & dof_mask) const
2497 {
2498  // enter the internal function with the respective block information set,
2499  // the actual implementation follows in the cm.templates.h file.
2501  local_dof_indices,
2502  sparsity_pattern,
2503  keep_constrained_entries,
2504  dof_mask,
2505  std::integral_constant<bool,
2507  SparsityPatternType>::value>());
2508 }
2509 
2510 
2511 
2512 template <typename number>
2514  const size_type & index,
2516  const number &inhomogeneity)
2517  : index(index)
2518  , entries(entries)
2519  , inhomogeneity(inhomogeneity)
2520 {}
2521 
2522 
2523 
2524 template <typename number>
2525 template <typename ConstraintLineType>
2527  const ConstraintLineType &other)
2528 {
2529  this->index = other.index;
2530 
2531  entries.clear();
2532  entries.insert(entries.begin(), other.entries.begin(), other.entries.end());
2533 
2534  this->inhomogeneity = other.inhomogeneity;
2535 }
2536 
2537 
2538 
2539 template <typename number>
2540 template <typename ConstraintLineType>
2543 operator=(const ConstraintLineType &other)
2544 {
2545  this->index = other.index;
2546 
2547  entries.clear();
2548  entries.insert(entries.begin(), other.entries.begin(), other.entries.end());
2549 
2550  this->inhomogeneity = other.inhomogeneity;
2551 
2552  return *this;
2553 }
2554 
2556 
2557 #endif
const types::global_dof_index invalid_size_type
Definition: types.h:205
void add_entry(const size_type constrained_dof_index, const size_type column, const number weight)
std::vector< ConstraintLine > lines
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
void set_zero(VectorType &vec) const
bool is_constrained(const size_type line_n) const
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:538
bool can_store_line(const size_type line_n) const
A class that provides a separate storage location on each thread that accesses the object...
Contents is actually a matrix.
types::global_dof_index size_type
Definition: cuda_kernels.h:45
typename internal::ProductTypeImpl< typename std::decay< T >::type, typename std::decay< U >::type >::type type
void copy_from(const AffineConstraints< other_number > &other)
static ::ExceptionBase & ExcColumnNotStoredHere(size_type arg1, size_type arg2)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
static ::ExceptionBase & ExcRowNotStoredHere(size_type arg1)
void distribute_local_to_global(const InVector &local_vector, const std::vector< size_type > &local_dof_indices, OutVector &global_vector) const
static ::ExceptionBase & ExcMatrixNotClosed()
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
static ::ExceptionBase & ExcMatrixIsClosed()
size_type size() const
Definition: index_set.h:1634
Matrix is diagonal.
void serialize(Archive &ar, const unsigned int)
size_type n_constraints() const
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:515
#define Assert(cond, exc)
Definition: exceptions.h:1465
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
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void make_sorted_row_list(const std::vector< size_type > &local_dof_indices, internal::AffineConstraints::GlobalRowsFromLocal< number > &global_rows) const
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1921
#define DeclException0(Exception0)
Definition: exceptions.h:470
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
ConstraintLine & operator=(const ConstraintLineType &other)
constexpr TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
friend class AffineConstraints
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
size_type calculate_line_index(const size_type line_n) const
std::vector< std::pair< size_type, number > > Entries
const IndexSet & get_local_lines() const
Threads::ThreadLocalStorage< internal::AffineConstraints::ScratchData< number > > scratch_data
void set_inhomogeneity(const size_type constrained_dof_index, const number value)
static ::ExceptionBase & ExcLineInexistant(int arg1, int arg2)
std::vector< size_type > lines_cache
boost::iterator_range< const_iterator > LineRange
const std::vector< std::pair< size_type, number > > * get_constraint_entries(const size_type line_n) const
bool is_inhomogeneously_constrained(const size_type index) const
unsigned int global_dof_index
Definition: types.h:76
void get_dof_values(const VectorType &global_vector, ForwardIteratorInd local_indices_begin, ForwardIteratorVec local_vector_begin, ForwardIteratorVec local_vector_end) const
number get_inhomogeneity(const size_type line_n) const
ConstraintLine(const size_type &index=numbers::invalid_dof_index, const Entries &entries={}, const number &inhomogeneity=0.0)
void compress() const
Definition: index_set.h:1642
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:584
bool is_element(const size_type index) const
Definition: index_set.h:1765
void add_line(const size_type line_n)
const types::global_dof_index invalid_dof_index
Definition: types.h:211
typename std::vector< ConstraintLine >::const_iterator const_iterator
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2022
static ::ExceptionBase & ExcEntryAlreadyExists(size_type arg1, size_type arg2, number arg3, number arg4)
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
static ::ExceptionBase & ExcInternalError()