16#ifndef dealii_matrix_free_operators_h
17#define dealii_matrix_free_operators_h
47 template <
typename VectorType>
48 std::enable_if_t<IsBlockVector<VectorType>::value,
unsigned int>
51 return vector.n_blocks();
54 template <
typename VectorType>
55 std::enable_if_t<!IsBlockVector<VectorType>::value,
unsigned int>
61 template <
typename VectorType>
62 std::enable_if_t<IsBlockVector<VectorType>::value,
63 typename VectorType::BlockType &>
64 subblock(VectorType &vector,
unsigned int block_no)
67 return vector.block(block_no);
64 subblock(VectorType &vector,
unsigned int block_no) {
…}
70 template <
typename VectorType>
71 std::enable_if_t<IsBlockVector<VectorType>::value,
72 const typename VectorType::BlockType &>
73 subblock(
const VectorType &vector,
unsigned int block_no)
76 return vector.block(block_no);
73 subblock(
const VectorType &vector,
unsigned int block_no) {
…}
79 template <
typename VectorType>
80 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType &>
86 template <
typename VectorType>
87 std::enable_if_t<!IsBlockVector<VectorType>::value,
const VectorType &>
88 subblock(
const VectorType &vector,
unsigned int)
88 subblock(
const VectorType &vector,
unsigned int) {
…}
93 template <
typename VectorType>
94 std::enable_if_t<IsBlockVector<VectorType>::value,
void>
97 vector.collect_sizes();
100 template <
typename VectorType>
101 std::enable_if_t<!IsBlockVector<VectorType>::value,
void>
42 namespace BlockHelper {
…}
185 typename VectorizedArrayType =
236 const std::vector<unsigned int> &selected_row_blocks =
237 std::vector<unsigned int>(),
238 const std::vector<unsigned int> &selected_column_blocks =
239 std::vector<unsigned int>());
258 const unsigned int level,
259 const std::vector<unsigned int> &selected_row_blocks =
260 std::vector<unsigned int>());
279 const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
280 const unsigned int level,
281 const std::vector<unsigned int> &selected_row_blocks =
282 std::vector<unsigned int>());
312 vmult(VectorType &dst,
const VectorType &src)
const;
318 Tvmult(VectorType &dst,
const VectorType &src)
const;
324 vmult_add(VectorType &dst,
const VectorType &src)
const;
337 el(
const unsigned int row,
const unsigned int col)
const;
370 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
376 const std::shared_ptr<DiagonalMatrix<VectorType>> &
382 const std::shared_ptr<DiagonalMatrix<VectorType>> &
392 const VectorType &src,
421 apply_add(VectorType &dst,
const VectorType &src)
const = 0;
434 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
470 mutable std::vector<std::vector<std::pair<value_type, value_type>>>
485 const VectorType &src,
497 const bool is_row)
const;
536 template <
typename OperatorType>
570 template <
typename VectorType>
572 vmult(VectorType &dst,
const VectorType &src)
const;
577 template <
typename VectorType>
579 Tvmult(VectorType &dst,
const VectorType &src)
const;
584 template <
typename VectorType>
618 int n_components = 1,
619 typename Number = double,
624 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
625 "Type of Number and of VectorizedArrayType do not match.");
637 VectorizedArrayType> &
fe_eval);
650 const unsigned int n_actual_components,
651 const VectorizedArrayType *in_array,
652 VectorizedArrayType *out_array)
const;
666 apply(
const VectorizedArrayType *in_array,
667 VectorizedArrayType *out_array)
const;
683 &inverse_dyadic_coefficients,
684 const VectorizedArrayType *in_array,
685 VectorizedArrayType *out_array)
const;
722 const VectorizedArrayType *in_array,
723 VectorizedArrayType *out_array)
const;
756 int n_q_points_1d = fe_degree + 1,
757 int n_components = 1,
759 typename VectorizedArrayType =
809 const std::shared_ptr<DiagonalMatrix<VectorType>> &
815 const std::shared_ptr<DiagonalMatrix<VectorType>> &
825 apply_add(VectorType &dst,
const VectorType &src)
const override;
834 const VectorType &src,
835 const std::pair<unsigned int, unsigned int> &cell_range)
const;
864 int n_q_points_1d = fe_degree + 1,
865 int n_components = 1,
867 typename VectorizedArrayType =
965 std::shared_ptr<Table<2, VectorizedArrayType>>
975 apply_add(VectorType &dst,
const VectorType &src)
const override;
984 const VectorType &src,
985 const std::pair<unsigned int, unsigned int> &cell_range)
const;
995 const std::pair<unsigned int, unsigned int> &cell_range)
const;
1000 template <
int n_components_compute>
1005 n_components_compute,
1007 VectorizedArrayType> &phi,
1008 const unsigned int cell)
const;
1024 typename VectorizedArrayType>
1029 VectorizedArrayType>::
1030 CellwiseInverseMassMatrix(
1035 VectorizedArrayType> &fe_eval)
1029 VectorizedArrayType>:: {
…}
1045 typename VectorizedArrayType>
1051 VectorizedArrayType>::
1052 fill_inverse_JxW_values(
1055 const unsigned int dofs_per_component_on_cell =
1058 Utilities::pow(fe_eval.get_shape_info().data.front().fe_degree + 1,
1062 inverse_jxw.
size() % dofs_per_component_on_cell == 0,
1064 "Expected diagonal to be a multiple of scalar dof per cells"));
1067 for (
unsigned int q = 0; q < dofs_per_component_on_cell; ++q)
1068 inverse_jxw[q] = 1. / fe_eval.JxW(q);
1070 for (
unsigned int q = dofs_per_component_on_cell; q < inverse_jxw.
size();)
1071 for (
unsigned int i = 0; i < dofs_per_component_on_cell; ++i, ++q)
1072 inverse_jxw[q] = inverse_jxw[i];
1051 VectorizedArrayType>:: {
…}
1081 typename VectorizedArrayType>
1088 VectorizedArrayType>::apply(
const VectorizedArrayType *in_array,
1089 VectorizedArrayType *out_array)
const
1093 template run<fe_degree>(n_components, fe_eval, in_array, out_array);
1096 n_components, fe_eval, in_array, out_array);
1088 VectorizedArrayType>::apply(
const VectorizedArrayType *in_array, {
…}
1105 typename VectorizedArrayType>
1111 VectorizedArrayType>::
1113 const unsigned int n_actual_components,
1114 const VectorizedArrayType *in_array,
1115 VectorizedArrayType *out_array)
const
1120 VectorizedArrayType>::template run<fe_degree>(n_actual_components,
1123 inverse_coefficients),
1129 n_actual_components,
1111 VectorizedArrayType>:: {
…}
1141 typename VectorizedArrayType>
1147 VectorizedArrayType>::
1149 &inverse_dyadic_coefficients,
1150 const VectorizedArrayType *in_array,
1151 VectorizedArrayType *out_array)
const
1153 const unsigned int unrolled_size =
1154 inverse_dyadic_coefficients.size() * (n_components * n_components);
1158 VectorizedArrayType>::
1159 template run<fe_degree>(n_components,
1162 &inverse_dyadic_coefficients[0][0][0],
1172 &inverse_dyadic_coefficients[0][0][0], unrolled_size),
1147 VectorizedArrayType>:: {
…}
1184 typename VectorizedArrayType>
1190 VectorizedArrayType>::
1191 transform_from_q_points_to_basis(
const unsigned int n_actual_components,
1192 const VectorizedArrayType *in_array,
1193 VectorizedArrayType *out_array)
const
1195 const auto n_q_points_1d = fe_eval.get_shape_info().data[0].n_q_points_1d;
1197 if (fe_degree > -1 && (fe_degree + 1 == n_q_points_1d))
1200 VectorizedArrayType>::template run<fe_degree,
1201 fe_degree + 1>(n_actual_components,
1190 VectorizedArrayType>:: {
…}
1216 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1219 , have_interface_matrices(false)
1224 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1231 for (
const unsigned int selected_row : selected_rows)
1232 total_size +=
data->get_vector_partitioner(selected_row)->size();
1238 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1245 for (
const unsigned int selected_column : selected_columns)
1246 total_size +=
data->get_vector_partitioner(selected_column)->size();
1252 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1257 inverse_diagonal_entries.reset();
1262 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1265 const unsigned int col)
const
1268 Assert(inverse_diagonal_entries.get() !=
nullptr &&
1269 inverse_diagonal_entries->m() > 0,
1271 return 1.0 / (*inverse_diagonal_entries)(row, row);
1276 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1279 VectorType &vec)
const
1285 const unsigned int index = selected_rows[i];
1287 .partitioners_are_compatible(
1288 *
data->get_dof_info(index).vector_partitioner))
1292 .partitioners_are_globally_compatible(
1293 *
data->get_dof_info(index).vector_partitioner),
1301 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1306 const std::vector<unsigned int> &given_row_selection,
1307 const std::vector<unsigned int> &given_column_selection)
1311 selected_rows.clear();
1312 selected_columns.clear();
1313 if (given_row_selection.empty())
1314 for (
unsigned int i = 0; i < data_->n_components(); ++i)
1315 selected_rows.push_back(i);
1318 for (
unsigned int i = 0; i < given_row_selection.size(); ++i)
1321 for (
unsigned int j = 0; j < given_row_selection.size(); ++j)
1323 Assert(given_row_selection[j] != given_row_selection[i],
1324 ExcMessage(
"Given row indices must be unique"));
1326 selected_rows.push_back(given_row_selection[i]);
1329 if (given_column_selection.empty())
1330 selected_columns = selected_rows;
1333 for (
unsigned int i = 0; i < given_column_selection.size(); ++i)
1336 for (
unsigned int j = 0; j < given_column_selection.size(); ++j)
1338 Assert(given_column_selection[j] != given_column_selection[i],
1339 ExcMessage(
"Given column indices must be unique"));
1341 selected_columns.push_back(given_column_selection[i]);
1345 edge_constrained_indices.clear();
1346 edge_constrained_indices.resize(selected_rows.size());
1347 edge_constrained_values.clear();
1348 edge_constrained_values.resize(selected_rows.size());
1349 have_interface_matrices =
false;
1354 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1360 const unsigned int level,
1361 const std::vector<unsigned int> &given_row_selection)
1363 std::vector<MGConstrainedDoFs> mg_constrained_dofs_vector(
1364 1, mg_constrained_dofs);
1365 initialize(data_, mg_constrained_dofs_vector,
level, given_row_selection);
1370 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1375 const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
1376 const unsigned int level,
1377 const std::vector<unsigned int> &given_row_selection)
1382 selected_rows.clear();
1383 selected_columns.clear();
1384 if (given_row_selection.empty())
1385 for (
unsigned int i = 0; i < data_->n_components(); ++i)
1386 selected_rows.push_back(i);
1389 for (
unsigned int i = 0; i < given_row_selection.size(); ++i)
1392 for (
unsigned int j = 0; j < given_row_selection.size(); ++j)
1394 Assert(given_row_selection[j] != given_row_selection[i],
1395 ExcMessage(
"Given row indices must be unique"));
1397 selected_rows.push_back(given_row_selection[i]);
1400 selected_columns = selected_rows;
1403 edge_constrained_indices.clear();
1404 edge_constrained_indices.resize(selected_rows.size());
1405 edge_constrained_values.clear();
1406 edge_constrained_values.resize(selected_rows.size());
1410 for (
unsigned int j = 0; j < selected_rows.size(); ++j)
1412 if (data_->n_cell_batches() > 0)
1418 const std::vector<types::global_dof_index> interface_indices =
1419 mg_constrained_dofs[j]
1420 .get_refinement_edge_indices(
level)
1421 .get_index_vector();
1422 edge_constrained_indices[j].clear();
1423 edge_constrained_indices[j].reserve(interface_indices.size());
1424 edge_constrained_values[j].resize(interface_indices.size());
1426 data->get_dof_handler(selected_rows[j]).locally_owned_mg_dofs(
level);
1427 for (
const auto interface_index : interface_indices)
1428 if (locally_owned.
is_element(interface_index))
1429 edge_constrained_indices[j].push_back(
1431 have_interface_matrices |=
1433 static_cast<unsigned int>(edge_constrained_indices[j].
size()),
1434 data->get_vector_partitioner()->get_mpi_communicator()) > 0;
1440 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1443 VectorType &dst)
const
1447 const std::vector<unsigned int> &constrained_dofs =
1448 data->get_constrained_dofs(selected_rows[j]);
1449 for (
const auto constrained_dof : constrained_dofs)
1451 for (
unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1453 edge_constrained_indices[j][i]) = 1.;
1459 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1462 const VectorType &src)
const
1467 vmult_add(dst, src);
1472 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1476 const VectorType &src)
const
1478 mult_add(dst, src,
false);
1483 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1487 const VectorType &src)
const
1489 mult_add(dst, src,
true);
1494 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1497 const VectorType &src,
1498 const bool is_row)
const
1504 const unsigned int mf_component =
1505 is_row ? selected_rows[i] : selected_columns[i];
1508 data->get_dof_info(mf_component).vector_partitioner.get())
1516 data->get_dof_info(mf_component)
1517 .vector_partitioner->locally_owned_size(),
1519 "The vector passed to the vmult() function does not have "
1520 "the correct size for compatibility with MatrixFree."));
1526 this->
data->initialize_dof_vector(
1530 .copy_locally_owned_data_from(copy_vec);
1536 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1540 const VectorType &src)
const
1544 adjust_ghost_range_if_necessary(src,
false);
1545 adjust_ghost_range_if_necessary(dst,
true);
1551 for (
unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1553 edge_constrained_values[j][i] = std::pair<Number, Number>(
1555 edge_constrained_indices[j][i]),
1557 edge_constrained_indices[j][i]));
1559 .local_element(edge_constrained_indices[j][i]) = 0.;
1566 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1570 const VectorType &src,
1576 preprocess_constraints(dst, src);
1578 Tapply_add(dst, src);
1580 apply_add(dst, src);
1581 postprocess_constraints(dst, src);
1586 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1590 const VectorType &src)
const
1594 const std::vector<unsigned int> &constrained_dofs =
1595 data->get_constrained_dofs(selected_rows[j]);
1596 for (
const auto constrained_dof : constrained_dofs)
1605 for (
unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1608 .local_element(edge_constrained_indices[j][i]) =
1609 edge_constrained_values[j][i].first;
1611 edge_constrained_indices[j][i]) =
1612 edge_constrained_values[j][i].second +
1613 edge_constrained_values[j][i].first;
1620 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1624 const VectorType &src)
const
1629 adjust_ghost_range_if_necessary(src,
false);
1630 adjust_ghost_range_if_necessary(dst,
true);
1634 if (!have_interface_matrices)
1640 for (
unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1642 edge_constrained_values[j][i] = std::pair<Number, Number>(
1644 edge_constrained_indices[j][i]),
1646 edge_constrained_indices[j][i]));
1648 .local_element(edge_constrained_indices[j][i]) = 0.;
1651 apply_add(dst, src);
1656 for (
unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1658 for (; c < edge_constrained_indices[j][i]; ++c)
1664 .local_element(edge_constrained_indices[j][i]) =
1665 edge_constrained_values[j][i].first;
1674 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1678 const VectorType &src)
const
1683 adjust_ghost_range_if_necessary(src,
false);
1684 adjust_ghost_range_if_necessary(dst,
true);
1688 if (!have_interface_matrices)
1691 VectorType src_cpy(src);
1695 for (
unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1697 for (; c < edge_constrained_indices[j][i]; ++c)
1705 apply_add(dst, src_cpy);
1708 for (
unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1710 edge_constrained_indices[j][i]) = 0.;
1715 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1719 const VectorType &src)
const
1724 Tvmult_add(dst, src);
1729 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1733 return inverse_diagonal_entries.get() !=
nullptr ?
1734 inverse_diagonal_entries->memory_consumption() :
1740 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1744 VectorizedArrayType>>
1752 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1753 const std::shared_ptr<DiagonalMatrix<VectorType>> &
1757 Assert(inverse_diagonal_entries.get() !=
nullptr &&
1758 inverse_diagonal_entries->m() > 0,
1760 return inverse_diagonal_entries;
1765 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1766 const std::shared_ptr<DiagonalMatrix<VectorType>> &
1769 Assert(diagonal_entries.get() !=
nullptr && diagonal_entries->m() > 0,
1771 return diagonal_entries;
1776 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1780 const VectorType &src)
const
1782 apply_add(dst, src);
1787 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1791 const VectorType &src,
1795 Assert(inverse_diagonal_entries.get() && inverse_diagonal_entries->m() > 0,
1797 inverse_diagonal_entries->vmult(dst, src);
1805 template <
typename OperatorType>
1808 , mf_base_operator(nullptr)
1813 template <
typename OperatorType>
1817 mf_base_operator =
nullptr;
1822 template <
typename OperatorType>
1826 mf_base_operator = &operator_in;
1831 template <
typename OperatorType>
1832 template <
typename VectorType>
1835 const VectorType &src)
const
1839 std::is_same_v<typename VectorType::value_type, value_type>,
1840 "The vector type must be based on the same value type as this "
1846 mf_base_operator->vmult_interface_down(dst, src);
1851 template <
typename OperatorType>
1852 template <
typename VectorType>
1855 const VectorType &src)
const
1859 std::is_same_v<typename VectorType::value_type, value_type>,
1860 "The vector type must be based on the same value type as this "
1866 mf_base_operator->vmult_interface_up(dst, src);
1871 template <
typename OperatorType>
1872 template <
typename VectorType>
1875 VectorType &vec)
const
1879 mf_base_operator->initialize_dof_vector(vec);
1890 typename VectorType,
1891 typename VectorizedArrayType>
1897 VectorizedArrayType>::MassOperator()
1898 :
Base<dim, VectorType, VectorizedArrayType>()
1903 "This class only supports the non-blocked vector variant of the Base "
1904 "operator because only a single FEEvaluation object is used in the "
1905 "apply function."));
1897 VectorizedArrayType>::MassOperator() {
…}
1914 typename VectorType,
1915 typename VectorizedArrayType>
1922 VectorizedArrayType>::compute_diagonal()
1926 Assert(this->selected_rows == this->selected_columns,
1927 ExcMessage(
"This function is only implemented for square (not "
1928 "rectangular) operators."));
1930 this->inverse_diagonal_entries =
1931 std::make_shared<DiagonalMatrix<VectorType>>();
1932 this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1933 VectorType &inverse_diagonal_vector =
1934 this->inverse_diagonal_entries->get_vector();
1935 VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1936 this->initialize_dof_vector(inverse_diagonal_vector);
1937 this->initialize_dof_vector(diagonal_vector);
1941 auto diagonal_evaluation = [](
auto &integrator) {
1943 for (
unsigned int q = 0; q < integrator.n_q_points; ++q)
1944 integrator.submit_value(integrator.get_value(q), q);
1955 VectorizedArrayType> &)>
1956 diagonal_evaluation_f(diagonal_evaluation);
1959 for (
unsigned int block_n = 0; block_n < this->selected_rows.size();
1964 diagonal_evaluation_f,
1965 this->selected_rows[block_n]);
1969 this->set_constrained_entries_to_one(diagonal_vector);
1971 inverse_diagonal_vector = diagonal_vector;
1973 for (
unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
1981 Assert(diagonal_vector.local_element(i) > Number(0),
1984 inverse_diagonal_vector.local_element(i) =
1985 1. / inverse_diagonal_vector.local_element(i);
1922 VectorizedArrayType>::compute_diagonal() {
…}
1997 typename VectorType,
1998 typename VectorizedArrayType>
2005 VectorizedArrayType>::compute_lumped_diagonal()
2011 Assert(this->selected_rows == this->selected_columns,
2012 ExcMessage(
"This function is only implemented for square (not "
2013 "rectangular) operators."));
2015 inverse_lumped_diagonal_entries =
2016 std::make_shared<DiagonalMatrix<VectorType>>();
2017 lumped_diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
2018 VectorType &inverse_lumped_diagonal_vector =
2019 inverse_lumped_diagonal_entries->get_vector();
2020 VectorType &lumped_diagonal_vector = lumped_diagonal_entries->get_vector();
2021 this->initialize_dof_vector(inverse_lumped_diagonal_vector);
2022 this->initialize_dof_vector(lumped_diagonal_vector);
2025 inverse_lumped_diagonal_vector = Number(1.);
2026 apply_add(lumped_diagonal_vector, inverse_lumped_diagonal_vector);
2027 this->set_constrained_entries_to_one(lumped_diagonal_vector);
2030 inverse_lumped_diagonal_vector.locally_owned_size();
2037 if (lumped_diagonal_vector.local_element(i) == Number(0.))
2038 inverse_lumped_diagonal_vector.local_element(i) = Number(1.);
2040 inverse_lumped_diagonal_vector.local_element(i) =
2041 Number(1.) / lumped_diagonal_vector.local_element(i);
2044 inverse_lumped_diagonal_vector.update_ghost_values();
2045 lumped_diagonal_vector.update_ghost_values();
2005 VectorizedArrayType>::compute_lumped_diagonal() {
…}
2054 typename VectorType,
2055 typename VectorizedArrayType>
2056 const std::shared_ptr<DiagonalMatrix<VectorType>> &
2062 VectorizedArrayType>::get_matrix_lumped_diagonal_inverse()
const
2064 Assert(inverse_lumped_diagonal_entries.get() !=
nullptr &&
2065 inverse_lumped_diagonal_entries->m() > 0,
2067 return inverse_lumped_diagonal_entries;
2062 VectorizedArrayType>::get_matrix_lumped_diagonal_inverse()
const {
…}
2076 typename VectorType,
2077 typename VectorizedArrayType>
2078 const std::shared_ptr<DiagonalMatrix<VectorType>> &
2084 VectorizedArrayType>::get_matrix_lumped_diagonal()
const
2086 Assert(lumped_diagonal_entries.get() !=
nullptr &&
2087 lumped_diagonal_entries->m() > 0,
2089 return lumped_diagonal_entries;
2084 VectorizedArrayType>::get_matrix_lumped_diagonal()
const {
…}
2098 typename VectorType,
2099 typename VectorizedArrayType>
2106 VectorizedArrayType>::apply_add(VectorType &dst,
2107 const VectorType &src)
const
2106 VectorizedArrayType>::apply_add(VectorType &dst, {
…}
2119 typename VectorType,
2120 typename VectorizedArrayType>
2127 VectorizedArrayType>::
2132 VectorizedArrayType> &
data,
2134 const VectorType &src,
2135 const std::pair<unsigned int, unsigned int> &cell_range)
const
2144 VectorizedArrayType>
2145 phi(
data, this->selected_rows[0]);
2146 for (
unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2149 phi.read_dof_values(src);
2151 for (
unsigned int q = 0; q < phi.n_q_points; ++q)
2152 phi.submit_value(phi.get_value(q), q);
2154 phi.distribute_local_to_global(dst);
2127 VectorizedArrayType>:: {
…}
2165 typename VectorType,
2166 typename VectorizedArrayType>
2172 VectorizedArrayType>::LaplaceOperator()
2173 :
Base<dim, VectorType, VectorizedArrayType>()
2172 VectorizedArrayType>::LaplaceOperator() {
…}
2182 typename VectorType,
2183 typename VectorizedArrayType>
2190 VectorizedArrayType>::clear()
2193 scalar_coefficient.reset();
2190 VectorizedArrayType>::clear() {
…}
2202 typename VectorType,
2203 typename VectorizedArrayType>
2210 VectorizedArrayType>::
2214 scalar_coefficient = scalar_coefficient_;
2210 VectorizedArrayType>:: {
…}
2223 typename VectorType,
2224 typename VectorizedArrayType>
2225 std::shared_ptr<Table<2, VectorizedArrayType>>
2231 VectorizedArrayType>::get_coefficient()
2234 return scalar_coefficient;
2231 VectorizedArrayType>::get_coefficient() {
…}
2243 typename VectorType,
2244 typename VectorizedArrayType>
2251 VectorizedArrayType>::compute_diagonal()
2258 this->inverse_diagonal_entries =
2259 std::make_shared<DiagonalMatrix<VectorType>>();
2260 this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
2261 VectorType &inverse_diagonal_vector =
2262 this->inverse_diagonal_entries->get_vector();
2263 VectorType &diagonal_vector = this->diagonal_entries->get_vector();
2264 this->initialize_dof_vector(inverse_diagonal_vector);
2265 this->initialize_dof_vector(diagonal_vector);
2271 this->set_constrained_entries_to_one(diagonal_vector);
2273 inverse_diagonal_vector = diagonal_vector;
2275 for (
unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
2277 if (
std::abs(inverse_diagonal_vector.local_element(i)) >
2278 std::sqrt(std::numeric_limits<Number>::epsilon()))
2279 inverse_diagonal_vector.local_element(i) =
2280 1. / inverse_diagonal_vector.local_element(i);
2282 inverse_diagonal_vector.local_element(i) = 1.;
2251 VectorizedArrayType>::compute_diagonal() {
…}
2293 typename VectorType,
2294 typename VectorizedArrayType>
2301 VectorizedArrayType>::apply_add(VectorType &dst,
2302 const VectorType &src)
const
2301 VectorizedArrayType>::apply_add(VectorType &dst, {
…}
2308 namespace Implementation
2310 template <
typename VectorizedArrayType>
2314 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2308 namespace Implementation {
…}
2328 typename VectorType,
2329 typename VectorizedArrayType>
2330 template <
int n_components_compute>
2332 LaplaceOperator<dim,
2337 VectorizedArrayType>::
2338 do_operation_on_cell(
2343 n_components_compute,
2345 VectorizedArrayType> &phi,
2346 const unsigned int cell)
const
2349 if (scalar_coefficient.get())
2351 Assert(scalar_coefficient->size(1) == 1 ||
2352 scalar_coefficient->size(1) == phi.n_q_points,
2353 ExcMessage(
"The number of columns in the coefficient table must "
2354 "be either 1 or the number of quadrature points " +
2355 std::to_string(phi.n_q_points) +
2356 ", but the given value was " +
2357 std::to_string(scalar_coefficient->size(1))));
2358 if (scalar_coefficient->size(1) == phi.n_q_points)
2359 for (
unsigned int q = 0; q < phi.n_q_points; ++q)
2362 (*scalar_coefficient)(cell, q)),
2363 ExcMessage(
"Coefficient must be non-negative"));
2364 phi.submit_gradient((*scalar_coefficient)(cell, q) *
2365 phi.get_gradient(q),
2371 ExcMessage(
"Coefficient must be non-negative"));
2372 const VectorizedArrayType coefficient =
2373 (*scalar_coefficient)(cell, 0);
2374 for (
unsigned int q = 0; q < phi.n_q_points; ++q)
2375 phi.submit_gradient(coefficient * phi.get_gradient(q), q);
2380 for (
unsigned int q = 0; q < phi.n_q_points; ++q)
2382 phi.submit_gradient(phi.get_gradient(q), q);
2337 VectorizedArrayType>:: {
…}
2394 typename VectorType,
2395 typename VectorizedArrayType>
2402 VectorizedArrayType>::
2407 VectorizedArrayType> &
data,
2409 const VectorType &src,
2410 const std::pair<unsigned int, unsigned int> &cell_range)
const
2419 VectorizedArrayType>
2420 phi(
data, this->selected_rows[0]);
2421 for (
unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2424 phi.read_dof_values(src);
2425 do_operation_on_cell(phi, cell);
2426 phi.distribute_local_to_global(dst);
2402 VectorizedArrayType>:: {
…}
2435 typename VectorType,
2436 typename VectorizedArrayType>
2443 VectorizedArrayType>::
2444 local_diagonal_cell(
2448 VectorizedArrayType> &
data,
2451 const std::pair<unsigned int, unsigned int> &cell_range)
const
2457 eval(
data, this->selected_rows[0]);
2463 VectorizedArrayType>
2464 eval_vector(
data, this->selected_rows[0]);
2465 for (
unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2468 eval_vector.reinit(cell);
2480 do_operation_on_cell(eval, cell);
2485 for (
unsigned int c = 0; c < n_components; ++c)
2487 .begin_dof_values()[i + c * eval_vector.dofs_per_component] =
2490 eval_vector.distribute_local_to_global(dst);
2443 VectorizedArrayType>:: {
…}
337 el(
const unsigned int row,
const unsigned int col)
const; {
…}
318 Tvmult(VectorType &dst,
const VectorType &src)
const; {
…}
312 vmult(VectorType &dst,
const VectorType &src)
const; {
…}
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
const Number * begin_dof_values() const
void reinit(const unsigned int cell_batch_index)
const unsigned int dofs_per_cell
size_type index_within_set(const size_type global_index) const
bool is_element(const size_type index) const
virtual ~Base() override=default
std::vector< unsigned int > selected_rows
void Tvmult_add(VectorType &dst, const VectorType &src) const
void set_constrained_entries_to_one(VectorType &dst) const
virtual void compute_diagonal()=0
void vmult_add(VectorType &dst, const VectorType &src) const
void vmult_interface_down(VectorType &dst, const VectorType &src) const
void preprocess_constraints(VectorType &dst, const VectorType &src) const
void mult_add(VectorType &dst, const VectorType &src, const bool transpose) const
std::vector< std::vector< unsigned int > > edge_constrained_indices
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal() const
void Tvmult(VectorType &dst, const VectorType &src) const
std::vector< std::vector< std::pair< value_type, value_type > > > edge_constrained_values
void initialize(std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > data, const std::vector< unsigned int > &selected_row_blocks=std::vector< unsigned int >(), const std::vector< unsigned int > &selected_column_blocks=std::vector< unsigned int >())
bool have_interface_matrices
std::vector< unsigned int > selected_columns
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > get_matrix_free() const
void initialize(std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > data_, const std::vector< MGConstrainedDoFs > &mg_constrained_dofs, const unsigned int level, const std::vector< unsigned int > &selected_row_blocks=std::vector< unsigned int >())
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal_inverse() const
void initialize(std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > data, const MGConstrainedDoFs &mg_constrained_dofs, const unsigned int level, const std::vector< unsigned int > &selected_row_blocks=std::vector< unsigned int >())
void initialize_dof_vector(VectorType &vec) const
std::shared_ptr< DiagonalMatrix< VectorType > > diagonal_entries
virtual void Tapply_add(VectorType &dst, const VectorType &src) const
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > data
void adjust_ghost_range_if_necessary(const VectorType &vec, const bool is_row) const
virtual std::size_t memory_consumption() const
std::shared_ptr< DiagonalMatrix< VectorType > > inverse_diagonal_entries
value_type el(const unsigned int row, const unsigned int col) const
virtual void apply_add(VectorType &dst, const VectorType &src) const =0
typename VectorType::size_type size_type
void precondition_Jacobi(VectorType &dst, const VectorType &src, const value_type omega) const
typename VectorType::value_type value_type
void vmult(VectorType &dst, const VectorType &src) const
void vmult_interface_up(VectorType &dst, const VectorType &src) const
void postprocess_constraints(VectorType &dst, const VectorType &src) const
const FEEvaluationBase< dim, n_components, Number, false, VectorizedArrayType > & fe_eval
void fill_inverse_JxW_values(AlignedVector< VectorizedArrayType > &inverse_jxw) const
void transform_from_q_points_to_basis(const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
void apply(const AlignedVector< VectorizedArrayType > &inverse_coefficient, const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
std::shared_ptr< Table< 2, VectorizedArrayType > > get_coefficient()
typename Base< dim, VectorType, VectorizedArrayType >::value_type value_type
virtual void compute_diagonal() override
std::shared_ptr< Table< 2, VectorizedArrayType > > scalar_coefficient
void local_diagonal_cell(const MatrixFree< dim, value_type, VectorizedArrayType > &data, VectorType &dst, const VectorType &, const std::pair< unsigned int, unsigned int > &cell_range) const
void local_apply_cell(const MatrixFree< dim, value_type, VectorizedArrayType > &data, VectorType &dst, const VectorType &src, const std::pair< unsigned int, unsigned int > &cell_range) const
virtual void apply_add(VectorType &dst, const VectorType &src) const override
typename Base< dim, VectorType, VectorizedArrayType >::size_type size_type
void set_coefficient(const std::shared_ptr< Table< 2, VectorizedArrayType > > &scalar_coefficient)
void do_operation_on_cell(FEEvaluation< dim, fe_degree, n_q_points_1d, n_components_compute, value_type, VectorizedArrayType > &phi, const unsigned int cell) const
virtual void clear() override
typename OperatorType::value_type value_type
void vmult(VectorType &dst, const VectorType &src) const
ObserverPointer< const OperatorType > mf_base_operator
void initialize(const OperatorType &operator_in)
void Tvmult(VectorType &dst, const VectorType &src) const
void initialize_dof_vector(VectorType &vec) const
typename OperatorType::size_type size_type
void compute_lumped_diagonal()
std::shared_ptr< DiagonalMatrix< VectorType > > lumped_diagonal_entries
typename Base< dim, VectorType, VectorizedArrayType >::size_type size_type
virtual void apply_add(VectorType &dst, const VectorType &src) const override
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_lumped_diagonal() const
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_lumped_diagonal_inverse() const
void local_apply_cell(const MatrixFree< dim, value_type, VectorizedArrayType > &data, VectorType &dst, const VectorType &src, const std::pair< unsigned int, unsigned int > &cell_range) const
virtual void compute_diagonal() override
typename Base< dim, VectorType, VectorizedArrayType >::value_type value_type
std::shared_ptr< DiagonalMatrix< VectorType > > inverse_lumped_diagonal_entries
#define DEAL_II_NAMESPACE_OPEN
constexpr bool running_in_debug_mode()
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
std::vector< index_type > data
types::global_dof_index locally_owned_size
std::enable_if_t< IsBlockVector< VectorType >::value, typename VectorType::BlockType & > subblock(VectorType &vector, unsigned int block_no)
std::enable_if_t< IsBlockVector< VectorType >::value, void > collect_sizes(VectorType &vector)
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
bool non_negative(const VectorizedArrayType &n)
T max(const T &t, const MPI_Comm mpi_communicator)
constexpr T pow(const T base, const int iexp)
constexpr unsigned int invalid_unsigned_int
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
static void apply(const unsigned int n_components, const FEEvaluationData< dim, Number, false > &fe_eval, const Number *in_array, Number *out_array)
static void transform_from_q_points_to_basis(const unsigned int n_components, const FEEvaluationData< dim, Number, false > &fe_eval, const Number *in_array, Number *out_array)