17 #ifndef dealii_matrix_free_operators_h
18 #define dealii_matrix_free_operators_h
48 template <
typename VectorType>
49 std::enable_if_t<IsBlockVector<VectorType>::value,
unsigned int>
52 return vector.n_blocks();
55 template <
typename VectorType>
56 std::enable_if_t<!IsBlockVector<VectorType>::value,
unsigned int>
62 template <
typename VectorType>
63 std::enable_if_t<IsBlockVector<VectorType>::value,
64 typename VectorType::BlockType &>
65 subblock(VectorType &vector,
unsigned int block_no)
68 return vector.block(block_no);
71 template <
typename VectorType>
72 std::enable_if_t<IsBlockVector<VectorType>::value,
73 const typename VectorType::BlockType &>
74 subblock(
const VectorType &vector,
unsigned int block_no)
77 return vector.block(block_no);
80 template <
typename VectorType>
81 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType &>
87 template <
typename VectorType>
88 std::enable_if_t<!IsBlockVector<VectorType>::value,
const VectorType &>
89 subblock(
const VectorType &vector,
unsigned int)
94 template <
typename VectorType>
95 std::enable_if_t<IsBlockVector<VectorType>::value,
void>
98 vector.collect_sizes();
101 template <
typename VectorType>
102 std::enable_if_t<!IsBlockVector<VectorType>::value,
void>
186 typename VectorizedArrayType =
237 const std::vector<unsigned int> &selected_row_blocks =
238 std::vector<unsigned int>(),
239 const std::vector<unsigned int> &selected_column_blocks =
240 std::vector<unsigned int>());
259 const unsigned int level,
260 const std::vector<unsigned int> &selected_row_blocks =
261 std::vector<unsigned int>());
280 const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
281 const unsigned int level,
282 const std::vector<unsigned int> & selected_row_blocks =
283 std::vector<unsigned int>());
313 vmult(VectorType &dst,
const VectorType &src)
const;
319 Tvmult(VectorType &dst,
const VectorType &src)
const;
325 vmult_add(VectorType &dst,
const VectorType &src)
const;
338 el(
const unsigned int row,
const unsigned int col)
const;
371 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
377 const std::shared_ptr<DiagonalMatrix<VectorType>> &
383 const std::shared_ptr<DiagonalMatrix<VectorType>> &
393 const VectorType &src,
422 apply_add(VectorType &dst,
const VectorType &src)
const = 0;
435 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
471 mutable std::vector<std::vector<std::pair<value_type, value_type>>>
486 const VectorType &src,
498 const bool is_row)
const;
537 template <
typename OperatorType>
571 template <
typename VectorType>
573 vmult(VectorType &dst,
const VectorType &src)
const;
578 template <
typename VectorType>
580 Tvmult(VectorType &dst,
const VectorType &src)
const;
585 template <
typename VectorType>
618 int n_components = 1,
624 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
625 "Type of Number and of VectorizedArrayType do not match.");
637 VectorizedArrayType> &
fe_eval);
649 const unsigned int n_actual_components,
650 const VectorizedArrayType * in_array,
651 VectorizedArrayType * out_array)
const;
665 apply(
const VectorizedArrayType *in_array,
666 VectorizedArrayType * out_array)
const;
703 const VectorizedArrayType *in_array,
704 VectorizedArrayType *out_array)
const;
737 int n_q_points_1d = fe_degree + 1,
738 int n_components = 1,
740 typename VectorizedArrayType =
789 const std::shared_ptr<DiagonalMatrix<VectorType>> &
795 const std::shared_ptr<DiagonalMatrix<VectorType>> &
805 apply_add(VectorType &dst,
const VectorType &src)
const override;
814 const VectorType & src,
815 const std::pair<unsigned int, unsigned int> & cell_range)
const;
844 int n_q_points_1d = fe_degree + 1,
845 int n_components = 1,
847 typename VectorizedArrayType =
945 std::shared_ptr<Table<2, VectorizedArrayType>>
955 apply_add(VectorType &dst,
const VectorType &src)
const override;
964 const VectorType & src,
965 const std::pair<unsigned int, unsigned int> & cell_range)
const;
975 const std::pair<unsigned int, unsigned int> &cell_range)
const;
980 template <
int n_components_compute>
985 n_components_compute,
987 VectorizedArrayType> &phi,
988 const unsigned int cell)
const;
1004 typename VectorizedArrayType>
1009 VectorizedArrayType>::
1010 CellwiseInverseMassMatrix(
1015 VectorizedArrayType> &fe_eval)
1025 typename VectorizedArrayType>
1031 VectorizedArrayType>::
1032 fill_inverse_JxW_values(
1035 constexpr
unsigned int dofs_per_component_on_cell =
1038 inverse_jxw.
size() % dofs_per_component_on_cell == 0,
1040 "Expected diagonal to be a multiple of scalar dof per cells"));
1043 for (
unsigned int q = 0; q < dofs_per_component_on_cell; ++q)
1044 inverse_jxw[q] = 1. / fe_eval.JxW(q);
1046 for (
unsigned int q = dofs_per_component_on_cell; q < inverse_jxw.
size();)
1047 for (
unsigned int i = 0; i < dofs_per_component_on_cell; ++i, ++q)
1048 inverse_jxw[q] = inverse_jxw[i];
1057 typename VectorizedArrayType>
1064 VectorizedArrayType>
::apply(
const VectorizedArrayType *in_array,
1065 VectorizedArrayType * out_array)
const
1069 template run<fe_degree>(n_components, fe_eval, in_array, out_array);
1072 n_components, fe_eval, in_array, out_array);
1081 typename VectorizedArrayType>
1089 const unsigned int n_actual_components,
1090 const VectorizedArrayType * in_array,
1091 VectorizedArrayType * out_array)
const
1093 const unsigned int given_degree =
1094 fe_eval.get_shape_info().data[0].fe_degree;
1097 VectorizedArrayType>::
1098 template run<fe_degree>(
1099 n_actual_components,
1100 fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
1101 inverse_coefficients,
1106 n_actual_components,
1108 fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
1109 inverse_coefficients,
1120 typename VectorizedArrayType>
1126 VectorizedArrayType>::
1127 transform_from_q_points_to_basis(
const unsigned int n_actual_components,
1128 const VectorizedArrayType *in_array,
1129 VectorizedArrayType * out_array)
const
1131 const auto n_q_points_1d = fe_eval.get_shape_info().data[0].n_q_points_1d;
1133 if (fe_degree > -1 && (fe_degree + 1 == n_q_points_1d))
1136 VectorizedArrayType>::template
run<fe_degree,
1137 fe_degree + 1>(n_actual_components,
1152 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1155 , have_interface_matrices(false)
1160 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1167 for (
const unsigned int selected_row : selected_rows)
1168 total_size += data->get_vector_partitioner(selected_row)->size();
1174 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1181 for (
const unsigned int selected_column : selected_columns)
1182 total_size += data->get_vector_partitioner(selected_column)->size();
1188 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1193 inverse_diagonal_entries.reset();
1198 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1201 const unsigned int col)
const
1205 Assert(inverse_diagonal_entries.get() !=
nullptr &&
1206 inverse_diagonal_entries->m() > 0,
1208 return 1.0 / (*inverse_diagonal_entries)(row, row);
1213 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1216 VectorType &vec)
const
1222 const unsigned int index = selected_rows[i];
1224 .partitioners_are_compatible(
1225 *data->get_dof_info(index).vector_partitioner))
1229 .partitioners_are_globally_compatible(
1230 *data->get_dof_info(index).vector_partitioner),
1238 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1243 const std::vector<unsigned int> &given_row_selection,
1244 const std::vector<unsigned int> &given_column_selection)
1248 selected_rows.clear();
1249 selected_columns.clear();
1250 if (given_row_selection.empty())
1251 for (
unsigned int i = 0; i < data_->n_components(); ++i)
1252 selected_rows.push_back(i);
1255 for (
unsigned int i = 0; i < given_row_selection.size(); ++i)
1258 for (
unsigned int j = 0; j < given_row_selection.size(); ++j)
1260 Assert(given_row_selection[j] != given_row_selection[i],
1261 ExcMessage(
"Given row indices must be unique"));
1263 selected_rows.push_back(given_row_selection[i]);
1266 if (given_column_selection.size() == 0)
1267 selected_columns = selected_rows;
1270 for (
unsigned int i = 0; i < given_column_selection.size(); ++i)
1273 for (
unsigned int j = 0; j < given_column_selection.size(); ++j)
1275 Assert(given_column_selection[j] != given_column_selection[i],
1276 ExcMessage(
"Given column indices must be unique"));
1278 selected_columns.push_back(given_column_selection[i]);
1282 edge_constrained_indices.clear();
1283 edge_constrained_indices.resize(selected_rows.size());
1284 edge_constrained_values.clear();
1285 edge_constrained_values.resize(selected_rows.size());
1286 have_interface_matrices =
false;
1291 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1297 const unsigned int level,
1298 const std::vector<unsigned int> &given_row_selection)
1300 std::vector<MGConstrainedDoFs> mg_constrained_dofs_vector(
1301 1, mg_constrained_dofs);
1302 initialize(data_, mg_constrained_dofs_vector,
level, given_row_selection);
1307 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1312 const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
1313 const unsigned int level,
1314 const std::vector<unsigned int> & given_row_selection)
1319 selected_rows.clear();
1320 selected_columns.clear();
1321 if (given_row_selection.empty())
1322 for (
unsigned int i = 0; i < data_->n_components(); ++i)
1323 selected_rows.push_back(i);
1326 for (
unsigned int i = 0; i < given_row_selection.size(); ++i)
1329 for (
unsigned int j = 0; j < given_row_selection.size(); ++j)
1331 Assert(given_row_selection[j] != given_row_selection[i],
1332 ExcMessage(
"Given row indices must be unique"));
1334 selected_rows.push_back(given_row_selection[i]);
1337 selected_columns = selected_rows;
1340 edge_constrained_indices.clear();
1341 edge_constrained_indices.resize(selected_rows.size());
1342 edge_constrained_values.clear();
1343 edge_constrained_values.resize(selected_rows.size());
1347 for (
unsigned int j = 0; j < selected_rows.size(); ++j)
1349 if (data_->n_cell_batches() > 0)
1355 const std::vector<types::global_dof_index> interface_indices =
1356 mg_constrained_dofs[j]
1357 .get_refinement_edge_indices(
level)
1358 .get_index_vector();
1359 edge_constrained_indices[j].clear();
1360 edge_constrained_indices[j].reserve(interface_indices.size());
1361 edge_constrained_values[j].resize(interface_indices.size());
1363 data->get_dof_handler(selected_rows[j]).locally_owned_mg_dofs(
level);
1364 for (
const auto interface_index : interface_indices)
1365 if (locally_owned.
is_element(interface_index))
1366 edge_constrained_indices[j].push_back(
1368 have_interface_matrices |=
1370 static_cast<unsigned int>(edge_constrained_indices[j].size()),
1371 data->get_vector_partitioner()->get_mpi_communicator()) > 0;
1377 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1380 VectorType &dst)
const
1384 const std::vector<unsigned int> &constrained_dofs =
1385 data->get_constrained_dofs(selected_rows[j]);
1386 for (
const auto constrained_dof : constrained_dofs)
1388 for (
unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1390 edge_constrained_indices[j][i]) = 1.;
1396 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1399 const VectorType &src)
const
1404 vmult_add(dst, src);
1409 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1413 const VectorType &src)
const
1415 mult_add(dst, src,
false);
1420 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1424 const VectorType &src)
const
1426 mult_add(dst, src,
true);
1431 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1434 const VectorType &src,
1435 const bool is_row)
const
1441 const unsigned int mf_component =
1442 is_row ? selected_rows[i] : selected_columns[i];
1445 data->get_dof_info(mf_component).vector_partitioner.get())
1452 ->locally_owned_size() ==
1453 data->get_dof_info(mf_component)
1454 .vector_partitioner->locally_owned_size(),
1456 "The vector passed to the vmult() function does not have "
1457 "the correct size for compatibility with MatrixFree."));
1463 this->data->initialize_dof_vector(
1467 .copy_locally_owned_data_from(copy_vec);
1473 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1477 const VectorType &src)
const
1481 adjust_ghost_range_if_necessary(src,
false);
1482 adjust_ghost_range_if_necessary(dst,
true);
1488 for (
unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1490 edge_constrained_values[j][i] = std::pair<Number, Number>(
1492 edge_constrained_indices[j][i]),
1494 edge_constrained_indices[j][i]));
1496 .local_element(edge_constrained_indices[j][i]) = 0.;
1503 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1507 const VectorType &src,
1513 preprocess_constraints(dst, src);
1515 Tapply_add(dst, src);
1517 apply_add(dst, src);
1518 postprocess_constraints(dst, src);
1523 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1527 const VectorType &src)
const
1531 const std::vector<unsigned int> &constrained_dofs =
1532 data->get_constrained_dofs(selected_rows[j]);
1533 for (
const auto constrained_dof : constrained_dofs)
1542 for (
unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1545 .local_element(edge_constrained_indices[j][i]) =
1546 edge_constrained_values[j][i].first;
1548 edge_constrained_indices[j][i]) =
1549 edge_constrained_values[j][i].second +
1550 edge_constrained_values[j][i].first;
1557 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1561 const VectorType &src)
const
1566 adjust_ghost_range_if_necessary(src,
false);
1567 adjust_ghost_range_if_necessary(dst,
true);
1571 if (!have_interface_matrices)
1577 for (
unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1579 edge_constrained_values[j][i] = std::pair<Number, Number>(
1581 edge_constrained_indices[j][i]),
1583 edge_constrained_indices[j][i]));
1585 .local_element(edge_constrained_indices[j][i]) = 0.;
1588 apply_add(dst, src);
1593 for (
unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1595 for (; c < edge_constrained_indices[j][i]; ++c)
1601 .local_element(edge_constrained_indices[j][i]) =
1602 edge_constrained_values[j][i].first;
1611 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1615 const VectorType &src)
const
1620 adjust_ghost_range_if_necessary(src,
false);
1621 adjust_ghost_range_if_necessary(dst,
true);
1625 if (!have_interface_matrices)
1628 VectorType src_cpy(src);
1632 for (
unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1634 for (; c < edge_constrained_indices[j][i]; ++c)
1642 apply_add(dst, src_cpy);
1645 for (
unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1647 edge_constrained_indices[j][i]) = 0.;
1652 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1656 const VectorType &src)
const
1661 Tvmult_add(dst, src);
1666 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1670 return inverse_diagonal_entries.get() !=
nullptr ?
1671 inverse_diagonal_entries->memory_consumption() :
1677 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1681 VectorizedArrayType>>
1689 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1690 const std::shared_ptr<DiagonalMatrix<VectorType>> &
1694 Assert(inverse_diagonal_entries.get() !=
nullptr &&
1695 inverse_diagonal_entries->m() > 0,
1697 return inverse_diagonal_entries;
1702 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1703 const std::shared_ptr<DiagonalMatrix<VectorType>> &
1706 Assert(diagonal_entries.get() !=
nullptr && diagonal_entries->m() > 0,
1708 return diagonal_entries;
1713 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1717 const VectorType &src)
const
1719 apply_add(dst, src);
1724 template <
int dim,
typename VectorType,
typename VectorizedArrayType>
1728 const VectorType & src,
1732 Assert(inverse_diagonal_entries.get() && inverse_diagonal_entries->m() > 0,
1734 inverse_diagonal_entries->vmult(dst, src);
1742 template <
typename OperatorType>
1745 , mf_base_operator(nullptr)
1750 template <
typename OperatorType>
1754 mf_base_operator =
nullptr;
1759 template <
typename OperatorType>
1763 mf_base_operator = &operator_in;
1768 template <
typename OperatorType>
1769 template <
typename VectorType>
1772 const VectorType &src)
const
1774 #ifndef DEAL_II_MSVC
1776 std::is_same<typename VectorType::value_type, value_type>::value,
1777 "The vector type must be based on the same value type as this "
1783 mf_base_operator->vmult_interface_down(dst, src);
1788 template <
typename OperatorType>
1789 template <
typename VectorType>
1792 const VectorType &src)
const
1794 #ifndef DEAL_II_MSVC
1796 std::is_same<typename VectorType::value_type, value_type>::value,
1797 "The vector type must be based on the same value type as this "
1803 mf_base_operator->vmult_interface_up(dst, src);
1808 template <
typename OperatorType>
1809 template <
typename VectorType>
1812 VectorType &vec)
const
1816 mf_base_operator->initialize_dof_vector(vec);
1827 typename VectorType,
1828 typename VectorizedArrayType>
1835 :
Base<dim, VectorType, VectorizedArrayType>()
1840 "This class only supports the non-blocked vector variant of the Base "
1841 "operator because only a single FEEvaluation object is used in the "
1842 "apply function."));
1851 typename VectorType,
1852 typename VectorizedArrayType>
1865 Assert(this->selected_rows == this->selected_columns,
1866 ExcMessage(
"This function is only implemented for square (not "
1867 "rectangular) operators."));
1869 this->inverse_diagonal_entries =
1870 std::make_shared<DiagonalMatrix<VectorType>>();
1871 this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1872 VectorType &inverse_diagonal_vector =
1873 this->inverse_diagonal_entries->get_vector();
1874 VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1875 this->initialize_dof_vector(inverse_diagonal_vector);
1876 this->initialize_dof_vector(diagonal_vector);
1880 auto diagonal_evaluation = [](
auto &integrator) {
1882 for (
unsigned int q = 0; q < integrator.n_q_points; ++q)
1883 integrator.submit_value(integrator.get_value(q), q);
1894 VectorizedArrayType> &)>
1895 diagonal_evaluation_f(diagonal_evaluation);
1898 for (
unsigned int block_n = 0; block_n < this->selected_rows.size();
1903 diagonal_evaluation_f,
1904 this->selected_rows[block_n]);
1908 this->set_constrained_entries_to_one(diagonal_vector);
1910 inverse_diagonal_vector = diagonal_vector;
1912 for (
unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
1915 Assert(diagonal_vector.local_element(i) > Number(0),
1917 inverse_diagonal_vector.local_element(i) =
1918 1. / inverse_diagonal_vector.local_element(i);
1930 typename VectorType,
1931 typename VectorizedArrayType>
1938 VectorizedArrayType>::compute_lumped_diagonal()
1944 Assert(this->selected_rows == this->selected_columns,
1945 ExcMessage(
"This function is only implemented for square (not "
1946 "rectangular) operators."));
1948 inverse_lumped_diagonal_entries =
1949 std::make_shared<DiagonalMatrix<VectorType>>();
1950 lumped_diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1951 VectorType &inverse_lumped_diagonal_vector =
1952 inverse_lumped_diagonal_entries->get_vector();
1953 VectorType &lumped_diagonal_vector = lumped_diagonal_entries->get_vector();
1954 this->initialize_dof_vector(inverse_lumped_diagonal_vector);
1955 this->initialize_dof_vector(lumped_diagonal_vector);
1958 inverse_lumped_diagonal_vector = Number(1.);
1959 apply_add(lumped_diagonal_vector, inverse_lumped_diagonal_vector);
1960 this->set_constrained_entries_to_one(lumped_diagonal_vector);
1963 inverse_lumped_diagonal_vector.locally_owned_size();
1968 for (
size_type i = 0; i < locally_owned_size; ++i)
1970 if (lumped_diagonal_vector.local_element(i) == Number(0.))
1971 inverse_lumped_diagonal_vector.local_element(i) = Number(1.);
1973 inverse_lumped_diagonal_vector.local_element(i) =
1974 Number(1.) / lumped_diagonal_vector.local_element(i);
1977 inverse_lumped_diagonal_vector.update_ghost_values();
1978 lumped_diagonal_vector.update_ghost_values();
1987 typename VectorType,
1988 typename VectorizedArrayType>
1989 const std::shared_ptr<DiagonalMatrix<VectorType>> &
1995 VectorizedArrayType>::get_matrix_lumped_diagonal_inverse()
const
1997 Assert(inverse_lumped_diagonal_entries.get() !=
nullptr &&
1998 inverse_lumped_diagonal_entries->m() > 0,
2000 return inverse_lumped_diagonal_entries;
2009 typename VectorType,
2010 typename VectorizedArrayType>
2011 const std::shared_ptr<DiagonalMatrix<VectorType>> &
2017 VectorizedArrayType>::get_matrix_lumped_diagonal()
const
2019 Assert(lumped_diagonal_entries.get() !=
nullptr &&
2020 lumped_diagonal_entries->m() > 0,
2022 return lumped_diagonal_entries;
2031 typename VectorType,
2032 typename VectorizedArrayType>
2039 VectorizedArrayType>::apply_add(VectorType & dst,
2040 const VectorType &src)
const
2052 typename VectorType,
2053 typename VectorizedArrayType>
2060 VectorizedArrayType>::
2065 VectorizedArrayType> & data,
2067 const VectorType & src,
2068 const std::pair<unsigned int, unsigned int> &cell_range)
const
2077 VectorizedArrayType>
2078 phi(data, this->selected_rows[0]);
2079 for (
unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2082 phi.read_dof_values(src);
2084 for (
unsigned int q = 0; q < phi.n_q_points; ++q)
2085 phi.submit_value(phi.get_value(q), q);
2087 phi.distribute_local_to_global(dst);
2098 typename VectorType,
2099 typename VectorizedArrayType>
2106 :
Base<dim, VectorType, VectorizedArrayType>()
2115 typename VectorType,
2116 typename VectorizedArrayType>
2123 VectorizedArrayType>::clear()
2126 scalar_coefficient.reset();
2135 typename VectorType,
2136 typename VectorizedArrayType>
2143 VectorizedArrayType>::
2147 scalar_coefficient = scalar_coefficient_;
2156 typename VectorType,
2157 typename VectorizedArrayType>
2158 std::shared_ptr<Table<2, VectorizedArrayType>>
2164 VectorizedArrayType>::get_coefficient()
2167 return scalar_coefficient;
2176 typename VectorType,
2177 typename VectorizedArrayType>
2191 this->inverse_diagonal_entries =
2192 std::make_shared<DiagonalMatrix<VectorType>>();
2193 this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
2194 VectorType &inverse_diagonal_vector =
2195 this->inverse_diagonal_entries->get_vector();
2196 VectorType &diagonal_vector = this->diagonal_entries->get_vector();
2197 this->initialize_dof_vector(inverse_diagonal_vector);
2198 this->initialize_dof_vector(diagonal_vector);
2204 this->set_constrained_entries_to_one(diagonal_vector);
2206 inverse_diagonal_vector = diagonal_vector;
2208 for (
unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
2210 if (std::abs(inverse_diagonal_vector.local_element(i)) >
2212 inverse_diagonal_vector.local_element(i) =
2213 1. / inverse_diagonal_vector.local_element(i);
2215 inverse_diagonal_vector.local_element(i) = 1.;
2226 typename VectorType,
2227 typename VectorizedArrayType>
2234 VectorizedArrayType>::apply_add(VectorType & dst,
2235 const VectorType &src)
const
2241 namespace Implementation
2243 template <
typename VectorizedArrayType>
2247 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2261 typename VectorType,
2262 typename VectorizedArrayType>
2263 template <
int n_components_compute>
2265 LaplaceOperator<dim,
2270 VectorizedArrayType>::
2271 do_operation_on_cell(
2276 n_components_compute,
2278 VectorizedArrayType> &phi,
2279 const unsigned int cell)
const
2282 if (scalar_coefficient.get())
2284 Assert(scalar_coefficient->size(1) == 1 ||
2285 scalar_coefficient->size(1) == phi.n_q_points,
2286 ExcMessage(
"The number of columns in the coefficient table must "
2287 "be either 1 or the number of quadrature points " +
2289 ", but the given value was " +
2291 if (scalar_coefficient->size(1) == phi.n_q_points)
2292 for (
unsigned int q = 0; q < phi.n_q_points; ++q)
2295 (*scalar_coefficient)(cell, q)),
2296 ExcMessage(
"Coefficient must be non-negative"));
2297 phi.submit_gradient((*scalar_coefficient)(cell, q) *
2298 phi.get_gradient(q),
2304 ExcMessage(
"Coefficient must be non-negative"));
2305 const VectorizedArrayType coefficient =
2306 (*scalar_coefficient)(cell, 0);
2307 for (
unsigned int q = 0; q < phi.n_q_points; ++q)
2308 phi.submit_gradient(coefficient * phi.get_gradient(q), q);
2313 for (
unsigned int q = 0; q < phi.n_q_points; ++q)
2315 phi.submit_gradient(phi.get_gradient(q), q);
2327 typename VectorType,
2328 typename VectorizedArrayType>
2335 VectorizedArrayType>::
2340 VectorizedArrayType> & data,
2342 const VectorType & src,
2343 const std::pair<unsigned int, unsigned int> &cell_range)
const
2352 VectorizedArrayType>
2353 phi(data, this->selected_rows[0]);
2354 for (
unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2357 phi.read_dof_values(src);
2358 do_operation_on_cell(phi, cell);
2359 phi.distribute_local_to_global(dst);
2368 typename VectorType,
2369 typename VectorizedArrayType>
2376 VectorizedArrayType>::
2377 local_diagonal_cell(
2381 VectorizedArrayType> &data,
2384 const std::pair<unsigned int, unsigned int> &cell_range)
const
2390 eval(data, this->selected_rows[0]);
2396 VectorizedArrayType>
2397 eval_vector(data, this->selected_rows[0]);
2398 for (
unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2401 eval_vector.reinit(cell);
2413 do_operation_on_cell(eval, cell);
2418 for (
unsigned int c = 0; c < n_components; ++c)
2420 .begin_dof_values()[i + c * eval_vector.dofs_per_component] =
2423 eval_vector.distribute_local_to_global(dst);
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
bool have_interface_matrices
std::vector< unsigned int > selected_columns
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > get_matrix_free() const
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal_inverse() const
void initialize_dof_vector(VectorType &vec) const
std::shared_ptr< DiagonalMatrix< VectorType > > diagonal_entries
virtual void Tapply_add(VectorType &dst, const VectorType &src) 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 >())
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
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 >())
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 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 >())
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
CellwiseInverseMassMatrix(const FEEvaluationBase< dim, n_components, Number, false, VectorizedArrayType > &fe_eval)
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 set_coefficient(const std::shared_ptr< Table< 2, VectorizedArrayType >> &scalar_coefficient)
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 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
SmartPointer< const OperatorType > mf_base_operator
void vmult(VectorType &dst, const VectorType &src) const
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
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
types::global_dof_index size_type
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
std::enable_if_t< IsBlockVector< VectorType >::value, void > collect_sizes(VectorType &vector)
std::enable_if_t< IsBlockVector< VectorType >::value, typename VectorType::BlockType & > subblock(VectorType &vector, unsigned int block_no)
bool non_negative(const VectorizedArrayType &n)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
T max(const T &t, const MPI_Comm &mpi_communicator)
constexpr T pow(const T base, const int iexp)
void run(const Iterator &begin, const typename identity< Iterator >::type &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
static const unsigned int invalid_unsigned_int
auto apply(F &&fn, Tuple &&t) -> decltype(apply_impl(std::forward< F >(fn), std::forward< Tuple >(t), std::make_index_sequence< std::tuple_size< typename std::remove_reference< Tuple >::type >::value >()))
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)