Reference documentation for deal.II version 9.4.1
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
operators.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2011 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
17#ifndef dealii_matrix_free_operators_h
18#define dealii_matrix_free_operators_h
19
20
21#include <deal.II/base/config.h>
22
23#include <deal.II/base/exceptions.h>
24#include <deal.II/base/subscriptor.h>
25#include <deal.II/base/vectorization.h>
26
27#include <deal.II/lac/diagonal_matrix.h>
28#include <deal.II/lac/la_parallel_vector.h>
29
30#include <deal.II/matrix_free/fe_evaluation.h>
31#include <deal.II/matrix_free/matrix_free.h>
32#include <deal.II/matrix_free/tools.h>
33
34#include <deal.II/multigrid/mg_constrained_dofs.h>
35
36
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 typename std::enable_if<IsBlockVector<VectorType>::value,
49 unsigned int>::type
50 n_blocks(const VectorType &vector)
51 {
52 return vector.n_blocks();
53 }
54
55 template <typename VectorType>
56 typename std::enable_if<!IsBlockVector<VectorType>::value,
57 unsigned int>::type
58 n_blocks(const VectorType &)
59 {
60 return 1;
61 }
62
63 template <typename VectorType>
64 typename std::enable_if<IsBlockVector<VectorType>::value,
65 typename VectorType::BlockType &>::type
66 subblock(VectorType &vector, unsigned int block_no)
67 {
68 return vector.block(block_no);
69 }
70
71 template <typename VectorType>
72 typename std::enable_if<IsBlockVector<VectorType>::value,
73 const typename VectorType::BlockType &>::type
74 subblock(const VectorType &vector, unsigned int block_no)
75 {
76 AssertIndexRange(block_no, vector.n_blocks());
77 return vector.block(block_no);
78 }
79
80 template <typename VectorType>
81 typename std::enable_if<!IsBlockVector<VectorType>::value,
82 VectorType &>::type
83 subblock(VectorType &vector, unsigned int)
84 {
85 return vector;
86 }
87
88 template <typename VectorType>
89 typename std::enable_if<!IsBlockVector<VectorType>::value,
90 const VectorType &>::type
91 subblock(const VectorType &vector, unsigned int)
92 {
93 return vector;
94 }
95
96 template <typename VectorType>
97 typename std::enable_if<IsBlockVector<VectorType>::value, void>::type
98 collect_sizes(VectorType &vector)
99 {
100 vector.collect_sizes();
101 }
102
103 template <typename VectorType>
104 typename std::enable_if<!IsBlockVector<VectorType>::value, void>::type
105 collect_sizes(const VectorType &)
106 {}
107 } // namespace BlockHelper
108
186 template <int dim,
187 typename VectorType = LinearAlgebra::distributed::Vector<double>,
188 typename VectorizedArrayType =
189 VectorizedArray<typename VectorType::value_type>>
190 class Base : public Subscriptor
191 {
192 public:
196 using value_type = typename VectorType::value_type;
197
201 using size_type = typename VectorType::size_type;
202
206 Base();
207
211 virtual ~Base() override = default;
212
217 virtual void
218 clear();
219
236 void
237 initialize(std::shared_ptr<
238 const MatrixFree<dim, value_type, VectorizedArrayType>> data,
239 const std::vector<unsigned int> &selected_row_blocks =
240 std::vector<unsigned int>(),
241 const std::vector<unsigned int> &selected_column_blocks =
242 std::vector<unsigned int>());
243
257 void
258 initialize(std::shared_ptr<
259 const MatrixFree<dim, value_type, VectorizedArrayType>> data,
260 const MGConstrainedDoFs & mg_constrained_dofs,
261 const unsigned int level,
262 const std::vector<unsigned int> &selected_row_blocks =
263 std::vector<unsigned int>());
264
279 void
280 initialize(std::shared_ptr<
281 const MatrixFree<dim, value_type, VectorizedArrayType>> data_,
282 const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
283 const unsigned int level,
284 const std::vector<unsigned int> & selected_row_blocks =
285 std::vector<unsigned int>());
286
290 size_type
291 m() const;
292
296 size_type
297 n() const;
298
302 void
303 vmult_interface_down(VectorType &dst, const VectorType &src) const;
304
308 void
309 vmult_interface_up(VectorType &dst, const VectorType &src) const;
310
314 void
315 vmult(VectorType &dst, const VectorType &src) const;
316
320 void
321 Tvmult(VectorType &dst, const VectorType &src) const;
322
326 void
327 vmult_add(VectorType &dst, const VectorType &src) const;
328
332 void
333 Tvmult_add(VectorType &dst, const VectorType &src) const;
334
339 value_type
340 el(const unsigned int row, const unsigned int col) const;
341
346 virtual std::size_t
347 memory_consumption() const;
348
352 void
353 initialize_dof_vector(VectorType &vec) const;
354
367 virtual void
368 compute_diagonal() = 0;
369
373 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
374 get_matrix_free() const;
375
379 const std::shared_ptr<DiagonalMatrix<VectorType>> &
380 get_matrix_diagonal_inverse() const;
381
385 const std::shared_ptr<DiagonalMatrix<VectorType>> &
386 get_matrix_diagonal() const;
387
393 void
394 precondition_Jacobi(VectorType & dst,
395 const VectorType &src,
396 const value_type omega) const;
397
398 protected:
403 void
404 preprocess_constraints(VectorType &dst, const VectorType &src) const;
405
410 void
411 postprocess_constraints(VectorType &dst, const VectorType &src) const;
412
417 void
418 set_constrained_entries_to_one(VectorType &dst) const;
419
423 virtual void
424 apply_add(VectorType &dst, const VectorType &src) const = 0;
425
431 virtual void
432 Tapply_add(VectorType &dst, const VectorType &src) const;
433
437 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
438 data;
439
444 std::shared_ptr<DiagonalMatrix<VectorType>> diagonal_entries;
445
450 std::shared_ptr<DiagonalMatrix<VectorType>> inverse_diagonal_entries;
451
456 std::vector<unsigned int> selected_rows;
457
462 std::vector<unsigned int> selected_columns;
463
464 private:
468 std::vector<std::vector<unsigned int>> edge_constrained_indices;
469
473 mutable std::vector<std::vector<std::pair<value_type, value_type>>>
474 edge_constrained_values;
475
480 bool have_interface_matrices;
481
486 void
487 mult_add(VectorType & dst,
488 const VectorType &src,
489 const bool transpose) const;
490
498 void
499 adjust_ghost_range_if_necessary(const VectorType &vec,
500 const bool is_row) const;
501 };
502
503
504
539 template <typename OperatorType>
540 class MGInterfaceOperator : public Subscriptor
541 {
542 public:
546 using value_type = typename OperatorType::value_type;
547
551 using size_type = typename OperatorType::size_type;
552
556 MGInterfaceOperator();
557
561 void
562 clear();
563
567 void
568 initialize(const OperatorType &operator_in);
569
573 template <typename VectorType>
574 void
575 vmult(VectorType &dst, const VectorType &src) const;
576
580 template <typename VectorType>
581 void
582 Tvmult(VectorType &dst, const VectorType &src) const;
583
587 template <typename VectorType>
588 void
589 initialize_dof_vector(VectorType &vec) const;
590
591
592 private:
596 SmartPointer<const OperatorType> mf_base_operator;
597 };
598
599
600
618 template <int dim,
619 int fe_degree,
620 int n_components = 1,
621 typename Number = double,
622 typename VectorizedArrayType = VectorizedArray<Number>>
623 class CellwiseInverseMassMatrix
624 {
625 static_assert(
626 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
627 "Type of Number and of VectorizedArrayType do not match.");
628
629 public:
634 CellwiseInverseMassMatrix(
635 const FEEvaluationBase<dim,
636 n_components,
637 Number,
638 false,
639 VectorizedArrayType> &fe_eval);
640
649 void
650 apply(const AlignedVector<VectorizedArrayType> &inverse_coefficient,
651 const unsigned int n_actual_components,
652 const VectorizedArrayType * in_array,
653 VectorizedArrayType * out_array) const;
654
666 void
667 apply(const VectorizedArrayType *in_array,
668 VectorizedArrayType * out_array) const;
669
703 void
704 transform_from_q_points_to_basis(const unsigned int n_actual_components,
705 const VectorizedArrayType *in_array,
706 VectorizedArrayType *out_array) const;
707
713 void
714 fill_inverse_JxW_values(
715 AlignedVector<VectorizedArrayType> &inverse_jxw) const;
716
717 private:
721 const FEEvaluationBase<dim,
722 n_components,
723 Number,
724 false,
725 VectorizedArrayType> &fe_eval;
726 };
727
728
729
737 template <int dim,
738 int fe_degree,
739 int n_q_points_1d = fe_degree + 1,
740 int n_components = 1,
741 typename VectorType = LinearAlgebra::distributed::Vector<double>,
742 typename VectorizedArrayType =
743 VectorizedArray<typename VectorType::value_type>>
744 class MassOperator : public Base<dim, VectorType, VectorizedArrayType>
745 {
746 public:
750 using value_type =
751 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
752
756 using size_type =
757 typename Base<dim, VectorType, VectorizedArrayType>::size_type;
758
762 MassOperator();
763
767 virtual void
768 compute_diagonal() override;
769
785 void
786 compute_lumped_diagonal();
787
791 const std::shared_ptr<DiagonalMatrix<VectorType>> &
792 get_matrix_lumped_diagonal() const;
793
797 const std::shared_ptr<DiagonalMatrix<VectorType>> &
798 get_matrix_lumped_diagonal_inverse() const;
799
800 private:
806 virtual void
807 apply_add(VectorType &dst, const VectorType &src) const override;
808
812 void
813 local_apply_cell(
814 const MatrixFree<dim, value_type, VectorizedArrayType> &data,
815 VectorType & dst,
816 const VectorType & src,
817 const std::pair<unsigned int, unsigned int> & cell_range) const;
818
823 std::shared_ptr<DiagonalMatrix<VectorType>> lumped_diagonal_entries;
824
829 std::shared_ptr<DiagonalMatrix<VectorType>> inverse_lumped_diagonal_entries;
830 };
831
832
833
844 template <int dim,
845 int fe_degree,
846 int n_q_points_1d = fe_degree + 1,
847 int n_components = 1,
848 typename VectorType = LinearAlgebra::distributed::Vector<double>,
849 typename VectorizedArrayType =
850 VectorizedArray<typename VectorType::value_type>>
851 class LaplaceOperator : public Base<dim, VectorType, VectorizedArrayType>
852 {
853 public:
857 using value_type =
858 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
859
863 using size_type =
864 typename Base<dim, VectorType, VectorizedArrayType>::size_type;
865
869 LaplaceOperator();
870
877 virtual void
878 compute_diagonal() override;
879
930 void
931 set_coefficient(
932 const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient);
933
938 virtual void
939 clear() override;
940
947 std::shared_ptr<Table<2, VectorizedArrayType>>
948 get_coefficient();
949
950 private:
956 virtual void
957 apply_add(VectorType &dst, const VectorType &src) const override;
958
962 void
963 local_apply_cell(
964 const MatrixFree<dim, value_type, VectorizedArrayType> &data,
965 VectorType & dst,
966 const VectorType & src,
967 const std::pair<unsigned int, unsigned int> & cell_range) const;
968
972 void
973 local_diagonal_cell(
974 const MatrixFree<dim, value_type, VectorizedArrayType> &data,
975 VectorType & dst,
976 const VectorType &,
977 const std::pair<unsigned int, unsigned int> &cell_range) const;
978
982 void
983 do_operation_on_cell(
984 FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, value_type>
985 & phi,
986 const unsigned int cell) const;
987
991 std::shared_ptr<Table<2, VectorizedArrayType>> scalar_coefficient;
992 };
993
994
995
996 // ------------------------------------ inline functions ---------------------
997
998 template <int dim,
999 int fe_degree,
1000 int n_components,
1001 typename Number,
1002 typename VectorizedArrayType>
1003 inline CellwiseInverseMassMatrix<dim,
1004 fe_degree,
1005 n_components,
1006 Number,
1007 VectorizedArrayType>::
1008 CellwiseInverseMassMatrix(
1009 const FEEvaluationBase<dim,
1010 n_components,
1011 Number,
1012 false,
1013 VectorizedArrayType> &fe_eval)
1014 : fe_eval(fe_eval)
1015 {}
1016
1017
1018
1019 template <int dim,
1020 int fe_degree,
1021 int n_components,
1022 typename Number,
1023 typename VectorizedArrayType>
1024 inline void
1025 CellwiseInverseMassMatrix<dim,
1026 fe_degree,
1027 n_components,
1028 Number,
1029 VectorizedArrayType>::
1030 fill_inverse_JxW_values(
1031 AlignedVector<VectorizedArrayType> &inverse_jxw) const
1032 {
1033 constexpr unsigned int dofs_per_component_on_cell =
1034 Utilities::pow(fe_degree + 1, dim);
1035 Assert(inverse_jxw.size() > 0 &&
1036 inverse_jxw.size() % dofs_per_component_on_cell == 0,
1037 ExcMessage(
1038 "Expected diagonal to be a multiple of scalar dof per cells"));
1039
1040 // compute values for the first component
1041 for (unsigned int q = 0; q < dofs_per_component_on_cell; ++q)
1042 inverse_jxw[q] = 1. / fe_eval.JxW(q);
1043 // copy values to rest of vector
1044 for (unsigned int q = dofs_per_component_on_cell; q < inverse_jxw.size();)
1045 for (unsigned int i = 0; i < dofs_per_component_on_cell; ++i, ++q)
1046 inverse_jxw[q] = inverse_jxw[i];
1047 }
1048
1049
1050
1051 template <int dim,
1052 int fe_degree,
1053 int n_components,
1054 typename Number,
1055 typename VectorizedArrayType>
1056 inline void
1057 CellwiseInverseMassMatrix<
1058 dim,
1059 fe_degree,
1060 n_components,
1061 Number,
1062 VectorizedArrayType>::apply(const VectorizedArrayType *in_array,
1063 VectorizedArrayType * out_array) const
1064 {
1065 if (fe_degree > -1)
1066 internal::CellwiseInverseMassMatrixImplBasic<dim, VectorizedArrayType>::
1067 template run<fe_degree>(n_components, fe_eval, in_array, out_array);
1068 else
1069 internal::CellwiseInverseMassFactory<dim, VectorizedArrayType>::apply(
1070 n_components, fe_eval, in_array, out_array);
1071 }
1072
1073
1074
1075 template <int dim,
1076 int fe_degree,
1077 int n_components,
1078 typename Number,
1079 typename VectorizedArrayType>
1080 inline void
1081 CellwiseInverseMassMatrix<dim,
1082 fe_degree,
1083 n_components,
1084 Number,
1085 VectorizedArrayType>::
1086 apply(const AlignedVector<VectorizedArrayType> &inverse_coefficients,
1087 const unsigned int n_actual_components,
1088 const VectorizedArrayType * in_array,
1089 VectorizedArrayType * out_array) const
1090 {
1091 const unsigned int given_degree =
1092 fe_eval.get_shape_info().data[0].fe_degree;
1093 if (fe_degree > -1)
1094 internal::CellwiseInverseMassMatrixImplFlexible<dim,
1095 VectorizedArrayType>::
1096 template run<fe_degree>(
1097 n_actual_components,
1098 fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
1099 inverse_coefficients,
1100 in_array,
1101 out_array);
1102 else
1103 internal::CellwiseInverseMassFactory<dim, VectorizedArrayType>::apply(
1104 n_actual_components,
1105 given_degree,
1106 fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
1107 inverse_coefficients,
1108 in_array,
1109 out_array);
1110 }
1111
1112
1113
1114 template <int dim,
1115 int fe_degree,
1116 int n_components,
1117 typename Number,
1118 typename VectorizedArrayType>
1119 inline void
1120 CellwiseInverseMassMatrix<dim,
1121 fe_degree,
1122 n_components,
1123 Number,
1124 VectorizedArrayType>::
1125 transform_from_q_points_to_basis(const unsigned int n_actual_components,
1126 const VectorizedArrayType *in_array,
1127 VectorizedArrayType * out_array) const
1128 {
1129 const auto n_q_points_1d = fe_eval.get_shape_info().data[0].n_q_points_1d;
1130
1131 if (fe_degree > -1 && (fe_degree + 1 == n_q_points_1d))
1132 internal::CellwiseInverseMassMatrixImplTransformFromQPoints<
1133 dim,
1134 VectorizedArrayType>::template run<fe_degree,
1135 fe_degree + 1>(n_actual_components,
1136 fe_eval,
1137 in_array,
1138 out_array);
1139 else
1140 internal::CellwiseInverseMassFactory<dim, VectorizedArrayType>::
1141 transform_from_q_points_to_basis(n_actual_components,
1142 fe_eval,
1143 in_array,
1144 out_array);
1145 }
1146
1147
1148
1149 //----------------- Base operator -----------------------------
1150 template <int dim, typename VectorType, typename VectorizedArrayType>
1151 Base<dim, VectorType, VectorizedArrayType>::Base()
1152 : Subscriptor()
1153 , have_interface_matrices(false)
1154 {}
1155
1156
1157
1158 template <int dim, typename VectorType, typename VectorizedArrayType>
1159 typename Base<dim, VectorType, VectorizedArrayType>::size_type
1160 Base<dim, VectorType, VectorizedArrayType>::m() const
1161 {
1162 Assert(data.get() != nullptr, ExcNotInitialized());
1163 typename Base<dim, VectorType, VectorizedArrayType>::size_type total_size =
1164 0;
1165 for (const unsigned int selected_row : selected_rows)
1166 total_size += data->get_vector_partitioner(selected_row)->size();
1167 return total_size;
1168 }
1169
1170
1171
1172 template <int dim, typename VectorType, typename VectorizedArrayType>
1173 typename Base<dim, VectorType, VectorizedArrayType>::size_type
1174 Base<dim, VectorType, VectorizedArrayType>::n() const
1175 {
1176 Assert(data.get() != nullptr, ExcNotInitialized());
1177 typename Base<dim, VectorType, VectorizedArrayType>::size_type total_size =
1178 0;
1179 for (const unsigned int selected_column : selected_columns)
1180 total_size += data->get_vector_partitioner(selected_column)->size();
1181 return total_size;
1182 }
1183
1184
1185
1186 template <int dim, typename VectorType, typename VectorizedArrayType>
1187 void
1188 Base<dim, VectorType, VectorizedArrayType>::clear()
1189 {
1190 data.reset();
1191 inverse_diagonal_entries.reset();
1192 }
1193
1194
1195
1196 template <int dim, typename VectorType, typename VectorizedArrayType>
1197 typename Base<dim, VectorType, VectorizedArrayType>::value_type
1198 Base<dim, VectorType, VectorizedArrayType>::el(const unsigned int row,
1199 const unsigned int col) const
1200 {
1201 (void)col;
1202 Assert(row == col, ExcNotImplemented());
1203 Assert(inverse_diagonal_entries.get() != nullptr &&
1204 inverse_diagonal_entries->m() > 0,
1205 ExcNotInitialized());
1206 return 1.0 / (*inverse_diagonal_entries)(row, row);
1207 }
1208
1209
1210
1211 template <int dim, typename VectorType, typename VectorizedArrayType>
1212 void
1213 Base<dim, VectorType, VectorizedArrayType>::initialize_dof_vector(
1214 VectorType &vec) const
1215 {
1216 Assert(data.get() != nullptr, ExcNotInitialized());
1217 AssertDimension(BlockHelper::n_blocks(vec), selected_rows.size());
1218 for (unsigned int i = 0; i < BlockHelper::n_blocks(vec); ++i)
1219 {
1220 const unsigned int index = selected_rows[i];
1221 if (!BlockHelper::subblock(vec, index)
1222 .partitioners_are_compatible(
1223 *data->get_dof_info(index).vector_partitioner))
1224 data->initialize_dof_vector(BlockHelper::subblock(vec, index), index);
1225
1226 Assert(BlockHelper::subblock(vec, index)
1227 .partitioners_are_globally_compatible(
1228 *data->get_dof_info(index).vector_partitioner),
1229 ExcInternalError());
1230 }
1231 BlockHelper::collect_sizes(vec);
1232 }
1233
1234
1235
1236 template <int dim, typename VectorType, typename VectorizedArrayType>
1237 void
1238 Base<dim, VectorType, VectorizedArrayType>::initialize(
1239 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
1240 data_,
1241 const std::vector<unsigned int> &given_row_selection,
1242 const std::vector<unsigned int> &given_column_selection)
1243 {
1244 data = data_;
1245
1246 selected_rows.clear();
1247 selected_columns.clear();
1248 if (given_row_selection.empty())
1249 for (unsigned int i = 0; i < data_->n_components(); ++i)
1250 selected_rows.push_back(i);
1251 else
1252 {
1253 for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1254 {
1255 AssertIndexRange(given_row_selection[i], data_->n_components());
1256 for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1257 if (j != i)
1258 Assert(given_row_selection[j] != given_row_selection[i],
1259 ExcMessage("Given row indices must be unique"));
1260
1261 selected_rows.push_back(given_row_selection[i]);
1262 }
1263 }
1264 if (given_column_selection.size() == 0)
1265 selected_columns = selected_rows;
1266 else
1267 {
1268 for (unsigned int i = 0; i < given_column_selection.size(); ++i)
1269 {
1270 AssertIndexRange(given_column_selection[i], data_->n_components());
1271 for (unsigned int j = 0; j < given_column_selection.size(); ++j)
1272 if (j != i)
1273 Assert(given_column_selection[j] != given_column_selection[i],
1274 ExcMessage("Given column indices must be unique"));
1275
1276 selected_columns.push_back(given_column_selection[i]);
1277 }
1278 }
1279
1280 edge_constrained_indices.clear();
1281 edge_constrained_indices.resize(selected_rows.size());
1282 edge_constrained_values.clear();
1283 edge_constrained_values.resize(selected_rows.size());
1284 have_interface_matrices = false;
1285 }
1286
1287
1288
1289 template <int dim, typename VectorType, typename VectorizedArrayType>
1290 void
1291 Base<dim, VectorType, VectorizedArrayType>::initialize(
1292 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
1293 data_,
1294 const MGConstrainedDoFs & mg_constrained_dofs,
1295 const unsigned int level,
1296 const std::vector<unsigned int> &given_row_selection)
1297 {
1298 std::vector<MGConstrainedDoFs> mg_constrained_dofs_vector(
1299 1, mg_constrained_dofs);
1300 initialize(data_, mg_constrained_dofs_vector, level, given_row_selection);
1301 }
1302
1303
1304
1305 template <int dim, typename VectorType, typename VectorizedArrayType>
1306 void
1307 Base<dim, VectorType, VectorizedArrayType>::initialize(
1308 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
1309 data_,
1310 const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
1311 const unsigned int level,
1312 const std::vector<unsigned int> & given_row_selection)
1313 {
1314 AssertThrow(level != numbers::invalid_unsigned_int,
1315 ExcMessage("level is not set"));
1316
1317 selected_rows.clear();
1318 selected_columns.clear();
1319 if (given_row_selection.empty())
1320 for (unsigned int i = 0; i < data_->n_components(); ++i)
1321 selected_rows.push_back(i);
1322 else
1323 {
1324 for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1325 {
1326 AssertIndexRange(given_row_selection[i], data_->n_components());
1327 for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1328 if (j != i)
1329 Assert(given_row_selection[j] != given_row_selection[i],
1330 ExcMessage("Given row indices must be unique"));
1331
1332 selected_rows.push_back(given_row_selection[i]);
1333 }
1334 }
1335 selected_columns = selected_rows;
1336
1337 AssertDimension(mg_constrained_dofs.size(), selected_rows.size());
1338 edge_constrained_indices.clear();
1339 edge_constrained_indices.resize(selected_rows.size());
1340 edge_constrained_values.clear();
1341 edge_constrained_values.resize(selected_rows.size());
1342
1343 data = data_;
1344
1345 for (unsigned int j = 0; j < selected_rows.size(); ++j)
1346 {
1347 if (data_->n_cell_batches() > 0)
1348 {
1349 AssertDimension(level, data_->get_cell_iterator(0, 0, j)->level());
1350 }
1351
1352 // setup edge_constrained indices
1353 std::vector<types::global_dof_index> interface_indices;
1354 mg_constrained_dofs[j]
1355 .get_refinement_edge_indices(level)
1356 .fill_index_vector(interface_indices);
1357 edge_constrained_indices[j].clear();
1358 edge_constrained_indices[j].reserve(interface_indices.size());
1359 edge_constrained_values[j].resize(interface_indices.size());
1360 const IndexSet &locally_owned =
1361 data->get_dof_handler(selected_rows[j]).locally_owned_mg_dofs(level);
1362 for (const auto interface_index : interface_indices)
1363 if (locally_owned.is_element(interface_index))
1364 edge_constrained_indices[j].push_back(
1365 locally_owned.index_within_set(interface_index));
1366 have_interface_matrices |=
1367 Utilities::MPI::max(
1368 static_cast<unsigned int>(edge_constrained_indices[j].size()),
1369 data->get_vector_partitioner()->get_mpi_communicator()) > 0;
1370 }
1371 }
1372
1373
1374
1375 template <int dim, typename VectorType, typename VectorizedArrayType>
1376 void
1377 Base<dim, VectorType, VectorizedArrayType>::set_constrained_entries_to_one(
1378 VectorType &dst) const
1379 {
1380 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1381 {
1382 const std::vector<unsigned int> &constrained_dofs =
1383 data->get_constrained_dofs(selected_rows[j]);
1384 for (const auto constrained_dof : constrained_dofs)
1385 BlockHelper::subblock(dst, j).local_element(constrained_dof) = 1.;
1386 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1387 BlockHelper::subblock(dst, j).local_element(
1388 edge_constrained_indices[j][i]) = 1.;
1389 }
1390 }
1391
1392
1393
1394 template <int dim, typename VectorType, typename VectorizedArrayType>
1395 void
1396 Base<dim, VectorType, VectorizedArrayType>::vmult(VectorType & dst,
1397 const VectorType &src) const
1398 {
1399 using Number =
1400 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1401 dst = Number(0.);
1402 vmult_add(dst, src);
1403 }
1404
1405
1406
1407 template <int dim, typename VectorType, typename VectorizedArrayType>
1408 void
1409 Base<dim, VectorType, VectorizedArrayType>::vmult_add(
1410 VectorType & dst,
1411 const VectorType &src) const
1412 {
1413 mult_add(dst, src, false);
1414 }
1415
1416
1417
1418 template <int dim, typename VectorType, typename VectorizedArrayType>
1419 void
1420 Base<dim, VectorType, VectorizedArrayType>::Tvmult_add(
1421 VectorType & dst,
1422 const VectorType &src) const
1423 {
1424 mult_add(dst, src, true);
1425 }
1426
1427
1428
1429 template <int dim, typename VectorType, typename VectorizedArrayType>
1430 void
1431 Base<dim, VectorType, VectorizedArrayType>::adjust_ghost_range_if_necessary(
1432 const VectorType &src,
1433 const bool is_row) const
1434 {
1435 using Number =
1436 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1437 for (unsigned int i = 0; i < BlockHelper::n_blocks(src); ++i)
1438 {
1439 const unsigned int mf_component =
1440 is_row ? selected_rows[i] : selected_columns[i];
1441 // If both vectors use the same partitioner -> done
1442 if (BlockHelper::subblock(src, i).get_partitioner().get() ==
1443 data->get_dof_info(mf_component).vector_partitioner.get())
1444 continue;
1445
1446 // If not, assert that the local ranges are the same and reset to the
1447 // current partitioner
1448 Assert(BlockHelper::subblock(src, i)
1449 .get_partitioner()
1450 ->locally_owned_size() ==
1451 data->get_dof_info(mf_component)
1452 .vector_partitioner->locally_owned_size(),
1453 ExcMessage(
1454 "The vector passed to the vmult() function does not have "
1455 "the correct size for compatibility with MatrixFree."));
1456
1457 // copy the vector content to a temporary vector so that it does not get
1458 // lost
1459 LinearAlgebra::distributed::Vector<Number> copy_vec(
1460 BlockHelper::subblock(src, i));
1461 this->data->initialize_dof_vector(
1462 BlockHelper::subblock(const_cast<VectorType &>(src), i),
1463 mf_component);
1464 BlockHelper::subblock(const_cast<VectorType &>(src), i)
1465 .copy_locally_owned_data_from(copy_vec);
1466 }
1467 }
1468
1469
1470
1471 template <int dim, typename VectorType, typename VectorizedArrayType>
1472 void
1473 Base<dim, VectorType, VectorizedArrayType>::preprocess_constraints(
1474 VectorType & dst,
1475 const VectorType &src) const
1476 {
1477 using Number =
1478 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1479 adjust_ghost_range_if_necessary(src, false);
1480 adjust_ghost_range_if_necessary(dst, true);
1481
1482 // set zero Dirichlet values on the input vector (and remember the src and
1483 // dst values because we need to reset them at the end)
1484 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1485 {
1486 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1487 {
1488 edge_constrained_values[j][i] = std::pair<Number, Number>(
1489 BlockHelper::subblock(src, j).local_element(
1490 edge_constrained_indices[j][i]),
1491 BlockHelper::subblock(dst, j).local_element(
1492 edge_constrained_indices[j][i]));
1493 BlockHelper::subblock(const_cast<VectorType &>(src), j)
1494 .local_element(edge_constrained_indices[j][i]) = 0.;
1495 }
1496 }
1497 }
1498
1499
1500
1501 template <int dim, typename VectorType, typename VectorizedArrayType>
1502 void
1503 Base<dim, VectorType, VectorizedArrayType>::mult_add(
1504 VectorType & dst,
1505 const VectorType &src,
1506 const bool transpose) const
1507 {
1508 AssertDimension(dst.size(), src.size());
1509 AssertDimension(BlockHelper::n_blocks(dst), BlockHelper::n_blocks(src));
1510 AssertDimension(BlockHelper::n_blocks(dst), selected_rows.size());
1511 preprocess_constraints(dst, src);
1512 if (transpose)
1513 Tapply_add(dst, src);
1514 else
1515 apply_add(dst, src);
1516 postprocess_constraints(dst, src);
1517 }
1518
1519
1520
1521 template <int dim, typename VectorType, typename VectorizedArrayType>
1522 void
1523 Base<dim, VectorType, VectorizedArrayType>::postprocess_constraints(
1524 VectorType & dst,
1525 const VectorType &src) const
1526 {
1527 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1528 {
1529 const std::vector<unsigned int> &constrained_dofs =
1530 data->get_constrained_dofs(selected_rows[j]);
1531 for (const auto constrained_dof : constrained_dofs)
1532 BlockHelper::subblock(dst, j).local_element(constrained_dof) +=
1533 BlockHelper::subblock(src, j).local_element(constrained_dof);
1534 }
1535
1536 // reset edge constrained values, multiply by unit matrix and add into
1537 // destination
1538 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1539 {
1540 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1541 {
1542 BlockHelper::subblock(const_cast<VectorType &>(src), j)
1543 .local_element(edge_constrained_indices[j][i]) =
1544 edge_constrained_values[j][i].first;
1545 BlockHelper::subblock(dst, j).local_element(
1546 edge_constrained_indices[j][i]) =
1547 edge_constrained_values[j][i].second +
1548 edge_constrained_values[j][i].first;
1549 }
1550 }
1551 }
1552
1553
1554
1555 template <int dim, typename VectorType, typename VectorizedArrayType>
1556 void
1557 Base<dim, VectorType, VectorizedArrayType>::vmult_interface_down(
1558 VectorType & dst,
1559 const VectorType &src) const
1560 {
1561 using Number =
1562 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1563 AssertDimension(dst.size(), src.size());
1564 adjust_ghost_range_if_necessary(src, false);
1565 adjust_ghost_range_if_necessary(dst, true);
1566
1567 dst = Number(0.);
1568
1569 if (!have_interface_matrices)
1570 return;
1571
1572 // set zero Dirichlet values on the input vector (and remember the src and
1573 // dst values because we need to reset them at the end)
1574 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1575 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1576 {
1577 edge_constrained_values[j][i] = std::pair<Number, Number>(
1578 BlockHelper::subblock(src, j).local_element(
1579 edge_constrained_indices[j][i]),
1580 BlockHelper::subblock(dst, j).local_element(
1581 edge_constrained_indices[j][i]));
1582 BlockHelper::subblock(const_cast<VectorType &>(src), j)
1583 .local_element(edge_constrained_indices[j][i]) = 0.;
1584 }
1585
1586 apply_add(dst, src);
1587
1588 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1589 {
1590 unsigned int c = 0;
1591 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1592 {
1593 for (; c < edge_constrained_indices[j][i]; ++c)
1594 BlockHelper::subblock(dst, j).local_element(c) = 0.;
1595 ++c;
1596
1597 // reset the src values
1598 BlockHelper::subblock(const_cast<VectorType &>(src), j)
1599 .local_element(edge_constrained_indices[j][i]) =
1600 edge_constrained_values[j][i].first;
1601 }
1602 for (; c < BlockHelper::subblock(dst, j).locally_owned_size(); ++c)
1603 BlockHelper::subblock(dst, j).local_element(c) = 0.;
1604 }
1605 }
1606
1607
1608
1609 template <int dim, typename VectorType, typename VectorizedArrayType>
1610 void
1611 Base<dim, VectorType, VectorizedArrayType>::vmult_interface_up(
1612 VectorType & dst,
1613 const VectorType &src) const
1614 {
1615 using Number =
1616 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1617 AssertDimension(dst.size(), src.size());
1618 adjust_ghost_range_if_necessary(src, false);
1619 adjust_ghost_range_if_necessary(dst, true);
1620
1621 dst = Number(0.);
1622
1623 if (!have_interface_matrices)
1624 return;
1625
1626 VectorType src_cpy(src);
1627 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1628 {
1629 unsigned int c = 0;
1630 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1631 {
1632 for (; c < edge_constrained_indices[j][i]; ++c)
1633 BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1634 ++c;
1635 }
1636 for (; c < BlockHelper::subblock(src_cpy, j).locally_owned_size(); ++c)
1637 BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1638 }
1639
1640 apply_add(dst, src_cpy);
1641
1642 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1643 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1644 BlockHelper::subblock(dst, j).local_element(
1645 edge_constrained_indices[j][i]) = 0.;
1646 }
1647
1648
1649
1650 template <int dim, typename VectorType, typename VectorizedArrayType>
1651 void
1652 Base<dim, VectorType, VectorizedArrayType>::Tvmult(
1653 VectorType & dst,
1654 const VectorType &src) const
1655 {
1656 using Number =
1657 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1658 dst = Number(0.);
1659 Tvmult_add(dst, src);
1660 }
1661
1662
1663
1664 template <int dim, typename VectorType, typename VectorizedArrayType>
1665 std::size_t
1666 Base<dim, VectorType, VectorizedArrayType>::memory_consumption() const
1667 {
1668 return inverse_diagonal_entries.get() != nullptr ?
1669 inverse_diagonal_entries->memory_consumption() :
1670 sizeof(*this);
1671 }
1672
1673
1674
1675 template <int dim, typename VectorType, typename VectorizedArrayType>
1676 std::shared_ptr<const MatrixFree<
1677 dim,
1678 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
1679 VectorizedArrayType>>
1680 Base<dim, VectorType, VectorizedArrayType>::get_matrix_free() const
1681 {
1682 return data;
1683 }
1684
1685
1686
1687 template <int dim, typename VectorType, typename VectorizedArrayType>
1688 const std::shared_ptr<DiagonalMatrix<VectorType>> &
1689 Base<dim, VectorType, VectorizedArrayType>::get_matrix_diagonal_inverse()
1690 const
1691 {
1692 Assert(inverse_diagonal_entries.get() != nullptr &&
1693 inverse_diagonal_entries->m() > 0,
1694 ExcNotInitialized());
1695 return inverse_diagonal_entries;
1696 }
1697
1698
1699
1700 template <int dim, typename VectorType, typename VectorizedArrayType>
1701 const std::shared_ptr<DiagonalMatrix<VectorType>> &
1702 Base<dim, VectorType, VectorizedArrayType>::get_matrix_diagonal() const
1703 {
1704 Assert(diagonal_entries.get() != nullptr && diagonal_entries->m() > 0,
1705 ExcNotInitialized());
1706 return diagonal_entries;
1707 }
1708
1709
1710
1711 template <int dim, typename VectorType, typename VectorizedArrayType>
1712 void
1713 Base<dim, VectorType, VectorizedArrayType>::Tapply_add(
1714 VectorType & dst,
1715 const VectorType &src) const
1716 {
1717 apply_add(dst, src);
1718 }
1719
1720
1721
1722 template <int dim, typename VectorType, typename VectorizedArrayType>
1723 void
1724 Base<dim, VectorType, VectorizedArrayType>::precondition_Jacobi(
1725 VectorType & dst,
1726 const VectorType & src,
1727 const typename Base<dim, VectorType, VectorizedArrayType>::value_type omega)
1728 const
1729 {
1730 Assert(inverse_diagonal_entries.get() && inverse_diagonal_entries->m() > 0,
1731 ExcNotInitialized());
1732 inverse_diagonal_entries->vmult(dst, src);
1733 dst *= omega;
1734 }
1735
1736
1737
1738 //------------------------- MGInterfaceOperator ------------------------------
1739
1740 template <typename OperatorType>
1741 MGInterfaceOperator<OperatorType>::MGInterfaceOperator()
1742 : Subscriptor()
1743 , mf_base_operator(nullptr)
1744 {}
1745
1746
1747
1748 template <typename OperatorType>
1749 void
1750 MGInterfaceOperator<OperatorType>::clear()
1751 {
1752 mf_base_operator = nullptr;
1753 }
1754
1755
1756
1757 template <typename OperatorType>
1758 void
1759 MGInterfaceOperator<OperatorType>::initialize(const OperatorType &operator_in)
1760 {
1761 mf_base_operator = &operator_in;
1762 }
1763
1764
1765
1766 template <typename OperatorType>
1767 template <typename VectorType>
1768 void
1769 MGInterfaceOperator<OperatorType>::vmult(VectorType & dst,
1770 const VectorType &src) const
1771 {
1772#ifndef DEAL_II_MSVC
1773 static_assert(
1774 std::is_same<typename VectorType::value_type, value_type>::value,
1775 "The vector type must be based on the same value type as this "
1776 "operator");
1777#endif
1778
1779 Assert(mf_base_operator != nullptr, ExcNotInitialized());
1780
1781 mf_base_operator->vmult_interface_down(dst, src);
1782 }
1783
1784
1785
1786 template <typename OperatorType>
1787 template <typename VectorType>
1788 void
1789 MGInterfaceOperator<OperatorType>::Tvmult(VectorType & dst,
1790 const VectorType &src) const
1791 {
1792#ifndef DEAL_II_MSVC
1793 static_assert(
1794 std::is_same<typename VectorType::value_type, value_type>::value,
1795 "The vector type must be based on the same value type as this "
1796 "operator");
1797#endif
1798
1799 Assert(mf_base_operator != nullptr, ExcNotInitialized());
1800
1801 mf_base_operator->vmult_interface_up(dst, src);
1802 }
1803
1804
1805
1806 template <typename OperatorType>
1807 template <typename VectorType>
1808 void
1809 MGInterfaceOperator<OperatorType>::initialize_dof_vector(
1810 VectorType &vec) const
1811 {
1812 Assert(mf_base_operator != nullptr, ExcNotInitialized());
1813
1814 mf_base_operator->initialize_dof_vector(vec);
1815 }
1816
1817
1818
1819 //-----------------------------MassOperator----------------------------------
1820
1821 template <int dim,
1822 int fe_degree,
1823 int n_q_points_1d,
1824 int n_components,
1825 typename VectorType,
1826 typename VectorizedArrayType>
1827 MassOperator<dim,
1828 fe_degree,
1829 n_q_points_1d,
1830 n_components,
1831 VectorType,
1832 VectorizedArrayType>::MassOperator()
1833 : Base<dim, VectorType, VectorizedArrayType>()
1834 {}
1835
1836
1837
1838 template <int dim,
1839 int fe_degree,
1840 int n_q_points_1d,
1841 int n_components,
1842 typename VectorType,
1843 typename VectorizedArrayType>
1844 void
1845 MassOperator<dim,
1846 fe_degree,
1847 n_q_points_1d,
1848 n_components,
1849 VectorType,
1850 VectorizedArrayType>::compute_diagonal()
1851 {
1852 using Number =
1853 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1854 Assert((Base<dim, VectorType, VectorizedArrayType>::data.get() != nullptr),
1855 ExcNotInitialized());
1856 Assert(this->selected_rows == this->selected_columns,
1857 ExcMessage("This function is only implemented for square (not "
1858 "rectangular) operators."));
1859
1860 this->inverse_diagonal_entries =
1861 std::make_shared<DiagonalMatrix<VectorType>>();
1862 this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1863 VectorType &inverse_diagonal_vector =
1864 this->inverse_diagonal_entries->get_vector();
1865 VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1866 this->initialize_dof_vector(inverse_diagonal_vector);
1867 this->initialize_dof_vector(diagonal_vector);
1868
1869 // Set up the action of the mass matrix in a way that's compatible with
1870 // MatrixFreeTools::compute_diagonal:
1871 auto diagonal_evaluation = [](auto &integrator) {
1872 integrator.evaluate(EvaluationFlags::values);
1873 for (unsigned int q = 0; q < integrator.n_q_points; ++q)
1874 integrator.submit_value(integrator.get_value(q), q);
1875 integrator.integrate(EvaluationFlags::values);
1876 };
1877
1878 std::function<void(
1879 FEEvaluation<
1880 dim,
1881 fe_degree,
1882 n_q_points_1d,
1883 n_components,
1884 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
1885 VectorizedArrayType> &)>
1886 diagonal_evaluation_f(diagonal_evaluation);
1887
1888 Assert(this->selected_rows.size() > 0, ExcInternalError());
1889 for (unsigned int block_n = 0; block_n < this->selected_rows.size();
1890 ++block_n)
1891 MatrixFreeTools::compute_diagonal(*this->data,
1892 BlockHelper::subblock(diagonal_vector,
1893 block_n),
1894 diagonal_evaluation_f,
1895 this->selected_rows[block_n]);
1896
1897 // Constrained entries will create zeros on the main diagonal, which we
1898 // don't want
1899 this->set_constrained_entries_to_one(diagonal_vector);
1900
1901 inverse_diagonal_vector = diagonal_vector;
1902
1903 for (unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
1904 ++i)
1905 {
1906 Assert(diagonal_vector.local_element(i) > Number(0),
1907 ExcInternalError());
1908 inverse_diagonal_vector.local_element(i) =
1909 1. / inverse_diagonal_vector.local_element(i);
1910 }
1911
1912 // We never need ghost values so don't update them
1913 }
1914
1915
1916
1917 template <int dim,
1918 int fe_degree,
1919 int n_q_points_1d,
1920 int n_components,
1921 typename VectorType,
1922 typename VectorizedArrayType>
1923 void
1924 MassOperator<dim,
1925 fe_degree,
1926 n_q_points_1d,
1927 n_components,
1928 VectorType,
1929 VectorizedArrayType>::compute_lumped_diagonal()
1930 {
1931 using Number =
1932 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1933 Assert((Base<dim, VectorType, VectorizedArrayType>::data.get() != nullptr),
1934 ExcNotInitialized());
1935 Assert(this->selected_rows == this->selected_columns,
1936 ExcMessage("This function is only implemented for square (not "
1937 "rectangular) operators."));
1938
1939 inverse_lumped_diagonal_entries =
1940 std::make_shared<DiagonalMatrix<VectorType>>();
1941 lumped_diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1942 VectorType &inverse_lumped_diagonal_vector =
1943 inverse_lumped_diagonal_entries->get_vector();
1944 VectorType &lumped_diagonal_vector = lumped_diagonal_entries->get_vector();
1945 this->initialize_dof_vector(inverse_lumped_diagonal_vector);
1946 this->initialize_dof_vector(lumped_diagonal_vector);
1947
1948 // Re-use the inverse_lumped_diagonal_vector as the vector of 1s
1949 inverse_lumped_diagonal_vector = Number(1.);
1950 apply_add(lumped_diagonal_vector, inverse_lumped_diagonal_vector);
1951 this->set_constrained_entries_to_one(lumped_diagonal_vector);
1952
1953 const size_type locally_owned_size =
1954 inverse_lumped_diagonal_vector.locally_owned_size();
1955 // A caller may request a lumped diagonal matrix when it doesn't make sense
1956 // (e.g., an element with negative-mean basis functions). Avoid division by
1957 // zero so we don't cause a floating point exception but permit negative
1958 // entries here.
1959 for (size_type i = 0; i < locally_owned_size; ++i)
1960 {
1961 if (lumped_diagonal_vector.local_element(i) == Number(0.))
1962 inverse_lumped_diagonal_vector.local_element(i) = Number(1.);
1963 else
1964 inverse_lumped_diagonal_vector.local_element(i) =
1965 Number(1.) / lumped_diagonal_vector.local_element(i);
1966 }
1967
1968 inverse_lumped_diagonal_vector.update_ghost_values();
1969 lumped_diagonal_vector.update_ghost_values();
1970 }
1971
1972
1973
1974 template <int dim,
1975 int fe_degree,
1976 int n_q_points_1d,
1977 int n_components,
1978 typename VectorType,
1979 typename VectorizedArrayType>
1980 const std::shared_ptr<DiagonalMatrix<VectorType>> &
1981 MassOperator<dim,
1982 fe_degree,
1983 n_q_points_1d,
1984 n_components,
1985 VectorType,
1986 VectorizedArrayType>::get_matrix_lumped_diagonal_inverse() const
1987 {
1988 Assert(inverse_lumped_diagonal_entries.get() != nullptr &&
1989 inverse_lumped_diagonal_entries->m() > 0,
1990 ExcNotInitialized());
1991 return inverse_lumped_diagonal_entries;
1992 }
1993
1994
1995
1996 template <int dim,
1997 int fe_degree,
1998 int n_q_points_1d,
1999 int n_components,
2000 typename VectorType,
2001 typename VectorizedArrayType>
2002 const std::shared_ptr<DiagonalMatrix<VectorType>> &
2003 MassOperator<dim,
2004 fe_degree,
2005 n_q_points_1d,
2006 n_components,
2007 VectorType,
2008 VectorizedArrayType>::get_matrix_lumped_diagonal() const
2009 {
2010 Assert(lumped_diagonal_entries.get() != nullptr &&
2011 lumped_diagonal_entries->m() > 0,
2012 ExcNotInitialized());
2013 return lumped_diagonal_entries;
2014 }
2015
2016
2017
2018 template <int dim,
2019 int fe_degree,
2020 int n_q_points_1d,
2021 int n_components,
2022 typename VectorType,
2023 typename VectorizedArrayType>
2024 void
2025 MassOperator<dim,
2026 fe_degree,
2027 n_q_points_1d,
2028 n_components,
2029 VectorType,
2030 VectorizedArrayType>::apply_add(VectorType & dst,
2031 const VectorType &src) const
2032 {
2033 Base<dim, VectorType, VectorizedArrayType>::data->cell_loop(
2034 &MassOperator::local_apply_cell, this, dst, src);
2035 }
2036
2037
2038
2039 template <int dim,
2040 int fe_degree,
2041 int n_q_points_1d,
2042 int n_components,
2043 typename VectorType,
2044 typename VectorizedArrayType>
2045 void
2046 MassOperator<dim,
2047 fe_degree,
2048 n_q_points_1d,
2049 n_components,
2050 VectorType,
2051 VectorizedArrayType>::
2052 local_apply_cell(
2053 const MatrixFree<
2054 dim,
2055 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2056 VectorizedArrayType> & data,
2057 VectorType & dst,
2058 const VectorType & src,
2059 const std::pair<unsigned int, unsigned int> &cell_range) const
2060 {
2061 using Number =
2062 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2063 FEEvaluation<dim,
2064 fe_degree,
2065 n_q_points_1d,
2066 n_components,
2067 Number,
2068 VectorizedArrayType>
2069 phi(data, this->selected_rows[0]);
2070 for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2071 {
2072 phi.reinit(cell);
2073 phi.read_dof_values(src);
2074 phi.evaluate(EvaluationFlags::values);
2075 for (unsigned int q = 0; q < phi.n_q_points; ++q)
2076 phi.submit_value(phi.get_value(q), q);
2077 phi.integrate(EvaluationFlags::values);
2078 phi.distribute_local_to_global(dst);
2079 }
2080 }
2081
2082
2083 //-----------------------------LaplaceOperator----------------------------------
2084
2085 template <int dim,
2086 int fe_degree,
2087 int n_q_points_1d,
2088 int n_components,
2089 typename VectorType,
2090 typename VectorizedArrayType>
2091 LaplaceOperator<dim,
2092 fe_degree,
2093 n_q_points_1d,
2094 n_components,
2095 VectorType,
2096 VectorizedArrayType>::LaplaceOperator()
2097 : Base<dim, VectorType, VectorizedArrayType>()
2098 {}
2099
2100
2101
2102 template <int dim,
2103 int fe_degree,
2104 int n_q_points_1d,
2105 int n_components,
2106 typename VectorType,
2107 typename VectorizedArrayType>
2108 void
2109 LaplaceOperator<dim,
2110 fe_degree,
2111 n_q_points_1d,
2112 n_components,
2113 VectorType,
2114 VectorizedArrayType>::clear()
2115 {
2116 Base<dim, VectorType, VectorizedArrayType>::clear();
2117 scalar_coefficient.reset();
2118 }
2119
2120
2121
2122 template <int dim,
2123 int fe_degree,
2124 int n_q_points_1d,
2125 int n_components,
2126 typename VectorType,
2127 typename VectorizedArrayType>
2128 void
2129 LaplaceOperator<dim,
2130 fe_degree,
2131 n_q_points_1d,
2132 n_components,
2133 VectorType,
2134 VectorizedArrayType>::
2135 set_coefficient(
2136 const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient_)
2137 {
2138 scalar_coefficient = scalar_coefficient_;
2139 }
2140
2141
2142
2143 template <int dim,
2144 int fe_degree,
2145 int n_q_points_1d,
2146 int n_components,
2147 typename VectorType,
2148 typename VectorizedArrayType>
2149 std::shared_ptr<Table<2, VectorizedArrayType>>
2150 LaplaceOperator<dim,
2151 fe_degree,
2152 n_q_points_1d,
2153 n_components,
2154 VectorType,
2155 VectorizedArrayType>::get_coefficient()
2156 {
2157 Assert(scalar_coefficient.get(), ExcNotInitialized());
2158 return scalar_coefficient;
2159 }
2160
2161
2162
2163 template <int dim,
2164 int fe_degree,
2165 int n_q_points_1d,
2166 int n_components,
2167 typename VectorType,
2168 typename VectorizedArrayType>
2169 void
2170 LaplaceOperator<dim,
2171 fe_degree,
2172 n_q_points_1d,
2173 n_components,
2174 VectorType,
2175 VectorizedArrayType>::compute_diagonal()
2176 {
2177 using Number =
2178 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2179 Assert((Base<dim, VectorType, VectorizedArrayType>::data.get() != nullptr),
2180 ExcNotInitialized());
2181
2182 this->inverse_diagonal_entries =
2183 std::make_shared<DiagonalMatrix<VectorType>>();
2184 this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
2185 VectorType &inverse_diagonal_vector =
2186 this->inverse_diagonal_entries->get_vector();
2187 VectorType &diagonal_vector = this->diagonal_entries->get_vector();
2188 this->initialize_dof_vector(inverse_diagonal_vector);
2189 this->initialize_dof_vector(diagonal_vector);
2190
2191 this->data->cell_loop(&LaplaceOperator::local_diagonal_cell,
2192 this,
2193 diagonal_vector,
2194 /*unused*/ diagonal_vector);
2195 this->set_constrained_entries_to_one(diagonal_vector);
2196
2197 inverse_diagonal_vector = diagonal_vector;
2198
2199 for (unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
2200 ++i)
2201 if (std::abs(inverse_diagonal_vector.local_element(i)) >
2202 std::sqrt(std::numeric_limits<Number>::epsilon()))
2203 inverse_diagonal_vector.local_element(i) =
2204 1. / inverse_diagonal_vector.local_element(i);
2205 else
2206 inverse_diagonal_vector.local_element(i) = 1.;
2207
2208 // We never need ghost values so don't update them
2209 }
2210
2211
2212
2213 template <int dim,
2214 int fe_degree,
2215 int n_q_points_1d,
2216 int n_components,
2217 typename VectorType,
2218 typename VectorizedArrayType>
2219 void
2220 LaplaceOperator<dim,
2221 fe_degree,
2222 n_q_points_1d,
2223 n_components,
2224 VectorType,
2225 VectorizedArrayType>::apply_add(VectorType & dst,
2226 const VectorType &src) const
2227 {
2228 Base<dim, VectorType, VectorizedArrayType>::data->cell_loop(
2229 &LaplaceOperator::local_apply_cell, this, dst, src);
2230 }
2231
2232 namespace Implementation
2233 {
2234 template <typename VectorizedArrayType>
2235 bool
2236 non_negative(const VectorizedArrayType &n)
2237 {
2238 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2239 if (n[v] < 0.)
2240 return false;
2241
2242 return true;
2243 }
2244 } // namespace Implementation
2245
2246
2247
2248 template <int dim,
2249 int fe_degree,
2250 int n_q_points_1d,
2251 int n_components,
2252 typename VectorType,
2253 typename VectorizedArrayType>
2254 void
2255 LaplaceOperator<dim,
2256 fe_degree,
2257 n_q_points_1d,
2258 n_components,
2259 VectorType,
2260 VectorizedArrayType>::
2261 do_operation_on_cell(
2262 FEEvaluation<
2263 dim,
2264 fe_degree,
2265 n_q_points_1d,
2266 n_components,
2267 typename Base<dim, VectorType, VectorizedArrayType>::value_type> &phi,
2268 const unsigned int cell) const
2269 {
2270 phi.evaluate(EvaluationFlags::gradients);
2271 if (scalar_coefficient.get())
2272 {
2273 Assert(scalar_coefficient->size(1) == 1 ||
2274 scalar_coefficient->size(1) == phi.n_q_points,
2275 ExcMessage("The number of columns in the coefficient table must "
2276 "be either 1 or the number of quadrature points " +
2277 std::to_string(phi.n_q_points) +
2278 ", but the given value was " +
2279 std::to_string(scalar_coefficient->size(1))));
2280 if (scalar_coefficient->size(1) == phi.n_q_points)
2281 for (unsigned int q = 0; q < phi.n_q_points; ++q)
2282 {
2283 Assert(Implementation::non_negative(
2284 (*scalar_coefficient)(cell, q)),
2285 ExcMessage("Coefficient must be non-negative"));
2286 phi.submit_gradient((*scalar_coefficient)(cell, q) *
2287 phi.get_gradient(q),
2288 q);
2289 }
2290 else
2291 {
2292 Assert(Implementation::non_negative((*scalar_coefficient)(cell, 0)),
2293 ExcMessage("Coefficient must be non-negative"));
2294 const VectorizedArrayType coefficient =
2295 (*scalar_coefficient)(cell, 0);
2296 for (unsigned int q = 0; q < phi.n_q_points; ++q)
2297 phi.submit_gradient(coefficient * phi.get_gradient(q), q);
2298 }
2299 }
2300 else
2301 {
2302 for (unsigned int q = 0; q < phi.n_q_points; ++q)
2303 {
2304 phi.submit_gradient(phi.get_gradient(q), q);
2305 }
2306 }
2307 phi.integrate(EvaluationFlags::gradients);
2308 }
2309
2310
2311
2312 template <int dim,
2313 int fe_degree,
2314 int n_q_points_1d,
2315 int n_components,
2316 typename VectorType,
2317 typename VectorizedArrayType>
2318 void
2319 LaplaceOperator<dim,
2320 fe_degree,
2321 n_q_points_1d,
2322 n_components,
2323 VectorType,
2324 VectorizedArrayType>::
2325 local_apply_cell(
2326 const MatrixFree<
2327 dim,
2328 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2329 VectorizedArrayType> & data,
2330 VectorType & dst,
2331 const VectorType & src,
2332 const std::pair<unsigned int, unsigned int> &cell_range) const
2333 {
2334 using Number =
2335 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2336 FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number> phi(
2337 data, this->selected_rows[0]);
2338 for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2339 {
2340 phi.reinit(cell);
2341 phi.read_dof_values(src);
2342 do_operation_on_cell(phi, cell);
2343 phi.distribute_local_to_global(dst);
2344 }
2345 }
2346
2347
2348 template <int dim,
2349 int fe_degree,
2350 int n_q_points_1d,
2351 int n_components,
2352 typename VectorType,
2353 typename VectorizedArrayType>
2354 void
2355 LaplaceOperator<dim,
2356 fe_degree,
2357 n_q_points_1d,
2358 n_components,
2359 VectorType,
2360 VectorizedArrayType>::
2361 local_diagonal_cell(
2362 const MatrixFree<
2363 dim,
2364 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2365 VectorizedArrayType> &data,
2366 VectorType & dst,
2367 const VectorType &,
2368 const std::pair<unsigned int, unsigned int> &cell_range) const
2369 {
2370 using Number =
2371 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2372
2373 FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number> phi(
2374 data, this->selected_rows[0]);
2375 for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2376 {
2377 phi.reinit(cell);
2378 VectorizedArrayType local_diagonal_vector[phi.static_dofs_per_cell];
2379 for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2380 {
2381 for (unsigned int j = 0; j < phi.dofs_per_component; ++j)
2382 phi.begin_dof_values()[j] = VectorizedArrayType();
2383 phi.begin_dof_values()[i] = 1.;
2384 do_operation_on_cell(phi, cell);
2385 local_diagonal_vector[i] = phi.begin_dof_values()[i];
2386 }
2387 for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2388 for (unsigned int c = 0; c < phi.n_components; ++c)
2389 phi.begin_dof_values()[i + c * phi.dofs_per_component] =
2390 local_diagonal_vector[i];
2391 phi.distribute_local_to_global(dst);
2392 }
2393 }
2394
2395
2396} // end of namespace MatrixFreeOperators
2397
2398
2399DEAL_II_NAMESPACE_CLOSE
2400
2401#endif
AlignedVector::size
size_type size() const
FEEvaluationBase::distribute_local_to_global
void distribute_local_to_global(VectorType &dst, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip()) const
FEEvaluationBase::read_dof_values
void read_dof_values(const VectorType &src, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip())
FEEvaluationData::begin_dof_values
const Number * begin_dof_values() const
FEEvaluation::dofs_per_component
const unsigned int dofs_per_component
Definition: fe_evaluation.h:2334
FEEvaluation::reinit
void reinit(const unsigned int cell_batch_index)
FEEvaluation::static_dofs_per_cell
static constexpr unsigned int static_dofs_per_cell
Definition: fe_evaluation.h:1986
FEEvaluation::n_components
static constexpr unsigned int n_components
Definition: fe_evaluation.h:1948
IndexSet::index_within_set
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1923
IndexSet::is_element
bool is_element(const size_type index) const
Definition: index_set.h:1767
MatrixFreeOperators::Base::~Base
virtual ~Base() override=default
MatrixFreeOperators::Base::selected_rows
std::vector< unsigned int > selected_rows
Definition: operators.h:456
MatrixFreeOperators::Base::Tvmult_add
void Tvmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1420
MatrixFreeOperators::Base::set_constrained_entries_to_one
void set_constrained_entries_to_one(VectorType &dst) const
Definition: operators.h:1377
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:1409
MatrixFreeOperators::Base::vmult_interface_down
void vmult_interface_down(VectorType &dst, const VectorType &src) const
Definition: operators.h:1557
MatrixFreeOperators::Base::n
size_type n() const
Definition: operators.h:1174
MatrixFreeOperators::Base::preprocess_constraints
void preprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1473
MatrixFreeOperators::Base::mult_add
void mult_add(VectorType &dst, const VectorType &src, const bool transpose) const
Definition: operators.h:1503
MatrixFreeOperators::Base::edge_constrained_indices
std::vector< std::vector< unsigned int > > edge_constrained_indices
Definition: operators.h:468
MatrixFreeOperators::Base::get_matrix_diagonal
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal() const
Definition: operators.h:1702
MatrixFreeOperators::Base::Tvmult
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1652
MatrixFreeOperators::Base::edge_constrained_values
std::vector< std::vector< std::pair< value_type, value_type > > > edge_constrained_values
Definition: operators.h:474
MatrixFreeOperators::Base::m
size_type m() const
Definition: operators.h:1160
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:1238
MatrixFreeOperators::Base::selected_columns
std::vector< unsigned int > selected_columns
Definition: operators.h:462
MatrixFreeOperators::Base::get_matrix_free
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > get_matrix_free() const
Definition: operators.h:1680
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:1307
MatrixFreeOperators::Base::get_matrix_diagonal_inverse
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal_inverse() const
Definition: operators.h:1689
MatrixFreeOperators::Base::clear
virtual void clear()
Definition: operators.h:1188
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:1291
MatrixFreeOperators::Base::initialize_dof_vector
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1213
MatrixFreeOperators::Base::diagonal_entries
std::shared_ptr< DiagonalMatrix< VectorType > > diagonal_entries
Definition: operators.h:444
MatrixFreeOperators::Base::Tapply_add
virtual void Tapply_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1713
MatrixFreeOperators::Base::data
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > data
Definition: operators.h:438
MatrixFreeOperators::Base::adjust_ghost_range_if_necessary
void adjust_ghost_range_if_necessary(const VectorType &vec, const bool is_row) const
Definition: operators.h:1431
MatrixFreeOperators::Base::memory_consumption
virtual std::size_t memory_consumption() const
Definition: operators.h:1666
MatrixFreeOperators::Base::inverse_diagonal_entries
std::shared_ptr< DiagonalMatrix< VectorType > > inverse_diagonal_entries
Definition: operators.h:450
MatrixFreeOperators::Base::el
value_type el(const unsigned int row, const unsigned int col) const
Definition: operators.h:1198
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:201
MatrixFreeOperators::Base::precondition_Jacobi
void precondition_Jacobi(VectorType &dst, const VectorType &src, const value_type omega) const
Definition: operators.h:1724
MatrixFreeOperators::Base::value_type
typename VectorType::value_type value_type
Definition: operators.h:196
MatrixFreeOperators::Base::vmult
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1396
MatrixFreeOperators::Base::vmult_interface_up
void vmult_interface_up(VectorType &dst, const VectorType &src) const
Definition: operators.h:1611
MatrixFreeOperators::Base::postprocess_constraints
void postprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1523
MatrixFreeOperators::CellwiseInverseMassMatrix::fe_eval
const FEEvaluationBase< dim, n_components, Number, false, VectorizedArrayType > & fe_eval
Definition: operators.h:725
MatrixFreeOperators::CellwiseInverseMassMatrix::fill_inverse_JxW_values
void fill_inverse_JxW_values(AlignedVector< VectorizedArrayType > &inverse_jxw) const
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:1125
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:1086
MatrixFreeOperators::LaplaceOperator::get_coefficient
std::shared_ptr< Table< 2, VectorizedArrayType > > get_coefficient()
Definition: operators.h:2155
MatrixFreeOperators::LaplaceOperator::value_type
typename Base< dim, VectorType, VectorizedArrayType >::value_type value_type
Definition: operators.h:858
MatrixFreeOperators::LaplaceOperator::compute_diagonal
virtual void compute_diagonal() override
Definition: operators.h:2175
MatrixFreeOperators::LaplaceOperator::scalar_coefficient
std::shared_ptr< Table< 2, VectorizedArrayType > > scalar_coefficient
Definition: operators.h:991
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:2361
MatrixFreeOperators::LaplaceOperator::do_operation_on_cell
void do_operation_on_cell(FEEvaluation< dim, fe_degree, n_q_points_1d, n_components, value_type > &phi, const unsigned int cell) const
Definition: operators.h:2261
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:2325
MatrixFreeOperators::LaplaceOperator::apply_add
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition: operators.h:2225
MatrixFreeOperators::LaplaceOperator::size_type
typename Base< dim, VectorType, VectorizedArrayType >::size_type size_type
Definition: operators.h:864
MatrixFreeOperators::LaplaceOperator::set_coefficient
void set_coefficient(const std::shared_ptr< Table< 2, VectorizedArrayType > > &scalar_coefficient)
Definition: operators.h:2135
MatrixFreeOperators::LaplaceOperator::clear
virtual void clear() override
Definition: operators.h:2114
MatrixFreeOperators::MGInterfaceOperator::value_type
typename OperatorType::value_type value_type
Definition: operators.h:546
MatrixFreeOperators::MGInterfaceOperator::mf_base_operator
SmartPointer< const OperatorType > mf_base_operator
Definition: operators.h:596
MatrixFreeOperators::MGInterfaceOperator::vmult
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1769
MatrixFreeOperators::MGInterfaceOperator::initialize
void initialize(const OperatorType &operator_in)
Definition: operators.h:1759
MatrixFreeOperators::MGInterfaceOperator::Tvmult
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1789
MatrixFreeOperators::MGInterfaceOperator::initialize_dof_vector
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1809
MatrixFreeOperators::MGInterfaceOperator::size_type
typename OperatorType::size_type size_type
Definition: operators.h:551
MatrixFreeOperators::MassOperator::lumped_diagonal_entries
std::shared_ptr< DiagonalMatrix< VectorType > > lumped_diagonal_entries
Definition: operators.h:823
MatrixFreeOperators::MassOperator::size_type
typename Base< dim, VectorType, VectorizedArrayType >::size_type size_type
Definition: operators.h:757
MatrixFreeOperators::MassOperator::apply_add
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition: operators.h:2030
MatrixFreeOperators::MassOperator::get_matrix_lumped_diagonal
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_lumped_diagonal() const
Definition: operators.h:2008
MatrixFreeOperators::MassOperator::get_matrix_lumped_diagonal_inverse
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_lumped_diagonal_inverse() const
Definition: operators.h:1986
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:2052
MatrixFreeOperators::MassOperator::compute_diagonal
virtual void compute_diagonal() override
Definition: operators.h:1850
MatrixFreeOperators::MassOperator::value_type
typename Base< dim, VectorType, VectorizedArrayType >::value_type value_type
Definition: operators.h:751
MatrixFreeOperators::MassOperator::inverse_lumped_diagonal_entries
std::shared_ptr< DiagonalMatrix< VectorType > > inverse_lumped_diagonal_entries
Definition: operators.h:829
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition: derivative_form.h:538
diagonal_matrix.h
fe_evaluation.h
level
unsigned int level
Definition: grid_out.cc:4606
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1473
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
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:1583
la_parallel_vector.h
matrix_free.h
mg_constrained_dofs.h
MatrixFreeOperators::BlockHelper::collect_sizes
std::enable_if< IsBlockVector< VectorType >::value, void >::type collect_sizes(VectorType &vector)
Definition: operators.h:98
MatrixFreeOperators::BlockHelper::n_blocks
std::enable_if< IsBlockVector< VectorType >::value, unsignedint >::type n_blocks(const VectorType &vector)
Definition: operators.h:50
MatrixFreeOperators::BlockHelper::subblock
std::enable_if< IsBlockVector< VectorType >::value, typenameVectorType::BlockType & >::type subblock(VectorType &vector, unsigned int block_no)
Definition: operators.h:66
MatrixFreeOperators::Implementation::non_negative
bool non_negative(const VectorizedArrayType &n)
Definition: operators.h:2236
MatrixFreeTools::compute_diagonal
void compute_diagonal(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, VectorType &diagonal_global, const std::function< void(FEEvaluation< dim, fe_degree, n_q_points_1d, n_components, Number, VectorizedArrayType > &)> &local_vmult, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
Utilities::MPI::max
T max(const T &t, const MPI_Comm &mpi_communicator)
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:462
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:201
std::sqrt
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5510
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5572
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.5