Reference documentation for deal.II version 9.2.0
\(\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 - 2020 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 =
190  VectorizedArray<typename VectorType::value_type>>
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<
239  const MatrixFree<dim, value_type, VectorizedArrayType>> data,
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<
260  const MatrixFree<dim, value_type, VectorizedArrayType>> data,
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<
282  const MatrixFree<dim, value_type, VectorizedArrayType>> data_,
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>>
435  data;
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>>>
471  edge_constrained_values;
472 
477  bool have_interface_matrices;
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>
539  class MGInterfaceOperator : public Subscriptor
540  {
541  public:
545  using value_type = typename OperatorType::value_type;
546 
550  using size_type = typename OperatorType::size_type;
551 
555  MGInterfaceOperator();
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:
595  SmartPointer<const OperatorType> mf_base_operator;
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>>
624  class CellwiseInverseMassMatrix
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:
635  CellwiseInverseMassMatrix(
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 =
746  VectorizedArray<typename VectorType::value_type>>
747  class MassOperator : public Base<dim, VectorType, VectorizedArrayType>
748  {
749  public:
753  using value_type =
754  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
755 
759  using size_type =
760  typename Base<dim, VectorType, VectorizedArrayType>::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(
788  const MatrixFree<dim, value_type, VectorizedArrayType> &data,
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 =
814  VectorizedArray<typename VectorType::value_type>>
815  class LaplaceOperator : public Base<dim, VectorType, VectorizedArrayType>
816  {
817  public:
821  using value_type =
822  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
823 
827  using size_type =
828  typename Base<dim, VectorType, VectorizedArrayType>::size_type;
829 
833  LaplaceOperator();
834 
841  virtual void
842  compute_diagonal() override;
843 
894  void
895  set_coefficient(
896  const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient);
897 
902  virtual void
903  clear() override;
904 
911  std::shared_ptr<Table<2, VectorizedArrayType>>
912  get_coefficient();
913 
914  private:
920  virtual void
921  apply_add(VectorType &dst, const VectorType &src) const override;
922 
926  void
927  local_apply_cell(
928  const MatrixFree<dim, value_type, VectorizedArrayType> &data,
929  VectorType & dst,
930  const VectorType & src,
931  const std::pair<unsigned int, unsigned int> & cell_range) const;
932 
936  void
937  local_diagonal_cell(
938  const MatrixFree<dim, value_type, VectorizedArrayType> &data,
939  VectorType & dst,
940  const VectorType &,
941  const std::pair<unsigned int, unsigned int> &cell_range) const;
942 
946  void
947  do_operation_on_cell(
948  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, value_type>
949  & phi,
950  const unsigned int cell) const;
951 
955  std::shared_ptr<Table<2, VectorizedArrayType>> scalar_coefficient;
956  };
957 
958 
959 
960  // ------------------------------------ inline functions ---------------------
961 
962  template <int dim,
963  int fe_degree,
964  int n_components,
965  typename Number,
966  typename VectorizedArrayType>
967  inline CellwiseInverseMassMatrix<dim,
968  fe_degree,
969  n_components,
970  Number,
971  VectorizedArrayType>::
972  CellwiseInverseMassMatrix(
973  const FEEvaluationBase<dim,
974  n_components,
975  Number,
976  false,
977  VectorizedArrayType> &fe_eval)
978  : fe_eval(fe_eval)
979  {}
980 
981 
982 
983  template <int dim,
984  int fe_degree,
985  int n_components,
986  typename Number,
987  typename VectorizedArrayType>
988  inline void
989  CellwiseInverseMassMatrix<dim,
990  fe_degree,
991  n_components,
992  Number,
993  VectorizedArrayType>::
994  fill_inverse_JxW_values(
995  AlignedVector<VectorizedArrayType> &inverse_jxw) const
996  {
997  constexpr unsigned int dofs_per_component_on_cell =
998  Utilities::pow(fe_degree + 1, dim);
999  Assert(inverse_jxw.size() > 0 &&
1000  inverse_jxw.size() % dofs_per_component_on_cell == 0,
1001  ExcMessage(
1002  "Expected diagonal to be a multiple of scalar dof per cells"));
1003 
1004  // temporarily reduce size of inverse_jxw to dofs_per_cell to get JxW values
1005  // from fe_eval (will not reallocate any memory)
1006  for (unsigned int q = 0; q < dofs_per_component_on_cell; ++q)
1007  inverse_jxw[q] = 1. / fe_eval.JxW(q);
1008  // copy values to rest of vector
1009  for (unsigned int q = dofs_per_component_on_cell; q < inverse_jxw.size();)
1010  for (unsigned int i = 0; i < dofs_per_component_on_cell; ++i, ++q)
1011  inverse_jxw[q] = inverse_jxw[i];
1012  }
1013 
1014 
1015 
1016  template <int dim,
1017  int fe_degree,
1018  int n_components,
1019  typename Number,
1020  typename VectorizedArrayType>
1021  inline void
1022  CellwiseInverseMassMatrix<
1023  dim,
1024  fe_degree,
1025  n_components,
1026  Number,
1027  VectorizedArrayType>::apply(const VectorizedArrayType *in_array,
1028  VectorizedArrayType * out_array) const
1029  {
1030  internal::CellwiseInverseMassMatrixImpl<
1031  dim,
1032  fe_degree,
1033  n_components,
1034  VectorizedArrayType>::apply(fe_eval, in_array, out_array);
1035  }
1036 
1037 
1038 
1039  template <int dim,
1040  int fe_degree,
1041  int n_components,
1042  typename Number,
1043  typename VectorizedArrayType>
1044  inline void
1045  CellwiseInverseMassMatrix<dim,
1046  fe_degree,
1047  n_components,
1048  Number,
1049  VectorizedArrayType>::
1050  apply(const AlignedVector<VectorizedArrayType> &inverse_coefficients,
1051  const unsigned int n_actual_components,
1052  const VectorizedArrayType * in_array,
1053  VectorizedArrayType * out_array) const
1054  {
1055  internal::CellwiseInverseMassMatrixImpl<dim,
1056  fe_degree,
1057  n_components,
1058  VectorizedArrayType>::
1059  apply(fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
1060  inverse_coefficients,
1061  n_actual_components,
1062  in_array,
1063  out_array);
1064  }
1065 
1066 
1067 
1068  template <int dim,
1069  int fe_degree,
1070  int n_components,
1071  typename Number,
1072  typename VectorizedArrayType>
1073  inline void
1074  CellwiseInverseMassMatrix<dim,
1075  fe_degree,
1076  n_components,
1077  Number,
1078  VectorizedArrayType>::
1079  transform_from_q_points_to_basis(const unsigned int n_actual_components,
1080  const VectorizedArrayType *in_array,
1081  VectorizedArrayType * out_array) const
1082  {
1083  internal::CellwiseInverseMassMatrixImpl<dim,
1084  fe_degree,
1085  n_components,
1086  VectorizedArrayType>::
1087  transform_from_q_points_to_basis(
1088  fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
1089  n_actual_components,
1090  in_array,
1091  out_array);
1092  }
1093 
1094 
1095 
1096  //----------------- Base operator -----------------------------
1097  template <int dim, typename VectorType, typename VectorizedArrayType>
1098  Base<dim, VectorType, VectorizedArrayType>::Base()
1099  : Subscriptor()
1100  , have_interface_matrices(false)
1101  {}
1102 
1103 
1104 
1105  template <int dim, typename VectorType, typename VectorizedArrayType>
1106  typename Base<dim, VectorType, VectorizedArrayType>::size_type
1107  Base<dim, VectorType, VectorizedArrayType>::m() const
1108  {
1109  Assert(data.get() != nullptr, ExcNotInitialized());
1110  typename Base<dim, VectorType, VectorizedArrayType>::size_type total_size =
1111  0;
1112  for (const unsigned int selected_row : selected_rows)
1113  total_size += data->get_vector_partitioner(selected_row)->size();
1114  return total_size;
1115  }
1116 
1117 
1118 
1119  template <int dim, typename VectorType, typename VectorizedArrayType>
1120  typename Base<dim, VectorType, VectorizedArrayType>::size_type
1121  Base<dim, VectorType, VectorizedArrayType>::n() const
1122  {
1123  Assert(data.get() != nullptr, ExcNotInitialized());
1124  typename Base<dim, VectorType, VectorizedArrayType>::size_type total_size =
1125  0;
1126  for (const unsigned int selected_column : selected_columns)
1127  total_size += data->get_vector_partitioner(selected_column)->size();
1128  return total_size;
1129  }
1130 
1131 
1132 
1133  template <int dim, typename VectorType, typename VectorizedArrayType>
1134  void
1135  Base<dim, VectorType, VectorizedArrayType>::clear()
1136  {
1137  data.reset();
1138  inverse_diagonal_entries.reset();
1139  }
1140 
1141 
1142 
1143  template <int dim, typename VectorType, typename VectorizedArrayType>
1144  typename Base<dim, VectorType, VectorizedArrayType>::value_type
1145  Base<dim, VectorType, VectorizedArrayType>::el(const unsigned int row,
1146  const unsigned int col) const
1147  {
1148  (void)col;
1149  Assert(row == col, ExcNotImplemented());
1150  Assert(inverse_diagonal_entries.get() != nullptr &&
1151  inverse_diagonal_entries->m() > 0,
1152  ExcNotInitialized());
1153  return 1.0 / (*inverse_diagonal_entries)(row, row);
1154  }
1155 
1156 
1157 
1158  template <int dim, typename VectorType, typename VectorizedArrayType>
1159  void
1160  Base<dim, VectorType, VectorizedArrayType>::initialize_dof_vector(
1161  VectorType &vec) const
1162  {
1163  Assert(data.get() != nullptr, ExcNotInitialized());
1164  AssertDimension(BlockHelper::n_blocks(vec), selected_rows.size());
1165  for (unsigned int i = 0; i < BlockHelper::n_blocks(vec); ++i)
1166  {
1167  const unsigned int index = selected_rows[i];
1168  if (!BlockHelper::subblock(vec, index)
1169  .partitioners_are_compatible(
1170  *data->get_dof_info(index).vector_partitioner))
1171  data->initialize_dof_vector(BlockHelper::subblock(vec, index), index);
1172 
1173  Assert(BlockHelper::subblock(vec, index)
1174  .partitioners_are_globally_compatible(
1175  *data->get_dof_info(index).vector_partitioner),
1176  ExcInternalError());
1177  }
1178  BlockHelper::collect_sizes(vec);
1179  }
1180 
1181 
1182 
1183  template <int dim, typename VectorType, typename VectorizedArrayType>
1184  void
1185  Base<dim, VectorType, VectorizedArrayType>::initialize(
1186  std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
1187  data_,
1188  const std::vector<unsigned int> &given_row_selection,
1189  const std::vector<unsigned int> &given_column_selection)
1190  {
1191  data = data_;
1192 
1193  selected_rows.clear();
1194  selected_columns.clear();
1195  if (given_row_selection.empty())
1196  for (unsigned int i = 0; i < data_->n_components(); ++i)
1197  selected_rows.push_back(i);
1198  else
1199  {
1200  for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1201  {
1202  AssertIndexRange(given_row_selection[i], data_->n_components());
1203  for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1204  if (j != i)
1205  Assert(given_row_selection[j] != given_row_selection[i],
1206  ExcMessage("Given row indices must be unique"));
1207 
1208  selected_rows.push_back(given_row_selection[i]);
1209  }
1210  }
1211  if (given_column_selection.size() == 0)
1212  selected_columns = selected_rows;
1213  else
1214  {
1215  for (unsigned int i = 0; i < given_column_selection.size(); ++i)
1216  {
1217  AssertIndexRange(given_column_selection[i], data_->n_components());
1218  for (unsigned int j = 0; j < given_column_selection.size(); ++j)
1219  if (j != i)
1220  Assert(given_column_selection[j] != given_column_selection[i],
1221  ExcMessage("Given column indices must be unique"));
1222 
1223  selected_columns.push_back(given_column_selection[i]);
1224  }
1225  }
1226 
1227  edge_constrained_indices.clear();
1228  edge_constrained_indices.resize(selected_rows.size());
1229  edge_constrained_values.clear();
1230  edge_constrained_values.resize(selected_rows.size());
1231  have_interface_matrices = false;
1232  }
1233 
1234 
1235 
1236  template <int dim, typename VectorType, typename VectorizedArrayType>
1237  void
1238  Base<dim, VectorType, VectorizedArrayType>::initialize(
1239  std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
1240  data_,
1241  const MGConstrainedDoFs & mg_constrained_dofs,
1242  const unsigned int level,
1243  const std::vector<unsigned int> &given_row_selection)
1244  {
1245  std::vector<MGConstrainedDoFs> mg_constrained_dofs_vector(
1246  1, mg_constrained_dofs);
1247  initialize(data_, mg_constrained_dofs_vector, level, given_row_selection);
1248  }
1249 
1250 
1251 
1252  template <int dim, typename VectorType, typename VectorizedArrayType>
1253  void
1254  Base<dim, VectorType, VectorizedArrayType>::initialize(
1255  std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
1256  data_,
1257  const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
1258  const unsigned int level,
1259  const std::vector<unsigned int> & given_row_selection)
1260  {
1261  AssertThrow(level != numbers::invalid_unsigned_int,
1262  ExcMessage("level is not set"));
1263 
1264  selected_rows.clear();
1265  selected_columns.clear();
1266  if (given_row_selection.empty())
1267  for (unsigned int i = 0; i < data_->n_components(); ++i)
1268  selected_rows.push_back(i);
1269  else
1270  {
1271  for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1272  {
1273  AssertIndexRange(given_row_selection[i], data_->n_components());
1274  for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1275  if (j != i)
1276  Assert(given_row_selection[j] != given_row_selection[i],
1277  ExcMessage("Given row indices must be unique"));
1278 
1279  selected_rows.push_back(given_row_selection[i]);
1280  }
1281  }
1282  selected_columns = selected_rows;
1283 
1284  AssertDimension(mg_constrained_dofs.size(), selected_rows.size());
1285  edge_constrained_indices.clear();
1286  edge_constrained_indices.resize(selected_rows.size());
1287  edge_constrained_values.clear();
1288  edge_constrained_values.resize(selected_rows.size());
1289 
1290  data = data_;
1291 
1292  for (unsigned int j = 0; j < selected_rows.size(); ++j)
1293  {
1294  if (data_->n_macro_cells() > 0)
1295  {
1296  AssertDimension(level, data_->get_cell_iterator(0, 0, j)->level());
1297  }
1298 
1299  // setup edge_constrained indices
1300  std::vector<types::global_dof_index> interface_indices;
1301  mg_constrained_dofs[j]
1302  .get_refinement_edge_indices(level)
1303  .fill_index_vector(interface_indices);
1304  edge_constrained_indices[j].clear();
1305  edge_constrained_indices[j].reserve(interface_indices.size());
1306  edge_constrained_values[j].resize(interface_indices.size());
1307  const IndexSet &locally_owned =
1308  data->get_dof_handler(selected_rows[j]).locally_owned_mg_dofs(level);
1309  for (const auto interface_index : interface_indices)
1310  if (locally_owned.is_element(interface_index))
1311  edge_constrained_indices[j].push_back(
1312  locally_owned.index_within_set(interface_index));
1313  have_interface_matrices |=
1314  Utilities::MPI::max(
1315  static_cast<unsigned int>(edge_constrained_indices[j].size()),
1316  data->get_vector_partitioner()->get_mpi_communicator()) > 0;
1317  }
1318  }
1319 
1320 
1321 
1322  template <int dim, typename VectorType, typename VectorizedArrayType>
1323  void
1324  Base<dim, VectorType, VectorizedArrayType>::set_constrained_entries_to_one(
1325  VectorType &dst) const
1326  {
1327  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1328  {
1329  const std::vector<unsigned int> &constrained_dofs =
1330  data->get_constrained_dofs(selected_rows[j]);
1331  for (const auto constrained_dof : constrained_dofs)
1332  BlockHelper::subblock(dst, j).local_element(constrained_dof) = 1.;
1333  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1334  BlockHelper::subblock(dst, j).local_element(
1335  edge_constrained_indices[j][i]) = 1.;
1336  }
1337  }
1338 
1339 
1340 
1341  template <int dim, typename VectorType, typename VectorizedArrayType>
1342  void
1343  Base<dim, VectorType, VectorizedArrayType>::vmult(VectorType & dst,
1344  const VectorType &src) const
1345  {
1346  using Number =
1347  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1348  dst = Number(0.);
1349  vmult_add(dst, src);
1350  }
1351 
1352 
1353 
1354  template <int dim, typename VectorType, typename VectorizedArrayType>
1355  void
1356  Base<dim, VectorType, VectorizedArrayType>::vmult_add(
1357  VectorType & dst,
1358  const VectorType &src) const
1359  {
1360  mult_add(dst, src, false);
1361  }
1362 
1363 
1364 
1365  template <int dim, typename VectorType, typename VectorizedArrayType>
1366  void
1367  Base<dim, VectorType, VectorizedArrayType>::Tvmult_add(
1368  VectorType & dst,
1369  const VectorType &src) const
1370  {
1371  mult_add(dst, src, true);
1372  }
1373 
1374 
1375 
1376  template <int dim, typename VectorType, typename VectorizedArrayType>
1377  void
1378  Base<dim, VectorType, VectorizedArrayType>::adjust_ghost_range_if_necessary(
1379  const VectorType &src,
1380  const bool is_row) const
1381  {
1382  using Number =
1383  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1384  for (unsigned int i = 0; i < BlockHelper::n_blocks(src); ++i)
1385  {
1386  const unsigned int mf_component =
1387  is_row ? selected_rows[i] : selected_columns[i];
1388  // If both vectors use the same partitioner -> done
1389  if (BlockHelper::subblock(src, i).get_partitioner().get() ==
1390  data->get_dof_info(mf_component).vector_partitioner.get())
1391  continue;
1392 
1393  // If not, assert that the local ranges are the same and reset to the
1394  // current partitioner
1395  Assert(
1396  BlockHelper::subblock(src, i).get_partitioner()->local_size() ==
1397  data->get_dof_info(mf_component).vector_partitioner->local_size(),
1398  ExcMessage("The vector passed to the vmult() function does not have "
1399  "the correct size for compatibility with MatrixFree."));
1400 
1401  // copy the vector content to a temporary vector so that it does not get
1402  // lost
1403  LinearAlgebra::distributed::Vector<Number> copy_vec(
1404  BlockHelper::subblock(src, i));
1405  BlockHelper::subblock(const_cast<VectorType &>(src), i)
1406  .reinit(data->get_dof_info(mf_component).vector_partitioner);
1407  BlockHelper::subblock(const_cast<VectorType &>(src), i)
1408  .copy_locally_owned_data_from(copy_vec);
1409  }
1410  }
1411 
1412 
1413 
1414  template <int dim, typename VectorType, typename VectorizedArrayType>
1415  void
1416  Base<dim, VectorType, VectorizedArrayType>::preprocess_constraints(
1417  VectorType & dst,
1418  const VectorType &src) const
1419  {
1420  using Number =
1421  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1422  adjust_ghost_range_if_necessary(src, false);
1423  adjust_ghost_range_if_necessary(dst, true);
1424 
1425  // set zero Dirichlet values on the input vector (and remember the src and
1426  // dst values because we need to reset them at the end)
1427  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1428  {
1429  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1430  {
1431  edge_constrained_values[j][i] = std::pair<Number, Number>(
1432  BlockHelper::subblock(src, j).local_element(
1433  edge_constrained_indices[j][i]),
1434  BlockHelper::subblock(dst, j).local_element(
1435  edge_constrained_indices[j][i]));
1436  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1437  .local_element(edge_constrained_indices[j][i]) = 0.;
1438  }
1439  }
1440  }
1441 
1442 
1443 
1444  template <int dim, typename VectorType, typename VectorizedArrayType>
1445  void
1446  Base<dim, VectorType, VectorizedArrayType>::mult_add(
1447  VectorType & dst,
1448  const VectorType &src,
1449  const bool transpose) const
1450  {
1451  AssertDimension(dst.size(), src.size());
1452  AssertDimension(BlockHelper::n_blocks(dst), BlockHelper::n_blocks(src));
1453  AssertDimension(BlockHelper::n_blocks(dst), selected_rows.size());
1454  preprocess_constraints(dst, src);
1455  if (transpose)
1456  Tapply_add(dst, src);
1457  else
1458  apply_add(dst, src);
1459  postprocess_constraints(dst, src);
1460  }
1461 
1462 
1463 
1464  template <int dim, typename VectorType, typename VectorizedArrayType>
1465  void
1466  Base<dim, VectorType, VectorizedArrayType>::postprocess_constraints(
1467  VectorType & dst,
1468  const VectorType &src) const
1469  {
1470  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1471  {
1472  const std::vector<unsigned int> &constrained_dofs =
1473  data->get_constrained_dofs(selected_rows[j]);
1474  for (const auto constrained_dof : constrained_dofs)
1475  BlockHelper::subblock(dst, j).local_element(constrained_dof) +=
1476  BlockHelper::subblock(src, j).local_element(constrained_dof);
1477  }
1478 
1479  // reset edge constrained values, multiply by unit matrix and add into
1480  // destination
1481  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1482  {
1483  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1484  {
1485  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1486  .local_element(edge_constrained_indices[j][i]) =
1487  edge_constrained_values[j][i].first;
1488  BlockHelper::subblock(dst, j).local_element(
1489  edge_constrained_indices[j][i]) =
1490  edge_constrained_values[j][i].second +
1491  edge_constrained_values[j][i].first;
1492  }
1493  }
1494  }
1495 
1496 
1497 
1498  template <int dim, typename VectorType, typename VectorizedArrayType>
1499  void
1500  Base<dim, VectorType, VectorizedArrayType>::vmult_interface_down(
1501  VectorType & dst,
1502  const VectorType &src) const
1503  {
1504  using Number =
1505  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1506  AssertDimension(dst.size(), src.size());
1507  adjust_ghost_range_if_necessary(src, false);
1508  adjust_ghost_range_if_necessary(dst, true);
1509 
1510  dst = Number(0.);
1511 
1512  if (!have_interface_matrices)
1513  return;
1514 
1515  // set zero Dirichlet values on the input vector (and remember the src and
1516  // dst values because we need to reset them at the end)
1517  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1518  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1519  {
1520  edge_constrained_values[j][i] = std::pair<Number, Number>(
1521  BlockHelper::subblock(src, j).local_element(
1522  edge_constrained_indices[j][i]),
1523  BlockHelper::subblock(dst, j).local_element(
1524  edge_constrained_indices[j][i]));
1525  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1526  .local_element(edge_constrained_indices[j][i]) = 0.;
1527  }
1528 
1529  apply_add(dst, src);
1530 
1531  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1532  {
1533  unsigned int c = 0;
1534  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1535  {
1536  for (; c < edge_constrained_indices[j][i]; ++c)
1537  BlockHelper::subblock(dst, j).local_element(c) = 0.;
1538  ++c;
1539 
1540  // reset the src values
1541  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1542  .local_element(edge_constrained_indices[j][i]) =
1543  edge_constrained_values[j][i].first;
1544  }
1545  for (; c < BlockHelper::subblock(dst, j).local_size(); ++c)
1546  BlockHelper::subblock(dst, j).local_element(c) = 0.;
1547  }
1548  }
1549 
1550 
1551 
1552  template <int dim, typename VectorType, typename VectorizedArrayType>
1553  void
1554  Base<dim, VectorType, VectorizedArrayType>::vmult_interface_up(
1555  VectorType & dst,
1556  const VectorType &src) const
1557  {
1558  using Number =
1559  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1560  AssertDimension(dst.size(), src.size());
1561  adjust_ghost_range_if_necessary(src, false);
1562  adjust_ghost_range_if_necessary(dst, true);
1563 
1564  dst = Number(0.);
1565 
1566  if (!have_interface_matrices)
1567  return;
1568 
1569  VectorType src_cpy(src);
1570  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1571  {
1572  unsigned int c = 0;
1573  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1574  {
1575  for (; c < edge_constrained_indices[j][i]; ++c)
1576  BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1577  ++c;
1578  }
1579  for (; c < BlockHelper::subblock(src_cpy, j).local_size(); ++c)
1580  BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1581  }
1582 
1583  apply_add(dst, src_cpy);
1584 
1585  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1586  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1587  BlockHelper::subblock(dst, j).local_element(
1588  edge_constrained_indices[j][i]) = 0.;
1589  }
1590 
1591 
1592 
1593  template <int dim, typename VectorType, typename VectorizedArrayType>
1594  void
1595  Base<dim, VectorType, VectorizedArrayType>::Tvmult(
1596  VectorType & dst,
1597  const VectorType &src) const
1598  {
1599  using Number =
1600  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1601  dst = Number(0.);
1602  Tvmult_add(dst, src);
1603  }
1604 
1605 
1606 
1607  template <int dim, typename VectorType, typename VectorizedArrayType>
1608  std::size_t
1609  Base<dim, VectorType, VectorizedArrayType>::memory_consumption() const
1610  {
1611  return inverse_diagonal_entries.get() != nullptr ?
1612  inverse_diagonal_entries->memory_consumption() :
1613  sizeof(*this);
1614  }
1615 
1616 
1617 
1618  template <int dim, typename VectorType, typename VectorizedArrayType>
1619  std::shared_ptr<const MatrixFree<
1620  dim,
1621  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
1622  VectorizedArrayType>>
1623  Base<dim, VectorType, VectorizedArrayType>::get_matrix_free() const
1624  {
1625  return data;
1626  }
1627 
1628 
1629 
1630  template <int dim, typename VectorType, typename VectorizedArrayType>
1631  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1632  Base<dim, VectorType, VectorizedArrayType>::get_matrix_diagonal_inverse()
1633  const
1634  {
1635  Assert(inverse_diagonal_entries.get() != nullptr &&
1636  inverse_diagonal_entries->m() > 0,
1637  ExcNotInitialized());
1638  return inverse_diagonal_entries;
1639  }
1640 
1641 
1642 
1643  template <int dim, typename VectorType, typename VectorizedArrayType>
1644  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1645  Base<dim, VectorType, VectorizedArrayType>::get_matrix_diagonal() const
1646  {
1647  Assert(diagonal_entries.get() != nullptr && diagonal_entries->m() > 0,
1648  ExcNotInitialized());
1649  return diagonal_entries;
1650  }
1651 
1652 
1653 
1654  template <int dim, typename VectorType, typename VectorizedArrayType>
1655  void
1656  Base<dim, VectorType, VectorizedArrayType>::Tapply_add(
1657  VectorType & dst,
1658  const VectorType &src) const
1659  {
1660  apply_add(dst, src);
1661  }
1662 
1663 
1664 
1665  template <int dim, typename VectorType, typename VectorizedArrayType>
1666  void
1667  Base<dim, VectorType, VectorizedArrayType>::precondition_Jacobi(
1668  VectorType & dst,
1669  const VectorType & src,
1670  const typename Base<dim, VectorType, VectorizedArrayType>::value_type omega)
1671  const
1672  {
1673  Assert(inverse_diagonal_entries.get() && inverse_diagonal_entries->m() > 0,
1674  ExcNotInitialized());
1675  inverse_diagonal_entries->vmult(dst, src);
1676  dst *= omega;
1677  }
1678 
1679 
1680 
1681  //------------------------- MGInterfaceOperator ------------------------------
1682 
1683  template <typename OperatorType>
1684  MGInterfaceOperator<OperatorType>::MGInterfaceOperator()
1685  : Subscriptor()
1686  , mf_base_operator(nullptr)
1687  {}
1688 
1689 
1690 
1691  template <typename OperatorType>
1692  void
1693  MGInterfaceOperator<OperatorType>::clear()
1694  {
1695  mf_base_operator = nullptr;
1696  }
1697 
1698 
1699 
1700  template <typename OperatorType>
1701  void
1702  MGInterfaceOperator<OperatorType>::initialize(const OperatorType &operator_in)
1703  {
1704  mf_base_operator = &operator_in;
1705  }
1706 
1707 
1708 
1709  template <typename OperatorType>
1710  template <typename VectorType>
1711  void
1712  MGInterfaceOperator<OperatorType>::vmult(VectorType & dst,
1713  const VectorType &src) const
1714  {
1715 #ifndef DEAL_II_MSVC
1716  static_assert(
1717  std::is_same<typename VectorType::value_type, value_type>::value,
1718  "The vector type must be based on the same value type as this "
1719  "operator");
1720 #endif
1721 
1722  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1723 
1724  mf_base_operator->vmult_interface_down(dst, src);
1725  }
1726 
1727 
1728 
1729  template <typename OperatorType>
1730  template <typename VectorType>
1731  void
1732  MGInterfaceOperator<OperatorType>::Tvmult(VectorType & dst,
1733  const VectorType &src) const
1734  {
1735 #ifndef DEAL_II_MSVC
1736  static_assert(
1737  std::is_same<typename VectorType::value_type, value_type>::value,
1738  "The vector type must be based on the same value type as this "
1739  "operator");
1740 #endif
1741 
1742  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1743 
1744  mf_base_operator->vmult_interface_up(dst, src);
1745  }
1746 
1747 
1748 
1749  template <typename OperatorType>
1750  template <typename VectorType>
1751  void
1752  MGInterfaceOperator<OperatorType>::initialize_dof_vector(
1753  VectorType &vec) const
1754  {
1755  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1756 
1757  mf_base_operator->initialize_dof_vector(vec);
1758  }
1759 
1760 
1761 
1762  //-----------------------------MassOperator----------------------------------
1763 
1764  template <int dim,
1765  int fe_degree,
1766  int n_q_points_1d,
1767  int n_components,
1768  typename VectorType,
1769  typename VectorizedArrayType>
1770  MassOperator<dim,
1771  fe_degree,
1772  n_q_points_1d,
1773  n_components,
1774  VectorType,
1775  VectorizedArrayType>::MassOperator()
1776  : Base<dim, VectorType, VectorizedArrayType>()
1777  {}
1778 
1779 
1780 
1781  template <int dim,
1782  int fe_degree,
1783  int n_q_points_1d,
1784  int n_components,
1785  typename VectorType,
1786  typename VectorizedArrayType>
1787  void
1788  MassOperator<dim,
1789  fe_degree,
1790  n_q_points_1d,
1791  n_components,
1792  VectorType,
1793  VectorizedArrayType>::compute_diagonal()
1794  {
1795  using Number =
1796  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1797  Assert((Base<dim, VectorType, VectorizedArrayType>::data.get() != nullptr),
1798  ExcNotInitialized());
1799 
1800  this->inverse_diagonal_entries =
1801  std::make_shared<DiagonalMatrix<VectorType>>();
1802  this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1803  VectorType &inverse_diagonal_vector =
1804  this->inverse_diagonal_entries->get_vector();
1805  VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1806  this->initialize_dof_vector(inverse_diagonal_vector);
1807  this->initialize_dof_vector(diagonal_vector);
1808  inverse_diagonal_vector = Number(1.);
1809  apply_add(diagonal_vector, inverse_diagonal_vector);
1810 
1811  this->set_constrained_entries_to_one(diagonal_vector);
1812  inverse_diagonal_vector = diagonal_vector;
1813 
1814  const unsigned int local_size = inverse_diagonal_vector.local_size();
1815  for (unsigned int i = 0; i < local_size; ++i)
1816  inverse_diagonal_vector.local_element(i) =
1817  Number(1.) / inverse_diagonal_vector.local_element(i);
1818 
1819  inverse_diagonal_vector.update_ghost_values();
1820  diagonal_vector.update_ghost_values();
1821  }
1822 
1823 
1824 
1825  template <int dim,
1826  int fe_degree,
1827  int n_q_points_1d,
1828  int n_components,
1829  typename VectorType,
1830  typename VectorizedArrayType>
1831  void
1832  MassOperator<dim,
1833  fe_degree,
1834  n_q_points_1d,
1835  n_components,
1836  VectorType,
1837  VectorizedArrayType>::apply_add(VectorType & dst,
1838  const VectorType &src) const
1839  {
1840  Base<dim, VectorType, VectorizedArrayType>::data->cell_loop(
1841  &MassOperator::local_apply_cell, this, dst, src);
1842  }
1843 
1844 
1845 
1846  template <int dim,
1847  int fe_degree,
1848  int n_q_points_1d,
1849  int n_components,
1850  typename VectorType,
1851  typename VectorizedArrayType>
1852  void
1853  MassOperator<dim,
1854  fe_degree,
1855  n_q_points_1d,
1856  n_components,
1857  VectorType,
1858  VectorizedArrayType>::
1859  local_apply_cell(
1860  const MatrixFree<
1861  dim,
1862  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
1863  VectorizedArrayType> & data,
1864  VectorType & dst,
1865  const VectorType & src,
1866  const std::pair<unsigned int, unsigned int> &cell_range) const
1867  {
1868  using Number =
1869  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1870  FEEvaluation<dim,
1871  fe_degree,
1872  n_q_points_1d,
1873  n_components,
1874  Number,
1875  VectorizedArrayType>
1876  phi(data, this->selected_rows[0]);
1877  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
1878  {
1879  phi.reinit(cell);
1880  phi.read_dof_values(src);
1881  phi.evaluate(true, false, false);
1882  for (unsigned int q = 0; q < phi.n_q_points; ++q)
1883  phi.submit_value(phi.get_value(q), q);
1884  phi.integrate(true, false);
1885  phi.distribute_local_to_global(dst);
1886  }
1887  }
1888 
1889 
1890  //-----------------------------LaplaceOperator----------------------------------
1891 
1892  template <int dim,
1893  int fe_degree,
1894  int n_q_points_1d,
1895  int n_components,
1896  typename VectorType,
1897  typename VectorizedArrayType>
1898  LaplaceOperator<dim,
1899  fe_degree,
1900  n_q_points_1d,
1901  n_components,
1902  VectorType,
1903  VectorizedArrayType>::LaplaceOperator()
1904  : Base<dim, VectorType, VectorizedArrayType>()
1905  {}
1906 
1907 
1908 
1909  template <int dim,
1910  int fe_degree,
1911  int n_q_points_1d,
1912  int n_components,
1913  typename VectorType,
1914  typename VectorizedArrayType>
1915  void
1916  LaplaceOperator<dim,
1917  fe_degree,
1918  n_q_points_1d,
1919  n_components,
1920  VectorType,
1921  VectorizedArrayType>::clear()
1922  {
1923  Base<dim, VectorType, VectorizedArrayType>::clear();
1924  scalar_coefficient.reset();
1925  }
1926 
1927 
1928 
1929  template <int dim,
1930  int fe_degree,
1931  int n_q_points_1d,
1932  int n_components,
1933  typename VectorType,
1934  typename VectorizedArrayType>
1935  void
1936  LaplaceOperator<dim,
1937  fe_degree,
1938  n_q_points_1d,
1939  n_components,
1940  VectorType,
1941  VectorizedArrayType>::
1942  set_coefficient(
1943  const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient_)
1944  {
1945  scalar_coefficient = scalar_coefficient_;
1946  }
1947 
1948 
1949 
1950  template <int dim,
1951  int fe_degree,
1952  int n_q_points_1d,
1953  int n_components,
1954  typename VectorType,
1955  typename VectorizedArrayType>
1956  std::shared_ptr<Table<2, VectorizedArrayType>>
1957  LaplaceOperator<dim,
1958  fe_degree,
1959  n_q_points_1d,
1960  n_components,
1961  VectorType,
1962  VectorizedArrayType>::get_coefficient()
1963  {
1964  Assert(scalar_coefficient.get(), ExcNotInitialized());
1965  return scalar_coefficient;
1966  }
1967 
1968 
1969 
1970  template <int dim,
1971  int fe_degree,
1972  int n_q_points_1d,
1973  int n_components,
1974  typename VectorType,
1975  typename VectorizedArrayType>
1976  void
1977  LaplaceOperator<dim,
1978  fe_degree,
1979  n_q_points_1d,
1980  n_components,
1981  VectorType,
1982  VectorizedArrayType>::compute_diagonal()
1983  {
1984  using Number =
1985  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1986  Assert((Base<dim, VectorType, VectorizedArrayType>::data.get() != nullptr),
1987  ExcNotInitialized());
1988 
1989  this->inverse_diagonal_entries =
1990  std::make_shared<DiagonalMatrix<VectorType>>();
1991  this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1992  VectorType &inverse_diagonal_vector =
1993  this->inverse_diagonal_entries->get_vector();
1994  VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1995  this->initialize_dof_vector(inverse_diagonal_vector);
1996  this->initialize_dof_vector(diagonal_vector);
1997 
1998  this->data->cell_loop(&LaplaceOperator::local_diagonal_cell,
1999  this,
2000  diagonal_vector,
2001  /*unused*/ diagonal_vector);
2002  this->set_constrained_entries_to_one(diagonal_vector);
2003 
2004  inverse_diagonal_vector = diagonal_vector;
2005 
2006  for (unsigned int i = 0; i < inverse_diagonal_vector.local_size(); ++i)
2007  if (std::abs(inverse_diagonal_vector.local_element(i)) >
2008  std::sqrt(std::numeric_limits<Number>::epsilon()))
2009  inverse_diagonal_vector.local_element(i) =
2010  1. / inverse_diagonal_vector.local_element(i);
2011  else
2012  inverse_diagonal_vector.local_element(i) = 1.;
2013 
2014  inverse_diagonal_vector.update_ghost_values();
2015  diagonal_vector.update_ghost_values();
2016  }
2017 
2018 
2019 
2020  template <int dim,
2021  int fe_degree,
2022  int n_q_points_1d,
2023  int n_components,
2024  typename VectorType,
2025  typename VectorizedArrayType>
2026  void
2027  LaplaceOperator<dim,
2028  fe_degree,
2029  n_q_points_1d,
2030  n_components,
2031  VectorType,
2032  VectorizedArrayType>::apply_add(VectorType & dst,
2033  const VectorType &src) const
2034  {
2035  Base<dim, VectorType, VectorizedArrayType>::data->cell_loop(
2036  &LaplaceOperator::local_apply_cell, this, dst, src);
2037  }
2038 
2039  namespace Implementation
2040  {
2041  template <typename VectorizedArrayType>
2042  bool
2043  non_negative(const VectorizedArrayType &n)
2044  {
2045  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2046  if (n[v] < 0.)
2047  return false;
2048 
2049  return true;
2050  }
2051  } // namespace Implementation
2052 
2053 
2054 
2055  template <int dim,
2056  int fe_degree,
2057  int n_q_points_1d,
2058  int n_components,
2059  typename VectorType,
2060  typename VectorizedArrayType>
2061  void
2062  LaplaceOperator<dim,
2063  fe_degree,
2064  n_q_points_1d,
2065  n_components,
2066  VectorType,
2067  VectorizedArrayType>::
2068  do_operation_on_cell(
2069  FEEvaluation<
2070  dim,
2071  fe_degree,
2072  n_q_points_1d,
2073  n_components,
2074  typename Base<dim, VectorType, VectorizedArrayType>::value_type> &phi,
2075  const unsigned int cell) const
2076  {
2077  phi.evaluate(false, true, false);
2078  if (scalar_coefficient.get())
2079  {
2080  Assert(scalar_coefficient->size(1) == 1 ||
2081  scalar_coefficient->size(1) == phi.n_q_points,
2082  ExcMessage("The number of columns in the coefficient table must "
2083  "be either 1 or the number of quadrature points " +
2084  std::to_string(phi.n_q_points) +
2085  ", but the given value was " +
2086  std::to_string(scalar_coefficient->size(1))));
2087  if (scalar_coefficient->size(1) == phi.n_q_points)
2088  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2089  {
2090  Assert(Implementation::non_negative(
2091  (*scalar_coefficient)(cell, q)),
2092  ExcMessage("Coefficient must be non-negative"));
2093  phi.submit_gradient((*scalar_coefficient)(cell, q) *
2094  phi.get_gradient(q),
2095  q);
2096  }
2097  else
2098  {
2099  Assert(Implementation::non_negative((*scalar_coefficient)(cell, 0)),
2100  ExcMessage("Coefficient must be non-negative"));
2101  const VectorizedArrayType coefficient =
2102  (*scalar_coefficient)(cell, 0);
2103  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2104  phi.submit_gradient(coefficient * phi.get_gradient(q), q);
2105  }
2106  }
2107  else
2108  {
2109  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2110  {
2111  phi.submit_gradient(phi.get_gradient(q), q);
2112  }
2113  }
2114  phi.integrate(false, true);
2115  }
2116 
2117 
2118 
2119  template <int dim,
2120  int fe_degree,
2121  int n_q_points_1d,
2122  int n_components,
2123  typename VectorType,
2124  typename VectorizedArrayType>
2125  void
2126  LaplaceOperator<dim,
2127  fe_degree,
2128  n_q_points_1d,
2129  n_components,
2130  VectorType,
2131  VectorizedArrayType>::
2132  local_apply_cell(
2133  const MatrixFree<
2134  dim,
2135  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2136  VectorizedArrayType> & data,
2137  VectorType & dst,
2138  const VectorType & src,
2139  const std::pair<unsigned int, unsigned int> &cell_range) const
2140  {
2141  using Number =
2142  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2143  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number> phi(
2144  data, this->selected_rows[0]);
2145  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2146  {
2147  phi.reinit(cell);
2148  phi.read_dof_values(src);
2149  do_operation_on_cell(phi, cell);
2150  phi.distribute_local_to_global(dst);
2151  }
2152  }
2153 
2154 
2155  template <int dim,
2156  int fe_degree,
2157  int n_q_points_1d,
2158  int n_components,
2159  typename VectorType,
2160  typename VectorizedArrayType>
2161  void
2162  LaplaceOperator<dim,
2163  fe_degree,
2164  n_q_points_1d,
2165  n_components,
2166  VectorType,
2167  VectorizedArrayType>::
2168  local_diagonal_cell(
2169  const MatrixFree<
2170  dim,
2171  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2172  VectorizedArrayType> &data,
2173  VectorType & dst,
2174  const VectorType &,
2175  const std::pair<unsigned int, unsigned int> &cell_range) const
2176  {
2177  using Number =
2178  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2179 
2180  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number> phi(
2181  data, this->selected_rows[0]);
2182  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2183  {
2184  phi.reinit(cell);
2185  VectorizedArrayType local_diagonal_vector[phi.static_dofs_per_cell];
2186  for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2187  {
2188  for (unsigned int j = 0; j < phi.dofs_per_component; ++j)
2189  phi.begin_dof_values()[j] = VectorizedArrayType();
2190  phi.begin_dof_values()[i] = 1.;
2191  do_operation_on_cell(phi, cell);
2192  local_diagonal_vector[i] = phi.begin_dof_values()[i];
2193  }
2194  for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2195  for (unsigned int c = 0; c < phi.n_components; ++c)
2196  phi.begin_dof_values()[i + c * phi.dofs_per_component] =
2197  local_diagonal_vector[i];
2198  phi.distribute_local_to_global(dst);
2199  }
2200  }
2201 
2202 
2203 } // end of namespace MatrixFreeOperators
2204 
2205 
2206 DEAL_II_NAMESPACE_CLOSE
2207 
2208 #endif
MatrixFreeOperators::MGInterfaceOperator::size_type
typename OperatorType::size_type size_type
Definition: operators.h:550
Physics::Elasticity::Kinematics::epsilon
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:476
MatrixFreeOperators::LaplaceOperator::do_operation_on_cell
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:2068
MatrixFreeOperators::MassOperator::MassOperator
MassOperator()
Definition: operators.h:1775
IndexSet
Definition: index_set.h:74
MatrixFreeOperators::Base::compute_diagonal
virtual void compute_diagonal()=0
MatrixFreeOperators::Base::have_interface_matrices
bool have_interface_matrices
Definition: operators.h:477
MatrixFreeOperators::BlockHelper::subblock
std::enable_if<!IsBlockVector< VectorType >::value, const VectorType & >::type subblock(const VectorType &vector, unsigned int)
Definition: operators.h:90
MatrixFreeOperators::MassOperator
Definition: operators.h:747
MatrixFreeOperators::Base::apply_add
virtual void apply_add(VectorType &dst, const VectorType &src) const =0
MatrixFreeOperators::MGInterfaceOperator
Definition: operators.h:539
FEEvaluation::reinit
void reinit(const unsigned int cell_batch_index)
LinearAlgebra::distributed::Vector< double >
vectorization.h
MGConstrainedDoFs
Definition: mg_constrained_dofs.h:45
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
MatrixFreeOperators::Base::initialize_dof_vector
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1160
FEEvaluation::n_components
static constexpr unsigned int n_components
Definition: fe_evaluation.h:2304
MatrixFreeOperators::MGInterfaceOperator::mf_base_operator
SmartPointer< const OperatorType > mf_base_operator
Definition: operators.h:595
MatrixFreeOperators::Base::Tvmult_add
void Tvmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1367
MatrixFreeOperators::MassOperator::size_type
typename Base< dim, VectorType, VectorizedArrayType >::size_type size_type
Definition: operators.h:760
MatrixFreeOperators::MGInterfaceOperator::MGInterfaceOperator
MGInterfaceOperator()
Definition: operators.h:1684
MatrixFreeOperators::MassOperator::local_apply_cell
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:1859
AlignedVector::size
size_type size() const
VectorType
MatrixFreeOperators::Base::get_matrix_free
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > get_matrix_free() const
Definition: operators.h:1623
MatrixFreeOperators::MGInterfaceOperator::vmult
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1712
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
Patterns::Tools::to_string
std::string to_string(const T &t)
Definition: patterns.h:2360
VectorizedArray
Definition: vectorization.h:395
MatrixFreeOperators::LaplaceOperator
Definition: operators.h:815
MatrixFreeOperators::MGInterfaceOperator::value_type
typename OperatorType::value_type value_type
Definition: operators.h:545
MatrixFreeOperators::CellwiseInverseMassMatrix::fe_eval
const FEEvaluationBase< dim, n_components, Number, false, VectorizedArrayType > & fe_eval
Definition: operators.h:726
MatrixFreeOperators::LaplaceOperator::set_coefficient
void set_coefficient(const std::shared_ptr< Table< 2, VectorizedArrayType >> &scalar_coefficient)
Definition: operators.h:1942
MatrixFreeOperators::Base::selected_columns
std::vector< unsigned int > selected_columns
Definition: operators.h:459
MatrixFreeOperators::Base::edge_constrained_indices
std::vector< std::vector< unsigned int > > edge_constrained_indices
Definition: operators.h:465
MatrixFree
Definition: matrix_free.h:117
MatrixFreeOperators::Base::mult_add
void mult_add(VectorType &dst, const VectorType &src, const bool transpose) const
Definition: operators.h:1446
matrix_free.h
MatrixFreeOperators::Base::clear
virtual void clear()
Definition: operators.h:1135
MatrixFreeOperators::MGInterfaceOperator::clear
void clear()
Definition: operators.h:1693
DoFTools::n_components
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
MatrixFreeOperators::Base::~Base
virtual ~Base() override=default
MatrixFreeOperators::BlockHelper::collect_sizes
std::enable_if<!IsBlockVector< VectorType >::value, void >::type collect_sizes(const VectorType &)
Definition: operators.h:104
MatrixFreeOperators
Definition: operators.h:39
MatrixFreeOperators::Base::vmult_interface_up
void vmult_interface_up(VectorType &dst, const VectorType &src) const
Definition: operators.h:1554
la_parallel_vector.h
MatrixFreeOperators::LaplaceOperator::local_apply_cell
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:2132
MatrixFreeOperators::Base::get_matrix_diagonal
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal() const
Definition: operators.h:1645
Table
Definition: table.h:699
MatrixFreeOperators::Base< dim, LinearAlgebra::distributed::Vector< double >, VectorizedArray< typename LinearAlgebra::distributed::Vector< double > ::value_type > >::value_type
typename LinearAlgebra::distributed::Vector< double > ::value_type value_type
Definition: operators.h:197
Subscriptor
Definition: subscriptor.h:62
MatrixFreeOperators::BlockHelper::n_blocks
std::enable_if< IsBlockVector< VectorType >::value, unsigned int >::type n_blocks(const VectorType &vector)
Definition: operators.h:49
MatrixFreeOperators::Base::Base
Base()
Definition: operators.h:1098
level
unsigned int level
Definition: grid_out.cc:4355
MatrixFreeOperators::Base::postprocess_constraints
void postprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1466
MatrixFreeOperators::LaplaceOperator::local_diagonal_cell
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:2168
MatrixFreeOperators::Base::inverse_diagonal_entries
std::shared_ptr< DiagonalMatrix< VectorType > > inverse_diagonal_entries
Definition: operators.h:447
MatrixFreeOperators::Base::vmult
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1343
MatrixFreeOperators::LaplaceOperator::size_type
typename Base< dim, VectorType, VectorizedArrayType >::size_type size_type
Definition: operators.h:828
MatrixFreeOperators::LaplaceOperator::LaplaceOperator
LaplaceOperator()
Definition: operators.h:1903
MatrixFreeOperators::CellwiseInverseMassMatrix::apply
void apply(const AlignedVector< VectorizedArrayType > &inverse_coefficient, const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition: operators.h:1050
MatrixFreeOperators::Base::n
size_type n() const
Definition: operators.h:1121
MatrixFreeOperators::MGInterfaceOperator::initialize
void initialize(const OperatorType &operator_in)
Definition: operators.h:1702
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
double
MatrixFreeOperators::Base::vmult_interface_down
void vmult_interface_down(VectorType &dst, const VectorType &src) const
Definition: operators.h:1500
MatrixFreeOperators::Base::initialize
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:1185
subscriptor.h
MatrixFreeOperators::Base::adjust_ghost_range_if_necessary
void adjust_ghost_range_if_necessary(const VectorType &vec, const bool is_row) const
Definition: operators.h:1378
MatrixFreeOperators::CellwiseInverseMassMatrix::fill_inverse_JxW_values
void fill_inverse_JxW_values(AlignedVector< VectorizedArrayType > &inverse_jxw) const
Definition: operators.h:994
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
MatrixFreeOperators::Base::edge_constrained_values
std::vector< std::vector< std::pair< value_type, value_type > > > edge_constrained_values
Definition: operators.h:471
transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition: derivative_form.h:470
internal::CellwiseInverseMassMatrixImpl
Definition: evaluation_kernels.h:3070
MemoryConsumption::memory_consumption
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Definition: memory_consumption.h:268
MatrixFreeOperators::CellwiseInverseMassMatrix
Definition: operators.h:624
AlignedVector< VectorizedArrayType >
std_cxx17::apply
auto apply(F &&fn, Tuple &&t) -> decltype(apply_impl(std::forward< F >(fn), std::forward< Tuple >(t), std_cxx14::make_index_sequence< std::tuple_size< typename std::remove_reference< Tuple >::type >::value >()))
Definition: tuple.h:40
IndexSet::is_element
bool is_element(const size_type index) const
Definition: index_set.h:1766
MatrixFreeOperators::LaplaceOperator::compute_diagonal
virtual void compute_diagonal() override
Definition: operators.h:1982
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
StandardExceptions::ExcNotInitialized
static ::ExceptionBase & ExcNotInitialized()
MatrixFreeOperators::Base::m
size_type m() const
Definition: operators.h:1107
FEEvaluation::dofs_per_component
const unsigned int dofs_per_component
Definition: fe_evaluation.h:2607
exceptions.h
MatrixFreeOperators::LaplaceOperator::apply_add
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition: operators.h:2032
mg_constrained_dofs.h
SmartPointer< const OperatorType >
MatrixFreeOperators::Base::selected_rows
std::vector< unsigned int > selected_rows
Definition: operators.h:453
value
static const bool value
Definition: dof_tools_constraints.cc:433
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
MatrixFreeOperators::Base::preprocess_constraints
void preprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1416
diagonal_matrix.h
LinearAlgebra::CUDAWrappers::kernel::size_type
types::global_dof_index size_type
Definition: cuda_kernels.h:45
MatrixFreeOperators::CellwiseInverseMassMatrix::CellwiseInverseMassMatrix
CellwiseInverseMassMatrix(const FEEvaluationBase< dim, n_components, Number, false, VectorizedArrayType > &fe_eval)
Definition: operators.h:972
MatrixFreeOperators::Base::set_constrained_entries_to_one
void set_constrained_entries_to_one(VectorType &dst) const
Definition: operators.h:1324
MatrixFreeOperators::LaplaceOperator::value_type
typename Base< dim, VectorType, VectorizedArrayType >::value_type value_type
Definition: operators.h:822
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
MatrixFreeOperators::Implementation::non_negative
bool non_negative(const VectorizedArrayType &n)
Definition: operators.h:2043
MatrixFreeOperators::MGInterfaceOperator::Tvmult
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1732
MatrixFreeOperators::LaplaceOperator::scalar_coefficient
std::shared_ptr< Table< 2, VectorizedArrayType > > scalar_coefficient
Definition: operators.h:955
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:191
MatrixFreeOperators::MGInterfaceOperator::initialize_dof_vector
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1752
MatrixFreeOperators::MassOperator::compute_diagonal
virtual void compute_diagonal() override
Definition: operators.h:1793
MatrixFreeOperators::CellwiseInverseMassMatrix::transform_from_q_points_to_basis
void transform_from_q_points_to_basis(const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition: operators.h:1079
MatrixFreeOperators::Base::get_matrix_diagonal_inverse
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal_inverse() const
Definition: operators.h:1632
MatrixFreeOperators::BlockHelper::collect_sizes
std::enable_if< IsBlockVector< VectorType >::value, void >::type collect_sizes(VectorType &vector)
Definition: operators.h:97
IndexSet::index_within_set
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1922
MatrixFreeOperators::LaplaceOperator::get_coefficient
std::shared_ptr< Table< 2, VectorizedArrayType > > get_coefficient()
Definition: operators.h:1962
MatrixFreeOperators::Base::vmult_add
void vmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1356
config.h
FEEvaluationBase
Definition: fe_evaluation.h:97
FEEvaluation::static_dofs_per_cell
static constexpr unsigned int static_dofs_per_cell
Definition: fe_evaluation.h:2342
MatrixFreeOperators::MassOperator::apply_add
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition: operators.h:1837
FEEvaluation
Definition: fe_evaluation.h:57
MatrixFreeOperators::Base::diagonal_entries
std::shared_ptr< DiagonalMatrix< VectorType > > diagonal_entries
Definition: operators.h:441
MatrixFreeOperators::LaplaceOperator::clear
virtual void clear() override
Definition: operators.h:1921
fe_evaluation.h
MatrixFreeOperators::BlockHelper::n_blocks
std::enable_if<!IsBlockVector< VectorType >::value, unsigned int >::type n_blocks(const VectorType &)
Definition: operators.h:57
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
MatrixFreeOperators::Base::precondition_Jacobi
void precondition_Jacobi(VectorType &dst, const VectorType &src, const value_type omega) const
Definition: operators.h:1667
MatrixFreeOperators::MassOperator::value_type
typename Base< dim, VectorType, VectorizedArrayType >::value_type value_type
Definition: operators.h:754
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
DerivativeForm::transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition: derivative_form.h:470
MatrixFreeOperators::Base::data
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > data
Definition: operators.h:435
MatrixFreeOperators::Base::Tvmult
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1595
MatrixFreeOperators::Base::el
value_type el(const unsigned int row, const unsigned int col) const
Definition: operators.h:1145
MatrixFreeOperators::Base
Definition: operators.h:191
Utilities::MPI::max
T max(const T &t, const MPI_Comm &mpi_communicator)
MatrixFreeOperators::Base::memory_consumption
virtual std::size_t memory_consumption() const
Definition: operators.h:1609
MatrixFreeOperators::BlockHelper::subblock
std::enable_if< IsBlockVector< VectorType >::value, typename VectorType::BlockType & >::type subblock(VectorType &vector, unsigned int block_no)
Definition: operators.h:65
MatrixFreeOperators::Base< dim, LinearAlgebra::distributed::Vector< double >, VectorizedArray< typename LinearAlgebra::distributed::Vector< double > ::value_type > >::size_type
typename LinearAlgebra::distributed::Vector< double > ::size_type size_type
Definition: operators.h:202
MatrixFreeOperators::Base::Tapply_add
virtual void Tapply_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1656
Generated by   doxygen 1.8.17