Reference documentation for deal.II version 9.5.0
\(\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// Copyright (C) 2011 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
17#ifndef dealii_matrix_free_operators_h
18#define dealii_matrix_free_operators_h
19
20
21#include <deal.II/base/config.h>
22
23#include <deal.II/base/exceptions.h>
24#include <deal.II/base/subscriptor.h>
25#include <deal.II/base/vectorization.h>
26
27#include <deal.II/lac/diagonal_matrix.h>
28#include <deal.II/lac/la_parallel_vector.h>
29
30#include <deal.II/matrix_free/fe_evaluation.h>
31#include <deal.II/matrix_free/matrix_free.h>
32#include <deal.II/matrix_free/tools.h>
33
34#include <deal.II/multigrid/mg_constrained_dofs.h>
35
36#include <limits>
37
38DEAL_II_NAMESPACE_OPEN
39
40
41namespace MatrixFreeOperators
42{
43 namespace BlockHelper
44 {
45 // workaround for unifying non-block vector and block vector implementations
46 // a non-block vector has one block and the only subblock is the vector
47 // itself
48 template <typename VectorType>
49 std::enable_if_t<IsBlockVector<VectorType>::value, unsigned int>
50 n_blocks(const VectorType &vector)
51 {
52 return vector.n_blocks();
53 }
54
55 template <typename VectorType>
56 std::enable_if_t<!IsBlockVector<VectorType>::value, unsigned int>
57 n_blocks(const VectorType &)
58 {
59 return 1;
60 }
61
62 template <typename VectorType>
63 std::enable_if_t<IsBlockVector<VectorType>::value,
64 typename VectorType::BlockType &>
65 subblock(VectorType &vector, unsigned int block_no)
66 {
67 AssertIndexRange(block_no, vector.n_blocks());
68 return vector.block(block_no);
69 }
70
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)
75 {
76 AssertIndexRange(block_no, vector.n_blocks());
77 return vector.block(block_no);
78 }
79
80 template <typename VectorType>
81 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType &>
82 subblock(VectorType &vector, unsigned int)
83 {
84 return vector;
85 }
86
87 template <typename VectorType>
88 std::enable_if_t<!IsBlockVector<VectorType>::value, const VectorType &>
89 subblock(const VectorType &vector, unsigned int)
90 {
91 return vector;
92 }
93
94 template <typename VectorType>
95 std::enable_if_t<IsBlockVector<VectorType>::value, void>
96 collect_sizes(VectorType &vector)
97 {
98 vector.collect_sizes();
99 }
100
101 template <typename VectorType>
102 std::enable_if_t<!IsBlockVector<VectorType>::value, void>
103 collect_sizes(const VectorType &)
104 {}
105 } // namespace BlockHelper
106
184 template <int dim,
185 typename VectorType = LinearAlgebra::distributed::Vector<double>,
186 typename VectorizedArrayType =
187 VectorizedArray<typename VectorType::value_type>>
188 class Base : public Subscriptor
189 {
190 public:
194 using value_type = typename VectorType::value_type;
195
199 using size_type = typename VectorType::size_type;
200
204 Base();
205
209 virtual ~Base() override = default;
210
215 virtual void
216 clear();
217
234 void
235 initialize(std::shared_ptr<
236 const MatrixFree<dim, value_type, VectorizedArrayType>> data,
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>());
241
255 void
256 initialize(std::shared_ptr<
257 const MatrixFree<dim, value_type, VectorizedArrayType>> data,
258 const MGConstrainedDoFs & mg_constrained_dofs,
259 const unsigned int level,
260 const std::vector<unsigned int> &selected_row_blocks =
261 std::vector<unsigned int>());
262
277 void
278 initialize(std::shared_ptr<
279 const MatrixFree<dim, value_type, VectorizedArrayType>> data_,
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>());
284
288 size_type
289 m() const;
290
294 size_type
295 n() const;
296
300 void
301 vmult_interface_down(VectorType &dst, const VectorType &src) const;
302
306 void
307 vmult_interface_up(VectorType &dst, const VectorType &src) const;
308
312 void
313 vmult(VectorType &dst, const VectorType &src) const;
314
318 void
319 Tvmult(VectorType &dst, const VectorType &src) const;
320
324 void
325 vmult_add(VectorType &dst, const VectorType &src) const;
326
330 void
331 Tvmult_add(VectorType &dst, const VectorType &src) const;
332
337 value_type
338 el(const unsigned int row, const unsigned int col) const;
339
344 virtual std::size_t
345 memory_consumption() const;
346
350 void
351 initialize_dof_vector(VectorType &vec) const;
352
365 virtual void
366 compute_diagonal() = 0;
367
371 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
372 get_matrix_free() const;
373
377 const std::shared_ptr<DiagonalMatrix<VectorType>> &
378 get_matrix_diagonal_inverse() const;
379
383 const std::shared_ptr<DiagonalMatrix<VectorType>> &
384 get_matrix_diagonal() const;
385
391 void
392 precondition_Jacobi(VectorType & dst,
393 const VectorType &src,
394 const value_type omega) const;
395
396 protected:
401 void
402 preprocess_constraints(VectorType &dst, const VectorType &src) const;
403
408 void
409 postprocess_constraints(VectorType &dst, const VectorType &src) const;
410
415 void
416 set_constrained_entries_to_one(VectorType &dst) const;
417
421 virtual void
422 apply_add(VectorType &dst, const VectorType &src) const = 0;
423
429 virtual void
430 Tapply_add(VectorType &dst, const VectorType &src) const;
431
435 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
436 data;
437
442 std::shared_ptr<DiagonalMatrix<VectorType>> diagonal_entries;
443
448 std::shared_ptr<DiagonalMatrix<VectorType>> inverse_diagonal_entries;
449
454 std::vector<unsigned int> selected_rows;
455
460 std::vector<unsigned int> selected_columns;
461
462 private:
466 std::vector<std::vector<unsigned int>> edge_constrained_indices;
467
471 mutable std::vector<std::vector<std::pair<value_type, value_type>>>
472 edge_constrained_values;
473
478 bool have_interface_matrices;
479
484 void
485 mult_add(VectorType & dst,
486 const VectorType &src,
487 const bool transpose) const;
488
496 void
497 adjust_ghost_range_if_necessary(const VectorType &vec,
498 const bool is_row) const;
499 };
500
501
502
537 template <typename OperatorType>
538 class MGInterfaceOperator : public Subscriptor
539 {
540 public:
544 using value_type = typename OperatorType::value_type;
545
549 using size_type = typename OperatorType::size_type;
550
554 MGInterfaceOperator();
555
559 void
560 clear();
561
565 void
566 initialize(const OperatorType &operator_in);
567
571 template <typename VectorType>
572 void
573 vmult(VectorType &dst, const VectorType &src) const;
574
578 template <typename VectorType>
579 void
580 Tvmult(VectorType &dst, const VectorType &src) const;
581
585 template <typename VectorType>
586 void
587 initialize_dof_vector(VectorType &vec) const;
588
589
590 private:
594 SmartPointer<const OperatorType> mf_base_operator;
595 };
596
597
598
617 template <int dim,
618 int fe_degree,
619 int n_components = 1,
620 typename Number = double,
621 typename VectorizedArrayType = VectorizedArray<Number>>
622 class CellwiseInverseMassMatrix
623 {
624 static_assert(
625 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
626 "Type of Number and of VectorizedArrayType do not match.");
627
628 public:
633 CellwiseInverseMassMatrix(
634 const FEEvaluationBase<dim,
635 n_components,
636 Number,
637 false,
638 VectorizedArrayType> &fe_eval);
639
649 void
650 apply(const AlignedVector<VectorizedArrayType> &inverse_coefficient,
651 const unsigned int n_actual_components,
652 const VectorizedArrayType * in_array,
653 VectorizedArrayType * out_array) const;
654
666 void
667 apply(const VectorizedArrayType *in_array,
668 VectorizedArrayType * out_array) const;
669
682 void
683 apply(const AlignedVector<Tensor<2, n_components, VectorizedArrayType>>
684 & inverse_dyadic_coefficients,
685 const VectorizedArrayType *in_array,
686 VectorizedArrayType * out_array) const;
687
721 void
722 transform_from_q_points_to_basis(const unsigned int n_actual_components,
723 const VectorizedArrayType *in_array,
724 VectorizedArrayType *out_array) const;
725
731 void
732 fill_inverse_JxW_values(
733 AlignedVector<VectorizedArrayType> &inverse_jxw) const;
734
735 private:
739 const FEEvaluationBase<dim,
740 n_components,
741 Number,
742 false,
743 VectorizedArrayType> &fe_eval;
744 };
745
746
747
755 template <int dim,
756 int fe_degree,
757 int n_q_points_1d = fe_degree + 1,
758 int n_components = 1,
759 typename VectorType = LinearAlgebra::distributed::Vector<double>,
760 typename VectorizedArrayType =
761 VectorizedArray<typename VectorType::value_type>>
762 class MassOperator : public Base<dim, VectorType, VectorizedArrayType>
763 {
764 public:
768 using value_type =
769 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
770
774 using size_type =
775 typename Base<dim, VectorType, VectorizedArrayType>::size_type;
776
780 MassOperator();
781
785 virtual void
786 compute_diagonal() override;
787
804 void
805 compute_lumped_diagonal();
806
810 const std::shared_ptr<DiagonalMatrix<VectorType>> &
811 get_matrix_lumped_diagonal() const;
812
816 const std::shared_ptr<DiagonalMatrix<VectorType>> &
817 get_matrix_lumped_diagonal_inverse() const;
818
819 private:
825 virtual void
826 apply_add(VectorType &dst, const VectorType &src) const override;
827
831 void
832 local_apply_cell(
833 const MatrixFree<dim, value_type, VectorizedArrayType> &data,
834 VectorType & dst,
835 const VectorType & src,
836 const std::pair<unsigned int, unsigned int> & cell_range) const;
837
842 std::shared_ptr<DiagonalMatrix<VectorType>> lumped_diagonal_entries;
843
848 std::shared_ptr<DiagonalMatrix<VectorType>> inverse_lumped_diagonal_entries;
849 };
850
851
852
863 template <int dim,
864 int fe_degree,
865 int n_q_points_1d = fe_degree + 1,
866 int n_components = 1,
867 typename VectorType = LinearAlgebra::distributed::Vector<double>,
868 typename VectorizedArrayType =
869 VectorizedArray<typename VectorType::value_type>>
870 class LaplaceOperator : public Base<dim, VectorType, VectorizedArrayType>
871 {
872 public:
876 using value_type =
877 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
878
882 using size_type =
883 typename Base<dim, VectorType, VectorizedArrayType>::size_type;
884
888 LaplaceOperator();
889
896 virtual void
897 compute_diagonal() override;
898
949 void
950 set_coefficient(
951 const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient);
952
957 virtual void
958 clear() override;
959
966 std::shared_ptr<Table<2, VectorizedArrayType>>
967 get_coefficient();
968
969 private:
975 virtual void
976 apply_add(VectorType &dst, const VectorType &src) const override;
977
981 void
982 local_apply_cell(
983 const MatrixFree<dim, value_type, VectorizedArrayType> &data,
984 VectorType & dst,
985 const VectorType & src,
986 const std::pair<unsigned int, unsigned int> & cell_range) const;
987
991 void
992 local_diagonal_cell(
993 const MatrixFree<dim, value_type, VectorizedArrayType> &data,
994 VectorType & dst,
995 const VectorType &,
996 const std::pair<unsigned int, unsigned int> &cell_range) const;
997
1001 template <int n_components_compute>
1002 void
1003 do_operation_on_cell(FEEvaluation<dim,
1004 fe_degree,
1005 n_q_points_1d,
1006 n_components_compute,
1007 value_type,
1008 VectorizedArrayType> &phi,
1009 const unsigned int cell) const;
1010
1014 std::shared_ptr<Table<2, VectorizedArrayType>> scalar_coefficient;
1015 };
1016
1017
1018
1019 // ------------------------------------ inline functions ---------------------
1020
1021 template <int dim,
1022 int fe_degree,
1023 int n_components,
1024 typename Number,
1025 typename VectorizedArrayType>
1026 inline CellwiseInverseMassMatrix<dim,
1027 fe_degree,
1028 n_components,
1029 Number,
1030 VectorizedArrayType>::
1031 CellwiseInverseMassMatrix(
1032 const FEEvaluationBase<dim,
1033 n_components,
1034 Number,
1035 false,
1036 VectorizedArrayType> &fe_eval)
1037 : fe_eval(fe_eval)
1038 {}
1039
1040
1041
1042 template <int dim,
1043 int fe_degree,
1044 int n_components,
1045 typename Number,
1046 typename VectorizedArrayType>
1047 inline void
1048 CellwiseInverseMassMatrix<dim,
1049 fe_degree,
1050 n_components,
1051 Number,
1052 VectorizedArrayType>::
1053 fill_inverse_JxW_values(
1054 AlignedVector<VectorizedArrayType> &inverse_jxw) const
1055 {
1056 const unsigned int dofs_per_component_on_cell =
1057 (fe_degree > -1) ?
1058 Utilities::pow(fe_degree + 1, dim) :
1059 Utilities::pow(fe_eval.get_shape_info().data.front().fe_degree + 1,
1060 dim);
1061
1062 Assert(inverse_jxw.size() > 0 &&
1063 inverse_jxw.size() % dofs_per_component_on_cell == 0,
1064 ExcMessage(
1065 "Expected diagonal to be a multiple of scalar dof per cells"));
1066
1067 // compute values for the first component
1068 for (unsigned int q = 0; q < dofs_per_component_on_cell; ++q)
1069 inverse_jxw[q] = 1. / fe_eval.JxW(q);
1070 // copy values to rest of vector
1071 for (unsigned int q = dofs_per_component_on_cell; q < inverse_jxw.size();)
1072 for (unsigned int i = 0; i < dofs_per_component_on_cell; ++i, ++q)
1073 inverse_jxw[q] = inverse_jxw[i];
1074 }
1075
1076
1077
1078 template <int dim,
1079 int fe_degree,
1080 int n_components,
1081 typename Number,
1082 typename VectorizedArrayType>
1083 inline void
1084 CellwiseInverseMassMatrix<
1085 dim,
1086 fe_degree,
1087 n_components,
1088 Number,
1089 VectorizedArrayType>::apply(const VectorizedArrayType *in_array,
1090 VectorizedArrayType * out_array) const
1091 {
1092 if (fe_degree > -1)
1093 internal::CellwiseInverseMassMatrixImplBasic<dim, VectorizedArrayType>::
1094 template run<fe_degree>(n_components, fe_eval, in_array, out_array);
1095 else
1096 internal::CellwiseInverseMassFactory<dim, VectorizedArrayType>::apply(
1097 n_components, fe_eval, in_array, out_array);
1098 }
1099
1100
1101
1102 template <int dim,
1103 int fe_degree,
1104 int n_components,
1105 typename Number,
1106 typename VectorizedArrayType>
1107 inline void
1108 CellwiseInverseMassMatrix<dim,
1109 fe_degree,
1110 n_components,
1111 Number,
1112 VectorizedArrayType>::
1113 apply(const AlignedVector<VectorizedArrayType> &inverse_coefficients,
1114 const unsigned int n_actual_components,
1115 const VectorizedArrayType * in_array,
1116 VectorizedArrayType * out_array) const
1117 {
1118 if (fe_degree > -1)
1119 internal::CellwiseInverseMassMatrixImplFlexible<
1120 dim,
1121 VectorizedArrayType>::template run<fe_degree>(n_actual_components,
1122 fe_eval,
1123 make_array_view(
1124 inverse_coefficients),
1125 false,
1126 in_array,
1127 out_array);
1128 else
1129 internal::CellwiseInverseMassFactory<dim, VectorizedArrayType>::apply(
1130 n_actual_components,
1131 fe_eval,
1132 make_array_view(inverse_coefficients),
1133 false,
1134 in_array,
1135 out_array);
1136 }
1137
1138 template <int dim,
1139 int fe_degree,
1140 int n_components,
1141 typename Number,
1142 typename VectorizedArrayType>
1143 inline void
1144 CellwiseInverseMassMatrix<dim,
1145 fe_degree,
1146 n_components,
1147 Number,
1148 VectorizedArrayType>::
1149 apply(const AlignedVector<Tensor<2, n_components, VectorizedArrayType>>
1150 & inverse_dyadic_coefficients,
1151 const VectorizedArrayType *in_array,
1152 VectorizedArrayType * out_array) const
1153 {
1154 const unsigned int unrolled_size =
1155 inverse_dyadic_coefficients.size() * (n_components * n_components);
1156
1157 if (fe_degree > -1)
1158 internal::CellwiseInverseMassMatrixImplFlexible<dim,
1159 VectorizedArrayType>::
1160 template run<fe_degree>(n_components,
1161 fe_eval,
1162 ArrayView<const VectorizedArrayType>(
1163 &inverse_dyadic_coefficients[0][0][0],
1164 unrolled_size),
1165 true,
1166 in_array,
1167 out_array);
1168 else
1169 internal::CellwiseInverseMassFactory<dim, VectorizedArrayType>::apply(
1170 n_components,
1171 fe_eval,
1172 ArrayView<const VectorizedArrayType>(
1173 &inverse_dyadic_coefficients[0][0][0], unrolled_size),
1174 true,
1175 in_array,
1176 out_array);
1177 }
1178
1179
1180
1181 template <int dim,
1182 int fe_degree,
1183 int n_components,
1184 typename Number,
1185 typename VectorizedArrayType>
1186 inline void
1187 CellwiseInverseMassMatrix<dim,
1188 fe_degree,
1189 n_components,
1190 Number,
1191 VectorizedArrayType>::
1192 transform_from_q_points_to_basis(const unsigned int n_actual_components,
1193 const VectorizedArrayType *in_array,
1194 VectorizedArrayType * out_array) const
1195 {
1196 const auto n_q_points_1d = fe_eval.get_shape_info().data[0].n_q_points_1d;
1197
1198 if (fe_degree > -1 && (fe_degree + 1 == n_q_points_1d))
1199 internal::CellwiseInverseMassMatrixImplTransformFromQPoints<
1200 dim,
1201 VectorizedArrayType>::template run<fe_degree,
1202 fe_degree + 1>(n_actual_components,
1203 fe_eval,
1204 in_array,
1205 out_array);
1206 else
1207 internal::CellwiseInverseMassFactory<dim, VectorizedArrayType>::
1208 transform_from_q_points_to_basis(n_actual_components,
1209 fe_eval,
1210 in_array,
1211 out_array);
1212 }
1213
1214
1215
1216 //----------------- Base operator -----------------------------
1217 template <int dim, typename VectorType, typename VectorizedArrayType>
1218 Base<dim, VectorType, VectorizedArrayType>::Base()
1219 : Subscriptor()
1220 , have_interface_matrices(false)
1221 {}
1222
1223
1224
1225 template <int dim, typename VectorType, typename VectorizedArrayType>
1226 typename Base<dim, VectorType, VectorizedArrayType>::size_type
1227 Base<dim, VectorType, VectorizedArrayType>::m() const
1228 {
1229 Assert(data.get() != nullptr, ExcNotInitialized());
1230 typename Base<dim, VectorType, VectorizedArrayType>::size_type total_size =
1231 0;
1232 for (const unsigned int selected_row : selected_rows)
1233 total_size += data->get_vector_partitioner(selected_row)->size();
1234 return total_size;
1235 }
1236
1237
1238
1239 template <int dim, typename VectorType, typename VectorizedArrayType>
1240 typename Base<dim, VectorType, VectorizedArrayType>::size_type
1241 Base<dim, VectorType, VectorizedArrayType>::n() const
1242 {
1243 Assert(data.get() != nullptr, ExcNotInitialized());
1244 typename Base<dim, VectorType, VectorizedArrayType>::size_type total_size =
1245 0;
1246 for (const unsigned int selected_column : selected_columns)
1247 total_size += data->get_vector_partitioner(selected_column)->size();
1248 return total_size;
1249 }
1250
1251
1252
1253 template <int dim, typename VectorType, typename VectorizedArrayType>
1254 void
1255 Base<dim, VectorType, VectorizedArrayType>::clear()
1256 {
1257 data.reset();
1258 inverse_diagonal_entries.reset();
1259 }
1260
1261
1262
1263 template <int dim, typename VectorType, typename VectorizedArrayType>
1264 typename Base<dim, VectorType, VectorizedArrayType>::value_type
1265 Base<dim, VectorType, VectorizedArrayType>::el(const unsigned int row,
1266 const unsigned int col) const
1267 {
1268 (void)col;
1269 Assert(row == col, ExcNotImplemented());
1270 Assert(inverse_diagonal_entries.get() != nullptr &&
1271 inverse_diagonal_entries->m() > 0,
1272 ExcNotInitialized());
1273 return 1.0 / (*inverse_diagonal_entries)(row, row);
1274 }
1275
1276
1277
1278 template <int dim, typename VectorType, typename VectorizedArrayType>
1279 void
1280 Base<dim, VectorType, VectorizedArrayType>::initialize_dof_vector(
1281 VectorType &vec) const
1282 {
1283 Assert(data.get() != nullptr, ExcNotInitialized());
1284 AssertDimension(BlockHelper::n_blocks(vec), selected_rows.size());
1285 for (unsigned int i = 0; i < BlockHelper::n_blocks(vec); ++i)
1286 {
1287 const unsigned int index = selected_rows[i];
1288 if (!BlockHelper::subblock(vec, index)
1289 .partitioners_are_compatible(
1290 *data->get_dof_info(index).vector_partitioner))
1291 data->initialize_dof_vector(BlockHelper::subblock(vec, index), index);
1292
1293 Assert(BlockHelper::subblock(vec, index)
1294 .partitioners_are_globally_compatible(
1295 *data->get_dof_info(index).vector_partitioner),
1296 ExcInternalError());
1297 }
1298 BlockHelper::collect_sizes(vec);
1299 }
1300
1301
1302
1303 template <int dim, typename VectorType, typename VectorizedArrayType>
1304 void
1305 Base<dim, VectorType, VectorizedArrayType>::initialize(
1306 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
1307 data_,
1308 const std::vector<unsigned int> &given_row_selection,
1309 const std::vector<unsigned int> &given_column_selection)
1310 {
1311 data = data_;
1312
1313 selected_rows.clear();
1314 selected_columns.clear();
1315 if (given_row_selection.empty())
1316 for (unsigned int i = 0; i < data_->n_components(); ++i)
1317 selected_rows.push_back(i);
1318 else
1319 {
1320 for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1321 {
1322 AssertIndexRange(given_row_selection[i], data_->n_components());
1323 for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1324 if (j != i)
1325 Assert(given_row_selection[j] != given_row_selection[i],
1326 ExcMessage("Given row indices must be unique"));
1327
1328 selected_rows.push_back(given_row_selection[i]);
1329 }
1330 }
1331 if (given_column_selection.size() == 0)
1332 selected_columns = selected_rows;
1333 else
1334 {
1335 for (unsigned int i = 0; i < given_column_selection.size(); ++i)
1336 {
1337 AssertIndexRange(given_column_selection[i], data_->n_components());
1338 for (unsigned int j = 0; j < given_column_selection.size(); ++j)
1339 if (j != i)
1340 Assert(given_column_selection[j] != given_column_selection[i],
1341 ExcMessage("Given column indices must be unique"));
1342
1343 selected_columns.push_back(given_column_selection[i]);
1344 }
1345 }
1346
1347 edge_constrained_indices.clear();
1348 edge_constrained_indices.resize(selected_rows.size());
1349 edge_constrained_values.clear();
1350 edge_constrained_values.resize(selected_rows.size());
1351 have_interface_matrices = false;
1352 }
1353
1354
1355
1356 template <int dim, typename VectorType, typename VectorizedArrayType>
1357 void
1358 Base<dim, VectorType, VectorizedArrayType>::initialize(
1359 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
1360 data_,
1361 const MGConstrainedDoFs & mg_constrained_dofs,
1362 const unsigned int level,
1363 const std::vector<unsigned int> &given_row_selection)
1364 {
1365 std::vector<MGConstrainedDoFs> mg_constrained_dofs_vector(
1366 1, mg_constrained_dofs);
1367 initialize(data_, mg_constrained_dofs_vector, level, given_row_selection);
1368 }
1369
1370
1371
1372 template <int dim, typename VectorType, typename VectorizedArrayType>
1373 void
1374 Base<dim, VectorType, VectorizedArrayType>::initialize(
1375 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
1376 data_,
1377 const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
1378 const unsigned int level,
1379 const std::vector<unsigned int> & given_row_selection)
1380 {
1381 AssertThrow(level != numbers::invalid_unsigned_int,
1382 ExcMessage("level is not set"));
1383
1384 selected_rows.clear();
1385 selected_columns.clear();
1386 if (given_row_selection.empty())
1387 for (unsigned int i = 0; i < data_->n_components(); ++i)
1388 selected_rows.push_back(i);
1389 else
1390 {
1391 for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1392 {
1393 AssertIndexRange(given_row_selection[i], data_->n_components());
1394 for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1395 if (j != i)
1396 Assert(given_row_selection[j] != given_row_selection[i],
1397 ExcMessage("Given row indices must be unique"));
1398
1399 selected_rows.push_back(given_row_selection[i]);
1400 }
1401 }
1402 selected_columns = selected_rows;
1403
1404 AssertDimension(mg_constrained_dofs.size(), selected_rows.size());
1405 edge_constrained_indices.clear();
1406 edge_constrained_indices.resize(selected_rows.size());
1407 edge_constrained_values.clear();
1408 edge_constrained_values.resize(selected_rows.size());
1409
1410 data = data_;
1411
1412 for (unsigned int j = 0; j < selected_rows.size(); ++j)
1413 {
1414 if (data_->n_cell_batches() > 0)
1415 {
1416 AssertDimension(level, data_->get_cell_iterator(0, 0, j)->level());
1417 }
1418
1419 // setup edge_constrained indices
1420 const std::vector<types::global_dof_index> interface_indices =
1421 mg_constrained_dofs[j]
1422 .get_refinement_edge_indices(level)
1423 .get_index_vector();
1424 edge_constrained_indices[j].clear();
1425 edge_constrained_indices[j].reserve(interface_indices.size());
1426 edge_constrained_values[j].resize(interface_indices.size());
1427 const IndexSet &locally_owned =
1428 data->get_dof_handler(selected_rows[j]).locally_owned_mg_dofs(level);
1429 for (const auto interface_index : interface_indices)
1430 if (locally_owned.is_element(interface_index))
1431 edge_constrained_indices[j].push_back(
1432 locally_owned.index_within_set(interface_index));
1433 have_interface_matrices |=
1434 Utilities::MPI::max(
1435 static_cast<unsigned int>(edge_constrained_indices[j].size()),
1436 data->get_vector_partitioner()->get_mpi_communicator()) > 0;
1437 }
1438 }
1439
1440
1441
1442 template <int dim, typename VectorType, typename VectorizedArrayType>
1443 void
1444 Base<dim, VectorType, VectorizedArrayType>::set_constrained_entries_to_one(
1445 VectorType &dst) const
1446 {
1447 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1448 {
1449 const std::vector<unsigned int> &constrained_dofs =
1450 data->get_constrained_dofs(selected_rows[j]);
1451 for (const auto constrained_dof : constrained_dofs)
1452 BlockHelper::subblock(dst, j).local_element(constrained_dof) = 1.;
1453 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1454 BlockHelper::subblock(dst, j).local_element(
1455 edge_constrained_indices[j][i]) = 1.;
1456 }
1457 }
1458
1459
1460
1461 template <int dim, typename VectorType, typename VectorizedArrayType>
1462 void
1463 Base<dim, VectorType, VectorizedArrayType>::vmult(VectorType & dst,
1464 const VectorType &src) const
1465 {
1466 using Number =
1467 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1468 dst = Number(0.);
1469 vmult_add(dst, src);
1470 }
1471
1472
1473
1474 template <int dim, typename VectorType, typename VectorizedArrayType>
1475 void
1476 Base<dim, VectorType, VectorizedArrayType>::vmult_add(
1477 VectorType & dst,
1478 const VectorType &src) const
1479 {
1480 mult_add(dst, src, false);
1481 }
1482
1483
1484
1485 template <int dim, typename VectorType, typename VectorizedArrayType>
1486 void
1487 Base<dim, VectorType, VectorizedArrayType>::Tvmult_add(
1488 VectorType & dst,
1489 const VectorType &src) const
1490 {
1491 mult_add(dst, src, true);
1492 }
1493
1494
1495
1496 template <int dim, typename VectorType, typename VectorizedArrayType>
1497 void
1498 Base<dim, VectorType, VectorizedArrayType>::adjust_ghost_range_if_necessary(
1499 const VectorType &src,
1500 const bool is_row) const
1501 {
1502 using Number =
1503 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1504 for (unsigned int i = 0; i < BlockHelper::n_blocks(src); ++i)
1505 {
1506 const unsigned int mf_component =
1507 is_row ? selected_rows[i] : selected_columns[i];
1508 // If both vectors use the same partitioner -> done
1509 if (BlockHelper::subblock(src, i).get_partitioner().get() ==
1510 data->get_dof_info(mf_component).vector_partitioner.get())
1511 continue;
1512
1513 // If not, assert that the local ranges are the same and reset to the
1514 // current partitioner
1515 Assert(BlockHelper::subblock(src, i)
1516 .get_partitioner()
1517 ->locally_owned_size() ==
1518 data->get_dof_info(mf_component)
1519 .vector_partitioner->locally_owned_size(),
1520 ExcMessage(
1521 "The vector passed to the vmult() function does not have "
1522 "the correct size for compatibility with MatrixFree."));
1523
1524 // copy the vector content to a temporary vector so that it does not get
1525 // lost
1526 LinearAlgebra::distributed::Vector<Number> copy_vec(
1527 BlockHelper::subblock(src, i));
1528 this->data->initialize_dof_vector(
1529 BlockHelper::subblock(const_cast<VectorType &>(src), i),
1530 mf_component);
1531 BlockHelper::subblock(const_cast<VectorType &>(src), i)
1532 .copy_locally_owned_data_from(copy_vec);
1533 }
1534 }
1535
1536
1537
1538 template <int dim, typename VectorType, typename VectorizedArrayType>
1539 void
1540 Base<dim, VectorType, VectorizedArrayType>::preprocess_constraints(
1541 VectorType & dst,
1542 const VectorType &src) const
1543 {
1544 using Number =
1545 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1546 adjust_ghost_range_if_necessary(src, false);
1547 adjust_ghost_range_if_necessary(dst, true);
1548
1549 // set zero Dirichlet values on the input vector (and remember the src and
1550 // dst values because we need to reset them at the end)
1551 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1552 {
1553 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1554 {
1555 edge_constrained_values[j][i] = std::pair<Number, Number>(
1556 BlockHelper::subblock(src, j).local_element(
1557 edge_constrained_indices[j][i]),
1558 BlockHelper::subblock(dst, j).local_element(
1559 edge_constrained_indices[j][i]));
1560 BlockHelper::subblock(const_cast<VectorType &>(src), j)
1561 .local_element(edge_constrained_indices[j][i]) = 0.;
1562 }
1563 }
1564 }
1565
1566
1567
1568 template <int dim, typename VectorType, typename VectorizedArrayType>
1569 void
1570 Base<dim, VectorType, VectorizedArrayType>::mult_add(
1571 VectorType & dst,
1572 const VectorType &src,
1573 const bool transpose) const
1574 {
1575 AssertDimension(dst.size(), src.size());
1576 AssertDimension(BlockHelper::n_blocks(dst), BlockHelper::n_blocks(src));
1577 AssertDimension(BlockHelper::n_blocks(dst), selected_rows.size());
1578 preprocess_constraints(dst, src);
1579 if (transpose)
1580 Tapply_add(dst, src);
1581 else
1582 apply_add(dst, src);
1583 postprocess_constraints(dst, src);
1584 }
1585
1586
1587
1588 template <int dim, typename VectorType, typename VectorizedArrayType>
1589 void
1590 Base<dim, VectorType, VectorizedArrayType>::postprocess_constraints(
1591 VectorType & dst,
1592 const VectorType &src) const
1593 {
1594 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1595 {
1596 const std::vector<unsigned int> &constrained_dofs =
1597 data->get_constrained_dofs(selected_rows[j]);
1598 for (const auto constrained_dof : constrained_dofs)
1599 BlockHelper::subblock(dst, j).local_element(constrained_dof) +=
1600 BlockHelper::subblock(src, j).local_element(constrained_dof);
1601 }
1602
1603 // reset edge constrained values, multiply by unit matrix and add into
1604 // destination
1605 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1606 {
1607 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1608 {
1609 BlockHelper::subblock(const_cast<VectorType &>(src), j)
1610 .local_element(edge_constrained_indices[j][i]) =
1611 edge_constrained_values[j][i].first;
1612 BlockHelper::subblock(dst, j).local_element(
1613 edge_constrained_indices[j][i]) =
1614 edge_constrained_values[j][i].second +
1615 edge_constrained_values[j][i].first;
1616 }
1617 }
1618 }
1619
1620
1621
1622 template <int dim, typename VectorType, typename VectorizedArrayType>
1623 void
1624 Base<dim, VectorType, VectorizedArrayType>::vmult_interface_down(
1625 VectorType & dst,
1626 const VectorType &src) const
1627 {
1628 using Number =
1629 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1630 AssertDimension(dst.size(), src.size());
1631 adjust_ghost_range_if_necessary(src, false);
1632 adjust_ghost_range_if_necessary(dst, true);
1633
1634 dst = Number(0.);
1635
1636 if (!have_interface_matrices)
1637 return;
1638
1639 // set zero Dirichlet values on the input vector (and remember the src and
1640 // dst values because we need to reset them at the end)
1641 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1642 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1643 {
1644 edge_constrained_values[j][i] = std::pair<Number, Number>(
1645 BlockHelper::subblock(src, j).local_element(
1646 edge_constrained_indices[j][i]),
1647 BlockHelper::subblock(dst, j).local_element(
1648 edge_constrained_indices[j][i]));
1649 BlockHelper::subblock(const_cast<VectorType &>(src), j)
1650 .local_element(edge_constrained_indices[j][i]) = 0.;
1651 }
1652
1653 apply_add(dst, src);
1654
1655 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1656 {
1657 unsigned int c = 0;
1658 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1659 {
1660 for (; c < edge_constrained_indices[j][i]; ++c)
1661 BlockHelper::subblock(dst, j).local_element(c) = 0.;
1662 ++c;
1663
1664 // reset the src values
1665 BlockHelper::subblock(const_cast<VectorType &>(src), j)
1666 .local_element(edge_constrained_indices[j][i]) =
1667 edge_constrained_values[j][i].first;
1668 }
1669 for (; c < BlockHelper::subblock(dst, j).locally_owned_size(); ++c)
1670 BlockHelper::subblock(dst, j).local_element(c) = 0.;
1671 }
1672 }
1673
1674
1675
1676 template <int dim, typename VectorType, typename VectorizedArrayType>
1677 void
1678 Base<dim, VectorType, VectorizedArrayType>::vmult_interface_up(
1679 VectorType & dst,
1680 const VectorType &src) const
1681 {
1682 using Number =
1683 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1684 AssertDimension(dst.size(), src.size());
1685 adjust_ghost_range_if_necessary(src, false);
1686 adjust_ghost_range_if_necessary(dst, true);
1687
1688 dst = Number(0.);
1689
1690 if (!have_interface_matrices)
1691 return;
1692
1693 VectorType src_cpy(src);
1694 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1695 {
1696 unsigned int c = 0;
1697 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1698 {
1699 for (; c < edge_constrained_indices[j][i]; ++c)
1700 BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1701 ++c;
1702 }
1703 for (; c < BlockHelper::subblock(src_cpy, j).locally_owned_size(); ++c)
1704 BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1705 }
1706
1707 apply_add(dst, src_cpy);
1708
1709 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1710 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1711 BlockHelper::subblock(dst, j).local_element(
1712 edge_constrained_indices[j][i]) = 0.;
1713 }
1714
1715
1716
1717 template <int dim, typename VectorType, typename VectorizedArrayType>
1718 void
1719 Base<dim, VectorType, VectorizedArrayType>::Tvmult(
1720 VectorType & dst,
1721 const VectorType &src) const
1722 {
1723 using Number =
1724 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1725 dst = Number(0.);
1726 Tvmult_add(dst, src);
1727 }
1728
1729
1730
1731 template <int dim, typename VectorType, typename VectorizedArrayType>
1732 std::size_t
1733 Base<dim, VectorType, VectorizedArrayType>::memory_consumption() const
1734 {
1735 return inverse_diagonal_entries.get() != nullptr ?
1736 inverse_diagonal_entries->memory_consumption() :
1737 sizeof(*this);
1738 }
1739
1740
1741
1742 template <int dim, typename VectorType, typename VectorizedArrayType>
1743 std::shared_ptr<const MatrixFree<
1744 dim,
1745 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
1746 VectorizedArrayType>>
1747 Base<dim, VectorType, VectorizedArrayType>::get_matrix_free() const
1748 {
1749 return data;
1750 }
1751
1752
1753
1754 template <int dim, typename VectorType, typename VectorizedArrayType>
1755 const std::shared_ptr<DiagonalMatrix<VectorType>> &
1756 Base<dim, VectorType, VectorizedArrayType>::get_matrix_diagonal_inverse()
1757 const
1758 {
1759 Assert(inverse_diagonal_entries.get() != nullptr &&
1760 inverse_diagonal_entries->m() > 0,
1761 ExcNotInitialized());
1762 return inverse_diagonal_entries;
1763 }
1764
1765
1766
1767 template <int dim, typename VectorType, typename VectorizedArrayType>
1768 const std::shared_ptr<DiagonalMatrix<VectorType>> &
1769 Base<dim, VectorType, VectorizedArrayType>::get_matrix_diagonal() const
1770 {
1771 Assert(diagonal_entries.get() != nullptr && diagonal_entries->m() > 0,
1772 ExcNotInitialized());
1773 return diagonal_entries;
1774 }
1775
1776
1777
1778 template <int dim, typename VectorType, typename VectorizedArrayType>
1779 void
1780 Base<dim, VectorType, VectorizedArrayType>::Tapply_add(
1781 VectorType & dst,
1782 const VectorType &src) const
1783 {
1784 apply_add(dst, src);
1785 }
1786
1787
1788
1789 template <int dim, typename VectorType, typename VectorizedArrayType>
1790 void
1791 Base<dim, VectorType, VectorizedArrayType>::precondition_Jacobi(
1792 VectorType & dst,
1793 const VectorType & src,
1794 const typename Base<dim, VectorType, VectorizedArrayType>::value_type omega)
1795 const
1796 {
1797 Assert(inverse_diagonal_entries.get() && inverse_diagonal_entries->m() > 0,
1798 ExcNotInitialized());
1799 inverse_diagonal_entries->vmult(dst, src);
1800 dst *= omega;
1801 }
1802
1803
1804
1805 //------------------------- MGInterfaceOperator ------------------------------
1806
1807 template <typename OperatorType>
1808 MGInterfaceOperator<OperatorType>::MGInterfaceOperator()
1809 : Subscriptor()
1810 , mf_base_operator(nullptr)
1811 {}
1812
1813
1814
1815 template <typename OperatorType>
1816 void
1817 MGInterfaceOperator<OperatorType>::clear()
1818 {
1819 mf_base_operator = nullptr;
1820 }
1821
1822
1823
1824 template <typename OperatorType>
1825 void
1826 MGInterfaceOperator<OperatorType>::initialize(const OperatorType &operator_in)
1827 {
1828 mf_base_operator = &operator_in;
1829 }
1830
1831
1832
1833 template <typename OperatorType>
1834 template <typename VectorType>
1835 void
1836 MGInterfaceOperator<OperatorType>::vmult(VectorType & dst,
1837 const VectorType &src) const
1838 {
1839#ifndef DEAL_II_MSVC
1840 static_assert(
1841 std::is_same<typename VectorType::value_type, value_type>::value,
1842 "The vector type must be based on the same value type as this "
1843 "operator");
1844#endif
1845
1846 Assert(mf_base_operator != nullptr, ExcNotInitialized());
1847
1848 mf_base_operator->vmult_interface_down(dst, src);
1849 }
1850
1851
1852
1853 template <typename OperatorType>
1854 template <typename VectorType>
1855 void
1856 MGInterfaceOperator<OperatorType>::Tvmult(VectorType & dst,
1857 const VectorType &src) const
1858 {
1859#ifndef DEAL_II_MSVC
1860 static_assert(
1861 std::is_same<typename VectorType::value_type, value_type>::value,
1862 "The vector type must be based on the same value type as this "
1863 "operator");
1864#endif
1865
1866 Assert(mf_base_operator != nullptr, ExcNotInitialized());
1867
1868 mf_base_operator->vmult_interface_up(dst, src);
1869 }
1870
1871
1872
1873 template <typename OperatorType>
1874 template <typename VectorType>
1875 void
1876 MGInterfaceOperator<OperatorType>::initialize_dof_vector(
1877 VectorType &vec) const
1878 {
1879 Assert(mf_base_operator != nullptr, ExcNotInitialized());
1880
1881 mf_base_operator->initialize_dof_vector(vec);
1882 }
1883
1884
1885
1886 //-----------------------------MassOperator----------------------------------
1887
1888 template <int dim,
1889 int fe_degree,
1890 int n_q_points_1d,
1891 int n_components,
1892 typename VectorType,
1893 typename VectorizedArrayType>
1894 MassOperator<dim,
1895 fe_degree,
1896 n_q_points_1d,
1897 n_components,
1898 VectorType,
1899 VectorizedArrayType>::MassOperator()
1900 : Base<dim, VectorType, VectorizedArrayType>()
1901 {
1902 AssertThrow(
1903 IsBlockVector<VectorType>::value == false,
1904 ExcNotImplemented(
1905 "This class only supports the non-blocked vector variant of the Base "
1906 "operator because only a single FEEvaluation object is used in the "
1907 "apply function."));
1908 }
1909
1910
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 void
1919 MassOperator<dim,
1920 fe_degree,
1921 n_q_points_1d,
1922 n_components,
1923 VectorType,
1924 VectorizedArrayType>::compute_diagonal()
1925 {
1926 Assert((Base<dim, VectorType, VectorizedArrayType>::data.get() != nullptr),
1927 ExcNotInitialized());
1928 Assert(this->selected_rows == this->selected_columns,
1929 ExcMessage("This function is only implemented for square (not "
1930 "rectangular) operators."));
1931
1932 this->inverse_diagonal_entries =
1933 std::make_shared<DiagonalMatrix<VectorType>>();
1934 this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1935 VectorType &inverse_diagonal_vector =
1936 this->inverse_diagonal_entries->get_vector();
1937 VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1938 this->initialize_dof_vector(inverse_diagonal_vector);
1939 this->initialize_dof_vector(diagonal_vector);
1940
1941 // Set up the action of the mass matrix in a way that's compatible with
1942 // MatrixFreeTools::compute_diagonal:
1943 auto diagonal_evaluation = [](auto &integrator) {
1944 integrator.evaluate(EvaluationFlags::values);
1945 for (unsigned int q = 0; q < integrator.n_q_points; ++q)
1946 integrator.submit_value(integrator.get_value(q), q);
1947 integrator.integrate(EvaluationFlags::values);
1948 };
1949
1950 std::function<void(
1951 FEEvaluation<
1952 dim,
1953 fe_degree,
1954 n_q_points_1d,
1955 n_components,
1956 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
1957 VectorizedArrayType> &)>
1958 diagonal_evaluation_f(diagonal_evaluation);
1959
1960 Assert(this->selected_rows.size() > 0, ExcInternalError());
1961 for (unsigned int block_n = 0; block_n < this->selected_rows.size();
1962 ++block_n)
1963 MatrixFreeTools::compute_diagonal(*this->data,
1964 BlockHelper::subblock(diagonal_vector,
1965 block_n),
1966 diagonal_evaluation_f,
1967 this->selected_rows[block_n]);
1968
1969 // Constrained entries will create zeros on the main diagonal, which we
1970 // don't want
1971 this->set_constrained_entries_to_one(diagonal_vector);
1972
1973 inverse_diagonal_vector = diagonal_vector;
1974
1975 for (unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
1976 ++i)
1977 {
1978#ifdef DEBUG
1979 // only define the type alias in debug mode to avoid a warning
1980 using Number =
1981 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1982 Assert(diagonal_vector.local_element(i) > Number(0),
1983 ExcInternalError());
1984#endif
1985 inverse_diagonal_vector.local_element(i) =
1986 1. / inverse_diagonal_vector.local_element(i);
1987 }
1988
1989 // We never need ghost values so don't update them
1990 }
1991
1992
1993
1994 template <int dim,
1995 int fe_degree,
1996 int n_q_points_1d,
1997 int n_components,
1998 typename VectorType,
1999 typename VectorizedArrayType>
2000 void
2001 MassOperator<dim,
2002 fe_degree,
2003 n_q_points_1d,
2004 n_components,
2005 VectorType,
2006 VectorizedArrayType>::compute_lumped_diagonal()
2007 {
2008 using Number =
2009 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2010 Assert((Base<dim, VectorType, VectorizedArrayType>::data.get() != nullptr),
2011 ExcNotInitialized());
2012 Assert(this->selected_rows == this->selected_columns,
2013 ExcMessage("This function is only implemented for square (not "
2014 "rectangular) operators."));
2015
2016 inverse_lumped_diagonal_entries =
2017 std::make_shared<DiagonalMatrix<VectorType>>();
2018 lumped_diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
2019 VectorType &inverse_lumped_diagonal_vector =
2020 inverse_lumped_diagonal_entries->get_vector();
2021 VectorType &lumped_diagonal_vector = lumped_diagonal_entries->get_vector();
2022 this->initialize_dof_vector(inverse_lumped_diagonal_vector);
2023 this->initialize_dof_vector(lumped_diagonal_vector);
2024
2025 // Re-use the inverse_lumped_diagonal_vector as the vector of 1s
2026 inverse_lumped_diagonal_vector = Number(1.);
2027 apply_add(lumped_diagonal_vector, inverse_lumped_diagonal_vector);
2028 this->set_constrained_entries_to_one(lumped_diagonal_vector);
2029
2030 const size_type locally_owned_size =
2031 inverse_lumped_diagonal_vector.locally_owned_size();
2032 // A caller may request a lumped diagonal matrix when it doesn't make sense
2033 // (e.g., an element with negative-mean basis functions). Avoid division by
2034 // zero so we don't cause a floating point exception but permit negative
2035 // entries here.
2036 for (size_type i = 0; i < locally_owned_size; ++i)
2037 {
2038 if (lumped_diagonal_vector.local_element(i) == Number(0.))
2039 inverse_lumped_diagonal_vector.local_element(i) = Number(1.);
2040 else
2041 inverse_lumped_diagonal_vector.local_element(i) =
2042 Number(1.) / lumped_diagonal_vector.local_element(i);
2043 }
2044
2045 inverse_lumped_diagonal_vector.update_ghost_values();
2046 lumped_diagonal_vector.update_ghost_values();
2047 }
2048
2049
2050
2051 template <int dim,
2052 int fe_degree,
2053 int n_q_points_1d,
2054 int n_components,
2055 typename VectorType,
2056 typename VectorizedArrayType>
2057 const std::shared_ptr<DiagonalMatrix<VectorType>> &
2058 MassOperator<dim,
2059 fe_degree,
2060 n_q_points_1d,
2061 n_components,
2062 VectorType,
2063 VectorizedArrayType>::get_matrix_lumped_diagonal_inverse() const
2064 {
2065 Assert(inverse_lumped_diagonal_entries.get() != nullptr &&
2066 inverse_lumped_diagonal_entries->m() > 0,
2067 ExcNotInitialized());
2068 return inverse_lumped_diagonal_entries;
2069 }
2070
2071
2072
2073 template <int dim,
2074 int fe_degree,
2075 int n_q_points_1d,
2076 int n_components,
2077 typename VectorType,
2078 typename VectorizedArrayType>
2079 const std::shared_ptr<DiagonalMatrix<VectorType>> &
2080 MassOperator<dim,
2081 fe_degree,
2082 n_q_points_1d,
2083 n_components,
2084 VectorType,
2085 VectorizedArrayType>::get_matrix_lumped_diagonal() const
2086 {
2087 Assert(lumped_diagonal_entries.get() != nullptr &&
2088 lumped_diagonal_entries->m() > 0,
2089 ExcNotInitialized());
2090 return lumped_diagonal_entries;
2091 }
2092
2093
2094
2095 template <int dim,
2096 int fe_degree,
2097 int n_q_points_1d,
2098 int n_components,
2099 typename VectorType,
2100 typename VectorizedArrayType>
2101 void
2102 MassOperator<dim,
2103 fe_degree,
2104 n_q_points_1d,
2105 n_components,
2106 VectorType,
2107 VectorizedArrayType>::apply_add(VectorType & dst,
2108 const VectorType &src) const
2109 {
2110 Base<dim, VectorType, VectorizedArrayType>::data->cell_loop(
2111 &MassOperator::local_apply_cell, this, dst, src);
2112 }
2113
2114
2115
2116 template <int dim,
2117 int fe_degree,
2118 int n_q_points_1d,
2119 int n_components,
2120 typename VectorType,
2121 typename VectorizedArrayType>
2122 void
2123 MassOperator<dim,
2124 fe_degree,
2125 n_q_points_1d,
2126 n_components,
2127 VectorType,
2128 VectorizedArrayType>::
2129 local_apply_cell(
2130 const MatrixFree<
2131 dim,
2132 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2133 VectorizedArrayType> & data,
2134 VectorType & dst,
2135 const VectorType & src,
2136 const std::pair<unsigned int, unsigned int> &cell_range) const
2137 {
2138 using Number =
2139 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2140 FEEvaluation<dim,
2141 fe_degree,
2142 n_q_points_1d,
2143 n_components,
2144 Number,
2145 VectorizedArrayType>
2146 phi(data, this->selected_rows[0]);
2147 for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2148 {
2149 phi.reinit(cell);
2150 phi.read_dof_values(src);
2151 phi.evaluate(EvaluationFlags::values);
2152 for (unsigned int q = 0; q < phi.n_q_points; ++q)
2153 phi.submit_value(phi.get_value(q), q);
2154 phi.integrate(EvaluationFlags::values);
2155 phi.distribute_local_to_global(dst);
2156 }
2157 }
2158
2159
2160 //-----------------------------LaplaceOperator----------------------------------
2161
2162 template <int dim,
2163 int fe_degree,
2164 int n_q_points_1d,
2165 int n_components,
2166 typename VectorType,
2167 typename VectorizedArrayType>
2168 LaplaceOperator<dim,
2169 fe_degree,
2170 n_q_points_1d,
2171 n_components,
2172 VectorType,
2173 VectorizedArrayType>::LaplaceOperator()
2174 : Base<dim, VectorType, VectorizedArrayType>()
2175 {}
2176
2177
2178
2179 template <int dim,
2180 int fe_degree,
2181 int n_q_points_1d,
2182 int n_components,
2183 typename VectorType,
2184 typename VectorizedArrayType>
2185 void
2186 LaplaceOperator<dim,
2187 fe_degree,
2188 n_q_points_1d,
2189 n_components,
2190 VectorType,
2191 VectorizedArrayType>::clear()
2192 {
2193 Base<dim, VectorType, VectorizedArrayType>::clear();
2194 scalar_coefficient.reset();
2195 }
2196
2197
2198
2199 template <int dim,
2200 int fe_degree,
2201 int n_q_points_1d,
2202 int n_components,
2203 typename VectorType,
2204 typename VectorizedArrayType>
2205 void
2206 LaplaceOperator<dim,
2207 fe_degree,
2208 n_q_points_1d,
2209 n_components,
2210 VectorType,
2211 VectorizedArrayType>::
2212 set_coefficient(
2213 const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient_)
2214 {
2215 scalar_coefficient = scalar_coefficient_;
2216 }
2217
2218
2219
2220 template <int dim,
2221 int fe_degree,
2222 int n_q_points_1d,
2223 int n_components,
2224 typename VectorType,
2225 typename VectorizedArrayType>
2226 std::shared_ptr<Table<2, VectorizedArrayType>>
2227 LaplaceOperator<dim,
2228 fe_degree,
2229 n_q_points_1d,
2230 n_components,
2231 VectorType,
2232 VectorizedArrayType>::get_coefficient()
2233 {
2234 Assert(scalar_coefficient.get(), ExcNotInitialized());
2235 return scalar_coefficient;
2236 }
2237
2238
2239
2240 template <int dim,
2241 int fe_degree,
2242 int n_q_points_1d,
2243 int n_components,
2244 typename VectorType,
2245 typename VectorizedArrayType>
2246 void
2247 LaplaceOperator<dim,
2248 fe_degree,
2249 n_q_points_1d,
2250 n_components,
2251 VectorType,
2252 VectorizedArrayType>::compute_diagonal()
2253 {
2254 using Number =
2255 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2256 Assert((Base<dim, VectorType, VectorizedArrayType>::data.get() != nullptr),
2257 ExcNotInitialized());
2258
2259 this->inverse_diagonal_entries =
2260 std::make_shared<DiagonalMatrix<VectorType>>();
2261 this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
2262 VectorType &inverse_diagonal_vector =
2263 this->inverse_diagonal_entries->get_vector();
2264 VectorType &diagonal_vector = this->diagonal_entries->get_vector();
2265 this->initialize_dof_vector(inverse_diagonal_vector);
2266 this->initialize_dof_vector(diagonal_vector);
2267
2268 this->data->cell_loop(&LaplaceOperator::local_diagonal_cell,
2269 this,
2270 diagonal_vector,
2271 /*unused*/ diagonal_vector);
2272 this->set_constrained_entries_to_one(diagonal_vector);
2273
2274 inverse_diagonal_vector = diagonal_vector;
2275
2276 for (unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
2277 ++i)
2278 if (std::abs(inverse_diagonal_vector.local_element(i)) >
2279 std::sqrt(std::numeric_limits<Number>::epsilon()))
2280 inverse_diagonal_vector.local_element(i) =
2281 1. / inverse_diagonal_vector.local_element(i);
2282 else
2283 inverse_diagonal_vector.local_element(i) = 1.;
2284
2285 // We never need ghost values so don't update them
2286 }
2287
2288
2289
2290 template <int dim,
2291 int fe_degree,
2292 int n_q_points_1d,
2293 int n_components,
2294 typename VectorType,
2295 typename VectorizedArrayType>
2296 void
2297 LaplaceOperator<dim,
2298 fe_degree,
2299 n_q_points_1d,
2300 n_components,
2301 VectorType,
2302 VectorizedArrayType>::apply_add(VectorType & dst,
2303 const VectorType &src) const
2304 {
2305 Base<dim, VectorType, VectorizedArrayType>::data->cell_loop(
2306 &LaplaceOperator::local_apply_cell, this, dst, src);
2307 }
2308
2309 namespace Implementation
2310 {
2311 template <typename VectorizedArrayType>
2312 bool
2313 non_negative(const VectorizedArrayType &n)
2314 {
2315 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2316 if (n[v] < 0.)
2317 return false;
2318
2319 return true;
2320 }
2321 } // namespace Implementation
2322
2323
2324
2325 template <int dim,
2326 int fe_degree,
2327 int n_q_points_1d,
2328 int n_components,
2329 typename VectorType,
2330 typename VectorizedArrayType>
2331 template <int n_components_compute>
2332 void
2333 LaplaceOperator<dim,
2334 fe_degree,
2335 n_q_points_1d,
2336 n_components,
2337 VectorType,
2338 VectorizedArrayType>::
2339 do_operation_on_cell(
2340 FEEvaluation<
2341 dim,
2342 fe_degree,
2343 n_q_points_1d,
2344 n_components_compute,
2345 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2346 VectorizedArrayType> &phi,
2347 const unsigned int cell) const
2348 {
2349 phi.evaluate(EvaluationFlags::gradients);
2350 if (scalar_coefficient.get())
2351 {
2352 Assert(scalar_coefficient->size(1) == 1 ||
2353 scalar_coefficient->size(1) == phi.n_q_points,
2354 ExcMessage("The number of columns in the coefficient table must "
2355 "be either 1 or the number of quadrature points " +
2356 std::to_string(phi.n_q_points) +
2357 ", but the given value was " +
2358 std::to_string(scalar_coefficient->size(1))));
2359 if (scalar_coefficient->size(1) == phi.n_q_points)
2360 for (unsigned int q = 0; q < phi.n_q_points; ++q)
2361 {
2362 Assert(Implementation::non_negative(
2363 (*scalar_coefficient)(cell, q)),
2364 ExcMessage("Coefficient must be non-negative"));
2365 phi.submit_gradient((*scalar_coefficient)(cell, q) *
2366 phi.get_gradient(q),
2367 q);
2368 }
2369 else
2370 {
2371 Assert(Implementation::non_negative((*scalar_coefficient)(cell, 0)),
2372 ExcMessage("Coefficient must be non-negative"));
2373 const VectorizedArrayType coefficient =
2374 (*scalar_coefficient)(cell, 0);
2375 for (unsigned int q = 0; q < phi.n_q_points; ++q)
2376 phi.submit_gradient(coefficient * phi.get_gradient(q), q);
2377 }
2378 }
2379 else
2380 {
2381 for (unsigned int q = 0; q < phi.n_q_points; ++q)
2382 {
2383 phi.submit_gradient(phi.get_gradient(q), q);
2384 }
2385 }
2386 phi.integrate(EvaluationFlags::gradients);
2387 }
2388
2389
2390
2391 template <int dim,
2392 int fe_degree,
2393 int n_q_points_1d,
2394 int n_components,
2395 typename VectorType,
2396 typename VectorizedArrayType>
2397 void
2398 LaplaceOperator<dim,
2399 fe_degree,
2400 n_q_points_1d,
2401 n_components,
2402 VectorType,
2403 VectorizedArrayType>::
2404 local_apply_cell(
2405 const MatrixFree<
2406 dim,
2407 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2408 VectorizedArrayType> & data,
2409 VectorType & dst,
2410 const VectorType & src,
2411 const std::pair<unsigned int, unsigned int> &cell_range) const
2412 {
2413 using Number =
2414 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2415 FEEvaluation<dim,
2416 fe_degree,
2417 n_q_points_1d,
2418 n_components,
2419 Number,
2420 VectorizedArrayType>
2421 phi(data, this->selected_rows[0]);
2422 for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2423 {
2424 phi.reinit(cell);
2425 phi.read_dof_values(src);
2426 do_operation_on_cell(phi, cell);
2427 phi.distribute_local_to_global(dst);
2428 }
2429 }
2430
2431
2432 template <int dim,
2433 int fe_degree,
2434 int n_q_points_1d,
2435 int n_components,
2436 typename VectorType,
2437 typename VectorizedArrayType>
2438 void
2439 LaplaceOperator<dim,
2440 fe_degree,
2441 n_q_points_1d,
2442 n_components,
2443 VectorType,
2444 VectorizedArrayType>::
2445 local_diagonal_cell(
2446 const MatrixFree<
2447 dim,
2448 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2449 VectorizedArrayType> &data,
2450 VectorType & dst,
2451 const VectorType &,
2452 const std::pair<unsigned int, unsigned int> &cell_range) const
2453 {
2454 using Number =
2455 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2456
2457 FEEvaluation<dim, fe_degree, n_q_points_1d, 1, Number, VectorizedArrayType>
2458 eval(data, this->selected_rows[0]);
2459 FEEvaluation<dim,
2460 fe_degree,
2461 n_q_points_1d,
2462 n_components,
2463 Number,
2464 VectorizedArrayType>
2465 eval_vector(data, this->selected_rows[0]);
2466 for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2467 {
2468 eval.reinit(cell);
2469 eval_vector.reinit(cell);
2470 // This function assumes that we have the same result on all
2471 // components, so we only need to go through the columns of one scalar
2472 // component, for which we have created a separate evaluator (attached
2473 // to the first component, but the component does not matter because
2474 // we only use the underlying integrals)
2475 for (unsigned int i = 0; i < eval.dofs_per_cell; ++i)
2476 {
2477 for (unsigned int j = 0; j < eval.dofs_per_cell; ++j)
2478 eval.begin_dof_values()[j] = VectorizedArrayType();
2479 eval.begin_dof_values()[i] = 1.;
2480
2481 do_operation_on_cell(eval, cell);
2482
2483 // We now pick up the value on the diagonal (row i) and broadcast
2484 // it to a second evaluator for all vector components, which we
2485 // will distribute to the result vector afterwards
2486 for (unsigned int c = 0; c < n_components; ++c)
2487 eval_vector
2488 .begin_dof_values()[i + c * eval_vector.dofs_per_component] =
2489 eval.begin_dof_values()[i];
2490 }
2491 eval_vector.distribute_local_to_global(dst);
2492 }
2493 }
2494
2495
2496} // end of namespace MatrixFreeOperators
2497
2498
2499DEAL_II_NAMESPACE_CLOSE
2500
2501#endif
make_array_view
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:704
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:2358
IndexSet::index_within_set
size_type index_within_set(const size_type global_index) const
Definition index_set.h:1883
IndexSet::is_element
bool is_element(const size_type index) const
Definition index_set.h:1776
MatrixFreeOperators::Base::~Base
virtual ~Base() override=default
MatrixFreeOperators::Base::selected_rows
std::vector< unsigned int > selected_rows
Definition operators.h:454
MatrixFreeOperators::Base::Tvmult_add
void Tvmult_add(VectorType &dst, const VectorType &src) const
Definition operators.h:1487
MatrixFreeOperators::Base::set_constrained_entries_to_one
void set_constrained_entries_to_one(VectorType &dst) const
Definition operators.h:1444
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:1476
MatrixFreeOperators::Base::vmult_interface_down
void vmult_interface_down(VectorType &dst, const VectorType &src) const
Definition operators.h:1624
MatrixFreeOperators::Base::n
size_type n() const
Definition operators.h:1241
MatrixFreeOperators::Base::preprocess_constraints
void preprocess_constraints(VectorType &dst, const VectorType &src) const
Definition operators.h:1540
MatrixFreeOperators::Base::mult_add
void mult_add(VectorType &dst, const VectorType &src, const bool transpose) const
Definition operators.h:1570
MatrixFreeOperators::Base::edge_constrained_indices
std::vector< std::vector< unsigned int > > edge_constrained_indices
Definition operators.h:466
MatrixFreeOperators::Base::get_matrix_diagonal
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal() const
Definition operators.h:1769
MatrixFreeOperators::Base::Tvmult
void Tvmult(VectorType &dst, const VectorType &src) const
Definition operators.h:1719
MatrixFreeOperators::Base::edge_constrained_values
std::vector< std::vector< std::pair< value_type, value_type > > > edge_constrained_values
Definition operators.h:472
MatrixFreeOperators::Base::m
size_type m() const
Definition operators.h:1227
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:1305
MatrixFreeOperators::Base::selected_columns
std::vector< unsigned int > selected_columns
Definition operators.h:460
MatrixFreeOperators::Base::get_matrix_free
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > get_matrix_free() const
Definition operators.h:1747
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:1374
MatrixFreeOperators::Base::get_matrix_diagonal_inverse
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal_inverse() const
Definition operators.h:1756
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:1358
MatrixFreeOperators::Base::initialize_dof_vector
void initialize_dof_vector(VectorType &vec) const
Definition operators.h:1280
MatrixFreeOperators::Base::diagonal_entries
std::shared_ptr< DiagonalMatrix< VectorType > > diagonal_entries
Definition operators.h:442
MatrixFreeOperators::Base::Tapply_add
virtual void Tapply_add(VectorType &dst, const VectorType &src) const
Definition operators.h:1780
MatrixFreeOperators::Base::data
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > data
Definition operators.h:436
MatrixFreeOperators::Base::adjust_ghost_range_if_necessary
void adjust_ghost_range_if_necessary(const VectorType &vec, const bool is_row) const
Definition operators.h:1498
MatrixFreeOperators::Base::memory_consumption
virtual std::size_t memory_consumption() const
Definition operators.h:1733
MatrixFreeOperators::Base::inverse_diagonal_entries
std::shared_ptr< DiagonalMatrix< VectorType > > inverse_diagonal_entries
Definition operators.h:448
MatrixFreeOperators::Base::el
value_type el(const unsigned int row, const unsigned int col) const
Definition operators.h:1265
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:199
MatrixFreeOperators::Base::precondition_Jacobi
void precondition_Jacobi(VectorType &dst, const VectorType &src, const value_type omega) const
Definition operators.h:1791
MatrixFreeOperators::Base::value_type
typename VectorType::value_type value_type
Definition operators.h:194
MatrixFreeOperators::Base::vmult
void vmult(VectorType &dst, const VectorType &src) const
Definition operators.h:1463
MatrixFreeOperators::Base::vmult_interface_up
void vmult_interface_up(VectorType &dst, const VectorType &src) const
Definition operators.h:1678
MatrixFreeOperators::Base::postprocess_constraints
void postprocess_constraints(VectorType &dst, const VectorType &src) const
Definition operators.h:1590
MatrixFreeOperators::CellwiseInverseMassMatrix::fe_eval
const FEEvaluationBase< dim, n_components, Number, false, VectorizedArrayType > & fe_eval
Definition operators.h:743
MatrixFreeOperators::CellwiseInverseMassMatrix::fill_inverse_JxW_values
void fill_inverse_JxW_values(AlignedVector< VectorizedArrayType > &inverse_jxw) const
Definition operators.h:1053
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:1192
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:1113
MatrixFreeOperators::LaplaceOperator::get_coefficient
std::shared_ptr< Table< 2, VectorizedArrayType > > get_coefficient()
Definition operators.h:2232
MatrixFreeOperators::LaplaceOperator::value_type
typename Base< dim, VectorType, VectorizedArrayType >::value_type value_type
Definition operators.h:877
MatrixFreeOperators::LaplaceOperator::compute_diagonal
virtual void compute_diagonal() override
Definition operators.h:2252
MatrixFreeOperators::LaplaceOperator::scalar_coefficient
std::shared_ptr< Table< 2, VectorizedArrayType > > scalar_coefficient
Definition operators.h:1014
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:2445
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:2404
MatrixFreeOperators::LaplaceOperator::apply_add
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition operators.h:2302
MatrixFreeOperators::LaplaceOperator::size_type
typename Base< dim, VectorType, VectorizedArrayType >::size_type size_type
Definition operators.h:883
MatrixFreeOperators::LaplaceOperator::set_coefficient
void set_coefficient(const std::shared_ptr< Table< 2, VectorizedArrayType > > &scalar_coefficient)
Definition operators.h:2212
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:2191
MatrixFreeOperators::MGInterfaceOperator::value_type
typename OperatorType::value_type value_type
Definition operators.h:544
MatrixFreeOperators::MGInterfaceOperator::mf_base_operator
SmartPointer< const OperatorType > mf_base_operator
Definition operators.h:594
MatrixFreeOperators::MGInterfaceOperator::vmult
void vmult(VectorType &dst, const VectorType &src) const
Definition operators.h:1836
MatrixFreeOperators::MGInterfaceOperator::initialize
void initialize(const OperatorType &operator_in)
Definition operators.h:1826
MatrixFreeOperators::MGInterfaceOperator::Tvmult
void Tvmult(VectorType &dst, const VectorType &src) const
Definition operators.h:1856
MatrixFreeOperators::MGInterfaceOperator::initialize_dof_vector
void initialize_dof_vector(VectorType &vec) const
Definition operators.h:1876
MatrixFreeOperators::MGInterfaceOperator::size_type
typename OperatorType::size_type size_type
Definition operators.h:549
MatrixFreeOperators::MassOperator::lumped_diagonal_entries
std::shared_ptr< DiagonalMatrix< VectorType > > lumped_diagonal_entries
Definition operators.h:842
MatrixFreeOperators::MassOperator::size_type
typename Base< dim, VectorType, VectorizedArrayType >::size_type size_type
Definition operators.h:775
MatrixFreeOperators::MassOperator::apply_add
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition operators.h:2107
MatrixFreeOperators::MassOperator::get_matrix_lumped_diagonal
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_lumped_diagonal() const
Definition operators.h:2085
MatrixFreeOperators::MassOperator::get_matrix_lumped_diagonal_inverse
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_lumped_diagonal_inverse() const
Definition operators.h:2063
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:2129
MatrixFreeOperators::MassOperator::compute_diagonal
virtual void compute_diagonal() override
Definition operators.h:1924
MatrixFreeOperators::MassOperator::value_type
typename Base< dim, VectorType, VectorizedArrayType >::value_type value_type
Definition operators.h:769
MatrixFreeOperators::MassOperator::inverse_lumped_diagonal_entries
std::shared_ptr< DiagonalMatrix< VectorType > > inverse_lumped_diagonal_entries
Definition operators.h:848
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition derivative_form.h:586
diagonal_matrix.h
fe_evaluation.h
level
unsigned int level
Definition grid_out.cc:4618
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exceptions.h:1614
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:1787
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exceptions.h:1855
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcNotInitialized
static ::ExceptionBase & ExcNotInitialized()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition exceptions.h:1703
la_parallel_vector.h
matrix_free.h
mg_constrained_dofs.h
MatrixFreeOperators::BlockHelper::subblock
std::enable_if_t< IsBlockVector< VectorType >::value, typename VectorType::BlockType & > subblock(VectorType &vector, unsigned int block_no)
Definition operators.h:65
MatrixFreeOperators::BlockHelper::collect_sizes
std::enable_if_t< IsBlockVector< VectorType >::value, void > collect_sizes(VectorType &vector)
Definition operators.h:96
MatrixFreeOperators::BlockHelper::n_blocks
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition operators.h:50
MatrixFreeOperators::Implementation::non_negative
bool non_negative(const VectorizedArrayType &n)
Definition operators.h:2313
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 > &)> &local_vmult, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_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:447
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition types.h:213
std::sqrt
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:5946
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6008
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.7