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);
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);
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)
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>
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)
1039 fe_eval.get_shape_info().n_q_points);
1048 typename VectorizedArrayType>
1054 VectorizedArrayType>::
1055 fill_inverse_JxW_values(
1058 const unsigned int dofs_per_component_on_cell =
1061 Utilities::pow(fe_eval.get_shape_info().data.front().fe_degree + 1,
1065 inverse_jxw.
size() % dofs_per_component_on_cell == 0,
1067 "Expected diagonal to be a multiple of scalar dof per cells"));
1070 for (
unsigned int q = 0; q < dofs_per_component_on_cell; ++q)
1071 inverse_jxw[q] = 1. / fe_eval.JxW(q);
1073 for (
unsigned int q = dofs_per_component_on_cell; q < inverse_jxw.
size();)
1074 for (
unsigned int i = 0; i < dofs_per_component_on_cell; ++i, ++q)
1075 inverse_jxw[q] = inverse_jxw[i];
1084 typename VectorizedArrayType>
1091 VectorizedArrayType>::apply(
const VectorizedArrayType *in_array,
1092 VectorizedArrayType *out_array)
const
1096 template run<fe_degree>(n_components, fe_eval, in_array, out_array);
1099 n_components, fe_eval, in_array, out_array);
1108 typename VectorizedArrayType>
1114 VectorizedArrayType>::
1116 const unsigned int n_actual_components,
1117 const VectorizedArrayType *in_array,
1118 VectorizedArrayType *out_array)
const
1123 VectorizedArrayType>::template run<fe_degree>(n_actual_components,
1126 inverse_coefficients),
1132 n_actual_components,
1144 typename VectorizedArrayType>
1150 VectorizedArrayType>::
1152 &inverse_dyadic_coefficients,
1153 const VectorizedArrayType *in_array,
1154 VectorizedArrayType *out_array)
const
1156 const unsigned int unrolled_size =
1157 inverse_dyadic_coefficients.size() * (n_components * n_components);
1161 VectorizedArrayType>::
1162 template run<fe_degree>(n_components,
1165 &inverse_dyadic_coefficients[0][0][0],
1175 &inverse_dyadic_coefficients[0][0][0], unrolled_size),
1187 typename VectorizedArrayType>
1193 VectorizedArrayType>::
1194 transform_from_q_points_to_basis(
const unsigned int n_actual_components,
1195 const VectorizedArrayType *in_array,
1196 VectorizedArrayType *out_array)
const
1198 const auto n_q_points_1d = fe_eval.get_shape_info().data[0].n_q_points_1d;
1200 if (fe_degree > -1 && (fe_degree + 1 == n_q_points_1d))
1203 VectorizedArrayType>::template run<fe_degree,
1204 fe_degree + 1>(n_actual_components,
1219 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1222 , have_interface_matrices(false)
1227 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1234 for (
const unsigned int selected_row : selected_rows)
1235 total_size +=
data->get_vector_partitioner(selected_row)->size();
1241 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1248 for (
const unsigned int selected_column : selected_columns)
1249 total_size +=
data->get_vector_partitioner(selected_column)->size();
1255 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1260 inverse_diagonal_entries.reset();
1265 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1268 const unsigned int col)
const
1271 Assert(inverse_diagonal_entries.get() !=
nullptr &&
1272 inverse_diagonal_entries->m() > 0,
1274 return 1.0 / (*inverse_diagonal_entries)(row, row);
1279 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1282 VectorType &vec)
const
1288 const unsigned int index = selected_rows[i];
1290 .partitioners_are_compatible(
1291 *
data->get_dof_info(index).vector_partitioner))
1295 .partitioners_are_globally_compatible(
1296 *
data->get_dof_info(index).vector_partitioner),
1304 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1309 const std::vector<unsigned int> &given_row_selection,
1310 const std::vector<unsigned int> &given_column_selection)
1314 selected_rows.clear();
1315 selected_columns.clear();
1316 if (given_row_selection.empty())
1317 for (
unsigned int i = 0; i < data_->n_components(); ++i)
1318 selected_rows.push_back(i);
1321 for (
unsigned int i = 0; i < given_row_selection.size(); ++i)
1324 for (
unsigned int j = 0; j < given_row_selection.size(); ++j)
1326 Assert(given_row_selection[j] != given_row_selection[i],
1327 ExcMessage(
"Given row indices must be unique"));
1329 selected_rows.push_back(given_row_selection[i]);
1332 if (given_column_selection.empty())
1333 selected_columns = selected_rows;
1336 for (
unsigned int i = 0; i < given_column_selection.size(); ++i)
1339 for (
unsigned int j = 0; j < given_column_selection.size(); ++j)
1341 Assert(given_column_selection[j] != given_column_selection[i],
1342 ExcMessage(
"Given column indices must be unique"));
1344 selected_columns.push_back(given_column_selection[i]);
1348 edge_constrained_indices.clear();
1349 edge_constrained_indices.resize(selected_rows.size());
1350 edge_constrained_values.clear();
1351 edge_constrained_values.resize(selected_rows.size());
1352 have_interface_matrices =
false;
1357 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1363 const unsigned int level,
1364 const std::vector<unsigned int> &given_row_selection)
1366 std::vector<MGConstrainedDoFs> mg_constrained_dofs_vector(
1367 1, mg_constrained_dofs);
1368 initialize(data_, mg_constrained_dofs_vector,
level, given_row_selection);
1373 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1378 const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
1379 const unsigned int level,
1380 const std::vector<unsigned int> &given_row_selection)
1385 selected_rows.clear();
1386 selected_columns.clear();
1387 if (given_row_selection.empty())
1388 for (
unsigned int i = 0; i < data_->n_components(); ++i)
1389 selected_rows.push_back(i);
1392 for (
unsigned int i = 0; i < given_row_selection.size(); ++i)
1395 for (
unsigned int j = 0; j < given_row_selection.size(); ++j)
1397 Assert(given_row_selection[j] != given_row_selection[i],
1398 ExcMessage(
"Given row indices must be unique"));
1400 selected_rows.push_back(given_row_selection[i]);
1403 selected_columns = selected_rows;
1406 edge_constrained_indices.clear();
1407 edge_constrained_indices.resize(selected_rows.size());
1408 edge_constrained_values.clear();
1409 edge_constrained_values.resize(selected_rows.size());
1413 for (
unsigned int j = 0; j < selected_rows.size(); ++j)
1415 if (data_->n_cell_batches() > 0)
1421 const std::vector<types::global_dof_index> interface_indices =
1422 mg_constrained_dofs[j]
1423 .get_refinement_edge_indices(
level)
1424 .get_index_vector();
1425 edge_constrained_indices[j].clear();
1426 edge_constrained_indices[j].reserve(interface_indices.size());
1427 edge_constrained_values[j].resize(interface_indices.size());
1429 data->get_dof_handler(selected_rows[j]).locally_owned_mg_dofs(
level);
1430 for (
const auto interface_index : interface_indices)
1431 if (locally_owned.
is_element(interface_index))
1432 edge_constrained_indices[j].push_back(
1434 have_interface_matrices |=
1436 static_cast<unsigned int>(edge_constrained_indices[j].
size()),
1437 data->get_vector_partitioner()->get_mpi_communicator()) > 0;
1443 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1446 VectorType &dst)
const
1450 const std::vector<unsigned int> &constrained_dofs =
1451 data->get_constrained_dofs(selected_rows[j]);
1452 for (
const auto constrained_dof : constrained_dofs)
1454 for (
unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1456 edge_constrained_indices[j][i]) = 1.;
1462 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1465 const VectorType &src)
const
1470 vmult_add(dst, src);
1475 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1479 const VectorType &src)
const
1481 mult_add(dst, src,
false);
1486 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1490 const VectorType &src)
const
1492 mult_add(dst, src,
true);
1497 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1500 const VectorType &src,
1501 const bool is_row)
const
1507 const unsigned int mf_component =
1508 is_row ? selected_rows[i] : selected_columns[i];
1511 data->get_dof_info(mf_component).vector_partitioner.get())
1519 data->get_dof_info(mf_component)
1520 .vector_partitioner->locally_owned_size(),
1522 "The vector passed to the vmult() function does not have "
1523 "the correct size for compatibility with MatrixFree."));
1529 this->
data->initialize_dof_vector(
1533 .copy_locally_owned_data_from(copy_vec);
1539 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1543 const VectorType &src)
const
1547 adjust_ghost_range_if_necessary(src,
false);
1548 adjust_ghost_range_if_necessary(dst,
true);
1554 for (
unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1556 edge_constrained_values[j][i] = std::pair<Number, Number>(
1558 edge_constrained_indices[j][i]),
1560 edge_constrained_indices[j][i]));
1562 .local_element(edge_constrained_indices[j][i]) = 0.;
1569 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1573 const VectorType &src,
1579 preprocess_constraints(dst, src);
1581 Tapply_add(dst, src);
1583 apply_add(dst, src);
1584 postprocess_constraints(dst, src);
1589 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1593 const VectorType &src)
const
1597 const std::vector<unsigned int> &constrained_dofs =
1598 data->get_constrained_dofs(selected_rows[j]);
1599 for (
const auto constrained_dof : constrained_dofs)
1608 for (
unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1611 .local_element(edge_constrained_indices[j][i]) =
1612 edge_constrained_values[j][i].first;
1614 edge_constrained_indices[j][i]) =
1615 edge_constrained_values[j][i].second +
1616 edge_constrained_values[j][i].first;
1623 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1627 const VectorType &src)
const
1632 adjust_ghost_range_if_necessary(src,
false);
1633 adjust_ghost_range_if_necessary(dst,
true);
1637 if (!have_interface_matrices)
1643 for (
unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1645 edge_constrained_values[j][i] = std::pair<Number, Number>(
1647 edge_constrained_indices[j][i]),
1649 edge_constrained_indices[j][i]));
1651 .local_element(edge_constrained_indices[j][i]) = 0.;
1654 apply_add(dst, src);
1659 for (
unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1661 for (; c < edge_constrained_indices[j][i]; ++c)
1667 .local_element(edge_constrained_indices[j][i]) =
1668 edge_constrained_values[j][i].first;
1677 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1681 const VectorType &src)
const
1686 adjust_ghost_range_if_necessary(src,
false);
1687 adjust_ghost_range_if_necessary(dst,
true);
1691 if (!have_interface_matrices)
1694 VectorType src_cpy(src);
1698 for (
unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1700 for (; c < edge_constrained_indices[j][i]; ++c)
1708 apply_add(dst, src_cpy);
1711 for (
unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1713 edge_constrained_indices[j][i]) = 0.;
1718 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1722 const VectorType &src)
const
1727 Tvmult_add(dst, src);
1732 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1736 return inverse_diagonal_entries.get() !=
nullptr ?
1737 inverse_diagonal_entries->memory_consumption() :
1743 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1747 VectorizedArrayType>>
1755 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1756 const std::shared_ptr<DiagonalMatrix<VectorType>> &
1760 Assert(inverse_diagonal_entries.get() !=
nullptr &&
1761 inverse_diagonal_entries->m() > 0,
1763 return inverse_diagonal_entries;
1768 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1769 const std::shared_ptr<DiagonalMatrix<VectorType>> &
1772 Assert(diagonal_entries.get() !=
nullptr && diagonal_entries->m() > 0,
1774 return diagonal_entries;
1779 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1783 const VectorType &src)
const
1785 apply_add(dst, src);
1790 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1794 const VectorType &src,
1798 Assert(inverse_diagonal_entries.get() && inverse_diagonal_entries->m() > 0,
1800 inverse_diagonal_entries->vmult(dst, src);
1808 template <
typename OperatorType>
1811 , mf_base_operator(nullptr)
1816 template <
typename OperatorType>
1820 mf_base_operator =
nullptr;
1825 template <
typename OperatorType>
1829 mf_base_operator = &operator_in;
1834 template <
typename OperatorType>
1835 template <
typename VectorType>
1838 const VectorType &src)
const
1842 std::is_same_v<typename VectorType::value_type, value_type>,
1843 "The vector type must be based on the same value type as this "
1849 mf_base_operator->vmult_interface_down(dst, src);
1854 template <
typename OperatorType>
1855 template <
typename VectorType>
1858 const VectorType &src)
const
1862 std::is_same_v<typename VectorType::value_type, value_type>,
1863 "The vector type must be based on the same value type as this "
1869 mf_base_operator->vmult_interface_up(dst, src);
1874 template <
typename OperatorType>
1875 template <
typename VectorType>
1878 VectorType &vec)
const
1882 mf_base_operator->initialize_dof_vector(vec);
1893 typename VectorType,
1894 typename VectorizedArrayType>
1900 VectorizedArrayType>::MassOperator()
1901 :
Base<dim, VectorType, VectorizedArrayType>()
1906 "This class only supports the non-blocked vector variant of the Base "
1907 "operator because only a single FEEvaluation object is used in the "
1908 "apply function."));
1917 typename VectorType,
1918 typename VectorizedArrayType>
1925 VectorizedArrayType>::compute_diagonal()
1929 Assert(this->selected_rows == this->selected_columns,
1930 ExcMessage(
"This function is only implemented for square (not "
1931 "rectangular) operators."));
1933 this->inverse_diagonal_entries =
1934 std::make_shared<DiagonalMatrix<VectorType>>();
1935 this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1936 VectorType &inverse_diagonal_vector =
1937 this->inverse_diagonal_entries->get_vector();
1938 VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1939 this->initialize_dof_vector(inverse_diagonal_vector);
1940 this->initialize_dof_vector(diagonal_vector);
1944 auto diagonal_evaluation = [](
auto &integrator) {
1946 for (
unsigned int q = 0; q < integrator.n_q_points; ++q)
1947 integrator.submit_value(integrator.get_value(q), q);
1958 VectorizedArrayType> &)>
1959 diagonal_evaluation_f(diagonal_evaluation);
1962 for (
unsigned int block_n = 0; block_n < this->selected_rows.size();
1967 diagonal_evaluation_f,
1968 this->selected_rows[block_n]);
1972 this->set_constrained_entries_to_one(diagonal_vector);
1974 inverse_diagonal_vector = diagonal_vector;
1976 for (
unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
1984 Assert(diagonal_vector.local_element(i) > Number(0),
1987 inverse_diagonal_vector.local_element(i) =
1988 1. / inverse_diagonal_vector.local_element(i);
2000 typename VectorType,
2001 typename VectorizedArrayType>
2008 VectorizedArrayType>::compute_lumped_diagonal()
2014 Assert(this->selected_rows == this->selected_columns,
2015 ExcMessage(
"This function is only implemented for square (not "
2016 "rectangular) operators."));
2018 inverse_lumped_diagonal_entries =
2019 std::make_shared<DiagonalMatrix<VectorType>>();
2020 lumped_diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
2021 VectorType &inverse_lumped_diagonal_vector =
2022 inverse_lumped_diagonal_entries->get_vector();
2023 VectorType &lumped_diagonal_vector = lumped_diagonal_entries->get_vector();
2024 this->initialize_dof_vector(inverse_lumped_diagonal_vector);
2025 this->initialize_dof_vector(lumped_diagonal_vector);
2028 inverse_lumped_diagonal_vector = Number(1.);
2029 apply_add(lumped_diagonal_vector, inverse_lumped_diagonal_vector);
2030 this->set_constrained_entries_to_one(lumped_diagonal_vector);
2033 inverse_lumped_diagonal_vector.locally_owned_size();
2040 if (lumped_diagonal_vector.local_element(i) == Number(0.))
2041 inverse_lumped_diagonal_vector.local_element(i) = Number(1.);
2043 inverse_lumped_diagonal_vector.local_element(i) =
2044 Number(1.) / lumped_diagonal_vector.local_element(i);
2047 inverse_lumped_diagonal_vector.update_ghost_values();
2048 lumped_diagonal_vector.update_ghost_values();
2057 typename VectorType,
2058 typename VectorizedArrayType>
2059 const std::shared_ptr<DiagonalMatrix<VectorType>> &
2065 VectorizedArrayType>::get_matrix_lumped_diagonal_inverse()
const
2067 Assert(inverse_lumped_diagonal_entries.get() !=
nullptr &&
2068 inverse_lumped_diagonal_entries->m() > 0,
2070 return inverse_lumped_diagonal_entries;
2079 typename VectorType,
2080 typename VectorizedArrayType>
2081 const std::shared_ptr<DiagonalMatrix<VectorType>> &
2087 VectorizedArrayType>::get_matrix_lumped_diagonal()
const
2089 Assert(lumped_diagonal_entries.get() !=
nullptr &&
2090 lumped_diagonal_entries->m() > 0,
2092 return lumped_diagonal_entries;
2101 typename VectorType,
2102 typename VectorizedArrayType>
2109 VectorizedArrayType>::apply_add(VectorType &dst,
2110 const VectorType &src)
const
2122 typename VectorType,
2123 typename VectorizedArrayType>
2130 VectorizedArrayType>::
2135 VectorizedArrayType> &
data,
2137 const VectorType &src,
2138 const std::pair<unsigned int, unsigned int> &cell_range)
const
2147 VectorizedArrayType>
2148 phi(
data, this->selected_rows[0]);
2149 for (
unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2152 phi.read_dof_values(src);
2154 for (
unsigned int q = 0; q < phi.n_q_points; ++q)
2155 phi.submit_value(phi.get_value(q), q);
2157 phi.distribute_local_to_global(dst);
2168 typename VectorType,
2169 typename VectorizedArrayType>
2175 VectorizedArrayType>::LaplaceOperator()
2176 :
Base<dim, VectorType, VectorizedArrayType>()
2185 typename VectorType,
2186 typename VectorizedArrayType>
2193 VectorizedArrayType>::clear()
2196 scalar_coefficient.reset();
2205 typename VectorType,
2206 typename VectorizedArrayType>
2213 VectorizedArrayType>::
2217 scalar_coefficient = scalar_coefficient_;
2226 typename VectorType,
2227 typename VectorizedArrayType>
2228 std::shared_ptr<Table<2, VectorizedArrayType>>
2234 VectorizedArrayType>::get_coefficient()
2237 return scalar_coefficient;
2246 typename VectorType,
2247 typename VectorizedArrayType>
2254 VectorizedArrayType>::compute_diagonal()
2261 this->inverse_diagonal_entries =
2262 std::make_shared<DiagonalMatrix<VectorType>>();
2263 this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
2264 VectorType &inverse_diagonal_vector =
2265 this->inverse_diagonal_entries->get_vector();
2266 VectorType &diagonal_vector = this->diagonal_entries->get_vector();
2267 this->initialize_dof_vector(inverse_diagonal_vector);
2268 this->initialize_dof_vector(diagonal_vector);
2274 this->set_constrained_entries_to_one(diagonal_vector);
2276 inverse_diagonal_vector = diagonal_vector;
2278 for (
unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
2280 if (
std::abs(inverse_diagonal_vector.local_element(i)) >
2281 std::sqrt(std::numeric_limits<Number>::epsilon()))
2282 inverse_diagonal_vector.local_element(i) =
2283 1. / inverse_diagonal_vector.local_element(i);
2285 inverse_diagonal_vector.local_element(i) = 1.;
2296 typename VectorType,
2297 typename VectorizedArrayType>
2304 VectorizedArrayType>::apply_add(VectorType &dst,
2305 const VectorType &src)
const
2311 namespace Implementation
2313 template <
typename VectorizedArrayType>
2317 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2331 typename VectorType,
2332 typename VectorizedArrayType>
2333 template <
int n_components_compute>
2335 LaplaceOperator<dim,
2340 VectorizedArrayType>::
2341 do_operation_on_cell(
2346 n_components_compute,
2348 VectorizedArrayType> &phi,
2349 const unsigned int cell)
const
2352 if (scalar_coefficient.get())
2354 Assert(scalar_coefficient->size(1) == 1 ||
2355 scalar_coefficient->size(1) == phi.n_q_points,
2356 ExcMessage(
"The number of columns in the coefficient table must "
2357 "be either 1 or the number of quadrature points " +
2358 std::to_string(phi.n_q_points) +
2359 ", but the given value was " +
2360 std::to_string(scalar_coefficient->size(1))));
2361 if (scalar_coefficient->size(1) == phi.n_q_points)
2362 for (
unsigned int q = 0; q < phi.n_q_points; ++q)
2365 (*scalar_coefficient)(cell, q)),
2366 ExcMessage(
"Coefficient must be non-negative"));
2367 phi.submit_gradient((*scalar_coefficient)(cell, q) *
2368 phi.get_gradient(q),
2374 ExcMessage(
"Coefficient must be non-negative"));
2375 const VectorizedArrayType coefficient =
2376 (*scalar_coefficient)(cell, 0);
2377 for (
unsigned int q = 0; q < phi.n_q_points; ++q)
2378 phi.submit_gradient(coefficient * phi.get_gradient(q), q);
2383 for (
unsigned int q = 0; q < phi.n_q_points; ++q)
2385 phi.submit_gradient(phi.get_gradient(q), q);
2397 typename VectorType,
2398 typename VectorizedArrayType>
2405 VectorizedArrayType>::
2410 VectorizedArrayType> &
data,
2412 const VectorType &src,
2413 const std::pair<unsigned int, unsigned int> &cell_range)
const
2422 VectorizedArrayType>
2423 phi(
data, this->selected_rows[0]);
2424 for (
unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2427 phi.read_dof_values(src);
2428 do_operation_on_cell(phi, cell);
2429 phi.distribute_local_to_global(dst);
2438 typename VectorType,
2439 typename VectorizedArrayType>
2446 VectorizedArrayType>::
2447 local_diagonal_cell(
2451 VectorizedArrayType> &
data,
2454 const std::pair<unsigned int, unsigned int> &cell_range)
const
2460 eval(
data, this->selected_rows[0]);
2466 VectorizedArrayType>
2467 eval_vector(
data, this->selected_rows[0]);
2468 for (
unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2471 eval_vector.reinit(cell);
2483 do_operation_on_cell(eval, cell);
2488 for (
unsigned int c = 0; c < n_components; ++c)
2490 .begin_dof_values()[i + c * eval_vector.dofs_per_component] =
2493 eval_vector.distribute_local_to_global(dst);
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)