Reference documentation for deal.II version GIT 6bdf9a4b6c 2022-08-12 19:20: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 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii_matrix_free_operators_h
18 #define dealii_matrix_free_operators_h
19 
20 
21 #include <deal.II/base/config.h>
22 
26 
29 
33 
35 
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 =
187  class Base : public Subscriptor
188  {
189  public:
193  using value_type = typename VectorType::value_type;
194 
199 
203  Base();
204 
208  virtual ~Base() override = default;
209 
214  virtual void
215  clear();
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 
287  size_type
288  m() const;
289 
293  size_type
294  n() const;
295 
299  void
300  vmult_interface_down(VectorType &dst, const VectorType &src) const;
301 
305  void
306  vmult_interface_up(VectorType &dst, const VectorType &src) const;
307 
311  void
312  vmult(VectorType &dst, const VectorType &src) const;
313 
317  void
318  Tvmult(VectorType &dst, const VectorType &src) const;
319 
323  void
324  vmult_add(VectorType &dst, const VectorType &src) const;
325 
329  void
330  Tvmult_add(VectorType &dst, const VectorType &src) const;
331 
336  value_type
337  el(const unsigned int row, const unsigned int col) const;
338 
343  virtual std::size_t
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 
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 
615  template <int dim,
616  int fe_degree,
617  int n_components = 1,
618  typename Number = double,
619  typename VectorizedArrayType = VectorizedArray<Number>>
621  {
622  static_assert(
623  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
624  "Type of Number and of VectorizedArrayType do not match.");
625 
626  public:
632  const FEEvaluationBase<dim,
633  n_components,
634  Number,
635  false,
636  VectorizedArrayType> &fe_eval);
637 
646  void
647  apply(const AlignedVector<VectorizedArrayType> &inverse_coefficient,
648  const unsigned int n_actual_components,
649  const VectorizedArrayType * in_array,
650  VectorizedArrayType * out_array) const;
651 
663  void
664  apply(const VectorizedArrayType *in_array,
665  VectorizedArrayType * out_array) const;
666 
700  void
701  transform_from_q_points_to_basis(const unsigned int n_actual_components,
702  const VectorizedArrayType *in_array,
703  VectorizedArrayType *out_array) const;
704 
710  void
712  AlignedVector<VectorizedArrayType> &inverse_jxw) const;
713 
714  private:
718  const FEEvaluationBase<dim,
719  n_components,
720  Number,
721  false,
722  VectorizedArrayType> &fe_eval;
723  };
724 
725 
726 
734  template <int dim,
735  int fe_degree,
736  int n_q_points_1d = fe_degree + 1,
737  int n_components = 1,
738  typename VectorType = LinearAlgebra::distributed::Vector<double>,
739  typename VectorizedArrayType =
741  class MassOperator : public Base<dim, VectorType, VectorizedArrayType>
742  {
743  public:
747  using value_type =
749 
753  using size_type =
755 
759  MassOperator();
760 
764  virtual void
765  compute_diagonal() override;
766 
782  void
784 
788  const std::shared_ptr<DiagonalMatrix<VectorType>> &
790 
794  const std::shared_ptr<DiagonalMatrix<VectorType>> &
796 
797  private:
803  virtual void
804  apply_add(VectorType &dst, const VectorType &src) const override;
805 
809  void
812  VectorType & dst,
813  const VectorType & src,
814  const std::pair<unsigned int, unsigned int> & cell_range) const;
815 
820  std::shared_ptr<DiagonalMatrix<VectorType>> lumped_diagonal_entries;
821 
826  std::shared_ptr<DiagonalMatrix<VectorType>> inverse_lumped_diagonal_entries;
827  };
828 
829 
830 
841  template <int dim,
842  int fe_degree,
843  int n_q_points_1d = fe_degree + 1,
844  int n_components = 1,
845  typename VectorType = LinearAlgebra::distributed::Vector<double>,
846  typename VectorizedArrayType =
848  class LaplaceOperator : public Base<dim, VectorType, VectorizedArrayType>
849  {
850  public:
854  using value_type =
856 
860  using size_type =
862 
866  LaplaceOperator();
867 
874  virtual void
875  compute_diagonal() override;
876 
927  void
929  const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient);
930 
935  virtual void
936  clear() override;
937 
944  std::shared_ptr<Table<2, VectorizedArrayType>>
945  get_coefficient();
946 
947  private:
953  virtual void
954  apply_add(VectorType &dst, const VectorType &src) const override;
955 
959  void
962  VectorType & dst,
963  const VectorType & src,
964  const std::pair<unsigned int, unsigned int> & cell_range) const;
965 
969  void
972  VectorType & dst,
973  const VectorType &,
974  const std::pair<unsigned int, unsigned int> &cell_range) const;
975 
979  void
982  & phi,
983  const unsigned int cell) const;
984 
988  std::shared_ptr<Table<2, VectorizedArrayType>> scalar_coefficient;
989  };
990 
991 
992 
993  // ------------------------------------ inline functions ---------------------
994 
995  template <int dim,
996  int fe_degree,
997  int n_components,
998  typename Number,
999  typename VectorizedArrayType>
1000  inline CellwiseInverseMassMatrix<dim,
1001  fe_degree,
1002  n_components,
1003  Number,
1004  VectorizedArrayType>::
1005  CellwiseInverseMassMatrix(
1006  const FEEvaluationBase<dim,
1007  n_components,
1008  Number,
1009  false,
1010  VectorizedArrayType> &fe_eval)
1011  : fe_eval(fe_eval)
1012  {}
1013 
1014 
1015 
1016  template <int dim,
1017  int fe_degree,
1018  int n_components,
1019  typename Number,
1020  typename VectorizedArrayType>
1021  inline void
1023  fe_degree,
1024  n_components,
1025  Number,
1026  VectorizedArrayType>::
1027  fill_inverse_JxW_values(
1028  AlignedVector<VectorizedArrayType> &inverse_jxw) const
1029  {
1030  constexpr unsigned int dofs_per_component_on_cell =
1031  Utilities::pow(fe_degree + 1, dim);
1032  Assert(inverse_jxw.size() > 0 &&
1033  inverse_jxw.size() % dofs_per_component_on_cell == 0,
1034  ExcMessage(
1035  "Expected diagonal to be a multiple of scalar dof per cells"));
1036 
1037  // compute values for the first component
1038  for (unsigned int q = 0; q < dofs_per_component_on_cell; ++q)
1039  inverse_jxw[q] = 1. / fe_eval.JxW(q);
1040  // copy values to rest of vector
1041  for (unsigned int q = dofs_per_component_on_cell; q < inverse_jxw.size();)
1042  for (unsigned int i = 0; i < dofs_per_component_on_cell; ++i, ++q)
1043  inverse_jxw[q] = inverse_jxw[i];
1044  }
1045 
1046 
1047 
1048  template <int dim,
1049  int fe_degree,
1050  int n_components,
1051  typename Number,
1052  typename VectorizedArrayType>
1053  inline void
1055  dim,
1056  fe_degree,
1057  n_components,
1058  Number,
1059  VectorizedArrayType>::apply(const VectorizedArrayType *in_array,
1060  VectorizedArrayType * out_array) const
1061  {
1062  if (fe_degree > -1)
1064  template run<fe_degree>(n_components, fe_eval, in_array, out_array);
1065  else
1067  n_components, fe_eval, in_array, out_array);
1068  }
1069 
1070 
1071 
1072  template <int dim,
1073  int fe_degree,
1074  int n_components,
1075  typename Number,
1076  typename VectorizedArrayType>
1077  inline void
1079  fe_degree,
1080  n_components,
1081  Number,
1082  VectorizedArrayType>::
1083  apply(const AlignedVector<VectorizedArrayType> &inverse_coefficients,
1084  const unsigned int n_actual_components,
1085  const VectorizedArrayType * in_array,
1086  VectorizedArrayType * out_array) const
1087  {
1088  const unsigned int given_degree =
1089  fe_eval.get_shape_info().data[0].fe_degree;
1090  if (fe_degree > -1)
1092  VectorizedArrayType>::
1093  template run<fe_degree>(
1094  n_actual_components,
1095  fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
1096  inverse_coefficients,
1097  in_array,
1098  out_array);
1099  else
1101  n_actual_components,
1102  given_degree,
1103  fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
1104  inverse_coefficients,
1105  in_array,
1106  out_array);
1107  }
1108 
1109 
1110 
1111  template <int dim,
1112  int fe_degree,
1113  int n_components,
1114  typename Number,
1115  typename VectorizedArrayType>
1116  inline void
1118  fe_degree,
1119  n_components,
1120  Number,
1121  VectorizedArrayType>::
1122  transform_from_q_points_to_basis(const unsigned int n_actual_components,
1123  const VectorizedArrayType *in_array,
1124  VectorizedArrayType * out_array) const
1125  {
1126  const auto n_q_points_1d = fe_eval.get_shape_info().data[0].n_q_points_1d;
1127 
1128  if (fe_degree > -1 && (fe_degree + 1 == n_q_points_1d))
1130  dim,
1131  VectorizedArrayType>::template run<fe_degree,
1132  fe_degree + 1>(n_actual_components,
1133  fe_eval,
1134  in_array,
1135  out_array);
1136  else
1138  transform_from_q_points_to_basis(n_actual_components,
1139  fe_eval,
1140  in_array,
1141  out_array);
1142  }
1143 
1144 
1145 
1146  //----------------- Base operator -----------------------------
1147  template <int dim, typename VectorType, typename VectorizedArrayType>
1149  : Subscriptor()
1150  , have_interface_matrices(false)
1151  {}
1152 
1153 
1154 
1155  template <int dim, typename VectorType, typename VectorizedArrayType>
1158  {
1159  Assert(data.get() != nullptr, ExcNotInitialized());
1161  0;
1162  for (const unsigned int selected_row : selected_rows)
1163  total_size += data->get_vector_partitioner(selected_row)->size();
1164  return total_size;
1165  }
1166 
1167 
1168 
1169  template <int dim, typename VectorType, typename VectorizedArrayType>
1172  {
1173  Assert(data.get() != nullptr, ExcNotInitialized());
1175  0;
1176  for (const unsigned int selected_column : selected_columns)
1177  total_size += data->get_vector_partitioner(selected_column)->size();
1178  return total_size;
1179  }
1180 
1181 
1182 
1183  template <int dim, typename VectorType, typename VectorizedArrayType>
1184  void
1186  {
1187  data.reset();
1188  inverse_diagonal_entries.reset();
1189  }
1190 
1191 
1192 
1193  template <int dim, typename VectorType, typename VectorizedArrayType>
1196  const unsigned int col) const
1197  {
1198  (void)col;
1199  Assert(row == col, ExcNotImplemented());
1200  Assert(inverse_diagonal_entries.get() != nullptr &&
1201  inverse_diagonal_entries->m() > 0,
1202  ExcNotInitialized());
1203  return 1.0 / (*inverse_diagonal_entries)(row, row);
1204  }
1205 
1206 
1207 
1208  template <int dim, typename VectorType, typename VectorizedArrayType>
1209  void
1211  VectorType &vec) const
1212  {
1213  Assert(data.get() != nullptr, ExcNotInitialized());
1214  AssertDimension(BlockHelper::n_blocks(vec), selected_rows.size());
1215  for (unsigned int i = 0; i < BlockHelper::n_blocks(vec); ++i)
1216  {
1217  const unsigned int index = selected_rows[i];
1218  if (!BlockHelper::subblock(vec, index)
1219  .partitioners_are_compatible(
1220  *data->get_dof_info(index).vector_partitioner))
1221  data->initialize_dof_vector(BlockHelper::subblock(vec, index), index);
1222 
1223  Assert(BlockHelper::subblock(vec, index)
1224  .partitioners_are_globally_compatible(
1225  *data->get_dof_info(index).vector_partitioner),
1226  ExcInternalError());
1227  }
1229  }
1230 
1231 
1232 
1233  template <int dim, typename VectorType, typename VectorizedArrayType>
1234  void
1237  data_,
1238  const std::vector<unsigned int> &given_row_selection,
1239  const std::vector<unsigned int> &given_column_selection)
1240  {
1241  data = data_;
1242 
1243  selected_rows.clear();
1244  selected_columns.clear();
1245  if (given_row_selection.empty())
1246  for (unsigned int i = 0; i < data_->n_components(); ++i)
1247  selected_rows.push_back(i);
1248  else
1249  {
1250  for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1251  {
1252  AssertIndexRange(given_row_selection[i], data_->n_components());
1253  for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1254  if (j != i)
1255  Assert(given_row_selection[j] != given_row_selection[i],
1256  ExcMessage("Given row indices must be unique"));
1257 
1258  selected_rows.push_back(given_row_selection[i]);
1259  }
1260  }
1261  if (given_column_selection.size() == 0)
1262  selected_columns = selected_rows;
1263  else
1264  {
1265  for (unsigned int i = 0; i < given_column_selection.size(); ++i)
1266  {
1267  AssertIndexRange(given_column_selection[i], data_->n_components());
1268  for (unsigned int j = 0; j < given_column_selection.size(); ++j)
1269  if (j != i)
1270  Assert(given_column_selection[j] != given_column_selection[i],
1271  ExcMessage("Given column indices must be unique"));
1272 
1273  selected_columns.push_back(given_column_selection[i]);
1274  }
1275  }
1276 
1277  edge_constrained_indices.clear();
1278  edge_constrained_indices.resize(selected_rows.size());
1279  edge_constrained_values.clear();
1280  edge_constrained_values.resize(selected_rows.size());
1281  have_interface_matrices = false;
1282  }
1283 
1284 
1285 
1286  template <int dim, typename VectorType, typename VectorizedArrayType>
1287  void
1290  data_,
1291  const MGConstrainedDoFs & mg_constrained_dofs,
1292  const unsigned int level,
1293  const std::vector<unsigned int> &given_row_selection)
1294  {
1295  std::vector<MGConstrainedDoFs> mg_constrained_dofs_vector(
1296  1, mg_constrained_dofs);
1297  initialize(data_, mg_constrained_dofs_vector, level, given_row_selection);
1298  }
1299 
1300 
1301 
1302  template <int dim, typename VectorType, typename VectorizedArrayType>
1303  void
1306  data_,
1307  const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
1308  const unsigned int level,
1309  const std::vector<unsigned int> & given_row_selection)
1310  {
1312  ExcMessage("level is not set"));
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  selected_columns = selected_rows;
1333 
1334  AssertDimension(mg_constrained_dofs.size(), selected_rows.size());
1335  edge_constrained_indices.clear();
1336  edge_constrained_indices.resize(selected_rows.size());
1337  edge_constrained_values.clear();
1338  edge_constrained_values.resize(selected_rows.size());
1339 
1340  data = data_;
1341 
1342  for (unsigned int j = 0; j < selected_rows.size(); ++j)
1343  {
1344  if (data_->n_cell_batches() > 0)
1345  {
1346  AssertDimension(level, data_->get_cell_iterator(0, 0, j)->level());
1347  }
1348 
1349  // setup edge_constrained indices
1350  std::vector<types::global_dof_index> interface_indices;
1351  mg_constrained_dofs[j]
1352  .get_refinement_edge_indices(level)
1353  .fill_index_vector(interface_indices);
1354  edge_constrained_indices[j].clear();
1355  edge_constrained_indices[j].reserve(interface_indices.size());
1356  edge_constrained_values[j].resize(interface_indices.size());
1357  const IndexSet &locally_owned =
1358  data->get_dof_handler(selected_rows[j]).locally_owned_mg_dofs(level);
1359  for (const auto interface_index : interface_indices)
1360  if (locally_owned.is_element(interface_index))
1361  edge_constrained_indices[j].push_back(
1362  locally_owned.index_within_set(interface_index));
1363  have_interface_matrices |=
1365  static_cast<unsigned int>(edge_constrained_indices[j].size()),
1366  data->get_vector_partitioner()->get_mpi_communicator()) > 0;
1367  }
1368  }
1369 
1370 
1371 
1372  template <int dim, typename VectorType, typename VectorizedArrayType>
1373  void
1375  VectorType &dst) const
1376  {
1377  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1378  {
1379  const std::vector<unsigned int> &constrained_dofs =
1380  data->get_constrained_dofs(selected_rows[j]);
1381  for (const auto constrained_dof : constrained_dofs)
1382  BlockHelper::subblock(dst, j).local_element(constrained_dof) = 1.;
1383  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1384  BlockHelper::subblock(dst, j).local_element(
1385  edge_constrained_indices[j][i]) = 1.;
1386  }
1387  }
1388 
1389 
1390 
1391  template <int dim, typename VectorType, typename VectorizedArrayType>
1392  void
1394  const VectorType &src) const
1395  {
1396  using Number =
1398  dst = Number(0.);
1399  vmult_add(dst, src);
1400  }
1401 
1402 
1403 
1404  template <int dim, typename VectorType, typename VectorizedArrayType>
1405  void
1407  VectorType & dst,
1408  const VectorType &src) const
1409  {
1410  mult_add(dst, src, false);
1411  }
1412 
1413 
1414 
1415  template <int dim, typename VectorType, typename VectorizedArrayType>
1416  void
1418  VectorType & dst,
1419  const VectorType &src) const
1420  {
1421  mult_add(dst, src, true);
1422  }
1423 
1424 
1425 
1426  template <int dim, typename VectorType, typename VectorizedArrayType>
1427  void
1429  const VectorType &src,
1430  const bool is_row) const
1431  {
1432  using Number =
1434  for (unsigned int i = 0; i < BlockHelper::n_blocks(src); ++i)
1435  {
1436  const unsigned int mf_component =
1437  is_row ? selected_rows[i] : selected_columns[i];
1438  // If both vectors use the same partitioner -> done
1439  if (BlockHelper::subblock(src, i).get_partitioner().get() ==
1440  data->get_dof_info(mf_component).vector_partitioner.get())
1441  continue;
1442 
1443  // If not, assert that the local ranges are the same and reset to the
1444  // current partitioner
1446  .get_partitioner()
1447  ->locally_owned_size() ==
1448  data->get_dof_info(mf_component)
1449  .vector_partitioner->locally_owned_size(),
1450  ExcMessage(
1451  "The vector passed to the vmult() function does not have "
1452  "the correct size for compatibility with MatrixFree."));
1453 
1454  // copy the vector content to a temporary vector so that it does not get
1455  // lost
1457  BlockHelper::subblock(src, i));
1458  this->data->initialize_dof_vector(
1459  BlockHelper::subblock(const_cast<VectorType &>(src), i),
1460  mf_component);
1461  BlockHelper::subblock(const_cast<VectorType &>(src), i)
1462  .copy_locally_owned_data_from(copy_vec);
1463  }
1464  }
1465 
1466 
1467 
1468  template <int dim, typename VectorType, typename VectorizedArrayType>
1469  void
1471  VectorType & dst,
1472  const VectorType &src) const
1473  {
1474  using Number =
1476  adjust_ghost_range_if_necessary(src, false);
1477  adjust_ghost_range_if_necessary(dst, true);
1478 
1479  // set zero Dirichlet values on the input vector (and remember the src and
1480  // dst values because we need to reset them at the end)
1481  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1482  {
1483  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1484  {
1485  edge_constrained_values[j][i] = std::pair<Number, Number>(
1486  BlockHelper::subblock(src, j).local_element(
1487  edge_constrained_indices[j][i]),
1488  BlockHelper::subblock(dst, j).local_element(
1489  edge_constrained_indices[j][i]));
1490  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1491  .local_element(edge_constrained_indices[j][i]) = 0.;
1492  }
1493  }
1494  }
1495 
1496 
1497 
1498  template <int dim, typename VectorType, typename VectorizedArrayType>
1499  void
1501  VectorType & dst,
1502  const VectorType &src,
1503  const bool transpose) const
1504  {
1505  AssertDimension(dst.size(), src.size());
1507  AssertDimension(BlockHelper::n_blocks(dst), selected_rows.size());
1508  preprocess_constraints(dst, src);
1509  if (transpose)
1510  Tapply_add(dst, src);
1511  else
1512  apply_add(dst, src);
1513  postprocess_constraints(dst, src);
1514  }
1515 
1516 
1517 
1518  template <int dim, typename VectorType, typename VectorizedArrayType>
1519  void
1521  VectorType & dst,
1522  const VectorType &src) const
1523  {
1524  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1525  {
1526  const std::vector<unsigned int> &constrained_dofs =
1527  data->get_constrained_dofs(selected_rows[j]);
1528  for (const auto constrained_dof : constrained_dofs)
1529  BlockHelper::subblock(dst, j).local_element(constrained_dof) +=
1530  BlockHelper::subblock(src, j).local_element(constrained_dof);
1531  }
1532 
1533  // reset edge constrained values, multiply by unit matrix and add into
1534  // destination
1535  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1536  {
1537  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1538  {
1539  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1540  .local_element(edge_constrained_indices[j][i]) =
1541  edge_constrained_values[j][i].first;
1542  BlockHelper::subblock(dst, j).local_element(
1543  edge_constrained_indices[j][i]) =
1544  edge_constrained_values[j][i].second +
1545  edge_constrained_values[j][i].first;
1546  }
1547  }
1548  }
1549 
1550 
1551 
1552  template <int dim, typename VectorType, typename VectorizedArrayType>
1553  void
1555  VectorType & dst,
1556  const VectorType &src) const
1557  {
1558  using Number =
1560  AssertDimension(dst.size(), src.size());
1561  adjust_ghost_range_if_necessary(src, false);
1562  adjust_ghost_range_if_necessary(dst, true);
1563 
1564  dst = Number(0.);
1565 
1566  if (!have_interface_matrices)
1567  return;
1568 
1569  // set zero Dirichlet values on the input vector (and remember the src and
1570  // dst values because we need to reset them at the end)
1571  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1572  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1573  {
1574  edge_constrained_values[j][i] = std::pair<Number, Number>(
1575  BlockHelper::subblock(src, j).local_element(
1576  edge_constrained_indices[j][i]),
1577  BlockHelper::subblock(dst, j).local_element(
1578  edge_constrained_indices[j][i]));
1579  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1580  .local_element(edge_constrained_indices[j][i]) = 0.;
1581  }
1582 
1583  apply_add(dst, src);
1584 
1585  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1586  {
1587  unsigned int c = 0;
1588  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1589  {
1590  for (; c < edge_constrained_indices[j][i]; ++c)
1591  BlockHelper::subblock(dst, j).local_element(c) = 0.;
1592  ++c;
1593 
1594  // reset the src values
1595  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1596  .local_element(edge_constrained_indices[j][i]) =
1597  edge_constrained_values[j][i].first;
1598  }
1599  for (; c < BlockHelper::subblock(dst, j).locally_owned_size(); ++c)
1600  BlockHelper::subblock(dst, j).local_element(c) = 0.;
1601  }
1602  }
1603 
1604 
1605 
1606  template <int dim, typename VectorType, typename VectorizedArrayType>
1607  void
1609  VectorType & dst,
1610  const VectorType &src) const
1611  {
1612  using Number =
1614  AssertDimension(dst.size(), src.size());
1615  adjust_ghost_range_if_necessary(src, false);
1616  adjust_ghost_range_if_necessary(dst, true);
1617 
1618  dst = Number(0.);
1619 
1620  if (!have_interface_matrices)
1621  return;
1622 
1623  VectorType src_cpy(src);
1624  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1625  {
1626  unsigned int c = 0;
1627  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1628  {
1629  for (; c < edge_constrained_indices[j][i]; ++c)
1630  BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1631  ++c;
1632  }
1633  for (; c < BlockHelper::subblock(src_cpy, j).locally_owned_size(); ++c)
1634  BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1635  }
1636 
1637  apply_add(dst, src_cpy);
1638 
1639  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1640  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1641  BlockHelper::subblock(dst, j).local_element(
1642  edge_constrained_indices[j][i]) = 0.;
1643  }
1644 
1645 
1646 
1647  template <int dim, typename VectorType, typename VectorizedArrayType>
1648  void
1650  VectorType & dst,
1651  const VectorType &src) const
1652  {
1653  using Number =
1655  dst = Number(0.);
1656  Tvmult_add(dst, src);
1657  }
1658 
1659 
1660 
1661  template <int dim, typename VectorType, typename VectorizedArrayType>
1662  std::size_t
1664  {
1665  return inverse_diagonal_entries.get() != nullptr ?
1666  inverse_diagonal_entries->memory_consumption() :
1667  sizeof(*this);
1668  }
1669 
1670 
1671 
1672  template <int dim, typename VectorType, typename VectorizedArrayType>
1673  std::shared_ptr<const MatrixFree<
1674  dim,
1676  VectorizedArrayType>>
1678  {
1679  return data;
1680  }
1681 
1682 
1683 
1684  template <int dim, typename VectorType, typename VectorizedArrayType>
1685  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1687  const
1688  {
1689  Assert(inverse_diagonal_entries.get() != nullptr &&
1690  inverse_diagonal_entries->m() > 0,
1691  ExcNotInitialized());
1692  return inverse_diagonal_entries;
1693  }
1694 
1695 
1696 
1697  template <int dim, typename VectorType, typename VectorizedArrayType>
1698  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1700  {
1701  Assert(diagonal_entries.get() != nullptr && diagonal_entries->m() > 0,
1702  ExcNotInitialized());
1703  return diagonal_entries;
1704  }
1705 
1706 
1707 
1708  template <int dim, typename VectorType, typename VectorizedArrayType>
1709  void
1711  VectorType & dst,
1712  const VectorType &src) const
1713  {
1714  apply_add(dst, src);
1715  }
1716 
1717 
1718 
1719  template <int dim, typename VectorType, typename VectorizedArrayType>
1720  void
1722  VectorType & dst,
1723  const VectorType & src,
1725  const
1726  {
1727  Assert(inverse_diagonal_entries.get() && inverse_diagonal_entries->m() > 0,
1728  ExcNotInitialized());
1729  inverse_diagonal_entries->vmult(dst, src);
1730  dst *= omega;
1731  }
1732 
1733 
1734 
1735  //------------------------- MGInterfaceOperator ------------------------------
1736 
1737  template <typename OperatorType>
1739  : Subscriptor()
1740  , mf_base_operator(nullptr)
1741  {}
1742 
1743 
1744 
1745  template <typename OperatorType>
1746  void
1748  {
1749  mf_base_operator = nullptr;
1750  }
1751 
1752 
1753 
1754  template <typename OperatorType>
1755  void
1756  MGInterfaceOperator<OperatorType>::initialize(const OperatorType &operator_in)
1757  {
1758  mf_base_operator = &operator_in;
1759  }
1760 
1761 
1762 
1763  template <typename OperatorType>
1764  template <typename VectorType>
1765  void
1767  const VectorType &src) const
1768  {
1769 #ifndef DEAL_II_MSVC
1770  static_assert(
1771  std::is_same<typename VectorType::value_type, value_type>::value,
1772  "The vector type must be based on the same value type as this "
1773  "operator");
1774 #endif
1775 
1776  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1777 
1778  mf_base_operator->vmult_interface_down(dst, src);
1779  }
1780 
1781 
1782 
1783  template <typename OperatorType>
1784  template <typename VectorType>
1785  void
1787  const VectorType &src) const
1788  {
1789 #ifndef DEAL_II_MSVC
1790  static_assert(
1791  std::is_same<typename VectorType::value_type, value_type>::value,
1792  "The vector type must be based on the same value type as this "
1793  "operator");
1794 #endif
1795 
1796  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1797 
1798  mf_base_operator->vmult_interface_up(dst, src);
1799  }
1800 
1801 
1802 
1803  template <typename OperatorType>
1804  template <typename VectorType>
1805  void
1807  VectorType &vec) const
1808  {
1809  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1810 
1811  mf_base_operator->initialize_dof_vector(vec);
1812  }
1813 
1814 
1815 
1816  //-----------------------------MassOperator----------------------------------
1817 
1818  template <int dim,
1819  int fe_degree,
1820  int n_q_points_1d,
1821  int n_components,
1822  typename VectorType,
1823  typename VectorizedArrayType>
1824  MassOperator<dim,
1825  fe_degree,
1826  n_q_points_1d,
1827  n_components,
1828  VectorType,
1829  VectorizedArrayType>::MassOperator()
1830  : Base<dim, VectorType, VectorizedArrayType>()
1831  {
1832  AssertThrow(
1835  "This class only supports the non-blocked vector variant of the Base "
1836  "operator because only a single FEEvaluation object is used in the "
1837  "apply function."));
1838  }
1839 
1840 
1841 
1842  template <int dim,
1843  int fe_degree,
1844  int n_q_points_1d,
1845  int n_components,
1846  typename VectorType,
1847  typename VectorizedArrayType>
1848  void
1849  MassOperator<dim,
1850  fe_degree,
1851  n_q_points_1d,
1852  n_components,
1853  VectorType,
1854  VectorizedArrayType>::compute_diagonal()
1855  {
1856  using Number =
1859  ExcNotInitialized());
1860  Assert(this->selected_rows == this->selected_columns,
1861  ExcMessage("This function is only implemented for square (not "
1862  "rectangular) operators."));
1863 
1864  this->inverse_diagonal_entries =
1865  std::make_shared<DiagonalMatrix<VectorType>>();
1866  this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1867  VectorType &inverse_diagonal_vector =
1868  this->inverse_diagonal_entries->get_vector();
1869  VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1870  this->initialize_dof_vector(inverse_diagonal_vector);
1871  this->initialize_dof_vector(diagonal_vector);
1872 
1873  // Set up the action of the mass matrix in a way that's compatible with
1874  // MatrixFreeTools::compute_diagonal:
1875  auto diagonal_evaluation = [](auto &integrator) {
1876  integrator.evaluate(EvaluationFlags::values);
1877  for (unsigned int q = 0; q < integrator.n_q_points; ++q)
1878  integrator.submit_value(integrator.get_value(q), q);
1879  integrator.integrate(EvaluationFlags::values);
1880  };
1881 
1882  std::function<void(
1883  FEEvaluation<
1884  dim,
1885  fe_degree,
1886  n_q_points_1d,
1887  n_components,
1889  VectorizedArrayType> &)>
1890  diagonal_evaluation_f(diagonal_evaluation);
1891 
1892  Assert(this->selected_rows.size() > 0, ExcInternalError());
1893  for (unsigned int block_n = 0; block_n < this->selected_rows.size();
1894  ++block_n)
1896  BlockHelper::subblock(diagonal_vector,
1897  block_n),
1898  diagonal_evaluation_f,
1899  this->selected_rows[block_n]);
1900 
1901  // Constrained entries will create zeros on the main diagonal, which we
1902  // don't want
1903  this->set_constrained_entries_to_one(diagonal_vector);
1904 
1905  inverse_diagonal_vector = diagonal_vector;
1906 
1907  for (unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
1908  ++i)
1909  {
1910  Assert(diagonal_vector.local_element(i) > Number(0),
1911  ExcInternalError());
1912  inverse_diagonal_vector.local_element(i) =
1913  1. / inverse_diagonal_vector.local_element(i);
1914  }
1915 
1916  // We never need ghost values so don't update them
1917  }
1918 
1919 
1920 
1921  template <int dim,
1922  int fe_degree,
1923  int n_q_points_1d,
1924  int n_components,
1925  typename VectorType,
1926  typename VectorizedArrayType>
1927  void
1928  MassOperator<dim,
1929  fe_degree,
1930  n_q_points_1d,
1931  n_components,
1932  VectorType,
1933  VectorizedArrayType>::compute_lumped_diagonal()
1934  {
1935  using Number =
1938  ExcNotInitialized());
1939  Assert(this->selected_rows == this->selected_columns,
1940  ExcMessage("This function is only implemented for square (not "
1941  "rectangular) operators."));
1942 
1943  inverse_lumped_diagonal_entries =
1944  std::make_shared<DiagonalMatrix<VectorType>>();
1945  lumped_diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1946  VectorType &inverse_lumped_diagonal_vector =
1947  inverse_lumped_diagonal_entries->get_vector();
1948  VectorType &lumped_diagonal_vector = lumped_diagonal_entries->get_vector();
1949  this->initialize_dof_vector(inverse_lumped_diagonal_vector);
1950  this->initialize_dof_vector(lumped_diagonal_vector);
1951 
1952  // Re-use the inverse_lumped_diagonal_vector as the vector of 1s
1953  inverse_lumped_diagonal_vector = Number(1.);
1954  apply_add(lumped_diagonal_vector, inverse_lumped_diagonal_vector);
1955  this->set_constrained_entries_to_one(lumped_diagonal_vector);
1956 
1957  const size_type locally_owned_size =
1958  inverse_lumped_diagonal_vector.locally_owned_size();
1959  // A caller may request a lumped diagonal matrix when it doesn't make sense
1960  // (e.g., an element with negative-mean basis functions). Avoid division by
1961  // zero so we don't cause a floating point exception but permit negative
1962  // entries here.
1963  for (size_type i = 0; i < locally_owned_size; ++i)
1964  {
1965  if (lumped_diagonal_vector.local_element(i) == Number(0.))
1966  inverse_lumped_diagonal_vector.local_element(i) = Number(1.);
1967  else
1968  inverse_lumped_diagonal_vector.local_element(i) =
1969  Number(1.) / lumped_diagonal_vector.local_element(i);
1970  }
1971 
1972  inverse_lumped_diagonal_vector.update_ghost_values();
1973  lumped_diagonal_vector.update_ghost_values();
1974  }
1975 
1976 
1977 
1978  template <int dim,
1979  int fe_degree,
1980  int n_q_points_1d,
1981  int n_components,
1982  typename VectorType,
1983  typename VectorizedArrayType>
1984  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1985  MassOperator<dim,
1986  fe_degree,
1987  n_q_points_1d,
1988  n_components,
1989  VectorType,
1990  VectorizedArrayType>::get_matrix_lumped_diagonal_inverse() const
1991  {
1992  Assert(inverse_lumped_diagonal_entries.get() != nullptr &&
1993  inverse_lumped_diagonal_entries->m() > 0,
1994  ExcNotInitialized());
1995  return inverse_lumped_diagonal_entries;
1996  }
1997 
1998 
1999 
2000  template <int dim,
2001  int fe_degree,
2002  int n_q_points_1d,
2003  int n_components,
2004  typename VectorType,
2005  typename VectorizedArrayType>
2006  const std::shared_ptr<DiagonalMatrix<VectorType>> &
2007  MassOperator<dim,
2008  fe_degree,
2009  n_q_points_1d,
2010  n_components,
2011  VectorType,
2012  VectorizedArrayType>::get_matrix_lumped_diagonal() const
2013  {
2014  Assert(lumped_diagonal_entries.get() != nullptr &&
2015  lumped_diagonal_entries->m() > 0,
2016  ExcNotInitialized());
2017  return lumped_diagonal_entries;
2018  }
2019 
2020 
2021 
2022  template <int dim,
2023  int fe_degree,
2024  int n_q_points_1d,
2025  int n_components,
2026  typename VectorType,
2027  typename VectorizedArrayType>
2028  void
2029  MassOperator<dim,
2030  fe_degree,
2031  n_q_points_1d,
2032  n_components,
2033  VectorType,
2034  VectorizedArrayType>::apply_add(VectorType & dst,
2035  const VectorType &src) const
2036  {
2038  &MassOperator::local_apply_cell, this, dst, src);
2039  }
2040 
2041 
2042 
2043  template <int dim,
2044  int fe_degree,
2045  int n_q_points_1d,
2046  int n_components,
2047  typename VectorType,
2048  typename VectorizedArrayType>
2049  void
2050  MassOperator<dim,
2051  fe_degree,
2052  n_q_points_1d,
2053  n_components,
2054  VectorType,
2055  VectorizedArrayType>::
2056  local_apply_cell(
2057  const MatrixFree<
2058  dim,
2060  VectorizedArrayType> & data,
2061  VectorType & dst,
2062  const VectorType & src,
2063  const std::pair<unsigned int, unsigned int> &cell_range) const
2064  {
2065  using Number =
2067  FEEvaluation<dim,
2068  fe_degree,
2069  n_q_points_1d,
2070  n_components,
2071  Number,
2072  VectorizedArrayType>
2073  phi(data, this->selected_rows[0]);
2074  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2075  {
2076  phi.reinit(cell);
2077  phi.read_dof_values(src);
2078  phi.evaluate(EvaluationFlags::values);
2079  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2080  phi.submit_value(phi.get_value(q), q);
2081  phi.integrate(EvaluationFlags::values);
2082  phi.distribute_local_to_global(dst);
2083  }
2084  }
2085 
2086 
2087  //-----------------------------LaplaceOperator----------------------------------
2088 
2089  template <int dim,
2090  int fe_degree,
2091  int n_q_points_1d,
2092  int n_components,
2093  typename VectorType,
2094  typename VectorizedArrayType>
2095  LaplaceOperator<dim,
2096  fe_degree,
2097  n_q_points_1d,
2098  n_components,
2099  VectorType,
2100  VectorizedArrayType>::LaplaceOperator()
2101  : Base<dim, VectorType, VectorizedArrayType>()
2102  {}
2103 
2104 
2105 
2106  template <int dim,
2107  int fe_degree,
2108  int n_q_points_1d,
2109  int n_components,
2110  typename VectorType,
2111  typename VectorizedArrayType>
2112  void
2113  LaplaceOperator<dim,
2114  fe_degree,
2115  n_q_points_1d,
2116  n_components,
2117  VectorType,
2118  VectorizedArrayType>::clear()
2119  {
2121  scalar_coefficient.reset();
2122  }
2123 
2124 
2125 
2126  template <int dim,
2127  int fe_degree,
2128  int n_q_points_1d,
2129  int n_components,
2130  typename VectorType,
2131  typename VectorizedArrayType>
2132  void
2133  LaplaceOperator<dim,
2134  fe_degree,
2135  n_q_points_1d,
2136  n_components,
2137  VectorType,
2138  VectorizedArrayType>::
2139  set_coefficient(
2140  const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient_)
2141  {
2142  scalar_coefficient = scalar_coefficient_;
2143  }
2144 
2145 
2146 
2147  template <int dim,
2148  int fe_degree,
2149  int n_q_points_1d,
2150  int n_components,
2151  typename VectorType,
2152  typename VectorizedArrayType>
2153  std::shared_ptr<Table<2, VectorizedArrayType>>
2154  LaplaceOperator<dim,
2155  fe_degree,
2156  n_q_points_1d,
2157  n_components,
2158  VectorType,
2159  VectorizedArrayType>::get_coefficient()
2160  {
2161  Assert(scalar_coefficient.get(), ExcNotInitialized());
2162  return scalar_coefficient;
2163  }
2164 
2165 
2166 
2167  template <int dim,
2168  int fe_degree,
2169  int n_q_points_1d,
2170  int n_components,
2171  typename VectorType,
2172  typename VectorizedArrayType>
2173  void
2174  LaplaceOperator<dim,
2175  fe_degree,
2176  n_q_points_1d,
2177  n_components,
2178  VectorType,
2179  VectorizedArrayType>::compute_diagonal()
2180  {
2181  using Number =
2184  ExcNotInitialized());
2185 
2186  this->inverse_diagonal_entries =
2187  std::make_shared<DiagonalMatrix<VectorType>>();
2188  this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
2189  VectorType &inverse_diagonal_vector =
2190  this->inverse_diagonal_entries->get_vector();
2191  VectorType &diagonal_vector = this->diagonal_entries->get_vector();
2192  this->initialize_dof_vector(inverse_diagonal_vector);
2193  this->initialize_dof_vector(diagonal_vector);
2194 
2195  this->data->cell_loop(&LaplaceOperator::local_diagonal_cell,
2196  this,
2197  diagonal_vector,
2198  /*unused*/ diagonal_vector);
2199  this->set_constrained_entries_to_one(diagonal_vector);
2200 
2201  inverse_diagonal_vector = diagonal_vector;
2202 
2203  for (unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
2204  ++i)
2205  if (std::abs(inverse_diagonal_vector.local_element(i)) >
2207  inverse_diagonal_vector.local_element(i) =
2208  1. / inverse_diagonal_vector.local_element(i);
2209  else
2210  inverse_diagonal_vector.local_element(i) = 1.;
2211 
2212  // We never need ghost values so don't update them
2213  }
2214 
2215 
2216 
2217  template <int dim,
2218  int fe_degree,
2219  int n_q_points_1d,
2220  int n_components,
2221  typename VectorType,
2222  typename VectorizedArrayType>
2223  void
2224  LaplaceOperator<dim,
2225  fe_degree,
2226  n_q_points_1d,
2227  n_components,
2228  VectorType,
2229  VectorizedArrayType>::apply_add(VectorType & dst,
2230  const VectorType &src) const
2231  {
2233  &LaplaceOperator::local_apply_cell, this, dst, src);
2234  }
2235 
2236  namespace Implementation
2237  {
2238  template <typename VectorizedArrayType>
2239  bool
2240  non_negative(const VectorizedArrayType &n)
2241  {
2242  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2243  if (n[v] < 0.)
2244  return false;
2245 
2246  return true;
2247  }
2248  } // namespace Implementation
2249 
2250 
2251 
2252  template <int dim,
2253  int fe_degree,
2254  int n_q_points_1d,
2255  int n_components,
2256  typename VectorType,
2257  typename VectorizedArrayType>
2258  void
2259  LaplaceOperator<dim,
2260  fe_degree,
2261  n_q_points_1d,
2262  n_components,
2263  VectorType,
2264  VectorizedArrayType>::
2265  do_operation_on_cell(
2266  FEEvaluation<
2267  dim,
2268  fe_degree,
2269  n_q_points_1d,
2270  n_components,
2272  const unsigned int cell) const
2273  {
2274  phi.evaluate(EvaluationFlags::gradients);
2275  if (scalar_coefficient.get())
2276  {
2277  Assert(scalar_coefficient->size(1) == 1 ||
2278  scalar_coefficient->size(1) == phi.n_q_points,
2279  ExcMessage("The number of columns in the coefficient table must "
2280  "be either 1 or the number of quadrature points " +
2281  std::to_string(phi.n_q_points) +
2282  ", but the given value was " +
2283  std::to_string(scalar_coefficient->size(1))));
2284  if (scalar_coefficient->size(1) == phi.n_q_points)
2285  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2286  {
2288  (*scalar_coefficient)(cell, q)),
2289  ExcMessage("Coefficient must be non-negative"));
2290  phi.submit_gradient((*scalar_coefficient)(cell, q) *
2291  phi.get_gradient(q),
2292  q);
2293  }
2294  else
2295  {
2296  Assert(Implementation::non_negative((*scalar_coefficient)(cell, 0)),
2297  ExcMessage("Coefficient must be non-negative"));
2298  const VectorizedArrayType coefficient =
2299  (*scalar_coefficient)(cell, 0);
2300  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2301  phi.submit_gradient(coefficient * phi.get_gradient(q), q);
2302  }
2303  }
2304  else
2305  {
2306  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2307  {
2308  phi.submit_gradient(phi.get_gradient(q), q);
2309  }
2310  }
2311  phi.integrate(EvaluationFlags::gradients);
2312  }
2313 
2314 
2315 
2316  template <int dim,
2317  int fe_degree,
2318  int n_q_points_1d,
2319  int n_components,
2320  typename VectorType,
2321  typename VectorizedArrayType>
2322  void
2323  LaplaceOperator<dim,
2324  fe_degree,
2325  n_q_points_1d,
2326  n_components,
2327  VectorType,
2328  VectorizedArrayType>::
2329  local_apply_cell(
2330  const MatrixFree<
2331  dim,
2333  VectorizedArrayType> & data,
2334  VectorType & dst,
2335  const VectorType & src,
2336  const std::pair<unsigned int, unsigned int> &cell_range) const
2337  {
2338  using Number =
2341  data, this->selected_rows[0]);
2342  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2343  {
2344  phi.reinit(cell);
2345  phi.read_dof_values(src);
2346  do_operation_on_cell(phi, cell);
2347  phi.distribute_local_to_global(dst);
2348  }
2349  }
2350 
2351 
2352  template <int dim,
2353  int fe_degree,
2354  int n_q_points_1d,
2355  int n_components,
2356  typename VectorType,
2357  typename VectorizedArrayType>
2358  void
2359  LaplaceOperator<dim,
2360  fe_degree,
2361  n_q_points_1d,
2362  n_components,
2363  VectorType,
2364  VectorizedArrayType>::
2365  local_diagonal_cell(
2366  const MatrixFree<
2367  dim,
2369  VectorizedArrayType> &data,
2370  VectorType & dst,
2371  const VectorType &,
2372  const std::pair<unsigned int, unsigned int> &cell_range) const
2373  {
2374  using Number =
2376 
2378  data, this->selected_rows[0]);
2379  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2380  {
2381  phi.reinit(cell);
2382  VectorizedArrayType local_diagonal_vector[phi.static_dofs_per_cell];
2383  for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2384  {
2385  for (unsigned int j = 0; j < phi.dofs_per_component; ++j)
2386  phi.begin_dof_values()[j] = VectorizedArrayType();
2387  phi.begin_dof_values()[i] = 1.;
2388  do_operation_on_cell(phi, cell);
2389  local_diagonal_vector[i] = phi.begin_dof_values()[i];
2390  }
2391  for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2392  for (unsigned int c = 0; c < phi.n_components; ++c)
2393  phi.begin_dof_values()[i + c * phi.dofs_per_component] =
2394  local_diagonal_vector[i];
2395  phi.distribute_local_to_global(dst);
2396  }
2397  }
2398 
2399 
2400 } // end of namespace MatrixFreeOperators
2401 
2402 
2404 
2405 #endif
size_type size() const
void distribute_local_to_global(VectorType &dst, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip()) const
void read_dof_values(const VectorType &src, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip())
const Number * begin_dof_values() const
const unsigned int dofs_per_component
void reinit(const unsigned int cell_batch_index)
static constexpr unsigned int static_dofs_per_cell
static constexpr unsigned int n_components
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1850
bool is_element(const size_type index) const
Definition: index_set.h:1743
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:1417
void set_constrained_entries_to_one(VectorType &dst) const
Definition: operators.h:1374
virtual void compute_diagonal()=0
void vmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1406
void vmult_interface_down(VectorType &dst, const VectorType &src) const
Definition: operators.h:1554
size_type n() const
Definition: operators.h:1171
void preprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1470
void mult_add(VectorType &dst, const VectorType &src, const bool transpose) const
Definition: operators.h:1500
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:1699
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1649
std::vector< std::vector< std::pair< value_type, value_type > > > edge_constrained_values
Definition: operators.h:471
size_type m() const
Definition: operators.h:1157
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:1677
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal_inverse() const
Definition: operators.h:1686
virtual void clear()
Definition: operators.h:1185
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1210
std::shared_ptr< DiagonalMatrix< VectorType > > diagonal_entries
Definition: operators.h:441
virtual void Tapply_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1710
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:1288
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:1428
virtual std::size_t memory_consumption() const
Definition: operators.h:1663
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:1304
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:1195
virtual void apply_add(VectorType &dst, const VectorType &src) const =0
typename VectorType::size_type size_type
Definition: operators.h:198
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:1235
void precondition_Jacobi(VectorType &dst, const VectorType &src, const value_type omega) const
Definition: operators.h:1721
typename VectorType::value_type value_type
Definition: operators.h:193
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1393
void vmult_interface_up(VectorType &dst, const VectorType &src) const
Definition: operators.h:1608
void postprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1520
const FEEvaluationBase< dim, n_components, Number, false, VectorizedArrayType > & fe_eval
Definition: operators.h:722
void fill_inverse_JxW_values(AlignedVector< VectorizedArrayType > &inverse_jxw) const
Definition: operators.h:1027
CellwiseInverseMassMatrix(const FEEvaluationBase< dim, n_components, Number, false, VectorizedArrayType > &fe_eval)
Definition: operators.h:1005
void transform_from_q_points_to_basis(const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition: operators.h:1122
void apply(const AlignedVector< VectorizedArrayType > &inverse_coefficient, const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition: operators.h:1083
std::shared_ptr< Table< 2, VectorizedArrayType > > get_coefficient()
Definition: operators.h:2159
typename Base< dim, VectorType, VectorizedArrayType >::value_type value_type
Definition: operators.h:855
virtual void compute_diagonal() override
Definition: operators.h:2179
std::shared_ptr< Table< 2, VectorizedArrayType > > scalar_coefficient
Definition: operators.h:988
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:2365
void do_operation_on_cell(FEEvaluation< dim, fe_degree, n_q_points_1d, n_components, value_type > &phi, const unsigned int cell) const
Definition: operators.h:2265
void set_coefficient(const std::shared_ptr< Table< 2, VectorizedArrayType >> &scalar_coefficient)
Definition: operators.h:2139
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:2329
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition: operators.h:2229
typename Base< dim, VectorType, VectorizedArrayType >::size_type size_type
Definition: operators.h:861
virtual void clear() override
Definition: operators.h:2118
typename OperatorType::value_type value_type
Definition: operators.h:543
SmartPointer< const OperatorType > mf_base_operator
Definition: operators.h:593
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1766
void initialize(const OperatorType &operator_in)
Definition: operators.h:1756
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1786
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1806
typename OperatorType::size_type size_type
Definition: operators.h:548
std::shared_ptr< DiagonalMatrix< VectorType > > lumped_diagonal_entries
Definition: operators.h:820
typename Base< dim, VectorType, VectorizedArrayType >::size_type size_type
Definition: operators.h:754
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition: operators.h:2034
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_lumped_diagonal() const
Definition: operators.h:2012
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_lumped_diagonal_inverse() const
Definition: operators.h:1990
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:2056
virtual void compute_diagonal() override
Definition: operators.h:1854
typename Base< dim, VectorType, VectorizedArrayType >::value_type value_type
Definition: operators.h:748
std::shared_ptr< DiagonalMatrix< VectorType > > inverse_lumped_diagonal_entries
Definition: operators.h:826
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
unsigned int level
Definition: grid_out.cc:4607
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
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:49
std::enable_if_t< IsBlockVector< VectorType >::value, void > collect_sizes(VectorType &vector)
Definition: operators.h:95
std::enable_if_t< IsBlockVector< VectorType >::value, typename VectorType::BlockType & > subblock(VectorType &vector, unsigned int block_no)
Definition: operators.h:64
bool non_negative(const VectorizedArrayType &n)
Definition: operators.h:2240
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:2393
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:462
void run(const Iterator &begin, const typename identity< Iterator >::type &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:474
static const unsigned int invalid_unsigned_int
Definition: types.h:201
auto apply(F &&fn, Tuple &&t) -> decltype(apply_impl(std::forward< F >(fn), std::forward< Tuple >(t), std::make_index_sequence< std::tuple_size< typename std::remove_reference< Tuple >::type >::value >()))
Definition: tuple.h:36
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)