Reference documentation for deal.II version Git 040c6ad7d4 2020-09-26 18:01:03 +0200
\(\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 - 2019 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 <utility>
35 #include <vector>
36 
38 
39 // Forward declarations
40 #ifndef DOXYGEN
41 template <typename>
42 class FullMatrix;
43 class SparsityPattern;
47 template <typename number>
48 class SparseMatrix;
49 template <typename number>
50 class BlockSparseMatrix;
51 
52 namespace internal
53 {
54  namespace AffineConstraints
55  {
57 
75  struct Distributing
76  {
77  Distributing(const size_type global_row = numbers::invalid_size_type,
78  const size_type local_row = numbers::invalid_size_type);
79 
80  Distributing(const Distributing &in);
81 
82  Distributing &
83  operator=(const Distributing &in);
84 
85  bool
86  operator<(const Distributing &in) const
87  {
88  return global_row < in.global_row;
89  }
90 
91  size_type global_row;
92  size_type local_row;
93  mutable size_type constraint_position;
94  };
95 
96 
97 
109  template <typename number>
110  struct DataCache
111  {
112  DataCache();
113 
114  void
115  reinit();
116 
117  size_type
118  insert_new_index(const std::pair<size_type, number> &pair);
119 
120  void
121  append_index(const size_type index,
122  const std::pair<size_type, number> &pair);
123 
124  size_type
125  get_size(const size_type index) const;
126 
127  const std::pair<size_type, number> *
128  get_entry(const size_type index) const;
129 
130  size_type row_length;
131 
132  std::vector<std::pair<size_type, number>> data;
133 
134  std::vector<size_type> individual_size;
135  };
136 
137 
138 
167  template <typename number>
168  class GlobalRowsFromLocal
169  {
170  public:
174  GlobalRowsFromLocal();
175 
176  void
177  reinit(const size_type n_local_rows);
178 
179  void
180  insert_index(const size_type global_row,
181  const size_type local_row,
182  const number constraint_value);
183  void
184  sort();
185 
186  void
187  print(std::ostream &os);
188 
192  size_type
193  size() const;
194 
199  size_type
200  size(const size_type counter_index) const;
201 
205  size_type
206  global_row(const size_type counter_index) const;
207 
211  size_type &
212  global_row(const size_type counter_index);
213 
219  size_type
220  local_row(const size_type counter_index) const;
221 
225  size_type &
226  local_row(const size_type counter_index);
227 
233  size_type
234  local_row(const size_type counter_index,
235  const size_type index_in_constraint) const;
236 
241  number
242  constraint_value(const size_type counter_index,
243  const size_type index_in_constraint) const;
244 
250  bool
251  have_indirect_rows() const;
252 
257  void
258  insert_constraint(const size_type constrained_local_dof);
259 
265  size_type
266  n_constraints() const;
267 
272  size_type
273  n_inhomogeneities() const;
274 
281  void
282  set_ith_constraint_inhomogeneous(const size_type i);
283 
288  size_type
289  constraint_origin(size_type i) const;
290 
296  std::vector<Distributing> total_row_indices;
297 
298  private:
302  DataCache<number> data_cache;
303 
308  size_type n_active_rows;
309 
314  size_type n_inhomogeneous_rows;
315  };
316 
317 
318 
336  template <typename number>
337  struct ScratchData
338  {
342  ScratchData()
343  : in_use(false)
344  {}
345 
349  ScratchData(const ScratchData &)
350  : in_use(false)
351  {}
352 
356  bool in_use;
357 
361  std::vector<size_type> columns;
362 
366  std::vector<number> values;
367 
371  std::vector<size_type> block_starts;
372 
376  std::vector<size_type> vector_indices;
377 
381  std::vector<number> vector_values;
382 
386  GlobalRowsFromLocal<number> global_rows;
387 
391  GlobalRowsFromLocal<number> global_columns;
392  };
393  } // namespace AffineConstraints
394 } // namespace internal
395 
396 namespace internal
397 {
398  namespace AffineConstraintsImplementation
399  {
400  template <class VectorType>
401  void
402  set_zero_all(const std::vector<types::global_dof_index> &cm,
403  VectorType & vec);
404 
405  template <class T>
406  void
407  set_zero_all(const std::vector<types::global_dof_index> &cm,
408  ::Vector<T> & vec);
409 
410  template <class T>
411  void
412  set_zero_all(const std::vector<types::global_dof_index> &cm,
413  ::BlockVector<T> & vec);
414  } // namespace AffineConstraintsImplementation
415 } // namespace internal
416 
417 
418 template <typename number>
419 class AffineConstraints;
420 #endif
421 
422 // TODO[WB]: We should have a function of the kind
423 // AffineConstraints::add_constraint (const size_type constrained_dof,
424 // const std::vector<std::pair<size_type, number> > &entries,
425 // const number inhomogeneity = 0);
426 // rather than building up constraints piecemeal through add_line/add_entry
427 // etc. This would also eliminate the possibility of accidentally changing
428 // existing constraints into something pointless, see the discussion on the
429 // mailing list on "Tiny bug in interpolate_boundary_values" in Sept. 2010.
430 
504 template <typename number = double>
506 {
507 public:
512 
519  {
525 
532 
538  right_object_wins
539  };
540 
558  explicit AffineConstraints(const IndexSet &local_constraints = IndexSet());
559 
563  explicit AffineConstraints(const AffineConstraints &affine_constraints);
564 
568  AffineConstraints(AffineConstraints &&affine_constraints) noexcept =
569  default; // NOLINT
570 
581  operator=(const AffineConstraints &) = delete;
582 
587  operator=(AffineConstraints &&affine_constraints) noexcept =
588  default; // NOLINT
589 
596  void
597  copy_from(const AffineConstraints &other);
598 
605  void
606  reinit(const IndexSet &local_constraints = IndexSet());
607 
614  bool
615  can_store_line(const size_type line_n) const;
616 
623  const IndexSet &
624  get_local_lines() const;
625 
645  void
646  add_selected_constraints(const AffineConstraints &constraints_in,
647  const IndexSet & filter);
648 
658  void
659  add_line(const size_type line_n);
660 
673  void
674  add_lines(const std::vector<bool> &lines);
675 
688  void
689  add_lines(const std::set<size_type> &lines);
690 
703  void
704  add_lines(const IndexSet &lines);
705 
717  void
718  add_entry(const size_type line_n, const size_type column, const number value);
719 
725  void
726  add_entries(const size_type line_n,
727  const std::vector<std::pair<size_type, number>> &col_val_pairs);
728 
735  void
736  set_inhomogeneity(const size_type line_n, const number value);
737 
758  void
759  close();
760 
785  void
786  merge(
787  const AffineConstraints & other_constraints,
788  const MergeConflictBehavior merge_conflict_behavior = no_conflicts_allowed,
789  const bool allow_different_local_lines = false);
790 
803  void
804  shift(const size_type offset);
805 
813  void
814  clear();
815 
828  size_type
829  n_constraints() const;
830 
840  bool
841  is_constrained(const size_type line_n) const;
842 
854  bool
855  is_identity_constrained(const size_type line_n) const;
856 
863  bool
864  are_identity_constrained(const size_type line_n_1,
865  const size_type line_n_2) const;
866 
877  size_type
878  max_constraint_indirections() const;
879 
884  bool
885  is_inhomogeneously_constrained(const size_type index) const;
886 
892  bool
893  has_inhomogeneities() const;
894 
899  const std::vector<std::pair<size_type, number>> *
900  get_constraint_entries(const size_type line_n) const;
901 
906  number
907  get_inhomogeneity(const size_type line_n) const;
908 
929  void
930  print(std::ostream &out) const;
931 
944  void
945  write_dot(std::ostream &) const;
946 
951  std::size_t
952  memory_consumption() const;
953 
960  void
961  resolve_indices(std::vector<types::global_dof_index> &indices) const;
962 
988  void
989  condense(SparsityPattern &sparsity) const;
990 
994  void
995  condense(BlockSparsityPattern &sparsity) const;
996 
1001  void
1002  condense(DynamicSparsityPattern &sparsity) const;
1003 
1008  void
1009  condense(BlockDynamicSparsityPattern &sparsity) const;
1010 
1018  void
1019  condense(SparseMatrix<number> &matrix) const;
1020 
1024  void
1025  condense(BlockSparseMatrix<number> &matrix) const;
1026 
1038  template <class VectorType>
1039  void
1040  condense(VectorType &vec) const;
1041 
1048  template <class VectorType>
1049  void
1050  condense(const VectorType &vec_ghosted, VectorType &output) const;
1051 
1064  template <class VectorType>
1065  void
1066  condense(SparseMatrix<number> &matrix, VectorType &vector) const;
1067 
1072  template <class BlockVectorType>
1073  void
1074  condense(BlockSparseMatrix<number> &matrix, BlockVectorType &vector) const;
1075 
1082  template <class VectorType>
1083  void
1084  set_zero(VectorType &vec) const;
1085 
1137  template <class InVector, class OutVector>
1138  void
1139  distribute_local_to_global(const InVector & local_vector,
1140  const std::vector<size_type> &local_dof_indices,
1141  OutVector & global_vector) const;
1142 
1186  template <typename VectorType>
1187  void
1188  distribute_local_to_global(const Vector<number> & local_vector,
1189  const std::vector<size_type> &local_dof_indices,
1190  VectorType & global_vector,
1191  const FullMatrix<number> & local_matrix) const;
1192 
1212  template <typename VectorType>
1213  void
1214  distribute_local_to_global(
1215  const Vector<number> & local_vector,
1216  const std::vector<size_type> &local_dof_indices_row,
1217  const std::vector<size_type> &local_dof_indices_col,
1218  VectorType & global_vector,
1219  const FullMatrix<number> & local_matrix,
1220  bool diagonal = false) const;
1221 
1225  template <class VectorType>
1226  void
1227  distribute_local_to_global(const size_type index,
1228  const number value,
1229  VectorType & global_vector) const;
1230 
1259  template <typename ForwardIteratorVec,
1260  typename ForwardIteratorInd,
1261  class VectorType>
1262  void
1263  distribute_local_to_global(ForwardIteratorVec local_vector_begin,
1264  ForwardIteratorVec local_vector_end,
1265  ForwardIteratorInd local_indices_begin,
1266  VectorType & global_vector) const;
1267 
1315  template <typename MatrixType>
1316  void
1317  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1318  const std::vector<size_type> &local_dof_indices,
1319  MatrixType & global_matrix) const;
1320 
1348  template <typename MatrixType>
1349  void
1350  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1351  const std::vector<size_type> &row_indices,
1352  const std::vector<size_type> &col_indices,
1353  MatrixType & global_matrix) const;
1354 
1371  template <typename MatrixType>
1372  void
1373  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1374  const std::vector<size_type> &row_indices,
1375  const AffineConstraints &column_affine_constraints,
1376  const std::vector<size_type> &column_indices,
1377  MatrixType & global_matrix) const;
1378 
1395  template <typename MatrixType, typename VectorType>
1396  void
1397  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1398  const Vector<number> & local_vector,
1399  const std::vector<size_type> &local_dof_indices,
1400  MatrixType & global_matrix,
1401  VectorType & global_vector,
1402  bool use_inhomogeneities_for_rhs = false) const;
1403 
1457  template <typename SparsityPatternType>
1458  void
1459  add_entries_local_to_global(
1460  const std::vector<size_type> &local_dof_indices,
1461  SparsityPatternType & sparsity_pattern,
1462  const bool keep_constrained_entries = true,
1463  const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1464 
1468  template <typename SparsityPatternType>
1469  void
1470  add_entries_local_to_global(
1471  const std::vector<size_type> &row_indices,
1472  const std::vector<size_type> &col_indices,
1473  SparsityPatternType & sparsity_pattern,
1474  const bool keep_constrained_entries = true,
1475  const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1476 
1496  template <typename ForwardIteratorVec,
1497  typename ForwardIteratorInd,
1498  class VectorType>
1499  void
1500  get_dof_values(const VectorType & global_vector,
1501  ForwardIteratorInd local_indices_begin,
1502  ForwardIteratorVec local_vector_begin,
1503  ForwardIteratorVec local_vector_end) const;
1504 
1526  template <class VectorType>
1527  void
1528  distribute(VectorType &vec) const;
1529 
1538  {
1543  using Entries = std::vector<std::pair<size_type, number>>;
1544 
1551 
1560 
1565 
1573  bool
1574  operator<(const ConstraintLine &) const;
1575 
1581  bool
1582  operator==(const ConstraintLine &) const;
1583 
1588  std::size_t
1589  memory_consumption() const;
1590 
1594  template <class Archive>
1595  void
1596  serialize(Archive &ar, const unsigned int)
1597  {
1598  ar &index &entries &inhomogeneity;
1599  }
1600  };
1601 
1605  using const_iterator = typename std::vector<ConstraintLine>::const_iterator;
1606 
1610  using LineRange = boost::iterator_range<const_iterator>;
1611 
1620  const LineRange
1621  get_lines() const;
1622 
1650  bool
1651  is_consistent_in_parallel(const std::vector<IndexSet> &locally_owned_dofs,
1652  const IndexSet & locally_active_dofs,
1653  const MPI_Comm mpi_communicator,
1654  const bool verbose = false) const;
1655 
1661  DeclException0(ExcMatrixIsClosed);
1667  DeclException0(ExcMatrixNotClosed);
1674  size_type,
1675  << "The specified line " << arg1 << " does not exist.");
1681  DeclException4(ExcEntryAlreadyExists,
1682  size_type,
1683  size_type,
1684  number,
1685  number,
1686  << "The entry for the indices " << arg1 << " and " << arg2
1687  << " already exists, but the values " << arg3 << " (old) and "
1688  << arg4 << " (new) differ "
1689  << "by " << (arg4 - arg3) << ".");
1695  DeclException2(ExcDoFConstrainedToConstrainedDoF,
1696  int,
1697  int,
1698  << "You tried to constrain DoF " << arg1 << " to DoF " << arg2
1699  << ", but that one is also constrained. This is not allowed!");
1705  DeclException1(ExcDoFIsConstrainedFromBothObjects,
1706  size_type,
1707  << "Degree of freedom " << arg1
1708  << " is constrained from both object in a merge operation.");
1714  DeclException1(ExcDoFIsConstrainedToConstrainedDoF,
1715  size_type,
1716  << "In the given argument a degree of freedom is constrained "
1717  << "to another DoF with number " << arg1
1718  << ", which however is constrained by this object. This is not"
1719  << " allowed.");
1725  DeclException1(ExcRowNotStoredHere,
1726  size_type,
1727  << "The index set given to this constraints object indicates "
1728  << "constraints for degree of freedom " << arg1
1729  << " should not be stored by this object, but a constraint "
1730  << "is being added.");
1731 
1737  DeclException2(ExcColumnNotStoredHere,
1738  size_type,
1739  size_type,
1740  << "The index set given to this constraints object indicates "
1741  << "constraints using degree of freedom " << arg2
1742  << " should not be stored by this object, but a constraint "
1743  << "for degree of freedom " << arg1 << " uses it.");
1744 
1750  DeclException2(ExcIncorrectConstraint,
1751  int,
1752  int,
1753  << "While distributing the constraint for DoF " << arg1
1754  << ", it turns out that one of the processors "
1755  << "who own the " << arg2 << " degrees of freedom that x_"
1756  << arg1 << " is constrained against does not know about "
1757  << "the constraint on x_" << arg1
1758  << ". Did you not initialize the AffineConstraints container "
1759  << "with the appropriate locally_relevant set so "
1760  << "that every processor who owns a DoF that constrains "
1761  << "another DoF also knows about this constraint?");
1762 
1763 private:
1775  std::vector<ConstraintLine> lines;
1776 
1809  std::vector<size_type> lines_cache;
1810 
1817 
1821  bool sorted;
1822 
1824  internal::AffineConstraints::ScratchData<number>>
1826 
1831  size_type
1832  calculate_line_index(const size_type line_n) const;
1833 
1838  template <typename MatrixType, typename VectorType>
1839  void
1840  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1841  const Vector<number> & local_vector,
1842  const std::vector<size_type> &local_dof_indices,
1843  MatrixType & global_matrix,
1844  VectorType & global_vector,
1845  bool use_inhomogeneities_for_rhs,
1846  std::integral_constant<bool, false>) const;
1847 
1852  template <typename MatrixType, typename VectorType>
1853  void
1854  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1855  const Vector<number> & local_vector,
1856  const std::vector<size_type> &local_dof_indices,
1857  MatrixType & global_matrix,
1858  VectorType & global_vector,
1859  bool use_inhomogeneities_for_rhs,
1860  std::integral_constant<bool, true>) const;
1861 
1866  template <typename SparsityPatternType>
1867  void
1868  add_entries_local_to_global(const std::vector<size_type> &local_dof_indices,
1869  SparsityPatternType & sparsity_pattern,
1870  const bool keep_constrained_entries,
1871  const Table<2, bool> &dof_mask,
1872  std::integral_constant<bool, false>) const;
1873 
1878  template <typename SparsityPatternType>
1879  void
1880  add_entries_local_to_global(const std::vector<size_type> &local_dof_indices,
1881  SparsityPatternType & sparsity_pattern,
1882  const bool keep_constrained_entries,
1883  const Table<2, bool> &dof_mask,
1884  std::integral_constant<bool, true>) const;
1885 
1893  void
1894  make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
1895  internal::AffineConstraints::GlobalRowsFromLocal<number>
1896  &global_rows) const;
1897 
1905  void
1906  make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
1907  std::vector<size_type> & active_dofs) const;
1908 
1912  template <typename MatrixScalar, typename VectorScalar>
1914  resolve_vector_entry(
1915  const size_type i,
1916  const internal::AffineConstraints::GlobalRowsFromLocal<number> &global_rows,
1917  const Vector<VectorScalar> & local_vector,
1918  const std::vector<size_type> & local_dof_indices,
1919  const FullMatrix<MatrixScalar> &local_matrix) const;
1920 };
1921 
1922 /* ---------------- template and inline functions ----------------- */
1923 
1924 template <typename number>
1926  const IndexSet &local_constraints)
1927  : lines()
1928  , local_lines(local_constraints)
1929  , sorted(false)
1930 {
1931  // make sure the IndexSet is compressed. Otherwise this can lead to crashes
1932  // that are hard to find (only happen in release mode).
1933  // see tests/mpi/affine_constraints_crash_01
1935 }
1936 
1937 template <typename number>
1939  const AffineConstraints &affine_constraints)
1940  : Subscriptor()
1941  , lines(affine_constraints.lines)
1942  , lines_cache(affine_constraints.lines_cache)
1943  , local_lines(affine_constraints.local_lines)
1944  , sorted(affine_constraints.sorted)
1945 {}
1946 
1947 template <typename number>
1948 inline void
1950 {
1951  Assert(sorted == false, ExcMatrixIsClosed());
1952 
1953  // the following can happen when we compute with distributed meshes and dof
1954  // handlers and we constrain a degree of freedom whose number we don't have
1955  // locally. if we don't abort here the program will try to allocate several
1956  // terabytes of memory to resize the various arrays below :-)
1958  const size_type line_index = calculate_line_index(line_n);
1959 
1960  // check whether line already exists; it may, in which case we can just quit
1961  if (is_constrained(line_n))
1962  return;
1963 
1964  // if necessary enlarge vector of existing entries for cache
1965  if (line_index >= lines_cache.size())
1966  lines_cache.resize(std::max(2 * static_cast<size_type>(lines_cache.size()),
1967  line_index + 1),
1969 
1970  // push a new line to the end of the list
1971  lines.emplace_back();
1972  lines.back().index = line_n;
1973  lines.back().inhomogeneity = 0.;
1974  lines_cache[line_index] = lines.size() - 1;
1975 }
1976 
1977 template <typename number>
1978 inline void
1980  const size_type column,
1981  const number value)
1982 {
1983  Assert(sorted == false, ExcMatrixIsClosed());
1984  Assert(line_n != column,
1985  ExcMessage("Can't constrain a degree of freedom to itself"));
1986 
1987  // Ensure that the current line is present in the cache:
1988  const size_type line_index = calculate_line_index(line_n);
1989  Assert(line_index < lines_cache.size(),
1990  ExcMessage("The current AffineConstraints does not contain the line "
1991  "for the current entry. Call AffineConstraints::add_line "
1992  "before calling this function."));
1993 
1994  // if in debug mode, check whether an entry for this column already exists
1995  // and if it's the same as the one entered at present
1996  //
1997  // in any case: exit the function if an entry for this column already
1998  // exists, since we don't want to enter it twice
2000  ExcInternalError());
2002  ExcColumnNotStoredHere(line_n, column));
2003  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
2004  Assert(line_ptr->index == line_n, ExcInternalError());
2005  for (const auto &p : line_ptr->entries)
2006  if (p.first == column)
2007  {
2008  Assert(std::abs(p.second - value) < 1.e-14,
2009  ExcEntryAlreadyExists(line_n, column, p.second, value));
2010  return;
2011  }
2012 
2013  line_ptr->entries.emplace_back(column, value);
2014 }
2015 
2016 template <typename number>
2017 inline void
2019  const number value)
2020 {
2021  const size_type line_index = calculate_line_index(line_n);
2022  Assert(line_index < lines_cache.size() &&
2023  lines_cache[line_index] != numbers::invalid_size_type,
2024  ExcMessage("call add_line() before calling set_inhomogeneity()"));
2025  Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
2026  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
2027  line_ptr->inhomogeneity = value;
2028 }
2029 
2030 template <typename number>
2031 template <class VectorType>
2032 inline void
2034 {
2035  // since lines is a private member, we cannot pass it to the functions
2036  // above. therefore, copy the content which is cheap
2037  std::vector<size_type> constrained_lines(lines.size());
2038  for (unsigned int i = 0; i < lines.size(); ++i)
2039  constrained_lines[i] = lines[i].index;
2040  internal::AffineConstraintsImplementation::set_zero_all(constrained_lines,
2041  vec);
2042 }
2043 
2044 template <typename number>
2047 {
2048  return lines.size();
2049 }
2050 
2051 template <typename number>
2052 inline bool
2054 {
2055  const size_type line_index = calculate_line_index(index);
2056  return ((line_index < lines_cache.size()) &&
2057  (lines_cache[line_index] != numbers::invalid_size_type));
2058 }
2059 
2060 template <typename number>
2061 inline bool
2063  const size_type line_n) const
2064 {
2065  // check whether the entry is constrained. could use is_constrained, but
2066  // that means computing the line index twice
2067  const size_type line_index = calculate_line_index(line_n);
2068  if (line_index >= lines_cache.size() ||
2069  lines_cache[line_index] == numbers::invalid_size_type)
2070  return false;
2071  else
2072  {
2073  Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
2074  return !(lines[lines_cache[line_index]].inhomogeneity == number(0.));
2075  }
2076 }
2077 
2078 template <typename number>
2079 inline const std::vector<std::pair<types::global_dof_index, number>> *
2081 {
2082  // check whether the entry is constrained. could use is_constrained, but
2083  // that means computing the line index twice
2084  const size_type line_index = calculate_line_index(line_n);
2085  if (line_index >= lines_cache.size() ||
2086  lines_cache[line_index] == numbers::invalid_size_type)
2087  return nullptr;
2088  else
2089  return &lines[lines_cache[line_index]].entries;
2090 }
2091 
2092 template <typename number>
2093 inline number
2095 {
2096  // check whether the entry is constrained. could use is_constrained, but
2097  // that means computing the line index twice
2098  const size_type line_index = calculate_line_index(line_n);
2099  if (line_index >= lines_cache.size() ||
2100  lines_cache[line_index] == numbers::invalid_size_type)
2101  return 0;
2102  else
2103  return lines[lines_cache[line_index]].inhomogeneity;
2104 }
2105 
2106 template <typename number>
2109 {
2110  // IndexSet is unused (serial case)
2111  if (!local_lines.size())
2112  return line_n;
2113 
2115 
2116  return local_lines.index_within_set(line_n);
2117 }
2118 
2119 template <typename number>
2120 inline bool
2122 {
2123  return !local_lines.size() || local_lines.is_element(line_n);
2124 }
2125 
2126 template <typename number>
2127 inline const IndexSet &
2129 {
2130  return local_lines;
2131 }
2132 
2133 template <typename number>
2134 template <class VectorType>
2135 inline void
2137  const size_type index,
2138  const number value,
2139  VectorType & global_vector) const
2140 {
2141  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2142 
2143  if (is_constrained(index) == false)
2144  global_vector(index) += value;
2145  else
2146  {
2147  const ConstraintLine &position =
2149  for (size_type j = 0; j < position.entries.size(); ++j)
2150  global_vector(position.entries[j].first) +=
2151  value * position.entries[j].second;
2152  }
2153 }
2154 
2155 template <typename number>
2156 template <typename ForwardIteratorVec,
2157  typename ForwardIteratorInd,
2158  class VectorType>
2159 inline void
2161  ForwardIteratorVec local_vector_begin,
2162  ForwardIteratorVec local_vector_end,
2163  ForwardIteratorInd local_indices_begin,
2164  VectorType & global_vector) const
2165 {
2166  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2167  for (; local_vector_begin != local_vector_end;
2168  ++local_vector_begin, ++local_indices_begin)
2169  {
2170  if (is_constrained(*local_indices_begin) == false)
2171  internal::ElementAccess<VectorType>::add(*local_vector_begin,
2172  *local_indices_begin,
2173  global_vector);
2174  else
2175  {
2176  const ConstraintLine &position =
2177  lines[lines_cache[calculate_line_index(*local_indices_begin)]];
2178  for (size_type j = 0; j < position.entries.size(); ++j)
2180  (*local_vector_begin) * position.entries[j].second,
2181  position.entries[j].first,
2182  global_vector);
2183  }
2184  }
2185 }
2186 
2187 template <typename number>
2188 template <class InVector, class OutVector>
2189 inline void
2191  const InVector & local_vector,
2192  const std::vector<size_type> &local_dof_indices,
2193  OutVector & global_vector) const
2194 {
2195  Assert(local_vector.size() == local_dof_indices.size(),
2196  ExcDimensionMismatch(local_vector.size(), local_dof_indices.size()));
2197  distribute_local_to_global(local_vector.begin(),
2198  local_vector.end(),
2199  local_dof_indices.begin(),
2200  global_vector);
2201 }
2202 
2203 template <typename number>
2204 template <typename ForwardIteratorVec,
2205  typename ForwardIteratorInd,
2206  class VectorType>
2207 inline void
2209  const VectorType & global_vector,
2210  ForwardIteratorInd local_indices_begin,
2211  ForwardIteratorVec local_vector_begin,
2212  ForwardIteratorVec local_vector_end) const
2213 {
2214  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2215  for (; local_vector_begin != local_vector_end;
2216  ++local_vector_begin, ++local_indices_begin)
2217  {
2218  if (is_constrained(*local_indices_begin) == false)
2219  *local_vector_begin = global_vector(*local_indices_begin);
2220  else
2221  {
2222  const ConstraintLine &position =
2223  lines[lines_cache[calculate_line_index(*local_indices_begin)]];
2224  typename VectorType::value_type value = position.inhomogeneity;
2225  for (size_type j = 0; j < position.entries.size(); ++j)
2226  value += (global_vector(position.entries[j].first) *
2227  position.entries[j].second);
2228  *local_vector_begin = value;
2229  }
2230  }
2231 }
2232 
2233 template <typename MatrixType>
2235 template <typename SparsityPatternType>
2237 template <typename number>
2239 
2257 template <typename MatrixType>
2259 {
2260 private:
2261  struct yes_type
2262  {
2263  char c[1];
2264  };
2265  struct no_type
2266  {
2267  char c[2];
2268  };
2269 
2275  template <typename T>
2276  static yes_type
2277  check_for_block_matrix(const BlockMatrixBase<T> *);
2278 
2283  template <typename T>
2284  static yes_type
2285  check_for_block_matrix(const BlockSparsityPatternBase<T> *);
2286 
2291  template <typename T>
2292  static yes_type
2293  check_for_block_matrix(const BlockSparseMatrixEZ<T> *);
2294 
2299  static no_type
2300  check_for_block_matrix(...);
2301 
2302 public:
2308  static const bool value =
2309  (sizeof(check_for_block_matrix(static_cast<MatrixType *>(nullptr))) ==
2310  sizeof(yes_type));
2311 };
2312 
2313 // instantiation of the static member
2314 template <typename MatrixType>
2316 
2317 template <typename number>
2318 template <typename MatrixType>
2319 inline void
2321  const FullMatrix<number> & local_matrix,
2322  const std::vector<size_type> &local_dof_indices,
2323  MatrixType & global_matrix) const
2324 {
2325  // create a dummy and hand on to the function actually implementing this
2326  // feature in the cm.templates.h file.
2327  Vector<typename MatrixType::value_type> dummy(0);
2329  local_matrix,
2330  dummy,
2331  local_dof_indices,
2332  global_matrix,
2333  dummy,
2334  false,
2335  std::integral_constant<bool, IsBlockMatrix<MatrixType>::value>());
2336 }
2337 
2338 template <typename number>
2339 template <typename MatrixType, typename VectorType>
2340 inline void
2342  const FullMatrix<number> & local_matrix,
2343  const Vector<number> & local_vector,
2344  const std::vector<size_type> &local_dof_indices,
2345  MatrixType & global_matrix,
2346  VectorType & global_vector,
2347  bool use_inhomogeneities_for_rhs) const
2348 {
2349  // enter the internal function with the respective block information set,
2350  // the actual implementation follows in the cm.templates.h file.
2352  local_matrix,
2353  local_vector,
2354  local_dof_indices,
2355  global_matrix,
2356  global_vector,
2357  use_inhomogeneities_for_rhs,
2358  std::integral_constant<bool, IsBlockMatrix<MatrixType>::value>());
2359 }
2360 
2361 template <typename number>
2362 template <typename SparsityPatternType>
2363 inline void
2365  const std::vector<size_type> &local_dof_indices,
2366  SparsityPatternType & sparsity_pattern,
2367  const bool keep_constrained_entries,
2368  const Table<2, bool> & dof_mask) const
2369 {
2370  // enter the internal function with the respective block information set,
2371  // the actual implementation follows in the cm.templates.h file.
2373  local_dof_indices,
2374  sparsity_pattern,
2375  keep_constrained_entries,
2376  dof_mask,
2377  std::integral_constant<bool, IsBlockMatrix<SparsityPatternType>::value>());
2378 }
2379 
2381 
2382 #endif
const types::global_dof_index invalid_size_type
Definition: types.h:205
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 add_entry(const size_type line_n, const size_type column, const number value)
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:1632
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:1411
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1919
#define DeclException0(Exception0)
Definition: exceptions.h:470
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
constexpr TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
AffineConstraints(const IndexSet &local_constraints=IndexSet())
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
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
void compress() const
Definition: index_set.h:1640
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:584
Definition: table.h:687
bool is_element(const size_type index) const
Definition: index_set.h:1763
void add_line(const size_type line_n)
void set_inhomogeneity(const size_type line_n, const number value)
typename std::vector< ConstraintLine >::const_iterator const_iterator
T max(const T &t, const MPI_Comm &mpi_communicator)
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:815
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()