Reference documentation for deal.II version GIT 29f9da0a34 2023-12-07 10:00:01+00:00
\(\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 - 2023 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 
29 #include <deal.II/lac/vector.h>
31 
32 #include <boost/range/iterator_range.hpp>
33 
34 #include <algorithm>
35 #include <set>
36 #include <type_traits>
37 #include <utility>
38 #include <vector>
39 
41 
42 // Forward declarations
43 #ifndef DOXYGEN
44 template <typename>
45 class FullMatrix;
46 class SparsityPattern;
50 template <typename number>
51 class SparseMatrix;
52 template <typename number>
53 class BlockSparseMatrix;
54 
55 namespace internal
56 {
57  namespace AffineConstraints
58  {
60 
78  struct Distributing
79  {
80  Distributing(const size_type global_row = numbers::invalid_size_type,
81  const size_type local_row = numbers::invalid_size_type);
82 
83  Distributing(const Distributing &in);
84 
85  Distributing &
86  operator=(const Distributing &in);
87 
88  bool
89  operator<(const Distributing &in) const
90  {
91  return global_row < in.global_row;
92  }
93 
94  size_type global_row;
95  size_type local_row;
96  mutable size_type constraint_position;
97  };
98 
99 
100 
112  template <typename number>
113  struct DataCache
114  {
115  DataCache();
116 
117  void
118  reinit();
119 
120  size_type
121  insert_new_index(const std::pair<size_type, number> &pair);
122 
123  void
124  append_index(const size_type index,
125  const std::pair<size_type, number> &pair);
126 
127  size_type
128  get_size(const size_type index) const;
129 
130  const std::pair<size_type, number> *
131  get_entry(const size_type index) const;
132 
133  size_type row_length;
134 
135  std::vector<std::pair<size_type, number>> data;
136 
137  std::vector<size_type> individual_size;
138  };
139 
140 
141 
170  template <typename number>
171  class GlobalRowsFromLocal
172  {
173  public:
177  GlobalRowsFromLocal();
178 
179  void
180  reinit(const size_type n_local_rows);
181 
182  void
183  insert_index(const size_type global_row,
184  const size_type local_row,
185  const number constraint_value);
186  void
187  sort();
188 
189  void
190  print(std::ostream &os);
191 
195  size_type
196  size() const;
197 
202  size_type
203  size(const size_type counter_index) const;
204 
208  size_type
209  global_row(const size_type counter_index) const;
210 
214  size_type &
215  global_row(const size_type counter_index);
216 
222  size_type
223  local_row(const size_type counter_index) const;
224 
228  size_type &
229  local_row(const size_type counter_index);
230 
236  size_type
237  local_row(const size_type counter_index,
238  const size_type index_in_constraint) const;
239 
244  number
245  constraint_value(const size_type counter_index,
246  const size_type index_in_constraint) const;
247 
254  bool
255  have_indirect_rows() const;
256 
261  void
262  insert_constraint(const size_type constrained_local_dof);
263 
269  size_type
270  n_constraints() const;
271 
276  size_type
277  n_inhomogeneities() const;
278 
285  void
286  set_ith_constraint_inhomogeneous(const size_type i);
287 
292  size_type
293  constraint_origin(size_type i) const;
294 
300  std::vector<Distributing> total_row_indices;
301 
302  private:
306  DataCache<number> data_cache;
307 
312  size_type n_active_rows;
313 
318  size_type n_inhomogeneous_rows;
319  };
320 
321 
322 
340  template <typename number>
341  struct ScratchData
342  {
346  ScratchData()
347  : in_use(false)
348  {}
349 
353  ScratchData(const ScratchData &)
354  : in_use(false)
355  {}
356 
360  bool in_use;
361 
365  std::vector<std::pair<size_type, size_type>> new_entries;
366 
370  std::vector<size_type> rows;
371 
375  std::vector<size_type> columns;
376 
380  std::vector<number> values;
381 
385  std::vector<size_type> block_starts;
386 
390  std::vector<size_type> vector_indices;
391 
395  std::vector<number> vector_values;
396 
400  GlobalRowsFromLocal<number> global_rows;
401 
405  GlobalRowsFromLocal<number> global_columns;
406  };
407  } // namespace AffineConstraints
408 } // namespace internal
409 
410 namespace internal
411 {
412  namespace AffineConstraintsImplementation
413  {
414  template <typename VectorType>
415  void
416  set_zero_all(const std::vector<types::global_dof_index> &cm,
417  VectorType &vec);
418 
419  template <class T>
420  void
421  set_zero_all(const std::vector<types::global_dof_index> &cm,
422  ::Vector<T> &vec);
423 
424  template <class T>
425  void
426  set_zero_all(const std::vector<types::global_dof_index> &cm,
427  ::BlockVector<T> &vec);
428  } // namespace AffineConstraintsImplementation
429 } // namespace internal
430 #endif
431 
432 
433 
506 template <typename number = double>
508 {
509 public:
514 
521  {
527 
534 
541  };
542 
561 
584  DEAL_II_DEPRECATED_EARLY_WITH_COMMENT(
585  "Use the constructor with two index set arguments.")
586  explicit AffineConstraints(const IndexSet &locally_stored_constraints);
587 
620  const IndexSet &locally_stored_constraints);
621 
625  explicit AffineConstraints(const AffineConstraints &affine_constraints);
626 
630  AffineConstraints(AffineConstraints &&affine_constraints) noexcept =
631  default; // NOLINT
632 
643  operator=(const AffineConstraints &) = delete;
644 
649  operator=(AffineConstraints &&affine_constraints) noexcept =
650  default; // NOLINT
651 
658  template <typename other_number>
659  void
660  copy_from(const AffineConstraints<other_number> &other);
661 
670  void
672 
681  DEAL_II_DEPRECATED_EARLY_WITH_COMMENT(
682  "Use the reinit() function with two index set arguments.")
683  void
684  reinit(const IndexSet &locally_stored_constraints);
685 
694  void
696  const IndexSet &locally_stored_constraints);
697 
706  bool
707  can_store_line(const size_type line_n) const;
708 
717  const IndexSet &
719 
726  const IndexSet &
728 
757  DEAL_II_DEPRECATED_EARLY_WITH_COMMENT("Use get_view() and merge() instead.")
758  void
760  const IndexSet &filter);
761 
812  void
814  const size_type constrained_dof,
815  const ArrayView<const std::pair<size_type, number>> &dependencies,
816  const number inhomogeneity = 0);
817 
835  void
836  constrain_dof_to_zero(const size_type constrained_dof);
837 
842  void
843  add_line(const size_type line_n);
844 
857  void
858  add_lines(const std::vector<bool> &lines);
859 
872  void
873  add_lines(const std::set<size_type> &lines);
874 
887  void
889 
905  void
906  add_entry(const size_type constrained_dof_index,
907  const size_type column,
908  const number weight);
909 
915  void
917  const size_type constrained_dof_index,
918  const std::vector<std::pair<size_type, number>> &col_weight_pairs);
919 
931  void
932  set_inhomogeneity(const size_type constrained_dof_index, const number value);
933 
955  void
956  close();
957 
963  bool
964  is_closed() const;
965 
971  bool
972  is_closed(const MPI_Comm comm) const;
973 
998  template <typename other_number>
999  void
1001  const AffineConstraints<other_number> &other_constraints,
1002  const MergeConflictBehavior merge_conflict_behavior = no_conflicts_allowed,
1003  const bool allow_different_local_lines = false);
1004 
1017  void
1018  shift(const size_type offset);
1019 
1069  get_view(const IndexSet &mask) const;
1070 
1078  void
1080 
1093  size_type
1094  n_constraints() const;
1095 
1100  size_type
1101  n_identities() const;
1102 
1108  size_type
1110 
1120  bool
1121  is_constrained(const size_type line_n) const;
1122 
1134  bool
1135  is_identity_constrained(const size_type line_n) const;
1136 
1143  bool
1145  const size_type line_n_2) const;
1146 
1157  size_type
1159 
1164  bool
1166 
1172  bool
1174 
1179  const std::vector<std::pair<size_type, number>> *
1180  get_constraint_entries(const size_type line_n) const;
1181 
1186  number
1187  get_inhomogeneity(const size_type line_n) const;
1188 
1209  void
1210  print(std::ostream &out) const;
1211 
1224  void
1225  write_dot(std::ostream &) const;
1226 
1231  std::size_t
1233 
1240  void
1241  resolve_indices(std::vector<types::global_dof_index> &indices) const;
1242 
1268  void
1269  condense(SparsityPattern &sparsity) const;
1270 
1274  void
1275  condense(BlockSparsityPattern &sparsity) const;
1276 
1281  void
1283 
1288  void
1290 
1298  void
1299  condense(SparseMatrix<number> &matrix) const;
1300 
1304  void
1306 
1318  template <typename VectorType>
1319  void
1320  condense(VectorType &vec) const;
1321 
1328  template <typename VectorType>
1329  void
1330  condense(const VectorType &vec_ghosted, VectorType &output) const;
1331 
1344  template <typename VectorType>
1345  void
1346  condense(SparseMatrix<number> &matrix, VectorType &vector) const;
1347 
1352  template <typename BlockVectorType>
1353  void
1354  condense(BlockSparseMatrix<number> &matrix, BlockVectorType &vector) const;
1355 
1362  template <typename VectorType>
1363  void
1364  set_zero(VectorType &vec) const;
1365 
1421  template <class InVector, class OutVector>
1422  void
1423  distribute_local_to_global(const InVector &local_vector,
1424  const std::vector<size_type> &local_dof_indices,
1425  OutVector &global_vector) const;
1426 
1474  template <typename VectorType>
1475  void
1476  distribute_local_to_global(const Vector<number> &local_vector,
1477  const std::vector<size_type> &local_dof_indices,
1478  VectorType &global_vector,
1479  const FullMatrix<number> &local_matrix) const;
1480 
1500  template <typename VectorType>
1501  void
1503  const Vector<number> &local_vector,
1504  const std::vector<size_type> &local_dof_indices_row,
1505  const std::vector<size_type> &local_dof_indices_col,
1506  VectorType &global_vector,
1507  const FullMatrix<number> &local_matrix,
1508  bool diagonal = false) const;
1509 
1513  template <typename VectorType>
1514  void
1516  const number value,
1517  VectorType &global_vector) const;
1518 
1551  template <typename ForwardIteratorVec,
1552  typename ForwardIteratorInd,
1553  typename VectorType>
1554  void
1555  distribute_local_to_global(ForwardIteratorVec local_vector_begin,
1556  ForwardIteratorVec local_vector_end,
1557  ForwardIteratorInd local_indices_begin,
1558  VectorType &global_vector) const;
1559 
1611  template <typename MatrixType>
1612  void
1613  distribute_local_to_global(const FullMatrix<number> &local_matrix,
1614  const std::vector<size_type> &local_dof_indices,
1615  MatrixType &global_matrix) const;
1616 
1644  template <typename MatrixType>
1645  void
1646  distribute_local_to_global(const FullMatrix<number> &local_matrix,
1647  const std::vector<size_type> &row_indices,
1648  const std::vector<size_type> &col_indices,
1649  MatrixType &global_matrix) const;
1650 
1667  template <typename MatrixType>
1668  void
1669  distribute_local_to_global(const FullMatrix<number> &local_matrix,
1670  const std::vector<size_type> &row_indices,
1671  const AffineConstraints &column_affine_constraints,
1672  const std::vector<size_type> &column_indices,
1673  MatrixType &global_matrix) const;
1674 
1695  template <typename MatrixType, typename VectorType>
1696  void
1697  distribute_local_to_global(const FullMatrix<number> &local_matrix,
1698  const Vector<number> &local_vector,
1699  const std::vector<size_type> &local_dof_indices,
1700  MatrixType &global_matrix,
1701  VectorType &global_vector,
1702  bool use_inhomogeneities_for_rhs = false) const;
1703 
1757  void
1759  const std::vector<size_type> &local_dof_indices,
1760  SparsityPatternBase &sparsity_pattern,
1761  const bool keep_constrained_entries = true,
1762  const Table<2, bool> &dof_mask = Table<2, bool>()) const;
1763 
1767  void
1769  const std::vector<size_type> &row_indices,
1770  const std::vector<size_type> &col_indices,
1771  SparsityPatternBase &sparsity_pattern,
1772  const bool keep_constrained_entries = true,
1773  const Table<2, bool> &dof_mask = Table<2, bool>()) const;
1774 
1779  void
1781  const std::vector<size_type> &row_indices,
1782  const AffineConstraints<number> &col_constraints,
1783  const std::vector<size_type> &col_indices,
1784  SparsityPatternBase &sparsity_pattern,
1785  const bool keep_constrained_entries = true,
1786  const Table<2, bool> &dof_mask = Table<2, bool>()) const;
1787 
1807  template <typename ForwardIteratorVec,
1808  typename ForwardIteratorInd,
1809  typename VectorType>
1810  void
1811  get_dof_values(const VectorType &global_vector,
1812  ForwardIteratorInd local_indices_begin,
1813  ForwardIteratorVec local_vector_begin,
1814  ForwardIteratorVec local_vector_end) const;
1815 
1837  template <typename VectorType>
1838  void
1839  distribute(VectorType &vec) const;
1840 
1849  {
1854  using Entries = std::vector<std::pair<size_type, number>>;
1855 
1862 
1871 
1876 
1880  ConstraintLine(const size_type &index = numbers::invalid_dof_index,
1881  const typename AffineConstraints<
1882  number>::ConstraintLine::Entries &entries = {},
1883  const number inhomogeneity = 0.0);
1884 
1888  ConstraintLine(const ConstraintLine &other) = default;
1889 
1893  ConstraintLine(ConstraintLine &&other) noexcept = default;
1894 
1898  ConstraintLine &
1899  operator=(const ConstraintLine &other) = default;
1900 
1904  ConstraintLine &
1905  operator=(ConstraintLine &&other) noexcept = default;
1906 
1911  std::size_t
1913 
1919  template <class Archive>
1920  void
1921  serialize(Archive &ar, const unsigned int)
1922  {
1923  ar &index &entries &inhomogeneity;
1924  }
1925 
1929  friend void
1931  {
1932  std::swap(l1.index, l2.index);
1933  std::swap(l1.entries, l2.entries);
1935  }
1936  };
1937 
1941  using const_iterator = typename std::vector<ConstraintLine>::const_iterator;
1942 
1946  using LineRange = boost::iterator_range<const_iterator>;
1947 
1956  LineRange
1957  get_lines() const;
1958 
1987  bool
1988  is_consistent_in_parallel(const std::vector<IndexSet> &locally_owned_dofs,
1989  const IndexSet &locally_active_dofs,
1990  const MPI_Comm mpi_communicator,
1991  const bool verbose = false) const;
1992 
1998  void
2000  const IndexSet &locally_relevant_dofs,
2001  const MPI_Comm mpi_communicator);
2002 
2010  "You are attempting an operation on an AffineConstraints object "
2011  "that requires the object to not be 'closed', i.e., for which you "
2012  "must not already have called the close() member function. But the "
2013  "object is already closed, and so the operation can not be "
2014  "performed.");
2022  "You are attempting an operation on an AffineConstraints object "
2023  "that requires the object to be 'closed', i.e., for which you "
2024  "needed to call the close() member function. But the object "
2025  "is not currently closed, and so the operation can not be "
2026  "performed.");
2033  size_type,
2034  << "The specified line " << arg1 << " does not exist.");
2041  size_type,
2042  size_type,
2043  number,
2044  number,
2045  << "The entry for the indices " << arg1 << " and " << arg2
2046  << " already exists, but the values " << arg3 << " (old) and "
2047  << arg4 << " (new) differ "
2048  << "by " << (arg4 - arg3) << '.');
2055  int,
2056  int,
2057  << "You tried to constrain DoF " << arg1 << " to DoF " << arg2
2058  << ", but that one is also constrained. This is not allowed!");
2065  size_type,
2066  << "Degree of freedom " << arg1
2067  << " is constrained from both object in a merge operation.");
2074  size_type,
2075  << "In the given argument a degree of freedom is constrained "
2076  << "to another DoF with number " << arg1
2077  << ", which however is constrained by this object. This is not"
2078  << " allowed.");
2085  size_type,
2086  << "The index set given to this constraints object indicates "
2087  << "constraints for degree of freedom " << arg1
2088  << " should not be stored by this object, but a constraint "
2089  << "is being added.");
2090 
2097  size_type,
2098  size_type,
2099  << "The index set given to this constraints object indicates "
2100  << "constraints using degree of freedom " << arg2
2101  << " should not be stored by this object, but a constraint "
2102  << "for degree of freedom " << arg1 << " uses it.");
2103 
2110  int,
2111  int,
2112  << "While distributing the constraint for DoF " << arg1
2113  << ", it turns out that one of the processors "
2114  << "who own the " << arg2 << " degrees of freedom that x_"
2115  << arg1 << " is constrained against does not know about "
2116  << "the constraint on x_" << arg1
2117  << ". Did you not initialize the AffineConstraints container "
2118  << "with the appropriate locally_relevant set so "
2119  << "that every processor who owns a DoF that constrains "
2120  << "another DoF also knows about this constraint?");
2121 
2122  template <typename>
2123  friend class AffineConstraints;
2124 
2125 private:
2137  std::vector<ConstraintLine> lines;
2138 
2171  std::vector<size_type> lines_cache;
2172 
2180 
2187 
2197 
2201  bool sorted;
2202 
2204  internal::AffineConstraints::ScratchData<number>>
2206 
2211  size_type
2212  calculate_line_index(const size_type line_n) const;
2213 
2218  template <typename MatrixType, typename VectorType>
2219  void
2221  const Vector<number> &local_vector,
2222  const std::vector<size_type> &local_dof_indices,
2223  MatrixType &global_matrix,
2224  VectorType &global_vector,
2225  const bool use_inhomogeneities_for_rhs,
2226  const std::bool_constant<false>) const;
2227 
2232  template <typename MatrixType, typename VectorType>
2233  void
2235  const Vector<number> &local_vector,
2236  const std::vector<size_type> &local_dof_indices,
2237  MatrixType &global_matrix,
2238  VectorType &global_vector,
2239  const bool use_inhomogeneities_for_rhs,
2240  const std::bool_constant<true>) const;
2241 
2249  void
2250  make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
2251  internal::AffineConstraints::GlobalRowsFromLocal<number>
2252  &global_rows) const;
2253 
2261  void
2262  make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
2263  std::vector<size_type> &active_dofs) const;
2264 
2268  template <typename MatrixScalar, typename VectorScalar>
2271  const size_type i,
2272  const internal::AffineConstraints::GlobalRowsFromLocal<number> &global_rows,
2273  const Vector<VectorScalar> &local_vector,
2274  const std::vector<size_type> &local_dof_indices,
2275  const FullMatrix<MatrixScalar> &local_matrix) const;
2276 };
2277 
2278 /* ---------------- template and inline functions ----------------- */
2279 
2280 template <typename number>
2282  : AffineConstraints<number>(IndexSet(), IndexSet())
2283 {}
2284 
2285 
2286 
2287 template <typename number>
2289  const IndexSet &locally_stored_constraints)
2290  : AffineConstraints<number>(locally_stored_constraints,
2291  locally_stored_constraints)
2292 {}
2293 
2294 
2295 
2296 template <typename number>
2299  const IndexSet &locally_stored_constraints)
2300  : lines()
2302  , local_lines(locally_stored_constraints)
2303  , sorted(false)
2304 {
2305  Assert(locally_owned_dofs.is_subset_of(locally_stored_constraints),
2306  ExcMessage("The set of locally stored constraints needs to be a "
2307  "superset of the locally owned DoFs."));
2308 
2309  // make sure the IndexSet is compressed. Otherwise this can lead to crashes
2310  // that are hard to find (only happen in release mode).
2311  // see tests/mpi/affine_constraints_crash_01
2313 }
2314 
2315 
2316 
2317 template <typename number>
2319  const AffineConstraints &affine_constraints)
2320  : Subscriptor()
2321  , lines(affine_constraints.lines)
2322  , lines_cache(affine_constraints.lines_cache)
2323  , locally_owned_dofs(affine_constraints.locally_owned_dofs)
2324  , local_lines(affine_constraints.local_lines)
2326  affine_constraints.needed_elements_for_distribute)
2327  , sorted(affine_constraints.sorted)
2328 {}
2329 
2330 
2331 
2332 template <typename number>
2333 inline void
2335 {
2336  Assert(sorted == false, ExcMatrixIsClosed());
2337 
2338  // the following can happen when we compute with distributed meshes and dof
2339  // handlers and we constrain a degree of freedom whose number we don't have
2340  // locally. if we don't abort here the program will try to allocate several
2341  // terabytes of memory to resize the various arrays below :-)
2343  const size_type line_index = calculate_line_index(line_n);
2344 
2345  // check whether line already exists; it may, in which case we can just quit
2346  if (is_constrained(line_n))
2347  return;
2348 
2349  // if necessary enlarge vector of existing entries for cache
2350  if (line_index >= lines_cache.size())
2351  lines_cache.resize(std::max(2 * static_cast<size_type>(lines_cache.size()),
2352  line_index + 1),
2354 
2355  // push a new line to the end of the list
2356  lines.emplace_back();
2357  lines.back().index = line_n;
2358  lines.back().inhomogeneity = 0.;
2359  lines_cache[line_index] = lines.size() - 1;
2360 }
2361 
2362 
2363 
2364 template <typename number>
2365 inline void
2367  const size_type column,
2368  const number weight)
2369 {
2370  Assert(sorted == false, ExcMatrixIsClosed());
2371  Assert(constrained_dof_index != column,
2372  ExcMessage("Can't constrain a degree of freedom to itself"));
2373 
2374  // Ensure that the current line is present in the cache:
2375  const size_type line_index = calculate_line_index(constrained_dof_index);
2376  Assert(line_index < lines_cache.size(),
2377  ExcMessage("The current AffineConstraints does not contain the line "
2378  "for the current entry. Call AffineConstraints::add_line "
2379  "before calling this function."));
2380 
2381  // if in debug mode, check whether an entry for this column already exists
2382  // and if it's the same as the one entered at present
2383  //
2384  // in any case: exit the function if an entry for this column already
2385  // exists, since we don't want to enter it twice
2387  ExcInternalError());
2389  ExcColumnNotStoredHere(constrained_dof_index, column));
2390  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
2391  Assert(line_ptr->index == constrained_dof_index, ExcInternalError());
2392  for (const auto &p : line_ptr->entries)
2393  if (p.first == column)
2394  {
2395  Assert(std::abs(p.second - weight) < 1.e-14,
2397  constrained_dof_index, column, p.second, weight));
2398  return;
2399  }
2400 
2401  line_ptr->entries.emplace_back(column, weight);
2402 }
2403 
2404 
2405 
2406 template <typename number>
2407 inline void
2409  const size_type constrained_dof_index,
2410  const number value)
2411 {
2412  const size_type line_index = calculate_line_index(constrained_dof_index);
2413  Assert(line_index < lines_cache.size() &&
2414  lines_cache[line_index] != numbers::invalid_size_type,
2415  ExcMessage("call add_line() before calling set_inhomogeneity()"));
2416  Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
2417  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
2418  line_ptr->inhomogeneity = value;
2419 }
2420 
2421 
2422 
2423 template <typename number>
2424 template <typename VectorType>
2425 inline void
2427 {
2428  // since lines is a private member, we cannot pass it to the functions
2429  // above. therefore, copy the content which is cheap
2430  std::vector<size_type> constrained_lines(lines.size());
2431  for (unsigned int i = 0; i < lines.size(); ++i)
2432  constrained_lines[i] = lines[i].index;
2433  internal::AffineConstraintsImplementation::set_zero_all(constrained_lines,
2434  vec);
2435 }
2436 
2437 template <typename number>
2440 {
2441  return lines.size();
2442 }
2443 
2444 template <typename number>
2447 {
2448  return std::count_if(lines.begin(),
2449  lines.end(),
2450  [](const ConstraintLine &line) {
2451  return (line.entries.size() == 1) &&
2452  (line.entries[0].second == number(1.));
2453  });
2454 }
2455 
2456 template <typename number>
2459 {
2460  return std::count_if(lines.begin(),
2461  lines.end(),
2462  [](const ConstraintLine &line) {
2463  return (line.inhomogeneity != number(0.));
2464  });
2465 }
2466 
2467 template <typename number>
2468 inline bool
2470 {
2471  if (lines.empty())
2472  return false;
2473 
2474  const size_type line_index = calculate_line_index(index);
2475  return ((line_index < lines_cache.size()) &&
2476  (lines_cache[line_index] != numbers::invalid_size_type));
2477 }
2478 
2479 template <typename number>
2480 inline bool
2482  const size_type line_n) const
2483 {
2484  // check whether the entry is constrained. could use is_constrained, but
2485  // that means computing the line index twice
2486  const size_type line_index = calculate_line_index(line_n);
2487  if (line_index >= lines_cache.size() ||
2488  lines_cache[line_index] == numbers::invalid_size_type)
2489  return false;
2490  else
2491  {
2492  Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
2493  return (lines[lines_cache[line_index]].inhomogeneity != number(0.));
2494  }
2495 }
2496 
2497 template <typename number>
2498 inline const std::vector<std::pair<types::global_dof_index, number>> *
2500 {
2501  if (lines.empty())
2502  return nullptr;
2503 
2504  // check whether the entry is constrained. could use is_constrained, but
2505  // that means computing the line index twice
2506  const size_type line_index = calculate_line_index(line_n);
2507  if (line_index >= lines_cache.size() ||
2508  lines_cache[line_index] == numbers::invalid_size_type)
2509  return nullptr;
2510  else
2511  return &lines[lines_cache[line_index]].entries;
2512 }
2513 
2514 template <typename number>
2515 inline number
2517 {
2518  // check whether the entry is constrained. could use is_constrained, but
2519  // that means computing the line index twice
2520  const size_type line_index = calculate_line_index(line_n);
2521  if (line_index >= lines_cache.size() ||
2522  lines_cache[line_index] == numbers::invalid_size_type)
2523  return 0;
2524  else
2525  return lines[lines_cache[line_index]].inhomogeneity;
2526 }
2527 
2528 template <typename number>
2531 {
2532  // IndexSet is unused (serial case)
2533  if (local_lines.size() == 0)
2534  return line_n;
2535 
2537 
2538  return local_lines.index_within_set(line_n);
2539 }
2540 
2541 
2542 
2543 template <typename number>
2544 inline bool
2546 {
2547  return local_lines.size() == 0 || local_lines.is_element(line_n);
2548 }
2549 
2550 
2551 
2552 template <typename number>
2553 inline const IndexSet &
2555 {
2556  return locally_owned_dofs;
2557 }
2558 
2559 
2560 
2561 template <typename number>
2562 inline const IndexSet &
2564 {
2565  return local_lines;
2566 }
2567 
2568 
2569 
2570 template <typename number>
2571 template <typename VectorType>
2572 inline void
2574  const size_type index,
2575  const number value,
2576  VectorType &global_vector) const
2577 {
2578  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2579 
2580  if (is_constrained(index) == false)
2581  global_vector(index) += value;
2582  else
2583  {
2584  const ConstraintLine &position =
2586  for (size_type j = 0; j < position.entries.size(); ++j)
2587  global_vector(position.entries[j].first) +=
2588  value * position.entries[j].second;
2589  }
2590 }
2591 
2592 template <typename number>
2593 template <typename ForwardIteratorVec,
2594  typename ForwardIteratorInd,
2595  typename VectorType>
2596 inline void
2598  ForwardIteratorVec local_vector_begin,
2599  ForwardIteratorVec local_vector_end,
2600  ForwardIteratorInd local_indices_begin,
2601  VectorType &global_vector) const
2602 {
2603  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2604  for (; local_vector_begin != local_vector_end;
2605  ++local_vector_begin, ++local_indices_begin)
2606  {
2607  if (is_constrained(*local_indices_begin) == false)
2608  internal::ElementAccess<VectorType>::add(*local_vector_begin,
2609  *local_indices_begin,
2610  global_vector);
2611  else
2612  {
2613  const ConstraintLine &position =
2614  lines[lines_cache[calculate_line_index(*local_indices_begin)]];
2615  for (size_type j = 0; j < position.entries.size(); ++j)
2617  (*local_vector_begin) * position.entries[j].second,
2618  position.entries[j].first,
2619  global_vector);
2620  }
2621  }
2622 }
2623 
2624 template <typename number>
2625 template <class InVector, class OutVector>
2626 inline void
2628  const InVector &local_vector,
2629  const std::vector<size_type> &local_dof_indices,
2630  OutVector &global_vector) const
2631 {
2632  Assert(local_vector.size() == local_dof_indices.size(),
2633  ExcDimensionMismatch(local_vector.size(), local_dof_indices.size()));
2634  distribute_local_to_global(local_vector.begin(),
2635  local_vector.end(),
2636  local_dof_indices.begin(),
2637  global_vector);
2638 }
2639 
2640 template <typename number>
2641 template <typename ForwardIteratorVec,
2642  typename ForwardIteratorInd,
2643  typename VectorType>
2644 inline void
2646  const VectorType &global_vector,
2647  ForwardIteratorInd local_indices_begin,
2648  ForwardIteratorVec local_vector_begin,
2649  ForwardIteratorVec local_vector_end) const
2650 {
2651  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2652  for (; local_vector_begin != local_vector_end;
2653  ++local_vector_begin, ++local_indices_begin)
2654  {
2655  if (is_constrained(*local_indices_begin) == false)
2656  *local_vector_begin = global_vector(*local_indices_begin);
2657  else
2658  {
2659  const ConstraintLine &position =
2660  lines[lines_cache[calculate_line_index(*local_indices_begin)]];
2661  typename VectorType::value_type value = position.inhomogeneity;
2662  for (size_type j = 0; j < position.entries.size(); ++j)
2663  value += (global_vector(position.entries[j].first) *
2664  position.entries[j].second);
2665  *local_vector_begin = value;
2666  }
2667  }
2668 }
2669 
2670 // Forward declarations
2671 #ifndef DOXYGEN
2672 template <typename MatrixType>
2673 class BlockMatrixBase;
2674 template <typename SparsityPatternType>
2676 template <typename number>
2677 class BlockSparseMatrixEZ;
2678 
2679 namespace internal
2680 {
2681  namespace AffineConstraints
2682  {
2700  template <typename MatrixType>
2701  struct IsBlockMatrix
2702  {
2703  private:
2709  template <typename T>
2710  static std::true_type
2711  check(const BlockMatrixBase<T> *);
2712 
2717  template <typename T>
2718  static std::true_type
2719  check(const BlockSparseMatrixEZ<T> *);
2720 
2725  static std::false_type
2726  check(...);
2727 
2728  public:
2735  static const bool value =
2736  std::is_same_v<decltype(check(std::declval<MatrixType *>())),
2737  std::true_type>;
2738  };
2739 
2740  // instantiation of the static member
2741  template <typename MatrixType>
2742  const bool IsBlockMatrix<MatrixType>::value;
2743 
2744  } // namespace AffineConstraints
2745 } // namespace internal
2746 #endif
2747 
2748 
2749 
2750 template <typename number>
2751 template <typename other_number>
2752 inline void
2754  const AffineConstraints<other_number> &other)
2755 {
2756  lines.clear();
2757  lines.reserve(other.lines.size());
2758 
2759  for (const auto l : other.lines)
2760  lines.emplace_back(l.index,
2761  typename ConstraintLine::Entries(l.entries.begin(),
2762  l.entries.end()),
2763  l.inhomogeneity);
2764 
2765  lines_cache = other.lines_cache;
2766  local_lines = other.local_lines;
2767  sorted = other.sorted;
2768 
2769  locally_owned_dofs = other.locally_owned_dofs;
2770  needed_elements_for_distribute = other.needed_elements_for_distribute;
2771 }
2772 
2773 
2774 template <typename number>
2775 template <typename other_number>
2776 void
2778  const AffineConstraints<other_number> &other_constraints,
2779  const MergeConflictBehavior merge_conflict_behavior,
2780  const bool allow_different_local_lines)
2781 {
2782  (void)allow_different_local_lines;
2783  Assert(allow_different_local_lines ||
2784  local_lines == other_constraints.local_lines,
2785  ExcMessage(
2786  "local_lines for this and the other objects are not the same "
2787  "although allow_different_local_lines is false."));
2788 
2789  // store the previous state with respect to sorting
2790  const bool object_was_sorted = sorted;
2791  sorted = false;
2792 
2793  // first action is to fold into the present object possible constraints
2794  // in the second object. we don't strictly need to do this any more since
2795  // the AffineConstraints container has learned to deal with chains of
2796  // constraints in the close() function, but we have traditionally done
2797  // this and it's not overly hard to do.
2798  //
2799  // for this, loop over all constraints and replace the constraint lines
2800  // with a new one where constraints are replaced if necessary.
2801  typename ConstraintLine::Entries tmp;
2802  for (ConstraintLine &line : lines)
2803  {
2804  tmp.clear();
2805  for (const std::pair<size_type, number> &entry : line.entries)
2806  {
2807  // if the present dof is not stored, or not constrained, or if we
2808  // won't take the constraint from the other object, then simply copy
2809  // it over
2810  if ((other_constraints.local_lines.size() != 0. &&
2811  other_constraints.local_lines.is_element(entry.first) ==
2812  false) ||
2813  other_constraints.is_constrained(entry.first) == false ||
2814  ((merge_conflict_behavior != right_object_wins) &&
2815  other_constraints.is_constrained(entry.first) &&
2816  this->is_constrained(entry.first)))
2817  tmp.push_back(entry);
2818  else
2819  // otherwise resolve further constraints by replacing the old
2820  // entry by a sequence of new entries taken from the other
2821  // object, but with multiplied weights
2822  {
2823  const auto *other_entries =
2824  other_constraints.get_constraint_entries(entry.first);
2825  Assert(other_entries != nullptr, ExcInternalError());
2826 
2827  const number weight = entry.second;
2828 
2829  for (const auto &other_entry : *other_entries)
2830  tmp.emplace_back(other_entry.first,
2831  other_entry.second * weight);
2832 
2833  line.inhomogeneity +=
2834  other_constraints.get_inhomogeneity(entry.first) * weight;
2835  }
2836  }
2837  // finally exchange old and newly resolved line
2838  line.entries.swap(tmp);
2839  }
2840 
2841  if (local_lines.size() != 0)
2842  local_lines.add_indices(other_constraints.local_lines);
2843 
2844  {
2845  // do not bother to resize the lines cache exactly since it is pretty
2846  // cheap to adjust it along the way.
2847  std::fill(lines_cache.begin(),
2848  lines_cache.end(),
2850 
2851  // reset lines_cache for our own constraints
2852  size_type index = 0;
2853  for (const ConstraintLine &line : lines)
2854  {
2855  const size_type local_line_no = calculate_line_index(line.index);
2856  if (local_line_no >= lines_cache.size())
2857  lines_cache.resize(local_line_no + 1, numbers::invalid_size_type);
2858  lines_cache[local_line_no] = index++;
2859  }
2860 
2861  // Add other_constraints to lines cache and our list of constraints
2862  for (const auto &line : other_constraints.lines)
2863  {
2864  const size_type local_line_no = calculate_line_index(line.index);
2865  if (local_line_no >= lines_cache.size())
2866  {
2867  lines_cache.resize(local_line_no + 1, numbers::invalid_size_type);
2868  lines.emplace_back(line.index,
2869  typename ConstraintLine::Entries(
2870  line.entries.begin(), line.entries.end()),
2871  line.inhomogeneity);
2872  lines_cache[local_line_no] = index++;
2873  }
2874  else if (lines_cache[local_line_no] == numbers::invalid_size_type)
2875  {
2876  // there are no constraints for that line yet
2877  lines.emplace_back(line.index,
2878  typename ConstraintLine::Entries(
2879  line.entries.begin(), line.entries.end()),
2880  line.inhomogeneity);
2881  AssertIndexRange(local_line_no, lines_cache.size());
2882  lines_cache[local_line_no] = index++;
2883  }
2884  else
2885  {
2886  // we already store that line
2887  switch (merge_conflict_behavior)
2888  {
2889  case no_conflicts_allowed:
2890  AssertThrow(false,
2891  ExcDoFIsConstrainedFromBothObjects(line.index));
2892  break;
2893 
2894  case left_object_wins:
2895  // ignore this constraint
2896  break;
2897 
2898  case right_object_wins:
2899  AssertIndexRange(local_line_no, lines_cache.size());
2900  lines[lines_cache[local_line_no]] = {
2901  line.index,
2902  typename ConstraintLine::Entries(line.entries.begin(),
2903  line.entries.end()),
2904  static_cast<number>(line.inhomogeneity)};
2905  break;
2906 
2907  default:
2908  Assert(false, ExcNotImplemented());
2909  }
2910  }
2911  }
2912 
2913  // check that we set the pointers correctly
2914  for (size_type i = 0; i < lines_cache.size(); ++i)
2915  if (lines_cache[i] != numbers::invalid_size_type)
2916  Assert(i == calculate_line_index(lines[lines_cache[i]].index),
2917  ExcInternalError());
2918  }
2919 
2920  // if the object was sorted before, then make sure it is so afterward as
2921  // well. otherwise leave everything in the unsorted state
2922  if (object_was_sorted == true)
2923  close();
2924 }
2925 
2926 
2927 
2928 template <typename number>
2929 template <typename MatrixType>
2930 inline void
2932  const FullMatrix<number> &local_matrix,
2933  const std::vector<size_type> &local_dof_indices,
2934  MatrixType &global_matrix) const
2935 {
2936  // create a dummy and hand on to the function actually implementing this
2937  // feature in the cm.templates.h file.
2939  distribute_local_to_global(
2940  local_matrix,
2941  dummy,
2942  local_dof_indices,
2943  global_matrix,
2944  dummy,
2945  false,
2946  std::integral_constant<
2947  bool,
2948  internal::AffineConstraints::IsBlockMatrix<MatrixType>::value>());
2949 }
2950 
2951 
2952 
2953 template <typename number>
2954 template <typename MatrixType, typename VectorType>
2955 inline void
2957  const FullMatrix<number> &local_matrix,
2958  const Vector<number> &local_vector,
2959  const std::vector<size_type> &local_dof_indices,
2960  MatrixType &global_matrix,
2961  VectorType &global_vector,
2962  bool use_inhomogeneities_for_rhs) const
2963 {
2964  // enter the internal function with the respective block information set,
2965  // the actual implementation follows in the cm.templates.h file.
2966  distribute_local_to_global(
2967  local_matrix,
2968  local_vector,
2969  local_dof_indices,
2970  global_matrix,
2971  global_vector,
2972  use_inhomogeneities_for_rhs,
2973  std::integral_constant<
2974  bool,
2975  internal::AffineConstraints::IsBlockMatrix<MatrixType>::value>());
2976 }
2977 
2978 
2979 
2980 template <typename number>
2982  const size_type &index,
2984  const number inhomogeneity)
2985  : index(index)
2986  , entries(entries)
2987  , inhomogeneity(inhomogeneity)
2988 {}
2989 
2990 
2991 
2993 
2994 #endif
bool is_closed() const
size_type n_inhomogeneities() const
types::global_dof_index size_type
bool has_inhomogeneities() const
LineRange get_lines() const
void add_line(const size_type line_n)
const IndexSet & get_locally_owned_indices() const
void add_selected_constraints(const AffineConstraints &constraints_in, const IndexSet &filter)
bool can_store_line(const size_type line_n) const
Threads::ThreadLocalStorage< internal::AffineConstraints::ScratchData< number > > scratch_data
friend class AffineConstraints
void merge(const AffineConstraints< other_number > &other_constraints, const MergeConflictBehavior merge_conflict_behavior=no_conflicts_allowed, const bool allow_different_local_lines=false)
void add_entry(const size_type constrained_dof_index, const size_type column, const number weight)
void make_consistent_in_parallel(const IndexSet &locally_owned_dofs, const IndexSet &locally_relevant_dofs, const MPI_Comm mpi_communicator)
number get_inhomogeneity(const size_type line_n) const
void make_sorted_row_list(const std::vector< size_type > &local_dof_indices, std::vector< size_type > &active_dofs) const
void get_dof_values(const VectorType &global_vector, ForwardIteratorInd local_indices_begin, ForwardIteratorVec local_vector_begin, ForwardIteratorVec local_vector_end) const
void distribute_local_to_global(const InVector &local_vector, const std::vector< size_type > &local_dof_indices, OutVector &global_vector) const
void make_sorted_row_list(const std::vector< size_type > &local_dof_indices, internal::AffineConstraints::GlobalRowsFromLocal< number > &global_rows) const
const IndexSet & get_local_lines() const
void add_entries_local_to_global(const std::vector< size_type > &local_dof_indices, SparsityPatternBase &sparsity_pattern, const bool keep_constrained_entries=true, const Table< 2, bool > &dof_mask=Table< 2, bool >()) const
void shift(const size_type offset)
size_type n_identities() const
std::size_t memory_consumption() const
typename std::vector< ConstraintLine >::const_iterator const_iterator
void set_inhomogeneity(const size_type constrained_dof_index, const number value)
void condense(SparsityPattern &sparsity) const
size_type max_constraint_indirections() const
void distribute(VectorType &vec) const
void resolve_indices(std::vector< types::global_dof_index > &indices) const
std::vector< ConstraintLine > lines
bool is_constrained(const size_type line_n) const
AffineConstraints get_view(const IndexSet &mask) const
IndexSet needed_elements_for_distribute
const std::vector< std::pair< size_type, number > > * get_constraint_entries(const size_type line_n) const
size_type calculate_line_index(const size_type line_n) const
void distribute_local_to_global(const FullMatrix< number > &local_matrix, const Vector< number > &local_vector, const std::vector< size_type > &local_dof_indices, MatrixType &global_matrix, VectorType &global_vector, const bool use_inhomogeneities_for_rhs, const std::bool_constant< true >) const
bool is_inhomogeneously_constrained(const size_type index) const
void copy_from(const AffineConstraints< other_number > &other)
std::vector< size_type > lines_cache
void constrain_dof_to_zero(const size_type constrained_dof)
void add_constraint(const size_type constrained_dof, const ArrayView< const std::pair< size_type, number >> &dependencies, const number inhomogeneity=0)
void distribute_local_to_global(const FullMatrix< number > &local_matrix, const Vector< number > &local_vector, const std::vector< size_type > &local_dof_indices, MatrixType &global_matrix, VectorType &global_vector, const bool use_inhomogeneities_for_rhs, const std::bool_constant< false >) const
bool is_consistent_in_parallel(const std::vector< IndexSet > &locally_owned_dofs, const IndexSet &locally_active_dofs, const MPI_Comm mpi_communicator, const bool verbose=false) const
bool is_identity_constrained(const size_type line_n) const
size_type n_constraints() const
void print(std::ostream &out) const
void write_dot(std::ostream &) const
void set_zero(VectorType &vec) const
void add_lines(const std::vector< bool > &lines)
boost::iterator_range< const_iterator > LineRange
bool are_identity_constrained(const size_type line_n_1, const size_type line_n_2) const
void add_entries(const size_type constrained_dof_index, const std::vector< std::pair< size_type, number >> &col_weight_pairs)
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
bool is_subset_of(const IndexSet &other) const
Definition: index_set.cc:704
size_type size() const
Definition: index_set.h:1761
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1986
bool is_element(const size_type index) const
Definition: index_set.h:1879
void compress() const
Definition: index_set.h:1769
A class that provides a separate storage location on each thread that accesses the object.
Definition: vector.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
__global__ void set(Number *val, const Number s, const size_type N)
static ::ExceptionBase & ExcLineInexistant(size_type arg1)
static ::ExceptionBase & ExcIncorrectConstraint(int arg1, int arg2)
static ::ExceptionBase & ExcInternalError()
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:586
static ::ExceptionBase & ExcColumnNotStoredHere(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcDoFConstrainedToConstrainedDoF(int arg1, int arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcDoFIsConstrainedToConstrainedDoF(size_type arg1)
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcDoFIsConstrainedFromBothObjects(size_type arg1)
static ::ExceptionBase & ExcRowNotStoredHere(size_type arg1)
static ::ExceptionBase & ExcMatrixNotClosed()
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:540
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:495
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:517
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcMatrixIsClosed()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcEntryAlreadyExists(size_type arg1, size_type arg2, number arg3, number arg4)
@ matrix
Contents is actually a matrix.
@ diagonal
Matrix is diagonal.
types::global_dof_index size_type
Definition: cuda_kernels.h:45
void swap(MemorySpaceData< T, MemorySpace > &u, MemorySpaceData< T, MemorySpace > &v)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
const types::global_dof_index invalid_dof_index
Definition: types.h:253
const types::global_dof_index invalid_size_type
Definition: types.h:234
Definition: types.h:33
unsigned int global_dof_index
Definition: types.h:82
void serialize(Archive &ar, const unsigned int)
friend void swap(ConstraintLine &l1, ConstraintLine &l2)
ConstraintLine(const ConstraintLine &other)=default
ConstraintLine(const size_type &index=numbers::invalid_dof_index, const typename AffineConstraints< number >::ConstraintLine::Entries &entries={}, const number inhomogeneity=0.0)
std::vector< std::pair< size_type, number > > Entries
ConstraintLine & operator=(ConstraintLine &&other) noexcept=default
ConstraintLine(ConstraintLine &&other) noexcept=default
std::size_t memory_consumption() const
ConstraintLine & operator=(const ConstraintLine &other)=default
typename internal::ProductTypeImpl< std::decay_t< T >, std::decay_t< U > >::type type
static void add(const typename VectorType::value_type value, const types::global_dof_index i, VectorType &V)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
const MPI_Comm comm