Reference documentation for deal.II version Git 24fafe1087 2021-11-29 17:11:51 +0100
\(\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 <algorithm>
34 #include <set>
35 #include <type_traits>
36 #include <utility>
37 #include <vector>
38 
40 
41 // Forward declarations
42 #ifndef DOXYGEN
43 template <typename>
44 class FullMatrix;
45 class SparsityPattern;
49 template <typename number>
50 class SparseMatrix;
51 template <typename number>
52 class BlockSparseMatrix;
53 
54 namespace internal
55 {
56  namespace AffineConstraints
57  {
59 
77  struct Distributing
78  {
79  Distributing(const size_type global_row = numbers::invalid_size_type,
80  const size_type local_row = numbers::invalid_size_type);
81 
82  Distributing(const Distributing &in);
83 
84  Distributing &
85  operator=(const Distributing &in);
86 
87  bool
88  operator<(const Distributing &in) const
89  {
90  return global_row < in.global_row;
91  }
92 
93  size_type global_row;
94  size_type local_row;
95  mutable size_type constraint_position;
96  };
97 
98 
99 
111  template <typename number>
112  struct DataCache
113  {
114  DataCache();
115 
116  void
117  reinit();
118 
119  size_type
120  insert_new_index(const std::pair<size_type, number> &pair);
121 
122  void
123  append_index(const size_type index,
124  const std::pair<size_type, number> &pair);
125 
126  size_type
127  get_size(const size_type index) const;
128 
129  const std::pair<size_type, number> *
130  get_entry(const size_type index) const;
131 
132  size_type row_length;
133 
134  std::vector<std::pair<size_type, number>> data;
135 
136  std::vector<size_type> individual_size;
137  };
138 
139 
140 
169  template <typename number>
170  class GlobalRowsFromLocal
171  {
172  public:
176  GlobalRowsFromLocal();
177 
178  void
179  reinit(const size_type n_local_rows);
180 
181  void
182  insert_index(const size_type global_row,
183  const size_type local_row,
184  const number constraint_value);
185  void
186  sort();
187 
188  void
189  print(std::ostream &os);
190 
194  size_type
195  size() const;
196 
201  size_type
202  size(const size_type counter_index) const;
203 
207  size_type
208  global_row(const size_type counter_index) const;
209 
213  size_type &
214  global_row(const size_type counter_index);
215 
221  size_type
222  local_row(const size_type counter_index) const;
223 
227  size_type &
228  local_row(const size_type counter_index);
229 
235  size_type
236  local_row(const size_type counter_index,
237  const size_type index_in_constraint) const;
238 
243  number
244  constraint_value(const size_type counter_index,
245  const size_type index_in_constraint) const;
246 
252  bool
253  have_indirect_rows() const;
254 
259  void
260  insert_constraint(const size_type constrained_local_dof);
261 
267  size_type
268  n_constraints() const;
269 
274  size_type
275  n_inhomogeneities() const;
276 
283  void
284  set_ith_constraint_inhomogeneous(const size_type i);
285 
290  size_type
291  constraint_origin(size_type i) const;
292 
298  std::vector<Distributing> total_row_indices;
299 
300  private:
304  DataCache<number> data_cache;
305 
310  size_type n_active_rows;
311 
316  size_type n_inhomogeneous_rows;
317  };
318 
319 
320 
338  template <typename number>
339  struct ScratchData
340  {
344  ScratchData()
345  : in_use(false)
346  {}
347 
351  ScratchData(const ScratchData &)
352  : in_use(false)
353  {}
354 
358  bool in_use;
359 
363  std::vector<size_type> columns;
364 
368  std::vector<number> values;
369 
373  std::vector<size_type> block_starts;
374 
378  std::vector<size_type> vector_indices;
379 
383  std::vector<number> vector_values;
384 
388  GlobalRowsFromLocal<number> global_rows;
389 
393  GlobalRowsFromLocal<number> global_columns;
394  };
395  } // namespace AffineConstraints
396 } // namespace internal
397 
398 namespace internal
399 {
400  namespace AffineConstraintsImplementation
401  {
402  template <class VectorType>
403  void
404  set_zero_all(const std::vector<types::global_dof_index> &cm,
405  VectorType & vec);
406 
407  template <class T>
408  void
409  set_zero_all(const std::vector<types::global_dof_index> &cm,
410  ::Vector<T> & vec);
411 
412  template <class T>
413  void
414  set_zero_all(const std::vector<types::global_dof_index> &cm,
415  ::BlockVector<T> & vec);
416  } // namespace AffineConstraintsImplementation
417 } // namespace internal
418 
419 
420 template <typename number>
421 class AffineConstraints;
422 #endif
423 
424 // TODO[WB]: We should have a function of the kind
425 // AffineConstraints::add_constraint (const size_type constrained_dof,
426 // const std::vector<std::pair<size_type, number> > &entries,
427 // const number inhomogeneity = 0);
428 // rather than building up constraints piecemeal through add_line/add_entry
429 // etc. This would also eliminate the possibility of accidentally changing
430 // existing constraints into something pointless, see the discussion on the
431 // mailing list on "Tiny bug in interpolate_boundary_values" in Sept. 2010.
432 
506 template <typename number = double>
508 {
509 public:
514 
521  {
527 
534 
540  right_object_wins
541  };
542 
560  explicit AffineConstraints(const IndexSet &local_constraints = IndexSet());
561 
565  explicit AffineConstraints(const AffineConstraints &affine_constraints);
566 
570  AffineConstraints(AffineConstraints &&affine_constraints) noexcept =
571  default; // NOLINT
572 
583  operator=(const AffineConstraints &) = delete;
584 
589  operator=(AffineConstraints &&affine_constraints) noexcept =
590  default; // NOLINT
591 
598  template <typename other_number>
599  void
600  copy_from(const AffineConstraints<other_number> &other);
601 
608  void
609  reinit(const IndexSet &local_constraints = IndexSet());
610 
617  bool
618  can_store_line(const size_type line_n) const;
619 
626  const IndexSet &
627  get_local_lines() const;
628 
648  void
649  add_selected_constraints(const AffineConstraints &constraints_in,
650  const IndexSet & filter);
651 
661  void
662  add_line(const size_type line_n);
663 
676  void
677  add_lines(const std::vector<bool> &lines);
678 
691  void
692  add_lines(const std::set<size_type> &lines);
693 
706  void
707  add_lines(const IndexSet &lines);
708 
724  void
725  add_entry(const size_type constrained_dof_index,
726  const size_type column,
727  const number weight);
728 
734  void
735  add_entries(
736  const size_type constrained_dof_index,
737  const std::vector<std::pair<size_type, number>> &col_weight_pairs);
738 
750  void
751  set_inhomogeneity(const size_type constrained_dof_index, const number value);
752 
774  void
775  close();
776 
801  void
802  merge(
803  const AffineConstraints & other_constraints,
804  const MergeConflictBehavior merge_conflict_behavior = no_conflicts_allowed,
805  const bool allow_different_local_lines = false);
806 
819  void
820  shift(const size_type offset);
821 
829  void
830  clear();
831 
844  size_type
845  n_constraints() const;
846 
851  size_type
852  n_identities() const;
853 
859  size_type
860  n_inhomogeneities() const;
861 
871  bool
872  is_constrained(const size_type line_n) const;
873 
885  bool
886  is_identity_constrained(const size_type line_n) const;
887 
894  bool
895  are_identity_constrained(const size_type line_n_1,
896  const size_type line_n_2) const;
897 
908  size_type
909  max_constraint_indirections() const;
910 
915  bool
916  is_inhomogeneously_constrained(const size_type index) const;
917 
923  bool
924  has_inhomogeneities() const;
925 
930  const std::vector<std::pair<size_type, number>> *
931  get_constraint_entries(const size_type line_n) const;
932 
937  number
938  get_inhomogeneity(const size_type line_n) const;
939 
960  void
961  print(std::ostream &out) const;
962 
975  void
976  write_dot(std::ostream &) const;
977 
982  std::size_t
983  memory_consumption() const;
984 
991  void
992  resolve_indices(std::vector<types::global_dof_index> &indices) const;
993 
1019  void
1020  condense(SparsityPattern &sparsity) const;
1021 
1025  void
1026  condense(BlockSparsityPattern &sparsity) const;
1027 
1032  void
1033  condense(DynamicSparsityPattern &sparsity) const;
1034 
1039  void
1040  condense(BlockDynamicSparsityPattern &sparsity) const;
1041 
1049  void
1050  condense(SparseMatrix<number> &matrix) const;
1051 
1055  void
1056  condense(BlockSparseMatrix<number> &matrix) const;
1057 
1069  template <class VectorType>
1070  void
1071  condense(VectorType &vec) const;
1072 
1079  template <class VectorType>
1080  void
1081  condense(const VectorType &vec_ghosted, VectorType &output) const;
1082 
1095  template <class VectorType>
1096  void
1097  condense(SparseMatrix<number> &matrix, VectorType &vector) const;
1098 
1103  template <class BlockVectorType>
1104  void
1105  condense(BlockSparseMatrix<number> &matrix, BlockVectorType &vector) const;
1106 
1113  template <class VectorType>
1114  void
1115  set_zero(VectorType &vec) const;
1116 
1172  template <class InVector, class OutVector>
1173  void
1174  distribute_local_to_global(const InVector & local_vector,
1175  const std::vector<size_type> &local_dof_indices,
1176  OutVector & global_vector) const;
1177 
1225  template <typename VectorType>
1226  void
1227  distribute_local_to_global(const Vector<number> & local_vector,
1228  const std::vector<size_type> &local_dof_indices,
1229  VectorType & global_vector,
1230  const FullMatrix<number> & local_matrix) const;
1231 
1251  template <typename VectorType>
1252  void
1253  distribute_local_to_global(
1254  const Vector<number> & local_vector,
1255  const std::vector<size_type> &local_dof_indices_row,
1256  const std::vector<size_type> &local_dof_indices_col,
1257  VectorType & global_vector,
1258  const FullMatrix<number> & local_matrix,
1259  bool diagonal = false) const;
1260 
1264  template <class VectorType>
1265  void
1266  distribute_local_to_global(const size_type index,
1267  const number value,
1268  VectorType & global_vector) const;
1269 
1302  template <typename ForwardIteratorVec,
1303  typename ForwardIteratorInd,
1304  class VectorType>
1305  void
1306  distribute_local_to_global(ForwardIteratorVec local_vector_begin,
1307  ForwardIteratorVec local_vector_end,
1308  ForwardIteratorInd local_indices_begin,
1309  VectorType & global_vector) const;
1310 
1362  template <typename MatrixType>
1363  void
1364  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1365  const std::vector<size_type> &local_dof_indices,
1366  MatrixType & global_matrix) const;
1367 
1395  template <typename MatrixType>
1396  void
1397  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1398  const std::vector<size_type> &row_indices,
1399  const std::vector<size_type> &col_indices,
1400  MatrixType & global_matrix) const;
1401 
1418  template <typename MatrixType>
1419  void
1420  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1421  const std::vector<size_type> &row_indices,
1422  const AffineConstraints &column_affine_constraints,
1423  const std::vector<size_type> &column_indices,
1424  MatrixType & global_matrix) const;
1425 
1446  template <typename MatrixType, typename VectorType>
1447  void
1448  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1449  const Vector<number> & local_vector,
1450  const std::vector<size_type> &local_dof_indices,
1451  MatrixType & global_matrix,
1452  VectorType & global_vector,
1453  bool use_inhomogeneities_for_rhs = false) const;
1454 
1508  template <typename SparsityPatternType>
1509  void
1510  add_entries_local_to_global(
1511  const std::vector<size_type> &local_dof_indices,
1512  SparsityPatternType & sparsity_pattern,
1513  const bool keep_constrained_entries = true,
1514  const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1515 
1519  template <typename SparsityPatternType>
1520  void
1521  add_entries_local_to_global(
1522  const std::vector<size_type> &row_indices,
1523  const std::vector<size_type> &col_indices,
1524  SparsityPatternType & sparsity_pattern,
1525  const bool keep_constrained_entries = true,
1526  const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1527 
1547  template <typename ForwardIteratorVec,
1548  typename ForwardIteratorInd,
1549  class VectorType>
1550  void
1551  get_dof_values(const VectorType & global_vector,
1552  ForwardIteratorInd local_indices_begin,
1553  ForwardIteratorVec local_vector_begin,
1554  ForwardIteratorVec local_vector_end) const;
1555 
1577  template <class VectorType>
1578  void
1579  distribute(VectorType &vec) const;
1580 
1589  {
1594  using Entries = std::vector<std::pair<size_type, number>>;
1595 
1602 
1611 
1616 
1621  const Entries & entries = {},
1622  const number & inhomogeneity = 0.0);
1623 
1627  template <typename ConstraintLineType>
1628  ConstraintLine(const ConstraintLineType &other);
1629 
1633  template <typename ConstraintLineType>
1634  ConstraintLine &
1635  operator=(const ConstraintLineType &other);
1636 
1644  bool
1645  operator<(const ConstraintLine &) const;
1646 
1652  bool
1653  operator==(const ConstraintLine &) const;
1654 
1659  std::size_t
1660  memory_consumption() const;
1661 
1667  template <class Archive>
1668  void
1669  serialize(Archive &ar, const unsigned int)
1670  {
1671  ar &index &entries &inhomogeneity;
1672  }
1673  };
1674 
1678  using const_iterator = typename std::vector<ConstraintLine>::const_iterator;
1679 
1683  using LineRange = boost::iterator_range<const_iterator>;
1684 
1693  const LineRange
1694  get_lines() const;
1695 
1724  bool
1725  is_consistent_in_parallel(const std::vector<IndexSet> &locally_owned_dofs,
1726  const IndexSet & locally_active_dofs,
1727  const MPI_Comm & mpi_communicator,
1728  const bool verbose = false) const;
1729 
1735  void
1736  make_consistent_in_parallel(const IndexSet &locally_owned_dofs,
1737  const IndexSet &locally_relevant_dofs,
1738  const MPI_Comm mpi_communicator);
1739 
1745  DeclException0(ExcMatrixIsClosed);
1751  DeclException0(ExcMatrixNotClosed);
1758  size_type,
1759  << "The specified line " << arg1 << " does not exist.");
1765  DeclException4(ExcEntryAlreadyExists,
1766  size_type,
1767  size_type,
1768  number,
1769  number,
1770  << "The entry for the indices " << arg1 << " and " << arg2
1771  << " already exists, but the values " << arg3 << " (old) and "
1772  << arg4 << " (new) differ "
1773  << "by " << (arg4 - arg3) << ".");
1779  DeclException2(ExcDoFConstrainedToConstrainedDoF,
1780  int,
1781  int,
1782  << "You tried to constrain DoF " << arg1 << " to DoF " << arg2
1783  << ", but that one is also constrained. This is not allowed!");
1789  DeclException1(ExcDoFIsConstrainedFromBothObjects,
1790  size_type,
1791  << "Degree of freedom " << arg1
1792  << " is constrained from both object in a merge operation.");
1798  DeclException1(ExcDoFIsConstrainedToConstrainedDoF,
1799  size_type,
1800  << "In the given argument a degree of freedom is constrained "
1801  << "to another DoF with number " << arg1
1802  << ", which however is constrained by this object. This is not"
1803  << " allowed.");
1809  DeclException1(ExcRowNotStoredHere,
1810  size_type,
1811  << "The index set given to this constraints object indicates "
1812  << "constraints for degree of freedom " << arg1
1813  << " should not be stored by this object, but a constraint "
1814  << "is being added.");
1815 
1821  DeclException2(ExcColumnNotStoredHere,
1822  size_type,
1823  size_type,
1824  << "The index set given to this constraints object indicates "
1825  << "constraints using degree of freedom " << arg2
1826  << " should not be stored by this object, but a constraint "
1827  << "for degree of freedom " << arg1 << " uses it.");
1828 
1834  DeclException2(ExcIncorrectConstraint,
1835  int,
1836  int,
1837  << "While distributing the constraint for DoF " << arg1
1838  << ", it turns out that one of the processors "
1839  << "who own the " << arg2 << " degrees of freedom that x_"
1840  << arg1 << " is constrained against does not know about "
1841  << "the constraint on x_" << arg1
1842  << ". Did you not initialize the AffineConstraints container "
1843  << "with the appropriate locally_relevant set so "
1844  << "that every processor who owns a DoF that constrains "
1845  << "another DoF also knows about this constraint?");
1846 
1847  template <typename>
1848  friend class AffineConstraints;
1849 
1850 private:
1862  std::vector<ConstraintLine> lines;
1863 
1896  std::vector<size_type> lines_cache;
1897 
1904 
1908  bool sorted;
1909 
1911  internal::AffineConstraints::ScratchData<number>>
1913 
1918  size_type
1919  calculate_line_index(const size_type line_n) const;
1920 
1925  template <typename MatrixType, typename VectorType>
1926  void
1927  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1928  const Vector<number> & local_vector,
1929  const std::vector<size_type> &local_dof_indices,
1930  MatrixType & global_matrix,
1931  VectorType & global_vector,
1932  const bool use_inhomogeneities_for_rhs,
1933  const std::integral_constant<bool, false>) const;
1934 
1939  template <typename MatrixType, typename VectorType>
1940  void
1941  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1942  const Vector<number> & local_vector,
1943  const std::vector<size_type> &local_dof_indices,
1944  MatrixType & global_matrix,
1945  VectorType & global_vector,
1946  const bool use_inhomogeneities_for_rhs,
1947  const std::integral_constant<bool, true>) const;
1948 
1953  template <typename SparsityPatternType>
1954  void
1955  add_entries_local_to_global(const std::vector<size_type> &local_dof_indices,
1956  SparsityPatternType & sparsity_pattern,
1957  const bool keep_constrained_entries,
1958  const Table<2, bool> &dof_mask,
1959  const std::integral_constant<bool, false>) const;
1960 
1965  template <typename SparsityPatternType>
1966  void
1967  add_entries_local_to_global(const std::vector<size_type> &local_dof_indices,
1968  SparsityPatternType & sparsity_pattern,
1969  const bool keep_constrained_entries,
1970  const Table<2, bool> &dof_mask,
1971  const std::integral_constant<bool, true>) const;
1972 
1980  void
1981  make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
1982  internal::AffineConstraints::GlobalRowsFromLocal<number>
1983  &global_rows) const;
1984 
1992  void
1993  make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
1994  std::vector<size_type> & active_dofs) const;
1995 
1999  template <typename MatrixScalar, typename VectorScalar>
2002  const size_type i,
2003  const internal::AffineConstraints::GlobalRowsFromLocal<number> &global_rows,
2004  const Vector<VectorScalar> & local_vector,
2005  const std::vector<size_type> & local_dof_indices,
2006  const FullMatrix<MatrixScalar> &local_matrix) const;
2007 };
2008 
2009 /* ---------------- template and inline functions ----------------- */
2010 
2011 template <typename number>
2013  const IndexSet &local_constraints)
2014  : lines()
2015  , local_lines(local_constraints)
2016  , sorted(false)
2017 {
2018  // make sure the IndexSet is compressed. Otherwise this can lead to crashes
2019  // that are hard to find (only happen in release mode).
2020  // see tests/mpi/affine_constraints_crash_01
2022 }
2023 
2024 template <typename number>
2026  const AffineConstraints &affine_constraints)
2027  : Subscriptor()
2028  , lines(affine_constraints.lines)
2029  , lines_cache(affine_constraints.lines_cache)
2030  , local_lines(affine_constraints.local_lines)
2031  , sorted(affine_constraints.sorted)
2032 {}
2033 
2034 template <typename number>
2035 inline void
2037 {
2038  Assert(sorted == false, ExcMatrixIsClosed());
2039 
2040  // the following can happen when we compute with distributed meshes and dof
2041  // handlers and we constrain a degree of freedom whose number we don't have
2042  // locally. if we don't abort here the program will try to allocate several
2043  // terabytes of memory to resize the various arrays below :-)
2045  const size_type line_index = calculate_line_index(line_n);
2046 
2047  // check whether line already exists; it may, in which case we can just quit
2048  if (is_constrained(line_n))
2049  return;
2050 
2051  // if necessary enlarge vector of existing entries for cache
2052  if (line_index >= lines_cache.size())
2053  lines_cache.resize(std::max(2 * static_cast<size_type>(lines_cache.size()),
2054  line_index + 1),
2056 
2057  // push a new line to the end of the list
2058  lines.emplace_back();
2059  lines.back().index = line_n;
2060  lines.back().inhomogeneity = 0.;
2061  lines_cache[line_index] = lines.size() - 1;
2062 }
2063 
2064 
2065 
2066 template <typename number>
2067 inline void
2069  const size_type column,
2070  const number weight)
2071 {
2072  Assert(sorted == false, ExcMatrixIsClosed());
2073  Assert(constrained_dof_index != column,
2074  ExcMessage("Can't constrain a degree of freedom to itself"));
2075 
2076  // Ensure that the current line is present in the cache:
2077  const size_type line_index = calculate_line_index(constrained_dof_index);
2078  Assert(line_index < lines_cache.size(),
2079  ExcMessage("The current AffineConstraints does not contain the line "
2080  "for the current entry. Call AffineConstraints::add_line "
2081  "before calling this function."));
2082 
2083  // if in debug mode, check whether an entry for this column already exists
2084  // and if it's the same as the one entered at present
2085  //
2086  // in any case: exit the function if an entry for this column already
2087  // exists, since we don't want to enter it twice
2089  ExcInternalError());
2091  ExcColumnNotStoredHere(constrained_dof_index, column));
2092  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
2093  Assert(line_ptr->index == constrained_dof_index, ExcInternalError());
2094  for (const auto &p : line_ptr->entries)
2095  if (p.first == column)
2096  {
2097  Assert(std::abs(p.second - weight) < 1.e-14,
2099  constrained_dof_index, column, p.second, weight));
2100  return;
2101  }
2102 
2103  line_ptr->entries.emplace_back(column, weight);
2104 }
2105 
2106 
2107 
2108 template <typename number>
2109 inline void
2111  const size_type constrained_dof_index,
2112  const number value)
2113 {
2114  const size_type line_index = calculate_line_index(constrained_dof_index);
2115  Assert(line_index < lines_cache.size() &&
2116  lines_cache[line_index] != numbers::invalid_size_type,
2117  ExcMessage("call add_line() before calling set_inhomogeneity()"));
2118  Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
2119  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
2120  line_ptr->inhomogeneity = value;
2121 }
2122 
2123 
2124 
2125 template <typename number>
2126 template <class VectorType>
2127 inline void
2129 {
2130  // since lines is a private member, we cannot pass it to the functions
2131  // above. therefore, copy the content which is cheap
2132  std::vector<size_type> constrained_lines(lines.size());
2133  for (unsigned int i = 0; i < lines.size(); ++i)
2134  constrained_lines[i] = lines[i].index;
2135  internal::AffineConstraintsImplementation::set_zero_all(constrained_lines,
2136  vec);
2137 }
2138 
2139 template <typename number>
2142 {
2143  return lines.size();
2144 }
2145 
2146 template <typename number>
2149 {
2150  return std::count_if(lines.begin(),
2151  lines.end(),
2152  [](const ConstraintLine &line) {
2153  return (line.entries.size() == 1) &&
2154  (line.entries[0].second == number(1.));
2155  });
2156 }
2157 
2158 template <typename number>
2161 {
2162  return std::count_if(lines.begin(),
2163  lines.end(),
2164  [](const ConstraintLine &line) {
2165  return (line.inhomogeneity != number(0.));
2166  });
2167 }
2168 
2169 template <typename number>
2170 inline bool
2172 {
2173  const size_type line_index = calculate_line_index(index);
2174  return ((line_index < lines_cache.size()) &&
2175  (lines_cache[line_index] != numbers::invalid_size_type));
2176 }
2177 
2178 template <typename number>
2179 inline bool
2181  const size_type line_n) const
2182 {
2183  // check whether the entry is constrained. could use is_constrained, but
2184  // that means computing the line index twice
2185  const size_type line_index = calculate_line_index(line_n);
2186  if (line_index >= lines_cache.size() ||
2187  lines_cache[line_index] == numbers::invalid_size_type)
2188  return false;
2189  else
2190  {
2191  Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
2192  return (lines[lines_cache[line_index]].inhomogeneity != number(0.));
2193  }
2194 }
2195 
2196 template <typename number>
2197 inline const std::vector<std::pair<types::global_dof_index, number>> *
2199 {
2200  // check whether the entry is constrained. could use is_constrained, but
2201  // that means computing the line index twice
2202  const size_type line_index = calculate_line_index(line_n);
2203  if (line_index >= lines_cache.size() ||
2204  lines_cache[line_index] == numbers::invalid_size_type)
2205  return nullptr;
2206  else
2207  return &lines[lines_cache[line_index]].entries;
2208 }
2209 
2210 template <typename number>
2211 inline number
2213 {
2214  // check whether the entry is constrained. could use is_constrained, but
2215  // that means computing the line index twice
2216  const size_type line_index = calculate_line_index(line_n);
2217  if (line_index >= lines_cache.size() ||
2218  lines_cache[line_index] == numbers::invalid_size_type)
2219  return 0;
2220  else
2221  return lines[lines_cache[line_index]].inhomogeneity;
2222 }
2223 
2224 template <typename number>
2227 {
2228  // IndexSet is unused (serial case)
2229  if (!local_lines.size())
2230  return line_n;
2231 
2233 
2234  return local_lines.index_within_set(line_n);
2235 }
2236 
2237 template <typename number>
2238 inline bool
2240 {
2241  return !local_lines.size() || local_lines.is_element(line_n);
2242 }
2243 
2244 template <typename number>
2245 inline const IndexSet &
2247 {
2248  return local_lines;
2249 }
2250 
2251 template <typename number>
2252 template <class VectorType>
2253 inline void
2255  const size_type index,
2256  const number value,
2257  VectorType & global_vector) const
2258 {
2259  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2260 
2261  if (is_constrained(index) == false)
2262  global_vector(index) += value;
2263  else
2264  {
2265  const ConstraintLine &position =
2267  for (size_type j = 0; j < position.entries.size(); ++j)
2268  global_vector(position.entries[j].first) +=
2269  value * position.entries[j].second;
2270  }
2271 }
2272 
2273 template <typename number>
2274 template <typename ForwardIteratorVec,
2275  typename ForwardIteratorInd,
2276  class VectorType>
2277 inline void
2279  ForwardIteratorVec local_vector_begin,
2280  ForwardIteratorVec local_vector_end,
2281  ForwardIteratorInd local_indices_begin,
2282  VectorType & global_vector) const
2283 {
2284  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2285  for (; local_vector_begin != local_vector_end;
2286  ++local_vector_begin, ++local_indices_begin)
2287  {
2288  if (is_constrained(*local_indices_begin) == false)
2289  internal::ElementAccess<VectorType>::add(*local_vector_begin,
2290  *local_indices_begin,
2291  global_vector);
2292  else
2293  {
2294  const ConstraintLine &position =
2295  lines[lines_cache[calculate_line_index(*local_indices_begin)]];
2296  for (size_type j = 0; j < position.entries.size(); ++j)
2298  (*local_vector_begin) * position.entries[j].second,
2299  position.entries[j].first,
2300  global_vector);
2301  }
2302  }
2303 }
2304 
2305 template <typename number>
2306 template <class InVector, class OutVector>
2307 inline void
2309  const InVector & local_vector,
2310  const std::vector<size_type> &local_dof_indices,
2311  OutVector & global_vector) const
2312 {
2313  Assert(local_vector.size() == local_dof_indices.size(),
2314  ExcDimensionMismatch(local_vector.size(), local_dof_indices.size()));
2315  distribute_local_to_global(local_vector.begin(),
2316  local_vector.end(),
2317  local_dof_indices.begin(),
2318  global_vector);
2319 }
2320 
2321 template <typename number>
2322 template <typename ForwardIteratorVec,
2323  typename ForwardIteratorInd,
2324  class VectorType>
2325 inline void
2327  const VectorType & global_vector,
2328  ForwardIteratorInd local_indices_begin,
2329  ForwardIteratorVec local_vector_begin,
2330  ForwardIteratorVec local_vector_end) const
2331 {
2332  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2333  for (; local_vector_begin != local_vector_end;
2334  ++local_vector_begin, ++local_indices_begin)
2335  {
2336  if (is_constrained(*local_indices_begin) == false)
2337  *local_vector_begin = global_vector(*local_indices_begin);
2338  else
2339  {
2340  const ConstraintLine &position =
2341  lines[lines_cache[calculate_line_index(*local_indices_begin)]];
2342  typename VectorType::value_type value = position.inhomogeneity;
2343  for (size_type j = 0; j < position.entries.size(); ++j)
2344  value += (global_vector(position.entries[j].first) *
2345  position.entries[j].second);
2346  *local_vector_begin = value;
2347  }
2348  }
2349 }
2350 
2351 template <typename MatrixType>
2353 template <typename SparsityPatternType>
2355 template <typename number>
2357 
2358 namespace internal
2359 {
2361  {
2379  template <typename MatrixType>
2381  {
2382  private:
2388  template <typename T>
2389  static std::true_type
2390  check(const BlockMatrixBase<T> *);
2391 
2396  template <typename T>
2397  static std::true_type
2398  check(const BlockSparseMatrixEZ<T> *);
2399 
2404  static std::false_type
2405  check(...);
2406 
2407  public:
2414  static const bool value =
2415  std::is_same<decltype(check(std::declval<MatrixType *>())),
2416  std::true_type>::value;
2417  };
2418 
2419  // instantiation of the static member
2420  template <typename MatrixType>
2422 
2423 
2432  template <typename MatrixType>
2434  {
2435  private:
2440  template <typename T>
2441  static std::true_type
2443 
2448  static std::false_type
2449  check(...);
2450 
2451  public:
2457  static const bool value =
2458  std::is_same<decltype(check(std::declval<MatrixType *>())),
2459  std::true_type>::value;
2460  };
2461 
2462  // instantiation of the static member
2463  template <typename MatrixType>
2465 
2466  } // namespace AffineConstraints
2467 } // namespace internal
2468 
2469 
2470 
2471 template <typename number>
2472 template <typename other_number>
2473 inline void
2475  const AffineConstraints<other_number> &other)
2476 {
2477  lines.clear();
2478  lines.insert(lines.begin(), other.lines.begin(), other.lines.end());
2479  lines_cache = other.lines_cache;
2480  local_lines = other.local_lines;
2481  sorted = other.sorted;
2482 }
2483 
2484 
2485 
2486 template <typename number>
2487 template <typename MatrixType>
2488 inline void
2490  const FullMatrix<number> & local_matrix,
2491  const std::vector<size_type> &local_dof_indices,
2492  MatrixType & global_matrix) const
2493 {
2494  // create a dummy and hand on to the function actually implementing this
2495  // feature in the cm.templates.h file.
2496  Vector<typename MatrixType::value_type> dummy(0);
2498  local_matrix,
2499  dummy,
2500  local_dof_indices,
2501  global_matrix,
2502  dummy,
2503  false,
2504  std::integral_constant<
2505  bool,
2507 }
2508 
2509 
2510 
2511 template <typename number>
2512 template <typename MatrixType, typename VectorType>
2513 inline void
2515  const FullMatrix<number> & local_matrix,
2516  const Vector<number> & local_vector,
2517  const std::vector<size_type> &local_dof_indices,
2518  MatrixType & global_matrix,
2519  VectorType & global_vector,
2520  bool use_inhomogeneities_for_rhs) const
2521 {
2522  // enter the internal function with the respective block information set,
2523  // the actual implementation follows in the cm.templates.h file.
2525  local_matrix,
2526  local_vector,
2527  local_dof_indices,
2528  global_matrix,
2529  global_vector,
2530  use_inhomogeneities_for_rhs,
2531  std::integral_constant<
2532  bool,
2534 }
2535 
2536 
2537 
2538 template <typename number>
2539 template <typename SparsityPatternType>
2540 inline void
2542  const std::vector<size_type> &local_dof_indices,
2543  SparsityPatternType & sparsity_pattern,
2544  const bool keep_constrained_entries,
2545  const Table<2, bool> & dof_mask) const
2546 {
2547  // enter the internal function with the respective block information set,
2548  // the actual implementation follows in the cm.templates.h file.
2550  local_dof_indices,
2551  sparsity_pattern,
2552  keep_constrained_entries,
2553  dof_mask,
2554  std::integral_constant<bool,
2556  SparsityPatternType>::value>());
2557 }
2558 
2559 
2560 
2561 template <typename number>
2563  const size_type & index,
2565  const number &inhomogeneity)
2566  : index(index)
2567  , entries(entries)
2568  , inhomogeneity(inhomogeneity)
2569 {}
2570 
2571 
2572 
2573 template <typename number>
2574 template <typename ConstraintLineType>
2576  const ConstraintLineType &other)
2577 {
2578  this->index = other.index;
2579 
2580  entries.clear();
2581  entries.insert(entries.begin(), other.entries.begin(), other.entries.end());
2582 
2583  this->inhomogeneity = other.inhomogeneity;
2584 }
2585 
2586 
2587 
2588 template <typename number>
2589 template <typename ConstraintLineType>
2592  const ConstraintLineType &other)
2593 {
2594  this->index = other.index;
2595 
2596  entries.clear();
2597  entries.insert(entries.begin(), other.entries.begin(), other.entries.end());
2598 
2599  this->inhomogeneity = other.inhomogeneity;
2600 
2601  return *this;
2602 }
2603 
2605 
2606 #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:532
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:1639
Matrix is diagonal.
void serialize(Archive &ar, const unsigned int)
size_type n_constraints() const
static ::ExceptionBase & ExcMessage(std::string arg1)
size_type n_inhomogeneities() const
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
bool check(const ConstraintKinds &kind_in, const unsigned int dim)
#define Assert(cond, exc)
Definition: exceptions.h:1461
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:1926
size_type n_identities() const
#define DeclException0(Exception0)
Definition: exceptions.h:464
#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:1647
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:578
bool is_element(const size_type index) const
Definition: index_set.h:1770
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:2013
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()