16 #ifndef dealii_affine_constraints_h
17 #define dealii_affine_constraints_h
30 #include <boost/range/iterator_range.hpp>
46 template <
typename number>
48 template <
typename number>
53 template <
typename number>
54 class GlobalRowsFromLocal;
59 namespace AffineConstraintsImplementation
61 template <
class VectorType>
63 set_zero_all(
const std::vector<types::global_dof_index> &cm,
68 set_zero_all(
const std::vector<types::global_dof_index> &cm,
73 set_zero_all(
const std::vector<types::global_dof_index> &cm,
79 template <
typename number>
179 template <
typename number =
double>
400 const std::vector<std::pair<size_type, number>> &col_val_pairs);
462 const bool allow_different_local_lines =
false);
572 const std::vector<std::pair<size_type, number>> *
603 print(std::ostream &out)
const;
711 template <
class VectorType>
721 template <
class VectorType>
737 template <
class VectorType>
745 template <
class BlockVectorType>
755 template <
class VectorType>
810 template <
class InVector,
class OutVector>
813 const std::vector<size_type> &local_dof_indices,
814 OutVector & global_vector)
const;
859 template <
typename VectorType>
862 const std::vector<size_type> &local_dof_indices,
885 template <
typename VectorType>
889 const std::vector<size_type> &local_dof_indices_row,
890 const std::vector<size_type> &local_dof_indices_col,
898 template <
class VectorType>
932 template <
typename ForwardIteratorVec,
933 typename ForwardIteratorInd,
937 ForwardIteratorVec local_vector_end,
938 ForwardIteratorInd local_indices_begin,
988 template <
typename MatrixType>
991 const std::vector<size_type> &local_dof_indices,
992 MatrixType & global_matrix)
const;
1021 template <
typename MatrixType>
1024 const std::vector<size_type> &row_indices,
1025 const std::vector<size_type> &col_indices,
1026 MatrixType & global_matrix)
const;
1044 template <
typename MatrixType>
1047 const std::vector<size_type> &row_indices,
1049 const std::vector<size_type> &column_indices,
1050 MatrixType & global_matrix)
const;
1068 template <
typename MatrixType,
typename VectorType>
1072 const std::vector<size_type> &local_dof_indices,
1073 MatrixType & global_matrix,
1075 bool use_inhomogeneities_for_rhs =
false)
const;
1130 template <
typename SparsityPatternType>
1133 const std::vector<size_type> &local_dof_indices,
1134 SparsityPatternType & sparsity_pattern,
1135 const bool keep_constrained_entries =
true,
1141 template <
typename SparsityPatternType>
1144 const std::vector<size_type> &row_indices,
1145 const std::vector<size_type> &col_indices,
1146 SparsityPatternType & sparsity_pattern,
1147 const bool keep_constrained_entries =
true,
1169 template <
typename ForwardIteratorVec,
1170 typename ForwardIteratorInd,
1174 ForwardIteratorInd local_indices_begin,
1175 ForwardIteratorVec local_vector_begin,
1176 ForwardIteratorVec local_vector_end)
const;
1199 template <
class VectorType>
1216 using Entries = std::vector<std::pair<size_type, number>>;
1267 template <
class Archive>
1325 const IndexSet & locally_active_dofs,
1327 const bool verbose =
false)
const;
1348 <<
"The specified line " << arg1 <<
" does not exist.");
1359 <<
"The entry for the indices " << arg1 <<
" and " << arg2
1360 <<
" already exists, but the values " << arg3 <<
" (old) and "
1361 << arg4 <<
" (new) differ "
1362 <<
"by " << (arg4 - arg3) <<
".");
1371 <<
"You tried to constrain DoF " << arg1 <<
" to DoF " << arg2
1372 <<
", but that one is also constrained. This is not allowed!");
1380 <<
"Degree of freedom " << arg1
1381 <<
" is constrained from both object in a merge operation.");
1389 <<
"In the given argument a degree of freedom is constrained "
1390 <<
"to another DoF with number " << arg1
1391 <<
", which however is constrained by this object. This is not"
1400 <<
"The index set given to this constraints object indicates "
1401 <<
"constraints for degree of freedom " << arg1
1402 <<
" should not be stored by this object, but a constraint "
1403 <<
"is being added.");
1413 <<
"The index set given to this constraints object indicates "
1414 <<
"constraints using degree of freedom " << arg2
1415 <<
" should not be stored by this object, but a constraint "
1416 <<
"for degree of freedom " << arg1 <<
" uses it.");
1426 <<
"While distributing the constraint for DoF " << arg1
1427 <<
", it turns out that one of the processors "
1428 <<
"who own the " << arg2 <<
" degrees of freedom that x_"
1429 << arg1 <<
" is constrained against does not know about "
1430 <<
"the constraint on x_" << arg1
1431 <<
". Did you not initialize the AffineConstraints container "
1432 <<
"with the appropriate locally_relevant set so "
1433 <<
"that every processor who owns a DoF that constrains "
1434 <<
"another DoF also knows about this constraint?");
1507 template <
typename MatrixType,
typename VectorType>
1511 const std::vector<size_type> &local_dof_indices,
1512 MatrixType & global_matrix,
1514 bool use_inhomogeneities_for_rhs,
1515 std::integral_constant<bool, false>)
const;
1521 template <
typename MatrixType,
typename VectorType>
1525 const std::vector<size_type> &local_dof_indices,
1526 MatrixType & global_matrix,
1528 bool use_inhomogeneities_for_rhs,
1529 std::integral_constant<bool, true>)
const;
1535 template <
typename SparsityPatternType>
1538 SparsityPatternType & sparsity_pattern,
1539 const bool keep_constrained_entries,
1541 std::integral_constant<bool, false>)
const;
1547 template <
typename SparsityPatternType>
1550 SparsityPatternType & sparsity_pattern,
1551 const bool keep_constrained_entries,
1553 std::integral_constant<bool, true>)
const;
1564 const std::vector<size_type> & local_dof_indices,
1565 internals::GlobalRowsFromLocal<number> &global_rows)
const;
1576 std::vector<size_type> & active_dofs)
const;
1581 template <
typename MatrixScalar,
typename VectorScalar>
1585 const internals::GlobalRowsFromLocal<number> &global_rows,
1586 const Vector<VectorScalar> & local_vector,
1587 const std::vector<size_type> & local_dof_indices,
1593 template <
typename number>
1597 , local_lines(local_constraints)
1606 template <
typename number>
1610 , lines(affine_constraints.lines)
1611 , lines_cache(affine_constraints.lines_cache)
1612 , local_lines(affine_constraints.local_lines)
1613 , sorted(affine_constraints.sorted)
1616 template <
typename number>
1620 Assert(sorted ==
false, ExcMatrixIsClosed());
1627 const size_type line_index = calculate_line_index(line_n);
1630 if (is_constrained(line_n))
1634 if (line_index >= lines_cache.size())
1640 lines.emplace_back();
1641 lines.back().index = line_n;
1642 lines.back().inhomogeneity = 0.;
1643 lines_cache[line_index] = lines.size() - 1;
1646 template <
typename number>
1652 Assert(sorted ==
false, ExcMatrixIsClosed());
1654 ExcMessage(
"Can't constrain a degree of freedom to itself"));
1657 const size_type line_index = calculate_line_index(line_n);
1658 Assert(line_index < lines_cache.size(),
1659 ExcMessage(
"The current AffineConstraints does not contain the line "
1660 "for the current entry. Call AffineConstraints::add_line "
1661 "before calling this function."));
1670 Assert(!local_lines.size() || local_lines.is_element(column),
1671 ExcColumnNotStoredHere(line_n, column));
1672 ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
1674 for (
const auto &p : line_ptr->entries)
1675 if (p.first == column)
1678 ExcEntryAlreadyExists(line_n, column, p.second,
value));
1682 line_ptr->entries.emplace_back(column,
value);
1685 template <
typename number>
1690 const size_type line_index = calculate_line_index(line_n);
1691 Assert(line_index < lines_cache.size() &&
1693 ExcMessage(
"call add_line() before calling set_inhomogeneity()"));
1695 ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
1696 line_ptr->inhomogeneity =
value;
1699 template <
typename number>
1700 template <
class VectorType>
1706 std::vector<size_type> constrained_lines(lines.size());
1707 for (
unsigned int i = 0; i < lines.size(); ++i)
1708 constrained_lines[i] = lines[i].index;
1709 internal::AffineConstraintsImplementation::set_zero_all(constrained_lines,
1713 template <
typename number>
1717 return lines.size();
1720 template <
typename number>
1724 const size_type line_index = calculate_line_index(index);
1725 return ((line_index < lines_cache.size()) &&
1729 template <
typename number>
1736 const size_type line_index = calculate_line_index(line_n);
1737 if (line_index >= lines_cache.size() ||
1743 return !(lines[lines_cache[line_index]].inhomogeneity == number(0.));
1747 template <
typename number>
1748 inline const std::vector<std::pair<types::global_dof_index, number>> *
1753 const size_type line_index = calculate_line_index(line_n);
1754 if (line_index >= lines_cache.size() ||
1758 return &lines[lines_cache[line_index]].entries;
1761 template <
typename number>
1767 const size_type line_index = calculate_line_index(line_n);
1768 if (line_index >= lines_cache.size() ||
1772 return lines[lines_cache[line_index]].inhomogeneity;
1775 template <
typename number>
1780 if (!local_lines.size())
1783 Assert(local_lines.is_element(line_n), ExcRowNotStoredHere(line_n));
1785 return local_lines.index_within_set(line_n);
1788 template <
typename number>
1792 return !local_lines.size() || local_lines.is_element(line_n);
1795 template <
typename number>
1802 template <
typename number>
1803 template <
class VectorType>
1810 Assert(lines.empty() || sorted ==
true, ExcMatrixNotClosed());
1812 if (is_constrained(index) ==
false)
1813 global_vector(index) +=
value;
1816 const ConstraintLine &position =
1817 lines[lines_cache[calculate_line_index(index)]];
1818 for (
size_type j = 0; j < position.entries.size(); ++j)
1819 global_vector(position.entries[j].first) +=
1820 value * position.entries[j].second;
1824 template <
typename number>
1825 template <
typename ForwardIteratorVec,
1826 typename ForwardIteratorInd,
1830 ForwardIteratorVec local_vector_begin,
1831 ForwardIteratorVec local_vector_end,
1832 ForwardIteratorInd local_indices_begin,
1835 Assert(lines.empty() || sorted ==
true, ExcMatrixNotClosed());
1836 for (; local_vector_begin != local_vector_end;
1837 ++local_vector_begin, ++local_indices_begin)
1839 if (is_constrained(*local_indices_begin) ==
false)
1841 *local_indices_begin,
1845 const ConstraintLine &position =
1846 lines[lines_cache[calculate_line_index(*local_indices_begin)]];
1847 for (
size_type j = 0; j < position.entries.size(); ++j)
1849 (*local_vector_begin) * position.entries[j].second,
1850 position.entries[j].first,
1856 template <
typename number>
1857 template <
class InVector,
class OutVector>
1860 const InVector & local_vector,
1861 const std::vector<size_type> &local_dof_indices,
1862 OutVector & global_vector)
const
1864 Assert(local_vector.size() == local_dof_indices.size(),
1866 distribute_local_to_global(local_vector.begin(),
1868 local_dof_indices.begin(),
1872 template <
typename number>
1873 template <
typename ForwardIteratorVec,
1874 typename ForwardIteratorInd,
1879 ForwardIteratorInd local_indices_begin,
1880 ForwardIteratorVec local_vector_begin,
1881 ForwardIteratorVec local_vector_end)
const
1883 Assert(lines.empty() || sorted ==
true, ExcMatrixNotClosed());
1884 for (; local_vector_begin != local_vector_end;
1885 ++local_vector_begin, ++local_indices_begin)
1887 if (is_constrained(*local_indices_begin) ==
false)
1888 *local_vector_begin = global_vector(*local_indices_begin);
1891 const ConstraintLine &position =
1892 lines[lines_cache[calculate_line_index(*local_indices_begin)]];
1893 typename VectorType::value_type
value = position.inhomogeneity;
1894 for (
size_type j = 0; j < position.entries.size(); ++j)
1895 value += (global_vector(position.entries[j].first) *
1896 position.entries[j].second);
1897 *local_vector_begin =
value;
1902 template <
typename MatrixType>
1904 template <
typename SparsityPatternType>
1906 template <
typename number>
1927 template <
typename MatrixType>
1945 template <
typename T>
1953 template <
typename T>
1961 template <
typename T>
1984 template <
typename MatrixType>
1987 template <
typename number>
1988 template <
typename MatrixType>
1992 const std::vector<size_type> &local_dof_indices,
1993 MatrixType & global_matrix)
const
1997 Vector<typename MatrixType::value_type>
dummy(0);
1998 distribute_local_to_global(
2008 template <
typename number>
2009 template <
typename MatrixType,
typename VectorType>
2014 const std::vector<size_type> &local_dof_indices,
2015 MatrixType & global_matrix,
2017 bool use_inhomogeneities_for_rhs)
const
2021 distribute_local_to_global(
2027 use_inhomogeneities_for_rhs,
2031 template <
typename number>
2032 template <
typename SparsityPatternType>
2035 const std::vector<size_type> &local_dof_indices,
2036 SparsityPatternType & sparsity_pattern,
2037 const bool keep_constrained_entries,
2042 add_entries_local_to_global(
2045 keep_constrained_entries,