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