Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
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
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
44template <typename>
45class FullMatrix;
46class SparsityPattern;
50template <typename number>
51class SparseMatrix;
52template <typename number>
54
55namespace 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
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
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
196 size() const;
197
203 size(const size_type counter_index) const;
204
209 global_row(const size_type counter_index) const;
210
214 size_type &
215 global_row(const size_type counter_index);
216
223 local_row(const size_type counter_index) const;
224
228 size_type &
229 local_row(const size_type counter_index);
230
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
253 bool
254 have_indirect_rows() const;
255
260 void
261 insert_constraint(const size_type constrained_local_dof);
262
269 n_constraints() const;
270
276 n_inhomogeneities() const;
277
284 void
285 set_ith_constraint_inhomogeneous(const size_type i);
286
292 constraint_origin(size_type i) const;
293
299 std::vector<Distributing> total_row_indices;
300
301 private:
305 DataCache<number> data_cache;
306
311 size_type n_active_rows;
312
317 size_type n_inhomogeneous_rows;
318 };
319
320
321
339 template <typename number>
340 struct ScratchData
341 {
345 ScratchData()
346 : in_use(false)
347 {}
348
352 ScratchData(const ScratchData &)
353 : in_use(false)
354 {}
355
359 bool in_use;
360
364 std::vector<std::pair<size_type, size_type>> new_entries;
365
369 std::vector<size_type> rows;
370
374 std::vector<size_type> columns;
375
379 std::vector<number> values;
380
384 std::vector<size_type> block_starts;
385
389 std::vector<size_type> vector_indices;
390
394 std::vector<number> vector_values;
395
399 GlobalRowsFromLocal<number> global_rows;
400
404 GlobalRowsFromLocal<number> global_columns;
405 };
406 } // namespace AffineConstraints
407} // namespace internal
408
409namespace internal
410{
411 namespace AffineConstraintsImplementation
412 {
413 template <class VectorType>
414 void
415 set_zero_all(const std::vector<types::global_dof_index> &cm,
416 VectorType & vec);
417
418 template <class T>
419 void
420 set_zero_all(const std::vector<types::global_dof_index> &cm,
421 ::Vector<T> & vec);
422
423 template <class T>
424 void
425 set_zero_all(const std::vector<types::global_dof_index> &cm,
426 ::BlockVector<T> & vec);
427 } // namespace AffineConstraintsImplementation
428} // namespace internal
429#endif
430
431// TODO[WB]: We should have a function of the kind
432// AffineConstraints::add_constraint (const size_type constrained_dof,
433// const std::vector<std::pair<size_type, number> > &entries,
434// const number inhomogeneity = 0);
435// rather than building up constraints piecemeal through add_line/add_entry
436// etc. This would also eliminate the possibility of accidentally changing
437// existing constraints into something pointless, see the discussion on the
438// mailing list on "Tiny bug in interpolate_boundary_values" in Sept. 2010.
439
513template <typename number = double>
515{
516public:
521
528 {
534
541
548 };
549
567 explicit AffineConstraints(const IndexSet &local_constraints = IndexSet());
568
572 explicit AffineConstraints(const AffineConstraints &affine_constraints);
573
577 AffineConstraints(AffineConstraints &&affine_constraints) noexcept =
578 default; // NOLINT
579
590 operator=(const AffineConstraints &) = delete;
591
596 operator=(AffineConstraints &&affine_constraints) noexcept =
597 default; // NOLINT
598
605 template <typename other_number>
606 void
608
615 void
616 reinit(const IndexSet &local_constraints = IndexSet());
617
624 bool
625 can_store_line(const size_type line_n) const;
626
633 const IndexSet &
635
655 void
657 const IndexSet & filter);
658
668 void
669 add_line(const size_type line_n);
670
683 void
684 add_lines(const std::vector<bool> &lines);
685
698 void
699 add_lines(const std::set<size_type> &lines);
700
713 void
715
731 void
732 add_entry(const size_type constrained_dof_index,
733 const size_type column,
734 const number weight);
735
741 void
743 const size_type constrained_dof_index,
744 const std::vector<std::pair<size_type, number>> &col_weight_pairs);
745
757 void
758 set_inhomogeneity(const size_type constrained_dof_index, const number value);
759
781 void
783
789 bool
790 is_closed() const;
791
797 bool
798 is_closed(const MPI_Comm comm) const;
799
824 void
826 const AffineConstraints & other_constraints,
827 const MergeConflictBehavior merge_conflict_behavior = no_conflicts_allowed,
828 const bool allow_different_local_lines = false);
829
842 void
843 shift(const size_type offset);
844
852 void
854
869
876
884
894 bool
895 is_constrained(const size_type line_n) const;
896
908 bool
910
917 bool
919 const size_type line_n_2) const;
920
933
938 bool
940
946 bool
948
953 const std::vector<std::pair<size_type, number>> *
955
960 number
961 get_inhomogeneity(const size_type line_n) const;
962
983 void
984 print(std::ostream &out) const;
985
998 void
999 write_dot(std::ostream &) const;
1000
1005 std::size_t
1007
1014 void
1015 resolve_indices(std::vector<types::global_dof_index> &indices) const;
1016
1042 void
1043 condense(SparsityPattern &sparsity) const;
1044
1048 void
1050
1055 void
1057
1062 void
1064
1072 void
1074
1078 void
1080
1092 template <class VectorType>
1093 void
1094 condense(VectorType &vec) const;
1095
1102 template <class VectorType>
1103 void
1104 condense(const VectorType &vec_ghosted, VectorType &output) const;
1105
1118 template <class VectorType>
1119 void
1120 condense(SparseMatrix<number> &matrix, VectorType &vector) const;
1121
1126 template <class BlockVectorType>
1127 void
1128 condense(BlockSparseMatrix<number> &matrix, BlockVectorType &vector) const;
1129
1136 template <class VectorType>
1137 void
1138 set_zero(VectorType &vec) const;
1139
1195 template <class InVector, class OutVector>
1196 void
1197 distribute_local_to_global(const InVector & local_vector,
1198 const std::vector<size_type> &local_dof_indices,
1199 OutVector & global_vector) const;
1200
1248 template <typename VectorType>
1249 void
1251 const std::vector<size_type> &local_dof_indices,
1252 VectorType & global_vector,
1253 const FullMatrix<number> & local_matrix) const;
1254
1274 template <typename VectorType>
1275 void
1277 const Vector<number> & local_vector,
1278 const std::vector<size_type> &local_dof_indices_row,
1279 const std::vector<size_type> &local_dof_indices_col,
1280 VectorType & global_vector,
1281 const FullMatrix<number> & local_matrix,
1282 bool diagonal = false) const;
1283
1287 template <class VectorType>
1288 void
1290 const number value,
1291 VectorType & global_vector) const;
1292
1325 template <typename ForwardIteratorVec,
1326 typename ForwardIteratorInd,
1327 class VectorType>
1328 void
1329 distribute_local_to_global(ForwardIteratorVec local_vector_begin,
1330 ForwardIteratorVec local_vector_end,
1331 ForwardIteratorInd local_indices_begin,
1332 VectorType & global_vector) const;
1333
1385 template <typename MatrixType>
1386 void
1388 const std::vector<size_type> &local_dof_indices,
1389 MatrixType & global_matrix) const;
1390
1418 template <typename MatrixType>
1419 void
1421 const std::vector<size_type> &row_indices,
1422 const std::vector<size_type> &col_indices,
1423 MatrixType & global_matrix) const;
1424
1441 template <typename MatrixType>
1442 void
1444 const std::vector<size_type> &row_indices,
1445 const AffineConstraints &column_affine_constraints,
1446 const std::vector<size_type> &column_indices,
1447 MatrixType & global_matrix) const;
1448
1469 template <typename MatrixType, typename VectorType>
1470 void
1472 const Vector<number> & local_vector,
1473 const std::vector<size_type> &local_dof_indices,
1474 MatrixType & global_matrix,
1475 VectorType & global_vector,
1476 bool use_inhomogeneities_for_rhs = false) const;
1477
1531 void
1533 const std::vector<size_type> &local_dof_indices,
1534 SparsityPatternBase & sparsity_pattern,
1535 const bool keep_constrained_entries = true,
1536 const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1537
1541 void
1543 const std::vector<size_type> &row_indices,
1544 const std::vector<size_type> &col_indices,
1545 SparsityPatternBase & sparsity_pattern,
1546 const bool keep_constrained_entries = true,
1547 const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1548
1553 void
1555 const std::vector<size_type> & row_indices,
1556 const AffineConstraints<number> &col_constraints,
1557 const std::vector<size_type> & col_indices,
1558 SparsityPatternBase & sparsity_pattern,
1559 const bool keep_constrained_entries = true,
1560 const Table<2, bool> & dof_mask = Table<2, bool>()) const;
1561
1581 template <typename ForwardIteratorVec,
1582 typename ForwardIteratorInd,
1583 class VectorType>
1584 void
1585 get_dof_values(const VectorType & global_vector,
1586 ForwardIteratorInd local_indices_begin,
1587 ForwardIteratorVec local_vector_begin,
1588 ForwardIteratorVec local_vector_end) const;
1589
1611 template <class VectorType>
1612 void
1613 distribute(VectorType &vec) const;
1614
1623 {
1628 using Entries = std::vector<std::pair<size_type, number>>;
1629
1636
1645
1650
1655 const Entries & entries = {},
1656 const number & inhomogeneity = 0.0);
1657
1661 template <typename ConstraintLineType>
1662 ConstraintLine(const ConstraintLineType &other);
1663
1667 template <typename ConstraintLineType>
1668 ConstraintLine &
1669 operator=(const ConstraintLineType &other);
1670
1678 bool
1679 operator<(const ConstraintLine &) const;
1680
1686 bool
1688
1693 std::size_t
1695
1701 template <class Archive>
1702 void
1703 serialize(Archive &ar, const unsigned int)
1704 {
1706 }
1707 };
1708
1712 using const_iterator = typename std::vector<ConstraintLine>::const_iterator;
1713
1717 using LineRange = boost::iterator_range<const_iterator>;
1718
1727 const LineRange
1728 get_lines() const;
1729
1758 bool
1759 is_consistent_in_parallel(const std::vector<IndexSet> &locally_owned_dofs,
1760 const IndexSet & locally_active_dofs,
1761 const MPI_Comm mpi_communicator,
1762 const bool verbose = false) const;
1763
1769 void
1770 make_consistent_in_parallel(const IndexSet &locally_owned_dofs,
1771 const IndexSet &locally_relevant_dofs,
1772 const MPI_Comm mpi_communicator);
1773
1792 size_type,
1793 << "The specified line " << arg1 << " does not exist.");
1800 size_type,
1801 size_type,
1802 number,
1803 number,
1804 << "The entry for the indices " << arg1 << " and " << arg2
1805 << " already exists, but the values " << arg3 << " (old) and "
1806 << arg4 << " (new) differ "
1807 << "by " << (arg4 - arg3) << '.');
1814 int,
1815 int,
1816 << "You tried to constrain DoF " << arg1 << " to DoF " << arg2
1817 << ", but that one is also constrained. This is not allowed!");
1824 size_type,
1825 << "Degree of freedom " << arg1
1826 << " is constrained from both object in a merge operation.");
1833 size_type,
1834 << "In the given argument a degree of freedom is constrained "
1835 << "to another DoF with number " << arg1
1836 << ", which however is constrained by this object. This is not"
1837 << " allowed.");
1844 size_type,
1845 << "The index set given to this constraints object indicates "
1846 << "constraints for degree of freedom " << arg1
1847 << " should not be stored by this object, but a constraint "
1848 << "is being added.");
1849
1856 size_type,
1857 size_type,
1858 << "The index set given to this constraints object indicates "
1859 << "constraints using degree of freedom " << arg2
1860 << " should not be stored by this object, but a constraint "
1861 << "for degree of freedom " << arg1 << " uses it.");
1862
1869 int,
1870 int,
1871 << "While distributing the constraint for DoF " << arg1
1872 << ", it turns out that one of the processors "
1873 << "who own the " << arg2 << " degrees of freedom that x_"
1874 << arg1 << " is constrained against does not know about "
1875 << "the constraint on x_" << arg1
1876 << ". Did you not initialize the AffineConstraints container "
1877 << "with the appropriate locally_relevant set so "
1878 << "that every processor who owns a DoF that constrains "
1879 << "another DoF also knows about this constraint?");
1880
1881 template <typename>
1882 friend class AffineConstraints;
1883
1884private:
1896 std::vector<ConstraintLine> lines;
1897
1930 std::vector<size_type> lines_cache;
1931
1938
1943
1945 internal::AffineConstraints::ScratchData<number>>
1947
1952 size_type
1953 calculate_line_index(const size_type line_n) const;
1954
1959 template <typename MatrixType, typename VectorType>
1960 void
1962 const Vector<number> & local_vector,
1963 const std::vector<size_type> &local_dof_indices,
1964 MatrixType & global_matrix,
1965 VectorType & global_vector,
1966 const bool use_inhomogeneities_for_rhs,
1967 const std::integral_constant<bool, false>) const;
1968
1973 template <typename MatrixType, typename VectorType>
1974 void
1976 const Vector<number> & local_vector,
1977 const std::vector<size_type> &local_dof_indices,
1978 MatrixType & global_matrix,
1979 VectorType & global_vector,
1980 const bool use_inhomogeneities_for_rhs,
1981 const std::integral_constant<bool, true>) const;
1982
1990 void
1991 make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
1992 internal::AffineConstraints::GlobalRowsFromLocal<number>
1993 &global_rows) const;
1994
2002 void
2003 make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
2004 std::vector<size_type> & active_dofs) const;
2005
2009 template <typename MatrixScalar, typename VectorScalar>
2012 const size_type i,
2013 const internal::AffineConstraints::GlobalRowsFromLocal<number> &global_rows,
2014 const Vector<VectorScalar> & local_vector,
2015 const std::vector<size_type> & local_dof_indices,
2016 const FullMatrix<MatrixScalar> &local_matrix) const;
2017};
2018
2019/* ---------------- template and inline functions ----------------- */
2020
2021template <typename number>
2023 const IndexSet &local_constraints)
2024 : lines()
2025 , local_lines(local_constraints)
2026 , sorted(false)
2027{
2028 // make sure the IndexSet is compressed. Otherwise this can lead to crashes
2029 // that are hard to find (only happen in release mode).
2030 // see tests/mpi/affine_constraints_crash_01
2032}
2033
2034template <typename number>
2036 const AffineConstraints &affine_constraints)
2037 : Subscriptor()
2038 , lines(affine_constraints.lines)
2039 , lines_cache(affine_constraints.lines_cache)
2040 , local_lines(affine_constraints.local_lines)
2041 , sorted(affine_constraints.sorted)
2042{}
2043
2044template <typename number>
2045inline void
2047{
2048 Assert(sorted == false, ExcMatrixIsClosed());
2049
2050 // the following can happen when we compute with distributed meshes and dof
2051 // handlers and we constrain a degree of freedom whose number we don't have
2052 // locally. if we don't abort here the program will try to allocate several
2053 // terabytes of memory to resize the various arrays below :-)
2055 const size_type line_index = calculate_line_index(line_n);
2056
2057 // check whether line already exists; it may, in which case we can just quit
2058 if (is_constrained(line_n))
2059 return;
2060
2061 // if necessary enlarge vector of existing entries for cache
2062 if (line_index >= lines_cache.size())
2063 lines_cache.resize(std::max(2 * static_cast<size_type>(lines_cache.size()),
2064 line_index + 1),
2066
2067 // push a new line to the end of the list
2068 lines.emplace_back();
2069 lines.back().index = line_n;
2070 lines.back().inhomogeneity = 0.;
2071 lines_cache[line_index] = lines.size() - 1;
2072}
2073
2074
2075
2076template <typename number>
2077inline void
2079 const size_type column,
2080 const number weight)
2081{
2082 Assert(sorted == false, ExcMatrixIsClosed());
2083 Assert(constrained_dof_index != column,
2084 ExcMessage("Can't constrain a degree of freedom to itself"));
2085
2086 // Ensure that the current line is present in the cache:
2087 const size_type line_index = calculate_line_index(constrained_dof_index);
2088 Assert(line_index < lines_cache.size(),
2089 ExcMessage("The current AffineConstraints does not contain the line "
2090 "for the current entry. Call AffineConstraints::add_line "
2091 "before calling this function."));
2092
2093 // if in debug mode, check whether an entry for this column already exists
2094 // and if it's the same as the one entered at present
2095 //
2096 // in any case: exit the function if an entry for this column already
2097 // exists, since we don't want to enter it twice
2101 ExcColumnNotStoredHere(constrained_dof_index, column));
2102 ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
2103 Assert(line_ptr->index == constrained_dof_index, ExcInternalError());
2104 for (const auto &p : line_ptr->entries)
2105 if (p.first == column)
2106 {
2107 Assert(std::abs(p.second - weight) < 1.e-14,
2109 constrained_dof_index, column, p.second, weight));
2110 return;
2111 }
2112
2113 line_ptr->entries.emplace_back(column, weight);
2114}
2115
2116
2117
2118template <typename number>
2119inline void
2121 const size_type constrained_dof_index,
2122 const number value)
2123{
2124 const size_type line_index = calculate_line_index(constrained_dof_index);
2125 Assert(line_index < lines_cache.size() &&
2127 ExcMessage("call add_line() before calling set_inhomogeneity()"));
2128 Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
2129 ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
2130 line_ptr->inhomogeneity = value;
2131}
2132
2133
2134
2135template <typename number>
2136template <class VectorType>
2137inline void
2139{
2140 // since lines is a private member, we cannot pass it to the functions
2141 // above. therefore, copy the content which is cheap
2142 std::vector<size_type> constrained_lines(lines.size());
2143 for (unsigned int i = 0; i < lines.size(); ++i)
2144 constrained_lines[i] = lines[i].index;
2145 internal::AffineConstraintsImplementation::set_zero_all(constrained_lines,
2146 vec);
2147}
2148
2149template <typename number>
2152{
2153 return lines.size();
2154}
2155
2156template <typename number>
2159{
2160 return std::count_if(lines.begin(),
2161 lines.end(),
2162 [](const ConstraintLine &line) {
2163 return (line.entries.size() == 1) &&
2164 (line.entries[0].second == number(1.));
2165 });
2166}
2167
2168template <typename number>
2171{
2172 return std::count_if(lines.begin(),
2173 lines.end(),
2174 [](const ConstraintLine &line) {
2175 return (line.inhomogeneity != number(0.));
2176 });
2177}
2178
2179template <typename number>
2180inline bool
2182{
2183 if (lines.empty())
2184 return false;
2185
2186 const size_type line_index = calculate_line_index(index);
2187 return ((line_index < lines_cache.size()) &&
2188 (lines_cache[line_index] != numbers::invalid_size_type));
2189}
2190
2191template <typename number>
2192inline bool
2194 const size_type line_n) const
2195{
2196 // check whether the entry is constrained. could use is_constrained, but
2197 // that means computing the line index twice
2198 const size_type line_index = calculate_line_index(line_n);
2199 if (line_index >= lines_cache.size() ||
2201 return false;
2202 else
2203 {
2204 Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
2205 return (lines[lines_cache[line_index]].inhomogeneity != number(0.));
2206 }
2207}
2208
2209template <typename number>
2210inline const std::vector<std::pair<types::global_dof_index, number>> *
2212{
2213 if (lines.empty())
2214 return nullptr;
2215
2216 // check whether the entry is constrained. could use is_constrained, but
2217 // that means computing the line index twice
2218 const size_type line_index = calculate_line_index(line_n);
2219 if (line_index >= lines_cache.size() ||
2221 return nullptr;
2222 else
2223 return &lines[lines_cache[line_index]].entries;
2224}
2225
2226template <typename number>
2227inline number
2229{
2230 // check whether the entry is constrained. could use is_constrained, but
2231 // that means computing the line index twice
2232 const size_type line_index = calculate_line_index(line_n);
2233 if (line_index >= lines_cache.size() ||
2235 return 0;
2236 else
2237 return lines[lines_cache[line_index]].inhomogeneity;
2238}
2239
2240template <typename number>
2243{
2244 // IndexSet is unused (serial case)
2245 if (local_lines.size() == 0)
2246 return line_n;
2247
2249
2250 return local_lines.index_within_set(line_n);
2251}
2252
2253template <typename number>
2254inline bool
2256{
2257 return local_lines.size() == 0 || local_lines.is_element(line_n);
2258}
2259
2260template <typename number>
2261inline const IndexSet &
2263{
2264 return local_lines;
2265}
2266
2267template <typename number>
2268template <class VectorType>
2269inline void
2271 const size_type index,
2272 const number value,
2273 VectorType & global_vector) const
2274{
2275 Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2276
2277 if (is_constrained(index) == false)
2278 global_vector(index) += value;
2279 else
2280 {
2281 const ConstraintLine &position =
2283 for (size_type j = 0; j < position.entries.size(); ++j)
2284 global_vector(position.entries[j].first) +=
2285 value * position.entries[j].second;
2286 }
2287}
2288
2289template <typename number>
2290template <typename ForwardIteratorVec,
2291 typename ForwardIteratorInd,
2292 class VectorType>
2293inline void
2295 ForwardIteratorVec local_vector_begin,
2296 ForwardIteratorVec local_vector_end,
2297 ForwardIteratorInd local_indices_begin,
2298 VectorType & global_vector) const
2299{
2300 Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2301 for (; local_vector_begin != local_vector_end;
2302 ++local_vector_begin, ++local_indices_begin)
2303 {
2304 if (is_constrained(*local_indices_begin) == false)
2305 internal::ElementAccess<VectorType>::add(*local_vector_begin,
2306 *local_indices_begin,
2307 global_vector);
2308 else
2309 {
2310 const ConstraintLine &position =
2311 lines[lines_cache[calculate_line_index(*local_indices_begin)]];
2312 for (size_type j = 0; j < position.entries.size(); ++j)
2314 (*local_vector_begin) * position.entries[j].second,
2315 position.entries[j].first,
2316 global_vector);
2317 }
2318 }
2319}
2320
2321template <typename number>
2322template <class InVector, class OutVector>
2323inline void
2325 const InVector & local_vector,
2326 const std::vector<size_type> &local_dof_indices,
2327 OutVector & global_vector) const
2328{
2329 Assert(local_vector.size() == local_dof_indices.size(),
2330 ExcDimensionMismatch(local_vector.size(), local_dof_indices.size()));
2331 distribute_local_to_global(local_vector.begin(),
2332 local_vector.end(),
2333 local_dof_indices.begin(),
2334 global_vector);
2335}
2336
2337template <typename number>
2338template <typename ForwardIteratorVec,
2339 typename ForwardIteratorInd,
2340 class VectorType>
2341inline void
2343 const VectorType & global_vector,
2344 ForwardIteratorInd local_indices_begin,
2345 ForwardIteratorVec local_vector_begin,
2346 ForwardIteratorVec local_vector_end) const
2347{
2348 Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2349 for (; local_vector_begin != local_vector_end;
2350 ++local_vector_begin, ++local_indices_begin)
2351 {
2352 if (is_constrained(*local_indices_begin) == false)
2353 *local_vector_begin = global_vector(*local_indices_begin);
2354 else
2355 {
2356 const ConstraintLine &position =
2357 lines[lines_cache[calculate_line_index(*local_indices_begin)]];
2358 typename VectorType::value_type value = position.inhomogeneity;
2359 for (size_type j = 0; j < position.entries.size(); ++j)
2360 value += (global_vector(position.entries[j].first) *
2361 position.entries[j].second);
2362 *local_vector_begin = value;
2363 }
2364 }
2365}
2366
2367// Forward declarations
2368#ifndef DOXYGEN
2369template <typename MatrixType>
2370class BlockMatrixBase;
2371template <typename SparsityPatternType>
2373template <typename number>
2375
2376namespace internal
2377{
2378 namespace AffineConstraints
2379 {
2397 template <typename MatrixType>
2398 struct IsBlockMatrix
2399 {
2400 private:
2406 template <typename T>
2407 static std::true_type
2408 check(const BlockMatrixBase<T> *);
2409
2414 template <typename T>
2415 static std::true_type
2416 check(const BlockSparseMatrixEZ<T> *);
2417
2422 static std::false_type
2423 check(...);
2424
2425 public:
2432 static const bool value =
2433 std::is_same<decltype(check(std::declval<MatrixType *>())),
2434 std::true_type>::value;
2435 };
2436
2437 // instantiation of the static member
2438 template <typename MatrixType>
2439 const bool IsBlockMatrix<MatrixType>::value;
2440
2441 } // namespace AffineConstraints
2442} // namespace internal
2443#endif
2444
2445
2446
2447template <typename number>
2448template <typename other_number>
2449inline void
2452{
2453 lines.clear();
2454 lines.insert(lines.begin(), other.lines.begin(), other.lines.end());
2455 lines_cache = other.lines_cache;
2456 local_lines = other.local_lines;
2457 sorted = other.sorted;
2458}
2459
2460
2461
2462template <typename number>
2463template <typename MatrixType>
2464inline void
2466 const FullMatrix<number> & local_matrix,
2467 const std::vector<size_type> &local_dof_indices,
2468 MatrixType & global_matrix) const
2469{
2470 // create a dummy and hand on to the function actually implementing this
2471 // feature in the cm.templates.h file.
2473 distribute_local_to_global(
2474 local_matrix,
2475 dummy,
2476 local_dof_indices,
2477 global_matrix,
2478 dummy,
2479 false,
2480 std::integral_constant<
2481 bool,
2482 internal::AffineConstraints::IsBlockMatrix<MatrixType>::value>());
2483}
2484
2485
2486
2487template <typename number>
2488template <typename MatrixType, typename VectorType>
2489inline void
2491 const FullMatrix<number> & local_matrix,
2492 const Vector<number> & local_vector,
2493 const std::vector<size_type> &local_dof_indices,
2494 MatrixType & global_matrix,
2495 VectorType & global_vector,
2496 bool use_inhomogeneities_for_rhs) const
2497{
2498 // enter the internal function with the respective block information set,
2499 // the actual implementation follows in the cm.templates.h file.
2500 distribute_local_to_global(
2501 local_matrix,
2502 local_vector,
2503 local_dof_indices,
2504 global_matrix,
2505 global_vector,
2506 use_inhomogeneities_for_rhs,
2507 std::integral_constant<
2508 bool,
2509 internal::AffineConstraints::IsBlockMatrix<MatrixType>::value>());
2510}
2511
2512
2513
2514template <typename number>
2516 const size_type & index,
2518 const number &inhomogeneity)
2519 : index(index)
2520 , entries(entries)
2521 , inhomogeneity(inhomogeneity)
2522{}
2523
2524
2525
2526template <typename number>
2527template <typename ConstraintLineType>
2529 const ConstraintLineType &other)
2530{
2531 this->index = other.index;
2532
2533 entries.clear();
2534 entries.insert(entries.begin(), other.entries.begin(), other.entries.end());
2535
2536 this->inhomogeneity = other.inhomogeneity;
2537}
2538
2539
2540
2541template <typename number>
2542template <typename ConstraintLineType>
2545 const ConstraintLineType &other)
2546{
2547 this->index = other.index;
2548
2549 entries.clear();
2550 entries.insert(entries.begin(), other.entries.begin(), other.entries.end());
2551
2552 this->inhomogeneity = other.inhomogeneity;
2553
2554 return *this;
2555}
2556
2558
2559#endif
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 add_entries_local_to_global(const std::vector< size_type > &row_indices, const AffineConstraints< number > &col_constraints, const std::vector< size_type > &col_indices, SparsityPatternBase &sparsity_pattern, const bool keep_constrained_entries=true, const Table< 2, bool > &dof_mask=Table< 2, bool >()) 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, 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 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
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
bool is_closed(const MPI_Comm comm) 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
void condense(SparsityPattern &sparsity) const
void distribute_local_to_global(const size_type index, const number value, VectorType &global_vector) const
size_type max_constraint_indirections() 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 add_entries_local_to_global(const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, SparsityPatternBase &sparsity_pattern, const bool keep_constrained_entries=true, const Table< 2, bool > &dof_mask=Table< 2, bool >()) 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
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
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 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:1661
size_type index_within_set(const size_type global_index) const
Definition index_set.h:1883
bool is_element(const size_type index) const
Definition index_set.h:1776
void compress() const
Definition index_set.h:1669
A class that provides a separate storage location on each thread that accesses the object.
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define DeclException0(Exception0)
Definition exceptions.h:465
static ::ExceptionBase & ExcMatrixNotClosed()
static ::ExceptionBase & ExcLineInexistant(size_type arg1)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition exceptions.h:579
static ::ExceptionBase & ExcDoFIsConstrainedFromBothObjects(size_type arg1)
static ::ExceptionBase & ExcDoFConstrainedToConstrainedDoF(int arg1, int arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcMatrixIsClosed()
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:533
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:510
static ::ExceptionBase & ExcEntryAlreadyExists(size_type arg1, size_type arg2, number arg3, number arg4)
static ::ExceptionBase & ExcMessage(std::string arg1)
types::global_dof_index size_type
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
const types::global_dof_index invalid_dof_index
Definition types.h:233
const types::global_dof_index invalid_size_type
Definition types.h:222
::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:82
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