Reference documentation for deal.II version Git 01d13113a1 2020-07-10 09:14:01 +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) = default; // NOLINT
569 
580  operator=(const AffineConstraints &) = delete;
581 
586  operator=(AffineConstraints &&affine_constraints) = default; // NOLINT
587 
594  void
595  copy_from(const AffineConstraints &other);
596 
603  void
604  reinit(const IndexSet &local_constraints = IndexSet());
605 
612  bool
613  can_store_line(const size_type line_n) const;
614 
621  const IndexSet &
622  get_local_lines() const;
623 
643  void
644  add_selected_constraints(const AffineConstraints &constraints_in,
645  const IndexSet & filter);
646 
656  void
657  add_line(const size_type line_n);
658 
671  void
672  add_lines(const std::vector<bool> &lines);
673 
686  void
687  add_lines(const std::set<size_type> &lines);
688 
701  void
702  add_lines(const IndexSet &lines);
703 
715  void
716  add_entry(const size_type line_n, const size_type column, const number value);
717 
723  void
724  add_entries(const size_type line_n,
725  const std::vector<std::pair<size_type, number>> &col_val_pairs);
726 
733  void
734  set_inhomogeneity(const size_type line_n, const number value);
735 
756  void
757  close();
758 
783  void
784  merge(
785  const AffineConstraints & other_constraints,
786  const MergeConflictBehavior merge_conflict_behavior = no_conflicts_allowed,
787  const bool allow_different_local_lines = false);
788 
801  void
802  shift(const size_type offset);
803 
811  void
812  clear();
813 
826  size_type
827  n_constraints() const;
828 
838  bool
839  is_constrained(const size_type line_n) const;
840 
852  bool
853  is_identity_constrained(const size_type line_n) const;
854 
861  bool
862  are_identity_constrained(const size_type line_n_1,
863  const size_type line_n_2) const;
864 
875  size_type
876  max_constraint_indirections() const;
877 
882  bool
883  is_inhomogeneously_constrained(const size_type index) const;
884 
890  bool
891  has_inhomogeneities() const;
892 
897  const std::vector<std::pair<size_type, number>> *
898  get_constraint_entries(const size_type line_n) const;
899 
904  number
905  get_inhomogeneity(const size_type line_n) const;
906 
927  void
928  print(std::ostream &out) const;
929 
942  void
943  write_dot(std::ostream &) const;
944 
949  std::size_t
950  memory_consumption() const;
951 
958  void
959  resolve_indices(std::vector<types::global_dof_index> &indices) const;
960 
986  void
987  condense(SparsityPattern &sparsity) const;
988 
992  void
993  condense(BlockSparsityPattern &sparsity) const;
994 
999  void
1000  condense(DynamicSparsityPattern &sparsity) const;
1001 
1006  void
1007  condense(BlockDynamicSparsityPattern &sparsity) const;
1008 
1016  void
1017  condense(SparseMatrix<number> &matrix) const;
1018 
1022  void
1023  condense(BlockSparseMatrix<number> &matrix) const;
1024 
1036  template <class VectorType>
1037  void
1038  condense(VectorType &vec) const;
1039 
1046  template <class VectorType>
1047  void
1048  condense(const VectorType &vec_ghosted, VectorType &output) const;
1049 
1062  template <class VectorType>
1063  void
1064  condense(SparseMatrix<number> &matrix, VectorType &vector) const;
1065 
1070  template <class BlockVectorType>
1071  void
1072  condense(BlockSparseMatrix<number> &matrix, BlockVectorType &vector) const;
1073 
1080  template <class VectorType>
1081  void
1082  set_zero(VectorType &vec) const;
1083 
1135  template <class InVector, class OutVector>
1136  void
1137  distribute_local_to_global(const InVector & local_vector,
1138  const std::vector<size_type> &local_dof_indices,
1139  OutVector & global_vector) const;
1140 
1184  template <typename VectorType>
1185  void
1186  distribute_local_to_global(const Vector<number> & local_vector,
1187  const std::vector<size_type> &local_dof_indices,
1188  VectorType & global_vector,
1189  const FullMatrix<number> & local_matrix) const;
1190 
1210  template <typename VectorType>
1211  void
1212  distribute_local_to_global(
1213  const Vector<number> & local_vector,
1214  const std::vector<size_type> &local_dof_indices_row,
1215  const std::vector<size_type> &local_dof_indices_col,
1216  VectorType & global_vector,
1217  const FullMatrix<number> & local_matrix,
1218  bool diagonal = false) const;
1219 
1223  template <class VectorType>
1224  void
1225  distribute_local_to_global(const size_type index,
1226  const number value,
1227  VectorType & global_vector) const;
1228 
1257  template <typename ForwardIteratorVec,
1258  typename ForwardIteratorInd,
1259  class VectorType>
1260  void
1261  distribute_local_to_global(ForwardIteratorVec local_vector_begin,
1262  ForwardIteratorVec local_vector_end,
1263  ForwardIteratorInd local_indices_begin,
1264  VectorType & global_vector) const;
1265 
1313  template <typename MatrixType>
1314  void
1315  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1316  const std::vector<size_type> &local_dof_indices,
1317  MatrixType & global_matrix) const;
1318 
1346  template <typename MatrixType>
1347  void
1348  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1349  const std::vector<size_type> &row_indices,
1350  const std::vector<size_type> &col_indices,
1351  MatrixType & global_matrix) const;
1352 
1369  template <typename MatrixType>
1370  void
1371  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1372  const std::vector<size_type> &row_indices,
1373  const AffineConstraints &column_affine_constraints,
1374  const std::vector<size_type> &column_indices,
1375  MatrixType & global_matrix) const;
1376 
1393  template <typename MatrixType, typename VectorType>
1394  void
1395  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1396  const Vector<number> & local_vector,
1397  const std::vector<size_type> &local_dof_indices,
1398  MatrixType & global_matrix,
1399  VectorType & global_vector,
1400  bool use_inhomogeneities_for_rhs = false) const;
1401 
1455  template <typename SparsityPatternType>
1456  void
1457  add_entries_local_to_global(
1458  const std::vector<size_type> &local_dof_indices,
1459  SparsityPatternType & sparsity_pattern,
1460  const bool keep_constrained_entries = true,
1461  const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1462 
1466  template <typename SparsityPatternType>
1467  void
1468  add_entries_local_to_global(
1469  const std::vector<size_type> &row_indices,
1470  const std::vector<size_type> &col_indices,
1471  SparsityPatternType & sparsity_pattern,
1472  const bool keep_constrained_entries = true,
1473  const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1474 
1494  template <typename ForwardIteratorVec,
1495  typename ForwardIteratorInd,
1496  class VectorType>
1497  void
1498  get_dof_values(const VectorType & global_vector,
1499  ForwardIteratorInd local_indices_begin,
1500  ForwardIteratorVec local_vector_begin,
1501  ForwardIteratorVec local_vector_end) const;
1502 
1524  template <class VectorType>
1525  void
1526  distribute(VectorType &vec) const;
1527 
1536  {
1541  using Entries = std::vector<std::pair<size_type, number>>;
1542 
1549 
1558 
1563 
1571  bool
1572  operator<(const ConstraintLine &) const;
1573 
1579  bool
1580  operator==(const ConstraintLine &) const;
1581 
1586  std::size_t
1587  memory_consumption() const;
1588 
1592  template <class Archive>
1593  void
1594  serialize(Archive &ar, const unsigned int)
1595  {
1596  ar &index &entries &inhomogeneity;
1597  }
1598  };
1599 
1603  using const_iterator = typename std::vector<ConstraintLine>::const_iterator;
1604 
1608  using LineRange = boost::iterator_range<const_iterator>;
1609 
1618  const LineRange
1619  get_lines() const;
1620 
1648  bool
1649  is_consistent_in_parallel(const std::vector<IndexSet> &locally_owned_dofs,
1650  const IndexSet & locally_active_dofs,
1651  const MPI_Comm mpi_communicator,
1652  const bool verbose = false) const;
1653 
1659  DeclException0(ExcMatrixIsClosed);
1665  DeclException0(ExcMatrixNotClosed);
1672  size_type,
1673  << "The specified line " << arg1 << " does not exist.");
1679  DeclException4(ExcEntryAlreadyExists,
1680  size_type,
1681  size_type,
1682  number,
1683  number,
1684  << "The entry for the indices " << arg1 << " and " << arg2
1685  << " already exists, but the values " << arg3 << " (old) and "
1686  << arg4 << " (new) differ "
1687  << "by " << (arg4 - arg3) << ".");
1693  DeclException2(ExcDoFConstrainedToConstrainedDoF,
1694  int,
1695  int,
1696  << "You tried to constrain DoF " << arg1 << " to DoF " << arg2
1697  << ", but that one is also constrained. This is not allowed!");
1703  DeclException1(ExcDoFIsConstrainedFromBothObjects,
1704  size_type,
1705  << "Degree of freedom " << arg1
1706  << " is constrained from both object in a merge operation.");
1712  DeclException1(ExcDoFIsConstrainedToConstrainedDoF,
1713  size_type,
1714  << "In the given argument a degree of freedom is constrained "
1715  << "to another DoF with number " << arg1
1716  << ", which however is constrained by this object. This is not"
1717  << " allowed.");
1723  DeclException1(ExcRowNotStoredHere,
1724  size_type,
1725  << "The index set given to this constraints object indicates "
1726  << "constraints for degree of freedom " << arg1
1727  << " should not be stored by this object, but a constraint "
1728  << "is being added.");
1729 
1735  DeclException2(ExcColumnNotStoredHere,
1736  size_type,
1737  size_type,
1738  << "The index set given to this constraints object indicates "
1739  << "constraints using degree of freedom " << arg2
1740  << " should not be stored by this object, but a constraint "
1741  << "for degree of freedom " << arg1 << " uses it.");
1742 
1748  DeclException2(ExcIncorrectConstraint,
1749  int,
1750  int,
1751  << "While distributing the constraint for DoF " << arg1
1752  << ", it turns out that one of the processors "
1753  << "who own the " << arg2 << " degrees of freedom that x_"
1754  << arg1 << " is constrained against does not know about "
1755  << "the constraint on x_" << arg1
1756  << ". Did you not initialize the AffineConstraints container "
1757  << "with the appropriate locally_relevant set so "
1758  << "that every processor who owns a DoF that constrains "
1759  << "another DoF also knows about this constraint?");
1760 
1761 private:
1773  std::vector<ConstraintLine> lines;
1774 
1807  std::vector<size_type> lines_cache;
1808 
1815 
1819  bool sorted;
1820 
1822  internal::AffineConstraints::ScratchData<number>>
1824 
1829  size_type
1830  calculate_line_index(const size_type line_n) const;
1831 
1836  template <typename MatrixType, typename VectorType>
1837  void
1838  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1839  const Vector<number> & local_vector,
1840  const std::vector<size_type> &local_dof_indices,
1841  MatrixType & global_matrix,
1842  VectorType & global_vector,
1843  bool use_inhomogeneities_for_rhs,
1844  std::integral_constant<bool, false>) const;
1845 
1850  template <typename MatrixType, typename VectorType>
1851  void
1852  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1853  const Vector<number> & local_vector,
1854  const std::vector<size_type> &local_dof_indices,
1855  MatrixType & global_matrix,
1856  VectorType & global_vector,
1857  bool use_inhomogeneities_for_rhs,
1858  std::integral_constant<bool, true>) const;
1859 
1864  template <typename SparsityPatternType>
1865  void
1866  add_entries_local_to_global(const std::vector<size_type> &local_dof_indices,
1867  SparsityPatternType & sparsity_pattern,
1868  const bool keep_constrained_entries,
1869  const Table<2, bool> &dof_mask,
1870  std::integral_constant<bool, false>) const;
1871 
1876  template <typename SparsityPatternType>
1877  void
1878  add_entries_local_to_global(const std::vector<size_type> &local_dof_indices,
1879  SparsityPatternType & sparsity_pattern,
1880  const bool keep_constrained_entries,
1881  const Table<2, bool> &dof_mask,
1882  std::integral_constant<bool, true>) const;
1883 
1891  void
1892  make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
1893  internal::AffineConstraints::GlobalRowsFromLocal<number>
1894  &global_rows) const;
1895 
1903  void
1904  make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
1905  std::vector<size_type> & active_dofs) const;
1906 
1910  template <typename MatrixScalar, typename VectorScalar>
1912  resolve_vector_entry(
1913  const size_type i,
1914  const internal::AffineConstraints::GlobalRowsFromLocal<number> &global_rows,
1915  const Vector<VectorScalar> & local_vector,
1916  const std::vector<size_type> & local_dof_indices,
1917  const FullMatrix<MatrixScalar> &local_matrix) const;
1918 };
1919 
1920 /* ---------------- template and inline functions ----------------- */
1921 
1922 template <typename number>
1924  const IndexSet &local_constraints)
1925  : lines()
1926  , local_lines(local_constraints)
1927  , sorted(false)
1928 {
1929  // make sure the IndexSet is compressed. Otherwise this can lead to crashes
1930  // that are hard to find (only happen in release mode).
1931  // see tests/mpi/affine_constraints_crash_01
1933 }
1934 
1935 template <typename number>
1937  const AffineConstraints &affine_constraints)
1938  : Subscriptor()
1939  , lines(affine_constraints.lines)
1940  , lines_cache(affine_constraints.lines_cache)
1941  , local_lines(affine_constraints.local_lines)
1942  , sorted(affine_constraints.sorted)
1943 {}
1944 
1945 template <typename number>
1946 inline void
1948 {
1949  Assert(sorted == false, ExcMatrixIsClosed());
1950 
1951  // the following can happen when we compute with distributed meshes and dof
1952  // handlers and we constrain a degree of freedom whose number we don't have
1953  // locally. if we don't abort here the program will try to allocate several
1954  // terabytes of memory to resize the various arrays below :-)
1956  const size_type line_index = calculate_line_index(line_n);
1957 
1958  // check whether line already exists; it may, in which case we can just quit
1959  if (is_constrained(line_n))
1960  return;
1961 
1962  // if necessary enlarge vector of existing entries for cache
1963  if (line_index >= lines_cache.size())
1964  lines_cache.resize(std::max(2 * static_cast<size_type>(lines_cache.size()),
1965  line_index + 1),
1967 
1968  // push a new line to the end of the list
1969  lines.emplace_back();
1970  lines.back().index = line_n;
1971  lines.back().inhomogeneity = 0.;
1972  lines_cache[line_index] = lines.size() - 1;
1973 }
1974 
1975 template <typename number>
1976 inline void
1978  const size_type column,
1979  const number value)
1980 {
1981  Assert(sorted == false, ExcMatrixIsClosed());
1982  Assert(line_n != column,
1983  ExcMessage("Can't constrain a degree of freedom to itself"));
1984 
1985  // Ensure that the current line is present in the cache:
1986  const size_type line_index = calculate_line_index(line_n);
1987  Assert(line_index < lines_cache.size(),
1988  ExcMessage("The current AffineConstraints does not contain the line "
1989  "for the current entry. Call AffineConstraints::add_line "
1990  "before calling this function."));
1991 
1992  // if in debug mode, check whether an entry for this column already exists
1993  // and if it's the same as the one entered at present
1994  //
1995  // in any case: exit the function if an entry for this column already
1996  // exists, since we don't want to enter it twice
1998  ExcInternalError());
2000  ExcColumnNotStoredHere(line_n, column));
2001  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
2002  Assert(line_ptr->index == line_n, ExcInternalError());
2003  for (const auto &p : line_ptr->entries)
2004  if (p.first == column)
2005  {
2006  Assert(std::abs(p.second - value) < 1.e-14,
2007  ExcEntryAlreadyExists(line_n, column, p.second, value));
2008  return;
2009  }
2010 
2011  line_ptr->entries.emplace_back(column, value);
2012 }
2013 
2014 template <typename number>
2015 inline void
2017  const number value)
2018 {
2019  const size_type line_index = calculate_line_index(line_n);
2020  Assert(line_index < lines_cache.size() &&
2021  lines_cache[line_index] != numbers::invalid_size_type,
2022  ExcMessage("call add_line() before calling set_inhomogeneity()"));
2023  Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
2024  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
2025  line_ptr->inhomogeneity = value;
2026 }
2027 
2028 template <typename number>
2029 template <class VectorType>
2030 inline void
2032 {
2033  // since lines is a private member, we cannot pass it to the functions
2034  // above. therefore, copy the content which is cheap
2035  std::vector<size_type> constrained_lines(lines.size());
2036  for (unsigned int i = 0; i < lines.size(); ++i)
2037  constrained_lines[i] = lines[i].index;
2038  internal::AffineConstraintsImplementation::set_zero_all(constrained_lines,
2039  vec);
2040 }
2041 
2042 template <typename number>
2045 {
2046  return lines.size();
2047 }
2048 
2049 template <typename number>
2050 inline bool
2052 {
2053  const size_type line_index = calculate_line_index(index);
2054  return ((line_index < lines_cache.size()) &&
2055  (lines_cache[line_index] != numbers::invalid_size_type));
2056 }
2057 
2058 template <typename number>
2059 inline bool
2061  const size_type line_n) const
2062 {
2063  // check whether the entry is constrained. could use is_constrained, but
2064  // that means computing the line index twice
2065  const size_type line_index = calculate_line_index(line_n);
2066  if (line_index >= lines_cache.size() ||
2067  lines_cache[line_index] == numbers::invalid_size_type)
2068  return false;
2069  else
2070  {
2071  Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
2072  return !(lines[lines_cache[line_index]].inhomogeneity == number(0.));
2073  }
2074 }
2075 
2076 template <typename number>
2077 inline const std::vector<std::pair<types::global_dof_index, number>> *
2079 {
2080  // check whether the entry is constrained. could use is_constrained, but
2081  // that means computing the line index twice
2082  const size_type line_index = calculate_line_index(line_n);
2083  if (line_index >= lines_cache.size() ||
2084  lines_cache[line_index] == numbers::invalid_size_type)
2085  return nullptr;
2086  else
2087  return &lines[lines_cache[line_index]].entries;
2088 }
2089 
2090 template <typename number>
2091 inline number
2093 {
2094  // check whether the entry is constrained. could use is_constrained, but
2095  // that means computing the line index twice
2096  const size_type line_index = calculate_line_index(line_n);
2097  if (line_index >= lines_cache.size() ||
2098  lines_cache[line_index] == numbers::invalid_size_type)
2099  return 0;
2100  else
2101  return lines[lines_cache[line_index]].inhomogeneity;
2102 }
2103 
2104 template <typename number>
2107 {
2108  // IndexSet is unused (serial case)
2109  if (!local_lines.size())
2110  return line_n;
2111 
2113 
2114  return local_lines.index_within_set(line_n);
2115 }
2116 
2117 template <typename number>
2118 inline bool
2120 {
2121  return !local_lines.size() || local_lines.is_element(line_n);
2122 }
2123 
2124 template <typename number>
2125 inline const IndexSet &
2127 {
2128  return local_lines;
2129 }
2130 
2131 template <typename number>
2132 template <class VectorType>
2133 inline void
2135  const size_type index,
2136  const number value,
2137  VectorType & global_vector) const
2138 {
2139  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2140 
2141  if (is_constrained(index) == false)
2142  global_vector(index) += value;
2143  else
2144  {
2145  const ConstraintLine &position =
2147  for (size_type j = 0; j < position.entries.size(); ++j)
2148  global_vector(position.entries[j].first) +=
2149  value * position.entries[j].second;
2150  }
2151 }
2152 
2153 template <typename number>
2154 template <typename ForwardIteratorVec,
2155  typename ForwardIteratorInd,
2156  class VectorType>
2157 inline void
2159  ForwardIteratorVec local_vector_begin,
2160  ForwardIteratorVec local_vector_end,
2161  ForwardIteratorInd local_indices_begin,
2162  VectorType & global_vector) const
2163 {
2164  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2165  for (; local_vector_begin != local_vector_end;
2166  ++local_vector_begin, ++local_indices_begin)
2167  {
2168  if (is_constrained(*local_indices_begin) == false)
2169  internal::ElementAccess<VectorType>::add(*local_vector_begin,
2170  *local_indices_begin,
2171  global_vector);
2172  else
2173  {
2174  const ConstraintLine &position =
2175  lines[lines_cache[calculate_line_index(*local_indices_begin)]];
2176  for (size_type j = 0; j < position.entries.size(); ++j)
2178  (*local_vector_begin) * position.entries[j].second,
2179  position.entries[j].first,
2180  global_vector);
2181  }
2182  }
2183 }
2184 
2185 template <typename number>
2186 template <class InVector, class OutVector>
2187 inline void
2189  const InVector & local_vector,
2190  const std::vector<size_type> &local_dof_indices,
2191  OutVector & global_vector) const
2192 {
2193  Assert(local_vector.size() == local_dof_indices.size(),
2194  ExcDimensionMismatch(local_vector.size(), local_dof_indices.size()));
2195  distribute_local_to_global(local_vector.begin(),
2196  local_vector.end(),
2197  local_dof_indices.begin(),
2198  global_vector);
2199 }
2200 
2201 template <typename number>
2202 template <typename ForwardIteratorVec,
2203  typename ForwardIteratorInd,
2204  class VectorType>
2205 inline void
2207  const VectorType & global_vector,
2208  ForwardIteratorInd local_indices_begin,
2209  ForwardIteratorVec local_vector_begin,
2210  ForwardIteratorVec local_vector_end) const
2211 {
2212  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2213  for (; local_vector_begin != local_vector_end;
2214  ++local_vector_begin, ++local_indices_begin)
2215  {
2216  if (is_constrained(*local_indices_begin) == false)
2217  *local_vector_begin = global_vector(*local_indices_begin);
2218  else
2219  {
2220  const ConstraintLine &position =
2221  lines[lines_cache[calculate_line_index(*local_indices_begin)]];
2222  typename VectorType::value_type value = position.inhomogeneity;
2223  for (size_type j = 0; j < position.entries.size(); ++j)
2224  value += (global_vector(position.entries[j].first) *
2225  position.entries[j].second);
2226  *local_vector_begin = value;
2227  }
2228  }
2229 }
2230 
2231 template <typename MatrixType>
2233 template <typename SparsityPatternType>
2235 template <typename number>
2237 
2255 template <typename MatrixType>
2257 {
2258 private:
2259  struct yes_type
2260  {
2261  char c[1];
2262  };
2263  struct no_type
2264  {
2265  char c[2];
2266  };
2267 
2273  template <typename T>
2274  static yes_type
2275  check_for_block_matrix(const BlockMatrixBase<T> *);
2276 
2281  template <typename T>
2282  static yes_type
2283  check_for_block_matrix(const BlockSparsityPatternBase<T> *);
2284 
2289  template <typename T>
2290  static yes_type
2291  check_for_block_matrix(const BlockSparseMatrixEZ<T> *);
2292 
2297  static no_type
2298  check_for_block_matrix(...);
2299 
2300 public:
2306  static const bool value =
2307  (sizeof(check_for_block_matrix(static_cast<MatrixType *>(nullptr))) ==
2308  sizeof(yes_type));
2309 };
2310 
2311 // instantiation of the static member
2312 template <typename MatrixType>
2314 
2315 template <typename number>
2316 template <typename MatrixType>
2317 inline void
2319  const FullMatrix<number> & local_matrix,
2320  const std::vector<size_type> &local_dof_indices,
2321  MatrixType & global_matrix) const
2322 {
2323  // create a dummy and hand on to the function actually implementing this
2324  // feature in the cm.templates.h file.
2325  Vector<typename MatrixType::value_type> dummy(0);
2327  local_matrix,
2328  dummy,
2329  local_dof_indices,
2330  global_matrix,
2331  dummy,
2332  false,
2333  std::integral_constant<bool, IsBlockMatrix<MatrixType>::value>());
2334 }
2335 
2336 template <typename number>
2337 template <typename MatrixType, typename VectorType>
2338 inline void
2340  const FullMatrix<number> & local_matrix,
2341  const Vector<number> & local_vector,
2342  const std::vector<size_type> &local_dof_indices,
2343  MatrixType & global_matrix,
2344  VectorType & global_vector,
2345  bool use_inhomogeneities_for_rhs) const
2346 {
2347  // enter the internal function with the respective block information set,
2348  // the actual implementation follows in the cm.templates.h file.
2350  local_matrix,
2351  local_vector,
2352  local_dof_indices,
2353  global_matrix,
2354  global_vector,
2355  use_inhomogeneities_for_rhs,
2356  std::integral_constant<bool, IsBlockMatrix<MatrixType>::value>());
2357 }
2358 
2359 template <typename number>
2360 template <typename SparsityPatternType>
2361 inline void
2363  const std::vector<size_type> &local_dof_indices,
2364  SparsityPatternType & sparsity_pattern,
2365  const bool keep_constrained_entries,
2366  const Table<2, bool> & dof_mask) const
2367 {
2368  // enter the internal function with the respective block information set,
2369  // the actual implementation follows in the cm.templates.h file.
2371  local_dof_indices,
2372  sparsity_pattern,
2373  keep_constrained_entries,
2374  dof_mask,
2375  std::integral_constant<bool, IsBlockMatrix<SparsityPatternType>::value>());
2376 }
2377 
2379 
2380 #endif
const types::global_dof_index invalid_size_type
Definition: types.h:200
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:1403
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:819
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()