Reference documentation for deal.II version Git 064a48c9cd 2021-01-23 10:50:23 -0500
\(\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 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_affine_constraints_h
17 #define dealii_affine_constraints_h
18 
19 #include <deal.II/base/config.h>
20 
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  void
598  copy_from(const AffineConstraints &other);
599 
606  void
607  reinit(const IndexSet &local_constraints = IndexSet());
608 
615  bool
616  can_store_line(const size_type line_n) const;
617 
624  const IndexSet &
625  get_local_lines() const;
626 
646  void
647  add_selected_constraints(const AffineConstraints &constraints_in,
648  const IndexSet & filter);
649 
659  void
660  add_line(const size_type line_n);
661 
674  void
675  add_lines(const std::vector<bool> &lines);
676 
689  void
690  add_lines(const std::set<size_type> &lines);
691 
704  void
705  add_lines(const IndexSet &lines);
706 
722  void
723  add_entry(const size_type constrained_dof_index,
724  const size_type column,
725  const number weight);
726 
732  void
733  add_entries(
734  const size_type constrained_dof_index,
735  const std::vector<std::pair<size_type, number>> &col_weight_pairs);
736 
748  void
749  set_inhomogeneity(const size_type constrained_dof_index, const number value);
750 
772  void
773  close();
774 
799  void
800  merge(
801  const AffineConstraints & other_constraints,
802  const MergeConflictBehavior merge_conflict_behavior = no_conflicts_allowed,
803  const bool allow_different_local_lines = false);
804 
817  void
818  shift(const size_type offset);
819 
827  void
828  clear();
829 
842  size_type
843  n_constraints() const;
844 
854  bool
855  is_constrained(const size_type line_n) const;
856 
868  bool
869  is_identity_constrained(const size_type line_n) const;
870 
877  bool
878  are_identity_constrained(const size_type line_n_1,
879  const size_type line_n_2) const;
880 
891  size_type
892  max_constraint_indirections() const;
893 
898  bool
899  is_inhomogeneously_constrained(const size_type index) const;
900 
906  bool
907  has_inhomogeneities() const;
908 
913  const std::vector<std::pair<size_type, number>> *
914  get_constraint_entries(const size_type line_n) const;
915 
920  number
921  get_inhomogeneity(const size_type line_n) const;
922 
943  void
944  print(std::ostream &out) const;
945 
958  void
959  write_dot(std::ostream &) const;
960 
965  std::size_t
966  memory_consumption() const;
967 
974  void
975  resolve_indices(std::vector<types::global_dof_index> &indices) const;
976 
1002  void
1003  condense(SparsityPattern &sparsity) const;
1004 
1008  void
1009  condense(BlockSparsityPattern &sparsity) const;
1010 
1015  void
1016  condense(DynamicSparsityPattern &sparsity) const;
1017 
1022  void
1023  condense(BlockDynamicSparsityPattern &sparsity) const;
1024 
1032  void
1033  condense(SparseMatrix<number> &matrix) const;
1034 
1038  void
1039  condense(BlockSparseMatrix<number> &matrix) const;
1040 
1052  template <class VectorType>
1053  void
1054  condense(VectorType &vec) const;
1055 
1062  template <class VectorType>
1063  void
1064  condense(const VectorType &vec_ghosted, VectorType &output) const;
1065 
1078  template <class VectorType>
1079  void
1080  condense(SparseMatrix<number> &matrix, VectorType &vector) const;
1081 
1086  template <class BlockVectorType>
1087  void
1088  condense(BlockSparseMatrix<number> &matrix, BlockVectorType &vector) const;
1089 
1096  template <class VectorType>
1097  void
1098  set_zero(VectorType &vec) const;
1099 
1155  template <class InVector, class OutVector>
1156  void
1157  distribute_local_to_global(const InVector & local_vector,
1158  const std::vector<size_type> &local_dof_indices,
1159  OutVector & global_vector) const;
1160 
1208  template <typename VectorType>
1209  void
1210  distribute_local_to_global(const Vector<number> & local_vector,
1211  const std::vector<size_type> &local_dof_indices,
1212  VectorType & global_vector,
1213  const FullMatrix<number> & local_matrix) const;
1214 
1234  template <typename VectorType>
1235  void
1236  distribute_local_to_global(
1237  const Vector<number> & local_vector,
1238  const std::vector<size_type> &local_dof_indices_row,
1239  const std::vector<size_type> &local_dof_indices_col,
1240  VectorType & global_vector,
1241  const FullMatrix<number> & local_matrix,
1242  bool diagonal = false) const;
1243 
1247  template <class VectorType>
1248  void
1249  distribute_local_to_global(const size_type index,
1250  const number value,
1251  VectorType & global_vector) const;
1252 
1285  template <typename ForwardIteratorVec,
1286  typename ForwardIteratorInd,
1287  class VectorType>
1288  void
1289  distribute_local_to_global(ForwardIteratorVec local_vector_begin,
1290  ForwardIteratorVec local_vector_end,
1291  ForwardIteratorInd local_indices_begin,
1292  VectorType & global_vector) const;
1293 
1345  template <typename MatrixType>
1346  void
1347  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1348  const std::vector<size_type> &local_dof_indices,
1349  MatrixType & global_matrix) const;
1350 
1378  template <typename MatrixType>
1379  void
1380  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1381  const std::vector<size_type> &row_indices,
1382  const std::vector<size_type> &col_indices,
1383  MatrixType & global_matrix) const;
1384 
1401  template <typename MatrixType>
1402  void
1403  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1404  const std::vector<size_type> &row_indices,
1405  const AffineConstraints &column_affine_constraints,
1406  const std::vector<size_type> &column_indices,
1407  MatrixType & global_matrix) const;
1408 
1429  template <typename MatrixType, typename VectorType>
1430  void
1431  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1432  const Vector<number> & local_vector,
1433  const std::vector<size_type> &local_dof_indices,
1434  MatrixType & global_matrix,
1435  VectorType & global_vector,
1436  bool use_inhomogeneities_for_rhs = false) const;
1437 
1491  template <typename SparsityPatternType>
1492  void
1493  add_entries_local_to_global(
1494  const std::vector<size_type> &local_dof_indices,
1495  SparsityPatternType & sparsity_pattern,
1496  const bool keep_constrained_entries = true,
1497  const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1498 
1502  template <typename SparsityPatternType>
1503  void
1504  add_entries_local_to_global(
1505  const std::vector<size_type> &row_indices,
1506  const std::vector<size_type> &col_indices,
1507  SparsityPatternType & sparsity_pattern,
1508  const bool keep_constrained_entries = true,
1509  const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1510 
1530  template <typename ForwardIteratorVec,
1531  typename ForwardIteratorInd,
1532  class VectorType>
1533  void
1534  get_dof_values(const VectorType & global_vector,
1535  ForwardIteratorInd local_indices_begin,
1536  ForwardIteratorVec local_vector_begin,
1537  ForwardIteratorVec local_vector_end) const;
1538 
1560  template <class VectorType>
1561  void
1562  distribute(VectorType &vec) const;
1563 
1572  {
1577  using Entries = std::vector<std::pair<size_type, number>>;
1578 
1585 
1594 
1599 
1607  bool
1608  operator<(const ConstraintLine &) const;
1609 
1615  bool
1616  operator==(const ConstraintLine &) const;
1617 
1622  std::size_t
1623  memory_consumption() const;
1624 
1630  template <class Archive>
1631  void
1632  serialize(Archive &ar, const unsigned int)
1633  {
1634  ar &index &entries &inhomogeneity;
1635  }
1636  };
1637 
1641  using const_iterator = typename std::vector<ConstraintLine>::const_iterator;
1642 
1646  using LineRange = boost::iterator_range<const_iterator>;
1647 
1656  const LineRange
1657  get_lines() const;
1658 
1687  bool
1688  is_consistent_in_parallel(const std::vector<IndexSet> &locally_owned_dofs,
1689  const IndexSet & locally_active_dofs,
1690  const MPI_Comm & mpi_communicator,
1691  const bool verbose = false) const;
1692 
1698  DeclException0(ExcMatrixIsClosed);
1704  DeclException0(ExcMatrixNotClosed);
1711  size_type,
1712  << "The specified line " << arg1 << " does not exist.");
1718  DeclException4(ExcEntryAlreadyExists,
1719  size_type,
1720  size_type,
1721  number,
1722  number,
1723  << "The entry for the indices " << arg1 << " and " << arg2
1724  << " already exists, but the values " << arg3 << " (old) and "
1725  << arg4 << " (new) differ "
1726  << "by " << (arg4 - arg3) << ".");
1732  DeclException2(ExcDoFConstrainedToConstrainedDoF,
1733  int,
1734  int,
1735  << "You tried to constrain DoF " << arg1 << " to DoF " << arg2
1736  << ", but that one is also constrained. This is not allowed!");
1742  DeclException1(ExcDoFIsConstrainedFromBothObjects,
1743  size_type,
1744  << "Degree of freedom " << arg1
1745  << " is constrained from both object in a merge operation.");
1751  DeclException1(ExcDoFIsConstrainedToConstrainedDoF,
1752  size_type,
1753  << "In the given argument a degree of freedom is constrained "
1754  << "to another DoF with number " << arg1
1755  << ", which however is constrained by this object. This is not"
1756  << " allowed.");
1762  DeclException1(ExcRowNotStoredHere,
1763  size_type,
1764  << "The index set given to this constraints object indicates "
1765  << "constraints for degree of freedom " << arg1
1766  << " should not be stored by this object, but a constraint "
1767  << "is being added.");
1768 
1774  DeclException2(ExcColumnNotStoredHere,
1775  size_type,
1776  size_type,
1777  << "The index set given to this constraints object indicates "
1778  << "constraints using degree of freedom " << arg2
1779  << " should not be stored by this object, but a constraint "
1780  << "for degree of freedom " << arg1 << " uses it.");
1781 
1787  DeclException2(ExcIncorrectConstraint,
1788  int,
1789  int,
1790  << "While distributing the constraint for DoF " << arg1
1791  << ", it turns out that one of the processors "
1792  << "who own the " << arg2 << " degrees of freedom that x_"
1793  << arg1 << " is constrained against does not know about "
1794  << "the constraint on x_" << arg1
1795  << ". Did you not initialize the AffineConstraints container "
1796  << "with the appropriate locally_relevant set so "
1797  << "that every processor who owns a DoF that constrains "
1798  << "another DoF also knows about this constraint?");
1799 
1800 private:
1812  std::vector<ConstraintLine> lines;
1813 
1846  std::vector<size_type> lines_cache;
1847 
1854 
1858  bool sorted;
1859 
1861  internal::AffineConstraints::ScratchData<number>>
1863 
1868  size_type
1869  calculate_line_index(const size_type line_n) const;
1870 
1875  template <typename MatrixType, typename VectorType>
1876  void
1877  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1878  const Vector<number> & local_vector,
1879  const std::vector<size_type> &local_dof_indices,
1880  MatrixType & global_matrix,
1881  VectorType & global_vector,
1882  const bool use_inhomogeneities_for_rhs,
1883  const std::integral_constant<bool, false>) const;
1884 
1889  template <typename MatrixType, typename VectorType>
1890  void
1891  distribute_local_to_global(const FullMatrix<number> & local_matrix,
1892  const Vector<number> & local_vector,
1893  const std::vector<size_type> &local_dof_indices,
1894  MatrixType & global_matrix,
1895  VectorType & global_vector,
1896  const bool use_inhomogeneities_for_rhs,
1897  const std::integral_constant<bool, true>) const;
1898 
1903  template <typename SparsityPatternType>
1904  void
1905  add_entries_local_to_global(const std::vector<size_type> &local_dof_indices,
1906  SparsityPatternType & sparsity_pattern,
1907  const bool keep_constrained_entries,
1908  const Table<2, bool> &dof_mask,
1909  const std::integral_constant<bool, false>) const;
1910 
1915  template <typename SparsityPatternType>
1916  void
1917  add_entries_local_to_global(const std::vector<size_type> &local_dof_indices,
1918  SparsityPatternType & sparsity_pattern,
1919  const bool keep_constrained_entries,
1920  const Table<2, bool> &dof_mask,
1921  const std::integral_constant<bool, true>) const;
1922 
1930  void
1931  make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
1932  internal::AffineConstraints::GlobalRowsFromLocal<number>
1933  &global_rows) const;
1934 
1942  void
1943  make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
1944  std::vector<size_type> & active_dofs) const;
1945 
1949  template <typename MatrixScalar, typename VectorScalar>
1951  resolve_vector_entry(
1952  const size_type i,
1953  const internal::AffineConstraints::GlobalRowsFromLocal<number> &global_rows,
1954  const Vector<VectorScalar> & local_vector,
1955  const std::vector<size_type> & local_dof_indices,
1956  const FullMatrix<MatrixScalar> &local_matrix) const;
1957 };
1958 
1959 /* ---------------- template and inline functions ----------------- */
1960 
1961 template <typename number>
1963  const IndexSet &local_constraints)
1964  : lines()
1965  , local_lines(local_constraints)
1966  , sorted(false)
1967 {
1968  // make sure the IndexSet is compressed. Otherwise this can lead to crashes
1969  // that are hard to find (only happen in release mode).
1970  // see tests/mpi/affine_constraints_crash_01
1972 }
1973 
1974 template <typename number>
1976  const AffineConstraints &affine_constraints)
1977  : Subscriptor()
1978  , lines(affine_constraints.lines)
1979  , lines_cache(affine_constraints.lines_cache)
1980  , local_lines(affine_constraints.local_lines)
1981  , sorted(affine_constraints.sorted)
1982 {}
1983 
1984 template <typename number>
1985 inline void
1987 {
1988  Assert(sorted == false, ExcMatrixIsClosed());
1989 
1990  // the following can happen when we compute with distributed meshes and dof
1991  // handlers and we constrain a degree of freedom whose number we don't have
1992  // locally. if we don't abort here the program will try to allocate several
1993  // terabytes of memory to resize the various arrays below :-)
1995  const size_type line_index = calculate_line_index(line_n);
1996 
1997  // check whether line already exists; it may, in which case we can just quit
1998  if (is_constrained(line_n))
1999  return;
2000 
2001  // if necessary enlarge vector of existing entries for cache
2002  if (line_index >= lines_cache.size())
2003  lines_cache.resize(std::max(2 * static_cast<size_type>(lines_cache.size()),
2004  line_index + 1),
2006 
2007  // push a new line to the end of the list
2008  lines.emplace_back();
2009  lines.back().index = line_n;
2010  lines.back().inhomogeneity = 0.;
2011  lines_cache[line_index] = lines.size() - 1;
2012 }
2013 
2014 
2015 
2016 template <typename number>
2017 inline void
2019  const size_type column,
2020  const number weight)
2021 {
2022  Assert(sorted == false, ExcMatrixIsClosed());
2023  Assert(constrained_dof_index != column,
2024  ExcMessage("Can't constrain a degree of freedom to itself"));
2025 
2026  // Ensure that the current line is present in the cache:
2027  const size_type line_index = calculate_line_index(constrained_dof_index);
2028  Assert(line_index < lines_cache.size(),
2029  ExcMessage("The current AffineConstraints does not contain the line "
2030  "for the current entry. Call AffineConstraints::add_line "
2031  "before calling this function."));
2032 
2033  // if in debug mode, check whether an entry for this column already exists
2034  // and if it's the same as the one entered at present
2035  //
2036  // in any case: exit the function if an entry for this column already
2037  // exists, since we don't want to enter it twice
2039  ExcInternalError());
2041  ExcColumnNotStoredHere(constrained_dof_index, column));
2042  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
2043  Assert(line_ptr->index == constrained_dof_index, ExcInternalError());
2044  for (const auto &p : line_ptr->entries)
2045  if (p.first == column)
2046  {
2047  Assert(std::abs(p.second - weight) < 1.e-14,
2049  constrained_dof_index, column, p.second, weight));
2050  return;
2051  }
2052 
2053  line_ptr->entries.emplace_back(column, weight);
2054 }
2055 
2056 
2057 
2058 template <typename number>
2059 inline void
2061  const size_type constrained_dof_index,
2062  const number value)
2063 {
2064  const size_type line_index = calculate_line_index(constrained_dof_index);
2065  Assert(line_index < lines_cache.size() &&
2066  lines_cache[line_index] != numbers::invalid_size_type,
2067  ExcMessage("call add_line() before calling set_inhomogeneity()"));
2068  Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
2069  ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
2070  line_ptr->inhomogeneity = value;
2071 }
2072 
2073 
2074 
2075 template <typename number>
2076 template <class VectorType>
2077 inline void
2079 {
2080  // since lines is a private member, we cannot pass it to the functions
2081  // above. therefore, copy the content which is cheap
2082  std::vector<size_type> constrained_lines(lines.size());
2083  for (unsigned int i = 0; i < lines.size(); ++i)
2084  constrained_lines[i] = lines[i].index;
2085  internal::AffineConstraintsImplementation::set_zero_all(constrained_lines,
2086  vec);
2087 }
2088 
2089 template <typename number>
2092 {
2093  return lines.size();
2094 }
2095 
2096 template <typename number>
2097 inline bool
2099 {
2100  const size_type line_index = calculate_line_index(index);
2101  return ((line_index < lines_cache.size()) &&
2102  (lines_cache[line_index] != numbers::invalid_size_type));
2103 }
2104 
2105 template <typename number>
2106 inline bool
2108  const size_type line_n) const
2109 {
2110  // check whether the entry is constrained. could use is_constrained, but
2111  // that means computing the line index twice
2112  const size_type line_index = calculate_line_index(line_n);
2113  if (line_index >= lines_cache.size() ||
2114  lines_cache[line_index] == numbers::invalid_size_type)
2115  return false;
2116  else
2117  {
2118  Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
2119  return !(lines[lines_cache[line_index]].inhomogeneity == number(0.));
2120  }
2121 }
2122 
2123 template <typename number>
2124 inline const std::vector<std::pair<types::global_dof_index, number>> *
2126 {
2127  // check whether the entry is constrained. could use is_constrained, but
2128  // that means computing the line index twice
2129  const size_type line_index = calculate_line_index(line_n);
2130  if (line_index >= lines_cache.size() ||
2131  lines_cache[line_index] == numbers::invalid_size_type)
2132  return nullptr;
2133  else
2134  return &lines[lines_cache[line_index]].entries;
2135 }
2136 
2137 template <typename number>
2138 inline number
2140 {
2141  // check whether the entry is constrained. could use is_constrained, but
2142  // that means computing the line index twice
2143  const size_type line_index = calculate_line_index(line_n);
2144  if (line_index >= lines_cache.size() ||
2145  lines_cache[line_index] == numbers::invalid_size_type)
2146  return 0;
2147  else
2148  return lines[lines_cache[line_index]].inhomogeneity;
2149 }
2150 
2151 template <typename number>
2154 {
2155  // IndexSet is unused (serial case)
2156  if (!local_lines.size())
2157  return line_n;
2158 
2160 
2161  return local_lines.index_within_set(line_n);
2162 }
2163 
2164 template <typename number>
2165 inline bool
2167 {
2168  return !local_lines.size() || local_lines.is_element(line_n);
2169 }
2170 
2171 template <typename number>
2172 inline const IndexSet &
2174 {
2175  return local_lines;
2176 }
2177 
2178 template <typename number>
2179 template <class VectorType>
2180 inline void
2182  const size_type index,
2183  const number value,
2184  VectorType & global_vector) const
2185 {
2186  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2187 
2188  if (is_constrained(index) == false)
2189  global_vector(index) += value;
2190  else
2191  {
2192  const ConstraintLine &position =
2194  for (size_type j = 0; j < position.entries.size(); ++j)
2195  global_vector(position.entries[j].first) +=
2196  value * position.entries[j].second;
2197  }
2198 }
2199 
2200 template <typename number>
2201 template <typename ForwardIteratorVec,
2202  typename ForwardIteratorInd,
2203  class VectorType>
2204 inline void
2206  ForwardIteratorVec local_vector_begin,
2207  ForwardIteratorVec local_vector_end,
2208  ForwardIteratorInd local_indices_begin,
2209  VectorType & global_vector) const
2210 {
2211  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2212  for (; local_vector_begin != local_vector_end;
2213  ++local_vector_begin, ++local_indices_begin)
2214  {
2215  if (is_constrained(*local_indices_begin) == false)
2216  internal::ElementAccess<VectorType>::add(*local_vector_begin,
2217  *local_indices_begin,
2218  global_vector);
2219  else
2220  {
2221  const ConstraintLine &position =
2222  lines[lines_cache[calculate_line_index(*local_indices_begin)]];
2223  for (size_type j = 0; j < position.entries.size(); ++j)
2225  (*local_vector_begin) * position.entries[j].second,
2226  position.entries[j].first,
2227  global_vector);
2228  }
2229  }
2230 }
2231 
2232 template <typename number>
2233 template <class InVector, class OutVector>
2234 inline void
2236  const InVector & local_vector,
2237  const std::vector<size_type> &local_dof_indices,
2238  OutVector & global_vector) const
2239 {
2240  Assert(local_vector.size() == local_dof_indices.size(),
2241  ExcDimensionMismatch(local_vector.size(), local_dof_indices.size()));
2242  distribute_local_to_global(local_vector.begin(),
2243  local_vector.end(),
2244  local_dof_indices.begin(),
2245  global_vector);
2246 }
2247 
2248 template <typename number>
2249 template <typename ForwardIteratorVec,
2250  typename ForwardIteratorInd,
2251  class VectorType>
2252 inline void
2254  const VectorType & global_vector,
2255  ForwardIteratorInd local_indices_begin,
2256  ForwardIteratorVec local_vector_begin,
2257  ForwardIteratorVec local_vector_end) const
2258 {
2259  Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2260  for (; local_vector_begin != local_vector_end;
2261  ++local_vector_begin, ++local_indices_begin)
2262  {
2263  if (is_constrained(*local_indices_begin) == false)
2264  *local_vector_begin = global_vector(*local_indices_begin);
2265  else
2266  {
2267  const ConstraintLine &position =
2268  lines[lines_cache[calculate_line_index(*local_indices_begin)]];
2269  typename VectorType::value_type value = position.inhomogeneity;
2270  for (size_type j = 0; j < position.entries.size(); ++j)
2271  value += (global_vector(position.entries[j].first) *
2272  position.entries[j].second);
2273  *local_vector_begin = value;
2274  }
2275  }
2276 }
2277 
2278 template <typename MatrixType>
2280 template <typename SparsityPatternType>
2282 template <typename number>
2284 
2285 namespace internal
2286 {
2288  {
2306  template <typename MatrixType>
2308  {
2309  private:
2315  template <typename T>
2316  static std::true_type
2317  check(const BlockMatrixBase<T> *);
2318 
2323  template <typename T>
2324  static std::true_type
2325  check(const BlockSparseMatrixEZ<T> *);
2326 
2331  static std::false_type
2332  check(...);
2333 
2334  public:
2341  static const bool value =
2342  std::is_same<decltype(check(std::declval<MatrixType *>())),
2343  std::true_type>::value;
2344  };
2345 
2346  // instantiation of the static member
2347  template <typename MatrixType>
2349 
2350 
2359  template <typename MatrixType>
2361  {
2362  private:
2367  template <typename T>
2368  static std::true_type
2369  check(const BlockSparsityPatternBase<T> *);
2370 
2375  static std::false_type
2376  check(...);
2377 
2378  public:
2384  static const bool value =
2385  std::is_same<decltype(check(std::declval<MatrixType *>())),
2386  std::true_type>::value;
2387  };
2388 
2389  // instantiation of the static member
2390  template <typename MatrixType>
2392 
2393  } // namespace AffineConstraints
2394 } // namespace internal
2395 
2396 
2397 template <typename number>
2398 template <typename MatrixType>
2399 inline void
2401  const FullMatrix<number> & local_matrix,
2402  const std::vector<size_type> &local_dof_indices,
2403  MatrixType & global_matrix) const
2404 {
2405  // create a dummy and hand on to the function actually implementing this
2406  // feature in the cm.templates.h file.
2407  Vector<typename MatrixType::value_type> dummy(0);
2409  local_matrix,
2410  dummy,
2411  local_dof_indices,
2412  global_matrix,
2413  dummy,
2414  false,
2415  std::integral_constant<
2416  bool,
2418 }
2419 
2420 
2421 
2422 template <typename number>
2423 template <typename MatrixType, typename VectorType>
2424 inline void
2426  const FullMatrix<number> & local_matrix,
2427  const Vector<number> & local_vector,
2428  const std::vector<size_type> &local_dof_indices,
2429  MatrixType & global_matrix,
2430  VectorType & global_vector,
2431  bool use_inhomogeneities_for_rhs) const
2432 {
2433  // enter the internal function with the respective block information set,
2434  // the actual implementation follows in the cm.templates.h file.
2436  local_matrix,
2437  local_vector,
2438  local_dof_indices,
2439  global_matrix,
2440  global_vector,
2441  use_inhomogeneities_for_rhs,
2442  std::integral_constant<
2443  bool,
2445 }
2446 
2447 
2448 
2449 template <typename number>
2450 template <typename SparsityPatternType>
2451 inline void
2453  const std::vector<size_type> &local_dof_indices,
2454  SparsityPatternType & sparsity_pattern,
2455  const bool keep_constrained_entries,
2456  const Table<2, bool> & dof_mask) const
2457 {
2458  // enter the internal function with the respective block information set,
2459  // the actual implementation follows in the cm.templates.h file.
2461  local_dof_indices,
2462  sparsity_pattern,
2463  keep_constrained_entries,
2464  dof_mask,
2465  std::integral_constant<bool,
2467  SparsityPatternType>::value>());
2468 }
2469 
2471 
2472 #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
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:1466
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1921
#define DeclException0(Exception0)
Definition: exceptions.h:470
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:380
constexpr TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
AffineConstraints(const IndexSet &local_constraints=IndexSet())
void add_entries_local_to_global(const std::vector< size_type > &local_dof_indices, SparsityPatternType &sparsity_pattern, const bool keep_constrained_entries=true, const Table< 2, bool > &dof_mask=Table< 2, bool >()) const
size_type calculate_line_index(const size_type line_n) const
std::vector< std::pair< size_type, number > > Entries
const IndexSet & get_local_lines() const
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
void compress() const
Definition: index_set.h:1642
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:379
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:584
Definition: table.h:688
bool is_element(const size_type index) const
Definition: index_set.h:1765
void add_line(const size_type line_n)
typename std::vector< ConstraintLine >::const_iterator const_iterator
T max(const T &t, const MPI_Comm &mpi_communicator)
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:931
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()