Reference documentation for deal.II version 9.4.0
\(\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 - 2022 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
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
43template <typename>
44class FullMatrix;
45class SparsityPattern;
49template <typename number>
50class SparseMatrix;
51template <typename number>
53
54namespace 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
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
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
195 size() const;
196
202 size(const size_type counter_index) const;
203
208 global_row(const size_type counter_index) const;
209
213 size_type &
214 global_row(const size_type counter_index);
215
222 local_row(const size_type counter_index) const;
223
227 size_type &
228 local_row(const size_type counter_index);
229
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
268 n_constraints() const;
269
275 n_inhomogeneities() const;
276
283 void
284 set_ith_constraint_inhomogeneous(const size_type i);
285
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
398namespace 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#endif
419
420// TODO[WB]: We should have a function of the kind
421// AffineConstraints::add_constraint (const size_type constrained_dof,
422// const std::vector<std::pair<size_type, number> > &entries,
423// const number inhomogeneity = 0);
424// rather than building up constraints piecemeal through add_line/add_entry
425// etc. This would also eliminate the possibility of accidentally changing
426// existing constraints into something pointless, see the discussion on the
427// mailing list on "Tiny bug in interpolate_boundary_values" in Sept. 2010.
428
502template <typename number = double>
504{
505public:
510
517 {
523
530
537 };
538
556 explicit AffineConstraints(const IndexSet &local_constraints = IndexSet());
557
561 explicit AffineConstraints(const AffineConstraints &affine_constraints);
562
566 AffineConstraints(AffineConstraints &&affine_constraints) noexcept =
567 default; // NOLINT
568
579 operator=(const AffineConstraints &) = delete;
580
585 operator=(AffineConstraints &&affine_constraints) noexcept =
586 default; // NOLINT
587
594 template <typename other_number>
595 void
597
604 void
605 reinit(const IndexSet &local_constraints = IndexSet());
606
613 bool
614 can_store_line(const size_type line_n) const;
615
622 const IndexSet &
624
644 void
646 const IndexSet & filter);
647
657 void
658 add_line(const size_type line_n);
659
672 void
673 add_lines(const std::vector<bool> &lines);
674
687 void
688 add_lines(const std::set<size_type> &lines);
689
702 void
704
720 void
721 add_entry(const size_type constrained_dof_index,
722 const size_type column,
723 const number weight);
724
730 void
732 const size_type constrained_dof_index,
733 const std::vector<std::pair<size_type, number>> &col_weight_pairs);
734
746 void
747 set_inhomogeneity(const size_type constrained_dof_index, const number value);
748
770 void
772
778 bool
779 is_closed() const;
780
786 bool
787 is_closed(const MPI_Comm &comm) const;
788
813 void
815 const AffineConstraints & other_constraints,
816 const MergeConflictBehavior merge_conflict_behavior = no_conflicts_allowed,
817 const bool allow_different_local_lines = false);
818
831 void
832 shift(const size_type offset);
833
841 void
843
858
865
873
883 bool
884 is_constrained(const size_type line_n) const;
885
897 bool
899
906 bool
908 const size_type line_n_2) const;
909
922
927 bool
929
935 bool
937
942 const std::vector<std::pair<size_type, number>> *
944
949 number
950 get_inhomogeneity(const size_type line_n) const;
951
972 void
973 print(std::ostream &out) const;
974
987 void
988 write_dot(std::ostream &) const;
989
994 std::size_t
996
1003 void
1004 resolve_indices(std::vector<types::global_dof_index> &indices) const;
1005
1031 void
1032 condense(SparsityPattern &sparsity) const;
1033
1037 void
1039
1044 void
1046
1051 void
1053
1061 void
1063
1067 void
1069
1081 template <class VectorType>
1082 void
1083 condense(VectorType &vec) const;
1084
1091 template <class VectorType>
1092 void
1093 condense(const VectorType &vec_ghosted, VectorType &output) const;
1094
1107 template <class VectorType>
1108 void
1109 condense(SparseMatrix<number> &matrix, VectorType &vector) const;
1110
1115 template <class BlockVectorType>
1116 void
1117 condense(BlockSparseMatrix<number> &matrix, BlockVectorType &vector) const;
1118
1125 template <class VectorType>
1126 void
1127 set_zero(VectorType &vec) const;
1128
1184 template <class InVector, class OutVector>
1185 void
1186 distribute_local_to_global(const InVector & local_vector,
1187 const std::vector<size_type> &local_dof_indices,
1188 OutVector & global_vector) const;
1189
1237 template <typename VectorType>
1238 void
1240 const std::vector<size_type> &local_dof_indices,
1241 VectorType & global_vector,
1242 const FullMatrix<number> & local_matrix) const;
1243
1263 template <typename VectorType>
1264 void
1266 const Vector<number> & local_vector,
1267 const std::vector<size_type> &local_dof_indices_row,
1268 const std::vector<size_type> &local_dof_indices_col,
1269 VectorType & global_vector,
1270 const FullMatrix<number> & local_matrix,
1271 bool diagonal = false) const;
1272
1276 template <class VectorType>
1277 void
1279 const number value,
1280 VectorType & global_vector) const;
1281
1314 template <typename ForwardIteratorVec,
1315 typename ForwardIteratorInd,
1316 class VectorType>
1317 void
1318 distribute_local_to_global(ForwardIteratorVec local_vector_begin,
1319 ForwardIteratorVec local_vector_end,
1320 ForwardIteratorInd local_indices_begin,
1321 VectorType & global_vector) const;
1322
1374 template <typename MatrixType>
1375 void
1377 const std::vector<size_type> &local_dof_indices,
1378 MatrixType & global_matrix) const;
1379
1407 template <typename MatrixType>
1408 void
1410 const std::vector<size_type> &row_indices,
1411 const std::vector<size_type> &col_indices,
1412 MatrixType & global_matrix) const;
1413
1430 template <typename MatrixType>
1431 void
1433 const std::vector<size_type> &row_indices,
1434 const AffineConstraints &column_affine_constraints,
1435 const std::vector<size_type> &column_indices,
1436 MatrixType & global_matrix) const;
1437
1458 template <typename MatrixType, typename VectorType>
1459 void
1461 const Vector<number> & local_vector,
1462 const std::vector<size_type> &local_dof_indices,
1463 MatrixType & global_matrix,
1464 VectorType & global_vector,
1465 bool use_inhomogeneities_for_rhs = false) const;
1466
1520 template <typename SparsityPatternType>
1521 void
1523 const std::vector<size_type> &local_dof_indices,
1524 SparsityPatternType & sparsity_pattern,
1525 const bool keep_constrained_entries = true,
1526 const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1527
1531 template <typename SparsityPatternType>
1532 void
1534 const std::vector<size_type> &row_indices,
1535 const std::vector<size_type> &col_indices,
1536 SparsityPatternType & sparsity_pattern,
1537 const bool keep_constrained_entries = true,
1538 const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1539
1544 template <typename SparsityPatternType>
1545 void
1547 const std::vector<size_type> & row_indices,
1548 const AffineConstraints<number> &col_constraints,
1549 const std::vector<size_type> & col_indices,
1550 SparsityPatternType & sparsity_pattern,
1551 const bool keep_constrained_entries = true,
1552 const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1553
1573 template <typename ForwardIteratorVec,
1574 typename ForwardIteratorInd,
1575 class VectorType>
1576 void
1577 get_dof_values(const VectorType & global_vector,
1578 ForwardIteratorInd local_indices_begin,
1579 ForwardIteratorVec local_vector_begin,
1580 ForwardIteratorVec local_vector_end) const;
1581
1603 template <class VectorType>
1604 void
1605 distribute(VectorType &vec) const;
1606
1615 {
1620 using Entries = std::vector<std::pair<size_type, number>>;
1621
1628
1637
1642
1647 const Entries & entries = {},
1648 const number & inhomogeneity = 0.0);
1649
1653 template <typename ConstraintLineType>
1654 ConstraintLine(const ConstraintLineType &other);
1655
1659 template <typename ConstraintLineType>
1661 operator=(const ConstraintLineType &other);
1662
1670 bool
1671 operator<(const ConstraintLine &) const;
1672
1678 bool
1680
1685 std::size_t
1687
1693 template <class Archive>
1694 void
1695 serialize(Archive &ar, const unsigned int)
1696 {
1698 }
1699 };
1700
1704 using const_iterator = typename std::vector<ConstraintLine>::const_iterator;
1705
1709 using LineRange = boost::iterator_range<const_iterator>;
1710
1719 const LineRange
1720 get_lines() const;
1721
1750 bool
1751 is_consistent_in_parallel(const std::vector<IndexSet> &locally_owned_dofs,
1752 const IndexSet & locally_active_dofs,
1753 const MPI_Comm & mpi_communicator,
1754 const bool verbose = false) const;
1755
1761 void
1762 make_consistent_in_parallel(const IndexSet &locally_owned_dofs,
1763 const IndexSet &locally_relevant_dofs,
1764 const MPI_Comm mpi_communicator);
1765
1784 size_type,
1785 << "The specified line " << arg1 << " does not exist.");
1792 size_type,
1793 size_type,
1794 number,
1795 number,
1796 << "The entry for the indices " << arg1 << " and " << arg2
1797 << " already exists, but the values " << arg3 << " (old) and "
1798 << arg4 << " (new) differ "
1799 << "by " << (arg4 - arg3) << '.');
1806 int,
1807 int,
1808 << "You tried to constrain DoF " << arg1 << " to DoF " << arg2
1809 << ", but that one is also constrained. This is not allowed!");
1816 size_type,
1817 << "Degree of freedom " << arg1
1818 << " is constrained from both object in a merge operation.");
1825 size_type,
1826 << "In the given argument a degree of freedom is constrained "
1827 << "to another DoF with number " << arg1
1828 << ", which however is constrained by this object. This is not"
1829 << " allowed.");
1836 size_type,
1837 << "The index set given to this constraints object indicates "
1838 << "constraints for degree of freedom " << arg1
1839 << " should not be stored by this object, but a constraint "
1840 << "is being added.");
1841
1848 size_type,
1849 size_type,
1850 << "The index set given to this constraints object indicates "
1851 << "constraints using degree of freedom " << arg2
1852 << " should not be stored by this object, but a constraint "
1853 << "for degree of freedom " << arg1 << " uses it.");
1854
1861 int,
1862 int,
1863 << "While distributing the constraint for DoF " << arg1
1864 << ", it turns out that one of the processors "
1865 << "who own the " << arg2 << " degrees of freedom that x_"
1866 << arg1 << " is constrained against does not know about "
1867 << "the constraint on x_" << arg1
1868 << ". Did you not initialize the AffineConstraints container "
1869 << "with the appropriate locally_relevant set so "
1870 << "that every processor who owns a DoF that constrains "
1871 << "another DoF also knows about this constraint?");
1872
1873 template <typename>
1874 friend class AffineConstraints;
1875
1876private:
1888 std::vector<ConstraintLine> lines;
1889
1922 std::vector<size_type> lines_cache;
1923
1930
1935
1937 internal::AffineConstraints::ScratchData<number>>
1939
1944 size_type
1945 calculate_line_index(const size_type line_n) const;
1946
1951 template <typename MatrixType, typename VectorType>
1952 void
1954 const Vector<number> & local_vector,
1955 const std::vector<size_type> &local_dof_indices,
1956 MatrixType & global_matrix,
1957 VectorType & global_vector,
1958 const bool use_inhomogeneities_for_rhs,
1959 const std::integral_constant<bool, false>) const;
1960
1965 template <typename MatrixType, typename VectorType>
1966 void
1968 const Vector<number> & local_vector,
1969 const std::vector<size_type> &local_dof_indices,
1970 MatrixType & global_matrix,
1971 VectorType & global_vector,
1972 const bool use_inhomogeneities_for_rhs,
1973 const std::integral_constant<bool, true>) const;
1974
1979 template <typename SparsityPatternType>
1980 void
1981 add_entries_local_to_global(const std::vector<size_type> &local_dof_indices,
1982 SparsityPatternType & sparsity_pattern,
1983 const bool keep_constrained_entries,
1984 const Table<2, bool> &dof_mask,
1985 const std::integral_constant<bool, false>) const;
1986
1991 template <typename SparsityPatternType>
1992 void
1993 add_entries_local_to_global(const std::vector<size_type> &local_dof_indices,
1994 SparsityPatternType & sparsity_pattern,
1995 const bool keep_constrained_entries,
1996 const Table<2, bool> &dof_mask,
1997 const std::integral_constant<bool, true>) const;
1998
2006 void
2007 make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
2008 internal::AffineConstraints::GlobalRowsFromLocal<number>
2009 &global_rows) const;
2010
2018 void
2019 make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
2020 std::vector<size_type> & active_dofs) const;
2021
2025 template <typename MatrixScalar, typename VectorScalar>
2028 const size_type i,
2029 const internal::AffineConstraints::GlobalRowsFromLocal<number> &global_rows,
2030 const Vector<VectorScalar> & local_vector,
2031 const std::vector<size_type> & local_dof_indices,
2032 const FullMatrix<MatrixScalar> &local_matrix) const;
2033};
2034
2035/* ---------------- template and inline functions ----------------- */
2036
2037template <typename number>
2039 const IndexSet &local_constraints)
2040 : lines()
2041 , local_lines(local_constraints)
2042 , sorted(false)
2043{
2044 // make sure the IndexSet is compressed. Otherwise this can lead to crashes
2045 // that are hard to find (only happen in release mode).
2046 // see tests/mpi/affine_constraints_crash_01
2048}
2049
2050template <typename number>
2052 const AffineConstraints &affine_constraints)
2053 : Subscriptor()
2054 , lines(affine_constraints.lines)
2055 , lines_cache(affine_constraints.lines_cache)
2056 , local_lines(affine_constraints.local_lines)
2057 , sorted(affine_constraints.sorted)
2058{}
2059
2060template <typename number>
2061inline void
2063{
2064 Assert(sorted == false, ExcMatrixIsClosed());
2065
2066 // the following can happen when we compute with distributed meshes and dof
2067 // handlers and we constrain a degree of freedom whose number we don't have
2068 // locally. if we don't abort here the program will try to allocate several
2069 // terabytes of memory to resize the various arrays below :-)
2071 const size_type line_index = calculate_line_index(line_n);
2072
2073 // check whether line already exists; it may, in which case we can just quit
2074 if (is_constrained(line_n))
2075 return;
2076
2077 // if necessary enlarge vector of existing entries for cache
2078 if (line_index >= lines_cache.size())
2079 lines_cache.resize(std::max(2 * static_cast<size_type>(lines_cache.size()),
2080 line_index + 1),
2082
2083 // push a new line to the end of the list
2084 lines.emplace_back();
2085 lines.back().index = line_n;
2086 lines.back().inhomogeneity = 0.;
2087 lines_cache[line_index] = lines.size() - 1;
2088}
2089
2090
2091
2092template <typename number>
2093inline void
2095 const size_type column,
2096 const number weight)
2097{
2098 Assert(sorted == false, ExcMatrixIsClosed());
2099 Assert(constrained_dof_index != column,
2100 ExcMessage("Can't constrain a degree of freedom to itself"));
2101
2102 // Ensure that the current line is present in the cache:
2103 const size_type line_index = calculate_line_index(constrained_dof_index);
2104 Assert(line_index < lines_cache.size(),
2105 ExcMessage("The current AffineConstraints does not contain the line "
2106 "for the current entry. Call AffineConstraints::add_line "
2107 "before calling this function."));
2108
2109 // if in debug mode, check whether an entry for this column already exists
2110 // and if it's the same as the one entered at present
2111 //
2112 // in any case: exit the function if an entry for this column already
2113 // exists, since we don't want to enter it twice
2117 ExcColumnNotStoredHere(constrained_dof_index, column));
2118 ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
2119 Assert(line_ptr->index == constrained_dof_index, ExcInternalError());
2120 for (const auto &p : line_ptr->entries)
2121 if (p.first == column)
2122 {
2123 Assert(std::abs(p.second - weight) < 1.e-14,
2125 constrained_dof_index, column, p.second, weight));
2126 return;
2127 }
2128
2129 line_ptr->entries.emplace_back(column, weight);
2130}
2131
2132
2133
2134template <typename number>
2135inline void
2137 const size_type constrained_dof_index,
2138 const number value)
2139{
2140 const size_type line_index = calculate_line_index(constrained_dof_index);
2141 Assert(line_index < lines_cache.size() &&
2143 ExcMessage("call add_line() before calling set_inhomogeneity()"));
2144 Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
2145 ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
2146 line_ptr->inhomogeneity = value;
2147}
2148
2149
2150
2151template <typename number>
2152template <class VectorType>
2153inline void
2155{
2156 // since lines is a private member, we cannot pass it to the functions
2157 // above. therefore, copy the content which is cheap
2158 std::vector<size_type> constrained_lines(lines.size());
2159 for (unsigned int i = 0; i < lines.size(); ++i)
2160 constrained_lines[i] = lines[i].index;
2161 internal::AffineConstraintsImplementation::set_zero_all(constrained_lines,
2162 vec);
2163}
2164
2165template <typename number>
2168{
2169 return lines.size();
2170}
2171
2172template <typename number>
2175{
2176 return std::count_if(lines.begin(),
2177 lines.end(),
2178 [](const ConstraintLine &line) {
2179 return (line.entries.size() == 1) &&
2180 (line.entries[0].second == number(1.));
2181 });
2182}
2183
2184template <typename number>
2187{
2188 return std::count_if(lines.begin(),
2189 lines.end(),
2190 [](const ConstraintLine &line) {
2191 return (line.inhomogeneity != number(0.));
2192 });
2193}
2194
2195template <typename number>
2196inline bool
2198{
2199 const size_type line_index = calculate_line_index(index);
2200 return ((line_index < lines_cache.size()) &&
2201 (lines_cache[line_index] != numbers::invalid_size_type));
2202}
2203
2204template <typename number>
2205inline bool
2207 const size_type line_n) const
2208{
2209 // check whether the entry is constrained. could use is_constrained, but
2210 // that means computing the line index twice
2211 const size_type line_index = calculate_line_index(line_n);
2212 if (line_index >= lines_cache.size() ||
2214 return false;
2215 else
2216 {
2217 Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
2218 return (lines[lines_cache[line_index]].inhomogeneity != number(0.));
2219 }
2220}
2221
2222template <typename number>
2223inline const std::vector<std::pair<types::global_dof_index, number>> *
2225{
2226 // check whether the entry is constrained. could use is_constrained, but
2227 // that means computing the line index twice
2228 const size_type line_index = calculate_line_index(line_n);
2229 if (line_index >= lines_cache.size() ||
2231 return nullptr;
2232 else
2233 return &lines[lines_cache[line_index]].entries;
2234}
2235
2236template <typename number>
2237inline number
2239{
2240 // check whether the entry is constrained. could use is_constrained, but
2241 // that means computing the line index twice
2242 const size_type line_index = calculate_line_index(line_n);
2243 if (line_index >= lines_cache.size() ||
2245 return 0;
2246 else
2247 return lines[lines_cache[line_index]].inhomogeneity;
2248}
2249
2250template <typename number>
2253{
2254 // IndexSet is unused (serial case)
2255 if (local_lines.size() == 0)
2256 return line_n;
2257
2259
2260 return local_lines.index_within_set(line_n);
2261}
2262
2263template <typename number>
2264inline bool
2266{
2267 return local_lines.size() == 0 || local_lines.is_element(line_n);
2268}
2269
2270template <typename number>
2271inline const IndexSet &
2273{
2274 return local_lines;
2275}
2276
2277template <typename number>
2278template <class VectorType>
2279inline void
2281 const size_type index,
2282 const number value,
2283 VectorType & global_vector) const
2284{
2285 Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2286
2287 if (is_constrained(index) == false)
2288 global_vector(index) += value;
2289 else
2290 {
2291 const ConstraintLine &position =
2293 for (size_type j = 0; j < position.entries.size(); ++j)
2294 global_vector(position.entries[j].first) +=
2295 value * position.entries[j].second;
2296 }
2297}
2298
2299template <typename number>
2300template <typename ForwardIteratorVec,
2301 typename ForwardIteratorInd,
2302 class VectorType>
2303inline void
2305 ForwardIteratorVec local_vector_begin,
2306 ForwardIteratorVec local_vector_end,
2307 ForwardIteratorInd local_indices_begin,
2308 VectorType & global_vector) const
2309{
2310 Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2311 for (; local_vector_begin != local_vector_end;
2312 ++local_vector_begin, ++local_indices_begin)
2313 {
2314 if (is_constrained(*local_indices_begin) == false)
2315 internal::ElementAccess<VectorType>::add(*local_vector_begin,
2316 *local_indices_begin,
2317 global_vector);
2318 else
2319 {
2320 const ConstraintLine &position =
2321 lines[lines_cache[calculate_line_index(*local_indices_begin)]];
2322 for (size_type j = 0; j < position.entries.size(); ++j)
2324 (*local_vector_begin) * position.entries[j].second,
2325 position.entries[j].first,
2326 global_vector);
2327 }
2328 }
2329}
2330
2331template <typename number>
2332template <class InVector, class OutVector>
2333inline void
2335 const InVector & local_vector,
2336 const std::vector<size_type> &local_dof_indices,
2337 OutVector & global_vector) const
2338{
2339 Assert(local_vector.size() == local_dof_indices.size(),
2340 ExcDimensionMismatch(local_vector.size(), local_dof_indices.size()));
2341 distribute_local_to_global(local_vector.begin(),
2342 local_vector.end(),
2343 local_dof_indices.begin(),
2344 global_vector);
2345}
2346
2347template <typename number>
2348template <typename ForwardIteratorVec,
2349 typename ForwardIteratorInd,
2350 class VectorType>
2351inline void
2353 const VectorType & global_vector,
2354 ForwardIteratorInd local_indices_begin,
2355 ForwardIteratorVec local_vector_begin,
2356 ForwardIteratorVec local_vector_end) const
2357{
2358 Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2359 for (; local_vector_begin != local_vector_end;
2360 ++local_vector_begin, ++local_indices_begin)
2361 {
2362 if (is_constrained(*local_indices_begin) == false)
2363 *local_vector_begin = global_vector(*local_indices_begin);
2364 else
2365 {
2366 const ConstraintLine &position =
2367 lines[lines_cache[calculate_line_index(*local_indices_begin)]];
2368 typename VectorType::value_type value = position.inhomogeneity;
2369 for (size_type j = 0; j < position.entries.size(); ++j)
2370 value += (global_vector(position.entries[j].first) *
2371 position.entries[j].second);
2372 *local_vector_begin = value;
2373 }
2374 }
2375}
2376
2377// Forward declarations
2378#ifndef DOXYGEN
2379template <typename MatrixType>
2380class BlockMatrixBase;
2381template <typename SparsityPatternType>
2383template <typename number>
2385
2386namespace internal
2387{
2388 namespace AffineConstraints
2389 {
2407 template <typename MatrixType>
2408 struct IsBlockMatrix
2409 {
2410 private:
2416 template <typename T>
2417 static std::true_type
2418 check(const BlockMatrixBase<T> *);
2419
2424 template <typename T>
2425 static std::true_type
2427
2432 static std::false_type
2433 check(...);
2434
2435 public:
2442 static const bool value =
2443 std::is_same<decltype(check(std::declval<MatrixType *>())),
2444 std::true_type>::value;
2445 };
2446
2447 // instantiation of the static member
2448 template <typename MatrixType>
2449 const bool IsBlockMatrix<MatrixType>::value;
2450
2451
2460 template <typename MatrixType>
2461 struct IsBlockSparsityPattern
2462 {
2463 private:
2468 template <typename T>
2469 static std::true_type
2471
2476 static std::false_type
2477 check(...);
2478
2479 public:
2485 static const bool value =
2486 std::is_same<decltype(check(std::declval<MatrixType *>())),
2487 std::true_type>::value;
2488 };
2489
2490 // instantiation of the static member
2491 template <typename MatrixType>
2492 const bool IsBlockSparsityPattern<MatrixType>::value;
2493
2494 } // namespace AffineConstraints
2495} // namespace internal
2496#endif
2497
2498
2499
2500template <typename number>
2501template <typename other_number>
2502inline void
2505{
2506 lines.clear();
2507 lines.insert(lines.begin(), other.lines.begin(), other.lines.end());
2508 lines_cache = other.lines_cache;
2509 local_lines = other.local_lines;
2510 sorted = other.sorted;
2511}
2512
2513
2514
2515template <typename number>
2516template <typename MatrixType>
2517inline void
2519 const FullMatrix<number> & local_matrix,
2520 const std::vector<size_type> &local_dof_indices,
2521 MatrixType & global_matrix) const
2522{
2523 // create a dummy and hand on to the function actually implementing this
2524 // feature in the cm.templates.h file.
2526 distribute_local_to_global(
2527 local_matrix,
2528 dummy,
2529 local_dof_indices,
2530 global_matrix,
2531 dummy,
2532 false,
2533 std::integral_constant<
2534 bool,
2535 internal::AffineConstraints::IsBlockMatrix<MatrixType>::value>());
2536}
2537
2538
2539
2540template <typename number>
2541template <typename MatrixType, typename VectorType>
2542inline void
2544 const FullMatrix<number> & local_matrix,
2545 const Vector<number> & local_vector,
2546 const std::vector<size_type> &local_dof_indices,
2547 MatrixType & global_matrix,
2548 VectorType & global_vector,
2549 bool use_inhomogeneities_for_rhs) const
2550{
2551 // enter the internal function with the respective block information set,
2552 // the actual implementation follows in the cm.templates.h file.
2553 distribute_local_to_global(
2554 local_matrix,
2555 local_vector,
2556 local_dof_indices,
2557 global_matrix,
2558 global_vector,
2559 use_inhomogeneities_for_rhs,
2560 std::integral_constant<
2561 bool,
2562 internal::AffineConstraints::IsBlockMatrix<MatrixType>::value>());
2563}
2564
2565
2566
2567template <typename number>
2568template <typename SparsityPatternType>
2569inline void
2571 const std::vector<size_type> &local_dof_indices,
2572 SparsityPatternType & sparsity_pattern,
2573 const bool keep_constrained_entries,
2574 const Table<2, bool> & dof_mask) const
2575{
2576 // enter the internal function with the respective block information set,
2577 // the actual implementation follows in the cm.templates.h file.
2578 add_entries_local_to_global(
2579 local_dof_indices,
2580 sparsity_pattern,
2581 keep_constrained_entries,
2582 dof_mask,
2583 std::integral_constant<bool,
2584 internal::AffineConstraints::IsBlockSparsityPattern<
2585 SparsityPatternType>::value>());
2586}
2587
2588
2589
2590template <typename number>
2592 const size_type & index,
2594 const number &inhomogeneity)
2595 : index(index)
2596 , entries(entries)
2597 , inhomogeneity(inhomogeneity)
2598{}
2599
2600
2601
2602template <typename number>
2603template <typename ConstraintLineType>
2605 const ConstraintLineType &other)
2606{
2607 this->index = other.index;
2608
2609 entries.clear();
2610 entries.insert(entries.begin(), other.entries.begin(), other.entries.end());
2611
2612 this->inhomogeneity = other.inhomogeneity;
2613}
2614
2615
2616
2617template <typename number>
2618template <typename ConstraintLineType>
2621 const ConstraintLineType &other)
2622{
2623 this->index = other.index;
2624
2625 entries.clear();
2626 entries.insert(entries.begin(), other.entries.begin(), other.entries.end());
2627
2628 this->inhomogeneity = other.inhomogeneity;
2629
2630 return *this;
2631}
2632
2634
2635#endif
bool is_closed(const MPI_Comm &comm) const
const LineRange get_lines() const
void add_entries(const size_type constrained_dof_index, const std::vector< std::pair< size_type, number > > &col_weight_pairs)
void condense(DynamicSparsityPattern &sparsity) const
void distribute_local_to_global(ForwardIteratorVec local_vector_begin, ForwardIteratorVec local_vector_end, ForwardIteratorInd local_indices_begin, VectorType &global_vector) const
bool is_closed() const
void merge(const AffineConstraints &other_constraints, const MergeConflictBehavior merge_conflict_behavior=no_conflicts_allowed, const bool allow_different_local_lines=false)
size_type n_inhomogeneities() const
bool has_inhomogeneities() const
void add_line(const size_type line_n)
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::integral_constant< bool, true >) const
void add_selected_constraints(const AffineConstraints &constraints_in, const IndexSet &filter)
void condense(VectorType &vec) const
bool can_store_line(const size_type line_n) const
AffineConstraints(AffineConstraints &&affine_constraints) noexcept=default
Threads::ThreadLocalStorage< internal::AffineConstraints::ScratchData< number > > scratch_data
friend class AffineConstraints
void add_entry(const size_type constrained_dof_index, const size_type column, const number weight)
void reinit(const IndexSet &local_constraints=IndexSet())
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 shift(const size_type offset)
size_type n_identities() const
std::size_t memory_consumption() const
void distribute_local_to_global(const Vector< number > &local_vector, const std::vector< size_type > &local_dof_indices_row, const std::vector< size_type > &local_dof_indices_col, VectorType &global_vector, const FullMatrix< number > &local_matrix, bool diagonal=false) const
typename std::vector< ConstraintLine >::const_iterator const_iterator
void condense(SparseMatrix< number > &matrix, VectorType &vector) const
void set_inhomogeneity(const size_type constrained_dof_index, const number value)
void distribute_local_to_global(const FullMatrix< number > &local_matrix, const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, MatrixType &global_matrix) const
void condense(BlockSparsityPattern &sparsity) 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::integral_constant< bool, 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
void condense(SparsityPattern &sparsity) const
void add_entries_local_to_global(const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, SparsityPatternType &sparsity_pattern, const bool keep_constrained_entries=true, const Table< 2, bool > &dof_mask=Table< 2, bool >()) const
void add_entries_local_to_global(const std::vector< size_type > &row_indices, const AffineConstraints< number > &col_constraints, const std::vector< size_type > &col_indices, SparsityPatternType &sparsity_pattern, const bool keep_constrained_entries=true, const Table< 2, bool > &dof_mask=Table< 2, bool >()) const
void distribute_local_to_global(const size_type index, const number value, VectorType &global_vector) const
size_type max_constraint_indirections() const
void add_entries_local_to_global(const std::vector< size_type > &local_dof_indices, SparsityPatternType &sparsity_pattern, const bool keep_constrained_entries, const Table< 2, bool > &dof_mask, const std::integral_constant< bool, false >) const
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
void add_lines(const std::set< size_type > &lines)
AffineConstraints(const AffineConstraints &affine_constraints)
void distribute(VectorType &vec) const
void add_lines(const IndexSet &lines)
void condense(BlockSparseMatrix< number > &matrix) 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
void condense(BlockDynamicSparsityPattern &sparsity) const
void distribute_local_to_global(const FullMatrix< number > &local_matrix, const std::vector< size_type > &local_dof_indices, MatrixType &global_matrix) const
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
AffineConstraints & operator=(AffineConstraints &&affine_constraints) noexcept=default
void condense(SparseMatrix< number > &matrix) const
void add_entries_local_to_global(const std::vector< size_type > &local_dof_indices, SparsityPatternType &sparsity_pattern, const bool keep_constrained_entries, const Table< 2, bool > &dof_mask, const std::integral_constant< bool, true >) const
bool is_inhomogeneously_constrained(const size_type index) 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, bool use_inhomogeneities_for_rhs=false) const
void condense(const VectorType &vec_ghosted, VectorType &output) const
void copy_from(const AffineConstraints< other_number > &other)
std::vector< size_type > lines_cache
AffineConstraints(const IndexSet &local_constraints=IndexSet())
void condense(BlockSparseMatrix< number > &matrix, BlockVectorType &vector) const
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
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 distribute_local_to_global(const Vector< number > &local_vector, const std::vector< size_type > &local_dof_indices, VectorType &global_vector, const FullMatrix< number > &local_matrix) const
void set_zero(VectorType &vec) const
void distribute_local_to_global(const FullMatrix< number > &local_matrix, const std::vector< size_type > &row_indices, const AffineConstraints &column_affine_constraints, const std::vector< size_type > &column_indices, MatrixType &global_matrix) 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
AffineConstraints & operator=(const AffineConstraints &)=delete
size_type size() const
Definition: index_set.h:1636
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1923
bool is_element(const size_type index) const
Definition: index_set.h:1767
void compress() const
Definition: index_set.h:1644
A class that provides a separate storage location on each thread that accesses the object.
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define DeclException0(Exception0)
Definition: exceptions.h:464
static ::ExceptionBase & ExcMatrixNotClosed()
static ::ExceptionBase & ExcLineInexistant(size_type arg1)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:578
static ::ExceptionBase & ExcDoFIsConstrainedFromBothObjects(size_type arg1)
static ::ExceptionBase & ExcDoFConstrainedToConstrainedDoF(int arg1, int arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcMatrixIsClosed()
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:532
static ::ExceptionBase & ExcDoFIsConstrainedToConstrainedDoF(size_type arg1)
static ::ExceptionBase & ExcColumnNotStoredHere(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcRowNotStoredHere(size_type arg1)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcIncorrectConstraint(int arg1, int arg2)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
static ::ExceptionBase & ExcEntryAlreadyExists(size_type arg1, size_type arg2, number arg3, number arg4)
static ::ExceptionBase & ExcMessage(std::string arg1)
@ matrix
Contents is actually a matrix.
@ diagonal
Matrix is diagonal.
types::global_dof_index size_type
Definition: cuda_kernels.h:45
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:216
const types::global_dof_index invalid_size_type
Definition: types.h:210
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
Definition: types.h:76
bool operator<(const ConstraintLine &) const
void serialize(Archive &ar, const unsigned int)
ConstraintLine & operator=(const ConstraintLineType &other)
ConstraintLine(const size_type &index=numbers::invalid_dof_index, const Entries &entries={}, const number &inhomogeneity=0.0)
std::vector< std::pair< size_type, number > > Entries
std::size_t memory_consumption() const
bool operator==(const ConstraintLine &) const
typename internal::ProductTypeImpl< typename std::decay< T >::type, typename std::decay< U >::type >::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