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