Reference documentation for deal.II version 9.1.1
\(\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 template <typename>
39 class FullMatrix;
40 class SparsityPattern;
44 template <typename number>
45 class SparseMatrix;
46 template <typename number>
48 
49 namespace internals
50 {
51  template <typename number>
52  class GlobalRowsFromLocal;
53 }
54 
55 
56 template <typename number>
57 class AffineConstraints;
58 
65 using ConstraintMatrix DEAL_II_DEPRECATED = AffineConstraints<double>;
66 // Note: Unfortunately, we cannot move this compatibility alias into
67 // constraint_matrix.h directly. This would break a lot of user projects
68 // that include constraint_matrix.h transitively due to various deal.II
69 // headers that include the file.
70 
71 
72 // TODO[WB]: We should have a function of the kind
73 // AffineConstraints::add_constraint (const size_type constrained_dof,
74 // const std::vector<std::pair<size_type, number> > &entries,
75 // const number inhomogeneity = 0);
76 // rather than building up constraints piecemeal through add_line/add_entry
77 // etc. This would also eliminate the possibility of accidentally changing
78 // existing constraints into something pointless, see the discussion on the
79 // mailing list on "Tiny bug in interpolate_boundary_values" in Sept. 2010.
80 
155 template <typename number = double>
156 class AffineConstraints : public Subscriptor
157 {
158 public:
163 
170  {
176 
183 
189  right_object_wins
190  };
191 
209  explicit AffineConstraints(const IndexSet &local_constraints = IndexSet());
210 
214  explicit AffineConstraints(const AffineConstraints &affine_constraints);
215 
219  AffineConstraints(AffineConstraints &&affine_constraints) = default; // NOLINT
220 
231  operator=(const AffineConstraints &) = delete;
232 
237  operator=(AffineConstraints &&affine_constraints) = default; // NOLINT
238 
245  void
246  copy_from(const AffineConstraints &other);
247 
254  void
255  reinit(const IndexSet &local_constraints = IndexSet());
256 
263  bool
264  can_store_line(const size_type line_n) const;
265 
272  const IndexSet &
273  get_local_lines() const;
274 
294  void
295  add_selected_constraints(const AffineConstraints &constraints_in,
296  const IndexSet & filter);
297 
307  void
308  add_line(const size_type line_n);
309 
322  void
323  add_lines(const std::vector<bool> &lines);
324 
337  void
338  add_lines(const std::set<size_type> &lines);
339 
352  void
353  add_lines(const IndexSet &lines);
354 
366  void
367  add_entry(const size_type line_n, const size_type column, const number value);
368 
374  void
375  add_entries(const size_type line_n,
376  const std::vector<std::pair<size_type, number>> &col_val_pairs);
377 
384  void
385  set_inhomogeneity(const size_type line_n, const number value);
386 
407  void
408  close();
409 
434  void
435  merge(
436  const AffineConstraints & other_constraints,
437  const MergeConflictBehavior merge_conflict_behavior = no_conflicts_allowed,
438  const bool allow_different_local_lines = false);
439 
452  void
453  shift(const size_type offset);
454 
462  void
463  clear();
464 
477  size_type
478  n_constraints() const;
479 
489  bool
490  is_constrained(const size_type line_n) const;
491 
503  bool
504  is_identity_constrained(const size_type line_n) const;
505 
512  bool
513  are_identity_constrained(const size_type line_n_1,
514  const size_type line_n_2) const;
515 
526  size_type
527  max_constraint_indirections() const;
528 
533  bool
534  is_inhomogeneously_constrained(const size_type index) const;
535 
541  bool
542  has_inhomogeneities() const;
543 
548  const std::vector<std::pair<size_type, number>> *
549  get_constraint_entries(const size_type line_n) const;
550 
555  number
556  get_inhomogeneity(const size_type line_n) const;
557 
578  void
579  print(std::ostream &out) const;
580 
593  void
594  write_dot(std::ostream &) const;
595 
600  std::size_t
601  memory_consumption() const;
602 
609  void
610  resolve_indices(std::vector<types::global_dof_index> &indices) const;
611 
637  void
638  condense(SparsityPattern &sparsity) const;
639 
643  void
644  condense(BlockSparsityPattern &sparsity) const;
645 
650  void
651  condense(DynamicSparsityPattern &sparsity) const;
652 
657  void
658  condense(BlockDynamicSparsityPattern &sparsity) const;
659 
667  void
668  condense(SparseMatrix<number> &matrix) const;
669 
673  void
674  condense(BlockSparseMatrix<number> &matrix) const;
675 
687  template <class VectorType>
688  void
689  condense(VectorType &vec) const;
690 
697  template <class VectorType>
698  void
699  condense(const VectorType &vec_ghosted, VectorType &output) const;
700 
713  template <class VectorType>
714  void
715  condense(SparseMatrix<number> &matrix, VectorType &vector) const;
716 
721  template <class BlockVectorType>
722  void
723  condense(BlockSparseMatrix<number> &matrix, BlockVectorType &vector) const;
724 
731  template <class VectorType>
732  void
733  set_zero(VectorType &vec) const;
734 
786  template <class InVector, class OutVector>
787  void
788  distribute_local_to_global(const InVector & local_vector,
789  const std::vector<size_type> &local_dof_indices,
790  OutVector & global_vector) const;
791 
835  template <typename VectorType>
836  void
837  distribute_local_to_global(const Vector<number> & local_vector,
838  const std::vector<size_type> &local_dof_indices,
839  VectorType & global_vector,
840  const FullMatrix<number> & local_matrix) const;
841 
861  template <typename VectorType>
862  void
863  distribute_local_to_global(
864  const Vector<number> & local_vector,
865  const std::vector<size_type> &local_dof_indices_row,
866  const std::vector<size_type> &local_dof_indices_col,
867  VectorType & global_vector,
868  const FullMatrix<number> & local_matrix,
869  bool diagonal = false) const;
870 
874  template <class VectorType>
875  void
876  distribute_local_to_global(const size_type index,
877  const number value,
878  VectorType & global_vector) const;
879 
908  template <typename ForwardIteratorVec,
909  typename ForwardIteratorInd,
910  class VectorType>
911  void
912  distribute_local_to_global(ForwardIteratorVec local_vector_begin,
913  ForwardIteratorVec local_vector_end,
914  ForwardIteratorInd local_indices_begin,
915  VectorType & global_vector) const;
916 
964  template <typename MatrixType>
965  void
966  distribute_local_to_global(const FullMatrix<number> & local_matrix,
967  const std::vector<size_type> &local_dof_indices,
968  MatrixType & global_matrix) const;
969 
997  template <typename MatrixType>
998  void
999  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1000  const std::vector<size_type> &row_indices,
1001  const std::vector<size_type> &col_indices,
1002  MatrixType & global_matrix) const;
1003 
1020  template <typename MatrixType>
1021  void
1022  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1023  const std::vector<size_type> &row_indices,
1024  const AffineConstraints &column_affine_constraints,
1025  const std::vector<size_type> &column_indices,
1026  MatrixType & global_matrix) const;
1027 
1044  template <typename MatrixType, typename VectorType>
1045  void
1046  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1047  const Vector<number> & local_vector,
1048  const std::vector<size_type> &local_dof_indices,
1049  MatrixType & global_matrix,
1050  VectorType & global_vector,
1051  bool use_inhomogeneities_for_rhs = false) const;
1052 
1106  template <typename SparsityPatternType>
1107  void
1108  add_entries_local_to_global(
1109  const std::vector<size_type> &local_dof_indices,
1110  SparsityPatternType & sparsity_pattern,
1111  const bool keep_constrained_entries = true,
1112  const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1113 
1117  template <typename SparsityPatternType>
1118  void
1119  add_entries_local_to_global(
1120  const std::vector<size_type> &row_indices,
1121  const std::vector<size_type> &col_indices,
1122  SparsityPatternType & sparsity_pattern,
1123  const bool keep_constrained_entries = true,
1124  const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1125 
1145  template <typename ForwardIteratorVec,
1146  typename ForwardIteratorInd,
1147  class VectorType>
1148  void
1149  get_dof_values(const VectorType & global_vector,
1150  ForwardIteratorInd local_indices_begin,
1151  ForwardIteratorVec local_vector_begin,
1152  ForwardIteratorVec local_vector_end) const;
1153 
1175  template <class VectorType>
1176  void
1177  distribute(VectorType &vec) const;
1178 
1187  {
1192  using Entries = std::vector<std::pair<size_type, number>>;
1193 
1200 
1209 
1214 
1222  bool
1223  operator<(const ConstraintLine &) const;
1224 
1230  bool
1231  operator==(const ConstraintLine &) const;
1232 
1237  std::size_t
1238  memory_consumption() const;
1239 
1243  template <class Archive>
1244  void
1245  serialize(Archive &ar, const unsigned int)
1246  {
1247  ar &index &entries &inhomogeneity;
1248  }
1249  };
1250 
1254  using const_iterator = typename std::vector<ConstraintLine>::const_iterator;
1255 
1259  using LineRange = boost::iterator_range<const_iterator>;
1260 
1269  const LineRange
1270  get_lines() const;
1271 
1299  bool
1300  is_consistent_in_parallel(const std::vector<IndexSet> &locally_owned_dofs,
1301  const IndexSet & locally_active_dofs,
1302  const MPI_Comm mpi_communicator,
1303  const bool verbose = false) const;
1304 
1310  DeclException0(ExcMatrixIsClosed);
1316  DeclException0(ExcMatrixNotClosed);
1323  size_type,
1324  << "The specified line " << arg1 << " does not exist.");
1330  DeclException4(ExcEntryAlreadyExists,
1331  size_type,
1332  size_type,
1333  number,
1334  number,
1335  << "The entry for the indices " << arg1 << " and " << arg2
1336  << " already exists, but the values " << arg3 << " (old) and "
1337  << arg4 << " (new) differ "
1338  << "by " << (arg4 - arg3) << ".");
1344  DeclException2(ExcDoFConstrainedToConstrainedDoF,
1345  int,
1346  int,
1347  << "You tried to constrain DoF " << arg1 << " to DoF " << arg2
1348  << ", but that one is also constrained. This is not allowed!");
1354  DeclException1(ExcDoFIsConstrainedFromBothObjects,
1355  size_type,
1356  << "Degree of freedom " << arg1
1357  << " is constrained from both object in a merge operation.");
1363  DeclException1(ExcDoFIsConstrainedToConstrainedDoF,
1364  size_type,
1365  << "In the given argument a degree of freedom is constrained "
1366  << "to another DoF with number " << arg1
1367  << ", which however is constrained by this object. This is not"
1368  << " allowed.");
1374  DeclException1(ExcRowNotStoredHere,
1375  size_type,
1376  << "The index set given to this constraints object indicates "
1377  << "constraints for degree of freedom " << arg1
1378  << " should not be stored by this object, but a constraint "
1379  << "is being added.");
1380 
1386  DeclException2(ExcColumnNotStoredHere,
1387  size_type,
1388  size_type,
1389  << "The index set given to this constraints object indicates "
1390  << "constraints using degree of freedom " << arg2
1391  << " should not be stored by this object, but a constraint "
1392  << "for degree of freedom " << arg1 << " uses it.");
1393 
1399  DeclException2(ExcIncorrectConstraint,
1400  int,
1401  int,
1402  << "While distributing the constraint for DoF " << arg1
1403  << ", it turns out that one of the processors "
1404  << "who own the " << arg2 << " degrees of freedom that x_"
1405  << arg1 << " is constrained against does not know about "
1406  << "the constraint on x_" << arg1
1407  << ". Did you not initialize the AffineConstraints container "
1408  << "with the appropriate locally_relevant set so "
1409  << "that every processor who owns a DoF that constrains "
1410  << "another DoF also knows about this constraint?");
1411 
1412 private:
1424  std::vector<ConstraintLine> lines;
1425 
1458  std::vector<size_type> lines_cache;
1459 
1466 
1470  bool sorted;
1471 
1476  size_type
1477  calculate_line_index(const size_type line_n) const;
1478 
1483  template <typename MatrixType, typename VectorType>
1484  void
1485  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1486  const Vector<number> & local_vector,
1487  const std::vector<size_type> &local_dof_indices,
1488  MatrixType & global_matrix,
1489  VectorType & global_vector,
1490  bool use_inhomogeneities_for_rhs,
1491  std::integral_constant<bool, false>) const;
1492 
1497  template <typename MatrixType, typename VectorType>
1498  void
1499  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1500  const Vector<number> & local_vector,
1501  const std::vector<size_type> &local_dof_indices,
1502  MatrixType & global_matrix,
1503  VectorType & global_vector,
1504  bool use_inhomogeneities_for_rhs,
1505  std::integral_constant<bool, true>) const;
1506 
1511  template <typename SparsityPatternType>
1512  void
1513  add_entries_local_to_global(const std::vector<size_type> &local_dof_indices,
1514  SparsityPatternType & sparsity_pattern,
1515  const bool keep_constrained_entries,
1516  const Table<2, bool> &dof_mask,
1517  std::integral_constant<bool, false>) const;
1518 
1523  template <typename SparsityPatternType>
1524  void
1525  add_entries_local_to_global(const std::vector<size_type> &local_dof_indices,
1526  SparsityPatternType & sparsity_pattern,
1527  const bool keep_constrained_entries,
1528  const Table<2, bool> &dof_mask,
1529  std::integral_constant<bool, true>) const;
1530 
1538  void
1539  make_sorted_row_list(
1540  const std::vector<size_type> & local_dof_indices,
1541  internals::GlobalRowsFromLocal<number> &global_rows) const;
1542 
1550  void
1551  make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
1552  std::vector<size_type> & active_dofs) const;
1553 
1557  template <typename MatrixScalar, typename VectorScalar>
1558  typename ProductType<VectorScalar, MatrixScalar>::type
1559  resolve_vector_entry(
1560  const size_type i,
1561  const internals::GlobalRowsFromLocal<number> &global_rows,
1562  const Vector<VectorScalar> & local_vector,
1563  const std::vector<size_type> & local_dof_indices,
1564  const FullMatrix<MatrixScalar> & local_matrix) const;
1565 };
1566 
1567 /* ---------------- template and inline functions ----------------- */
1568 
1569 template <typename number>
1571  const IndexSet &local_constraints)
1572  : lines()
1573  , local_lines(local_constraints)
1574  , sorted(false)
1575 {
1576  // make sure the IndexSet is compressed. Otherwise this can lead to crashes
1577  // that are hard to find (only happen in release mode).
1578  // see tests/mpi/affine_constraints_crash_01
1580 }
1581 
1582 template <typename number>
1584  const AffineConstraints &affine_constraints)
1585  : Subscriptor()
1586  , lines(affine_constraints.lines)
1587  , lines_cache(affine_constraints.lines_cache)
1588  , local_lines(affine_constraints.local_lines)
1589  , sorted(affine_constraints.sorted)
1590 {}
1591 
1592 template <typename number>
1593 inline void
1595 {
1596  Assert(sorted == false, ExcMatrixIsClosed());
1597 
1598  // the following can happen when we compute with distributed meshes and dof
1599  // handlers and we constrain a degree of freedom whose number we don't have
1600  // locally. if we don't abort here the program will try to allocate several
1601  // terabytes of memory to resize the various arrays below :-)
1603  const size_type line_index = calculate_line_index(line_n);
1604 
1605  // check whether line already exists; it may, in which case we can just quit
1606  if (is_constrained(line_n))
1607  return;
1608 
1609  // if necessary enlarge vector of existing entries for cache
1610  if (line_index >= lines_cache.size())
1611  lines_cache.resize(std::max(2 * static_cast<size_type>(lines_cache.size()),
1612  line_index + 1),
1614 
1615  // push a new line to the end of the list
1616  lines.emplace_back();
1617  lines.back().index = line_n;
1618  lines.back().inhomogeneity = 0.;
1619  lines_cache[line_index] = lines.size() - 1;
1620 }
1621 
1622 template <typename number>
1623 inline void
1625  const size_type column,
1626  const number value)
1627 {
1628  Assert(sorted == false, ExcMatrixIsClosed());
1629  Assert(line_n != column,
1630  ExcMessage("Can't constrain a degree of freedom to itself"));
1631 
1632  // Ensure that the current line is present in the cache:
1633  const size_type line_index = calculate_line_index(line_n);
1634  Assert(line_index < lines_cache.size(),
1635  ExcMessage("The current AffineConstraints does not contain the line "
1636  "for the current entry. Call AffineConstraints::add_line "
1637  "before calling this function."));
1638 
1639  // if in debug mode, check whether an entry for this column already exists
1640  // and if it's the same as the one entered at present
1641  //
1642  // in any case: exit the function if an entry for this column already
1643  // exists, since we don't want to enter it twice
1645  ExcInternalError());
1647  ExcColumnNotStoredHere(line_n, column));
1648  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
1649  Assert(line_ptr->index == line_n, ExcInternalError());
1650  for (const auto &p : line_ptr->entries)
1651  if (p.first == column)
1652  {
1653  Assert(std::abs(p.second - value) < 1.e-14,
1654  ExcEntryAlreadyExists(line_n, column, p.second, value));
1655  return;
1656  }
1657 
1658  line_ptr->entries.emplace_back(column, value);
1659 }
1660 
1661 template <typename number>
1662 inline void
1664  const number value)
1665 {
1666  const size_type line_index = calculate_line_index(line_n);
1667  Assert(line_index < lines_cache.size() &&
1668  lines_cache[line_index] != numbers::invalid_size_type,
1669  ExcMessage("call add_line() before calling set_inhomogeneity()"));
1670  Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
1671  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
1672  line_ptr->inhomogeneity = value;
1673 }
1674 
1675 template <typename number>
1678 {
1679  return lines.size();
1680 }
1681 
1682 template <typename number>
1683 inline bool
1685 {
1686  const size_type line_index = calculate_line_index(index);
1687  return ((line_index < lines_cache.size()) &&
1688  (lines_cache[line_index] != numbers::invalid_size_type));
1689 }
1690 
1691 template <typename number>
1692 inline bool
1694  const size_type line_n) const
1695 {
1696  // check whether the entry is constrained. could use is_constrained, but
1697  // that means computing the line index twice
1698  const size_type line_index = calculate_line_index(line_n);
1699  if (line_index >= lines_cache.size() ||
1700  lines_cache[line_index] == numbers::invalid_size_type)
1701  return false;
1702  else
1703  {
1704  Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
1705  return !(lines[lines_cache[line_index]].inhomogeneity == number(0.));
1706  }
1707 }
1708 
1709 template <typename number>
1710 inline const std::vector<std::pair<types::global_dof_index, number>> *
1712 {
1713  // check whether the entry is constrained. could use is_constrained, but
1714  // that means computing the line index twice
1715  const size_type line_index = calculate_line_index(line_n);
1716  if (line_index >= lines_cache.size() ||
1717  lines_cache[line_index] == numbers::invalid_size_type)
1718  return nullptr;
1719  else
1720  return &lines[lines_cache[line_index]].entries;
1721 }
1722 
1723 template <typename number>
1724 inline number
1726 {
1727  // check whether the entry is constrained. could use is_constrained, but
1728  // that means computing the line index twice
1729  const size_type line_index = calculate_line_index(line_n);
1730  if (line_index >= lines_cache.size() ||
1731  lines_cache[line_index] == numbers::invalid_size_type)
1732  return 0;
1733  else
1734  return lines[lines_cache[line_index]].inhomogeneity;
1735 }
1736 
1737 template <typename number>
1740 {
1741  // IndexSet is unused (serial case)
1742  if (!local_lines.size())
1743  return line_n;
1744 
1746 
1747  return local_lines.index_within_set(line_n);
1748 }
1749 
1750 template <typename number>
1751 inline bool
1753 {
1754  return !local_lines.size() || local_lines.is_element(line_n);
1755 }
1756 
1757 template <typename number>
1758 inline const IndexSet &
1760 {
1761  return local_lines;
1762 }
1763 
1764 template <typename number>
1765 template <class VectorType>
1766 inline void
1768  const size_type index,
1769  const number value,
1770  VectorType & global_vector) const
1771 {
1772  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
1773 
1774  if (is_constrained(index) == false)
1775  global_vector(index) += value;
1776  else
1777  {
1778  const ConstraintLine &position =
1780  for (size_type j = 0; j < position.entries.size(); ++j)
1781  global_vector(position.entries[j].first) +=
1782  value * position.entries[j].second;
1783  }
1784 }
1785 
1786 template <typename number>
1787 template <typename ForwardIteratorVec,
1788  typename ForwardIteratorInd,
1789  class VectorType>
1790 inline void
1792  ForwardIteratorVec local_vector_begin,
1793  ForwardIteratorVec local_vector_end,
1794  ForwardIteratorInd local_indices_begin,
1795  VectorType & global_vector) const
1796 {
1797  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
1798  for (; local_vector_begin != local_vector_end;
1799  ++local_vector_begin, ++local_indices_begin)
1800  {
1801  if (is_constrained(*local_indices_begin) == false)
1802  internal::ElementAccess<VectorType>::add(*local_vector_begin,
1803  *local_indices_begin,
1804  global_vector);
1805  else
1806  {
1807  const ConstraintLine &position =
1808  lines[lines_cache[calculate_line_index(*local_indices_begin)]];
1809  for (size_type j = 0; j < position.entries.size(); ++j)
1810  internal::ElementAccess<VectorType>::add(
1811  (*local_vector_begin) * position.entries[j].second,
1812  position.entries[j].first,
1813  global_vector);
1814  }
1815  }
1816 }
1817 
1818 template <typename number>
1819 template <class InVector, class OutVector>
1820 inline void
1822  const InVector & local_vector,
1823  const std::vector<size_type> &local_dof_indices,
1824  OutVector & global_vector) const
1825 {
1826  Assert(local_vector.size() == local_dof_indices.size(),
1827  ExcDimensionMismatch(local_vector.size(), local_dof_indices.size()));
1828  distribute_local_to_global(local_vector.begin(),
1829  local_vector.end(),
1830  local_dof_indices.begin(),
1831  global_vector);
1832 }
1833 
1834 template <typename number>
1835 template <typename ForwardIteratorVec,
1836  typename ForwardIteratorInd,
1837  class VectorType>
1838 inline void
1840  const VectorType & global_vector,
1841  ForwardIteratorInd local_indices_begin,
1842  ForwardIteratorVec local_vector_begin,
1843  ForwardIteratorVec local_vector_end) const
1844 {
1845  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
1846  for (; local_vector_begin != local_vector_end;
1847  ++local_vector_begin, ++local_indices_begin)
1848  {
1849  if (is_constrained(*local_indices_begin) == false)
1850  *local_vector_begin = global_vector(*local_indices_begin);
1851  else
1852  {
1853  const ConstraintLine &position =
1854  lines[lines_cache[calculate_line_index(*local_indices_begin)]];
1855  typename VectorType::value_type value = position.inhomogeneity;
1856  for (size_type j = 0; j < position.entries.size(); ++j)
1857  value += (global_vector(position.entries[j].first) *
1858  position.entries[j].second);
1859  *local_vector_begin = value;
1860  }
1861  }
1862 }
1863 
1864 template <typename MatrixType>
1866 template <typename SparsityPatternType>
1868 template <typename number>
1870 
1889 template <typename MatrixType>
1891 {
1892 private:
1893  struct yes_type
1894  {
1895  char c[1];
1896  };
1897  struct no_type
1898  {
1899  char c[2];
1900  };
1901 
1907  template <typename T>
1908  static yes_type
1909  check_for_block_matrix(const BlockMatrixBase<T> *);
1910 
1915  template <typename T>
1916  static yes_type
1917  check_for_block_matrix(const BlockSparsityPatternBase<T> *);
1918 
1923  template <typename T>
1924  static yes_type
1925  check_for_block_matrix(const BlockSparseMatrixEZ<T> *);
1926 
1931  static no_type
1932  check_for_block_matrix(...);
1933 
1934 public:
1940  static const bool value =
1941  (sizeof(check_for_block_matrix(static_cast<MatrixType *>(nullptr))) ==
1942  sizeof(yes_type));
1943 };
1944 
1945 // instantiation of the static member
1946 template <typename MatrixType>
1948 
1949 template <typename number>
1950 template <typename MatrixType>
1951 inline void
1953  const FullMatrix<number> & local_matrix,
1954  const std::vector<size_type> &local_dof_indices,
1955  MatrixType & global_matrix) const
1956 {
1957  // create a dummy and hand on to the function actually implementing this
1958  // feature in the cm.templates.h file.
1959  Vector<typename MatrixType::value_type> dummy(0);
1961  local_matrix,
1962  dummy,
1963  local_dof_indices,
1964  global_matrix,
1965  dummy,
1966  false,
1967  std::integral_constant<bool, IsBlockMatrix<MatrixType>::value>());
1968 }
1969 
1970 template <typename number>
1971 template <typename MatrixType, typename VectorType>
1972 inline void
1974  const FullMatrix<number> & local_matrix,
1975  const Vector<number> & local_vector,
1976  const std::vector<size_type> &local_dof_indices,
1977  MatrixType & global_matrix,
1978  VectorType & global_vector,
1979  bool use_inhomogeneities_for_rhs) const
1980 {
1981  // enter the internal function with the respective block information set,
1982  // the actual implementation follows in the cm.templates.h file.
1984  local_matrix,
1985  local_vector,
1986  local_dof_indices,
1987  global_matrix,
1988  global_vector,
1989  use_inhomogeneities_for_rhs,
1990  std::integral_constant<bool, IsBlockMatrix<MatrixType>::value>());
1991 }
1992 
1993 template <typename number>
1994 template <typename SparsityPatternType>
1995 inline void
1997  const std::vector<size_type> &local_dof_indices,
1998  SparsityPatternType & sparsity_pattern,
1999  const bool keep_constrained_entries,
2000  const Table<2, bool> & dof_mask) const
2001 {
2002  // enter the internal function with the respective block information set,
2003  // the actual implementation follows in the cm.templates.h file.
2005  local_dof_indices,
2006  sparsity_pattern,
2007  keep_constrained_entries,
2008  dof_mask,
2009  std::integral_constant<bool, IsBlockMatrix<SparsityPatternType>::value>());
2010 }
2011 
2012 DEAL_II_NAMESPACE_CLOSE
2013 
2014 #endif
const types::global_dof_index invalid_size_type
Definition: types.h:182
std::vector< ConstraintLine > lines
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)
TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
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:1600
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:1407
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:1821
#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:1608
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:587
Definition: table.h:37
bool is_element(const size_type index) const
Definition: index_set.h:1665
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()