Reference documentation for deal.II version 9.3.3
\(\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\}}\)
operators.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2011 - 2021 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
33#include <deal.II/multigrid/mg_constrained_dofs.h>
34
35
36DEAL_II_NAMESPACE_OPEN
37
38
39namespace MatrixFreeOperators
40{
41 namespace BlockHelper
42 {
43 // workaroud for unifying non-block vector and block vector implementations
44 // a non-block vector has one block and the only subblock is the vector
45 // itself
46 template <typename VectorType>
47 typename std::enable_if<IsBlockVector<VectorType>::value,
48 unsigned int>::type
49 n_blocks(const VectorType &vector)
50 {
51 return vector.n_blocks();
52 }
53
54 template <typename VectorType>
55 typename std::enable_if<!IsBlockVector<VectorType>::value,
56 unsigned int>::type
57 n_blocks(const VectorType &)
58 {
59 return 1;
60 }
61
62 template <typename VectorType>
63 typename std::enable_if<IsBlockVector<VectorType>::value,
64 typename VectorType::BlockType &>::type
65 subblock(VectorType &vector, unsigned int block_no)
66 {
67 return vector.block(block_no);
68 }
69
70 template <typename VectorType>
71 typename std::enable_if<IsBlockVector<VectorType>::value,
72 const typename VectorType::BlockType &>::type
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 typename std::enable_if<!IsBlockVector<VectorType>::value,
81 VectorType &>::type
82 subblock(VectorType &vector, unsigned int)
83 {
84 return vector;
85 }
86
87 template <typename VectorType>
88 typename std::enable_if<!IsBlockVector<VectorType>::value,
89 const VectorType &>::type
90 subblock(const VectorType &vector, unsigned int)
91 {
92 return vector;
93 }
94
95 template <typename VectorType>
96 typename std::enable_if<IsBlockVector<VectorType>::value, void>::type
97 collect_sizes(VectorType &vector)
98 {
99 vector.collect_sizes();
100 }
101
102 template <typename VectorType>
103 typename std::enable_if<!IsBlockVector<VectorType>::value, void>::type
104 collect_sizes(const VectorType &)
105 {}
106 } // namespace BlockHelper
107
185 template <int dim,
186 typename VectorType = LinearAlgebra::distributed::Vector<double>,
187 typename VectorizedArrayType =
188 VectorizedArray<typename VectorType::value_type>>
189 class Base : public Subscriptor
190 {
191 public:
195 using value_type = typename VectorType::value_type;
196
200 using size_type = typename VectorType::size_type;
201
205 Base();
206
210 virtual ~Base() override = default;
211
216 virtual void
217 clear();
218
235 void
236 initialize(std::shared_ptr<
237 const MatrixFree<dim, value_type, VectorizedArrayType>> data,
238 const std::vector<unsigned int> &selected_row_blocks =
239 std::vector<unsigned int>(),
240 const std::vector<unsigned int> &selected_column_blocks =
241 std::vector<unsigned int>());
242
256 void
257 initialize(std::shared_ptr<
258 const MatrixFree<dim, value_type, VectorizedArrayType>> data,
259 const MGConstrainedDoFs & mg_constrained_dofs,
260 const unsigned int level,
261 const std::vector<unsigned int> &selected_row_blocks =
262 std::vector<unsigned int>());
263
278 void
279 initialize(std::shared_ptr<
280 const MatrixFree<dim, value_type, VectorizedArrayType>> data_,
281 const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
282 const unsigned int level,
283 const std::vector<unsigned int> & selected_row_blocks =
284 std::vector<unsigned int>());
285
289 size_type
290 m() const;
291
295 size_type
296 n() const;
297
301 void
302 vmult_interface_down(VectorType &dst, const VectorType &src) const;
303
307 void
308 vmult_interface_up(VectorType &dst, const VectorType &src) const;
309
313 void
314 vmult(VectorType &dst, const VectorType &src) const;
315
319 void
320 Tvmult(VectorType &dst, const VectorType &src) const;
321
325 void
326 vmult_add(VectorType &dst, const VectorType &src) const;
327
331 void
332 Tvmult_add(VectorType &dst, const VectorType &src) const;
333
338 value_type
339 el(const unsigned int row, const unsigned int col) const;
340
345 virtual std::size_t
346 memory_consumption() const;
347
351 void
352 initialize_dof_vector(VectorType &vec) const;
353
361 virtual void
362 compute_diagonal() = 0;
363
367 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
368 get_matrix_free() const;
369
373 const std::shared_ptr<DiagonalMatrix<VectorType>> &
374 get_matrix_diagonal_inverse() const;
375
379 const std::shared_ptr<DiagonalMatrix<VectorType>> &
380 get_matrix_diagonal() const;
381
382
388 void
389 precondition_Jacobi(VectorType & dst,
390 const VectorType &src,
391 const value_type omega) const;
392
393 protected:
398 void
399 preprocess_constraints(VectorType &dst, const VectorType &src) const;
400
405 void
406 postprocess_constraints(VectorType &dst, const VectorType &src) const;
407
412 void
413 set_constrained_entries_to_one(VectorType &dst) const;
414
418 virtual void
419 apply_add(VectorType &dst, const VectorType &src) const = 0;
420
426 virtual void
427 Tapply_add(VectorType &dst, const VectorType &src) const;
428
432 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
433 data;
434
439 std::shared_ptr<DiagonalMatrix<VectorType>> diagonal_entries;
440
445 std::shared_ptr<DiagonalMatrix<VectorType>> inverse_diagonal_entries;
446
451 std::vector<unsigned int> selected_rows;
452
457 std::vector<unsigned int> selected_columns;
458
459 private:
463 std::vector<std::vector<unsigned int>> edge_constrained_indices;
464
468 mutable std::vector<std::vector<std::pair<value_type, value_type>>>
469 edge_constrained_values;
470
475 bool have_interface_matrices;
476
481 void
482 mult_add(VectorType & dst,
483 const VectorType &src,
484 const bool transpose) const;
485
493 void
494 adjust_ghost_range_if_necessary(const VectorType &vec,
495 const bool is_row) const;
496 };
497
498
499
534 template <typename OperatorType>
535 class MGInterfaceOperator : public Subscriptor
536 {
537 public:
541 using value_type = typename OperatorType::value_type;
542
546 using size_type = typename OperatorType::size_type;
547
551 MGInterfaceOperator();
552
556 void
557 clear();
558
562 void
563 initialize(const OperatorType &operator_in);
564
568 template <typename VectorType>
569 void
570 vmult(VectorType &dst, const VectorType &src) const;
571
575 template <typename VectorType>
576 void
577 Tvmult(VectorType &dst, const VectorType &src) const;
578
582 template <typename VectorType>
583 void
584 initialize_dof_vector(VectorType &vec) const;
585
586
587 private:
591 SmartPointer<const OperatorType> mf_base_operator;
592 };
593
594
595
613 template <int dim,
614 int fe_degree,
615 int n_components = 1,
616 typename Number = double,
617 typename VectorizedArrayType = VectorizedArray<Number>>
618 class CellwiseInverseMassMatrix
619 {
620 static_assert(
621 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
622 "Type of Number and of VectorizedArrayType do not match.");
623
624 public:
629 CellwiseInverseMassMatrix(
630 const FEEvaluationBase<dim,
631 n_components,
632 Number,
633 false,
634 VectorizedArrayType> &fe_eval);
635
644 void
645 apply(const AlignedVector<VectorizedArrayType> &inverse_coefficient,
646 const unsigned int n_actual_components,
647 const VectorizedArrayType * in_array,
648 VectorizedArrayType * out_array) const;
649
661 void
662 apply(const VectorizedArrayType *in_array,
663 VectorizedArrayType * out_array) const;
664
698 void
699 transform_from_q_points_to_basis(const unsigned int n_actual_components,
700 const VectorizedArrayType *in_array,
701 VectorizedArrayType *out_array) const;
702
708 void
709 fill_inverse_JxW_values(
710 AlignedVector<VectorizedArrayType> &inverse_jxw) const;
711
712 private:
716 const FEEvaluationBase<dim,
717 n_components,
718 Number,
719 false,
720 VectorizedArrayType> &fe_eval;
721 };
722
723
724
732 template <int dim,
733 int fe_degree,
734 int n_q_points_1d = fe_degree + 1,
735 int n_components = 1,
736 typename VectorType = LinearAlgebra::distributed::Vector<double>,
737 typename VectorizedArrayType =
738 VectorizedArray<typename VectorType::value_type>>
739 class MassOperator : public Base<dim, VectorType, VectorizedArrayType>
740 {
741 public:
745 using value_type =
746 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
747
751 using size_type =
752 typename Base<dim, VectorType, VectorizedArrayType>::size_type;
753
757 MassOperator();
758
763 virtual void
764 compute_diagonal() override;
765
766 private:
772 virtual void
773 apply_add(VectorType &dst, const VectorType &src) const override;
774
778 void
779 local_apply_cell(
780 const MatrixFree<dim, value_type, VectorizedArrayType> &data,
781 VectorType & dst,
782 const VectorType & src,
783 const std::pair<unsigned int, unsigned int> & cell_range) const;
784 };
785
786
787
798 template <int dim,
799 int fe_degree,
800 int n_q_points_1d = fe_degree + 1,
801 int n_components = 1,
802 typename VectorType = LinearAlgebra::distributed::Vector<double>,
803 typename VectorizedArrayType =
804 VectorizedArray<typename VectorType::value_type>>
805 class LaplaceOperator : public Base<dim, VectorType, VectorizedArrayType>
806 {
807 public:
811 using value_type =
812 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
813
817 using size_type =
818 typename Base<dim, VectorType, VectorizedArrayType>::size_type;
819
823 LaplaceOperator();
824
831 virtual void
832 compute_diagonal() override;
833
884 void
885 set_coefficient(
886 const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient);
887
892 virtual void
893 clear() override;
894
901 std::shared_ptr<Table<2, VectorizedArrayType>>
902 get_coefficient();
903
904 private:
910 virtual void
911 apply_add(VectorType &dst, const VectorType &src) const override;
912
916 void
917 local_apply_cell(
918 const MatrixFree<dim, value_type, VectorizedArrayType> &data,
919 VectorType & dst,
920 const VectorType & src,
921 const std::pair<unsigned int, unsigned int> & cell_range) const;
922
926 void
927 local_diagonal_cell(
928 const MatrixFree<dim, value_type, VectorizedArrayType> &data,
929 VectorType & dst,
930 const VectorType &,
931 const std::pair<unsigned int, unsigned int> &cell_range) const;
932
936 void
937 do_operation_on_cell(
938 FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, value_type>
939 & phi,
940 const unsigned int cell) const;
941
945 std::shared_ptr<Table<2, VectorizedArrayType>> scalar_coefficient;
946 };
947
948
949
950 // ------------------------------------ inline functions ---------------------
951
952 template <int dim,
953 int fe_degree,
954 int n_components,
955 typename Number,
956 typename VectorizedArrayType>
957 inline CellwiseInverseMassMatrix<dim,
958 fe_degree,
959 n_components,
960 Number,
961 VectorizedArrayType>::
962 CellwiseInverseMassMatrix(
963 const FEEvaluationBase<dim,
964 n_components,
965 Number,
966 false,
967 VectorizedArrayType> &fe_eval)
968 : fe_eval(fe_eval)
969 {}
970
971
972
973 template <int dim,
974 int fe_degree,
975 int n_components,
976 typename Number,
977 typename VectorizedArrayType>
978 inline void
979 CellwiseInverseMassMatrix<dim,
980 fe_degree,
981 n_components,
982 Number,
983 VectorizedArrayType>::
984 fill_inverse_JxW_values(
985 AlignedVector<VectorizedArrayType> &inverse_jxw) const
986 {
987 constexpr unsigned int dofs_per_component_on_cell =
988 Utilities::pow(fe_degree + 1, dim);
989 Assert(inverse_jxw.size() > 0 &&
990 inverse_jxw.size() % dofs_per_component_on_cell == 0,
991 ExcMessage(
992 "Expected diagonal to be a multiple of scalar dof per cells"));
993
994 // compute values for the first component
995 for (unsigned int q = 0; q < dofs_per_component_on_cell; ++q)
996 inverse_jxw[q] = 1. / fe_eval.JxW(q);
997 // copy values to rest of vector
998 for (unsigned int q = dofs_per_component_on_cell; q < inverse_jxw.size();)
999 for (unsigned int i = 0; i < dofs_per_component_on_cell; ++i, ++q)
1000 inverse_jxw[q] = inverse_jxw[i];
1001 }
1002
1003
1004
1005 template <int dim,
1006 int fe_degree,
1007 int n_components,
1008 typename Number,
1009 typename VectorizedArrayType>
1010 inline void
1011 CellwiseInverseMassMatrix<
1012 dim,
1013 fe_degree,
1014 n_components,
1015 Number,
1016 VectorizedArrayType>::apply(const VectorizedArrayType *in_array,
1017 VectorizedArrayType * out_array) const
1018 {
1019 internal::CellwiseInverseMassMatrixImplBasic<dim, VectorizedArrayType>::
1020 template run<fe_degree>(n_components, fe_eval, in_array, out_array);
1021 }
1022
1023
1024
1025 template <int dim,
1026 int fe_degree,
1027 int n_components,
1028 typename Number,
1029 typename VectorizedArrayType>
1030 inline void
1031 CellwiseInverseMassMatrix<dim,
1032 fe_degree,
1033 n_components,
1034 Number,
1035 VectorizedArrayType>::
1036 apply(const AlignedVector<VectorizedArrayType> &inverse_coefficients,
1037 const unsigned int n_actual_components,
1038 const VectorizedArrayType * in_array,
1039 VectorizedArrayType * out_array) const
1040 {
1041 internal::CellwiseInverseMassMatrixImplFlexible<dim, VectorizedArrayType>::
1042 template run<fe_degree>(
1043 n_actual_components,
1044 fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
1045 inverse_coefficients,
1046 in_array,
1047 out_array);
1048 }
1049
1050
1051
1052 template <int dim,
1053 int fe_degree,
1054 int n_components,
1055 typename Number,
1056 typename VectorizedArrayType>
1057 inline void
1058 CellwiseInverseMassMatrix<dim,
1059 fe_degree,
1060 n_components,
1061 Number,
1062 VectorizedArrayType>::
1063 transform_from_q_points_to_basis(const unsigned int n_actual_components,
1064 const VectorizedArrayType *in_array,
1065 VectorizedArrayType * out_array) const
1066 {
1067 internal::CellwiseInverseMassMatrixImplTransformFromQPoints<
1068 dim,
1069 VectorizedArrayType>::template run<fe_degree>(n_actual_components,
1070 fe_eval.get_shape_info()
1071 .data.front()
1072 .inverse_shape_values_eo,
1073 in_array,
1074 out_array);
1075 }
1076
1077
1078
1079 //----------------- Base operator -----------------------------
1080 template <int dim, typename VectorType, typename VectorizedArrayType>
1081 Base<dim, VectorType, VectorizedArrayType>::Base()
1082 : Subscriptor()
1083 , have_interface_matrices(false)
1084 {}
1085
1086
1087
1088 template <int dim, typename VectorType, typename VectorizedArrayType>
1089 typename Base<dim, VectorType, VectorizedArrayType>::size_type
1090 Base<dim, VectorType, VectorizedArrayType>::m() const
1091 {
1092 Assert(data.get() != nullptr, ExcNotInitialized());
1093 typename Base<dim, VectorType, VectorizedArrayType>::size_type total_size =
1094 0;
1095 for (const unsigned int selected_row : selected_rows)
1096 total_size += data->get_vector_partitioner(selected_row)->size();
1097 return total_size;
1098 }
1099
1100
1101
1102 template <int dim, typename VectorType, typename VectorizedArrayType>
1103 typename Base<dim, VectorType, VectorizedArrayType>::size_type
1104 Base<dim, VectorType, VectorizedArrayType>::n() const
1105 {
1106 Assert(data.get() != nullptr, ExcNotInitialized());
1107 typename Base<dim, VectorType, VectorizedArrayType>::size_type total_size =
1108 0;
1109 for (const unsigned int selected_column : selected_columns)
1110 total_size += data->get_vector_partitioner(selected_column)->size();
1111 return total_size;
1112 }
1113
1114
1115
1116 template <int dim, typename VectorType, typename VectorizedArrayType>
1117 void
1118 Base<dim, VectorType, VectorizedArrayType>::clear()
1119 {
1120 data.reset();
1121 inverse_diagonal_entries.reset();
1122 }
1123
1124
1125
1126 template <int dim, typename VectorType, typename VectorizedArrayType>
1127 typename Base<dim, VectorType, VectorizedArrayType>::value_type
1128 Base<dim, VectorType, VectorizedArrayType>::el(const unsigned int row,
1129 const unsigned int col) const
1130 {
1131 (void)col;
1132 Assert(row == col, ExcNotImplemented());
1133 Assert(inverse_diagonal_entries.get() != nullptr &&
1134 inverse_diagonal_entries->m() > 0,
1135 ExcNotInitialized());
1136 return 1.0 / (*inverse_diagonal_entries)(row, row);
1137 }
1138
1139
1140
1141 template <int dim, typename VectorType, typename VectorizedArrayType>
1142 void
1143 Base<dim, VectorType, VectorizedArrayType>::initialize_dof_vector(
1144 VectorType &vec) const
1145 {
1146 Assert(data.get() != nullptr, ExcNotInitialized());
1147 AssertDimension(BlockHelper::n_blocks(vec), selected_rows.size());
1148 for (unsigned int i = 0; i < BlockHelper::n_blocks(vec); ++i)
1149 {
1150 const unsigned int index = selected_rows[i];
1151 if (!BlockHelper::subblock(vec, index)
1152 .partitioners_are_compatible(
1153 *data->get_dof_info(index).vector_partitioner))
1154 data->initialize_dof_vector(BlockHelper::subblock(vec, index), index);
1155
1156 Assert(BlockHelper::subblock(vec, index)
1157 .partitioners_are_globally_compatible(
1158 *data->get_dof_info(index).vector_partitioner),
1159 ExcInternalError());
1160 }
1161 BlockHelper::collect_sizes(vec);
1162 }
1163
1164
1165
1166 template <int dim, typename VectorType, typename VectorizedArrayType>
1167 void
1168 Base<dim, VectorType, VectorizedArrayType>::initialize(
1169 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
1170 data_,
1171 const std::vector<unsigned int> &given_row_selection,
1172 const std::vector<unsigned int> &given_column_selection)
1173 {
1174 data = data_;
1175
1176 selected_rows.clear();
1177 selected_columns.clear();
1178 if (given_row_selection.empty())
1179 for (unsigned int i = 0; i < data_->n_components(); ++i)
1180 selected_rows.push_back(i);
1181 else
1182 {
1183 for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1184 {
1185 AssertIndexRange(given_row_selection[i], data_->n_components());
1186 for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1187 if (j != i)
1188 Assert(given_row_selection[j] != given_row_selection[i],
1189 ExcMessage("Given row indices must be unique"));
1190
1191 selected_rows.push_back(given_row_selection[i]);
1192 }
1193 }
1194 if (given_column_selection.size() == 0)
1195 selected_columns = selected_rows;
1196 else
1197 {
1198 for (unsigned int i = 0; i < given_column_selection.size(); ++i)
1199 {
1200 AssertIndexRange(given_column_selection[i], data_->n_components());
1201 for (unsigned int j = 0; j < given_column_selection.size(); ++j)
1202 if (j != i)
1203 Assert(given_column_selection[j] != given_column_selection[i],
1204 ExcMessage("Given column indices must be unique"));
1205
1206 selected_columns.push_back(given_column_selection[i]);
1207 }
1208 }
1209
1210 edge_constrained_indices.clear();
1211 edge_constrained_indices.resize(selected_rows.size());
1212 edge_constrained_values.clear();
1213 edge_constrained_values.resize(selected_rows.size());
1214 have_interface_matrices = false;
1215 }
1216
1217
1218
1219 template <int dim, typename VectorType, typename VectorizedArrayType>
1220 void
1221 Base<dim, VectorType, VectorizedArrayType>::initialize(
1222 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
1223 data_,
1224 const MGConstrainedDoFs & mg_constrained_dofs,
1225 const unsigned int level,
1226 const std::vector<unsigned int> &given_row_selection)
1227 {
1228 std::vector<MGConstrainedDoFs> mg_constrained_dofs_vector(
1229 1, mg_constrained_dofs);
1230 initialize(data_, mg_constrained_dofs_vector, level, given_row_selection);
1231 }
1232
1233
1234
1235 template <int dim, typename VectorType, typename VectorizedArrayType>
1236 void
1237 Base<dim, VectorType, VectorizedArrayType>::initialize(
1238 std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
1239 data_,
1240 const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
1241 const unsigned int level,
1242 const std::vector<unsigned int> & given_row_selection)
1243 {
1244 AssertThrow(level != numbers::invalid_unsigned_int,
1245 ExcMessage("level is not set"));
1246
1247 selected_rows.clear();
1248 selected_columns.clear();
1249 if (given_row_selection.empty())
1250 for (unsigned int i = 0; i < data_->n_components(); ++i)
1251 selected_rows.push_back(i);
1252 else
1253 {
1254 for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1255 {
1256 AssertIndexRange(given_row_selection[i], data_->n_components());
1257 for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1258 if (j != i)
1259 Assert(given_row_selection[j] != given_row_selection[i],
1260 ExcMessage("Given row indices must be unique"));
1261
1262 selected_rows.push_back(given_row_selection[i]);
1263 }
1264 }
1265 selected_columns = selected_rows;
1266
1267 AssertDimension(mg_constrained_dofs.size(), selected_rows.size());
1268 edge_constrained_indices.clear();
1269 edge_constrained_indices.resize(selected_rows.size());
1270 edge_constrained_values.clear();
1271 edge_constrained_values.resize(selected_rows.size());
1272
1273 data = data_;
1274
1275 for (unsigned int j = 0; j < selected_rows.size(); ++j)
1276 {
1277 if (data_->n_cell_batches() > 0)
1278 {
1279 AssertDimension(level, data_->get_cell_iterator(0, 0, j)->level());
1280 }
1281
1282 // setup edge_constrained indices
1283 std::vector<types::global_dof_index> interface_indices;
1284 mg_constrained_dofs[j]
1285 .get_refinement_edge_indices(level)
1286 .fill_index_vector(interface_indices);
1287 edge_constrained_indices[j].clear();
1288 edge_constrained_indices[j].reserve(interface_indices.size());
1289 edge_constrained_values[j].resize(interface_indices.size());
1290 const IndexSet &locally_owned =
1291 data->get_dof_handler(selected_rows[j]).locally_owned_mg_dofs(level);
1292 for (const auto interface_index : interface_indices)
1293 if (locally_owned.is_element(interface_index))
1294 edge_constrained_indices[j].push_back(
1295 locally_owned.index_within_set(interface_index));
1296 have_interface_matrices |=
1297 Utilities::MPI::max(
1298 static_cast<unsigned int>(edge_constrained_indices[j].size()),
1299 data->get_vector_partitioner()->get_mpi_communicator()) > 0;
1300 }
1301 }
1302
1303
1304
1305 template <int dim, typename VectorType, typename VectorizedArrayType>
1306 void
1307 Base<dim, VectorType, VectorizedArrayType>::set_constrained_entries_to_one(
1308 VectorType &dst) const
1309 {
1310 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1311 {
1312 const std::vector<unsigned int> &constrained_dofs =
1313 data->get_constrained_dofs(selected_rows[j]);
1314 for (const auto constrained_dof : constrained_dofs)
1315 BlockHelper::subblock(dst, j).local_element(constrained_dof) = 1.;
1316 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1317 BlockHelper::subblock(dst, j).local_element(
1318 edge_constrained_indices[j][i]) = 1.;
1319 }
1320 }
1321
1322
1323
1324 template <int dim, typename VectorType, typename VectorizedArrayType>
1325 void
1326 Base<dim, VectorType, VectorizedArrayType>::vmult(VectorType & dst,
1327 const VectorType &src) const
1328 {
1329 using Number =
1330 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1331 dst = Number(0.);
1332 vmult_add(dst, src);
1333 }
1334
1335
1336
1337 template <int dim, typename VectorType, typename VectorizedArrayType>
1338 void
1339 Base<dim, VectorType, VectorizedArrayType>::vmult_add(
1340 VectorType & dst,
1341 const VectorType &src) const
1342 {
1343 mult_add(dst, src, false);
1344 }
1345
1346
1347
1348 template <int dim, typename VectorType, typename VectorizedArrayType>
1349 void
1350 Base<dim, VectorType, VectorizedArrayType>::Tvmult_add(
1351 VectorType & dst,
1352 const VectorType &src) const
1353 {
1354 mult_add(dst, src, true);
1355 }
1356
1357
1358
1359 template <int dim, typename VectorType, typename VectorizedArrayType>
1360 void
1361 Base<dim, VectorType, VectorizedArrayType>::adjust_ghost_range_if_necessary(
1362 const VectorType &src,
1363 const bool is_row) const
1364 {
1365 using Number =
1366 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1367 for (unsigned int i = 0; i < BlockHelper::n_blocks(src); ++i)
1368 {
1369 const unsigned int mf_component =
1370 is_row ? selected_rows[i] : selected_columns[i];
1371 // If both vectors use the same partitioner -> done
1372 if (BlockHelper::subblock(src, i).get_partitioner().get() ==
1373 data->get_dof_info(mf_component).vector_partitioner.get())
1374 continue;
1375
1376 // If not, assert that the local ranges are the same and reset to the
1377 // current partitioner
1378 Assert(BlockHelper::subblock(src, i)
1379 .get_partitioner()
1380 ->locally_owned_size() ==
1381 data->get_dof_info(mf_component)
1382 .vector_partitioner->locally_owned_size(),
1383 ExcMessage(
1384 "The vector passed to the vmult() function does not have "
1385 "the correct size for compatibility with MatrixFree."));
1386
1387 // copy the vector content to a temporary vector so that it does not get
1388 // lost
1389 LinearAlgebra::distributed::Vector<Number> copy_vec(
1390 BlockHelper::subblock(src, i));
1391 this->data->initialize_dof_vector(
1392 BlockHelper::subblock(const_cast<VectorType &>(src), i),
1393 mf_component);
1394 BlockHelper::subblock(const_cast<VectorType &>(src), i)
1395 .copy_locally_owned_data_from(copy_vec);
1396 }
1397 }
1398
1399
1400
1401 template <int dim, typename VectorType, typename VectorizedArrayType>
1402 void
1403 Base<dim, VectorType, VectorizedArrayType>::preprocess_constraints(
1404 VectorType & dst,
1405 const VectorType &src) const
1406 {
1407 using Number =
1408 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1409 adjust_ghost_range_if_necessary(src, false);
1410 adjust_ghost_range_if_necessary(dst, true);
1411
1412 // set zero Dirichlet values on the input vector (and remember the src and
1413 // dst values because we need to reset them at the end)
1414 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1415 {
1416 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1417 {
1418 edge_constrained_values[j][i] = std::pair<Number, Number>(
1419 BlockHelper::subblock(src, j).local_element(
1420 edge_constrained_indices[j][i]),
1421 BlockHelper::subblock(dst, j).local_element(
1422 edge_constrained_indices[j][i]));
1423 BlockHelper::subblock(const_cast<VectorType &>(src), j)
1424 .local_element(edge_constrained_indices[j][i]) = 0.;
1425 }
1426 }
1427 }
1428
1429
1430
1431 template <int dim, typename VectorType, typename VectorizedArrayType>
1432 void
1433 Base<dim, VectorType, VectorizedArrayType>::mult_add(
1434 VectorType & dst,
1435 const VectorType &src,
1436 const bool transpose) const
1437 {
1438 AssertDimension(dst.size(), src.size());
1439 AssertDimension(BlockHelper::n_blocks(dst), BlockHelper::n_blocks(src));
1440 AssertDimension(BlockHelper::n_blocks(dst), selected_rows.size());
1441 preprocess_constraints(dst, src);
1442 if (transpose)
1443 Tapply_add(dst, src);
1444 else
1445 apply_add(dst, src);
1446 postprocess_constraints(dst, src);
1447 }
1448
1449
1450
1451 template <int dim, typename VectorType, typename VectorizedArrayType>
1452 void
1453 Base<dim, VectorType, VectorizedArrayType>::postprocess_constraints(
1454 VectorType & dst,
1455 const VectorType &src) const
1456 {
1457 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1458 {
1459 const std::vector<unsigned int> &constrained_dofs =
1460 data->get_constrained_dofs(selected_rows[j]);
1461 for (const auto constrained_dof : constrained_dofs)
1462 BlockHelper::subblock(dst, j).local_element(constrained_dof) +=
1463 BlockHelper::subblock(src, j).local_element(constrained_dof);
1464 }
1465
1466 // reset edge constrained values, multiply by unit matrix and add into
1467 // destination
1468 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1469 {
1470 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1471 {
1472 BlockHelper::subblock(const_cast<VectorType &>(src), j)
1473 .local_element(edge_constrained_indices[j][i]) =
1474 edge_constrained_values[j][i].first;
1475 BlockHelper::subblock(dst, j).local_element(
1476 edge_constrained_indices[j][i]) =
1477 edge_constrained_values[j][i].second +
1478 edge_constrained_values[j][i].first;
1479 }
1480 }
1481 }
1482
1483
1484
1485 template <int dim, typename VectorType, typename VectorizedArrayType>
1486 void
1487 Base<dim, VectorType, VectorizedArrayType>::vmult_interface_down(
1488 VectorType & dst,
1489 const VectorType &src) const
1490 {
1491 using Number =
1492 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1493 AssertDimension(dst.size(), src.size());
1494 adjust_ghost_range_if_necessary(src, false);
1495 adjust_ghost_range_if_necessary(dst, true);
1496
1497 dst = Number(0.);
1498
1499 if (!have_interface_matrices)
1500 return;
1501
1502 // set zero Dirichlet values on the input vector (and remember the src and
1503 // dst values because we need to reset them at the end)
1504 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1505 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1506 {
1507 edge_constrained_values[j][i] = std::pair<Number, Number>(
1508 BlockHelper::subblock(src, j).local_element(
1509 edge_constrained_indices[j][i]),
1510 BlockHelper::subblock(dst, j).local_element(
1511 edge_constrained_indices[j][i]));
1512 BlockHelper::subblock(const_cast<VectorType &>(src), j)
1513 .local_element(edge_constrained_indices[j][i]) = 0.;
1514 }
1515
1516 apply_add(dst, src);
1517
1518 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1519 {
1520 unsigned int c = 0;
1521 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1522 {
1523 for (; c < edge_constrained_indices[j][i]; ++c)
1524 BlockHelper::subblock(dst, j).local_element(c) = 0.;
1525 ++c;
1526
1527 // reset the src values
1528 BlockHelper::subblock(const_cast<VectorType &>(src), j)
1529 .local_element(edge_constrained_indices[j][i]) =
1530 edge_constrained_values[j][i].first;
1531 }
1532 for (; c < BlockHelper::subblock(dst, j).locally_owned_size(); ++c)
1533 BlockHelper::subblock(dst, j).local_element(c) = 0.;
1534 }
1535 }
1536
1537
1538
1539 template <int dim, typename VectorType, typename VectorizedArrayType>
1540 void
1541 Base<dim, VectorType, VectorizedArrayType>::vmult_interface_up(
1542 VectorType & dst,
1543 const VectorType &src) const
1544 {
1545 using Number =
1546 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1547 AssertDimension(dst.size(), src.size());
1548 adjust_ghost_range_if_necessary(src, false);
1549 adjust_ghost_range_if_necessary(dst, true);
1550
1551 dst = Number(0.);
1552
1553 if (!have_interface_matrices)
1554 return;
1555
1556 VectorType src_cpy(src);
1557 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1558 {
1559 unsigned int c = 0;
1560 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1561 {
1562 for (; c < edge_constrained_indices[j][i]; ++c)
1563 BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1564 ++c;
1565 }
1566 for (; c < BlockHelper::subblock(src_cpy, j).locally_owned_size(); ++c)
1567 BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1568 }
1569
1570 apply_add(dst, src_cpy);
1571
1572 for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1573 for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1574 BlockHelper::subblock(dst, j).local_element(
1575 edge_constrained_indices[j][i]) = 0.;
1576 }
1577
1578
1579
1580 template <int dim, typename VectorType, typename VectorizedArrayType>
1581 void
1582 Base<dim, VectorType, VectorizedArrayType>::Tvmult(
1583 VectorType & dst,
1584 const VectorType &src) const
1585 {
1586 using Number =
1587 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1588 dst = Number(0.);
1589 Tvmult_add(dst, src);
1590 }
1591
1592
1593
1594 template <int dim, typename VectorType, typename VectorizedArrayType>
1595 std::size_t
1596 Base<dim, VectorType, VectorizedArrayType>::memory_consumption() const
1597 {
1598 return inverse_diagonal_entries.get() != nullptr ?
1599 inverse_diagonal_entries->memory_consumption() :
1600 sizeof(*this);
1601 }
1602
1603
1604
1605 template <int dim, typename VectorType, typename VectorizedArrayType>
1606 std::shared_ptr<const MatrixFree<
1607 dim,
1608 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
1609 VectorizedArrayType>>
1610 Base<dim, VectorType, VectorizedArrayType>::get_matrix_free() const
1611 {
1612 return data;
1613 }
1614
1615
1616
1617 template <int dim, typename VectorType, typename VectorizedArrayType>
1618 const std::shared_ptr<DiagonalMatrix<VectorType>> &
1619 Base<dim, VectorType, VectorizedArrayType>::get_matrix_diagonal_inverse()
1620 const
1621 {
1622 Assert(inverse_diagonal_entries.get() != nullptr &&
1623 inverse_diagonal_entries->m() > 0,
1624 ExcNotInitialized());
1625 return inverse_diagonal_entries;
1626 }
1627
1628
1629
1630 template <int dim, typename VectorType, typename VectorizedArrayType>
1631 const std::shared_ptr<DiagonalMatrix<VectorType>> &
1632 Base<dim, VectorType, VectorizedArrayType>::get_matrix_diagonal() const
1633 {
1634 Assert(diagonal_entries.get() != nullptr && diagonal_entries->m() > 0,
1635 ExcNotInitialized());
1636 return diagonal_entries;
1637 }
1638
1639
1640
1641 template <int dim, typename VectorType, typename VectorizedArrayType>
1642 void
1643 Base<dim, VectorType, VectorizedArrayType>::Tapply_add(
1644 VectorType & dst,
1645 const VectorType &src) const
1646 {
1647 apply_add(dst, src);
1648 }
1649
1650
1651
1652 template <int dim, typename VectorType, typename VectorizedArrayType>
1653 void
1654 Base<dim, VectorType, VectorizedArrayType>::precondition_Jacobi(
1655 VectorType & dst,
1656 const VectorType & src,
1657 const typename Base<dim, VectorType, VectorizedArrayType>::value_type omega)
1658 const
1659 {
1660 Assert(inverse_diagonal_entries.get() && inverse_diagonal_entries->m() > 0,
1661 ExcNotInitialized());
1662 inverse_diagonal_entries->vmult(dst, src);
1663 dst *= omega;
1664 }
1665
1666
1667
1668 //------------------------- MGInterfaceOperator ------------------------------
1669
1670 template <typename OperatorType>
1671 MGInterfaceOperator<OperatorType>::MGInterfaceOperator()
1672 : Subscriptor()
1673 , mf_base_operator(nullptr)
1674 {}
1675
1676
1677
1678 template <typename OperatorType>
1679 void
1680 MGInterfaceOperator<OperatorType>::clear()
1681 {
1682 mf_base_operator = nullptr;
1683 }
1684
1685
1686
1687 template <typename OperatorType>
1688 void
1689 MGInterfaceOperator<OperatorType>::initialize(const OperatorType &operator_in)
1690 {
1691 mf_base_operator = &operator_in;
1692 }
1693
1694
1695
1696 template <typename OperatorType>
1697 template <typename VectorType>
1698 void
1699 MGInterfaceOperator<OperatorType>::vmult(VectorType & dst,
1700 const VectorType &src) const
1701 {
1702#ifndef DEAL_II_MSVC
1703 static_assert(
1704 std::is_same<typename VectorType::value_type, value_type>::value,
1705 "The vector type must be based on the same value type as this "
1706 "operator");
1707#endif
1708
1709 Assert(mf_base_operator != nullptr, ExcNotInitialized());
1710
1711 mf_base_operator->vmult_interface_down(dst, src);
1712 }
1713
1714
1715
1716 template <typename OperatorType>
1717 template <typename VectorType>
1718 void
1719 MGInterfaceOperator<OperatorType>::Tvmult(VectorType & dst,
1720 const VectorType &src) const
1721 {
1722#ifndef DEAL_II_MSVC
1723 static_assert(
1724 std::is_same<typename VectorType::value_type, value_type>::value,
1725 "The vector type must be based on the same value type as this "
1726 "operator");
1727#endif
1728
1729 Assert(mf_base_operator != nullptr, ExcNotInitialized());
1730
1731 mf_base_operator->vmult_interface_up(dst, src);
1732 }
1733
1734
1735
1736 template <typename OperatorType>
1737 template <typename VectorType>
1738 void
1739 MGInterfaceOperator<OperatorType>::initialize_dof_vector(
1740 VectorType &vec) const
1741 {
1742 Assert(mf_base_operator != nullptr, ExcNotInitialized());
1743
1744 mf_base_operator->initialize_dof_vector(vec);
1745 }
1746
1747
1748
1749 //-----------------------------MassOperator----------------------------------
1750
1751 template <int dim,
1752 int fe_degree,
1753 int n_q_points_1d,
1754 int n_components,
1755 typename VectorType,
1756 typename VectorizedArrayType>
1757 MassOperator<dim,
1758 fe_degree,
1759 n_q_points_1d,
1760 n_components,
1761 VectorType,
1762 VectorizedArrayType>::MassOperator()
1763 : Base<dim, VectorType, VectorizedArrayType>()
1764 {}
1765
1766
1767
1768 template <int dim,
1769 int fe_degree,
1770 int n_q_points_1d,
1771 int n_components,
1772 typename VectorType,
1773 typename VectorizedArrayType>
1774 void
1775 MassOperator<dim,
1776 fe_degree,
1777 n_q_points_1d,
1778 n_components,
1779 VectorType,
1780 VectorizedArrayType>::compute_diagonal()
1781 {
1782 using Number =
1783 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1784 Assert((Base<dim, VectorType, VectorizedArrayType>::data.get() != nullptr),
1785 ExcNotInitialized());
1786
1787 this->inverse_diagonal_entries =
1788 std::make_shared<DiagonalMatrix<VectorType>>();
1789 this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1790 VectorType &inverse_diagonal_vector =
1791 this->inverse_diagonal_entries->get_vector();
1792 VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1793 this->initialize_dof_vector(inverse_diagonal_vector);
1794 this->initialize_dof_vector(diagonal_vector);
1795 inverse_diagonal_vector = Number(1.);
1796 apply_add(diagonal_vector, inverse_diagonal_vector);
1797
1798 this->set_constrained_entries_to_one(diagonal_vector);
1799 inverse_diagonal_vector = diagonal_vector;
1800
1801 const unsigned int locally_owned_size =
1802 inverse_diagonal_vector.locally_owned_size();
1803 for (unsigned int i = 0; i < locally_owned_size; ++i)
1804 inverse_diagonal_vector.local_element(i) =
1805 Number(1.) / inverse_diagonal_vector.local_element(i);
1806
1807 inverse_diagonal_vector.update_ghost_values();
1808 diagonal_vector.update_ghost_values();
1809 }
1810
1811
1812
1813 template <int dim,
1814 int fe_degree,
1815 int n_q_points_1d,
1816 int n_components,
1817 typename VectorType,
1818 typename VectorizedArrayType>
1819 void
1820 MassOperator<dim,
1821 fe_degree,
1822 n_q_points_1d,
1823 n_components,
1824 VectorType,
1825 VectorizedArrayType>::apply_add(VectorType & dst,
1826 const VectorType &src) const
1827 {
1828 Base<dim, VectorType, VectorizedArrayType>::data->cell_loop(
1829 &MassOperator::local_apply_cell, this, dst, src);
1830 }
1831
1832
1833
1834 template <int dim,
1835 int fe_degree,
1836 int n_q_points_1d,
1837 int n_components,
1838 typename VectorType,
1839 typename VectorizedArrayType>
1840 void
1841 MassOperator<dim,
1842 fe_degree,
1843 n_q_points_1d,
1844 n_components,
1845 VectorType,
1846 VectorizedArrayType>::
1847 local_apply_cell(
1848 const MatrixFree<
1849 dim,
1850 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
1851 VectorizedArrayType> & data,
1852 VectorType & dst,
1853 const VectorType & src,
1854 const std::pair<unsigned int, unsigned int> &cell_range) const
1855 {
1856 using Number =
1857 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1858 FEEvaluation<dim,
1859 fe_degree,
1860 n_q_points_1d,
1861 n_components,
1862 Number,
1863 VectorizedArrayType>
1864 phi(data, this->selected_rows[0]);
1865 for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
1866 {
1867 phi.reinit(cell);
1868 phi.read_dof_values(src);
1869 phi.evaluate(EvaluationFlags::values);
1870 for (unsigned int q = 0; q < phi.n_q_points; ++q)
1871 phi.submit_value(phi.get_value(q), q);
1872 phi.integrate(EvaluationFlags::values);
1873 phi.distribute_local_to_global(dst);
1874 }
1875 }
1876
1877
1878 //-----------------------------LaplaceOperator----------------------------------
1879
1880 template <int dim,
1881 int fe_degree,
1882 int n_q_points_1d,
1883 int n_components,
1884 typename VectorType,
1885 typename VectorizedArrayType>
1886 LaplaceOperator<dim,
1887 fe_degree,
1888 n_q_points_1d,
1889 n_components,
1890 VectorType,
1891 VectorizedArrayType>::LaplaceOperator()
1892 : Base<dim, VectorType, VectorizedArrayType>()
1893 {}
1894
1895
1896
1897 template <int dim,
1898 int fe_degree,
1899 int n_q_points_1d,
1900 int n_components,
1901 typename VectorType,
1902 typename VectorizedArrayType>
1903 void
1904 LaplaceOperator<dim,
1905 fe_degree,
1906 n_q_points_1d,
1907 n_components,
1908 VectorType,
1909 VectorizedArrayType>::clear()
1910 {
1911 Base<dim, VectorType, VectorizedArrayType>::clear();
1912 scalar_coefficient.reset();
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 LaplaceOperator<dim,
1925 fe_degree,
1926 n_q_points_1d,
1927 n_components,
1928 VectorType,
1929 VectorizedArrayType>::
1930 set_coefficient(
1931 const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient_)
1932 {
1933 scalar_coefficient = scalar_coefficient_;
1934 }
1935
1936
1937
1938 template <int dim,
1939 int fe_degree,
1940 int n_q_points_1d,
1941 int n_components,
1942 typename VectorType,
1943 typename VectorizedArrayType>
1944 std::shared_ptr<Table<2, VectorizedArrayType>>
1945 LaplaceOperator<dim,
1946 fe_degree,
1947 n_q_points_1d,
1948 n_components,
1949 VectorType,
1950 VectorizedArrayType>::get_coefficient()
1951 {
1952 Assert(scalar_coefficient.get(), ExcNotInitialized());
1953 return scalar_coefficient;
1954 }
1955
1956
1957
1958 template <int dim,
1959 int fe_degree,
1960 int n_q_points_1d,
1961 int n_components,
1962 typename VectorType,
1963 typename VectorizedArrayType>
1964 void
1965 LaplaceOperator<dim,
1966 fe_degree,
1967 n_q_points_1d,
1968 n_components,
1969 VectorType,
1970 VectorizedArrayType>::compute_diagonal()
1971 {
1972 using Number =
1973 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1974 Assert((Base<dim, VectorType, VectorizedArrayType>::data.get() != nullptr),
1975 ExcNotInitialized());
1976
1977 this->inverse_diagonal_entries =
1978 std::make_shared<DiagonalMatrix<VectorType>>();
1979 this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1980 VectorType &inverse_diagonal_vector =
1981 this->inverse_diagonal_entries->get_vector();
1982 VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1983 this->initialize_dof_vector(inverse_diagonal_vector);
1984 this->initialize_dof_vector(diagonal_vector);
1985
1986 this->data->cell_loop(&LaplaceOperator::local_diagonal_cell,
1987 this,
1988 diagonal_vector,
1989 /*unused*/ diagonal_vector);
1990 this->set_constrained_entries_to_one(diagonal_vector);
1991
1992 inverse_diagonal_vector = diagonal_vector;
1993
1994 for (unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
1995 ++i)
1996 if (std::abs(inverse_diagonal_vector.local_element(i)) >
1997 std::sqrt(std::numeric_limits<Number>::epsilon()))
1998 inverse_diagonal_vector.local_element(i) =
1999 1. / inverse_diagonal_vector.local_element(i);
2000 else
2001 inverse_diagonal_vector.local_element(i) = 1.;
2002
2003 inverse_diagonal_vector.update_ghost_values();
2004 diagonal_vector.update_ghost_values();
2005 }
2006
2007
2008
2009 template <int dim,
2010 int fe_degree,
2011 int n_q_points_1d,
2012 int n_components,
2013 typename VectorType,
2014 typename VectorizedArrayType>
2015 void
2016 LaplaceOperator<dim,
2017 fe_degree,
2018 n_q_points_1d,
2019 n_components,
2020 VectorType,
2021 VectorizedArrayType>::apply_add(VectorType & dst,
2022 const VectorType &src) const
2023 {
2024 Base<dim, VectorType, VectorizedArrayType>::data->cell_loop(
2025 &LaplaceOperator::local_apply_cell, this, dst, src);
2026 }
2027
2028 namespace Implementation
2029 {
2030 template <typename VectorizedArrayType>
2031 bool
2032 non_negative(const VectorizedArrayType &n)
2033 {
2034 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2035 if (n[v] < 0.)
2036 return false;
2037
2038 return true;
2039 }
2040 } // namespace Implementation
2041
2042
2043
2044 template <int dim,
2045 int fe_degree,
2046 int n_q_points_1d,
2047 int n_components,
2048 typename VectorType,
2049 typename VectorizedArrayType>
2050 void
2051 LaplaceOperator<dim,
2052 fe_degree,
2053 n_q_points_1d,
2054 n_components,
2055 VectorType,
2056 VectorizedArrayType>::
2057 do_operation_on_cell(
2058 FEEvaluation<
2059 dim,
2060 fe_degree,
2061 n_q_points_1d,
2062 n_components,
2063 typename Base<dim, VectorType, VectorizedArrayType>::value_type> &phi,
2064 const unsigned int cell) const
2065 {
2066 phi.evaluate(EvaluationFlags::gradients);
2067 if (scalar_coefficient.get())
2068 {
2069 Assert(scalar_coefficient->size(1) == 1 ||
2070 scalar_coefficient->size(1) == phi.n_q_points,
2071 ExcMessage("The number of columns in the coefficient table must "
2072 "be either 1 or the number of quadrature points " +
2073 std::to_string(phi.n_q_points) +
2074 ", but the given value was " +
2075 std::to_string(scalar_coefficient->size(1))));
2076 if (scalar_coefficient->size(1) == phi.n_q_points)
2077 for (unsigned int q = 0; q < phi.n_q_points; ++q)
2078 {
2079 Assert(Implementation::non_negative(
2080 (*scalar_coefficient)(cell, q)),
2081 ExcMessage("Coefficient must be non-negative"));
2082 phi.submit_gradient((*scalar_coefficient)(cell, q) *
2083 phi.get_gradient(q),
2084 q);
2085 }
2086 else
2087 {
2088 Assert(Implementation::non_negative((*scalar_coefficient)(cell, 0)),
2089 ExcMessage("Coefficient must be non-negative"));
2090 const VectorizedArrayType coefficient =
2091 (*scalar_coefficient)(cell, 0);
2092 for (unsigned int q = 0; q < phi.n_q_points; ++q)
2093 phi.submit_gradient(coefficient * phi.get_gradient(q), q);
2094 }
2095 }
2096 else
2097 {
2098 for (unsigned int q = 0; q < phi.n_q_points; ++q)
2099 {
2100 phi.submit_gradient(phi.get_gradient(q), q);
2101 }
2102 }
2103 phi.integrate(EvaluationFlags::gradients);
2104 }
2105
2106
2107
2108 template <int dim,
2109 int fe_degree,
2110 int n_q_points_1d,
2111 int n_components,
2112 typename VectorType,
2113 typename VectorizedArrayType>
2114 void
2115 LaplaceOperator<dim,
2116 fe_degree,
2117 n_q_points_1d,
2118 n_components,
2119 VectorType,
2120 VectorizedArrayType>::
2121 local_apply_cell(
2122 const MatrixFree<
2123 dim,
2124 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2125 VectorizedArrayType> & data,
2126 VectorType & dst,
2127 const VectorType & src,
2128 const std::pair<unsigned int, unsigned int> &cell_range) const
2129 {
2130 using Number =
2131 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2132 FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number> phi(
2133 data, this->selected_rows[0]);
2134 for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2135 {
2136 phi.reinit(cell);
2137 phi.read_dof_values(src);
2138 do_operation_on_cell(phi, cell);
2139 phi.distribute_local_to_global(dst);
2140 }
2141 }
2142
2143
2144 template <int dim,
2145 int fe_degree,
2146 int n_q_points_1d,
2147 int n_components,
2148 typename VectorType,
2149 typename VectorizedArrayType>
2150 void
2151 LaplaceOperator<dim,
2152 fe_degree,
2153 n_q_points_1d,
2154 n_components,
2155 VectorType,
2156 VectorizedArrayType>::
2157 local_diagonal_cell(
2158 const MatrixFree<
2159 dim,
2160 typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2161 VectorizedArrayType> &data,
2162 VectorType & dst,
2163 const VectorType &,
2164 const std::pair<unsigned int, unsigned int> &cell_range) const
2165 {
2166 using Number =
2167 typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2168
2169 FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number> phi(
2170 data, this->selected_rows[0]);
2171 for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2172 {
2173 phi.reinit(cell);
2174 VectorizedArrayType local_diagonal_vector[phi.static_dofs_per_cell];
2175 for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2176 {
2177 for (unsigned int j = 0; j < phi.dofs_per_component; ++j)
2178 phi.begin_dof_values()[j] = VectorizedArrayType();
2179 phi.begin_dof_values()[i] = 1.;
2180 do_operation_on_cell(phi, cell);
2181 local_diagonal_vector[i] = phi.begin_dof_values()[i];
2182 }
2183 for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2184 for (unsigned int c = 0; c < phi.n_components; ++c)
2185 phi.begin_dof_values()[i + c * phi.dofs_per_component] =
2186 local_diagonal_vector[i];
2187 phi.distribute_local_to_global(dst);
2188 }
2189 }
2190
2191
2192} // end of namespace MatrixFreeOperators
2193
2194
2195DEAL_II_NAMESPACE_CLOSE
2196
2197#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::begin_dof_values
const VectorizedArrayType * begin_dof_values() const
FEEvaluationBase::read_dof_values
void read_dof_values(const VectorType &src, const unsigned int first_index=0)
FEEvaluation::dofs_per_component
const unsigned int dofs_per_component
Definition: fe_evaluation.h:2934
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:2597
FEEvaluation::n_components
static constexpr unsigned int n_components
Definition: fe_evaluation.h:2559
IndexSet::index_within_set
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1921
IndexSet::is_element
bool is_element(const size_type index) const
Definition: index_set.h:1765
MatrixFreeOperators::Base::~Base
virtual ~Base() override=default
MatrixFreeOperators::Base::selected_rows
std::vector< unsigned int > selected_rows
Definition: operators.h:451
MatrixFreeOperators::Base::Tvmult_add
void Tvmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1350
MatrixFreeOperators::Base::set_constrained_entries_to_one
void set_constrained_entries_to_one(VectorType &dst) const
Definition: operators.h:1307
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:1339
MatrixFreeOperators::Base::vmult_interface_down
void vmult_interface_down(VectorType &dst, const VectorType &src) const
Definition: operators.h:1487
MatrixFreeOperators::Base::n
size_type n() const
Definition: operators.h:1104
MatrixFreeOperators::Base::preprocess_constraints
void preprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1403
MatrixFreeOperators::Base::mult_add
void mult_add(VectorType &dst, const VectorType &src, const bool transpose) const
Definition: operators.h:1433
MatrixFreeOperators::Base::edge_constrained_indices
std::vector< std::vector< unsigned int > > edge_constrained_indices
Definition: operators.h:463
MatrixFreeOperators::Base::get_matrix_diagonal
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal() const
Definition: operators.h:1632
MatrixFreeOperators::Base::Tvmult
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1582
MatrixFreeOperators::Base::edge_constrained_values
std::vector< std::vector< std::pair< value_type, value_type > > > edge_constrained_values
Definition: operators.h:469
MatrixFreeOperators::Base::m
size_type m() const
Definition: operators.h:1090
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:1168
MatrixFreeOperators::Base::selected_columns
std::vector< unsigned int > selected_columns
Definition: operators.h:457
MatrixFreeOperators::Base::get_matrix_free
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > get_matrix_free() const
Definition: operators.h:1610
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:1237
MatrixFreeOperators::Base::get_matrix_diagonal_inverse
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal_inverse() const
Definition: operators.h:1619
MatrixFreeOperators::Base::clear
virtual void clear()
Definition: operators.h:1118
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:1221
MatrixFreeOperators::Base::initialize_dof_vector
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1143
MatrixFreeOperators::Base::diagonal_entries
std::shared_ptr< DiagonalMatrix< VectorType > > diagonal_entries
Definition: operators.h:439
MatrixFreeOperators::Base::Tapply_add
virtual void Tapply_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1643
MatrixFreeOperators::Base::data
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > data
Definition: operators.h:433
MatrixFreeOperators::Base::adjust_ghost_range_if_necessary
void adjust_ghost_range_if_necessary(const VectorType &vec, const bool is_row) const
Definition: operators.h:1361
MatrixFreeOperators::Base::memory_consumption
virtual std::size_t memory_consumption() const
Definition: operators.h:1596
MatrixFreeOperators::Base::inverse_diagonal_entries
std::shared_ptr< DiagonalMatrix< VectorType > > inverse_diagonal_entries
Definition: operators.h:445
MatrixFreeOperators::Base::el
value_type el(const unsigned int row, const unsigned int col) const
Definition: operators.h:1128
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:200
MatrixFreeOperators::Base::precondition_Jacobi
void precondition_Jacobi(VectorType &dst, const VectorType &src, const value_type omega) const
Definition: operators.h:1654
MatrixFreeOperators::Base::value_type
typename VectorType::value_type value_type
Definition: operators.h:195
MatrixFreeOperators::Base::vmult
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1326
MatrixFreeOperators::Base::vmult_interface_up
void vmult_interface_up(VectorType &dst, const VectorType &src) const
Definition: operators.h:1541
MatrixFreeOperators::Base::postprocess_constraints
void postprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1453
MatrixFreeOperators::CellwiseInverseMassMatrix::fe_eval
const FEEvaluationBase< dim, n_components, Number, false, VectorizedArrayType > & fe_eval
Definition: operators.h:720
MatrixFreeOperators::CellwiseInverseMassMatrix::fill_inverse_JxW_values
void fill_inverse_JxW_values(AlignedVector< VectorizedArrayType > &inverse_jxw) const
Definition: operators.h:984
MatrixFreeOperators::CellwiseInverseMassMatrix::CellwiseInverseMassMatrix
CellwiseInverseMassMatrix(const FEEvaluationBase< dim, n_components, Number, false, VectorizedArrayType > &fe_eval)
Definition: operators.h:962
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:1063
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:1036
MatrixFreeOperators::LaplaceOperator::get_coefficient
std::shared_ptr< Table< 2, VectorizedArrayType > > get_coefficient()
Definition: operators.h:1950
MatrixFreeOperators::LaplaceOperator::value_type
typename Base< dim, VectorType, VectorizedArrayType >::value_type value_type
Definition: operators.h:812
MatrixFreeOperators::LaplaceOperator::compute_diagonal
virtual void compute_diagonal() override
Definition: operators.h:1970
MatrixFreeOperators::LaplaceOperator::scalar_coefficient
std::shared_ptr< Table< 2, VectorizedArrayType > > scalar_coefficient
Definition: operators.h:945
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:2157
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:2057
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:2121
MatrixFreeOperators::LaplaceOperator::apply_add
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition: operators.h:2021
MatrixFreeOperators::LaplaceOperator::size_type
typename Base< dim, VectorType, VectorizedArrayType >::size_type size_type
Definition: operators.h:818
MatrixFreeOperators::LaplaceOperator::set_coefficient
void set_coefficient(const std::shared_ptr< Table< 2, VectorizedArrayType > > &scalar_coefficient)
Definition: operators.h:1930
MatrixFreeOperators::LaplaceOperator::clear
virtual void clear() override
Definition: operators.h:1909
MatrixFreeOperators::MGInterfaceOperator::value_type
typename OperatorType::value_type value_type
Definition: operators.h:541
MatrixFreeOperators::MGInterfaceOperator::mf_base_operator
SmartPointer< const OperatorType > mf_base_operator
Definition: operators.h:591
MatrixFreeOperators::MGInterfaceOperator::vmult
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1699
MatrixFreeOperators::MGInterfaceOperator::initialize
void initialize(const OperatorType &operator_in)
Definition: operators.h:1689
MatrixFreeOperators::MGInterfaceOperator::Tvmult
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1719
MatrixFreeOperators::MGInterfaceOperator::initialize_dof_vector
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1739
MatrixFreeOperators::MGInterfaceOperator::size_type
typename OperatorType::size_type size_type
Definition: operators.h:546
MatrixFreeOperators::MassOperator::size_type
typename Base< dim, VectorType, VectorizedArrayType >::size_type size_type
Definition: operators.h:752
MatrixFreeOperators::MassOperator::apply_add
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition: operators.h:1825
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:1847
MatrixFreeOperators::MassOperator::compute_diagonal
virtual void compute_diagonal() override
Definition: operators.h:1780
MatrixFreeOperators::MassOperator::value_type
typename Base< dim, VectorType, VectorizedArrayType >::value_type value_type
Definition: operators.h:746
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition: derivative_form.h:536
diagonal_matrix.h
fe_evaluation.h
level
unsigned int level
Definition: grid_out.cc:4590
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1465
Patterns::Tools::to_string
std::string to_string(const T &t)
Definition: patterns.h:2329
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
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:1575
la_parallel_vector.h
matrix_free.h
mg_constrained_dofs.h
LinearAlgebra::CUDAWrappers::kernel::size_type
types::global_dof_index size_type
Definition: cuda_kernels.h:45
MatrixFreeOperators::BlockHelper::collect_sizes
std::enable_if< IsBlockVector< VectorType >::value, void >::type collect_sizes(VectorType &vector)
Definition: operators.h:97
MatrixFreeOperators::BlockHelper::n_blocks
std::enable_if< IsBlockVector< VectorType >::value, unsignedint >::type n_blocks(const VectorType &vector)
Definition: operators.h:49
MatrixFreeOperators::BlockHelper::subblock
std::enable_if< IsBlockVector< VectorType >::value, typenameVectorType::BlockType & >::type subblock(VectorType &vector, unsigned int block_no)
Definition: operators.h:65
MatrixFreeOperators::Implementation::non_negative
bool non_negative(const VectorizedArrayType &n)
Definition: operators.h:2032
MatrixFreeTools::compute_diagonal
void compute_diagonal(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, LinearAlgebra::distributed::Vector< Number > &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)
Physics::Elasticity::Kinematics::epsilon
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
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:461
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:196
std::sqrt
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5406
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5468
Generated by doxygen 1.9.4