Reference documentation for deal.II version 9.6.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// 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/exceptions.h>
23#include <deal.II/base/subscriptor.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 Subscriptor
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 Subscriptor
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 SmartPointer<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 : Subscriptor()
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 (void)col;
1268 Assert(row == col, ExcNotImplemented());
1269 Assert(inverse_diagonal_entries.get() != nullptr &&
1270 inverse_diagonal_entries->m() > 0,
1271 ExcNotInitialized());
1272 return 1.0 / (*inverse_diagonal_entries)(row, row);
1273 }
1274
1275
1276
1277 template <int dim, typename VectorType, typename VectorizedArrayType>
1278 void
1279 Base<dim, VectorType, VectorizedArrayType>::initialize_dof_vector(
1280 VectorType &vec) const
1281 {
1282 Assert(data.get() != nullptr, ExcNotInitialized());
1283 AssertDimension(BlockHelper::n_blocks(vec), selected_rows.size());
1284 for (unsigned int i = 0; i < BlockHelper::n_blocks(vec); ++i)
1285 {
1286 const unsigned int index = selected_rows[i];
1287 if (!BlockHelper::subblock(vec, index)
1288 .partitioners_are_compatible(
1289 *data->get_dof_info(index).vector_partitioner))
1290 data->initialize_dof_vector(BlockHelper::subblock(vec, index), index);
1291
1292 Assert(BlockHelper::subblock(vec, index)
1293 .partitioners_are_globally_compatible(
1294 *data->get_dof_info(index).vector_partitioner),
1295 ExcInternalError());
1296 }
1297 BlockHelper::collect_sizes(vec);
1298 }
1299
1300
1301
1302 template <int dim, typename VectorType, typename VectorizedArrayType>
1303 void
1304 Base<dim, VectorType, VectorizedArrayType>::initialize(
1305 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
1306 data_,
1307 const std::vector<unsigned int> &given_row_selection,
1308 const std::vector<unsigned int> &given_column_selection)
1309 {
1310 data = data_;
1311
1312 selected_rows.clear();
1313 selected_columns.clear();
1314 if (given_row_selection.empty())
1315 for (unsigned int i = 0; i < data_->n_components(); ++i)
1316 selected_rows.push_back(i);
1317 else
1318 {
1319 for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1320 {
1321 AssertIndexRange(given_row_selection[i], data_->n_components());
1322 for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1323 if (j != i)
1324 Assert(given_row_selection[j] != given_row_selection[i],
1325 ExcMessage("Given row indices must be unique"));
1326
1327 selected_rows.push_back(given_row_selection[i]);
1328 }
1329 }
1330 if (given_column_selection.empty())
1331 selected_columns = selected_rows;
1332 else
1333 {
1334 for (unsigned int i = 0; i < given_column_selection.size(); ++i)
1335 {
1336 AssertIndexRange(given_column_selection[i], data_->n_components());
1337 for (unsigned int j = 0; j < given_column_selection.size(); ++j)
1338 if (j != i)
1339 Assert(given_column_selection[j] != given_column_selection[i],
1340 ExcMessage("Given column indices must be unique"));
1341
1342 selected_columns.push_back(given_column_selection[i]);
1343 }
1344 }
1345
1346 edge_constrained_indices.clear();
1347 edge_constrained_indices.resize(selected_rows.size());
1348 edge_constrained_values.clear();
1349 edge_constrained_values.resize(selected_rows.size());
1350 have_interface_matrices = false;
1351 }
1352
1353
1354
1355 template <int dim, typename VectorType, typename VectorizedArrayType>
1356 void
1357 Base<dim, VectorType, VectorizedArrayType>::initialize(
1358 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
1359 data_,
1360 const MGConstrainedDoFs &mg_constrained_dofs,
1361 const unsigned int level,
1362 const std::vector<unsigned int> &given_row_selection)
1363 {
1364 std::vector<MGConstrainedDoFs> mg_constrained_dofs_vector(
1365 1, mg_constrained_dofs);
1366 initialize(data_, mg_constrained_dofs_vector, level, given_row_selection);
1367 }
1368
1369
1370
1371 template <int dim, typename VectorType, typename VectorizedArrayType>
1372 void
1373 Base<dim, VectorType, VectorizedArrayType>::initialize(
1374 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
1375 data_,
1376 const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
1377 const unsigned int level,
1378 const std::vector<unsigned int> &given_row_selection)
1379 {
1380 AssertThrow(level != numbers::invalid_unsigned_int,
1381 ExcMessage("level is not set"));
1382
1383 selected_rows.clear();
1384 selected_columns.clear();
1385 if (given_row_selection.empty())
1386 for (unsigned int i = 0; i < data_->n_components(); ++i)
1387 selected_rows.push_back(i);
1388 else
1389 {
1390 for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1391 {
1392 AssertIndexRange(given_row_selection[i], data_->n_components());
1393 for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1394 if (j != i)
1395 Assert(given_row_selection[j] != given_row_selection[i],
1396 ExcMessage("Given row indices must be unique"));
1397
1398 selected_rows.push_back(given_row_selection[i]);
1399 }
1400 }
1401 selected_columns = selected_rows;
1402
1403 AssertDimension(mg_constrained_dofs.size(), selected_rows.size());
1404 edge_constrained_indices.clear();
1405 edge_constrained_indices.resize(selected_rows.size());
1406 edge_constrained_values.clear();
1407 edge_constrained_values.resize(selected_rows.size());
1408
1409 data = data_;
1410
1411 for (unsigned int j = 0; j < selected_rows.size(); ++j)
1412 {
1413 if (data_->n_cell_batches() > 0)
1414 {
1415 AssertDimension(level, data_->get_cell_iterator(0, 0, j)->level());
1416 }
1417
1418 // setup edge_constrained indices
1419 const std::vector<types::global_dof_index> interface_indices =
1420 mg_constrained_dofs[j]
1421 .get_refinement_edge_indices(level)
1422 .get_index_vector();
1423 edge_constrained_indices[j].clear();
1424 edge_constrained_indices[j].reserve(interface_indices.size());
1425 edge_constrained_values[j].resize(interface_indices.size());
1426 const IndexSet &locally_owned =
1427 data->get_dof_handler(selected_rows[j]).locally_owned_mg_dofs(level);
1428 for (const auto interface_index : interface_indices)
1429 if (locally_owned.is_element(interface_index))
1430 edge_constrained_indices[j].push_back(
1431 locally_owned.index_within_set(interface_index));
1432 have_interface_matrices |=
1433 Utilities::MPI::max(
1434 static_cast<unsigned int>(edge_constrained_indices[j].size()),
1435 data->get_vector_partitioner()->get_mpi_communicator()) > 0;
1436 }
1437 }
1438
1439
1440
1441 template <int dim, typename VectorType, typename VectorizedArrayType>
1442 void
1443 Base<dim, VectorType, VectorizedArrayType>::set_constrained_entries_to_one(
1444 VectorType &dst) const
1445 {
1446 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1447 {
1448 const std::vector<unsigned int> &constrained_dofs =
1449 data->get_constrained_dofs(selected_rows[j]);
1450 for (const auto constrained_dof : constrained_dofs)
1451 BlockHelper::subblock(dst, j).local_element(constrained_dof) = 1.;
1452 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1453 BlockHelper::subblock(dst, j).local_element(
1454 edge_constrained_indices[j][i]) = 1.;
1455 }
1456 }
1457
1458
1459
1460 template <int dim, typename VectorType, typename VectorizedArrayType>
1461 void
1462 Base<dim, VectorType, VectorizedArrayType>::vmult(VectorType &dst,
1463 const VectorType &src) const
1464 {
1465 using Number =
1466 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1467 dst = Number(0.);
1468 vmult_add(dst, src);
1469 }
1470
1471
1472
1473 template <int dim, typename VectorType, typename VectorizedArrayType>
1474 void
1475 Base<dim, VectorType, VectorizedArrayType>::vmult_add(
1476 VectorType &dst,
1477 const VectorType &src) const
1478 {
1479 mult_add(dst, src, false);
1480 }
1481
1482
1483
1484 template <int dim, typename VectorType, typename VectorizedArrayType>
1485 void
1486 Base<dim, VectorType, VectorizedArrayType>::Tvmult_add(
1487 VectorType &dst,
1488 const VectorType &src) const
1489 {
1490 mult_add(dst, src, true);
1491 }
1492
1493
1494
1495 template <int dim, typename VectorType, typename VectorizedArrayType>
1496 void
1497 Base<dim, VectorType, VectorizedArrayType>::adjust_ghost_range_if_necessary(
1498 const VectorType &src,
1499 const bool is_row) const
1500 {
1501 using Number =
1502 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1503 for (unsigned int i = 0; i < BlockHelper::n_blocks(src); ++i)
1504 {
1505 const unsigned int mf_component =
1506 is_row ? selected_rows[i] : selected_columns[i];
1507 // If both vectors use the same partitioner -> done
1508 if (BlockHelper::subblock(src, i).get_partitioner().get() ==
1509 data->get_dof_info(mf_component).vector_partitioner.get())
1510 continue;
1511
1512 // If not, assert that the local ranges are the same and reset to the
1513 // current partitioner
1514 Assert(BlockHelper::subblock(src, i)
1515 .get_partitioner()
1516 ->locally_owned_size() ==
1517 data->get_dof_info(mf_component)
1518 .vector_partitioner->locally_owned_size(),
1519 ExcMessage(
1520 "The vector passed to the vmult() function does not have "
1521 "the correct size for compatibility with MatrixFree."));
1522
1523 // copy the vector content to a temporary vector so that it does not get
1524 // lost
1525 LinearAlgebra::distributed::Vector<Number> copy_vec(
1526 BlockHelper::subblock(src, i));
1527 this->data->initialize_dof_vector(
1528 BlockHelper::subblock(const_cast<VectorType &>(src), i),
1529 mf_component);
1530 BlockHelper::subblock(const_cast<VectorType &>(src), i)
1531 .copy_locally_owned_data_from(copy_vec);
1532 }
1533 }
1534
1535
1536
1537 template <int dim, typename VectorType, typename VectorizedArrayType>
1538 void
1539 Base<dim, VectorType, VectorizedArrayType>::preprocess_constraints(
1540 VectorType &dst,
1541 const VectorType &src) const
1542 {
1543 using Number =
1544 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1545 adjust_ghost_range_if_necessary(src, false);
1546 adjust_ghost_range_if_necessary(dst, true);
1547
1548 // set zero Dirichlet values on the input vector (and remember the src and
1549 // dst values because we need to reset them at the end)
1550 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1551 {
1552 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1553 {
1554 edge_constrained_values[j][i] = std::pair<Number, Number>(
1555 BlockHelper::subblock(src, j).local_element(
1556 edge_constrained_indices[j][i]),
1557 BlockHelper::subblock(dst, j).local_element(
1558 edge_constrained_indices[j][i]));
1559 BlockHelper::subblock(const_cast<VectorType &>(src), j)
1560 .local_element(edge_constrained_indices[j][i]) = 0.;
1561 }
1562 }
1563 }
1564
1565
1566
1567 template <int dim, typename VectorType, typename VectorizedArrayType>
1568 void
1569 Base<dim, VectorType, VectorizedArrayType>::mult_add(
1570 VectorType &dst,
1571 const VectorType &src,
1572 const bool transpose) const
1573 {
1574 AssertDimension(dst.size(), src.size());
1575 AssertDimension(BlockHelper::n_blocks(dst), BlockHelper::n_blocks(src));
1576 AssertDimension(BlockHelper::n_blocks(dst), selected_rows.size());
1577 preprocess_constraints(dst, src);
1578 if (transpose)
1579 Tapply_add(dst, src);
1580 else
1581 apply_add(dst, src);
1582 postprocess_constraints(dst, src);
1583 }
1584
1585
1586
1587 template <int dim, typename VectorType, typename VectorizedArrayType>
1588 void
1589 Base<dim, VectorType, VectorizedArrayType>::postprocess_constraints(
1590 VectorType &dst,
1591 const VectorType &src) const
1592 {
1593 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1594 {
1595 const std::vector<unsigned int> &constrained_dofs =
1596 data->get_constrained_dofs(selected_rows[j]);
1597 for (const auto constrained_dof : constrained_dofs)
1598 BlockHelper::subblock(dst, j).local_element(constrained_dof) +=
1599 BlockHelper::subblock(src, j).local_element(constrained_dof);
1600 }
1601
1602 // reset edge constrained values, multiply by unit matrix and add into
1603 // destination
1604 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1605 {
1606 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1607 {
1608 BlockHelper::subblock(const_cast<VectorType &>(src), j)
1609 .local_element(edge_constrained_indices[j][i]) =
1610 edge_constrained_values[j][i].first;
1611 BlockHelper::subblock(dst, j).local_element(
1612 edge_constrained_indices[j][i]) =
1613 edge_constrained_values[j][i].second +
1614 edge_constrained_values[j][i].first;
1615 }
1616 }
1617 }
1618
1619
1620
1621 template <int dim, typename VectorType, typename VectorizedArrayType>
1622 void
1623 Base<dim, VectorType, VectorizedArrayType>::vmult_interface_down(
1624 VectorType &dst,
1625 const VectorType &src) const
1626 {
1627 using Number =
1628 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1629 AssertDimension(dst.size(), src.size());
1630 adjust_ghost_range_if_necessary(src, false);
1631 adjust_ghost_range_if_necessary(dst, true);
1632
1633 dst = Number(0.);
1634
1635 if (!have_interface_matrices)
1636 return;
1637
1638 // set zero Dirichlet values on the input vector (and remember the src and
1639 // dst values because we need to reset them at the end)
1640 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1641 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1642 {
1643 edge_constrained_values[j][i] = std::pair<Number, Number>(
1644 BlockHelper::subblock(src, j).local_element(
1645 edge_constrained_indices[j][i]),
1646 BlockHelper::subblock(dst, j).local_element(
1647 edge_constrained_indices[j][i]));
1648 BlockHelper::subblock(const_cast<VectorType &>(src), j)
1649 .local_element(edge_constrained_indices[j][i]) = 0.;
1650 }
1651
1652 apply_add(dst, src);
1653
1654 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1655 {
1656 unsigned int c = 0;
1657 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1658 {
1659 for (; c < edge_constrained_indices[j][i]; ++c)
1660 BlockHelper::subblock(dst, j).local_element(c) = 0.;
1661 ++c;
1662
1663 // reset the src values
1664 BlockHelper::subblock(const_cast<VectorType &>(src), j)
1665 .local_element(edge_constrained_indices[j][i]) =
1666 edge_constrained_values[j][i].first;
1667 }
1668 for (; c < BlockHelper::subblock(dst, j).locally_owned_size(); ++c)
1669 BlockHelper::subblock(dst, j).local_element(c) = 0.;
1670 }
1671 }
1672
1673
1674
1675 template <int dim, typename VectorType, typename VectorizedArrayType>
1676 void
1677 Base<dim, VectorType, VectorizedArrayType>::vmult_interface_up(
1678 VectorType &dst,
1679 const VectorType &src) const
1680 {
1681 using Number =
1682 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1683 AssertDimension(dst.size(), src.size());
1684 adjust_ghost_range_if_necessary(src, false);
1685 adjust_ghost_range_if_necessary(dst, true);
1686
1687 dst = Number(0.);
1688
1689 if (!have_interface_matrices)
1690 return;
1691
1692 VectorType src_cpy(src);
1693 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1694 {
1695 unsigned int c = 0;
1696 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1697 {
1698 for (; c < edge_constrained_indices[j][i]; ++c)
1699 BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1700 ++c;
1701 }
1702 for (; c < BlockHelper::subblock(src_cpy, j).locally_owned_size(); ++c)
1703 BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1704 }
1705
1706 apply_add(dst, src_cpy);
1707
1708 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1709 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1710 BlockHelper::subblock(dst, j).local_element(
1711 edge_constrained_indices[j][i]) = 0.;
1712 }
1713
1714
1715
1716 template <int dim, typename VectorType, typename VectorizedArrayType>
1717 void
1718 Base<dim, VectorType, VectorizedArrayType>::Tvmult(
1719 VectorType &dst,
1720 const VectorType &src) const
1721 {
1722 using Number =
1723 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1724 dst = Number(0.);
1725 Tvmult_add(dst, src);
1726 }
1727
1728
1729
1730 template <int dim, typename VectorType, typename VectorizedArrayType>
1731 std::size_t
1732 Base<dim, VectorType, VectorizedArrayType>::memory_consumption() const
1733 {
1734 return inverse_diagonal_entries.get() != nullptr ?
1735 inverse_diagonal_entries->memory_consumption() :
1736 sizeof(*this);
1737 }
1738
1739
1740
1741 template <int dim, typename VectorType, typename VectorizedArrayType>
1742 std::shared_ptr<const MatrixFree<
1743 dim,
1744 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
1745 VectorizedArrayType>>
1750
1751
1752
1753 template <int dim, typename VectorType, typename VectorizedArrayType>
1754 const std::shared_ptr<DiagonalMatrix<VectorType>> &
1755 Base<dim, VectorType, VectorizedArrayType>::get_matrix_diagonal_inverse()
1756 const
1757 {
1758 Assert(inverse_diagonal_entries.get() != nullptr &&
1759 inverse_diagonal_entries->m() > 0,
1760 ExcNotInitialized());
1761 return inverse_diagonal_entries;
1762 }
1763
1764
1765
1766 template <int dim, typename VectorType, typename VectorizedArrayType>
1767 const std::shared_ptr<DiagonalMatrix<VectorType>> &
1768 Base<dim, VectorType, VectorizedArrayType>::get_matrix_diagonal() const
1769 {
1770 Assert(diagonal_entries.get() != nullptr && diagonal_entries->m() > 0,
1771 ExcNotInitialized());
1772 return diagonal_entries;
1773 }
1774
1775
1776
1777 template <int dim, typename VectorType, typename VectorizedArrayType>
1778 void
1779 Base<dim, VectorType, VectorizedArrayType>::Tapply_add(
1780 VectorType &dst,
1781 const VectorType &src) const
1782 {
1783 apply_add(dst, src);
1784 }
1785
1786
1787
1788 template <int dim, typename VectorType, typename VectorizedArrayType>
1789 void
1790 Base<dim, VectorType, VectorizedArrayType>::precondition_Jacobi(
1791 VectorType &dst,
1792 const VectorType &src,
1793 const typename Base<dim, VectorType, VectorizedArrayType>::value_type omega)
1794 const
1795 {
1796 Assert(inverse_diagonal_entries.get() && inverse_diagonal_entries->m() > 0,
1797 ExcNotInitialized());
1798 inverse_diagonal_entries->vmult(dst, src);
1799 dst *= omega;
1800 }
1801
1802
1803
1804 //------------------------- MGInterfaceOperator ------------------------------
1805
1806 template <typename OperatorType>
1807 MGInterfaceOperator<OperatorType>::MGInterfaceOperator()
1808 : Subscriptor()
1809 , mf_base_operator(nullptr)
1810 {}
1811
1812
1813
1814 template <typename OperatorType>
1815 void
1816 MGInterfaceOperator<OperatorType>::clear()
1817 {
1818 mf_base_operator = nullptr;
1819 }
1820
1821
1822
1823 template <typename OperatorType>
1824 void
1825 MGInterfaceOperator<OperatorType>::initialize(const OperatorType &operator_in)
1826 {
1827 mf_base_operator = &operator_in;
1828 }
1829
1830
1831
1832 template <typename OperatorType>
1833 template <typename VectorType>
1834 void
1835 MGInterfaceOperator<OperatorType>::vmult(VectorType &dst,
1836 const VectorType &src) const
1837 {
1838#ifndef DEAL_II_MSVC
1839 static_assert(
1840 std::is_same_v<typename VectorType::value_type, value_type>,
1841 "The vector type must be based on the same value type as this "
1842 "operator");
1843#endif
1844
1845 Assert(mf_base_operator != nullptr, ExcNotInitialized());
1846
1847 mf_base_operator->vmult_interface_down(dst, src);
1848 }
1849
1850
1851
1852 template <typename OperatorType>
1853 template <typename VectorType>
1854 void
1855 MGInterfaceOperator<OperatorType>::Tvmult(VectorType &dst,
1856 const VectorType &src) const
1857 {
1858#ifndef DEAL_II_MSVC
1859 static_assert(
1860 std::is_same_v<typename VectorType::value_type, value_type>,
1861 "The vector type must be based on the same value type as this "
1862 "operator");
1863#endif
1864
1865 Assert(mf_base_operator != nullptr, ExcNotInitialized());
1866
1867 mf_base_operator->vmult_interface_up(dst, src);
1868 }
1869
1870
1871
1872 template <typename OperatorType>
1873 template <typename VectorType>
1874 void
1875 MGInterfaceOperator<OperatorType>::initialize_dof_vector(
1876 VectorType &vec) const
1877 {
1878 Assert(mf_base_operator != nullptr, ExcNotInitialized());
1879
1880 mf_base_operator->initialize_dof_vector(vec);
1881 }
1882
1883
1884
1885 //-----------------------------MassOperator----------------------------------
1886
1887 template <int dim,
1888 int fe_degree,
1889 int n_q_points_1d,
1890 int n_components,
1891 typename VectorType,
1892 typename VectorizedArrayType>
1893 MassOperator<dim,
1894 fe_degree,
1895 n_q_points_1d,
1896 n_components,
1897 VectorType,
1898 VectorizedArrayType>::MassOperator()
1899 : Base<dim, VectorType, VectorizedArrayType>()
1900 {
1901 AssertThrow(
1902 IsBlockVector<VectorType>::value == false,
1903 ExcNotImplemented(
1904 "This class only supports the non-blocked vector variant of the Base "
1905 "operator because only a single FEEvaluation object is used in the "
1906 "apply function."));
1907 }
1908
1909
1910
1911 template <int dim,
1912 int fe_degree,
1913 int n_q_points_1d,
1914 int n_components,
1915 typename VectorType,
1916 typename VectorizedArrayType>
1917 void
1918 MassOperator<dim,
1919 fe_degree,
1920 n_q_points_1d,
1921 n_components,
1922 VectorType,
1923 VectorizedArrayType>::compute_diagonal()
1924 {
1925 Assert((Base<dim, VectorType, VectorizedArrayType>::data.get() != nullptr),
1926 ExcNotInitialized());
1927 Assert(this->selected_rows == this->selected_columns,
1928 ExcMessage("This function is only implemented for square (not "
1929 "rectangular) operators."));
1930
1931 this->inverse_diagonal_entries =
1932 std::make_shared<DiagonalMatrix<VectorType>>();
1933 this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1934 VectorType &inverse_diagonal_vector =
1935 this->inverse_diagonal_entries->get_vector();
1936 VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1937 this->initialize_dof_vector(inverse_diagonal_vector);
1938 this->initialize_dof_vector(diagonal_vector);
1939
1940 // Set up the action of the mass matrix in a way that's compatible with
1941 // MatrixFreeTools::compute_diagonal:
1942 auto diagonal_evaluation = [](auto &integrator) {
1943 integrator.evaluate(EvaluationFlags::values);
1944 for (unsigned int q = 0; q < integrator.n_q_points; ++q)
1945 integrator.submit_value(integrator.get_value(q), q);
1946 integrator.integrate(EvaluationFlags::values);
1947 };
1948
1949 std::function<void(
1950 FEEvaluation<
1951 dim,
1952 fe_degree,
1953 n_q_points_1d,
1954 n_components,
1955 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
1956 VectorizedArrayType> &)>
1957 diagonal_evaluation_f(diagonal_evaluation);
1958
1959 Assert(this->selected_rows.size() > 0, ExcInternalError());
1960 for (unsigned int block_n = 0; block_n < this->selected_rows.size();
1961 ++block_n)
1962 MatrixFreeTools::compute_diagonal(*this->data,
1963 BlockHelper::subblock(diagonal_vector,
1964 block_n),
1965 diagonal_evaluation_f,
1966 this->selected_rows[block_n]);
1967
1968 // Constrained entries will create zeros on the main diagonal, which we
1969 // don't want
1970 this->set_constrained_entries_to_one(diagonal_vector);
1971
1972 inverse_diagonal_vector = diagonal_vector;
1973
1974 for (unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
1975 ++i)
1976 {
1977#ifdef DEBUG
1978 // only define the type alias in debug mode to avoid a warning
1979 using Number =
1980 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1981 Assert(diagonal_vector.local_element(i) > Number(0),
1982 ExcInternalError());
1983#endif
1984 inverse_diagonal_vector.local_element(i) =
1985 1. / inverse_diagonal_vector.local_element(i);
1986 }
1987
1988 // We never need ghost values so don't update them
1989 }
1990
1991
1992
1993 template <int dim,
1994 int fe_degree,
1995 int n_q_points_1d,
1996 int n_components,
1997 typename VectorType,
1998 typename VectorizedArrayType>
1999 void
2000 MassOperator<dim,
2001 fe_degree,
2002 n_q_points_1d,
2003 n_components,
2004 VectorType,
2005 VectorizedArrayType>::compute_lumped_diagonal()
2006 {
2007 using Number =
2008 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2009 Assert((Base<dim, VectorType, VectorizedArrayType>::data.get() != nullptr),
2010 ExcNotInitialized());
2011 Assert(this->selected_rows == this->selected_columns,
2012 ExcMessage("This function is only implemented for square (not "
2013 "rectangular) operators."));
2014
2015 inverse_lumped_diagonal_entries =
2016 std::make_shared<DiagonalMatrix<VectorType>>();
2017 lumped_diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
2018 VectorType &inverse_lumped_diagonal_vector =
2019 inverse_lumped_diagonal_entries->get_vector();
2020 VectorType &lumped_diagonal_vector = lumped_diagonal_entries->get_vector();
2021 this->initialize_dof_vector(inverse_lumped_diagonal_vector);
2022 this->initialize_dof_vector(lumped_diagonal_vector);
2023
2024 // Re-use the inverse_lumped_diagonal_vector as the vector of 1s
2025 inverse_lumped_diagonal_vector = Number(1.);
2026 apply_add(lumped_diagonal_vector, inverse_lumped_diagonal_vector);
2027 this->set_constrained_entries_to_one(lumped_diagonal_vector);
2028
2029 const size_type locally_owned_size =
2030 inverse_lumped_diagonal_vector.locally_owned_size();
2031 // A caller may request a lumped diagonal matrix when it doesn't make sense
2032 // (e.g., an element with negative-mean basis functions). Avoid division by
2033 // zero so we don't cause a floating point exception but permit negative
2034 // entries here.
2035 for (size_type i = 0; i < locally_owned_size; ++i)
2036 {
2037 if (lumped_diagonal_vector.local_element(i) == Number(0.))
2038 inverse_lumped_diagonal_vector.local_element(i) = Number(1.);
2039 else
2040 inverse_lumped_diagonal_vector.local_element(i) =
2041 Number(1.) / lumped_diagonal_vector.local_element(i);
2042 }
2043
2044 inverse_lumped_diagonal_vector.update_ghost_values();
2045 lumped_diagonal_vector.update_ghost_values();
2046 }
2047
2048
2049
2050 template <int dim,
2051 int fe_degree,
2052 int n_q_points_1d,
2053 int n_components,
2054 typename VectorType,
2055 typename VectorizedArrayType>
2056 const std::shared_ptr<DiagonalMatrix<VectorType>> &
2057 MassOperator<dim,
2058 fe_degree,
2059 n_q_points_1d,
2060 n_components,
2061 VectorType,
2062 VectorizedArrayType>::get_matrix_lumped_diagonal_inverse() const
2063 {
2064 Assert(inverse_lumped_diagonal_entries.get() != nullptr &&
2065 inverse_lumped_diagonal_entries->m() > 0,
2066 ExcNotInitialized());
2067 return inverse_lumped_diagonal_entries;
2068 }
2069
2070
2071
2072 template <int dim,
2073 int fe_degree,
2074 int n_q_points_1d,
2075 int n_components,
2076 typename VectorType,
2077 typename VectorizedArrayType>
2078 const std::shared_ptr<DiagonalMatrix<VectorType>> &
2079 MassOperator<dim,
2080 fe_degree,
2081 n_q_points_1d,
2082 n_components,
2083 VectorType,
2084 VectorizedArrayType>::get_matrix_lumped_diagonal() const
2085 {
2086 Assert(lumped_diagonal_entries.get() != nullptr &&
2087 lumped_diagonal_entries->m() > 0,
2088 ExcNotInitialized());
2089 return lumped_diagonal_entries;
2090 }
2091
2092
2093
2094 template <int dim,
2095 int fe_degree,
2096 int n_q_points_1d,
2097 int n_components,
2098 typename VectorType,
2099 typename VectorizedArrayType>
2100 void
2101 MassOperator<dim,
2102 fe_degree,
2103 n_q_points_1d,
2104 n_components,
2105 VectorType,
2106 VectorizedArrayType>::apply_add(VectorType &dst,
2107 const VectorType &src) const
2108 {
2109 Base<dim, VectorType, VectorizedArrayType>::data->cell_loop(
2110 &MassOperator::local_apply_cell, this, dst, src);
2111 }
2112
2113
2114
2115 template <int dim,
2116 int fe_degree,
2117 int n_q_points_1d,
2118 int n_components,
2119 typename VectorType,
2120 typename VectorizedArrayType>
2121 void
2122 MassOperator<dim,
2123 fe_degree,
2124 n_q_points_1d,
2125 n_components,
2126 VectorType,
2127 VectorizedArrayType>::
2128 local_apply_cell(
2129 const MatrixFree<
2130 dim,
2131 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2132 VectorizedArrayType> &data,
2133 VectorType &dst,
2134 const VectorType &src,
2135 const std::pair<unsigned int, unsigned int> &cell_range) const
2136 {
2137 using Number =
2138 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2139 FEEvaluation<dim,
2140 fe_degree,
2141 n_q_points_1d,
2142 n_components,
2143 Number,
2144 VectorizedArrayType>
2145 phi(data, this->selected_rows[0]);
2146 for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2147 {
2148 phi.reinit(cell);
2149 phi.read_dof_values(src);
2150 phi.evaluate(EvaluationFlags::values);
2151 for (unsigned int q = 0; q < phi.n_q_points; ++q)
2152 phi.submit_value(phi.get_value(q), q);
2153 phi.integrate(EvaluationFlags::values);
2154 phi.distribute_local_to_global(dst);
2155 }
2156 }
2157
2158
2159 //-----------------------------LaplaceOperator----------------------------------
2160
2161 template <int dim,
2162 int fe_degree,
2163 int n_q_points_1d,
2164 int n_components,
2165 typename VectorType,
2166 typename VectorizedArrayType>
2167 LaplaceOperator<dim,
2168 fe_degree,
2169 n_q_points_1d,
2170 n_components,
2171 VectorType,
2172 VectorizedArrayType>::LaplaceOperator()
2173 : Base<dim, VectorType, VectorizedArrayType>()
2174 {}
2175
2176
2177
2178 template <int dim,
2179 int fe_degree,
2180 int n_q_points_1d,
2181 int n_components,
2182 typename VectorType,
2183 typename VectorizedArrayType>
2184 void
2185 LaplaceOperator<dim,
2186 fe_degree,
2187 n_q_points_1d,
2188 n_components,
2189 VectorType,
2190 VectorizedArrayType>::clear()
2191 {
2192 Base<dim, VectorType, VectorizedArrayType>::clear();
2193 scalar_coefficient.reset();
2194 }
2195
2196
2197
2198 template <int dim,
2199 int fe_degree,
2200 int n_q_points_1d,
2201 int n_components,
2202 typename VectorType,
2203 typename VectorizedArrayType>
2204 void
2205 LaplaceOperator<dim,
2206 fe_degree,
2207 n_q_points_1d,
2208 n_components,
2209 VectorType,
2210 VectorizedArrayType>::
2211 set_coefficient(
2212 const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient_)
2213 {
2214 scalar_coefficient = scalar_coefficient_;
2215 }
2216
2217
2218
2219 template <int dim,
2220 int fe_degree,
2221 int n_q_points_1d,
2222 int n_components,
2223 typename VectorType,
2224 typename VectorizedArrayType>
2225 std::shared_ptr<Table<2, VectorizedArrayType>>
2226 LaplaceOperator<dim,
2227 fe_degree,
2228 n_q_points_1d,
2229 n_components,
2230 VectorType,
2231 VectorizedArrayType>::get_coefficient()
2232 {
2233 Assert(scalar_coefficient.get(), ExcNotInitialized());
2234 return scalar_coefficient;
2235 }
2236
2237
2238
2239 template <int dim,
2240 int fe_degree,
2241 int n_q_points_1d,
2242 int n_components,
2243 typename VectorType,
2244 typename VectorizedArrayType>
2245 void
2246 LaplaceOperator<dim,
2247 fe_degree,
2248 n_q_points_1d,
2249 n_components,
2250 VectorType,
2251 VectorizedArrayType>::compute_diagonal()
2252 {
2253 using Number =
2254 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2255 Assert((Base<dim, VectorType, VectorizedArrayType>::data.get() != nullptr),
2256 ExcNotInitialized());
2257
2258 this->inverse_diagonal_entries =
2259 std::make_shared<DiagonalMatrix<VectorType>>();
2260 this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
2261 VectorType &inverse_diagonal_vector =
2262 this->inverse_diagonal_entries->get_vector();
2263 VectorType &diagonal_vector = this->diagonal_entries->get_vector();
2264 this->initialize_dof_vector(inverse_diagonal_vector);
2265 this->initialize_dof_vector(diagonal_vector);
2266
2267 this->data->cell_loop(&LaplaceOperator::local_diagonal_cell,
2268 this,
2269 diagonal_vector,
2270 /*unused*/ diagonal_vector);
2271 this->set_constrained_entries_to_one(diagonal_vector);
2272
2273 inverse_diagonal_vector = diagonal_vector;
2274
2275 for (unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
2276 ++i)
2277 if (std::abs(inverse_diagonal_vector.local_element(i)) >
2278 std::sqrt(std::numeric_limits<Number>::epsilon()))
2279 inverse_diagonal_vector.local_element(i) =
2280 1. / inverse_diagonal_vector.local_element(i);
2281 else
2282 inverse_diagonal_vector.local_element(i) = 1.;
2283
2284 // We never need ghost values so don't update them
2285 }
2286
2287
2288
2289 template <int dim,
2290 int fe_degree,
2291 int n_q_points_1d,
2292 int n_components,
2293 typename VectorType,
2294 typename VectorizedArrayType>
2295 void
2296 LaplaceOperator<dim,
2297 fe_degree,
2298 n_q_points_1d,
2299 n_components,
2300 VectorType,
2301 VectorizedArrayType>::apply_add(VectorType &dst,
2302 const VectorType &src) const
2303 {
2304 Base<dim, VectorType, VectorizedArrayType>::data->cell_loop(
2305 &LaplaceOperator::local_apply_cell, this, dst, src);
2306 }
2307
2308 namespace Implementation
2309 {
2310 template <typename VectorizedArrayType>
2311 bool
2312 non_negative(const VectorizedArrayType &n)
2313 {
2314 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2315 if (n[v] < 0.)
2316 return false;
2317
2318 return true;
2319 }
2320 } // namespace Implementation
2321
2322
2323
2324 template <int dim,
2325 int fe_degree,
2326 int n_q_points_1d,
2327 int n_components,
2328 typename VectorType,
2329 typename VectorizedArrayType>
2330 template <int n_components_compute>
2331 void
2332 LaplaceOperator<dim,
2333 fe_degree,
2334 n_q_points_1d,
2335 n_components,
2336 VectorType,
2337 VectorizedArrayType>::
2338 do_operation_on_cell(
2339 FEEvaluation<
2340 dim,
2341 fe_degree,
2342 n_q_points_1d,
2343 n_components_compute,
2344 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2345 VectorizedArrayType> &phi,
2346 const unsigned int cell) const
2347 {
2348 phi.evaluate(EvaluationFlags::gradients);
2349 if (scalar_coefficient.get())
2350 {
2351 Assert(scalar_coefficient->size(1) == 1 ||
2352 scalar_coefficient->size(1) == phi.n_q_points,
2353 ExcMessage("The number of columns in the coefficient table must "
2354 "be either 1 or the number of quadrature points " +
2355 std::to_string(phi.n_q_points) +
2356 ", but the given value was " +
2357 std::to_string(scalar_coefficient->size(1))));
2358 if (scalar_coefficient->size(1) == phi.n_q_points)
2359 for (unsigned int q = 0; q < phi.n_q_points; ++q)
2360 {
2361 Assert(Implementation::non_negative(
2362 (*scalar_coefficient)(cell, q)),
2363 ExcMessage("Coefficient must be non-negative"));
2364 phi.submit_gradient((*scalar_coefficient)(cell, q) *
2365 phi.get_gradient(q),
2366 q);
2367 }
2368 else
2369 {
2370 Assert(Implementation::non_negative((*scalar_coefficient)(cell, 0)),
2371 ExcMessage("Coefficient must be non-negative"));
2372 const VectorizedArrayType coefficient =
2373 (*scalar_coefficient)(cell, 0);
2374 for (unsigned int q = 0; q < phi.n_q_points; ++q)
2375 phi.submit_gradient(coefficient * phi.get_gradient(q), q);
2376 }
2377 }
2378 else
2379 {
2380 for (unsigned int q = 0; q < phi.n_q_points; ++q)
2381 {
2382 phi.submit_gradient(phi.get_gradient(q), q);
2383 }
2384 }
2385 phi.integrate(EvaluationFlags::gradients);
2386 }
2387
2388
2389
2390 template <int dim,
2391 int fe_degree,
2392 int n_q_points_1d,
2393 int n_components,
2394 typename VectorType,
2395 typename VectorizedArrayType>
2396 void
2397 LaplaceOperator<dim,
2398 fe_degree,
2399 n_q_points_1d,
2400 n_components,
2401 VectorType,
2402 VectorizedArrayType>::
2403 local_apply_cell(
2404 const MatrixFree<
2405 dim,
2406 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2407 VectorizedArrayType> &data,
2408 VectorType &dst,
2409 const VectorType &src,
2410 const std::pair<unsigned int, unsigned int> &cell_range) const
2411 {
2412 using Number =
2413 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2414 FEEvaluation<dim,
2415 fe_degree,
2416 n_q_points_1d,
2417 n_components,
2418 Number,
2419 VectorizedArrayType>
2420 phi(data, this->selected_rows[0]);
2421 for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2422 {
2423 phi.reinit(cell);
2424 phi.read_dof_values(src);
2425 do_operation_on_cell(phi, cell);
2426 phi.distribute_local_to_global(dst);
2427 }
2428 }
2429
2430
2431 template <int dim,
2432 int fe_degree,
2433 int n_q_points_1d,
2434 int n_components,
2435 typename VectorType,
2436 typename VectorizedArrayType>
2437 void
2438 LaplaceOperator<dim,
2439 fe_degree,
2440 n_q_points_1d,
2441 n_components,
2442 VectorType,
2443 VectorizedArrayType>::
2444 local_diagonal_cell(
2445 const MatrixFree<
2446 dim,
2447 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2448 VectorizedArrayType> &data,
2449 VectorType &dst,
2450 const VectorType &,
2451 const std::pair<unsigned int, unsigned int> &cell_range) const
2452 {
2453 using Number =
2454 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2455
2456 FEEvaluation<dim, fe_degree, n_q_points_1d, 1, Number, VectorizedArrayType>
2457 eval(data, this->selected_rows[0]);
2458 FEEvaluation<dim,
2459 fe_degree,
2460 n_q_points_1d,
2461 n_components,
2462 Number,
2463 VectorizedArrayType>
2464 eval_vector(data, this->selected_rows[0]);
2465 for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2466 {
2467 eval.reinit(cell);
2468 eval_vector.reinit(cell);
2469 // This function assumes that we have the same result on all
2470 // components, so we only need to go through the columns of one scalar
2471 // component, for which we have created a separate evaluator (attached
2472 // to the first component, but the component does not matter because
2473 // we only use the underlying integrals)
2474 for (unsigned int i = 0; i < eval.dofs_per_cell; ++i)
2475 {
2476 for (unsigned int j = 0; j < eval.dofs_per_cell; ++j)
2477 eval.begin_dof_values()[j] = VectorizedArrayType();
2478 eval.begin_dof_values()[i] = 1.;
2479
2480 do_operation_on_cell(eval, cell);
2481
2482 // We now pick up the value on the diagonal (row i) and broadcast
2483 // it to a second evaluator for all vector components, which we
2484 // will distribute to the result vector afterwards
2485 for (unsigned int c = 0; c < n_components; ++c)
2486 eval_vector
2487 .begin_dof_values()[i + c * eval_vector.dofs_per_component] =
2488 eval.begin_dof_values()[i];
2489 }
2490 eval_vector.distribute_local_to_global(dst);
2491 }
2492 }
2493
2494
2495} // end of namespace MatrixFreeOperators
2496
2497
2498DEAL_II_NAMESPACE_CLOSE
2499
2500#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:1749
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:1486
MatrixFreeOperators::Base::set_constrained_entries_to_one
void set_constrained_entries_to_one(VectorType &dst) const
Definition operators.h:1443
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:1475
MatrixFreeOperators::Base::vmult_interface_down
void vmult_interface_down(VectorType &dst, const VectorType &src) const
Definition operators.h:1623
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:1539
MatrixFreeOperators::Base::mult_add
void mult_add(VectorType &dst, const VectorType &src, const bool transpose) const
Definition operators.h:1569
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:1768
MatrixFreeOperators::Base::Tvmult
void Tvmult(VectorType &dst, const VectorType &src) const
Definition operators.h:1718
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:1304
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:1746
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:1373
MatrixFreeOperators::Base::get_matrix_diagonal_inverse
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal_inverse() const
Definition operators.h:1755
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:1357
MatrixFreeOperators::Base::initialize_dof_vector
void initialize_dof_vector(VectorType &vec) const
Definition operators.h:1279
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:1779
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:1497
MatrixFreeOperators::Base::memory_consumption
virtual std::size_t memory_consumption() const
Definition operators.h:1732
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:1790
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:1462
MatrixFreeOperators::Base::vmult_interface_up
void vmult_interface_up(VectorType &dst, const VectorType &src) const
Definition operators.h:1677
MatrixFreeOperators::Base::postprocess_constraints
void postprocess_constraints(VectorType &dst, const VectorType &src) const
Definition operators.h:1589
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::CellwiseInverseMassMatrix
CellwiseInverseMassMatrix(const FEEvaluationBase< dim, n_components, Number, false, VectorizedArrayType > &fe_eval)
Definition operators.h:1030
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::value_type
typename Base< dim, VectorType, VectorizedArrayType >::value_type value_type
Definition operators.h:875
MatrixFreeOperators::LaplaceOperator::get_coefficient
std::shared_ptr< Table< 2, VectorizedArrayType > > get_coefficient()
Definition operators.h:2231
MatrixFreeOperators::LaplaceOperator::compute_diagonal
virtual void compute_diagonal() override
Definition operators.h:2251
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:2444
MatrixFreeOperators::LaplaceOperator::size_type
typename Base< dim, VectorType, VectorizedArrayType >::size_type size_type
Definition operators.h:881
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:2403
MatrixFreeOperators::LaplaceOperator::apply_add
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition operators.h:2301
MatrixFreeOperators::LaplaceOperator::set_coefficient
void set_coefficient(const std::shared_ptr< Table< 2, VectorizedArrayType > > &scalar_coefficient)
Definition operators.h:2211
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:2190
MatrixFreeOperators::MGInterfaceOperator::value_type
typename OperatorType::value_type value_type
Definition operators.h:543
MatrixFreeOperators::MGInterfaceOperator::mf_base_operator
SmartPointer< const OperatorType > mf_base_operator
Definition operators.h:593
MatrixFreeOperators::MGInterfaceOperator::vmult
void vmult(VectorType &dst, const VectorType &src) const
Definition operators.h:1835
MatrixFreeOperators::MGInterfaceOperator::initialize
void initialize(const OperatorType &operator_in)
Definition operators.h:1825
MatrixFreeOperators::MGInterfaceOperator::Tvmult
void Tvmult(VectorType &dst, const VectorType &src) const
Definition operators.h:1855
MatrixFreeOperators::MGInterfaceOperator::initialize_dof_vector
void initialize_dof_vector(VectorType &vec) const
Definition operators.h:1875
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::apply_add
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition operators.h:2106
MatrixFreeOperators::MassOperator::size_type
typename Base< dim, VectorType, VectorizedArrayType >::size_type size_type
Definition operators.h:773
MatrixFreeOperators::MassOperator::get_matrix_lumped_diagonal
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_lumped_diagonal() const
Definition operators.h:2084
MatrixFreeOperators::MassOperator::get_matrix_lumped_diagonal_inverse
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_lumped_diagonal_inverse() const
Definition operators.h:2062
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:2128
MatrixFreeOperators::MassOperator::compute_diagonal
virtual void compute_diagonal() override
Definition operators.h:1923
MatrixFreeOperators::MassOperator::value_type
typename Base< dim, VectorType, VectorizedArrayType >::value_type value_type
Definition operators.h:767
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:503
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition derivative_form.h:610
diagonal_matrix.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:1638
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:1985
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exceptions.h:2053
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:1739
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: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:2312
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)
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.11.0