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