Reference documentation for deal.II version Git cf38156a17 2020-04-01 20:09:03 -0400
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
operators.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2011 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii_matrix_free_operators_h
18 #define dealii_matrix_free_operators_h
19 
20 
21 #include <deal.II/base/config.h>
22 
23 #include <deal.II/base/exceptions.h>
24 #include <deal.II/base/subscriptor.h>
25 #include <deal.II/base/vectorization.h>
26 
27 #include <deal.II/lac/diagonal_matrix.h>
28 #include <deal.II/lac/la_parallel_vector.h>
29 
30 #include <deal.II/matrix_free/fe_evaluation.h>
31 #include <deal.II/matrix_free/matrix_free.h>
32 
33 #include <deal.II/multigrid/mg_constrained_dofs.h>
34 
35 
36 DEAL_II_NAMESPACE_OPEN
37 
38 
39 namespace MatrixFreeOperators
40 {
41  namespace BlockHelper
42  {
43  // workaroud for unifying non-block vector and block vector implementations
44  // a non-block vector has one block and the only subblock is the vector
45  // itself
46  template <typename VectorType>
47  typename std::enable_if<IsBlockVector<VectorType>::value,
48  unsigned int>::type
49  n_blocks(const VectorType &vector)
50  {
51  return vector.n_blocks();
52  }
53 
54  template <typename VectorType>
55  typename std::enable_if<!IsBlockVector<VectorType>::value,
56  unsigned int>::type
57  n_blocks(const VectorType &)
58  {
59  return 1;
60  }
61 
62  template <typename VectorType>
63  typename std::enable_if<IsBlockVector<VectorType>::value,
64  typename VectorType::BlockType &>::type
65  subblock(VectorType &vector, unsigned int block_no)
66  {
67  return vector.block(block_no);
68  }
69 
70  template <typename VectorType>
71  typename std::enable_if<IsBlockVector<VectorType>::value,
72  const typename VectorType::BlockType &>::type
73  subblock(const VectorType &vector, unsigned int block_no)
74  {
75  AssertIndexRange(block_no, vector.n_blocks());
76  return vector.block(block_no);
77  }
78 
79  template <typename VectorType>
80  typename std::enable_if<!IsBlockVector<VectorType>::value,
81  VectorType &>::type
82  subblock(VectorType &vector, unsigned int)
83  {
84  return vector;
85  }
86 
87  template <typename VectorType>
88  typename std::enable_if<!IsBlockVector<VectorType>::value,
89  const VectorType &>::type
90  subblock(const VectorType &vector, unsigned int)
91  {
92  return vector;
93  }
94 
95  template <typename VectorType>
96  typename std::enable_if<IsBlockVector<VectorType>::value, void>::type
97  collect_sizes(VectorType &vector)
98  {
99  vector.collect_sizes();
100  }
101 
102  template <typename VectorType>
103  typename std::enable_if<!IsBlockVector<VectorType>::value, void>::type
104  collect_sizes(const VectorType &)
105  {}
106  } // namespace BlockHelper
107 
187  template <int dim,
188  typename VectorType = LinearAlgebra::distributed::Vector<double>,
189  typename VectorizedArrayType =
191  class Base : public Subscriptor
192  {
193  public:
197  using value_type = typename VectorType::value_type;
198 
202  using size_type = typename VectorType::size_type;
203 
207  Base();
208 
212  virtual ~Base() override = default;
213 
218  virtual void
219  clear();
220 
237  void
238  initialize(std::shared_ptr<
240  const std::vector<unsigned int> &selected_row_blocks =
241  std::vector<unsigned int>(),
242  const std::vector<unsigned int> &selected_column_blocks =
243  std::vector<unsigned int>());
244 
258  void
259  initialize(std::shared_ptr<
261  const MGConstrainedDoFs & mg_constrained_dofs,
262  const unsigned int level,
263  const std::vector<unsigned int> &selected_row_blocks =
264  std::vector<unsigned int>());
265 
280  void
281  initialize(std::shared_ptr<
283  const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
284  const unsigned int level,
285  const std::vector<unsigned int> & selected_row_blocks =
286  std::vector<unsigned int>());
287 
291  size_type
292  m() const;
293 
297  size_type
298  n() const;
299 
303  void
304  vmult_interface_down(VectorType &dst, const VectorType &src) const;
305 
309  void
310  vmult_interface_up(VectorType &dst, const VectorType &src) const;
311 
315  void
316  vmult(VectorType &dst, const VectorType &src) const;
317 
321  void
322  Tvmult(VectorType &dst, const VectorType &src) const;
323 
327  void
328  vmult_add(VectorType &dst, const VectorType &src) const;
329 
333  void
334  Tvmult_add(VectorType &dst, const VectorType &src) const;
335 
340  value_type
341  el(const unsigned int row, const unsigned int col) const;
342 
347  virtual std::size_t
348  memory_consumption() const;
349 
353  void
354  initialize_dof_vector(VectorType &vec) const;
355 
363  virtual void
364  compute_diagonal() = 0;
365 
369  std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
370  get_matrix_free() const;
371 
375  const std::shared_ptr<DiagonalMatrix<VectorType>> &
376  get_matrix_diagonal_inverse() const;
377 
381  const std::shared_ptr<DiagonalMatrix<VectorType>> &
382  get_matrix_diagonal() const;
383 
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 
538  template <typename OperatorType>
540  {
541  public:
545  using value_type = typename OperatorType::value_type;
546 
550  using size_type = typename OperatorType::size_type;
551 
556 
560  void
561  clear();
562 
566  void
567  initialize(const OperatorType &operator_in);
568 
572  template <typename VectorType>
573  void
574  vmult(VectorType &dst, const VectorType &src) const;
575 
579  template <typename VectorType>
580  void
581  Tvmult(VectorType &dst, const VectorType &src) const;
582 
586  template <typename VectorType>
587  void
588  initialize_dof_vector(VectorType &vec) const;
589 
590 
591  private:
596  };
597 
598 
599 
619  template <int dim,
620  int fe_degree,
621  int n_components = 1,
622  typename Number = double,
623  typename VectorizedArrayType = VectorizedArray<Number>>
625  {
626  static_assert(
627  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
628  "Type of Number and of VectorizedArrayType do not match.");
629 
630  public:
636  const FEEvaluationBase<dim,
637  n_components,
638  Number,
639  false,
640  VectorizedArrayType> &fe_eval);
641 
650  void
651  apply(const AlignedVector<VectorizedArrayType> &inverse_coefficient,
652  const unsigned int n_actual_components,
653  const VectorizedArrayType * in_array,
654  VectorizedArrayType * out_array) const;
655 
667  void
668  apply(const VectorizedArrayType *in_array,
669  VectorizedArrayType * out_array) const;
670 
704  void
705  transform_from_q_points_to_basis(const unsigned int n_actual_components,
706  const VectorizedArrayType *in_array,
707  VectorizedArrayType *out_array) const;
708 
714  void
715  fill_inverse_JxW_values(
716  AlignedVector<VectorizedArrayType> &inverse_jxw) const;
717 
718  private:
722  const FEEvaluationBase<dim,
723  n_components,
724  Number,
725  false,
726  VectorizedArrayType> &fe_eval;
727  };
728 
729 
730 
740  template <int dim,
741  int fe_degree,
742  int n_q_points_1d = fe_degree + 1,
743  int n_components = 1,
744  typename VectorType = LinearAlgebra::distributed::Vector<double>,
745  typename VectorizedArrayType =
747  class MassOperator : public Base<dim, VectorType, VectorizedArrayType>
748  {
749  public:
753  using value_type =
755 
759  using size_type =
761 
765  MassOperator();
766 
771  virtual void
772  compute_diagonal() override;
773 
774  private:
780  virtual void
781  apply_add(VectorType &dst, const VectorType &src) const override;
782 
786  void
787  local_apply_cell(
789  VectorType & dst,
790  const VectorType & src,
791  const std::pair<unsigned int, unsigned int> & cell_range) const;
792  };
793 
794 
795 
808  template <int dim,
809  int fe_degree,
810  int n_q_points_1d = fe_degree + 1,
811  int n_components = 1,
812  typename VectorType = LinearAlgebra::distributed::Vector<double>,
813  typename VectorizedArrayType =
815  class LaplaceOperator : public Base<dim, VectorType, VectorizedArrayType>
816  {
817  public:
821  using value_type =
823 
827  using size_type =
829 
833  LaplaceOperator();
834 
841  virtual void
842  compute_diagonal();
843 
890  void
891  set_coefficient(
892  const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient);
893 
894  virtual void
895  clear();
896 
903  std::shared_ptr<Table<2, VectorizedArrayType>>
904  get_coefficient();
905 
906  private:
912  virtual void
913  apply_add(VectorType &dst, const VectorType &src) const;
914 
918  void
919  local_apply_cell(
921  VectorType & dst,
922  const VectorType & src,
923  const std::pair<unsigned int, unsigned int> & cell_range) const;
924 
928  void
929  local_diagonal_cell(
931  VectorType & dst,
932  const VectorType &,
933  const std::pair<unsigned int, unsigned int> &cell_range) const;
934 
938  void
939  do_operation_on_cell(
941  & phi,
942  const unsigned int cell) const;
943 
947  std::shared_ptr<Table<2, VectorizedArrayType>> scalar_coefficient;
948  };
949 
950 
951 
952  // ------------------------------------ inline functions ---------------------
953 
954  template <int dim,
955  int fe_degree,
956  int n_components,
957  typename Number,
958  typename VectorizedArrayType>
959  inline CellwiseInverseMassMatrix<dim,
960  fe_degree,
961  n_components,
962  Number,
963  VectorizedArrayType>::
964  CellwiseInverseMassMatrix(
965  const FEEvaluationBase<dim,
966  n_components,
967  Number,
968  false,
969  VectorizedArrayType> &fe_eval)
970  : fe_eval(fe_eval)
971  {}
972 
973 
974 
975  template <int dim,
976  int fe_degree,
977  int n_components,
978  typename Number,
979  typename VectorizedArrayType>
980  inline void
982  fe_degree,
983  n_components,
984  Number,
985  VectorizedArrayType>::
987  AlignedVector<VectorizedArrayType> &inverse_jxw) const
988  {
989  constexpr unsigned int dofs_per_component_on_cell =
990  Utilities::pow(fe_degree + 1, dim);
991  Assert(inverse_jxw.size() > 0 &&
992  inverse_jxw.size() % dofs_per_component_on_cell == 0,
993  ExcMessage(
994  "Expected diagonal to be a multiple of scalar dof per cells"));
995 
996  // temporarily reduce size of inverse_jxw to dofs_per_cell to get JxW values
997  // from fe_eval (will not reallocate any memory)
998  for (unsigned int q = 0; q < dofs_per_component_on_cell; ++q)
999  inverse_jxw[q] = 1. / fe_eval.JxW(q);
1000  // copy values to rest of vector
1001  for (unsigned int q = dofs_per_component_on_cell; q < inverse_jxw.size();)
1002  for (unsigned int i = 0; i < dofs_per_component_on_cell; ++i, ++q)
1003  inverse_jxw[q] = inverse_jxw[i];
1004  }
1005 
1006 
1007 
1008  template <int dim,
1009  int fe_degree,
1010  int n_components,
1011  typename Number,
1012  typename VectorizedArrayType>
1013  inline void
1015  dim,
1016  fe_degree,
1017  n_components,
1018  Number,
1019  VectorizedArrayType>::apply(const VectorizedArrayType *in_array,
1020  VectorizedArrayType * out_array) const
1021  {
1023  dim,
1024  fe_degree,
1025  n_components,
1026  VectorizedArrayType>::apply(fe_eval, in_array, out_array);
1027  }
1028 
1029 
1030 
1031  template <int dim,
1032  int fe_degree,
1033  int n_components,
1034  typename Number,
1035  typename VectorizedArrayType>
1036  inline void
1038  fe_degree,
1039  n_components,
1040  Number,
1041  VectorizedArrayType>::
1042  apply(const AlignedVector<VectorizedArrayType> &inverse_coefficients,
1043  const unsigned int n_actual_components,
1044  const VectorizedArrayType * in_array,
1045  VectorizedArrayType * out_array) const
1046  {
1048  fe_degree,
1049  n_components,
1050  VectorizedArrayType>::
1051  apply(fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
1052  inverse_coefficients,
1053  n_actual_components,
1054  in_array,
1055  out_array);
1056  }
1057 
1058 
1059 
1060  template <int dim,
1061  int fe_degree,
1062  int n_components,
1063  typename Number,
1064  typename VectorizedArrayType>
1065  inline void
1067  fe_degree,
1068  n_components,
1069  Number,
1070  VectorizedArrayType>::
1071  transform_from_q_points_to_basis(const unsigned int n_actual_components,
1072  const VectorizedArrayType *in_array,
1073  VectorizedArrayType * out_array) const
1074  {
1076  fe_degree,
1077  n_components,
1078  VectorizedArrayType>::
1080  fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
1081  n_actual_components,
1082  in_array,
1083  out_array);
1084  }
1085 
1086 
1087 
1088  //----------------- Base operator -----------------------------
1089  template <int dim, typename VectorType, typename VectorizedArrayType>
1091  : Subscriptor()
1092  , have_interface_matrices(false)
1093  {}
1094 
1095 
1096 
1097  template <int dim, typename VectorType, typename VectorizedArrayType>
1100  {
1101  Assert(data.get() != nullptr, ExcNotInitialized());
1103  0;
1104  for (unsigned int i = 0; i < selected_rows.size(); ++i)
1105  total_size += data->get_vector_partitioner(selected_rows[i])->size();
1106  return total_size;
1107  }
1108 
1109 
1110 
1111  template <int dim, typename VectorType, typename VectorizedArrayType>
1114  {
1115  Assert(data.get() != nullptr, ExcNotInitialized());
1117  0;
1118  for (unsigned int i = 0; i < selected_columns.size(); ++i)
1119  total_size += data->get_vector_partitioner(selected_columns[i])->size();
1120  return total_size;
1121  }
1122 
1123 
1124 
1125  template <int dim, typename VectorType, typename VectorizedArrayType>
1126  void
1128  {
1129  data.reset();
1130  inverse_diagonal_entries.reset();
1131  }
1132 
1133 
1134 
1135  template <int dim, typename VectorType, typename VectorizedArrayType>
1138  const unsigned int col) const
1139  {
1140  (void)col;
1141  Assert(row == col, ExcNotImplemented());
1142  Assert(inverse_diagonal_entries.get() != nullptr &&
1143  inverse_diagonal_entries->m() > 0,
1144  ExcNotInitialized());
1145  return 1.0 / (*inverse_diagonal_entries)(row, row);
1146  }
1147 
1148 
1149 
1150  template <int dim, typename VectorType, typename VectorizedArrayType>
1151  void
1153  VectorType &vec) const
1154  {
1155  Assert(data.get() != nullptr, ExcNotInitialized());
1156  AssertDimension(BlockHelper::n_blocks(vec), selected_rows.size());
1157  for (unsigned int i = 0; i < BlockHelper::n_blocks(vec); ++i)
1158  {
1159  const unsigned int index = selected_rows[i];
1160  if (!BlockHelper::subblock(vec, index)
1161  .partitioners_are_compatible(
1162  *data->get_dof_info(index).vector_partitioner))
1163  data->initialize_dof_vector(BlockHelper::subblock(vec, index), index);
1164 
1165  Assert(BlockHelper::subblock(vec, index)
1166  .partitioners_are_globally_compatible(
1167  *data->get_dof_info(index).vector_partitioner),
1168  ExcInternalError());
1169  }
1170  BlockHelper::collect_sizes(vec);
1171  }
1172 
1173 
1174 
1175  template <int dim, typename VectorType, typename VectorizedArrayType>
1176  void
1179  data_,
1180  const std::vector<unsigned int> &given_row_selection,
1181  const std::vector<unsigned int> &given_column_selection)
1182  {
1183  data = data_;
1184 
1185  selected_rows.clear();
1186  selected_columns.clear();
1187  if (given_row_selection.empty())
1188  for (unsigned int i = 0; i < data_->n_components(); ++i)
1189  selected_rows.push_back(i);
1190  else
1191  {
1192  for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1193  {
1194  AssertIndexRange(given_row_selection[i], data_->n_components());
1195  for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1196  if (j != i)
1197  Assert(given_row_selection[j] != given_row_selection[i],
1198  ExcMessage("Given row indices must be unique"));
1199 
1200  selected_rows.push_back(given_row_selection[i]);
1201  }
1202  }
1203  if (given_column_selection.size() == 0)
1205  else
1206  {
1207  for (unsigned int i = 0; i < given_column_selection.size(); ++i)
1208  {
1209  AssertIndexRange(given_column_selection[i], data_->n_components());
1210  for (unsigned int j = 0; j < given_column_selection.size(); ++j)
1211  if (j != i)
1212  Assert(given_column_selection[j] != given_column_selection[i],
1213  ExcMessage("Given column indices must be unique"));
1214 
1215  selected_columns.push_back(given_column_selection[i]);
1216  }
1217  }
1218 
1219  edge_constrained_indices.clear();
1220  edge_constrained_indices.resize(selected_rows.size());
1221  edge_constrained_values.clear();
1222  edge_constrained_values.resize(selected_rows.size());
1223  have_interface_matrices = false;
1224  }
1225 
1226 
1227 
1228  template <int dim, typename VectorType, typename VectorizedArrayType>
1229  void
1232  data_,
1233  const MGConstrainedDoFs & mg_constrained_dofs,
1234  const unsigned int level,
1235  const std::vector<unsigned int> &given_row_selection)
1236  {
1237  std::vector<MGConstrainedDoFs> mg_constrained_dofs_vector(
1238  1, mg_constrained_dofs);
1239  initialize(data_, mg_constrained_dofs_vector, level, given_row_selection);
1240  }
1241 
1242 
1243 
1244  template <int dim, typename VectorType, typename VectorizedArrayType>
1245  void
1248  data_,
1249  const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
1250  const unsigned int level,
1251  const std::vector<unsigned int> & given_row_selection)
1252  {
1254  ExcMessage("level is not set"));
1255 
1256  selected_rows.clear();
1257  selected_columns.clear();
1258  if (given_row_selection.empty())
1259  for (unsigned int i = 0; i < data_->n_components(); ++i)
1260  selected_rows.push_back(i);
1261  else
1262  {
1263  for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1264  {
1265  AssertIndexRange(given_row_selection[i], data_->n_components());
1266  for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1267  if (j != i)
1268  Assert(given_row_selection[j] != given_row_selection[i],
1269  ExcMessage("Given row indices must be unique"));
1270 
1271  selected_rows.push_back(given_row_selection[i]);
1272  }
1273  }
1275 
1276  AssertDimension(mg_constrained_dofs.size(), selected_rows.size());
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 
1282  data = data_;
1283 
1284  for (unsigned int j = 0; j < selected_rows.size(); ++j)
1285  {
1286  if (data_->n_macro_cells() > 0)
1287  {
1288  AssertDimension(level, data_->get_cell_iterator(0, 0, j)->level());
1289  }
1290 
1291  // setup edge_constrained indices
1292  std::vector<types::global_dof_index> interface_indices;
1293  mg_constrained_dofs[j]
1294  .get_refinement_edge_indices(level)
1295  .fill_index_vector(interface_indices);
1296  edge_constrained_indices[j].clear();
1297  edge_constrained_indices[j].reserve(interface_indices.size());
1298  edge_constrained_values[j].resize(interface_indices.size());
1299  const IndexSet &locally_owned =
1300  data->get_dof_handler(selected_rows[j]).locally_owned_mg_dofs(level);
1301  for (const auto interface_index : interface_indices)
1302  if (locally_owned.is_element(interface_index))
1303  edge_constrained_indices[j].push_back(
1304  locally_owned.index_within_set(interface_index));
1307  static_cast<unsigned int>(edge_constrained_indices[j].size()),
1308  data->get_vector_partitioner()->get_mpi_communicator()) > 0;
1309  }
1310  }
1311 
1312 
1313 
1314  template <int dim, typename VectorType, typename VectorizedArrayType>
1315  void
1317  VectorType &dst) const
1318  {
1319  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1320  {
1321  const std::vector<unsigned int> &constrained_dofs =
1322  data->get_constrained_dofs(selected_rows[j]);
1323  for (const auto constrained_dof : constrained_dofs)
1324  BlockHelper::subblock(dst, j).local_element(constrained_dof) = 1.;
1325  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1326  BlockHelper::subblock(dst, j).local_element(
1327  edge_constrained_indices[j][i]) = 1.;
1328  }
1329  }
1330 
1331 
1332 
1333  template <int dim, typename VectorType, typename VectorizedArrayType>
1334  void
1336  const VectorType &src) const
1337  {
1338  using Number =
1340  dst = Number(0.);
1341  vmult_add(dst, src);
1342  }
1343 
1344 
1345 
1346  template <int dim, typename VectorType, typename VectorizedArrayType>
1347  void
1349  VectorType & dst,
1350  const VectorType &src) const
1351  {
1352  mult_add(dst, src, false);
1353  }
1354 
1355 
1356 
1357  template <int dim, typename VectorType, typename VectorizedArrayType>
1358  void
1360  VectorType & dst,
1361  const VectorType &src) const
1362  {
1363  mult_add(dst, src, true);
1364  }
1365 
1366 
1367 
1368  template <int dim, typename VectorType, typename VectorizedArrayType>
1369  void
1371  const VectorType &src,
1372  const bool is_row) const
1373  {
1374  using Number =
1376  for (unsigned int i = 0; i < BlockHelper::n_blocks(src); ++i)
1377  {
1378  const unsigned int mf_component =
1379  is_row ? selected_rows[i] : selected_columns[i];
1380  // If both vectors use the same partitioner -> done
1381  if (BlockHelper::subblock(src, i).get_partitioner().get() ==
1382  data->get_dof_info(mf_component).vector_partitioner.get())
1383  continue;
1384 
1385  // If not, assert that the local ranges are the same and reset to the
1386  // current partitioner
1387  Assert(
1388  BlockHelper::subblock(src, i).get_partitioner()->local_size() ==
1389  data->get_dof_info(mf_component).vector_partitioner->local_size(),
1390  ExcMessage("The vector passed to the vmult() function does not have "
1391  "the correct size for compatibility with MatrixFree."));
1392 
1393  // copy the vector content to a temporary vector so that it does not get
1394  // lost
1396  BlockHelper::subblock(src, i));
1397  BlockHelper::subblock(const_cast<VectorType &>(src), i)
1398  .reinit(data->get_dof_info(mf_component).vector_partitioner);
1399  BlockHelper::subblock(const_cast<VectorType &>(src), i)
1400  .copy_locally_owned_data_from(copy_vec);
1401  }
1402  }
1403 
1404 
1405 
1406  template <int dim, typename VectorType, typename VectorizedArrayType>
1407  void
1409  VectorType & dst,
1410  const VectorType &src) const
1411  {
1412  using Number =
1414  adjust_ghost_range_if_necessary(src, false);
1416 
1417  // set zero Dirichlet values on the input vector (and remember the src and
1418  // dst values because we need to reset them at the end)
1419  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1420  {
1421  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1422  {
1423  edge_constrained_values[j][i] = std::pair<Number, Number>(
1424  BlockHelper::subblock(src, j).local_element(
1425  edge_constrained_indices[j][i]),
1426  BlockHelper::subblock(dst, j).local_element(
1427  edge_constrained_indices[j][i]));
1428  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1429  .local_element(edge_constrained_indices[j][i]) = 0.;
1430  }
1431  }
1432  }
1433 
1434 
1435 
1436  template <int dim, typename VectorType, typename VectorizedArrayType>
1437  void
1439  VectorType & dst,
1440  const VectorType &src,
1441  const bool transpose) const
1442  {
1443  AssertDimension(dst.size(), src.size());
1444  AssertDimension(BlockHelper::n_blocks(dst), BlockHelper::n_blocks(src));
1445  AssertDimension(BlockHelper::n_blocks(dst), selected_rows.size());
1446  preprocess_constraints(dst, src);
1447  if (transpose)
1448  Tapply_add(dst, src);
1449  else
1450  apply_add(dst, src);
1451  postprocess_constraints(dst, src);
1452  }
1453 
1454 
1455 
1456  template <int dim, typename VectorType, typename VectorizedArrayType>
1457  void
1459  VectorType & dst,
1460  const VectorType &src) const
1461  {
1462  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1463  {
1464  const std::vector<unsigned int> &constrained_dofs =
1465  data->get_constrained_dofs(selected_rows[j]);
1466  for (const auto constrained_dof : constrained_dofs)
1467  BlockHelper::subblock(dst, j).local_element(constrained_dof) +=
1468  BlockHelper::subblock(src, j).local_element(constrained_dof);
1469  }
1470 
1471  // reset edge constrained values, multiply by unit matrix and add into
1472  // destination
1473  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1474  {
1475  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1476  {
1477  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1478  .local_element(edge_constrained_indices[j][i]) =
1479  edge_constrained_values[j][i].first;
1480  BlockHelper::subblock(dst, j).local_element(
1481  edge_constrained_indices[j][i]) =
1482  edge_constrained_values[j][i].second +
1483  edge_constrained_values[j][i].first;
1484  }
1485  }
1486  }
1487 
1488 
1489 
1490  template <int dim, typename VectorType, typename VectorizedArrayType>
1491  void
1493  VectorType & dst,
1494  const VectorType &src) const
1495  {
1496  using Number =
1498  AssertDimension(dst.size(), src.size());
1499  adjust_ghost_range_if_necessary(src, false);
1501 
1502  dst = Number(0.);
1503 
1505  return;
1506 
1507  // set zero Dirichlet values on the input vector (and remember the src and
1508  // dst values because we need to reset them at the end)
1509  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1510  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1511  {
1512  edge_constrained_values[j][i] = std::pair<Number, Number>(
1513  BlockHelper::subblock(src, j).local_element(
1514  edge_constrained_indices[j][i]),
1515  BlockHelper::subblock(dst, j).local_element(
1516  edge_constrained_indices[j][i]));
1517  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1518  .local_element(edge_constrained_indices[j][i]) = 0.;
1519  }
1520 
1521  apply_add(dst, src);
1522 
1523  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1524  {
1525  unsigned int c = 0;
1526  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1527  {
1528  for (; c < edge_constrained_indices[j][i]; ++c)
1529  BlockHelper::subblock(dst, j).local_element(c) = 0.;
1530  ++c;
1531 
1532  // reset the src values
1533  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1534  .local_element(edge_constrained_indices[j][i]) =
1535  edge_constrained_values[j][i].first;
1536  }
1537  for (; c < BlockHelper::subblock(dst, j).local_size(); ++c)
1538  BlockHelper::subblock(dst, j).local_element(c) = 0.;
1539  }
1540  }
1541 
1542 
1543 
1544  template <int dim, typename VectorType, typename VectorizedArrayType>
1545  void
1547  VectorType & dst,
1548  const VectorType &src) const
1549  {
1550  using Number =
1552  AssertDimension(dst.size(), src.size());
1553  adjust_ghost_range_if_necessary(src, false);
1555 
1556  dst = Number(0.);
1557 
1559  return;
1560 
1561  VectorType src_cpy(src);
1562  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1563  {
1564  unsigned int c = 0;
1565  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1566  {
1567  for (; c < edge_constrained_indices[j][i]; ++c)
1568  BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1569  ++c;
1570  }
1571  for (; c < BlockHelper::subblock(src_cpy, j).local_size(); ++c)
1572  BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1573  }
1574 
1575  apply_add(dst, src_cpy);
1576 
1577  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1578  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1579  BlockHelper::subblock(dst, j).local_element(
1580  edge_constrained_indices[j][i]) = 0.;
1581  }
1582 
1583 
1584 
1585  template <int dim, typename VectorType, typename VectorizedArrayType>
1586  void
1588  VectorType & dst,
1589  const VectorType &src) const
1590  {
1591  using Number =
1593  dst = Number(0.);
1594  Tvmult_add(dst, src);
1595  }
1596 
1597 
1598 
1599  template <int dim, typename VectorType, typename VectorizedArrayType>
1600  std::size_t
1602  {
1603  return inverse_diagonal_entries.get() != nullptr ?
1604  inverse_diagonal_entries->memory_consumption() :
1605  sizeof(*this);
1606  }
1607 
1608 
1609 
1610  template <int dim, typename VectorType, typename VectorizedArrayType>
1611  std::shared_ptr<const MatrixFree<
1612  dim,
1614  VectorizedArrayType>>
1616  {
1617  return data;
1618  }
1619 
1620 
1621 
1622  template <int dim, typename VectorType, typename VectorizedArrayType>
1623  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1625  const
1626  {
1627  Assert(inverse_diagonal_entries.get() != nullptr &&
1628  inverse_diagonal_entries->m() > 0,
1629  ExcNotInitialized());
1630  return inverse_diagonal_entries;
1631  }
1632 
1633 
1634 
1635  template <int dim, typename VectorType, typename VectorizedArrayType>
1636  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1638  {
1639  Assert(diagonal_entries.get() != nullptr && diagonal_entries->m() > 0,
1640  ExcNotInitialized());
1641  return diagonal_entries;
1642  }
1643 
1644 
1645 
1646  template <int dim, typename VectorType, typename VectorizedArrayType>
1647  void
1649  VectorType & dst,
1650  const VectorType &src) const
1651  {
1652  apply_add(dst, src);
1653  }
1654 
1655 
1656 
1657  template <int dim, typename VectorType, typename VectorizedArrayType>
1658  void
1660  VectorType & dst,
1661  const VectorType & src,
1662  const typename Base<dim, VectorType, VectorizedArrayType>::value_type omega)
1663  const
1664  {
1666  ExcNotInitialized());
1667  inverse_diagonal_entries->vmult(dst, src);
1668  dst *= omega;
1669  }
1670 
1671 
1672 
1673  //------------------------- MGInterfaceOperator ------------------------------
1674 
1675  template <typename OperatorType>
1677  : Subscriptor()
1678  , mf_base_operator(nullptr)
1679  {}
1680 
1681 
1682 
1683  template <typename OperatorType>
1684  void
1686  {
1687  mf_base_operator = nullptr;
1688  }
1689 
1690 
1691 
1692  template <typename OperatorType>
1693  void
1694  MGInterfaceOperator<OperatorType>::initialize(const OperatorType &operator_in)
1695  {
1696  mf_base_operator = &operator_in;
1697  }
1698 
1699 
1700 
1701  template <typename OperatorType>
1702  template <typename VectorType>
1703  void
1705  const VectorType &src) const
1706  {
1707 #ifndef DEAL_II_MSVC
1708  static_assert(
1709  std::is_same<typename VectorType::value_type, value_type>::value,
1710  "The vector type must be based on the same value type as this "
1711  "operator");
1712 #endif
1713 
1714  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1715 
1716  mf_base_operator->vmult_interface_down(dst, src);
1717  }
1718 
1719 
1720 
1721  template <typename OperatorType>
1722  template <typename VectorType>
1723  void
1725  const VectorType &src) const
1726  {
1727 #ifndef DEAL_II_MSVC
1728  static_assert(
1729  std::is_same<typename VectorType::value_type, value_type>::value,
1730  "The vector type must be based on the same value type as this "
1731  "operator");
1732 #endif
1733 
1734  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1735 
1736  mf_base_operator->vmult_interface_up(dst, src);
1737  }
1738 
1739 
1740 
1741  template <typename OperatorType>
1742  template <typename VectorType>
1743  void
1745  VectorType &vec) const
1746  {
1747  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1748 
1749  mf_base_operator->initialize_dof_vector(vec);
1750  }
1751 
1752 
1753 
1754  //-----------------------------MassOperator----------------------------------
1755 
1756  template <int dim,
1757  int fe_degree,
1758  int n_q_points_1d,
1759  int n_components,
1760  typename VectorType,
1761  typename VectorizedArrayType>
1762  MassOperator<dim,
1763  fe_degree,
1764  n_q_points_1d,
1765  n_components,
1766  VectorType,
1767  VectorizedArrayType>::MassOperator()
1768  : Base<dim, VectorType, VectorizedArrayType>()
1769  {}
1770 
1771 
1772 
1773  template <int dim,
1774  int fe_degree,
1775  int n_q_points_1d,
1776  int n_components,
1777  typename VectorType,
1778  typename VectorizedArrayType>
1779  void
1780  MassOperator<dim,
1781  fe_degree,
1782  n_q_points_1d,
1783  n_components,
1784  VectorType,
1785  VectorizedArrayType>::compute_diagonal()
1786  {
1787  using Number =
1788  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1790  ExcNotInitialized());
1791 
1792  this->inverse_diagonal_entries =
1793  std::make_shared<DiagonalMatrix<VectorType>>();
1794  this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1795  VectorType &inverse_diagonal_vector =
1796  this->inverse_diagonal_entries->get_vector();
1797  VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1798  this->initialize_dof_vector(inverse_diagonal_vector);
1799  this->initialize_dof_vector(diagonal_vector);
1800  inverse_diagonal_vector = Number(1.);
1801  apply_add(diagonal_vector, inverse_diagonal_vector);
1802 
1803  this->set_constrained_entries_to_one(diagonal_vector);
1804  inverse_diagonal_vector = diagonal_vector;
1805 
1806  const unsigned int local_size = inverse_diagonal_vector.local_size();
1807  for (unsigned int i = 0; i < local_size; ++i)
1808  inverse_diagonal_vector.local_element(i) =
1809  Number(1.) / inverse_diagonal_vector.local_element(i);
1810 
1811  inverse_diagonal_vector.update_ghost_values();
1812  diagonal_vector.update_ghost_values();
1813  }
1814 
1815 
1816 
1817  template <int dim,
1818  int fe_degree,
1819  int n_q_points_1d,
1820  int n_components,
1821  typename VectorType,
1822  typename VectorizedArrayType>
1823  void
1824  MassOperator<dim,
1825  fe_degree,
1826  n_q_points_1d,
1827  n_components,
1828  VectorType,
1829  VectorizedArrayType>::apply_add(VectorType & dst,
1830  const VectorType &src) const
1831  {
1833  &MassOperator::local_apply_cell, this, dst, src);
1834  }
1835 
1836 
1837 
1838  template <int dim,
1839  int fe_degree,
1840  int n_q_points_1d,
1841  int n_components,
1842  typename VectorType,
1843  typename VectorizedArrayType>
1844  void
1845  MassOperator<dim,
1846  fe_degree,
1847  n_q_points_1d,
1848  n_components,
1849  VectorType,
1850  VectorizedArrayType>::
1852  const MatrixFree<
1853  dim,
1854  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
1855  VectorizedArrayType> & data,
1856  VectorType & dst,
1857  const VectorType & src,
1858  const std::pair<unsigned int, unsigned int> &cell_range) const
1859  {
1860  using Number =
1861  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1862  FEEvaluation<dim,
1863  fe_degree,
1864  n_q_points_1d,
1865  n_components,
1866  Number,
1867  VectorizedArrayType>
1868  phi(data, this->selected_rows[0]);
1869  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
1870  {
1871  phi.reinit(cell);
1872  phi.read_dof_values(src);
1873  phi.evaluate(true, false, false);
1874  for (unsigned int q = 0; q < phi.n_q_points; ++q)
1875  phi.submit_value(phi.get_value(q), q);
1876  phi.integrate(true, false);
1877  phi.distribute_local_to_global(dst);
1878  }
1879  }
1880 
1881 
1882  //-----------------------------LaplaceOperator----------------------------------
1883 
1884  template <int dim,
1885  int fe_degree,
1886  int n_q_points_1d,
1887  int n_components,
1888  typename VectorType,
1889  typename VectorizedArrayType>
1890  LaplaceOperator<dim,
1891  fe_degree,
1892  n_q_points_1d,
1893  n_components,
1894  VectorType,
1895  VectorizedArrayType>::LaplaceOperator()
1896  : Base<dim, VectorType, VectorizedArrayType>()
1897  {}
1898 
1899 
1900 
1901  template <int dim,
1902  int fe_degree,
1903  int n_q_points_1d,
1904  int n_components,
1905  typename VectorType,
1906  typename VectorizedArrayType>
1907  void
1908  LaplaceOperator<dim,
1909  fe_degree,
1910  n_q_points_1d,
1911  n_components,
1912  VectorType,
1913  VectorizedArrayType>::clear()
1914  {
1916  scalar_coefficient.reset();
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  LaplaceOperator<dim,
1929  fe_degree,
1930  n_q_points_1d,
1931  n_components,
1932  VectorType,
1933  VectorizedArrayType>::
1935  const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient_)
1936  {
1937  scalar_coefficient = scalar_coefficient_;
1938  }
1939 
1940 
1941 
1942  template <int dim,
1943  int fe_degree,
1944  int n_q_points_1d,
1945  int n_components,
1946  typename VectorType,
1947  typename VectorizedArrayType>
1948  std::shared_ptr<Table<2, VectorizedArrayType>>
1949  LaplaceOperator<dim,
1950  fe_degree,
1951  n_q_points_1d,
1952  n_components,
1953  VectorType,
1954  VectorizedArrayType>::get_coefficient()
1955  {
1957  return scalar_coefficient;
1958  }
1959 
1960 
1961 
1962  template <int dim,
1963  int fe_degree,
1964  int n_q_points_1d,
1965  int n_components,
1966  typename VectorType,
1967  typename VectorizedArrayType>
1968  void
1969  LaplaceOperator<dim,
1970  fe_degree,
1971  n_q_points_1d,
1972  n_components,
1973  VectorType,
1974  VectorizedArrayType>::compute_diagonal()
1975  {
1976  using Number =
1977  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1979  ExcNotInitialized());
1980 
1981  this->inverse_diagonal_entries =
1982  std::make_shared<DiagonalMatrix<VectorType>>();
1983  this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1984  VectorType &inverse_diagonal_vector =
1985  this->inverse_diagonal_entries->get_vector();
1986  VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1987  this->initialize_dof_vector(inverse_diagonal_vector);
1988  this->initialize_dof_vector(diagonal_vector);
1989 
1990  this->data->cell_loop(&LaplaceOperator::local_diagonal_cell,
1991  this,
1992  diagonal_vector,
1993  /*unused*/ diagonal_vector);
1994  this->set_constrained_entries_to_one(diagonal_vector);
1995 
1996  inverse_diagonal_vector = diagonal_vector;
1997 
1998  for (unsigned int i = 0; i < inverse_diagonal_vector.local_size(); ++i)
1999  if (std::abs(inverse_diagonal_vector.local_element(i)) >
2000  std::sqrt(std::numeric_limits<Number>::epsilon()))
2001  inverse_diagonal_vector.local_element(i) =
2002  1. / inverse_diagonal_vector.local_element(i);
2003  else
2004  inverse_diagonal_vector.local_element(i) = 1.;
2005 
2006  inverse_diagonal_vector.update_ghost_values();
2007  diagonal_vector.update_ghost_values();
2008  }
2009 
2010 
2011 
2012  template <int dim,
2013  int fe_degree,
2014  int n_q_points_1d,
2015  int n_components,
2016  typename VectorType,
2017  typename VectorizedArrayType>
2018  void
2019  LaplaceOperator<dim,
2020  fe_degree,
2021  n_q_points_1d,
2022  n_components,
2023  VectorType,
2024  VectorizedArrayType>::apply_add(VectorType & dst,
2025  const VectorType &src) const
2026  {
2028  &LaplaceOperator::local_apply_cell, this, dst, src);
2029  }
2030 
2031  namespace Implementation
2032  {
2033  template <typename VectorizedArrayType>
2034  bool
2035  non_negative(const VectorizedArrayType &n)
2036  {
2037  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2038  if (n[v] < 0.)
2039  return false;
2040 
2041  return true;
2042  }
2043  } // namespace Implementation
2044 
2045 
2046 
2047  template <int dim,
2048  int fe_degree,
2049  int n_q_points_1d,
2050  int n_components,
2051  typename VectorType,
2052  typename VectorizedArrayType>
2053  void
2054  LaplaceOperator<dim,
2055  fe_degree,
2056  n_q_points_1d,
2057  n_components,
2058  VectorType,
2059  VectorizedArrayType>::
2061  FEEvaluation<
2062  dim,
2063  fe_degree,
2064  n_q_points_1d,
2065  n_components,
2066  typename Base<dim, VectorType, VectorizedArrayType>::value_type> &phi,
2067  const unsigned int cell) const
2068  {
2069  phi.evaluate(false, true, false);
2070  if (scalar_coefficient.get())
2071  {
2072  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2073  {
2074  Assert(Implementation::non_negative((*scalar_coefficient)(cell, q)),
2075  ExcMessage("Coefficient must be non-negative"));
2076  phi.submit_gradient((*scalar_coefficient)(cell, q) *
2077  phi.get_gradient(q),
2078  q);
2079  }
2080  }
2081  else
2082  {
2083  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2084  {
2085  phi.submit_gradient(phi.get_gradient(q), q);
2086  }
2087  }
2088  phi.integrate(false, true);
2089  }
2090 
2091 
2092 
2093  template <int dim,
2094  int fe_degree,
2095  int n_q_points_1d,
2096  int n_components,
2097  typename VectorType,
2098  typename VectorizedArrayType>
2099  void
2100  LaplaceOperator<dim,
2101  fe_degree,
2102  n_q_points_1d,
2103  n_components,
2104  VectorType,
2105  VectorizedArrayType>::
2107  const MatrixFree<
2108  dim,
2109  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2110  VectorizedArrayType> & data,
2111  VectorType & dst,
2112  const VectorType & src,
2113  const std::pair<unsigned int, unsigned int> &cell_range) const
2114  {
2115  using Number =
2116  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2118  data, this->selected_rows[0]);
2119  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2120  {
2121  phi.reinit(cell);
2122  phi.read_dof_values(src);
2123  do_operation_on_cell(phi, cell);
2124  phi.distribute_local_to_global(dst);
2125  }
2126  }
2127 
2128 
2129  template <int dim,
2130  int fe_degree,
2131  int n_q_points_1d,
2132  int n_components,
2133  typename VectorType,
2134  typename VectorizedArrayType>
2135  void
2136  LaplaceOperator<dim,
2137  fe_degree,
2138  n_q_points_1d,
2139  n_components,
2140  VectorType,
2141  VectorizedArrayType>::
2143  const MatrixFree<
2144  dim,
2145  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2146  VectorizedArrayType> &data,
2147  VectorType & dst,
2148  const VectorType &,
2149  const std::pair<unsigned int, unsigned int> &cell_range) const
2150  {
2151  using Number =
2152  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2153 
2155  data, this->selected_rows[0]);
2156  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2157  {
2158  phi.reinit(cell);
2159  VectorizedArrayType local_diagonal_vector[phi.static_dofs_per_cell];
2160  for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2161  {
2162  for (unsigned int j = 0; j < phi.dofs_per_component; ++j)
2163  phi.begin_dof_values()[j] = VectorizedArrayType();
2164  phi.begin_dof_values()[i] = 1.;
2165  do_operation_on_cell(phi, cell);
2166  local_diagonal_vector[i] = phi.begin_dof_values()[i];
2167  }
2168  for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2169  for (unsigned int c = 0; c < phi.n_components; ++c)
2170  phi.begin_dof_values()[i + c * phi.dofs_per_component] =
2171  local_diagonal_vector[i];
2172  phi.distribute_local_to_global(dst);
2173  }
2174  }
2175 
2176 
2177 } // end of namespace MatrixFreeOperators
2178 
2179 
2180 DEAL_II_NAMESPACE_CLOSE
2181 
2182 #endif
typename OperatorType::value_type value_type
Definition: operators.h:545
VectorizedArrayType JxW(const unsigned int q_index) const
virtual void apply_add(VectorType &dst, const VectorType &src) const =0
void fill_inverse_JxW_values(AlignedVector< VectorizedArrayType > &inverse_jxw) const
Definition: operators.h:986
static const unsigned int invalid_unsigned_int
Definition: types.h:190
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal_inverse() const
Definition: operators.h:1624
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
virtual std::size_t memory_consumption() const
Definition: operators.h:1601
const FEEvaluationBase< dim, n_components, Number, false, VectorizedArrayType > & fe_eval
Definition: operators.h:726
size_type n() const
Definition: operators.h:1113
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:2142
void precondition_Jacobi(VectorType &dst, const VectorType &src, const value_type omega) const
Definition: operators.h:1659
std::shared_ptr< DiagonalMatrix< VectorType > > diagonal_entries
Definition: operators.h:441
void read_dof_values(const VectorType &src, const unsigned int first_index=0)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
std::vector< std::vector< std::pair< value_type, value_type > > > edge_constrained_values
Definition: operators.h:471
static ::ExceptionBase & ExcNotInitialized()
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > data
Definition: operators.h:435
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1587
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
constexpr SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
typename VectorType::value_type value_type
Definition: operators.h:197
typename OperatorType::size_type size_type
Definition: operators.h:550
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:488
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > get_matrix_free() const
Definition: operators.h:1615
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1744
void vmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1348
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal() const
Definition: operators.h:1637
static ::ExceptionBase & ExcMessage(std::string arg1)
void vmult_interface_up(VectorType &dst, const VectorType &src) const
Definition: operators.h:1546
#define Assert(cond, exc)
Definition: exceptions.h:1419
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:2106
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1922
void postprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1458
virtual void Tapply_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1648
std::shared_ptr< Table< 2, VectorizedArrayType > > scalar_coefficient
Definition: operators.h:947
const unsigned int dofs_per_component
std::vector< unsigned int > selected_rows
Definition: operators.h:453
void reinit(const unsigned int cell_batch_index)
virtual void clear()
Definition: operators.h:1127
void set_constrained_entries_to_one(VectorType &dst) const
Definition: operators.h:1316
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1152
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1704
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition: operators.h:1829
void Tvmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1359
static constexpr unsigned int n_components
virtual void apply_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:2024
size_type m() const
Definition: operators.h:1099
static constexpr unsigned int static_dofs_per_cell
typename VectorType::size_type size_type
Definition: operators.h:202
void mult_add(VectorType &dst, const VectorType &src, const bool transpose) const
Definition: operators.h:1438
std::shared_ptr< Table< 2, VectorizedArrayType > > get_coefficient()
Definition: operators.h:1954
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:2060
const VectorizedArrayType * begin_dof_values() const
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info() const
std::vector< unsigned int > selected_columns
Definition: operators.h:459
void vmult_interface_down(VectorType &dst, const VectorType &src) const
Definition: operators.h:1492
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
size_type size() const
static ::ExceptionBase & ExcNotImplemented()
value_type el(const unsigned int row, const unsigned int col) const
Definition: operators.h:1137
Definition: table.h:699
std::vector< UnivariateShapeData< Number > > data
Definition: shape_info.h:355
bool is_element(const size_type index) const
Definition: index_set.h:1766
void initialize(const OperatorType &operator_in)
Definition: operators.h:1694
SmartPointer< const OperatorType > mf_base_operator
Definition: operators.h:595
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1724
void apply(const AlignedVector< VectorizedArrayType > &inverse_coefficient, const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition: operators.h:1042
std::shared_ptr< DiagonalMatrix< VectorType > > inverse_diagonal_entries
Definition: operators.h:447
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1335
virtual void compute_diagonal() override
Definition: operators.h:1785
T max(const T &t, const MPI_Comm &mpi_communicator)
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:1851
std::vector< std::vector< unsigned int > > edge_constrained_indices
Definition: operators.h:465
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:1177
void preprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1408
void transform_from_q_points_to_basis(const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition: operators.h:1071
void set_coefficient(const std::shared_ptr< Table< 2, VectorizedArrayType >> &scalar_coefficient)
Definition: operators.h:1934
static ::ExceptionBase & ExcInternalError()
void adjust_ghost_range_if_necessary(const VectorType &vec, const bool is_row) const
Definition: operators.h:1370