Reference documentation for deal.II version GIT relicensing-1062-gc06da148b8 2024-07-15 19:20:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
affine_constraints.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1998 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_affine_constraints_h
16#define dealii_affine_constraints_h
17
18#include <deal.II/base/config.h>
19
23#include <deal.II/base/table.h>
26
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
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 <typename 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
432
505template <typename number = double>
507{
508public:
513
541
560
583 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT(
584 "Use the constructor with two index set arguments.")
585 explicit AffineConstraints(const IndexSet &locally_stored_constraints);
586
619 const IndexSet &locally_stored_constraints);
620
624 explicit AffineConstraints(const AffineConstraints &affine_constraints);
625
629 AffineConstraints(AffineConstraints &&affine_constraints) noexcept =
630 default; // NOLINT
631
642 operator=(const AffineConstraints &) = delete;
643
648 operator=(AffineConstraints &&affine_constraints) noexcept =
649 default; // NOLINT
650
657 template <typename other_number>
658 void
659 copy_from(const AffineConstraints<other_number> &other);
660
669 void
671
680 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT(
681 "Use the reinit() function with two index set arguments.")
682 void
683 reinit(const IndexSet &locally_stored_constraints);
684
693 void
695 const IndexSet &locally_stored_constraints);
696
705 bool
706 can_store_line(const size_type line_n) const;
707
716 const IndexSet &
718
725 const IndexSet &
727
756 DEAL_II_DEPRECATED_EARLY_WITH_COMMENT("Use get_view() and merge() instead.")
757 void
759 const IndexSet &filter);
760
811 void
813 const size_type constrained_dof,
814 const ArrayView<const std::pair<size_type, number>> &dependencies,
815 const number inhomogeneity = 0);
816
834 void
835 constrain_dof_to_zero(const size_type constrained_dof);
836
841 void
842 add_line(const size_type line_n);
843
856 void
857 add_lines(const std::vector<bool> &lines);
858
871 void
873
886 void
888
904 void
905 add_entry(const size_type constrained_dof_index,
906 const size_type column,
907 const number weight);
908
914 void
916 const size_type constrained_dof_index,
917 const std::vector<std::pair<size_type, number>> &col_weight_pairs);
918
930 void
931 set_inhomogeneity(const size_type constrained_dof_index, const number value);
932
954 void
956
962 bool
963 is_closed() const;
964
970 bool
971 is_closed(const MPI_Comm comm) const;
972
997 template <typename other_number>
998 void
1000 const AffineConstraints<other_number> &other_constraints,
1001 const MergeConflictBehavior merge_conflict_behavior = no_conflicts_allowed,
1002 const bool allow_different_local_lines = false);
1003
1016 void
1017 shift(const size_type offset);
1018
1068 get_view(const IndexSet &mask) const;
1069
1077 void
1079
1092 size_type
1094
1099 size_type
1101
1107 size_type
1109
1119 bool
1120 is_constrained(const size_type line_n) const;
1121
1133 bool
1135
1142 bool
1144 const size_type line_n_2) const;
1145
1156 size_type
1158
1163 bool
1165
1171 bool
1173
1178 const std::vector<std::pair<size_type, number>> *
1180
1185 number
1186 get_inhomogeneity(const size_type line_n) const;
1187
1208 void
1209 print(std::ostream &out) const;
1210
1223 void
1224 write_dot(std::ostream &) const;
1225
1230 std::size_t
1232
1239 void
1240 resolve_indices(std::vector<types::global_dof_index> &indices) const;
1241
1267 void
1268 condense(SparsityPattern &sparsity) const;
1269
1273 void
1275
1280 void
1282
1287 void
1289
1297 void
1298 condense(SparseMatrix<number> &matrix) const;
1299
1303 void
1304 condense(BlockSparseMatrix<number> &matrix) const;
1305
1317 template <typename VectorType>
1318 void
1319 condense(VectorType &vec) const;
1320
1327 template <typename VectorType>
1328 void
1329 condense(const VectorType &vec_ghosted, VectorType &output) const;
1330
1343 template <typename VectorType>
1344 void
1345 condense(SparseMatrix<number> &matrix, VectorType &vector) const;
1346
1351 template <typename BlockVectorType>
1352 void
1353 condense(BlockSparseMatrix<number> &matrix, BlockVectorType &vector) const;
1354
1361 template <typename VectorType>
1362 void
1363 set_zero(VectorType &vec) const;
1364
1420 template <class InVector, class OutVector>
1421 void
1422 distribute_local_to_global(const InVector &local_vector,
1423 const std::vector<size_type> &local_dof_indices,
1424 OutVector &global_vector) const;
1425
1473 template <typename VectorType>
1474 void
1475 distribute_local_to_global(const Vector<number> &local_vector,
1476 const std::vector<size_type> &local_dof_indices,
1477 VectorType &global_vector,
1478 const FullMatrix<number> &local_matrix) const;
1479
1499 template <typename VectorType>
1500 void
1502 const Vector<number> &local_vector,
1503 const std::vector<size_type> &local_dof_indices_row,
1504 const std::vector<size_type> &local_dof_indices_col,
1505 VectorType &global_vector,
1506 const FullMatrix<number> &local_matrix,
1507 bool diagonal = false) const;
1508
1512 template <typename VectorType>
1513 void
1515 const number value,
1516 VectorType &global_vector) const;
1517
1550 template <typename ForwardIteratorVec,
1551 typename ForwardIteratorInd,
1552 typename VectorType>
1553 void
1554 distribute_local_to_global(ForwardIteratorVec local_vector_begin,
1555 ForwardIteratorVec local_vector_end,
1556 ForwardIteratorInd local_indices_begin,
1557 VectorType &global_vector) const;
1558
1610 template <typename MatrixType>
1611 void
1612 distribute_local_to_global(const FullMatrix<number> &local_matrix,
1613 const std::vector<size_type> &local_dof_indices,
1614 MatrixType &global_matrix) const;
1615
1643 template <typename MatrixType>
1644 void
1645 distribute_local_to_global(const FullMatrix<number> &local_matrix,
1646 const std::vector<size_type> &row_indices,
1647 const std::vector<size_type> &col_indices,
1648 MatrixType &global_matrix) const;
1649
1666 template <typename MatrixType>
1667 void
1668 distribute_local_to_global(const FullMatrix<number> &local_matrix,
1669 const std::vector<size_type> &row_indices,
1670 const AffineConstraints &column_affine_constraints,
1671 const std::vector<size_type> &column_indices,
1672 MatrixType &global_matrix) const;
1673
1694 template <typename MatrixType, typename VectorType>
1695 void
1696 distribute_local_to_global(const FullMatrix<number> &local_matrix,
1697 const Vector<number> &local_vector,
1698 const std::vector<size_type> &local_dof_indices,
1699 MatrixType &global_matrix,
1700 VectorType &global_vector,
1701 bool use_inhomogeneities_for_rhs = false) const;
1702
1756 void
1758 const std::vector<size_type> &local_dof_indices,
1759 SparsityPatternBase &sparsity_pattern,
1760 const bool keep_constrained_entries = true,
1761 const Table<2, bool> &dof_mask = Table<2, bool>()) const;
1762
1766 void
1768 const std::vector<size_type> &row_indices,
1769 const std::vector<size_type> &col_indices,
1770 SparsityPatternBase &sparsity_pattern,
1771 const bool keep_constrained_entries = true,
1772 const Table<2, bool> &dof_mask = Table<2, bool>()) const;
1773
1778 void
1780 const std::vector<size_type> &row_indices,
1781 const AffineConstraints<number> &col_constraints,
1782 const std::vector<size_type> &col_indices,
1783 SparsityPatternBase &sparsity_pattern,
1784 const bool keep_constrained_entries = true,
1785 const Table<2, bool> &dof_mask = Table<2, bool>()) const;
1786
1806 template <typename ForwardIteratorVec,
1807 typename ForwardIteratorInd,
1808 typename VectorType>
1809 void
1810 get_dof_values(const VectorType &global_vector,
1811 ForwardIteratorInd local_indices_begin,
1812 ForwardIteratorVec local_vector_begin,
1813 ForwardIteratorVec local_vector_end) const;
1814
1836 template <typename VectorType>
1837 void
1838 distribute(VectorType &vec) const;
1839
1848 {
1853 using Entries = std::vector<std::pair<size_type, number>>;
1854
1861
1870
1875
1879 ConstraintLine(const size_type &index = numbers::invalid_dof_index,
1880 const typename AffineConstraints<
1881 number>::ConstraintLine::Entries &entries = {},
1882 const number inhomogeneity = 0.0);
1883
1887 ConstraintLine(const ConstraintLine &other) = default;
1888
1892 ConstraintLine(ConstraintLine &&other) noexcept = default;
1893
1898 operator=(const ConstraintLine &other) = default;
1899
1904 operator=(ConstraintLine &&other) noexcept = default;
1905
1910 std::size_t
1912
1918 template <class Archive>
1919 void
1920 serialize(Archive &ar, const unsigned int)
1921 {
1922 ar &index &entries &inhomogeneity;
1923 }
1924
1928 friend void
1930 {
1931 std::swap(l1.index, l2.index);
1932 std::swap(l1.entries, l2.entries);
1933 std::swap(l1.inhomogeneity, l2.inhomogeneity);
1934 }
1935 };
1936
1940 using const_iterator = typename std::vector<ConstraintLine>::const_iterator;
1941
1945 using LineRange = boost::iterator_range<const_iterator>;
1946
1955 LineRange
1956 get_lines() const;
1957
1986 bool
1987 is_consistent_in_parallel(const std::vector<IndexSet> &locally_owned_dofs,
1988 const IndexSet &locally_active_dofs,
1989 const MPI_Comm mpi_communicator,
1990 const bool verbose = false) const;
1991
2022 void
2024 const IndexSet &constraints_to_make_consistent,
2025 const MPI_Comm mpi_communicator);
2026
2034 "You are attempting an operation on an AffineConstraints object "
2035 "that requires the object to not be 'closed', i.e., for which you "
2036 "must not already have called the close() member function. But the "
2037 "object is already closed, and so the operation can not be "
2038 "performed.");
2046 "You are attempting an operation on an AffineConstraints object "
2047 "that requires the object to be 'closed', i.e., for which you "
2048 "needed to call the close() member function. But the object "
2049 "is not currently closed, and so the operation can not be "
2050 "performed.");
2057 size_type,
2058 << "The specified line " << arg1 << " does not exist.");
2065 size_type,
2066 size_type,
2067 number,
2068 number,
2069 << "The entry for the indices " << arg1 << " and " << arg2
2070 << " already exists, but the values " << arg3 << " (old) and "
2071 << arg4 << " (new) differ "
2072 << "by " << (arg4 - arg3) << '.');
2079 int,
2080 int,
2081 << "You tried to constrain DoF " << arg1 << " to DoF " << arg2
2082 << ", but that one is also constrained. This is not allowed!");
2089 size_type,
2090 << "Degree of freedom " << arg1
2091 << " is constrained from both object in a merge operation.");
2098 size_type,
2099 << "In the given argument a degree of freedom is constrained "
2100 << "to another DoF with number " << arg1
2101 << ", which however is constrained by this object. This is not"
2102 << " allowed.");
2109 size_type,
2110 << "The index set given to this constraints object indicates "
2111 << "constraints for degree of freedom " << arg1
2112 << " should not be stored by this object, but a constraint "
2113 << "is being added.");
2114
2121 size_type,
2122 size_type,
2123 << "The index set given to this constraints object indicates "
2124 << "constraints using degree of freedom " << arg2
2125 << " should not be stored by this object, but a constraint "
2126 << "for degree of freedom " << arg1 << " uses it.");
2127
2134 int,
2135 int,
2136 << "While distributing the constraint for DoF " << arg1
2137 << ", it turns out that one of the processors "
2138 << "who own the " << arg2 << " degrees of freedom that x_"
2139 << arg1 << " is constrained against does not know about "
2140 << "the constraint on x_" << arg1
2141 << ". Did you not initialize the AffineConstraints container "
2142 << "with the appropriate locally_relevant set so "
2143 << "that every processor who owns a DoF that constrains "
2144 << "another DoF also knows about this constraint?");
2145
2146 template <typename>
2147 friend class AffineConstraints;
2148
2149private:
2161 std::vector<ConstraintLine> lines;
2162
2195 std::vector<size_type> lines_cache;
2196
2204
2211
2221
2226
2228 internal::AffineConstraints::ScratchData<number>>
2230
2235 size_type
2236 calculate_line_index(const size_type line_n) const;
2237
2242 template <typename MatrixType, typename VectorType>
2243 void
2245 const Vector<number> &local_vector,
2246 const std::vector<size_type> &local_dof_indices,
2247 MatrixType &global_matrix,
2248 VectorType &global_vector,
2249 const bool use_inhomogeneities_for_rhs,
2250 const std::bool_constant<false>) const;
2251
2256 template <typename MatrixType, typename VectorType>
2257 void
2259 const Vector<number> &local_vector,
2260 const std::vector<size_type> &local_dof_indices,
2261 MatrixType &global_matrix,
2262 VectorType &global_vector,
2263 const bool use_inhomogeneities_for_rhs,
2264 const std::bool_constant<true>) const;
2265
2273 void
2274 make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
2275 internal::AffineConstraints::GlobalRowsFromLocal<number>
2276 &global_rows) const;
2277
2285 void
2286 make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
2287 std::vector<size_type> &active_dofs) const;
2288
2292 template <typename MatrixScalar, typename VectorScalar>
2295 const size_type i,
2296 const internal::AffineConstraints::GlobalRowsFromLocal<number> &global_rows,
2297 const Vector<VectorScalar> &local_vector,
2298 const std::vector<size_type> &local_dof_indices,
2299 const FullMatrix<MatrixScalar> &local_matrix) const;
2300};
2301
2302/* ---------------- template and inline functions ----------------- */
2303
2304template <typename number>
2308
2309
2310
2311template <typename number>
2313 const IndexSet &locally_stored_constraints)
2314 : AffineConstraints<number>(locally_stored_constraints,
2315 locally_stored_constraints)
2316{}
2317
2318
2319
2320template <typename number>
2323 const IndexSet &locally_stored_constraints)
2324 : lines()
2326 , local_lines(locally_stored_constraints)
2327 , sorted(false)
2328{
2329 Assert(locally_owned_dofs.is_subset_of(locally_stored_constraints),
2330 ExcMessage("The set of locally stored constraints needs to be a "
2331 "superset of the locally owned DoFs."));
2332
2333 // make sure the IndexSet is compressed. Otherwise this can lead to crashes
2334 // that are hard to find (only happen in release mode).
2335 // see tests/mpi/affine_constraints_crash_01
2337}
2338
2339
2340
2341template <typename number>
2343 const AffineConstraints &affine_constraints)
2344 : Subscriptor()
2345 , lines(affine_constraints.lines)
2346 , lines_cache(affine_constraints.lines_cache)
2347 , locally_owned_dofs(affine_constraints.locally_owned_dofs)
2348 , local_lines(affine_constraints.local_lines)
2350 affine_constraints.needed_elements_for_distribute)
2351 , sorted(affine_constraints.sorted)
2352{}
2353
2354
2355
2356template <typename number>
2357inline void
2359{
2360 Assert(sorted == false, ExcMatrixIsClosed());
2361
2362 // the following can happen when we compute with distributed meshes and dof
2363 // handlers and we constrain a degree of freedom whose number we don't have
2364 // locally. if we don't abort here the program will try to allocate several
2365 // terabytes of memory to resize the various arrays below :-)
2367 const size_type line_index = calculate_line_index(line_n);
2368
2369 // check whether line already exists; it may, in which case we can just quit
2370 if (is_constrained(line_n))
2371 return;
2372
2373 // if necessary enlarge vector of existing entries for cache
2374 if (line_index >= lines_cache.size())
2375 lines_cache.resize(std::max(2 * static_cast<size_type>(lines_cache.size()),
2376 line_index + 1),
2378
2379 // push a new line to the end of the list
2380 lines.emplace_back();
2381 lines.back().index = line_n;
2382 lines.back().inhomogeneity = 0.;
2383 lines_cache[line_index] = lines.size() - 1;
2384}
2385
2386
2387
2388template <typename number>
2389inline void
2391 const size_type column,
2392 const number weight)
2393{
2394 Assert(sorted == false, ExcMatrixIsClosed());
2395 Assert(constrained_dof_index != column,
2396 ExcMessage("Can't constrain a degree of freedom to itself"));
2397
2398 // Ensure that the current line is present in the cache:
2399 const size_type line_index = calculate_line_index(constrained_dof_index);
2400 Assert(line_index < lines_cache.size(),
2401 ExcMessage("The current AffineConstraints does not contain the line "
2402 "for the current entry. Call AffineConstraints::add_line "
2403 "before calling this function."));
2404
2405 // if in debug mode, check whether an entry for this column already exists
2406 // and if it's the same as the one entered at present
2407 //
2408 // in any case: exit the function if an entry for this column already
2409 // exists, since we don't want to enter it twice
2413 ExcColumnNotStoredHere(constrained_dof_index, column));
2414 ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
2415 Assert(line_ptr->index == constrained_dof_index, ExcInternalError());
2416 for (const auto &p : line_ptr->entries)
2417 if (p.first == column)
2418 {
2419 Assert(std::abs(p.second - weight) < 1.e-14,
2421 constrained_dof_index, column, p.second, weight));
2422 return;
2423 }
2424
2425 line_ptr->entries.emplace_back(column, weight);
2426}
2427
2428
2429
2430template <typename number>
2431inline void
2433 const size_type constrained_dof_index,
2434 const number value)
2435{
2436 const size_type line_index = calculate_line_index(constrained_dof_index);
2437 Assert(line_index < lines_cache.size() &&
2439 ExcMessage("call add_line() before calling set_inhomogeneity()"));
2440 Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
2441 ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
2442 line_ptr->inhomogeneity = value;
2443}
2444
2445
2446
2447template <typename number>
2448template <typename VectorType>
2449inline void
2451{
2452 // since lines is a private member, we cannot pass it to the functions
2453 // above. therefore, copy the content which is cheap
2454 std::vector<size_type> constrained_lines(lines.size());
2455 for (unsigned int i = 0; i < lines.size(); ++i)
2456 constrained_lines[i] = lines[i].index;
2457 internal::AffineConstraintsImplementation::set_zero_all(constrained_lines,
2458 vec);
2459}
2460
2461template <typename number>
2464{
2465 return lines.size();
2466}
2467
2468template <typename number>
2471{
2472 return std::count_if(lines.begin(),
2473 lines.end(),
2474 [](const ConstraintLine &line) {
2475 return (line.entries.size() == 1) &&
2476 (line.entries[0].second == number(1.));
2477 });
2478}
2479
2480template <typename number>
2483{
2484 return std::count_if(lines.begin(),
2485 lines.end(),
2486 [](const ConstraintLine &line) {
2487 return (line.inhomogeneity != number(0.));
2488 });
2489}
2490
2491template <typename number>
2492inline bool
2494{
2495 if (lines.empty())
2496 return false;
2497
2498 const size_type line_index = calculate_line_index(index);
2499 return ((line_index < lines_cache.size()) &&
2500 (lines_cache[line_index] != numbers::invalid_size_type));
2501}
2502
2503template <typename number>
2504inline bool
2506 const size_type line_n) const
2507{
2508 // check whether the entry is constrained. could use is_constrained, but
2509 // that means computing the line index twice
2510 const size_type line_index = calculate_line_index(line_n);
2511 if (line_index >= lines_cache.size() ||
2513 return false;
2514 else
2515 {
2516 Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
2517 return (lines[lines_cache[line_index]].inhomogeneity != number(0.));
2518 }
2519}
2520
2521template <typename number>
2522inline const std::vector<std::pair<types::global_dof_index, number>> *
2524{
2525 if (lines.empty())
2526 return nullptr;
2527
2528 // check whether the entry is constrained. could use is_constrained, but
2529 // that means computing the line index twice
2530 const size_type line_index = calculate_line_index(line_n);
2531 if (line_index >= lines_cache.size() ||
2533 return nullptr;
2534 else
2535 return &lines[lines_cache[line_index]].entries;
2536}
2537
2538template <typename number>
2539inline number
2541{
2542 // check whether the entry is constrained. could use is_constrained, but
2543 // that means computing the line index twice
2544 const size_type line_index = calculate_line_index(line_n);
2545 if (line_index >= lines_cache.size() ||
2547 return 0;
2548 else
2549 return lines[lines_cache[line_index]].inhomogeneity;
2550}
2551
2552template <typename number>
2555{
2556 // IndexSet is unused (serial case)
2557 if (local_lines.size() == 0)
2558 return line_n;
2559
2561
2562 return local_lines.index_within_set(line_n);
2563}
2564
2565
2566
2567template <typename number>
2568inline bool
2570{
2571 return local_lines.size() == 0 || local_lines.is_element(line_n);
2572}
2573
2574
2575
2576template <typename number>
2577inline const IndexSet &
2582
2583
2584
2585template <typename number>
2586inline const IndexSet &
2591
2592
2593
2594template <typename number>
2595template <typename VectorType>
2596inline void
2598 const size_type index,
2599 const number value,
2600 VectorType &global_vector) const
2601{
2602 Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2603
2604 if (is_constrained(index) == false)
2605 global_vector(index) += value;
2606 else
2607 {
2608 const ConstraintLine &position =
2610 for (size_type j = 0; j < position.entries.size(); ++j)
2611 global_vector(position.entries[j].first) +=
2612 value * position.entries[j].second;
2613 }
2614}
2615
2616template <typename number>
2617template <typename ForwardIteratorVec,
2618 typename ForwardIteratorInd,
2619 typename VectorType>
2620inline void
2622 ForwardIteratorVec local_vector_begin,
2623 ForwardIteratorVec local_vector_end,
2624 ForwardIteratorInd local_indices_begin,
2625 VectorType &global_vector) const
2626{
2627 Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2628 for (; local_vector_begin != local_vector_end;
2629 ++local_vector_begin, ++local_indices_begin)
2630 {
2631 if (is_constrained(*local_indices_begin) == false)
2632 internal::ElementAccess<VectorType>::add(*local_vector_begin,
2633 *local_indices_begin,
2634 global_vector);
2635 else
2636 {
2637 const ConstraintLine &position =
2638 lines[lines_cache[calculate_line_index(*local_indices_begin)]];
2639 for (size_type j = 0; j < position.entries.size(); ++j)
2641 (*local_vector_begin) * position.entries[j].second,
2642 position.entries[j].first,
2643 global_vector);
2644 }
2645 }
2646}
2647
2648template <typename number>
2649template <class InVector, class OutVector>
2650inline void
2652 const InVector &local_vector,
2653 const std::vector<size_type> &local_dof_indices,
2654 OutVector &global_vector) const
2655{
2656 Assert(global_vector.has_ghost_elements() == false, ExcGhostsPresent());
2657 Assert(local_vector.size() == local_dof_indices.size(),
2658 ExcDimensionMismatch(local_vector.size(), local_dof_indices.size()));
2659 distribute_local_to_global(local_vector.begin(),
2660 local_vector.end(),
2661 local_dof_indices.begin(),
2662 global_vector);
2663}
2664
2665template <typename number>
2666template <typename ForwardIteratorVec,
2667 typename ForwardIteratorInd,
2668 typename VectorType>
2669inline void
2671 const VectorType &global_vector,
2672 ForwardIteratorInd local_indices_begin,
2673 ForwardIteratorVec local_vector_begin,
2674 ForwardIteratorVec local_vector_end) const
2675{
2676 Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
2677 for (; local_vector_begin != local_vector_end;
2678 ++local_vector_begin, ++local_indices_begin)
2679 {
2680 if (is_constrained(*local_indices_begin) == false)
2681 *local_vector_begin = global_vector(*local_indices_begin);
2682 else
2683 {
2684 const ConstraintLine &position =
2685 lines[lines_cache[calculate_line_index(*local_indices_begin)]];
2686 typename VectorType::value_type value = position.inhomogeneity;
2687 for (size_type j = 0; j < position.entries.size(); ++j)
2688 value += (global_vector(position.entries[j].first) *
2689 position.entries[j].second);
2690 *local_vector_begin = value;
2691 }
2692 }
2693}
2694
2695// Forward declarations
2696#ifndef DOXYGEN
2697template <typename MatrixType>
2698class BlockMatrixBase;
2699template <typename SparsityPatternType>
2701template <typename number>
2703
2704namespace internal
2705{
2706 namespace AffineConstraints
2707 {
2725 template <typename MatrixType>
2726 struct IsBlockMatrix
2727 {
2728 private:
2734 template <typename T>
2735 static std::true_type
2736 check(const BlockMatrixBase<T> *);
2737
2742 template <typename T>
2743 static std::true_type
2744 check(const BlockSparseMatrixEZ<T> *);
2745
2750 static std::false_type
2751 check(...);
2752
2753 public:
2760 static const bool value =
2761 std::is_same_v<decltype(check(std::declval<MatrixType *>())),
2762 std::true_type>;
2763 };
2764
2765 // instantiation of the static member
2766 template <typename MatrixType>
2767 const bool IsBlockMatrix<MatrixType>::value;
2768
2769 } // namespace AffineConstraints
2770} // namespace internal
2771#endif
2772
2773
2774
2775template <typename number>
2776template <typename other_number>
2777inline void
2780{
2781 lines.clear();
2782 lines.reserve(other.lines.size());
2783
2784 for (const auto &l : other.lines)
2785 lines.emplace_back(l.index,
2786 typename ConstraintLine::Entries(l.entries.begin(),
2787 l.entries.end()),
2788 l.inhomogeneity);
2789
2790 lines_cache = other.lines_cache;
2791 local_lines = other.local_lines;
2792 sorted = other.sorted;
2793
2794 locally_owned_dofs = other.locally_owned_dofs;
2795 needed_elements_for_distribute = other.needed_elements_for_distribute;
2796}
2797
2798
2799template <typename number>
2800template <typename other_number>
2801void
2803 const AffineConstraints<other_number> &other_constraints,
2804 const MergeConflictBehavior merge_conflict_behavior,
2805 const bool allow_different_local_lines)
2806{
2807 (void)allow_different_local_lines;
2808 Assert(allow_different_local_lines ||
2809 local_lines == other_constraints.local_lines,
2810 ExcMessage(
2811 "local_lines for this and the other objects are not the same "
2812 "although allow_different_local_lines is false."));
2813
2814 // store the previous state with respect to sorting
2815 const bool object_was_sorted = sorted;
2816 sorted = false;
2817
2818 // first action is to fold into the present object possible constraints
2819 // in the second object. we don't strictly need to do this any more since
2820 // the AffineConstraints container has learned to deal with chains of
2821 // constraints in the close() function, but we have traditionally done
2822 // this and it's not overly hard to do.
2823 //
2824 // for this, loop over all constraints and replace the constraint lines
2825 // with a new one where constraints are replaced if necessary.
2826 typename ConstraintLine::Entries tmp;
2827 for (ConstraintLine &line : lines)
2828 {
2829 tmp.clear();
2830 for (const std::pair<size_type, number> &entry : line.entries)
2831 {
2832 // if the present dof is not stored, or not constrained, or if we
2833 // won't take the constraint from the other object, then simply copy
2834 // it over
2835 if ((other_constraints.local_lines.size() != 0. &&
2836 other_constraints.local_lines.is_element(entry.first) ==
2837 false) ||
2838 other_constraints.is_constrained(entry.first) == false ||
2839 ((merge_conflict_behavior != right_object_wins) &&
2840 other_constraints.is_constrained(entry.first) &&
2841 this->is_constrained(entry.first)))
2842 tmp.push_back(entry);
2843 else
2844 // otherwise resolve further constraints by replacing the old
2845 // entry by a sequence of new entries taken from the other
2846 // object, but with multiplied weights
2847 {
2848 const auto *other_entries =
2849 other_constraints.get_constraint_entries(entry.first);
2850 Assert(other_entries != nullptr, ExcInternalError());
2851
2852 const number weight = entry.second;
2853
2854 for (const auto &other_entry : *other_entries)
2855 tmp.emplace_back(other_entry.first,
2856 other_entry.second * weight);
2857
2858 line.inhomogeneity +=
2859 other_constraints.get_inhomogeneity(entry.first) * weight;
2860 }
2861 }
2862 // finally exchange old and newly resolved line
2863 line.entries.swap(tmp);
2864 }
2865
2866 if (local_lines.size() != 0)
2867 local_lines.add_indices(other_constraints.local_lines);
2868
2869 {
2870 // do not bother to resize the lines cache exactly since it is pretty
2871 // cheap to adjust it along the way.
2872 std::fill(lines_cache.begin(),
2873 lines_cache.end(),
2875
2876 // reset lines_cache for our own constraints
2877 size_type index = 0;
2878 for (const ConstraintLine &line : lines)
2879 {
2880 const size_type local_line_no = calculate_line_index(line.index);
2881 if (local_line_no >= lines_cache.size())
2882 lines_cache.resize(local_line_no + 1, numbers::invalid_size_type);
2883 lines_cache[local_line_no] = index++;
2884 }
2885
2886 // Add other_constraints to lines cache and our list of constraints
2887 for (const auto &line : other_constraints.lines)
2888 {
2889 const size_type local_line_no = calculate_line_index(line.index);
2890 if (local_line_no >= lines_cache.size())
2891 {
2892 lines_cache.resize(local_line_no + 1, numbers::invalid_size_type);
2893 lines.emplace_back(line.index,
2894 typename ConstraintLine::Entries(
2895 line.entries.begin(), line.entries.end()),
2896 line.inhomogeneity);
2897 lines_cache[local_line_no] = index++;
2898 }
2899 else if (lines_cache[local_line_no] == numbers::invalid_size_type)
2900 {
2901 // there are no constraints for that line yet
2902 lines.emplace_back(line.index,
2903 typename ConstraintLine::Entries(
2904 line.entries.begin(), line.entries.end()),
2905 line.inhomogeneity);
2906 AssertIndexRange(local_line_no, lines_cache.size());
2907 lines_cache[local_line_no] = index++;
2908 }
2909 else
2910 {
2911 // we already store that line
2912 switch (merge_conflict_behavior)
2913 {
2914 case no_conflicts_allowed:
2915 AssertThrow(false,
2916 ExcDoFIsConstrainedFromBothObjects(line.index));
2917 break;
2918
2919 case left_object_wins:
2920 // ignore this constraint
2921 break;
2922
2923 case right_object_wins:
2924 AssertIndexRange(local_line_no, lines_cache.size());
2925 lines[lines_cache[local_line_no]] = {
2926 line.index,
2927 typename ConstraintLine::Entries(line.entries.begin(),
2928 line.entries.end()),
2929 static_cast<number>(line.inhomogeneity)};
2930 break;
2931
2932 default:
2934 }
2935 }
2936 }
2937
2938 // check that we set the pointers correctly
2939 for (size_type i = 0; i < lines_cache.size(); ++i)
2940 if (lines_cache[i] != numbers::invalid_size_type)
2941 Assert(i == calculate_line_index(lines[lines_cache[i]].index),
2943 }
2944
2945 // if the object was sorted before, then make sure it is so afterward as
2946 // well. otherwise leave everything in the unsorted state
2947 if (object_was_sorted == true)
2948 close();
2949}
2950
2951
2952
2953template <typename number>
2954template <typename MatrixType>
2955inline void
2957 const FullMatrix<number> &local_matrix,
2958 const std::vector<size_type> &local_dof_indices,
2959 MatrixType &global_matrix) const
2960{
2961 // create a dummy and hand on to the function actually implementing this
2962 // feature in the cm.templates.h file.
2964 distribute_local_to_global(
2965 local_matrix,
2966 dummy,
2967 local_dof_indices,
2968 global_matrix,
2969 dummy,
2970 false,
2971 std::integral_constant<
2972 bool,
2973 internal::AffineConstraints::IsBlockMatrix<MatrixType>::value>());
2974}
2975
2976
2977
2978template <typename number>
2979template <typename MatrixType, typename VectorType>
2980inline void
2982 const FullMatrix<number> &local_matrix,
2983 const Vector<number> &local_vector,
2984 const std::vector<size_type> &local_dof_indices,
2985 MatrixType &global_matrix,
2986 VectorType &global_vector,
2987 bool use_inhomogeneities_for_rhs) const
2988{
2989 // enter the internal function with the respective block information set,
2990 // the actual implementation follows in the cm.templates.h file.
2991 distribute_local_to_global(
2992 local_matrix,
2993 local_vector,
2994 local_dof_indices,
2995 global_matrix,
2996 global_vector,
2997 use_inhomogeneities_for_rhs,
2998 std::integral_constant<
2999 bool,
3000 internal::AffineConstraints::IsBlockMatrix<MatrixType>::value>());
3001}
3002
3003
3004
3005template <typename number>
3007 const size_type &index,
3009 const number inhomogeneity)
3010 : index(index)
3011 , entries(entries)
3012 , inhomogeneity(inhomogeneity)
3013{}
3014
3015
3016
3018
3019#endif
void add_entries(const size_type constrained_dof_index, const std::vector< std::pair< size_type, number > > &col_weight_pairs)
bool is_closed() const
size_type n_inhomogeneities() const
bool has_inhomogeneities() const
LineRange get_lines() const
void add_line(const size_type line_n)
const IndexSet & get_locally_owned_indices() const
void add_selected_constraints(const AffineConstraints &constraints_in, const IndexSet &filter)
bool can_store_line(const size_type line_n) const
void add_constraint(const size_type constrained_dof, const ArrayView< const std::pair< size_type, number > > &dependencies, const number inhomogeneity=0)
Threads::ThreadLocalStorage< internal::AffineConstraints::ScratchData< number > > scratch_data
friend class AffineConstraints
void merge(const AffineConstraints< other_number > &other_constraints, const MergeConflictBehavior merge_conflict_behavior=no_conflicts_allowed, const bool allow_different_local_lines=false)
void add_entry(const size_type constrained_dof_index, const size_type column, const number weight)
number get_inhomogeneity(const size_type line_n) const
void make_sorted_row_list(const std::vector< size_type > &local_dof_indices, std::vector< size_type > &active_dofs) const
void get_dof_values(const VectorType &global_vector, ForwardIteratorInd local_indices_begin, ForwardIteratorVec local_vector_begin, ForwardIteratorVec local_vector_end) const
void distribute_local_to_global(const InVector &local_vector, const std::vector< size_type > &local_dof_indices, OutVector &global_vector) const
void make_sorted_row_list(const std::vector< size_type > &local_dof_indices, internal::AffineConstraints::GlobalRowsFromLocal< number > &global_rows) const
const IndexSet & get_local_lines() const
void add_entries_local_to_global(const std::vector< size_type > &local_dof_indices, SparsityPatternBase &sparsity_pattern, const bool keep_constrained_entries=true, const Table< 2, bool > &dof_mask=Table< 2, bool >()) const
void shift(const size_type offset)
size_type n_identities() const
std::size_t memory_consumption() const
typename std::vector< ConstraintLine >::const_iterator const_iterator
void set_inhomogeneity(const size_type constrained_dof_index, const number value)
void condense(SparsityPattern &sparsity) const
size_type max_constraint_indirections() const
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 distribute(VectorType &vec) const
void resolve_indices(std::vector< types::global_dof_index > &indices) const
std::vector< ConstraintLine > lines
bool is_constrained(const size_type line_n) const
AffineConstraints get_view(const IndexSet &mask) const
IndexSet needed_elements_for_distribute
void make_consistent_in_parallel(const IndexSet &locally_owned_dofs, const IndexSet &constraints_to_make_consistent, const MPI_Comm mpi_communicator)
const std::vector< std::pair< size_type, number > > * get_constraint_entries(const size_type line_n) const
size_type calculate_line_index(const size_type line_n) const
void distribute_local_to_global(const FullMatrix< number > &local_matrix, const Vector< number > &local_vector, const std::vector< size_type > &local_dof_indices, MatrixType &global_matrix, VectorType &global_vector, const bool use_inhomogeneities_for_rhs, const std::bool_constant< true >) const
bool is_inhomogeneously_constrained(const size_type index) const
void copy_from(const AffineConstraints< other_number > &other)
std::vector< size_type > lines_cache
void constrain_dof_to_zero(const size_type constrained_dof)
void distribute_local_to_global(const FullMatrix< number > &local_matrix, const Vector< number > &local_vector, const std::vector< size_type > &local_dof_indices, MatrixType &global_matrix, VectorType &global_vector, const bool use_inhomogeneities_for_rhs, const std::bool_constant< false >) const
bool is_consistent_in_parallel(const std::vector< IndexSet > &locally_owned_dofs, const IndexSet &locally_active_dofs, const MPI_Comm mpi_communicator, const bool verbose=false) const
bool is_identity_constrained(const size_type line_n) const
size_type n_constraints() const
void print(std::ostream &out) const
void write_dot(std::ostream &) const
void set_zero(VectorType &vec) const
void add_lines(const std::vector< bool > &lines)
boost::iterator_range< const_iterator > LineRange
bool are_identity_constrained(const size_type line_n_1, const size_type line_n_2) const
bool is_subset_of(const IndexSet &other) const
Definition index_set.cc:704
size_type size() const
Definition index_set.h:1766
size_type index_within_set(const size_type global_index) const
Definition index_set.h:1991
bool is_element(const size_type index) const
Definition index_set.h:1884
void compress() const
Definition index_set.h:1774
A class that provides a separate storage location on each thread that accesses the object.
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcMatrixNotClosed()
static ::ExceptionBase & ExcGhostsPresent()
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition exceptions.h:585
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:539
static ::ExceptionBase & ExcDoFIsConstrainedToConstrainedDoF(size_type arg1)
static ::ExceptionBase & ExcColumnNotStoredHere(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcRowNotStoredHere(size_type arg1)
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:494
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:516
static ::ExceptionBase & ExcEntryAlreadyExists(size_type arg1, size_type arg2, number arg3, number arg4)
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcLineInexistent(size_type arg1)
#define AssertThrow(cond, exc)
#define DEAL_II_NOT_IMPLEMENTED()
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:252
const types::global_dof_index invalid_size_type
Definition types.h:233
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition types.h:32
unsigned int global_dof_index
Definition types.h:81
*braid_SplitCommworld & comm
void serialize(Archive &ar, const unsigned int)
ConstraintLine & operator=(const ConstraintLine &other)=default
ConstraintLine(const ConstraintLine &other)=default
ConstraintLine(const size_type &index=numbers::invalid_dof_index, const typename AffineConstraints< number >::ConstraintLine::Entries &entries={}, const number inhomogeneity=0.0)
std::vector< std::pair< size_type, number > > Entries
ConstraintLine(ConstraintLine &&other) noexcept=default
ConstraintLine & operator=(ConstraintLine &&other) noexcept=default
std::size_t memory_consumption() const
friend void swap(ConstraintLine &l1, ConstraintLine &l2) noexcept
typename internal::ProductTypeImpl< std::decay_t< T >, std::decay_t< U > >::type type
static void add(const typename VectorType::value_type value, const types::global_dof_index i, VectorType &V)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)