Reference documentation for deal.II version Git 8236ec9b21 2019-11-20 14:24:58 -0500
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
affine_constraints.h
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 
21 #include <deal.II/base/exceptions.h>
22 #include <deal.II/base/index_set.h>
23 #include <deal.II/base/subscriptor.h>
24 #include <deal.II/base/table.h>
25 #include <deal.II/base/template_constraints.h>
26 
27 #include <deal.II/lac/vector.h>
28 #include <deal.II/lac/vector_element_access.h>
29 
30 #include <boost/range/iterator_range.hpp>
31 
32 #include <set>
33 #include <utility>
34 #include <vector>
35 
36 DEAL_II_NAMESPACE_OPEN
37 
38 // Forward declarations
39 #ifndef DOXYGEN
40 template <typename>
41 class FullMatrix;
42 class SparsityPattern;
46 template <typename number>
47 class SparseMatrix;
48 template <typename number>
49 class BlockSparseMatrix;
50 
51 namespace internals
52 {
53  template <typename number>
54  class GlobalRowsFromLocal;
55 }
56 
57 namespace internal
58 {
59  namespace AffineConstraintsImplementation
60  {
61  template <class VectorType>
62  void
63  set_zero_all(const std::vector<types::global_dof_index> &cm,
64  VectorType & vec);
65 
66  template <class T>
67  void
68  set_zero_all(const std::vector<types::global_dof_index> &cm,
69  ::Vector<T> & vec);
70 
71  template <class T>
72  void
73  set_zero_all(const std::vector<types::global_dof_index> &cm,
74  ::BlockVector<T> & vec);
75  } // namespace AffineConstraintsImplementation
76 } // namespace internal
77 
78 
79 template <typename number>
80 class AffineConstraints;
81 #endif
82 
89 using ConstraintMatrix DEAL_II_DEPRECATED = AffineConstraints<double>;
90 // Note: Unfortunately, we cannot move this compatibility alias into
91 // constraint_matrix.h directly. This would break a lot of user projects
92 // that include constraint_matrix.h transitively due to various deal.II
93 // headers that include the file.
94 
95 
96 // TODO[WB]: We should have a function of the kind
97 // AffineConstraints::add_constraint (const size_type constrained_dof,
98 // const std::vector<std::pair<size_type, number> > &entries,
99 // const number inhomogeneity = 0);
100 // rather than building up constraints piecemeal through add_line/add_entry
101 // etc. This would also eliminate the possibility of accidentally changing
102 // existing constraints into something pointless, see the discussion on the
103 // mailing list on "Tiny bug in interpolate_boundary_values" in Sept. 2010.
104 
179 template <typename number = double>
181 {
182 public:
187 
194  {
200 
207 
213  right_object_wins
214  };
215 
233  explicit AffineConstraints(const IndexSet &local_constraints = IndexSet());
234 
238  explicit AffineConstraints(const AffineConstraints &affine_constraints);
239 
243  AffineConstraints(AffineConstraints &&affine_constraints) = default; // NOLINT
244 
255  operator=(const AffineConstraints &) = delete;
256 
261  operator=(AffineConstraints &&affine_constraints) = default; // NOLINT
262 
269  void
270  copy_from(const AffineConstraints &other);
271 
278  void
279  reinit(const IndexSet &local_constraints = IndexSet());
280 
287  bool
288  can_store_line(const size_type line_n) const;
289 
296  const IndexSet &
297  get_local_lines() const;
298 
318  void
319  add_selected_constraints(const AffineConstraints &constraints_in,
320  const IndexSet & filter);
321 
331  void
332  add_line(const size_type line_n);
333 
346  void
347  add_lines(const std::vector<bool> &lines);
348 
361  void
362  add_lines(const std::set<size_type> &lines);
363 
376  void
377  add_lines(const IndexSet &lines);
378 
390  void
391  add_entry(const size_type line_n, const size_type column, const number value);
392 
398  void
399  add_entries(const size_type line_n,
400  const std::vector<std::pair<size_type, number>> &col_val_pairs);
401 
408  void
409  set_inhomogeneity(const size_type line_n, const number value);
410 
431  void
432  close();
433 
458  void
459  merge(
460  const AffineConstraints & other_constraints,
461  const MergeConflictBehavior merge_conflict_behavior = no_conflicts_allowed,
462  const bool allow_different_local_lines = false);
463 
476  void
477  shift(const size_type offset);
478 
486  void
487  clear();
488 
501  size_type
502  n_constraints() const;
503 
513  bool
514  is_constrained(const size_type line_n) const;
515 
527  bool
528  is_identity_constrained(const size_type line_n) const;
529 
536  bool
537  are_identity_constrained(const size_type line_n_1,
538  const size_type line_n_2) const;
539 
550  size_type
551  max_constraint_indirections() const;
552 
557  bool
558  is_inhomogeneously_constrained(const size_type index) const;
559 
565  bool
566  has_inhomogeneities() const;
567 
572  const std::vector<std::pair<size_type, number>> *
573  get_constraint_entries(const size_type line_n) const;
574 
579  number
580  get_inhomogeneity(const size_type line_n) const;
581 
602  void
603  print(std::ostream &out) const;
604 
617  void
618  write_dot(std::ostream &) const;
619 
624  std::size_t
625  memory_consumption() const;
626 
633  void
634  resolve_indices(std::vector<types::global_dof_index> &indices) const;
635 
661  void
662  condense(SparsityPattern &sparsity) const;
663 
667  void
668  condense(BlockSparsityPattern &sparsity) const;
669 
674  void
675  condense(DynamicSparsityPattern &sparsity) const;
676 
681  void
682  condense(BlockDynamicSparsityPattern &sparsity) const;
683 
691  void
692  condense(SparseMatrix<number> &matrix) const;
693 
697  void
698  condense(BlockSparseMatrix<number> &matrix) const;
699 
711  template <class VectorType>
712  void
713  condense(VectorType &vec) const;
714 
721  template <class VectorType>
722  void
723  condense(const VectorType &vec_ghosted, VectorType &output) const;
724 
737  template <class VectorType>
738  void
739  condense(SparseMatrix<number> &matrix, VectorType &vector) const;
740 
745  template <class BlockVectorType>
746  void
747  condense(BlockSparseMatrix<number> &matrix, BlockVectorType &vector) const;
748 
755  template <class VectorType>
756  void
757  set_zero(VectorType &vec) const;
758 
810  template <class InVector, class OutVector>
811  void
812  distribute_local_to_global(const InVector & local_vector,
813  const std::vector<size_type> &local_dof_indices,
814  OutVector & global_vector) const;
815 
859  template <typename VectorType>
860  void
861  distribute_local_to_global(const Vector<number> & local_vector,
862  const std::vector<size_type> &local_dof_indices,
863  VectorType & global_vector,
864  const FullMatrix<number> & local_matrix) const;
865 
885  template <typename VectorType>
886  void
887  distribute_local_to_global(
888  const Vector<number> & local_vector,
889  const std::vector<size_type> &local_dof_indices_row,
890  const std::vector<size_type> &local_dof_indices_col,
891  VectorType & global_vector,
892  const FullMatrix<number> & local_matrix,
893  bool diagonal = false) const;
894 
898  template <class VectorType>
899  void
900  distribute_local_to_global(const size_type index,
901  const number value,
902  VectorType & global_vector) const;
903 
932  template <typename ForwardIteratorVec,
933  typename ForwardIteratorInd,
934  class VectorType>
935  void
936  distribute_local_to_global(ForwardIteratorVec local_vector_begin,
937  ForwardIteratorVec local_vector_end,
938  ForwardIteratorInd local_indices_begin,
939  VectorType & global_vector) const;
940 
988  template <typename MatrixType>
989  void
990  distribute_local_to_global(const FullMatrix<number> & local_matrix,
991  const std::vector<size_type> &local_dof_indices,
992  MatrixType & global_matrix) const;
993 
1021  template <typename MatrixType>
1022  void
1023  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1024  const std::vector<size_type> &row_indices,
1025  const std::vector<size_type> &col_indices,
1026  MatrixType & global_matrix) const;
1027 
1044  template <typename MatrixType>
1045  void
1046  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1047  const std::vector<size_type> &row_indices,
1048  const AffineConstraints &column_affine_constraints,
1049  const std::vector<size_type> &column_indices,
1050  MatrixType & global_matrix) const;
1051 
1068  template <typename MatrixType, typename VectorType>
1069  void
1070  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1071  const Vector<number> & local_vector,
1072  const std::vector<size_type> &local_dof_indices,
1073  MatrixType & global_matrix,
1074  VectorType & global_vector,
1075  bool use_inhomogeneities_for_rhs = false) const;
1076 
1130  template <typename SparsityPatternType>
1131  void
1132  add_entries_local_to_global(
1133  const std::vector<size_type> &local_dof_indices,
1134  SparsityPatternType & sparsity_pattern,
1135  const bool keep_constrained_entries = true,
1136  const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1137 
1141  template <typename SparsityPatternType>
1142  void
1143  add_entries_local_to_global(
1144  const std::vector<size_type> &row_indices,
1145  const std::vector<size_type> &col_indices,
1146  SparsityPatternType & sparsity_pattern,
1147  const bool keep_constrained_entries = true,
1148  const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1149 
1169  template <typename ForwardIteratorVec,
1170  typename ForwardIteratorInd,
1171  class VectorType>
1172  void
1173  get_dof_values(const VectorType & global_vector,
1174  ForwardIteratorInd local_indices_begin,
1175  ForwardIteratorVec local_vector_begin,
1176  ForwardIteratorVec local_vector_end) const;
1177 
1199  template <class VectorType>
1200  void
1201  distribute(VectorType &vec) const;
1202 
1211  {
1216  using Entries = std::vector<std::pair<size_type, number>>;
1217 
1224 
1233 
1238 
1246  bool
1247  operator<(const ConstraintLine &) const;
1248 
1254  bool
1255  operator==(const ConstraintLine &) const;
1256 
1261  std::size_t
1262  memory_consumption() const;
1263 
1267  template <class Archive>
1268  void
1269  serialize(Archive &ar, const unsigned int)
1270  {
1271  ar &index &entries &inhomogeneity;
1272  }
1273  };
1274 
1278  using const_iterator = typename std::vector<ConstraintLine>::const_iterator;
1279 
1283  using LineRange = boost::iterator_range<const_iterator>;
1284 
1293  const LineRange
1294  get_lines() const;
1295 
1323  bool
1324  is_consistent_in_parallel(const std::vector<IndexSet> &locally_owned_dofs,
1325  const IndexSet & locally_active_dofs,
1326  const MPI_Comm mpi_communicator,
1327  const bool verbose = false) const;
1328 
1334  DeclException0(ExcMatrixIsClosed);
1340  DeclException0(ExcMatrixNotClosed);
1347  size_type,
1348  << "The specified line " << arg1 << " does not exist.");
1354  DeclException4(ExcEntryAlreadyExists,
1355  size_type,
1356  size_type,
1357  number,
1358  number,
1359  << "The entry for the indices " << arg1 << " and " << arg2
1360  << " already exists, but the values " << arg3 << " (old) and "
1361  << arg4 << " (new) differ "
1362  << "by " << (arg4 - arg3) << ".");
1368  DeclException2(ExcDoFConstrainedToConstrainedDoF,
1369  int,
1370  int,
1371  << "You tried to constrain DoF " << arg1 << " to DoF " << arg2
1372  << ", but that one is also constrained. This is not allowed!");
1378  DeclException1(ExcDoFIsConstrainedFromBothObjects,
1379  size_type,
1380  << "Degree of freedom " << arg1
1381  << " is constrained from both object in a merge operation.");
1387  DeclException1(ExcDoFIsConstrainedToConstrainedDoF,
1388  size_type,
1389  << "In the given argument a degree of freedom is constrained "
1390  << "to another DoF with number " << arg1
1391  << ", which however is constrained by this object. This is not"
1392  << " allowed.");
1398  DeclException1(ExcRowNotStoredHere,
1399  size_type,
1400  << "The index set given to this constraints object indicates "
1401  << "constraints for degree of freedom " << arg1
1402  << " should not be stored by this object, but a constraint "
1403  << "is being added.");
1404 
1410  DeclException2(ExcColumnNotStoredHere,
1411  size_type,
1412  size_type,
1413  << "The index set given to this constraints object indicates "
1414  << "constraints using degree of freedom " << arg2
1415  << " should not be stored by this object, but a constraint "
1416  << "for degree of freedom " << arg1 << " uses it.");
1417 
1423  DeclException2(ExcIncorrectConstraint,
1424  int,
1425  int,
1426  << "While distributing the constraint for DoF " << arg1
1427  << ", it turns out that one of the processors "
1428  << "who own the " << arg2 << " degrees of freedom that x_"
1429  << arg1 << " is constrained against does not know about "
1430  << "the constraint on x_" << arg1
1431  << ". Did you not initialize the AffineConstraints container "
1432  << "with the appropriate locally_relevant set so "
1433  << "that every processor who owns a DoF that constrains "
1434  << "another DoF also knows about this constraint?");
1435 
1436 private:
1448  std::vector<ConstraintLine> lines;
1449 
1482  std::vector<size_type> lines_cache;
1483 
1490 
1494  bool sorted;
1495 
1500  size_type
1501  calculate_line_index(const size_type line_n) const;
1502 
1507  template <typename MatrixType, typename VectorType>
1508  void
1509  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1510  const Vector<number> & local_vector,
1511  const std::vector<size_type> &local_dof_indices,
1512  MatrixType & global_matrix,
1513  VectorType & global_vector,
1514  bool use_inhomogeneities_for_rhs,
1515  std::integral_constant<bool, false>) const;
1516 
1521  template <typename MatrixType, typename VectorType>
1522  void
1523  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1524  const Vector<number> & local_vector,
1525  const std::vector<size_type> &local_dof_indices,
1526  MatrixType & global_matrix,
1527  VectorType & global_vector,
1528  bool use_inhomogeneities_for_rhs,
1529  std::integral_constant<bool, true>) const;
1530 
1535  template <typename SparsityPatternType>
1536  void
1537  add_entries_local_to_global(const std::vector<size_type> &local_dof_indices,
1538  SparsityPatternType & sparsity_pattern,
1539  const bool keep_constrained_entries,
1540  const Table<2, bool> &dof_mask,
1541  std::integral_constant<bool, false>) const;
1542 
1547  template <typename SparsityPatternType>
1548  void
1549  add_entries_local_to_global(const std::vector<size_type> &local_dof_indices,
1550  SparsityPatternType & sparsity_pattern,
1551  const bool keep_constrained_entries,
1552  const Table<2, bool> &dof_mask,
1553  std::integral_constant<bool, true>) const;
1554 
1562  void
1563  make_sorted_row_list(
1564  const std::vector<size_type> & local_dof_indices,
1565  internals::GlobalRowsFromLocal<number> &global_rows) const;
1566 
1574  void
1575  make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
1576  std::vector<size_type> & active_dofs) const;
1577 
1581  template <typename MatrixScalar, typename VectorScalar>
1582  typename ProductType<VectorScalar, MatrixScalar>::type
1583  resolve_vector_entry(
1584  const size_type i,
1585  const internals::GlobalRowsFromLocal<number> &global_rows,
1586  const Vector<VectorScalar> & local_vector,
1587  const std::vector<size_type> & local_dof_indices,
1588  const FullMatrix<MatrixScalar> & local_matrix) const;
1589 };
1590 
1591 /* ---------------- template and inline functions ----------------- */
1592 
1593 template <typename number>
1595  const IndexSet &local_constraints)
1596  : lines()
1597  , local_lines(local_constraints)
1598  , sorted(false)
1599 {
1600  // make sure the IndexSet is compressed. Otherwise this can lead to crashes
1601  // that are hard to find (only happen in release mode).
1602  // see tests/mpi/affine_constraints_crash_01
1604 }
1605 
1606 template <typename number>
1608  const AffineConstraints &affine_constraints)
1609  : Subscriptor()
1610  , lines(affine_constraints.lines)
1611  , lines_cache(affine_constraints.lines_cache)
1612  , local_lines(affine_constraints.local_lines)
1613  , sorted(affine_constraints.sorted)
1614 {}
1615 
1616 template <typename number>
1617 inline void
1619 {
1620  Assert(sorted == false, ExcMatrixIsClosed());
1621 
1622  // the following can happen when we compute with distributed meshes and dof
1623  // handlers and we constrain a degree of freedom whose number we don't have
1624  // locally. if we don't abort here the program will try to allocate several
1625  // terabytes of memory to resize the various arrays below :-)
1627  const size_type line_index = calculate_line_index(line_n);
1628 
1629  // check whether line already exists; it may, in which case we can just quit
1630  if (is_constrained(line_n))
1631  return;
1632 
1633  // if necessary enlarge vector of existing entries for cache
1634  if (line_index >= lines_cache.size())
1635  lines_cache.resize(std::max(2 * static_cast<size_type>(lines_cache.size()),
1636  line_index + 1),
1638 
1639  // push a new line to the end of the list
1640  lines.emplace_back();
1641  lines.back().index = line_n;
1642  lines.back().inhomogeneity = 0.;
1643  lines_cache[line_index] = lines.size() - 1;
1644 }
1645 
1646 template <typename number>
1647 inline void
1649  const size_type column,
1650  const number value)
1651 {
1652  Assert(sorted == false, ExcMatrixIsClosed());
1653  Assert(line_n != column,
1654  ExcMessage("Can't constrain a degree of freedom to itself"));
1655 
1656  // Ensure that the current line is present in the cache:
1657  const size_type line_index = calculate_line_index(line_n);
1658  Assert(line_index < lines_cache.size(),
1659  ExcMessage("The current AffineConstraints does not contain the line "
1660  "for the current entry. Call AffineConstraints::add_line "
1661  "before calling this function."));
1662 
1663  // if in debug mode, check whether an entry for this column already exists
1664  // and if it's the same as the one entered at present
1665  //
1666  // in any case: exit the function if an entry for this column already
1667  // exists, since we don't want to enter it twice
1669  ExcInternalError());
1671  ExcColumnNotStoredHere(line_n, column));
1672  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
1673  Assert(line_ptr->index == line_n, ExcInternalError());
1674  for (const auto &p : line_ptr->entries)
1675  if (p.first == column)
1676  {
1677  Assert(std::abs(p.second - value) < 1.e-14,
1678  ExcEntryAlreadyExists(line_n, column, p.second, value));
1679  return;
1680  }
1681 
1682  line_ptr->entries.emplace_back(column, value);
1683 }
1684 
1685 template <typename number>
1686 inline void
1688  const number value)
1689 {
1690  const size_type line_index = calculate_line_index(line_n);
1691  Assert(line_index < lines_cache.size() &&
1692  lines_cache[line_index] != numbers::invalid_size_type,
1693  ExcMessage("call add_line() before calling set_inhomogeneity()"));
1694  Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
1695  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
1696  line_ptr->inhomogeneity = value;
1697 }
1698 
1699 template <typename number>
1700 template <class VectorType>
1701 inline void
1703 {
1704  // since lines is a private member, we cannot pass it to the functions
1705  // above. therefore, copy the content which is cheap
1706  std::vector<size_type> constrained_lines(lines.size());
1707  for (unsigned int i = 0; i < lines.size(); ++i)
1708  constrained_lines[i] = lines[i].index;
1709  internal::AffineConstraintsImplementation::set_zero_all(constrained_lines,
1710  vec);
1711 }
1712 
1713 template <typename number>
1716 {
1717  return lines.size();
1718 }
1719 
1720 template <typename number>
1721 inline bool
1723 {
1724  const size_type line_index = calculate_line_index(index);
1725  return ((line_index < lines_cache.size()) &&
1726  (lines_cache[line_index] != numbers::invalid_size_type));
1727 }
1728 
1729 template <typename number>
1730 inline bool
1732  const size_type line_n) const
1733 {
1734  // check whether the entry is constrained. could use is_constrained, but
1735  // that means computing the line index twice
1736  const size_type line_index = calculate_line_index(line_n);
1737  if (line_index >= lines_cache.size() ||
1738  lines_cache[line_index] == numbers::invalid_size_type)
1739  return false;
1740  else
1741  {
1742  Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
1743  return !(lines[lines_cache[line_index]].inhomogeneity == number(0.));
1744  }
1745 }
1746 
1747 template <typename number>
1748 inline const std::vector<std::pair<types::global_dof_index, number>> *
1750 {
1751  // check whether the entry is constrained. could use is_constrained, but
1752  // that means computing the line index twice
1753  const size_type line_index = calculate_line_index(line_n);
1754  if (line_index >= lines_cache.size() ||
1755  lines_cache[line_index] == numbers::invalid_size_type)
1756  return nullptr;
1757  else
1758  return &lines[lines_cache[line_index]].entries;
1759 }
1760 
1761 template <typename number>
1762 inline number
1764 {
1765  // check whether the entry is constrained. could use is_constrained, but
1766  // that means computing the line index twice
1767  const size_type line_index = calculate_line_index(line_n);
1768  if (line_index >= lines_cache.size() ||
1769  lines_cache[line_index] == numbers::invalid_size_type)
1770  return 0;
1771  else
1772  return lines[lines_cache[line_index]].inhomogeneity;
1773 }
1774 
1775 template <typename number>
1778 {
1779  // IndexSet is unused (serial case)
1780  if (!local_lines.size())
1781  return line_n;
1782 
1784 
1785  return local_lines.index_within_set(line_n);
1786 }
1787 
1788 template <typename number>
1789 inline bool
1791 {
1792  return !local_lines.size() || local_lines.is_element(line_n);
1793 }
1794 
1795 template <typename number>
1796 inline const IndexSet &
1798 {
1799  return local_lines;
1800 }
1801 
1802 template <typename number>
1803 template <class VectorType>
1804 inline void
1806  const size_type index,
1807  const number value,
1808  VectorType & global_vector) const
1809 {
1810  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
1811 
1812  if (is_constrained(index) == false)
1813  global_vector(index) += value;
1814  else
1815  {
1816  const ConstraintLine &position =
1818  for (size_type j = 0; j < position.entries.size(); ++j)
1819  global_vector(position.entries[j].first) +=
1820  value * position.entries[j].second;
1821  }
1822 }
1823 
1824 template <typename number>
1825 template <typename ForwardIteratorVec,
1826  typename ForwardIteratorInd,
1827  class VectorType>
1828 inline void
1830  ForwardIteratorVec local_vector_begin,
1831  ForwardIteratorVec local_vector_end,
1832  ForwardIteratorInd local_indices_begin,
1833  VectorType & global_vector) const
1834 {
1835  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
1836  for (; local_vector_begin != local_vector_end;
1837  ++local_vector_begin, ++local_indices_begin)
1838  {
1839  if (is_constrained(*local_indices_begin) == false)
1840  internal::ElementAccess<VectorType>::add(*local_vector_begin,
1841  *local_indices_begin,
1842  global_vector);
1843  else
1844  {
1845  const ConstraintLine &position =
1846  lines[lines_cache[calculate_line_index(*local_indices_begin)]];
1847  for (size_type j = 0; j < position.entries.size(); ++j)
1848  internal::ElementAccess<VectorType>::add(
1849  (*local_vector_begin) * position.entries[j].second,
1850  position.entries[j].first,
1851  global_vector);
1852  }
1853  }
1854 }
1855 
1856 template <typename number>
1857 template <class InVector, class OutVector>
1858 inline void
1860  const InVector & local_vector,
1861  const std::vector<size_type> &local_dof_indices,
1862  OutVector & global_vector) const
1863 {
1864  Assert(local_vector.size() == local_dof_indices.size(),
1865  ExcDimensionMismatch(local_vector.size(), local_dof_indices.size()));
1866  distribute_local_to_global(local_vector.begin(),
1867  local_vector.end(),
1868  local_dof_indices.begin(),
1869  global_vector);
1870 }
1871 
1872 template <typename number>
1873 template <typename ForwardIteratorVec,
1874  typename ForwardIteratorInd,
1875  class VectorType>
1876 inline void
1878  const VectorType & global_vector,
1879  ForwardIteratorInd local_indices_begin,
1880  ForwardIteratorVec local_vector_begin,
1881  ForwardIteratorVec local_vector_end) const
1882 {
1883  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
1884  for (; local_vector_begin != local_vector_end;
1885  ++local_vector_begin, ++local_indices_begin)
1886  {
1887  if (is_constrained(*local_indices_begin) == false)
1888  *local_vector_begin = global_vector(*local_indices_begin);
1889  else
1890  {
1891  const ConstraintLine &position =
1892  lines[lines_cache[calculate_line_index(*local_indices_begin)]];
1893  typename VectorType::value_type value = position.inhomogeneity;
1894  for (size_type j = 0; j < position.entries.size(); ++j)
1895  value += (global_vector(position.entries[j].first) *
1896  position.entries[j].second);
1897  *local_vector_begin = value;
1898  }
1899  }
1900 }
1901 
1902 template <typename MatrixType>
1904 template <typename SparsityPatternType>
1906 template <typename number>
1908 
1927 template <typename MatrixType>
1929 {
1930 private:
1931  struct yes_type
1932  {
1933  char c[1];
1934  };
1935  struct no_type
1936  {
1937  char c[2];
1938  };
1939 
1945  template <typename T>
1946  static yes_type
1947  check_for_block_matrix(const BlockMatrixBase<T> *);
1948 
1953  template <typename T>
1954  static yes_type
1955  check_for_block_matrix(const BlockSparsityPatternBase<T> *);
1956 
1961  template <typename T>
1962  static yes_type
1963  check_for_block_matrix(const BlockSparseMatrixEZ<T> *);
1964 
1969  static no_type
1970  check_for_block_matrix(...);
1971 
1972 public:
1978  static const bool value =
1979  (sizeof(check_for_block_matrix(static_cast<MatrixType *>(nullptr))) ==
1980  sizeof(yes_type));
1981 };
1982 
1983 // instantiation of the static member
1984 template <typename MatrixType>
1986 
1987 template <typename number>
1988 template <typename MatrixType>
1989 inline void
1991  const FullMatrix<number> & local_matrix,
1992  const std::vector<size_type> &local_dof_indices,
1993  MatrixType & global_matrix) const
1994 {
1995  // create a dummy and hand on to the function actually implementing this
1996  // feature in the cm.templates.h file.
1997  Vector<typename MatrixType::value_type> dummy(0);
1999  local_matrix,
2000  dummy,
2001  local_dof_indices,
2002  global_matrix,
2003  dummy,
2004  false,
2005  std::integral_constant<bool, IsBlockMatrix<MatrixType>::value>());
2006 }
2007 
2008 template <typename number>
2009 template <typename MatrixType, typename VectorType>
2010 inline void
2012  const FullMatrix<number> & local_matrix,
2013  const Vector<number> & local_vector,
2014  const std::vector<size_type> &local_dof_indices,
2015  MatrixType & global_matrix,
2016  VectorType & global_vector,
2017  bool use_inhomogeneities_for_rhs) const
2018 {
2019  // enter the internal function with the respective block information set,
2020  // the actual implementation follows in the cm.templates.h file.
2022  local_matrix,
2023  local_vector,
2024  local_dof_indices,
2025  global_matrix,
2026  global_vector,
2027  use_inhomogeneities_for_rhs,
2028  std::integral_constant<bool, IsBlockMatrix<MatrixType>::value>());
2029 }
2030 
2031 template <typename number>
2032 template <typename SparsityPatternType>
2033 inline void
2035  const std::vector<size_type> &local_dof_indices,
2036  SparsityPatternType & sparsity_pattern,
2037  const bool keep_constrained_entries,
2038  const Table<2, bool> & dof_mask) const
2039 {
2040  // enter the internal function with the respective block information set,
2041  // the actual implementation follows in the cm.templates.h file.
2043  local_dof_indices,
2044  sparsity_pattern,
2045  keep_constrained_entries,
2046  dof_mask,
2047  std::integral_constant<bool, IsBlockMatrix<SparsityPatternType>::value>());
2048 }
2049 
2050 DEAL_II_NAMESPACE_CLOSE
2051 
2052 #endif
const types::global_dof_index invalid_size_type
Definition: types.h:196
std::vector< ConstraintLine > lines
void set_zero(VectorType &vec) const
bool is_constrained(const size_type line_n) const
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
bool can_store_line(const size_type line_n) const
void add_entry(const size_type line_n, const size_type column, const number value)
static ::ExceptionBase & ExcColumnNotStoredHere(size_type arg1, size_type arg2)
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()
static ::ExceptionBase & ExcMatrixIsClosed()
size_type size() const
Definition: index_set.h:1607
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:518
#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:1895
#define DeclException0(Exception0)
Definition: exceptions.h:473
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
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:89
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:1615
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:587
Definition: table.h:699
bool is_element(const size_type index) const
Definition: index_set.h:1739
DEAL_II_CONSTEXPR TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
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
static ::ExceptionBase & ExcEntryAlreadyExists(size_type arg1, size_type arg2, number arg3, number arg4)
static ::ExceptionBase & ExcInternalError()