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