deal.II version GIT relicensing-6809-ge913b9bb34 2026-09-25 17:20:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
operators.h
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2014 - 2026 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13
14#ifndef dealii_matrix_free_operators_h
15#define dealii_matrix_free_operators_h
16
17
18#include <deal.II/base/config.h>
19
20#include <deal.II/base/enable_observer_pointer.h>
21#include <deal.II/base/exceptions.h>
22#include <deal.II/base/vectorization.h>
23
24#include <deal.II/lac/diagonal_matrix.h>
25#include <deal.II/lac/la_parallel_vector.h>
26
27#include <deal.II/matrix_free/fe_evaluation.h>
28#include <deal.II/matrix_free/matrix_free.h>
29#include <deal.II/matrix_free/tools.h>
30
31#include <deal.II/multigrid/mg_constrained_dofs.h>
32
33#include <limits>
34
35DEAL_II_NAMESPACE_OPEN
36
37
38namespace MatrixFreeOperators
39{
40 namespace BlockHelper
41 {
42 // workaround for unifying non-block vector and block vector implementations
43 // a non-block vector has one block and the only subblock is the vector
44 // itself
45 template <typename VectorType>
46 std::enable_if_t<IsBlockVector<VectorType>::value, unsigned int>
47 n_blocks(const VectorType &vector)
48 {
49 return vector.n_blocks();
50 }
51
52 template <typename VectorType>
53 std::enable_if_t<!IsBlockVector<VectorType>::value, unsigned int>
54 n_blocks(const VectorType &)
55 {
56 return 1;
57 }
58
59 template <typename VectorType>
60 std::enable_if_t<IsBlockVector<VectorType>::value,
61 typename VectorType::BlockType &>
62 subblock(VectorType &vector, unsigned int block_no)
63 {
64 AssertIndexRange(block_no, vector.n_blocks());
65 return vector.block(block_no);
66 }
67
68 template <typename VectorType>
69 std::enable_if_t<IsBlockVector<VectorType>::value,
70 const typename VectorType::BlockType &>
71 subblock(const VectorType &vector, unsigned int block_no)
72 {
73 AssertIndexRange(block_no, vector.n_blocks());
74 return vector.block(block_no);
75 }
76
77 template <typename VectorType>
78 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType &>
79 subblock(VectorType &vector, unsigned int)
80 {
81 return vector;
82 }
83
84 template <typename VectorType>
85 std::enable_if_t<!IsBlockVector<VectorType>::value, const VectorType &>
86 subblock(const VectorType &vector, unsigned int)
87 {
88 return vector;
89 }
90
91 template <typename VectorType>
92 std::enable_if_t<IsBlockVector<VectorType>::value, void>
93 collect_sizes(VectorType &vector)
94 {
95 vector.collect_sizes();
96 }
97
98 template <typename VectorType>
99 std::enable_if_t<!IsBlockVector<VectorType>::value, void>
100 collect_sizes(const VectorType &)
101 {}
102 } // namespace BlockHelper
103
181 template <int dim,
182 typename VectorType = LinearAlgebra::distributed::Vector<double>,
183 typename VectorizedArrayType =
184 VectorizedArray<typename VectorType::value_type>>
185 class Base : public EnableObserverPointer
186 {
187 public:
191 using value_type = typename VectorType::value_type;
192
196 using size_type = typename VectorType::size_type;
197
201 Base();
202
206 virtual ~Base() override = default;
207
212 virtual void
213 clear();
214
231 void
232 initialize(std::shared_ptr<
233 const MatrixFree<dim, value_type, VectorizedArrayType>> data,
234 const std::vector<unsigned int> &selected_row_blocks =
235 std::vector<unsigned int>(),
236 const std::vector<unsigned int> &selected_column_blocks =
237 std::vector<unsigned int>());
238
252 void
253 initialize(std::shared_ptr<
254 const MatrixFree<dim, value_type, VectorizedArrayType>> data,
255 const MGConstrainedDoFs &mg_constrained_dofs,
256 const unsigned int level,
257 const std::vector<unsigned int> &selected_row_blocks =
258 std::vector<unsigned int>());
259
274 void
275 initialize(std::shared_ptr<
276 const MatrixFree<dim, value_type, VectorizedArrayType>> data_,
277 const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
278 const unsigned int level,
279 const std::vector<unsigned int> &selected_row_blocks =
280 std::vector<unsigned int>());
281
285 size_type
286 m() const;
287
291 size_type
292 n() const;
293
297 void
298 vmult_interface_down(VectorType &dst, const VectorType &src) const;
299
303 void
304 vmult_interface_up(VectorType &dst, const VectorType &src) const;
305
309 void
310 vmult(VectorType &dst, const VectorType &src) const;
311
315 void
316 Tvmult(VectorType &dst, const VectorType &src) const;
317
321 void
322 vmult_add(VectorType &dst, const VectorType &src) const;
323
327 void
328 Tvmult_add(VectorType &dst, const VectorType &src) const;
329
334 value_type
335 el(const unsigned int row, const unsigned int col) const;
336
341 virtual std::size_t
342 memory_consumption() const;
343
347 void
348 initialize_dof_vector(VectorType &vec) const;
349
362 virtual void
363 compute_diagonal() = 0;
364
368 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
369 get_matrix_free() const;
370
374 const std::shared_ptr<DiagonalMatrix<VectorType>> &
375 get_matrix_diagonal_inverse() const;
376
380 const std::shared_ptr<DiagonalMatrix<VectorType>> &
381 get_matrix_diagonal() const;
382
388 void
389 precondition_Jacobi(VectorType &dst,
390 const VectorType &src,
391 const value_type omega) const;
392
393 protected:
398 void
399 preprocess_constraints(VectorType &dst, const VectorType &src) const;
400
405 void
406 postprocess_constraints(VectorType &dst, const VectorType &src) const;
407
412 void
413 set_constrained_entries_to_one(VectorType &dst) const;
414
418 virtual void
419 apply_add(VectorType &dst, const VectorType &src) const = 0;
420
426 virtual void
427 Tapply_add(VectorType &dst, const VectorType &src) const;
428
432 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
433 data;
434
439 std::shared_ptr<DiagonalMatrix<VectorType>> diagonal_entries;
440
445 std::shared_ptr<DiagonalMatrix<VectorType>> inverse_diagonal_entries;
446
451 std::vector<unsigned int> selected_rows;
452
457 std::vector<unsigned int> selected_columns;
458
459 private:
463 std::vector<std::vector<unsigned int>> edge_constrained_indices;
464
468 mutable std::vector<std::vector<std::pair<value_type, value_type>>>
469 edge_constrained_values;
470
475 bool have_interface_matrices;
476
481 void
482 mult_add(VectorType &dst,
483 const VectorType &src,
484 const bool transpose) const;
485
493 void
494 adjust_ghost_range_if_necessary(const VectorType &vec,
495 const bool is_row) const;
496 };
497
498
499
534 template <typename OperatorType>
535 class MGInterfaceOperator : public EnableObserverPointer
536 {
537 public:
541 using value_type = typename OperatorType::value_type;
542
546 using size_type = typename OperatorType::size_type;
547
551 MGInterfaceOperator();
552
556 void
557 clear();
558
562 void
563 initialize(const OperatorType &operator_in);
564
568 template <typename VectorType>
569 void
570 vmult(VectorType &dst, const VectorType &src) const;
571
575 template <typename VectorType>
576 void
577 Tvmult(VectorType &dst, const VectorType &src) const;
578
582 template <typename VectorType>
583 void
584 initialize_dof_vector(VectorType &vec) const;
585
586
587 private:
591 ObserverPointer<const OperatorType> mf_base_operator;
592 };
593
594
595
614 template <int dim,
615 int fe_degree,
616 int n_components = 1,
617 typename Number = double,
618 typename VectorizedArrayType = VectorizedArray<Number>>
619 class CellwiseInverseMassMatrix
620 {
621 static_assert(
622 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
623 "Type of Number and of VectorizedArrayType do not match.");
624
625 public:
630 CellwiseInverseMassMatrix(
631 const FEEvaluationBase<dim,
632 n_components,
633 Number,
634 false,
635 VectorizedArrayType> &fe_eval);
636
646 void
647 apply(const AlignedVector<VectorizedArrayType> &inverse_coefficient,
648 const unsigned int n_actual_components,
649 const VectorizedArrayType *in_array,
650 VectorizedArrayType *out_array) const;
651
663 void
664 apply(const VectorizedArrayType *in_array,
665 VectorizedArrayType *out_array) const;
666
679 void
680 apply(const AlignedVector<Tensor<2, n_components, VectorizedArrayType>>
681 &inverse_dyadic_coefficients,
682 const VectorizedArrayType *in_array,
683 VectorizedArrayType *out_array) const;
684
718 void
719 transform_from_q_points_to_basis(const unsigned int n_actual_components,
720 const VectorizedArrayType *in_array,
721 VectorizedArrayType *out_array) const;
722
728 void
729 fill_inverse_JxW_values(
730 AlignedVector<VectorizedArrayType> &inverse_jxw) const;
731
732 private:
736 const FEEvaluationBase<dim,
737 n_components,
738 Number,
739 false,
740 VectorizedArrayType> &fe_eval;
741 };
742
743
744
752 template <int dim,
753 int fe_degree,
754 int n_q_points_1d = fe_degree + 1,
755 int n_components = 1,
756 typename VectorType = LinearAlgebra::distributed::Vector<double>,
757 typename VectorizedArrayType =
758 VectorizedArray<typename VectorType::value_type>>
759 class MassOperator : public Base<dim, VectorType, VectorizedArrayType>
760 {
761 public:
765 using value_type =
766 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
767
771 using size_type =
772 typename Base<dim, VectorType, VectorizedArrayType>::size_type;
773
777 MassOperator();
778
782 virtual void
783 compute_diagonal() override;
784
801 void
802 compute_lumped_diagonal();
803
807 const std::shared_ptr<DiagonalMatrix<VectorType>> &
808 get_matrix_lumped_diagonal() const;
809
813 const std::shared_ptr<DiagonalMatrix<VectorType>> &
814 get_matrix_lumped_diagonal_inverse() const;
815
816 private:
822 virtual void
823 apply_add(VectorType &dst, const VectorType &src) const override;
824
828 void
829 local_apply_cell(
830 const MatrixFree<dim, value_type, VectorizedArrayType> &data,
831 VectorType &dst,
832 const VectorType &src,
833 const std::pair<unsigned int, unsigned int> &cell_range) const;
834
839 std::shared_ptr<DiagonalMatrix<VectorType>> lumped_diagonal_entries;
840
845 std::shared_ptr<DiagonalMatrix<VectorType>> inverse_lumped_diagonal_entries;
846 };
847
848
849
885 template <int dim,
886 int fe_degree,
887 int n_q_points_1d = fe_degree + 1,
888 int n_components = 1,
889 typename VectorType = LinearAlgebra::distributed::Vector<double>,
890 typename VectorizedArrayType =
891 VectorizedArray<typename VectorType::value_type>>
892 class LaplaceOperator : public Base<dim, VectorType, VectorizedArrayType>
893 {
894 public:
898 using value_type =
899 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
900
904 using size_type =
905 typename Base<dim, VectorType, VectorizedArrayType>::size_type;
906
910 LaplaceOperator();
911
918 virtual void
919 compute_diagonal() override;
920
971 void
972 set_coefficient(
973 const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient);
974
979 virtual void
980 clear() override;
981
988 std::shared_ptr<Table<2, VectorizedArrayType>>
989 get_coefficient();
990
991 private:
997 virtual void
998 apply_add(VectorType &dst, const VectorType &src) const override;
999
1003 void
1004 local_apply_cell(
1005 const MatrixFree<dim, value_type, VectorizedArrayType> &data,
1006 VectorType &dst,
1007 const VectorType &src,
1008 const std::pair<unsigned int, unsigned int> &cell_range) const;
1009
1013 void
1014 local_diagonal_cell(
1015 const MatrixFree<dim, value_type, VectorizedArrayType> &data,
1016 VectorType &dst,
1017 const VectorType &,
1018 const std::pair<unsigned int, unsigned int> &cell_range) const;
1019
1023 template <int n_components_compute>
1024 void
1025 do_operation_on_cell(FEEvaluation<dim,
1026 fe_degree,
1027 n_q_points_1d,
1028 n_components_compute,
1029 value_type,
1030 VectorizedArrayType> &phi,
1031 const unsigned int cell) const;
1032
1036 std::shared_ptr<Table<2, VectorizedArrayType>> scalar_coefficient;
1037 };
1038
1039
1040
1041 // ------------------------------------ inline functions ---------------------
1042
1043 template <int dim,
1044 int fe_degree,
1045 int n_components,
1046 typename Number,
1047 typename VectorizedArrayType>
1048 inline CellwiseInverseMassMatrix<dim,
1049 fe_degree,
1050 n_components,
1051 Number,
1052 VectorizedArrayType>::
1053 CellwiseInverseMassMatrix(
1054 const FEEvaluationBase<dim,
1055 n_components,
1056 Number,
1057 false,
1058 VectorizedArrayType> &fe_eval)
1059 : fe_eval(fe_eval)
1060 {
1061 AssertDimension(fe_eval.get_shape_info().dofs_per_component_on_cell,
1062 fe_eval.get_shape_info().n_q_points);
1063 }
1064
1065
1066
1067 template <int dim,
1068 int fe_degree,
1069 int n_components,
1070 typename Number,
1071 typename VectorizedArrayType>
1072 inline void
1073 CellwiseInverseMassMatrix<dim,
1074 fe_degree,
1075 n_components,
1076 Number,
1077 VectorizedArrayType>::
1078 fill_inverse_JxW_values(
1079 AlignedVector<VectorizedArrayType> &inverse_jxw) const
1080 {
1081 const unsigned int dofs_per_component_on_cell =
1082 (fe_degree > -1) ?
1083 Utilities::pow(fe_degree + 1, dim) :
1084 Utilities::pow(fe_eval.get_shape_info().data.front().fe_degree + 1,
1085 dim);
1086
1087 Assert(inverse_jxw.size() > 0 &&
1088 inverse_jxw.size() % dofs_per_component_on_cell == 0,
1089 ExcMessage(
1090 "Expected diagonal to be a multiple of scalar dof per cells"));
1091
1092 // compute values for the first component
1093 for (unsigned int q = 0; q < dofs_per_component_on_cell; ++q)
1094 inverse_jxw[q] = 1. / fe_eval.JxW(q);
1095 // copy values to rest of vector
1096 for (unsigned int q = dofs_per_component_on_cell; q < inverse_jxw.size();)
1097 for (unsigned int i = 0; i < dofs_per_component_on_cell; ++i, ++q)
1098 inverse_jxw[q] = inverse_jxw[i];
1099 }
1100
1101
1102
1103 template <int dim,
1104 int fe_degree,
1105 int n_components,
1106 typename Number,
1107 typename VectorizedArrayType>
1108 inline void
1109 CellwiseInverseMassMatrix<
1110 dim,
1111 fe_degree,
1112 n_components,
1113 Number,
1114 VectorizedArrayType>::apply(const VectorizedArrayType *in_array,
1115 VectorizedArrayType *out_array) const
1116 {
1117 if (fe_degree > -1)
1118 internal::CellwiseInverseMassMatrixImplBasic<dim, VectorizedArrayType>::
1119 template run<fe_degree>(n_components, fe_eval, in_array, out_array);
1120 else
1121 internal::CellwiseInverseMassFactory<dim, VectorizedArrayType>::apply(
1122 n_components, fe_eval, in_array, out_array);
1123 }
1124
1125
1126
1127 template <int dim,
1128 int fe_degree,
1129 int n_components,
1130 typename Number,
1131 typename VectorizedArrayType>
1132 inline void
1133 CellwiseInverseMassMatrix<dim,
1134 fe_degree,
1135 n_components,
1136 Number,
1137 VectorizedArrayType>::
1138 apply(const AlignedVector<VectorizedArrayType> &inverse_coefficients,
1139 const unsigned int n_actual_components,
1140 const VectorizedArrayType *in_array,
1141 VectorizedArrayType *out_array) const
1142 {
1143 if (fe_degree > -1)
1144 internal::CellwiseInverseMassMatrixImplFlexible<
1145 dim,
1146 VectorizedArrayType>::template run<fe_degree>(n_actual_components,
1147 fe_eval,
1148 make_array_view(
1149 inverse_coefficients),
1150 false,
1151 in_array,
1152 out_array);
1153 else
1154 internal::CellwiseInverseMassFactory<dim, VectorizedArrayType>::apply(
1155 n_actual_components,
1156 fe_eval,
1157 make_array_view(inverse_coefficients),
1158 false,
1159 in_array,
1160 out_array);
1161 }
1162
1163 template <int dim,
1164 int fe_degree,
1165 int n_components,
1166 typename Number,
1167 typename VectorizedArrayType>
1168 inline void
1169 CellwiseInverseMassMatrix<dim,
1170 fe_degree,
1171 n_components,
1172 Number,
1173 VectorizedArrayType>::
1174 apply(const AlignedVector<Tensor<2, n_components, VectorizedArrayType>>
1175 &inverse_dyadic_coefficients,
1176 const VectorizedArrayType *in_array,
1177 VectorizedArrayType *out_array) const
1178 {
1179 const unsigned int unrolled_size =
1180 inverse_dyadic_coefficients.size() * (n_components * n_components);
1181
1182 if (fe_degree > -1)
1183 internal::CellwiseInverseMassMatrixImplFlexible<dim,
1184 VectorizedArrayType>::
1185 template run<fe_degree>(n_components,
1186 fe_eval,
1187 ArrayView<const VectorizedArrayType>(
1188 &inverse_dyadic_coefficients[0][0][0],
1189 unrolled_size),
1190 true,
1191 in_array,
1192 out_array);
1193 else
1194 internal::CellwiseInverseMassFactory<dim, VectorizedArrayType>::apply(
1195 n_components,
1196 fe_eval,
1197 ArrayView<const VectorizedArrayType>(
1198 &inverse_dyadic_coefficients[0][0][0], unrolled_size),
1199 true,
1200 in_array,
1201 out_array);
1202 }
1203
1204
1205
1206 template <int dim,
1207 int fe_degree,
1208 int n_components,
1209 typename Number,
1210 typename VectorizedArrayType>
1211 inline void
1212 CellwiseInverseMassMatrix<dim,
1213 fe_degree,
1214 n_components,
1215 Number,
1216 VectorizedArrayType>::
1217 transform_from_q_points_to_basis(const unsigned int n_actual_components,
1218 const VectorizedArrayType *in_array,
1219 VectorizedArrayType *out_array) const
1220 {
1221 const auto n_q_points_1d = fe_eval.get_shape_info().data[0].n_q_points_1d;
1222
1223 if (fe_degree > -1 && (fe_degree + 1 == n_q_points_1d))
1224 internal::CellwiseInverseMassMatrixImplTransformFromQPoints<
1225 dim,
1226 VectorizedArrayType>::template run<fe_degree,
1227 fe_degree + 1>(n_actual_components,
1228 fe_eval,
1229 in_array,
1230 out_array);
1231 else
1232 internal::CellwiseInverseMassFactory<dim, VectorizedArrayType>::
1233 transform_from_q_points_to_basis(n_actual_components,
1234 fe_eval,
1235 in_array,
1236 out_array);
1237 }
1238
1239
1240
1241 //----------------- Base operator -----------------------------
1242 template <int dim, typename VectorType, typename VectorizedArrayType>
1243 Base<dim, VectorType, VectorizedArrayType>::Base()
1244 : EnableObserverPointer()
1245 , have_interface_matrices(false)
1246 {}
1247
1248
1249
1250 template <int dim, typename VectorType, typename VectorizedArrayType>
1251 typename Base<dim, VectorType, VectorizedArrayType>::size_type
1252 Base<dim, VectorType, VectorizedArrayType>::m() const
1253 {
1254 Assert(data.get() != nullptr, ExcNotInitialized());
1255 typename Base<dim, VectorType, VectorizedArrayType>::size_type total_size =
1256 0;
1257 for (const unsigned int selected_row : selected_rows)
1258 total_size += data->get_vector_partitioner(selected_row)->size();
1259 return total_size;
1260 }
1261
1262
1263
1264 template <int dim, typename VectorType, typename VectorizedArrayType>
1265 typename Base<dim, VectorType, VectorizedArrayType>::size_type
1266 Base<dim, VectorType, VectorizedArrayType>::n() const
1267 {
1268 Assert(data.get() != nullptr, ExcNotInitialized());
1269 typename Base<dim, VectorType, VectorizedArrayType>::size_type total_size =
1270 0;
1271 for (const unsigned int selected_column : selected_columns)
1272 total_size += data->get_vector_partitioner(selected_column)->size();
1273 return total_size;
1274 }
1275
1276
1277
1278 template <int dim, typename VectorType, typename VectorizedArrayType>
1279 void
1280 Base<dim, VectorType, VectorizedArrayType>::clear()
1281 {
1282 data.reset();
1283 inverse_diagonal_entries.reset();
1284 }
1285
1286
1287
1288 template <int dim, typename VectorType, typename VectorizedArrayType>
1289 typename Base<dim, VectorType, VectorizedArrayType>::value_type
1290 Base<dim, VectorType, VectorizedArrayType>::el(const unsigned int row,
1291 const unsigned int col) const
1292 {
1293 Assert(row == col, ExcNotImplemented());
1294 Assert(inverse_diagonal_entries.get() != nullptr &&
1295 inverse_diagonal_entries->m() > 0,
1296 ExcNotInitialized());
1297 return 1.0 / (*inverse_diagonal_entries)(row, row);
1298 }
1299
1300
1301
1302 template <int dim, typename VectorType, typename VectorizedArrayType>
1303 void
1304 Base<dim, VectorType, VectorizedArrayType>::initialize_dof_vector(
1305 VectorType &vec) const
1306 {
1307 Assert(data.get() != nullptr, ExcNotInitialized());
1308 AssertDimension(BlockHelper::n_blocks(vec), selected_rows.size());
1309 for (unsigned int i = 0; i < BlockHelper::n_blocks(vec); ++i)
1310 {
1311 const unsigned int index = selected_rows[i];
1312 if (!BlockHelper::subblock(vec, index)
1313 .partitioners_are_compatible(
1314 *data->get_dof_info(index).vector_partitioner))
1315 data->initialize_dof_vector(BlockHelper::subblock(vec, index), index);
1316
1317 Assert(BlockHelper::subblock(vec, index)
1318 .partitioners_are_globally_compatible(
1319 *data->get_dof_info(index).vector_partitioner),
1320 ExcInternalError());
1321 }
1322 BlockHelper::collect_sizes(vec);
1323 }
1324
1325
1326
1327 template <int dim, typename VectorType, typename VectorizedArrayType>
1328 void
1329 Base<dim, VectorType, VectorizedArrayType>::initialize(
1330 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
1331 data_,
1332 const std::vector<unsigned int> &given_row_selection,
1333 const std::vector<unsigned int> &given_column_selection)
1334 {
1335 data = data_;
1336
1337 selected_rows.clear();
1338 selected_columns.clear();
1339 if (given_row_selection.empty())
1340 for (unsigned int i = 0; i < data_->n_components(); ++i)
1341 selected_rows.push_back(i);
1342 else
1343 {
1344 for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1345 {
1346 AssertIndexRange(given_row_selection[i], data_->n_components());
1347 for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1348 if (j != i)
1349 Assert(given_row_selection[j] != given_row_selection[i],
1350 ExcMessage("Given row indices must be unique"));
1351
1352 selected_rows.push_back(given_row_selection[i]);
1353 }
1354 }
1355 if (given_column_selection.empty())
1356 selected_columns = selected_rows;
1357 else
1358 {
1359 for (unsigned int i = 0; i < given_column_selection.size(); ++i)
1360 {
1361 AssertIndexRange(given_column_selection[i], data_->n_components());
1362 for (unsigned int j = 0; j < given_column_selection.size(); ++j)
1363 if (j != i)
1364 Assert(given_column_selection[j] != given_column_selection[i],
1365 ExcMessage("Given column indices must be unique"));
1366
1367 selected_columns.push_back(given_column_selection[i]);
1368 }
1369 }
1370
1371 edge_constrained_indices.clear();
1372 edge_constrained_indices.resize(selected_rows.size());
1373 edge_constrained_values.clear();
1374 edge_constrained_values.resize(selected_rows.size());
1375 have_interface_matrices = false;
1376 }
1377
1378
1379
1380 template <int dim, typename VectorType, typename VectorizedArrayType>
1381 void
1382 Base<dim, VectorType, VectorizedArrayType>::initialize(
1383 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
1384 data_,
1385 const MGConstrainedDoFs &mg_constrained_dofs,
1386 const unsigned int level,
1387 const std::vector<unsigned int> &given_row_selection)
1388 {
1389 std::vector<MGConstrainedDoFs> mg_constrained_dofs_vector(
1390 1, mg_constrained_dofs);
1391 initialize(data_, mg_constrained_dofs_vector, level, given_row_selection);
1392 }
1393
1394
1395
1396 template <int dim, typename VectorType, typename VectorizedArrayType>
1397 void
1398 Base<dim, VectorType, VectorizedArrayType>::initialize(
1399 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
1400 data_,
1401 const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
1402 const unsigned int level,
1403 const std::vector<unsigned int> &given_row_selection)
1404 {
1405 AssertThrow(level != numbers::invalid_unsigned_int,
1406 ExcMessage("level is not set"));
1407
1408 selected_rows.clear();
1409 selected_columns.clear();
1410 if (given_row_selection.empty())
1411 for (unsigned int i = 0; i < data_->n_components(); ++i)
1412 selected_rows.push_back(i);
1413 else
1414 {
1415 for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1416 {
1417 AssertIndexRange(given_row_selection[i], data_->n_components());
1418 for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1419 if (j != i)
1420 Assert(given_row_selection[j] != given_row_selection[i],
1421 ExcMessage("Given row indices must be unique"));
1422
1423 selected_rows.push_back(given_row_selection[i]);
1424 }
1425 }
1426 selected_columns = selected_rows;
1427
1428 AssertDimension(mg_constrained_dofs.size(), selected_rows.size());
1429 edge_constrained_indices.clear();
1430 edge_constrained_indices.resize(selected_rows.size());
1431 edge_constrained_values.clear();
1432 edge_constrained_values.resize(selected_rows.size());
1433
1434 data = data_;
1435
1436 for (unsigned int j = 0; j < selected_rows.size(); ++j)
1437 {
1438 if (data_->n_cell_batches() > 0)
1439 {
1440 AssertDimension(level, data_->get_cell_iterator(0, 0, j)->level());
1441 }
1442
1443 // setup edge_constrained indices
1444 const std::vector<types::global_dof_index> interface_indices =
1445 mg_constrained_dofs[j]
1446 .get_refinement_edge_indices(level)
1447 .get_index_vector();
1448 edge_constrained_indices[j].clear();
1449 edge_constrained_indices[j].reserve(interface_indices.size());
1450 edge_constrained_values[j].resize(interface_indices.size());
1451 const IndexSet &locally_owned =
1452 data->get_dof_handler(selected_rows[j]).locally_owned_mg_dofs(level);
1453 for (const auto interface_index : interface_indices)
1454 if (locally_owned.is_element(interface_index))
1455 edge_constrained_indices[j].push_back(
1456 locally_owned.index_within_set(interface_index));
1457 have_interface_matrices |=
1458 Utilities::MPI::max(
1459 static_cast<unsigned int>(edge_constrained_indices[j].size()),
1460 data->get_vector_partitioner()->get_mpi_communicator()) > 0;
1461 }
1462 }
1463
1464
1465
1466 template <int dim, typename VectorType, typename VectorizedArrayType>
1467 void
1468 Base<dim, VectorType, VectorizedArrayType>::set_constrained_entries_to_one(
1469 VectorType &dst) const
1470 {
1471 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1472 {
1473 const std::vector<unsigned int> &constrained_dofs =
1474 data->get_constrained_dofs(selected_rows[j]);
1475 for (const auto constrained_dof : constrained_dofs)
1476 BlockHelper::subblock(dst, j).local_element(constrained_dof) = 1.;
1477 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1478 BlockHelper::subblock(dst, j).local_element(
1479 edge_constrained_indices[j][i]) = 1.;
1480 }
1481 }
1482
1483
1484
1485 template <int dim, typename VectorType, typename VectorizedArrayType>
1486 void
1487 Base<dim, VectorType, VectorizedArrayType>::vmult(VectorType &dst,
1488 const VectorType &src) const
1489 {
1490 using Number =
1491 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1492 dst = Number(0.);
1493 vmult_add(dst, src);
1494 }
1495
1496
1497
1498 template <int dim, typename VectorType, typename VectorizedArrayType>
1499 void
1500 Base<dim, VectorType, VectorizedArrayType>::vmult_add(
1501 VectorType &dst,
1502 const VectorType &src) const
1503 {
1504 mult_add(dst, src, false);
1505 }
1506
1507
1508
1509 template <int dim, typename VectorType, typename VectorizedArrayType>
1510 void
1511 Base<dim, VectorType, VectorizedArrayType>::Tvmult_add(
1512 VectorType &dst,
1513 const VectorType &src) const
1514 {
1515 mult_add(dst, src, true);
1516 }
1517
1518
1519
1520 template <int dim, typename VectorType, typename VectorizedArrayType>
1521 void
1522 Base<dim, VectorType, VectorizedArrayType>::adjust_ghost_range_if_necessary(
1523 const VectorType &src,
1524 const bool is_row) const
1525 {
1526 using Number =
1527 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1528 for (unsigned int i = 0; i < BlockHelper::n_blocks(src); ++i)
1529 {
1530 const unsigned int mf_component =
1531 is_row ? selected_rows[i] : selected_columns[i];
1532 // If both vectors use the same partitioner -> done
1533 if (BlockHelper::subblock(src, i).get_partitioner().get() ==
1534 data->get_dof_info(mf_component).vector_partitioner.get())
1535 continue;
1536
1537 // If not, assert that the local ranges are the same and reset to the
1538 // current partitioner
1539 Assert(BlockHelper::subblock(src, i)
1540 .get_partitioner()
1541 ->locally_owned_size() ==
1542 data->get_dof_info(mf_component)
1543 .vector_partitioner->locally_owned_size(),
1544 ExcMessage(
1545 "The vector passed to the vmult() function does not have "
1546 "the correct size for compatibility with MatrixFree."));
1547
1548 // copy the vector content to a temporary vector so that it does not get
1549 // lost
1550 LinearAlgebra::distributed::Vector<Number> copy_vec(
1551 BlockHelper::subblock(src, i));
1552 this->data->initialize_dof_vector(
1553 BlockHelper::subblock(const_cast<VectorType &>(src), i),
1554 mf_component);
1555 BlockHelper::subblock(const_cast<VectorType &>(src), i)
1556 .copy_locally_owned_data_from(copy_vec);
1557 }
1558 }
1559
1560
1561
1562 template <int dim, typename VectorType, typename VectorizedArrayType>
1563 void
1564 Base<dim, VectorType, VectorizedArrayType>::preprocess_constraints(
1565 VectorType &dst,
1566 const VectorType &src) const
1567 {
1568 using Number =
1569 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1570 adjust_ghost_range_if_necessary(src, false);
1571 adjust_ghost_range_if_necessary(dst, true);
1572
1573 // set zero Dirichlet values on the input vector (and remember the src and
1574 // dst values because we need to reset them at the end)
1575 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1576 {
1577 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1578 {
1579 edge_constrained_values[j][i] = std::pair<Number, Number>(
1580 BlockHelper::subblock(src, j).local_element(
1581 edge_constrained_indices[j][i]),
1582 BlockHelper::subblock(dst, j).local_element(
1583 edge_constrained_indices[j][i]));
1584 BlockHelper::subblock(const_cast<VectorType &>(src), j)
1585 .local_element(edge_constrained_indices[j][i]) = 0.;
1586 }
1587 }
1588 }
1589
1590
1591
1592 template <int dim, typename VectorType, typename VectorizedArrayType>
1593 void
1594 Base<dim, VectorType, VectorizedArrayType>::mult_add(
1595 VectorType &dst,
1596 const VectorType &src,
1597 const bool transpose) const
1598 {
1599 AssertDimension(dst.size(), src.size());
1600 AssertDimension(BlockHelper::n_blocks(dst), BlockHelper::n_blocks(src));
1601 AssertDimension(BlockHelper::n_blocks(dst), selected_rows.size());
1602 preprocess_constraints(dst, src);
1603 if (transpose)
1604 Tapply_add(dst, src);
1605 else
1606 apply_add(dst, src);
1607 postprocess_constraints(dst, src);
1608 }
1609
1610
1611
1612 template <int dim, typename VectorType, typename VectorizedArrayType>
1613 void
1614 Base<dim, VectorType, VectorizedArrayType>::postprocess_constraints(
1615 VectorType &dst,
1616 const VectorType &src) const
1617 {
1618 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1619 {
1620 const std::vector<unsigned int> &constrained_dofs =
1621 data->get_constrained_dofs(selected_rows[j]);
1622 for (const auto constrained_dof : constrained_dofs)
1623 BlockHelper::subblock(dst, j).local_element(constrained_dof) +=
1624 BlockHelper::subblock(src, j).local_element(constrained_dof);
1625 }
1626
1627 // reset edge constrained values, multiply by unit matrix and add into
1628 // destination
1629 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1630 {
1631 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1632 {
1633 BlockHelper::subblock(const_cast<VectorType &>(src), j)
1634 .local_element(edge_constrained_indices[j][i]) =
1635 edge_constrained_values[j][i].first;
1636 BlockHelper::subblock(dst, j).local_element(
1637 edge_constrained_indices[j][i]) =
1638 edge_constrained_values[j][i].second +
1639 edge_constrained_values[j][i].first;
1640 }
1641 }
1642 }
1643
1644
1645
1646 template <int dim, typename VectorType, typename VectorizedArrayType>
1647 void
1648 Base<dim, VectorType, VectorizedArrayType>::vmult_interface_down(
1649 VectorType &dst,
1650 const VectorType &src) const
1651 {
1652 using Number =
1653 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1654 AssertDimension(dst.size(), src.size());
1655 adjust_ghost_range_if_necessary(src, false);
1656 adjust_ghost_range_if_necessary(dst, true);
1657
1658 dst = Number(0.);
1659
1660 if (!have_interface_matrices)
1661 return;
1662
1663 // set zero Dirichlet values on the input vector (and remember the src and
1664 // dst values because we need to reset them at the end)
1665 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1666 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1667 {
1668 edge_constrained_values[j][i] = std::pair<Number, Number>(
1669 BlockHelper::subblock(src, j).local_element(
1670 edge_constrained_indices[j][i]),
1671 BlockHelper::subblock(dst, j).local_element(
1672 edge_constrained_indices[j][i]));
1673 BlockHelper::subblock(const_cast<VectorType &>(src), j)
1674 .local_element(edge_constrained_indices[j][i]) = 0.;
1675 }
1676
1677 apply_add(dst, src);
1678
1679 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1680 {
1681 unsigned int c = 0;
1682 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1683 {
1684 for (; c < edge_constrained_indices[j][i]; ++c)
1685 BlockHelper::subblock(dst, j).local_element(c) = 0.;
1686 ++c;
1687
1688 // reset the src values
1689 BlockHelper::subblock(const_cast<VectorType &>(src), j)
1690 .local_element(edge_constrained_indices[j][i]) =
1691 edge_constrained_values[j][i].first;
1692 }
1693 for (; c < BlockHelper::subblock(dst, j).locally_owned_size(); ++c)
1694 BlockHelper::subblock(dst, j).local_element(c) = 0.;
1695 }
1696 }
1697
1698
1699
1700 template <int dim, typename VectorType, typename VectorizedArrayType>
1701 void
1702 Base<dim, VectorType, VectorizedArrayType>::vmult_interface_up(
1703 VectorType &dst,
1704 const VectorType &src) const
1705 {
1706 using Number =
1707 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1708 AssertDimension(dst.size(), src.size());
1709 adjust_ghost_range_if_necessary(src, false);
1710 adjust_ghost_range_if_necessary(dst, true);
1711
1712 dst = Number(0.);
1713
1714 if (!have_interface_matrices)
1715 return;
1716
1717 VectorType src_cpy(src);
1718 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1719 {
1720 unsigned int c = 0;
1721 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1722 {
1723 for (; c < edge_constrained_indices[j][i]; ++c)
1724 BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1725 ++c;
1726 }
1727 for (; c < BlockHelper::subblock(src_cpy, j).locally_owned_size(); ++c)
1728 BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1729 }
1730
1731 apply_add(dst, src_cpy);
1732
1733 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1734 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1735 BlockHelper::subblock(dst, j).local_element(
1736 edge_constrained_indices[j][i]) = 0.;
1737 }
1738
1739
1740
1741 template <int dim, typename VectorType, typename VectorizedArrayType>
1742 void
1743 Base<dim, VectorType, VectorizedArrayType>::Tvmult(
1744 VectorType &dst,
1745 const VectorType &src) const
1746 {
1747 using Number =
1748 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1749 dst = Number(0.);
1750 Tvmult_add(dst, src);
1751 }
1752
1753
1754
1755 template <int dim, typename VectorType, typename VectorizedArrayType>
1756 std::size_t
1757 Base<dim, VectorType, VectorizedArrayType>::memory_consumption() const
1758 {
1759 return inverse_diagonal_entries.get() != nullptr ?
1760 inverse_diagonal_entries->memory_consumption() :
1761 sizeof(*this);
1762 }
1763
1764
1765
1766 template <int dim, typename VectorType, typename VectorizedArrayType>
1767 std::shared_ptr<const MatrixFree<
1768 dim,
1769 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
1770 VectorizedArrayType>>
1775
1776
1777
1778 template <int dim, typename VectorType, typename VectorizedArrayType>
1779 const std::shared_ptr<DiagonalMatrix<VectorType>> &
1780 Base<dim, VectorType, VectorizedArrayType>::get_matrix_diagonal_inverse()
1781 const
1782 {
1783 Assert(inverse_diagonal_entries.get() != nullptr &&
1784 inverse_diagonal_entries->m() > 0,
1785 ExcNotInitialized());
1786 return inverse_diagonal_entries;
1787 }
1788
1789
1790
1791 template <int dim, typename VectorType, typename VectorizedArrayType>
1792 const std::shared_ptr<DiagonalMatrix<VectorType>> &
1793 Base<dim, VectorType, VectorizedArrayType>::get_matrix_diagonal() const
1794 {
1795 Assert(diagonal_entries.get() != nullptr && diagonal_entries->m() > 0,
1796 ExcNotInitialized());
1797 return diagonal_entries;
1798 }
1799
1800
1801
1802 template <int dim, typename VectorType, typename VectorizedArrayType>
1803 void
1804 Base<dim, VectorType, VectorizedArrayType>::Tapply_add(
1805 VectorType &dst,
1806 const VectorType &src) const
1807 {
1808 apply_add(dst, src);
1809 }
1810
1811
1812
1813 template <int dim, typename VectorType, typename VectorizedArrayType>
1814 void
1815 Base<dim, VectorType, VectorizedArrayType>::precondition_Jacobi(
1816 VectorType &dst,
1817 const VectorType &src,
1818 const typename Base<dim, VectorType, VectorizedArrayType>::value_type omega)
1819 const
1820 {
1821 Assert(inverse_diagonal_entries.get() && inverse_diagonal_entries->m() > 0,
1822 ExcNotInitialized());
1823 inverse_diagonal_entries->vmult(dst, src);
1824 dst *= omega;
1825 }
1826
1827
1828
1829 //------------------------- MGInterfaceOperator ------------------------------
1830
1831 template <typename OperatorType>
1832 MGInterfaceOperator<OperatorType>::MGInterfaceOperator()
1833 : EnableObserverPointer()
1834 , mf_base_operator(nullptr)
1835 {}
1836
1837
1838
1839 template <typename OperatorType>
1840 void
1841 MGInterfaceOperator<OperatorType>::clear()
1842 {
1843 mf_base_operator = nullptr;
1844 }
1845
1846
1847
1848 template <typename OperatorType>
1849 void
1850 MGInterfaceOperator<OperatorType>::initialize(const OperatorType &operator_in)
1851 {
1852 mf_base_operator = &operator_in;
1853 }
1854
1855
1856
1857 template <typename OperatorType>
1858 template <typename VectorType>
1859 void
1860 MGInterfaceOperator<OperatorType>::vmult(VectorType &dst,
1861 const VectorType &src) const
1862 {
1863#ifndef DEAL_II_MSVC
1864 static_assert(
1865 std::is_same_v<typename VectorType::value_type, value_type>,
1866 "The vector type must be based on the same value type as this "
1867 "operator");
1868#endif
1869
1870 Assert(mf_base_operator != nullptr, ExcNotInitialized());
1871
1872 mf_base_operator->vmult_interface_down(dst, src);
1873 }
1874
1875
1876
1877 template <typename OperatorType>
1878 template <typename VectorType>
1879 void
1880 MGInterfaceOperator<OperatorType>::Tvmult(VectorType &dst,
1881 const VectorType &src) const
1882 {
1883#ifndef DEAL_II_MSVC
1884 static_assert(
1885 std::is_same_v<typename VectorType::value_type, value_type>,
1886 "The vector type must be based on the same value type as this "
1887 "operator");
1888#endif
1889
1890 Assert(mf_base_operator != nullptr, ExcNotInitialized());
1891
1892 mf_base_operator->vmult_interface_up(dst, src);
1893 }
1894
1895
1896
1897 template <typename OperatorType>
1898 template <typename VectorType>
1899 void
1900 MGInterfaceOperator<OperatorType>::initialize_dof_vector(
1901 VectorType &vec) const
1902 {
1903 Assert(mf_base_operator != nullptr, ExcNotInitialized());
1904
1905 mf_base_operator->initialize_dof_vector(vec);
1906 }
1907
1908
1909
1910 //-----------------------------MassOperator----------------------------------
1911
1912 template <int dim,
1913 int fe_degree,
1914 int n_q_points_1d,
1915 int n_components,
1916 typename VectorType,
1917 typename VectorizedArrayType>
1918 MassOperator<dim,
1919 fe_degree,
1920 n_q_points_1d,
1921 n_components,
1922 VectorType,
1923 VectorizedArrayType>::MassOperator()
1924 : Base<dim, VectorType, VectorizedArrayType>()
1925 {
1926 AssertThrow(
1927 IsBlockVector<VectorType>::value == false,
1928 ExcNotImplemented(
1929 "This class only supports the non-blocked vector variant of the Base "
1930 "operator because only a single FEEvaluation object is used in the "
1931 "apply function."));
1932 }
1933
1934
1935
1936 template <int dim,
1937 int fe_degree,
1938 int n_q_points_1d,
1939 int n_components,
1940 typename VectorType,
1941 typename VectorizedArrayType>
1942 void
1943 MassOperator<dim,
1944 fe_degree,
1945 n_q_points_1d,
1946 n_components,
1947 VectorType,
1948 VectorizedArrayType>::compute_diagonal()
1949 {
1950 Assert((Base<dim, VectorType, VectorizedArrayType>::data.get() != nullptr),
1951 ExcNotInitialized());
1952 Assert(this->selected_rows == this->selected_columns,
1953 ExcMessage("This function is only implemented for square (not "
1954 "rectangular) operators."));
1955
1956 this->inverse_diagonal_entries =
1957 std::make_shared<DiagonalMatrix<VectorType>>();
1958 this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1959 VectorType &inverse_diagonal_vector =
1960 this->inverse_diagonal_entries->get_vector();
1961 VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1962 this->initialize_dof_vector(inverse_diagonal_vector);
1963 this->initialize_dof_vector(diagonal_vector);
1964
1965 // Set up the action of the mass matrix in a way that's compatible with
1966 // MatrixFreeTools::compute_diagonal:
1967 auto diagonal_evaluation = [](auto &integrator) {
1968 integrator.evaluate(EvaluationFlags::values);
1969 for (unsigned int q = 0; q < integrator.n_q_points; ++q)
1970 integrator.submit_value(integrator.get_value(q), q);
1971 integrator.integrate(EvaluationFlags::values);
1972 };
1973
1974 std::function<void(
1975 FEEvaluation<
1976 dim,
1977 fe_degree,
1978 n_q_points_1d,
1979 n_components,
1980 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
1981 VectorizedArrayType> &)>
1982 diagonal_evaluation_f(diagonal_evaluation);
1983
1984 Assert(this->selected_rows.size() > 0, ExcInternalError());
1985 for (unsigned int block_n = 0; block_n < this->selected_rows.size();
1986 ++block_n)
1987 MatrixFreeTools::compute_diagonal(*this->data,
1988 BlockHelper::subblock(diagonal_vector,
1989 block_n),
1990 diagonal_evaluation_f,
1991 this->selected_rows[block_n]);
1992
1993 // Constrained entries will create zeros on the main diagonal, which we
1994 // don't want
1995 this->set_constrained_entries_to_one(diagonal_vector);
1996
1997 inverse_diagonal_vector = diagonal_vector;
1998
1999 for (unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
2000 ++i)
2001 {
2002 if constexpr (running_in_debug_mode())
2003 {
2004 // only define the type alias in debug mode to avoid a warning
2005 using Number =
2006 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2007 Assert(diagonal_vector.local_element(i) > Number(0),
2008 ExcInternalError());
2009 }
2010 inverse_diagonal_vector.local_element(i) =
2011 1. / inverse_diagonal_vector.local_element(i);
2012 }
2013
2014 // We never need ghost values so don't update them
2015 }
2016
2017
2018
2019 template <int dim,
2020 int fe_degree,
2021 int n_q_points_1d,
2022 int n_components,
2023 typename VectorType,
2024 typename VectorizedArrayType>
2025 void
2026 MassOperator<dim,
2027 fe_degree,
2028 n_q_points_1d,
2029 n_components,
2030 VectorType,
2031 VectorizedArrayType>::compute_lumped_diagonal()
2032 {
2033 using Number =
2034 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2035 Assert((Base<dim, VectorType, VectorizedArrayType>::data.get() != nullptr),
2036 ExcNotInitialized());
2037 Assert(this->selected_rows == this->selected_columns,
2038 ExcMessage("This function is only implemented for square (not "
2039 "rectangular) operators."));
2040
2041 inverse_lumped_diagonal_entries =
2042 std::make_shared<DiagonalMatrix<VectorType>>();
2043 lumped_diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
2044 VectorType &inverse_lumped_diagonal_vector =
2045 inverse_lumped_diagonal_entries->get_vector();
2046 VectorType &lumped_diagonal_vector = lumped_diagonal_entries->get_vector();
2047 this->initialize_dof_vector(inverse_lumped_diagonal_vector);
2048 this->initialize_dof_vector(lumped_diagonal_vector);
2049
2050 // Re-use the inverse_lumped_diagonal_vector as the vector of 1s
2051 inverse_lumped_diagonal_vector = Number(1.);
2052 apply_add(lumped_diagonal_vector, inverse_lumped_diagonal_vector);
2053 this->set_constrained_entries_to_one(lumped_diagonal_vector);
2054
2055 const size_type locally_owned_size =
2056 inverse_lumped_diagonal_vector.locally_owned_size();
2057 // A caller may request a lumped diagonal matrix when it doesn't make sense
2058 // (e.g., an element with negative-mean basis functions). Avoid division by
2059 // zero so we don't cause a floating point exception but permit negative
2060 // entries here.
2061 for (size_type i = 0; i < locally_owned_size; ++i)
2062 {
2063 if (lumped_diagonal_vector.local_element(i) == Number(0.))
2064 inverse_lumped_diagonal_vector.local_element(i) = Number(1.);
2065 else
2066 inverse_lumped_diagonal_vector.local_element(i) =
2067 Number(1.) / lumped_diagonal_vector.local_element(i);
2068 }
2069
2070 inverse_lumped_diagonal_vector.update_ghost_values();
2071 lumped_diagonal_vector.update_ghost_values();
2072 }
2073
2074
2075
2076 template <int dim,
2077 int fe_degree,
2078 int n_q_points_1d,
2079 int n_components,
2080 typename VectorType,
2081 typename VectorizedArrayType>
2082 const std::shared_ptr<DiagonalMatrix<VectorType>> &
2083 MassOperator<dim,
2084 fe_degree,
2085 n_q_points_1d,
2086 n_components,
2087 VectorType,
2088 VectorizedArrayType>::get_matrix_lumped_diagonal_inverse() const
2089 {
2090 Assert(inverse_lumped_diagonal_entries.get() != nullptr &&
2091 inverse_lumped_diagonal_entries->m() > 0,
2092 ExcNotInitialized());
2093 return inverse_lumped_diagonal_entries;
2094 }
2095
2096
2097
2098 template <int dim,
2099 int fe_degree,
2100 int n_q_points_1d,
2101 int n_components,
2102 typename VectorType,
2103 typename VectorizedArrayType>
2104 const std::shared_ptr<DiagonalMatrix<VectorType>> &
2105 MassOperator<dim,
2106 fe_degree,
2107 n_q_points_1d,
2108 n_components,
2109 VectorType,
2110 VectorizedArrayType>::get_matrix_lumped_diagonal() const
2111 {
2112 Assert(lumped_diagonal_entries.get() != nullptr &&
2113 lumped_diagonal_entries->m() > 0,
2114 ExcNotInitialized());
2115 return lumped_diagonal_entries;
2116 }
2117
2118
2119
2120 template <int dim,
2121 int fe_degree,
2122 int n_q_points_1d,
2123 int n_components,
2124 typename VectorType,
2125 typename VectorizedArrayType>
2126 void
2127 MassOperator<dim,
2128 fe_degree,
2129 n_q_points_1d,
2130 n_components,
2131 VectorType,
2132 VectorizedArrayType>::apply_add(VectorType &dst,
2133 const VectorType &src) const
2134 {
2135 Base<dim, VectorType, VectorizedArrayType>::data->cell_loop(
2136 &MassOperator::local_apply_cell, this, dst, src);
2137 }
2138
2139
2140
2141 template <int dim,
2142 int fe_degree,
2143 int n_q_points_1d,
2144 int n_components,
2145 typename VectorType,
2146 typename VectorizedArrayType>
2147 void
2148 MassOperator<dim,
2149 fe_degree,
2150 n_q_points_1d,
2151 n_components,
2152 VectorType,
2153 VectorizedArrayType>::
2154 local_apply_cell(
2155 const MatrixFree<
2156 dim,
2157 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2158 VectorizedArrayType> &data,
2159 VectorType &dst,
2160 const VectorType &src,
2161 const std::pair<unsigned int, unsigned int> &cell_range) const
2162 {
2163 using Number =
2164 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2165 FEEvaluation<dim,
2166 fe_degree,
2167 n_q_points_1d,
2168 n_components,
2169 Number,
2170 VectorizedArrayType>
2171 phi(data, this->selected_rows[0]);
2172 for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2173 {
2174 phi.reinit(cell);
2175 phi.read_dof_values(src);
2176 phi.evaluate(EvaluationFlags::values);
2177 for (unsigned int q = 0; q < phi.n_q_points; ++q)
2178 phi.submit_value(phi.get_value(q), q);
2179 phi.integrate(EvaluationFlags::values);
2180 phi.distribute_local_to_global(dst);
2181 }
2182 }
2183
2184
2185 //-----------------------------LaplaceOperator----------------------------------
2186
2187 template <int dim,
2188 int fe_degree,
2189 int n_q_points_1d,
2190 int n_components,
2191 typename VectorType,
2192 typename VectorizedArrayType>
2193 LaplaceOperator<dim,
2194 fe_degree,
2195 n_q_points_1d,
2196 n_components,
2197 VectorType,
2198 VectorizedArrayType>::LaplaceOperator()
2199 : Base<dim, VectorType, VectorizedArrayType>()
2200 {}
2201
2202
2203
2204 template <int dim,
2205 int fe_degree,
2206 int n_q_points_1d,
2207 int n_components,
2208 typename VectorType,
2209 typename VectorizedArrayType>
2210 void
2211 LaplaceOperator<dim,
2212 fe_degree,
2213 n_q_points_1d,
2214 n_components,
2215 VectorType,
2216 VectorizedArrayType>::clear()
2217 {
2218 Base<dim, VectorType, VectorizedArrayType>::clear();
2219 scalar_coefficient.reset();
2220 }
2221
2222
2223
2224 template <int dim,
2225 int fe_degree,
2226 int n_q_points_1d,
2227 int n_components,
2228 typename VectorType,
2229 typename VectorizedArrayType>
2230 void
2231 LaplaceOperator<dim,
2232 fe_degree,
2233 n_q_points_1d,
2234 n_components,
2235 VectorType,
2236 VectorizedArrayType>::
2237 set_coefficient(
2238 const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient_)
2239 {
2240 scalar_coefficient = scalar_coefficient_;
2241 }
2242
2243
2244
2245 template <int dim,
2246 int fe_degree,
2247 int n_q_points_1d,
2248 int n_components,
2249 typename VectorType,
2250 typename VectorizedArrayType>
2251 std::shared_ptr<Table<2, VectorizedArrayType>>
2252 LaplaceOperator<dim,
2253 fe_degree,
2254 n_q_points_1d,
2255 n_components,
2256 VectorType,
2257 VectorizedArrayType>::get_coefficient()
2258 {
2259 Assert(scalar_coefficient.get(), ExcNotInitialized());
2260 return scalar_coefficient;
2261 }
2262
2263
2264
2265 template <int dim,
2266 int fe_degree,
2267 int n_q_points_1d,
2268 int n_components,
2269 typename VectorType,
2270 typename VectorizedArrayType>
2271 void
2272 LaplaceOperator<dim,
2273 fe_degree,
2274 n_q_points_1d,
2275 n_components,
2276 VectorType,
2277 VectorizedArrayType>::compute_diagonal()
2278 {
2279 using Number =
2280 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2281 Assert((Base<dim, VectorType, VectorizedArrayType>::data.get() != nullptr),
2282 ExcNotInitialized());
2283
2284 this->inverse_diagonal_entries =
2285 std::make_shared<DiagonalMatrix<VectorType>>();
2286 this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
2287 VectorType &inverse_diagonal_vector =
2288 this->inverse_diagonal_entries->get_vector();
2289 VectorType &diagonal_vector = this->diagonal_entries->get_vector();
2290 this->initialize_dof_vector(inverse_diagonal_vector);
2291 this->initialize_dof_vector(diagonal_vector);
2292
2293 this->data->cell_loop(&LaplaceOperator::local_diagonal_cell,
2294 this,
2295 diagonal_vector,
2296 /*unused*/ diagonal_vector);
2297 this->set_constrained_entries_to_one(diagonal_vector);
2298
2299 inverse_diagonal_vector = diagonal_vector;
2300
2301 for (unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
2302 ++i)
2303 if (std::abs(inverse_diagonal_vector.local_element(i)) >
2304 std::sqrt(std::numeric_limits<Number>::epsilon()))
2305 inverse_diagonal_vector.local_element(i) =
2306 1. / inverse_diagonal_vector.local_element(i);
2307 else
2308 inverse_diagonal_vector.local_element(i) = 1.;
2309
2310 // We never need ghost values so don't update them
2311 }
2312
2313
2314
2315 template <int dim,
2316 int fe_degree,
2317 int n_q_points_1d,
2318 int n_components,
2319 typename VectorType,
2320 typename VectorizedArrayType>
2321 void
2322 LaplaceOperator<dim,
2323 fe_degree,
2324 n_q_points_1d,
2325 n_components,
2326 VectorType,
2327 VectorizedArrayType>::apply_add(VectorType &dst,
2328 const VectorType &src) const
2329 {
2330 Base<dim, VectorType, VectorizedArrayType>::data->cell_loop(
2331 &LaplaceOperator::local_apply_cell, this, dst, src);
2332 }
2333
2334 namespace Implementation
2335 {
2336 template <typename VectorizedArrayType>
2337 bool
2338 non_negative(const VectorizedArrayType &n)
2339 {
2340 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2341 if (n[v] < 0.)
2342 return false;
2343
2344 return true;
2345 }
2346 } // namespace Implementation
2347
2348
2349
2350 template <int dim,
2351 int fe_degree,
2352 int n_q_points_1d,
2353 int n_components,
2354 typename VectorType,
2355 typename VectorizedArrayType>
2356 template <int n_components_compute>
2357 void
2358 LaplaceOperator<dim,
2359 fe_degree,
2360 n_q_points_1d,
2361 n_components,
2362 VectorType,
2363 VectorizedArrayType>::
2364 do_operation_on_cell(
2365 FEEvaluation<
2366 dim,
2367 fe_degree,
2368 n_q_points_1d,
2369 n_components_compute,
2370 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2371 VectorizedArrayType> &phi,
2372 const unsigned int cell) const
2373 {
2374 phi.evaluate(EvaluationFlags::gradients);
2375 if (scalar_coefficient.get())
2376 {
2377 Assert(scalar_coefficient->size(1) == 1 ||
2378 scalar_coefficient->size(1) == phi.n_q_points,
2379 ExcMessage("The number of columns in the coefficient table must "
2380 "be either 1 or the number of quadrature points " +
2381 std::to_string(phi.n_q_points) +
2382 ", but the given value was " +
2383 std::to_string(scalar_coefficient->size(1))));
2384 if (scalar_coefficient->size(1) == phi.n_q_points)
2385 for (unsigned int q = 0; q < phi.n_q_points; ++q)
2386 {
2387 Assert(Implementation::non_negative(
2388 (*scalar_coefficient)(cell, q)),
2389 ExcMessage("Coefficient must be non-negative"));
2390 phi.submit_gradient((*scalar_coefficient)(cell, q) *
2391 phi.get_gradient(q),
2392 q);
2393 }
2394 else
2395 {
2396 Assert(Implementation::non_negative((*scalar_coefficient)(cell, 0)),
2397 ExcMessage("Coefficient must be non-negative"));
2398 const VectorizedArrayType coefficient =
2399 (*scalar_coefficient)(cell, 0);
2400 for (unsigned int q = 0; q < phi.n_q_points; ++q)
2401 phi.submit_gradient(coefficient * phi.get_gradient(q), q);
2402 }
2403 }
2404 else
2405 {
2406 for (unsigned int q = 0; q < phi.n_q_points; ++q)
2407 {
2408 phi.submit_gradient(phi.get_gradient(q), q);
2409 }
2410 }
2411 phi.integrate(EvaluationFlags::gradients);
2412 }
2413
2414
2415
2416 template <int dim,
2417 int fe_degree,
2418 int n_q_points_1d,
2419 int n_components,
2420 typename VectorType,
2421 typename VectorizedArrayType>
2422 void
2423 LaplaceOperator<dim,
2424 fe_degree,
2425 n_q_points_1d,
2426 n_components,
2427 VectorType,
2428 VectorizedArrayType>::
2429 local_apply_cell(
2430 const MatrixFree<
2431 dim,
2432 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2433 VectorizedArrayType> &data,
2434 VectorType &dst,
2435 const VectorType &src,
2436 const std::pair<unsigned int, unsigned int> &cell_range) const
2437 {
2438 using Number =
2439 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2440 FEEvaluation<dim,
2441 fe_degree,
2442 n_q_points_1d,
2443 n_components,
2444 Number,
2445 VectorizedArrayType>
2446 phi(data, this->selected_rows[0]);
2447 for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2448 {
2449 phi.reinit(cell);
2450 phi.read_dof_values(src);
2451 do_operation_on_cell(phi, cell);
2452 phi.distribute_local_to_global(dst);
2453 }
2454 }
2455
2456
2457 template <int dim,
2458 int fe_degree,
2459 int n_q_points_1d,
2460 int n_components,
2461 typename VectorType,
2462 typename VectorizedArrayType>
2463 void
2464 LaplaceOperator<dim,
2465 fe_degree,
2466 n_q_points_1d,
2467 n_components,
2468 VectorType,
2469 VectorizedArrayType>::
2470 local_diagonal_cell(
2471 const MatrixFree<
2472 dim,
2473 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2474 VectorizedArrayType> &data,
2475 VectorType &dst,
2476 const VectorType &,
2477 const std::pair<unsigned int, unsigned int> &cell_range) const
2478 {
2479 using Number =
2480 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2481
2482 FEEvaluation<dim, fe_degree, n_q_points_1d, 1, Number, VectorizedArrayType>
2483 eval(data, this->selected_rows[0]);
2484 FEEvaluation<dim,
2485 fe_degree,
2486 n_q_points_1d,
2487 n_components,
2488 Number,
2489 VectorizedArrayType>
2490 eval_vector(data, this->selected_rows[0]);
2491 for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2492 {
2493 eval.reinit(cell);
2494 eval_vector.reinit(cell);
2495 // This function assumes that we have the same result on all
2496 // components, so we only need to go through the columns of one scalar
2497 // component, for which we have created a separate evaluator (attached
2498 // to the first component, but the component does not matter because
2499 // we only use the underlying integrals)
2500 for (unsigned int i = 0; i < eval.dofs_per_cell; ++i)
2501 {
2502 for (unsigned int j = 0; j < eval.dofs_per_cell; ++j)
2503 eval.begin_dof_values()[j] = VectorizedArrayType();
2504 eval.begin_dof_values()[i] = 1.;
2505
2506 do_operation_on_cell(eval, cell);
2507
2508 // We now pick up the value on the diagonal (row i) and broadcast
2509 // it to a second evaluator for all vector components, which we
2510 // will distribute to the result vector afterwards
2511 for (unsigned int c = 0; c < n_components; ++c)
2512 eval_vector
2513 .begin_dof_values()[i + c * eval_vector.dofs_per_component] =
2514 eval.begin_dof_values()[i];
2515 }
2516 eval_vector.distribute_local_to_global(dst);
2517 }
2518 }
2519
2520
2521} // end of namespace MatrixFreeOperators
2522
2523
2524DEAL_II_NAMESPACE_CLOSE
2525
2526#endif
push_back
***mech_lbc_system increment_interpolation_handlers push_back(scale_z_handler)
make_array_view
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:1002
AlignedVector::size
size_type size() const
FEEvaluationData::begin_dof_values
const Number * begin_dof_values() const
FEEvaluation::reinit
void reinit(const unsigned int cell_batch_index)
FEEvaluation::dofs_per_cell
const unsigned int dofs_per_cell
Definition fe_evaluation.h:1774
IndexSet::index_within_set
size_type index_within_set(const size_type global_index) const
Definition index_set.h:1977
IndexSet::is_element
bool is_element(const size_type index) const
Definition index_set.h:1877
MatrixFreeOperators::Base::~Base
virtual ~Base() override=default
MatrixFreeOperators::Base::selected_rows
std::vector< unsigned int > selected_rows
Definition operators.h:451
MatrixFreeOperators::Base::Tvmult_add
void Tvmult_add(VectorType &dst, const VectorType &src) const
Definition operators.h:1511
MatrixFreeOperators::Base::set_constrained_entries_to_one
void set_constrained_entries_to_one(VectorType &dst) const
Definition operators.h:1468
MatrixFreeOperators::Base::compute_diagonal
virtual void compute_diagonal()=0
MatrixFreeOperators::Base::vmult_add
void vmult_add(VectorType &dst, const VectorType &src) const
Definition operators.h:1500
MatrixFreeOperators::Base::vmult_interface_down
void vmult_interface_down(VectorType &dst, const VectorType &src) const
Definition operators.h:1648
MatrixFreeOperators::Base::n
size_type n() const
Definition operators.h:1266
MatrixFreeOperators::Base::preprocess_constraints
void preprocess_constraints(VectorType &dst, const VectorType &src) const
Definition operators.h:1564
MatrixFreeOperators::Base::mult_add
void mult_add(VectorType &dst, const VectorType &src, const bool transpose) const
Definition operators.h:1594
MatrixFreeOperators::Base::edge_constrained_indices
std::vector< std::vector< unsigned int > > edge_constrained_indices
Definition operators.h:463
MatrixFreeOperators::Base::get_matrix_diagonal
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal() const
Definition operators.h:1793
MatrixFreeOperators::Base::Tvmult
void Tvmult(VectorType &dst, const VectorType &src) const
Definition operators.h:1743
MatrixFreeOperators::Base::edge_constrained_values
std::vector< std::vector< std::pair< value_type, value_type > > > edge_constrained_values
Definition operators.h:469
MatrixFreeOperators::Base::m
size_type m() const
Definition operators.h:1252
MatrixFreeOperators::Base::initialize
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 >())
Definition operators.h:1329
MatrixFreeOperators::Base::selected_columns
std::vector< unsigned int > selected_columns
Definition operators.h:457
MatrixFreeOperators::Base::get_matrix_free
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > get_matrix_free() const
Definition operators.h:1771
MatrixFreeOperators::Base::initialize
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 >())
Definition operators.h:1398
MatrixFreeOperators::Base::get_matrix_diagonal_inverse
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal_inverse() const
Definition operators.h:1780
MatrixFreeOperators::Base::initialize
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 >())
Definition operators.h:1382
MatrixFreeOperators::Base::initialize_dof_vector
void initialize_dof_vector(VectorType &vec) const
Definition operators.h:1304
MatrixFreeOperators::Base::diagonal_entries
std::shared_ptr< DiagonalMatrix< VectorType > > diagonal_entries
Definition operators.h:439
MatrixFreeOperators::Base::Tapply_add
virtual void Tapply_add(VectorType &dst, const VectorType &src) const
Definition operators.h:1804
MatrixFreeOperators::Base::data
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > data
Definition operators.h:433
MatrixFreeOperators::Base::adjust_ghost_range_if_necessary
void adjust_ghost_range_if_necessary(const VectorType &vec, const bool is_row) const
Definition operators.h:1522
MatrixFreeOperators::Base::memory_consumption
virtual std::size_t memory_consumption() const
Definition operators.h:1757
MatrixFreeOperators::Base::inverse_diagonal_entries
std::shared_ptr< DiagonalMatrix< VectorType > > inverse_diagonal_entries
Definition operators.h:445
MatrixFreeOperators::Base::el
value_type el(const unsigned int row, const unsigned int col) const
Definition operators.h:1290
MatrixFreeOperators::Base::apply_add
virtual void apply_add(VectorType &dst, const VectorType &src) const =0
MatrixFreeOperators::Base::size_type
typename VectorType::size_type size_type
Definition operators.h:196
MatrixFreeOperators::Base::precondition_Jacobi
void precondition_Jacobi(VectorType &dst, const VectorType &src, const value_type omega) const
Definition operators.h:1815
MatrixFreeOperators::Base::value_type
typename VectorType::value_type value_type
Definition operators.h:191
MatrixFreeOperators::Base::vmult
void vmult(VectorType &dst, const VectorType &src) const
Definition operators.h:1487
MatrixFreeOperators::Base::vmult_interface_up
void vmult_interface_up(VectorType &dst, const VectorType &src) const
Definition operators.h:1702
MatrixFreeOperators::Base::postprocess_constraints
void postprocess_constraints(VectorType &dst, const VectorType &src) const
Definition operators.h:1614
MatrixFreeOperators::CellwiseInverseMassMatrix::fe_eval
const FEEvaluationBase< dim, n_components, Number, false, VectorizedArrayType > & fe_eval
Definition operators.h:740
MatrixFreeOperators::CellwiseInverseMassMatrix::fill_inverse_JxW_values
void fill_inverse_JxW_values(AlignedVector< VectorizedArrayType > &inverse_jxw) const
Definition operators.h:1078
MatrixFreeOperators::CellwiseInverseMassMatrix::transform_from_q_points_to_basis
void transform_from_q_points_to_basis(const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition operators.h:1217
MatrixFreeOperators::CellwiseInverseMassMatrix::apply
void apply(const AlignedVector< VectorizedArrayType > &inverse_coefficient, const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition operators.h:1138
MatrixFreeOperators::LaplaceOperator::get_coefficient
std::shared_ptr< Table< 2, VectorizedArrayType > > get_coefficient()
Definition operators.h:2257
MatrixFreeOperators::LaplaceOperator::value_type
typename Base< dim, VectorType, VectorizedArrayType >::value_type value_type
Definition operators.h:899
MatrixFreeOperators::LaplaceOperator::compute_diagonal
virtual void compute_diagonal() override
Definition operators.h:2277
MatrixFreeOperators::LaplaceOperator::scalar_coefficient
std::shared_ptr< Table< 2, VectorizedArrayType > > scalar_coefficient
Definition operators.h:1036
MatrixFreeOperators::LaplaceOperator::local_diagonal_cell
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
Definition operators.h:2470
MatrixFreeOperators::LaplaceOperator::local_apply_cell
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
Definition operators.h:2429
MatrixFreeOperators::LaplaceOperator::apply_add
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition operators.h:2327
MatrixFreeOperators::LaplaceOperator::size_type
typename Base< dim, VectorType, VectorizedArrayType >::size_type size_type
Definition operators.h:905
MatrixFreeOperators::LaplaceOperator::set_coefficient
void set_coefficient(const std::shared_ptr< Table< 2, VectorizedArrayType > > &scalar_coefficient)
Definition operators.h:2237
MatrixFreeOperators::LaplaceOperator::do_operation_on_cell
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
MatrixFreeOperators::LaplaceOperator::clear
virtual void clear() override
Definition operators.h:2216
MatrixFreeOperators::MGInterfaceOperator::value_type
typename OperatorType::value_type value_type
Definition operators.h:541
MatrixFreeOperators::MGInterfaceOperator::vmult
void vmult(VectorType &dst, const VectorType &src) const
Definition operators.h:1860
MatrixFreeOperators::MGInterfaceOperator::mf_base_operator
ObserverPointer< const OperatorType > mf_base_operator
Definition operators.h:591
MatrixFreeOperators::MGInterfaceOperator::initialize
void initialize(const OperatorType &operator_in)
Definition operators.h:1850
MatrixFreeOperators::MGInterfaceOperator::Tvmult
void Tvmult(VectorType &dst, const VectorType &src) const
Definition operators.h:1880
MatrixFreeOperators::MGInterfaceOperator::initialize_dof_vector
void initialize_dof_vector(VectorType &vec) const
Definition operators.h:1900
MatrixFreeOperators::MGInterfaceOperator::size_type
typename OperatorType::size_type size_type
Definition operators.h:546
MatrixFreeOperators::MassOperator::lumped_diagonal_entries
std::shared_ptr< DiagonalMatrix< VectorType > > lumped_diagonal_entries
Definition operators.h:839
MatrixFreeOperators::MassOperator::size_type
typename Base< dim, VectorType, VectorizedArrayType >::size_type size_type
Definition operators.h:772
MatrixFreeOperators::MassOperator::apply_add
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition operators.h:2132
MatrixFreeOperators::MassOperator::get_matrix_lumped_diagonal
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_lumped_diagonal() const
Definition operators.h:2110
MatrixFreeOperators::MassOperator::get_matrix_lumped_diagonal_inverse
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_lumped_diagonal_inverse() const
Definition operators.h:2088
MatrixFreeOperators::MassOperator::local_apply_cell
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
Definition operators.h:2154
MatrixFreeOperators::MassOperator::compute_diagonal
virtual void compute_diagonal() override
Definition operators.h:1948
MatrixFreeOperators::MassOperator::value_type
typename Base< dim, VectorType, VectorizedArrayType >::value_type value_type
Definition operators.h:766
MatrixFreeOperators::MassOperator::inverse_lumped_diagonal_entries
std::shared_ptr< DiagonalMatrix< VectorType > > inverse_lumped_diagonal_entries
Definition operators.h:845
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
running_in_debug_mode
constexpr bool running_in_debug_mode()
Definition config.h:76
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition derivative_form.h:629
diagonal_matrix.h
enable_observer_pointer.h
fe_evaluation.h
level
unsigned int level
Definition grid_out.cc:4642
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exception_macros.h:543
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exception_macros.h:829
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exception_macros.h:901
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcNotInitialized
static ::ExceptionBase & ExcNotInitialized()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition exception_macros.h:647
la_parallel_vector.h
matrix_free.h
mg_constrained_dofs.h
data
std::vector< index_type > data
Definition mpi.cc:734
size
std::size_t size
Definition mpi.cc:733
locally_owned_size
types::global_dof_index locally_owned_size
Definition mpi.cc:821
MatrixFreeOperators::BlockHelper::subblock
std::enable_if_t< IsBlockVector< VectorType >::value, typename VectorType::BlockType & > subblock(VectorType &vector, unsigned int block_no)
Definition operators.h:62
MatrixFreeOperators::BlockHelper::collect_sizes
std::enable_if_t< IsBlockVector< VectorType >::value, void > collect_sizes(VectorType &vector)
Definition operators.h:93
MatrixFreeOperators::BlockHelper::n_blocks
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition operators.h:47
MatrixFreeOperators::Implementation::non_negative
bool non_negative(const VectorizedArrayType &n)
Definition operators.h:2338
MatrixFreeTools::compute_diagonal
void compute_diagonal(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, VectorType &diagonal_global, const std::function< void(FEEvaluation< dim, fe_degree, n_q_points_1d, n_components, Number, VectorizedArrayType > &)> &cell_operation, const unsigned int dof_handler_index=0, const unsigned int quadrature_index=0, const unsigned int first_selected_component=0, const unsigned int first_vector_component=0)
Utilities::MPI::max
T max(const T &t, const MPI_Comm mpi_communicator)
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
numbers::invalid_unsigned_int
constexpr unsigned int invalid_unsigned_int
Definition types.h:228
std::sqrt
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6505
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6563
internal::CellwiseInverseMassFactory::apply
static void apply(const unsigned int n_components, const FEEvaluationData< dim, Number, false > &fe_eval, const Number *in_array, Number *out_array)
internal::CellwiseInverseMassFactory::transform_from_q_points_to_basis
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)
Generated by doxygen 1.9.8