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