Reference documentation for deal.II version Git 9c5841c5ee 2019-11-11 14:05:52 -0700
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
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/exceptions.h>
22 #include <deal.II/base/subscriptor.h>
23 #include <deal.II/base/vectorization.h>
24 
25 #include <deal.II/lac/diagonal_matrix.h>
26 #include <deal.II/lac/la_parallel_vector.h>
27 
28 #include <deal.II/matrix_free/fe_evaluation.h>
29 #include <deal.II/matrix_free/matrix_free.h>
30 
31 #include <deal.II/multigrid/mg_constrained_dofs.h>
32 
33 
34 DEAL_II_NAMESPACE_OPEN
35 
36 
37 namespace MatrixFreeOperators
38 {
39  namespace BlockHelper
40  {
41  // workaroud for unifying non-block vector and block vector implementations
42  // a non-block vector has one block and the only subblock is the vector
43  // itself
44  template <typename VectorType>
45  typename std::enable_if<IsBlockVector<VectorType>::value,
46  unsigned int>::type
47  n_blocks(const VectorType &vector)
48  {
49  return vector.n_blocks();
50  }
51 
52  template <typename VectorType>
53  typename std::enable_if<!IsBlockVector<VectorType>::value,
54  unsigned int>::type
55  n_blocks(const VectorType &)
56  {
57  return 1;
58  }
59 
60  template <typename VectorType>
61  typename std::enable_if<IsBlockVector<VectorType>::value,
62  typename VectorType::BlockType &>::type
63  subblock(VectorType &vector, unsigned int block_no)
64  {
65  return vector.block(block_no);
66  }
67 
68  template <typename VectorType>
69  typename std::enable_if<IsBlockVector<VectorType>::value,
70  const typename VectorType::BlockType &>::type
71  subblock(const VectorType &vector, unsigned int block_no)
72  {
73  AssertIndexRange(block_no, vector.n_blocks());
74  return vector.block(block_no);
75  }
76 
77  template <typename VectorType>
78  typename std::enable_if<!IsBlockVector<VectorType>::value,
79  VectorType &>::type
80  subblock(VectorType &vector, unsigned int)
81  {
82  return vector;
83  }
84 
85  template <typename VectorType>
86  typename std::enable_if<!IsBlockVector<VectorType>::value,
87  const VectorType &>::type
88  subblock(const VectorType &vector, unsigned int)
89  {
90  return vector;
91  }
92 
93  template <typename VectorType>
94  typename std::enable_if<IsBlockVector<VectorType>::value, void>::type
95  collect_sizes(VectorType &vector)
96  {
97  vector.collect_sizes();
98  }
99 
100  template <typename VectorType>
101  typename std::enable_if<!IsBlockVector<VectorType>::value, void>::type
102  collect_sizes(const VectorType &)
103  {}
104  } // namespace BlockHelper
105 
185  template <int dim,
186  typename VectorType = LinearAlgebra::distributed::Vector<double>,
187  typename VectorizedArrayType =
189  class Base : public Subscriptor
190  {
191  public:
195  using value_type = typename VectorType::value_type;
196 
200  using size_type = typename VectorType::size_type;
201 
205  Base();
206 
210  virtual ~Base() override = default;
211 
216  virtual void
217  clear();
218 
235  void
236  initialize(std::shared_ptr<
238  const std::vector<unsigned int> &selected_row_blocks =
239  std::vector<unsigned int>(),
240  const std::vector<unsigned int> &selected_column_blocks =
241  std::vector<unsigned int>());
242 
256  void
257  initialize(std::shared_ptr<
259  const MGConstrainedDoFs & mg_constrained_dofs,
260  const unsigned int level,
261  const std::vector<unsigned int> &selected_row_blocks =
262  std::vector<unsigned int>());
263 
278  void
279  initialize(std::shared_ptr<
281  const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
282  const unsigned int level,
283  const std::vector<unsigned int> & selected_row_blocks =
284  std::vector<unsigned int>());
285 
289  size_type
290  m() const;
291 
295  size_type
296  n() const;
297 
301  void
302  vmult_interface_down(VectorType &dst, const VectorType &src) const;
303 
307  void
308  vmult_interface_up(VectorType &dst, const VectorType &src) const;
309 
313  void
314  vmult(VectorType &dst, const VectorType &src) const;
315 
319  void
320  Tvmult(VectorType &dst, const VectorType &src) const;
321 
325  void
326  vmult_add(VectorType &dst, const VectorType &src) const;
327 
331  void
332  Tvmult_add(VectorType &dst, const VectorType &src) const;
333 
338  value_type
339  el(const unsigned int row, const unsigned int col) const;
340 
345  virtual std::size_t
346  memory_consumption() const;
347 
351  void
352  initialize_dof_vector(VectorType &vec) const;
353 
361  virtual void
362  compute_diagonal() = 0;
363 
367  std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
368  get_matrix_free() const;
369 
373  const std::shared_ptr<DiagonalMatrix<VectorType>> &
374  get_matrix_diagonal_inverse() const;
375 
379  const std::shared_ptr<DiagonalMatrix<VectorType>> &
380  get_matrix_diagonal() const;
381 
382 
388  void
389  precondition_Jacobi(VectorType & dst,
390  const VectorType &src,
391  const value_type omega) const;
392 
393  protected:
398  void
399  preprocess_constraints(VectorType &dst, const VectorType &src) const;
400 
405  void
406  postprocess_constraints(VectorType &dst, const VectorType &src) const;
407 
412  void
413  set_constrained_entries_to_one(VectorType &dst) const;
414 
418  virtual void
419  apply_add(VectorType &dst, const VectorType &src) const = 0;
420 
426  virtual void
427  Tapply_add(VectorType &dst, const VectorType &src) const;
428 
432  std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
434 
439  std::shared_ptr<DiagonalMatrix<VectorType>> diagonal_entries;
440 
445  std::shared_ptr<DiagonalMatrix<VectorType>> inverse_diagonal_entries;
446 
451  std::vector<unsigned int> selected_rows;
452 
457  std::vector<unsigned int> selected_columns;
458 
459  private:
463  std::vector<std::vector<unsigned int>> edge_constrained_indices;
464 
468  mutable std::vector<std::vector<std::pair<value_type, value_type>>>
470 
476 
481  void
482  mult_add(VectorType & dst,
483  const VectorType &src,
484  const bool transpose) const;
485 
493  void
494  adjust_ghost_range_if_necessary(const VectorType &vec,
495  const bool is_row) const;
496  };
497 
498 
499 
536  template <typename OperatorType>
538  {
539  public:
543  using value_type = typename OperatorType::value_type;
544 
548  using size_type = typename OperatorType::size_type;
549 
554 
558  void
559  clear();
560 
564  void
565  initialize(const OperatorType &operator_in);
566 
570  template <typename VectorType>
571  void
572  vmult(VectorType &dst, const VectorType &src) const;
573 
577  template <typename VectorType>
578  void
579  Tvmult(VectorType &dst, const VectorType &src) const;
580 
584  template <typename VectorType>
585  void
586  initialize_dof_vector(VectorType &vec) const;
587 
588 
589  private:
594  };
595 
596 
597 
617  template <int dim,
618  int fe_degree,
619  int n_components = 1,
620  typename Number = double,
621  typename VectorizedArrayType = VectorizedArray<Number>>
623  {
624  static_assert(
625  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
626  "Type of Number and of VectorizedArrayType do not match.");
627 
628  public:
634  const FEEvaluationBase<dim,
635  n_components,
636  Number,
637  false,
638  VectorizedArrayType> &fe_eval);
639 
648  void
649  apply(const AlignedVector<VectorizedArrayType> &inverse_coefficient,
650  const unsigned int n_actual_components,
651  const VectorizedArrayType * in_array,
652  VectorizedArrayType * out_array) const;
653 
687  void
688  transform_from_q_points_to_basis(const unsigned int n_actual_components,
689  const VectorizedArrayType *in_array,
690  VectorizedArrayType *out_array) const;
691 
697  void
698  fill_inverse_JxW_values(
699  AlignedVector<VectorizedArrayType> &inverse_jxw) const;
700 
701  private:
705  const FEEvaluationBase<dim,
706  n_components,
707  Number,
708  false,
709  VectorizedArrayType> &fe_eval;
710 
715  };
716 
717 
718 
728  template <int dim,
729  int fe_degree,
730  int n_q_points_1d = fe_degree + 1,
731  int n_components = 1,
732  typename VectorType = LinearAlgebra::distributed::Vector<double>,
733  typename VectorizedArrayType =
735  class MassOperator : public Base<dim, VectorType, VectorizedArrayType>
736  {
737  public:
741  using value_type =
743 
747  using size_type =
749 
753  MassOperator();
754 
759  virtual void
760  compute_diagonal() override;
761 
762  private:
768  virtual void
769  apply_add(VectorType &dst, const VectorType &src) const override;
770 
774  void
775  local_apply_cell(
777  VectorType & dst,
778  const VectorType & src,
779  const std::pair<unsigned int, unsigned int> & cell_range) const;
780  };
781 
782 
783 
796  template <int dim,
797  int fe_degree,
798  int n_q_points_1d = fe_degree + 1,
799  int n_components = 1,
800  typename VectorType = LinearAlgebra::distributed::Vector<double>,
801  typename VectorizedArrayType =
803  class LaplaceOperator : public Base<dim, VectorType, VectorizedArrayType>
804  {
805  public:
809  using value_type =
811 
815  using size_type =
817 
821  LaplaceOperator();
822 
829  virtual void
830  compute_diagonal();
831 
878  void
879  set_coefficient(
880  const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient);
881 
882  virtual void
883  clear();
884 
891  std::shared_ptr<Table<2, VectorizedArrayType>>
892  get_coefficient();
893 
894  private:
900  virtual void
901  apply_add(VectorType &dst, const VectorType &src) const;
902 
906  void
907  local_apply_cell(
909  VectorType & dst,
910  const VectorType & src,
911  const std::pair<unsigned int, unsigned int> & cell_range) const;
912 
916  void
917  local_diagonal_cell(
919  VectorType & dst,
920  const VectorType &,
921  const std::pair<unsigned int, unsigned int> &cell_range) const;
922 
926  void
927  do_operation_on_cell(
929  & phi,
930  const unsigned int cell) const;
931 
935  std::shared_ptr<Table<2, VectorizedArrayType>> scalar_coefficient;
936  };
937 
938 
939 
940  // ------------------------------------ inline functions ---------------------
941 
942  template <int dim,
943  int fe_degree,
944  int n_components,
945  typename Number,
946  typename VectorizedArrayType>
947  inline CellwiseInverseMassMatrix<dim,
948  fe_degree,
949  n_components,
950  Number,
951  VectorizedArrayType>::
952  CellwiseInverseMassMatrix(
953  const FEEvaluationBase<dim,
954  n_components,
955  Number,
956  false,
957  VectorizedArrayType> &fe_eval)
958  : fe_eval(fe_eval)
959  {
960  FullMatrix<double> shapes_1d(fe_degree + 1, fe_degree + 1);
961  for (unsigned int i = 0, c = 0; i < shapes_1d.m(); ++i)
962  for (unsigned int j = 0; j < shapes_1d.n(); ++j, ++c)
963  shapes_1d(i, j) = fe_eval.get_shape_info().shape_values[c][0];
964  shapes_1d.gauss_jordan();
965  const unsigned int stride = (fe_degree + 2) / 2;
966  inverse_shape.resize(stride * (fe_degree + 1));
967  for (unsigned int i = 0; i < stride; ++i)
968  for (unsigned int q = 0; q < (fe_degree + 2) / 2; ++q)
969  {
970  inverse_shape[i * stride + q] =
971  0.5 * (shapes_1d(i, q) + shapes_1d(i, fe_degree - q));
972  inverse_shape[(fe_degree - i) * stride + q] =
973  0.5 * (shapes_1d(i, q) - shapes_1d(i, fe_degree - q));
974  }
975  if (fe_degree % 2 == 0)
976  for (unsigned int q = 0; q < (fe_degree + 2) / 2; ++q)
977  inverse_shape[fe_degree / 2 * stride + q] = shapes_1d(fe_degree / 2, q);
978  }
979 
980 
981 
982  template <int dim,
983  int fe_degree,
984  int n_components,
985  typename Number,
986  typename VectorizedArrayType>
987  inline void
989  fe_degree,
990  n_components,
991  Number,
992  VectorizedArrayType>::
994  AlignedVector<VectorizedArrayType> &inverse_jxw) const
995  {
996  constexpr unsigned int dofs_per_component_on_cell =
997  Utilities::pow(fe_degree + 1, dim);
998  Assert(inverse_jxw.size() > 0 &&
999  inverse_jxw.size() % dofs_per_component_on_cell == 0,
1000  ExcMessage(
1001  "Expected diagonal to be a multiple of scalar dof per cells"));
1002 
1003  // temporarily reduce size of inverse_jxw to dofs_per_cell to get JxW values
1004  // from fe_eval (will not reallocate any memory)
1005  for (unsigned int q = 0; q < dofs_per_component_on_cell; ++q)
1006  inverse_jxw[q] = 1. / fe_eval.JxW(q);
1007  // copy values to rest of vector
1008  for (unsigned int q = dofs_per_component_on_cell; q < inverse_jxw.size();)
1009  for (unsigned int i = 0; i < dofs_per_component_on_cell; ++i, ++q)
1010  inverse_jxw[q] = inverse_jxw[i];
1011  }
1012 
1013 
1014 
1015  template <int dim,
1016  int fe_degree,
1017  int n_components,
1018  typename Number,
1019  typename VectorizedArrayType>
1020  inline void
1022  fe_degree,
1023  n_components,
1024  Number,
1025  VectorizedArrayType>::
1026  apply(const AlignedVector<VectorizedArrayType> &inverse_coefficients,
1027  const unsigned int n_actual_components,
1028  const VectorizedArrayType * in_array,
1029  VectorizedArrayType * out_array) const
1030  {
1031  internal::CellwiseInverseMassMatrixImpl<
1032  dim,
1033  fe_degree,
1034  n_components,
1035  VectorizedArrayType>::apply(inverse_shape,
1036  inverse_coefficients,
1037  n_actual_components,
1038  in_array,
1039  out_array);
1040  }
1041 
1042 
1043 
1044  template <int dim,
1045  int fe_degree,
1046  int n_components,
1047  typename Number,
1048  typename VectorizedArrayType>
1049  inline void
1051  fe_degree,
1052  n_components,
1053  Number,
1054  VectorizedArrayType>::
1055  transform_from_q_points_to_basis(const unsigned int n_actual_components,
1056  const VectorizedArrayType *in_array,
1057  VectorizedArrayType * out_array) const
1058  {
1059  internal::CellwiseInverseMassMatrixImpl<dim,
1060  fe_degree,
1061  n_components,
1062  VectorizedArrayType>::
1064  n_actual_components,
1065  in_array,
1066  out_array);
1067  }
1068 
1069 
1070 
1071  //----------------- Base operator -----------------------------
1072  template <int dim, typename VectorType, typename VectorizedArrayType>
1074  : Subscriptor()
1075  , have_interface_matrices(false)
1076  {}
1077 
1078 
1079 
1080  template <int dim, typename VectorType, typename VectorizedArrayType>
1083  {
1084  Assert(data.get() != nullptr, ExcNotInitialized());
1086  0;
1087  for (unsigned int i = 0; i < selected_rows.size(); ++i)
1088  total_size += data->get_vector_partitioner(selected_rows[i])->size();
1089  return total_size;
1090  }
1091 
1092 
1093 
1094  template <int dim, typename VectorType, typename VectorizedArrayType>
1097  {
1098  Assert(data.get() != nullptr, ExcNotInitialized());
1100  0;
1101  for (unsigned int i = 0; i < selected_columns.size(); ++i)
1102  total_size += data->get_vector_partitioner(selected_columns[i])->size();
1103  return total_size;
1104  }
1105 
1106 
1107 
1108  template <int dim, typename VectorType, typename VectorizedArrayType>
1109  void
1111  {
1112  data.reset();
1113  inverse_diagonal_entries.reset();
1114  }
1115 
1116 
1117 
1118  template <int dim, typename VectorType, typename VectorizedArrayType>
1121  const unsigned int col) const
1122  {
1123  (void)col;
1124  Assert(row == col, ExcNotImplemented());
1125  Assert(inverse_diagonal_entries.get() != nullptr &&
1126  inverse_diagonal_entries->m() > 0,
1127  ExcNotInitialized());
1128  return 1.0 / (*inverse_diagonal_entries)(row, row);
1129  }
1130 
1131 
1132 
1133  template <int dim, typename VectorType, typename VectorizedArrayType>
1134  void
1136  VectorType &vec) const
1137  {
1138  Assert(data.get() != nullptr, ExcNotInitialized());
1139  AssertDimension(BlockHelper::n_blocks(vec), selected_rows.size());
1140  for (unsigned int i = 0; i < BlockHelper::n_blocks(vec); ++i)
1141  {
1142  const unsigned int index = selected_rows[i];
1143  if (!BlockHelper::subblock(vec, index)
1144  .partitioners_are_compatible(
1145  *data->get_dof_info(index).vector_partitioner))
1146  data->initialize_dof_vector(BlockHelper::subblock(vec, index), index);
1147 
1148  Assert(BlockHelper::subblock(vec, index)
1149  .partitioners_are_globally_compatible(
1150  *data->get_dof_info(index).vector_partitioner),
1151  ExcInternalError());
1152  }
1153  BlockHelper::collect_sizes(vec);
1154  }
1155 
1156 
1157 
1158  template <int dim, typename VectorType, typename VectorizedArrayType>
1159  void
1162  data_,
1163  const std::vector<unsigned int> &given_row_selection,
1164  const std::vector<unsigned int> &given_column_selection)
1165  {
1166  data = data_;
1167 
1168  selected_rows.clear();
1169  selected_columns.clear();
1170  if (given_row_selection.empty())
1171  for (unsigned int i = 0; i < data_->n_components(); ++i)
1172  selected_rows.push_back(i);
1173  else
1174  {
1175  for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1176  {
1177  AssertIndexRange(given_row_selection[i], data_->n_components());
1178  for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1179  if (j != i)
1180  Assert(given_row_selection[j] != given_row_selection[i],
1181  ExcMessage("Given row indices must be unique"));
1182 
1183  selected_rows.push_back(given_row_selection[i]);
1184  }
1185  }
1186  if (given_column_selection.size() == 0)
1188  else
1189  {
1190  for (unsigned int i = 0; i < given_column_selection.size(); ++i)
1191  {
1192  AssertIndexRange(given_column_selection[i], data_->n_components());
1193  for (unsigned int j = 0; j < given_column_selection.size(); ++j)
1194  if (j != i)
1195  Assert(given_column_selection[j] != given_column_selection[i],
1196  ExcMessage("Given column indices must be unique"));
1197 
1198  selected_columns.push_back(given_column_selection[i]);
1199  }
1200  }
1201 
1202  edge_constrained_indices.clear();
1203  edge_constrained_indices.resize(selected_rows.size());
1204  edge_constrained_values.clear();
1205  edge_constrained_values.resize(selected_rows.size());
1206  have_interface_matrices = false;
1207  }
1208 
1209 
1210 
1211  template <int dim, typename VectorType, typename VectorizedArrayType>
1212  void
1215  data_,
1216  const MGConstrainedDoFs & mg_constrained_dofs,
1217  const unsigned int level,
1218  const std::vector<unsigned int> &given_row_selection)
1219  {
1220  std::vector<MGConstrainedDoFs> mg_constrained_dofs_vector(
1221  1, mg_constrained_dofs);
1222  initialize(data_, mg_constrained_dofs_vector, level, given_row_selection);
1223  }
1224 
1225 
1226 
1227  template <int dim, typename VectorType, typename VectorizedArrayType>
1228  void
1231  data_,
1232  const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
1233  const unsigned int level,
1234  const std::vector<unsigned int> & given_row_selection)
1235  {
1237  ExcMessage("level is not set"));
1238 
1239  selected_rows.clear();
1240  selected_columns.clear();
1241  if (given_row_selection.empty())
1242  for (unsigned int i = 0; i < data_->n_components(); ++i)
1243  selected_rows.push_back(i);
1244  else
1245  {
1246  for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1247  {
1248  AssertIndexRange(given_row_selection[i], data_->n_components());
1249  for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1250  if (j != i)
1251  Assert(given_row_selection[j] != given_row_selection[i],
1252  ExcMessage("Given row indices must be unique"));
1253 
1254  selected_rows.push_back(given_row_selection[i]);
1255  }
1256  }
1258 
1259  AssertDimension(mg_constrained_dofs.size(), selected_rows.size());
1260  edge_constrained_indices.clear();
1261  edge_constrained_indices.resize(selected_rows.size());
1262  edge_constrained_values.clear();
1263  edge_constrained_values.resize(selected_rows.size());
1264 
1265  data = data_;
1266 
1267  for (unsigned int j = 0; j < selected_rows.size(); ++j)
1268  {
1269  if (data_->n_macro_cells() > 0)
1270  {
1271  AssertDimension(level, data_->get_cell_iterator(0, 0, j)->level());
1272  }
1273 
1274  // setup edge_constrained indices
1275  std::vector<types::global_dof_index> interface_indices;
1276  mg_constrained_dofs[j]
1277  .get_refinement_edge_indices(level)
1278  .fill_index_vector(interface_indices);
1279  edge_constrained_indices[j].clear();
1280  edge_constrained_indices[j].reserve(interface_indices.size());
1281  edge_constrained_values[j].resize(interface_indices.size());
1282  const IndexSet &locally_owned =
1283  data->get_dof_handler(selected_rows[j]).locally_owned_mg_dofs(level);
1284  for (const auto interface_index : interface_indices)
1285  if (locally_owned.is_element(interface_index))
1286  edge_constrained_indices[j].push_back(
1287  locally_owned.index_within_set(interface_index));
1290  static_cast<unsigned int>(edge_constrained_indices[j].size()),
1291  data->get_vector_partitioner()->get_mpi_communicator()) > 0;
1292  }
1293  }
1294 
1295 
1296 
1297  template <int dim, typename VectorType, typename VectorizedArrayType>
1298  void
1300  VectorType &dst) const
1301  {
1302  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1303  {
1304  const std::vector<unsigned int> &constrained_dofs =
1305  data->get_constrained_dofs(selected_rows[j]);
1306  for (const auto constrained_dof : constrained_dofs)
1307  BlockHelper::subblock(dst, j).local_element(constrained_dof) = 1.;
1308  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1309  BlockHelper::subblock(dst, j).local_element(
1310  edge_constrained_indices[j][i]) = 1.;
1311  }
1312  }
1313 
1314 
1315 
1316  template <int dim, typename VectorType, typename VectorizedArrayType>
1317  void
1319  const VectorType &src) const
1320  {
1321  using Number =
1323  dst = Number(0.);
1324  vmult_add(dst, src);
1325  }
1326 
1327 
1328 
1329  template <int dim, typename VectorType, typename VectorizedArrayType>
1330  void
1332  VectorType & dst,
1333  const VectorType &src) const
1334  {
1335  mult_add(dst, src, false);
1336  }
1337 
1338 
1339 
1340  template <int dim, typename VectorType, typename VectorizedArrayType>
1341  void
1343  VectorType & dst,
1344  const VectorType &src) const
1345  {
1346  mult_add(dst, src, true);
1347  }
1348 
1349 
1350 
1351  template <int dim, typename VectorType, typename VectorizedArrayType>
1352  void
1354  const VectorType &src,
1355  const bool is_row) const
1356  {
1357  using Number =
1359  for (unsigned int i = 0; i < BlockHelper::n_blocks(src); ++i)
1360  {
1361  const unsigned int mf_component =
1362  is_row ? selected_rows[i] : selected_columns[i];
1363  // If both vectors use the same partitioner -> done
1364  if (BlockHelper::subblock(src, i).get_partitioner().get() ==
1365  data->get_dof_info(mf_component).vector_partitioner.get())
1366  continue;
1367 
1368  // If not, assert that the local ranges are the same and reset to the
1369  // current partitioner
1370  Assert(
1371  BlockHelper::subblock(src, i).get_partitioner()->local_size() ==
1372  data->get_dof_info(mf_component).vector_partitioner->local_size(),
1373  ExcMessage("The vector passed to the vmult() function does not have "
1374  "the correct size for compatibility with MatrixFree."));
1375 
1376  // copy the vector content to a temporary vector so that it does not get
1377  // lost
1379  BlockHelper::subblock(src, i));
1380  BlockHelper::subblock(const_cast<VectorType &>(src), i)
1381  .reinit(data->get_dof_info(mf_component).vector_partitioner);
1382  BlockHelper::subblock(const_cast<VectorType &>(src), i)
1383  .copy_locally_owned_data_from(copy_vec);
1384  }
1385  }
1386 
1387 
1388 
1389  template <int dim, typename VectorType, typename VectorizedArrayType>
1390  void
1392  VectorType & dst,
1393  const VectorType &src) const
1394  {
1395  using Number =
1397  adjust_ghost_range_if_necessary(src, false);
1399 
1400  // set zero Dirichlet values on the input vector (and remember the src and
1401  // dst values because we need to reset them at the end)
1402  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1403  {
1404  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1405  {
1406  edge_constrained_values[j][i] = std::pair<Number, Number>(
1407  BlockHelper::subblock(src, j).local_element(
1408  edge_constrained_indices[j][i]),
1409  BlockHelper::subblock(dst, j).local_element(
1410  edge_constrained_indices[j][i]));
1411  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1412  .local_element(edge_constrained_indices[j][i]) = 0.;
1413  }
1414  }
1415  }
1416 
1417 
1418 
1419  template <int dim, typename VectorType, typename VectorizedArrayType>
1420  void
1422  VectorType & dst,
1423  const VectorType &src,
1424  const bool transpose) const
1425  {
1426  AssertDimension(dst.size(), src.size());
1427  AssertDimension(BlockHelper::n_blocks(dst), BlockHelper::n_blocks(src));
1428  AssertDimension(BlockHelper::n_blocks(dst), selected_rows.size());
1429  preprocess_constraints(dst, src);
1430  if (transpose)
1431  Tapply_add(dst, src);
1432  else
1433  apply_add(dst, src);
1434  postprocess_constraints(dst, src);
1435  }
1436 
1437 
1438 
1439  template <int dim, typename VectorType, typename VectorizedArrayType>
1440  void
1442  VectorType & dst,
1443  const VectorType &src) const
1444  {
1445  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1446  {
1447  const std::vector<unsigned int> &constrained_dofs =
1448  data->get_constrained_dofs(selected_rows[j]);
1449  for (const auto constrained_dof : constrained_dofs)
1450  BlockHelper::subblock(dst, j).local_element(constrained_dof) +=
1451  BlockHelper::subblock(src, j).local_element(constrained_dof);
1452  }
1453 
1454  // reset edge constrained values, multiply by unit matrix and add into
1455  // destination
1456  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1457  {
1458  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1459  {
1460  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1461  .local_element(edge_constrained_indices[j][i]) =
1462  edge_constrained_values[j][i].first;
1463  BlockHelper::subblock(dst, j).local_element(
1464  edge_constrained_indices[j][i]) =
1465  edge_constrained_values[j][i].second +
1466  edge_constrained_values[j][i].first;
1467  }
1468  }
1469  }
1470 
1471 
1472 
1473  template <int dim, typename VectorType, typename VectorizedArrayType>
1474  void
1476  VectorType & dst,
1477  const VectorType &src) const
1478  {
1479  using Number =
1481  AssertDimension(dst.size(), src.size());
1482  adjust_ghost_range_if_necessary(src, false);
1484 
1485  dst = Number(0.);
1486 
1488  return;
1489 
1490  // set zero Dirichlet values on the input vector (and remember the src and
1491  // dst values because we need to reset them at the end)
1492  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1493  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1494  {
1495  edge_constrained_values[j][i] = std::pair<Number, Number>(
1496  BlockHelper::subblock(src, j).local_element(
1497  edge_constrained_indices[j][i]),
1498  BlockHelper::subblock(dst, j).local_element(
1499  edge_constrained_indices[j][i]));
1500  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1501  .local_element(edge_constrained_indices[j][i]) = 0.;
1502  }
1503 
1504  apply_add(dst, src);
1505 
1506  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1507  {
1508  unsigned int c = 0;
1509  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1510  {
1511  for (; c < edge_constrained_indices[j][i]; ++c)
1512  BlockHelper::subblock(dst, j).local_element(c) = 0.;
1513  ++c;
1514 
1515  // reset the src values
1516  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1517  .local_element(edge_constrained_indices[j][i]) =
1518  edge_constrained_values[j][i].first;
1519  }
1520  for (; c < BlockHelper::subblock(dst, j).local_size(); ++c)
1521  BlockHelper::subblock(dst, j).local_element(c) = 0.;
1522  }
1523  }
1524 
1525 
1526 
1527  template <int dim, typename VectorType, typename VectorizedArrayType>
1528  void
1530  VectorType & dst,
1531  const VectorType &src) const
1532  {
1533  using Number =
1535  AssertDimension(dst.size(), src.size());
1536  adjust_ghost_range_if_necessary(src, false);
1538 
1539  dst = Number(0.);
1540 
1542  return;
1543 
1544  VectorType src_cpy(src);
1545  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1546  {
1547  unsigned int c = 0;
1548  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1549  {
1550  for (; c < edge_constrained_indices[j][i]; ++c)
1551  BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1552  ++c;
1553  }
1554  for (; c < BlockHelper::subblock(src_cpy, j).local_size(); ++c)
1555  BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1556  }
1557 
1558  apply_add(dst, src_cpy);
1559 
1560  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1561  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1562  BlockHelper::subblock(dst, j).local_element(
1563  edge_constrained_indices[j][i]) = 0.;
1564  }
1565 
1566 
1567 
1568  template <int dim, typename VectorType, typename VectorizedArrayType>
1569  void
1571  VectorType & dst,
1572  const VectorType &src) const
1573  {
1574  using Number =
1576  dst = Number(0.);
1577  Tvmult_add(dst, src);
1578  }
1579 
1580 
1581 
1582  template <int dim, typename VectorType, typename VectorizedArrayType>
1583  std::size_t
1585  {
1586  return inverse_diagonal_entries.get() != nullptr ?
1587  inverse_diagonal_entries->memory_consumption() :
1588  sizeof(*this);
1589  }
1590 
1591 
1592 
1593  template <int dim, typename VectorType, typename VectorizedArrayType>
1594  std::shared_ptr<const MatrixFree<
1595  dim,
1597  VectorizedArrayType>>
1599  {
1600  return data;
1601  }
1602 
1603 
1604 
1605  template <int dim, typename VectorType, typename VectorizedArrayType>
1606  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1608  const
1609  {
1610  Assert(inverse_diagonal_entries.get() != nullptr &&
1611  inverse_diagonal_entries->m() > 0,
1612  ExcNotInitialized());
1613  return inverse_diagonal_entries;
1614  }
1615 
1616 
1617 
1618  template <int dim, typename VectorType, typename VectorizedArrayType>
1619  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1621  {
1622  Assert(diagonal_entries.get() != nullptr && diagonal_entries->m() > 0,
1623  ExcNotInitialized());
1624  return diagonal_entries;
1625  }
1626 
1627 
1628 
1629  template <int dim, typename VectorType, typename VectorizedArrayType>
1630  void
1632  VectorType & dst,
1633  const VectorType &src) const
1634  {
1635  apply_add(dst, src);
1636  }
1637 
1638 
1639 
1640  template <int dim, typename VectorType, typename VectorizedArrayType>
1641  void
1643  VectorType & dst,
1644  const VectorType & src,
1645  const typename Base<dim, VectorType, VectorizedArrayType>::value_type omega)
1646  const
1647  {
1649  ExcNotInitialized());
1650  inverse_diagonal_entries->vmult(dst, src);
1651  dst *= omega;
1652  }
1653 
1654 
1655 
1656  //------------------------- MGInterfaceOperator ------------------------------
1657 
1658  template <typename OperatorType>
1660  : Subscriptor()
1661  , mf_base_operator(nullptr)
1662  {}
1663 
1664 
1665 
1666  template <typename OperatorType>
1667  void
1669  {
1670  mf_base_operator = nullptr;
1671  }
1672 
1673 
1674 
1675  template <typename OperatorType>
1676  void
1677  MGInterfaceOperator<OperatorType>::initialize(const OperatorType &operator_in)
1678  {
1679  mf_base_operator = &operator_in;
1680  }
1681 
1682 
1683 
1684  template <typename OperatorType>
1685  template <typename VectorType>
1686  void
1688  const VectorType &src) const
1689  {
1690 #ifndef DEAL_II_MSVC
1691  static_assert(
1692  std::is_same<typename VectorType::value_type, value_type>::value,
1693  "The vector type must be based on the same value type as this "
1694  "operator");
1695 #endif
1696 
1697  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1698 
1699  mf_base_operator->vmult_interface_down(dst, src);
1700  }
1701 
1702 
1703 
1704  template <typename OperatorType>
1705  template <typename VectorType>
1706  void
1708  const VectorType &src) const
1709  {
1710 #ifndef DEAL_II_MSVC
1711  static_assert(
1712  std::is_same<typename VectorType::value_type, value_type>::value,
1713  "The vector type must be based on the same value type as this "
1714  "operator");
1715 #endif
1716 
1717  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1718 
1719  mf_base_operator->vmult_interface_up(dst, src);
1720  }
1721 
1722 
1723 
1724  template <typename OperatorType>
1725  template <typename VectorType>
1726  void
1728  VectorType &vec) const
1729  {
1730  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1731 
1732  mf_base_operator->initialize_dof_vector(vec);
1733  }
1734 
1735 
1736 
1737  //-----------------------------MassOperator----------------------------------
1738 
1739  template <int dim,
1740  int fe_degree,
1741  int n_q_points_1d,
1742  int n_components,
1743  typename VectorType,
1744  typename VectorizedArrayType>
1745  MassOperator<dim,
1746  fe_degree,
1747  n_q_points_1d,
1748  n_components,
1749  VectorType,
1750  VectorizedArrayType>::MassOperator()
1751  : Base<dim, VectorType, VectorizedArrayType>()
1752  {}
1753 
1754 
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  void
1763  MassOperator<dim,
1764  fe_degree,
1765  n_q_points_1d,
1766  n_components,
1767  VectorType,
1768  VectorizedArrayType>::compute_diagonal()
1769  {
1770  using Number =
1771  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1773  ExcNotInitialized());
1774 
1776  this->diagonal_entries.reset(new DiagonalMatrix<VectorType>());
1777  VectorType &inverse_diagonal_vector =
1778  this->inverse_diagonal_entries->get_vector();
1779  VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1780  this->initialize_dof_vector(inverse_diagonal_vector);
1781  this->initialize_dof_vector(diagonal_vector);
1782  inverse_diagonal_vector = Number(1.);
1783  apply_add(diagonal_vector, inverse_diagonal_vector);
1784 
1785  this->set_constrained_entries_to_one(diagonal_vector);
1786  inverse_diagonal_vector = diagonal_vector;
1787 
1788  const unsigned int local_size = inverse_diagonal_vector.local_size();
1789  for (unsigned int i = 0; i < local_size; ++i)
1790  inverse_diagonal_vector.local_element(i) =
1791  Number(1.) / inverse_diagonal_vector.local_element(i);
1792 
1793  inverse_diagonal_vector.update_ghost_values();
1794  diagonal_vector.update_ghost_values();
1795  }
1796 
1797 
1798 
1799  template <int dim,
1800  int fe_degree,
1801  int n_q_points_1d,
1802  int n_components,
1803  typename VectorType,
1804  typename VectorizedArrayType>
1805  void
1806  MassOperator<dim,
1807  fe_degree,
1808  n_q_points_1d,
1809  n_components,
1810  VectorType,
1811  VectorizedArrayType>::apply_add(VectorType & dst,
1812  const VectorType &src) const
1813  {
1815  &MassOperator::local_apply_cell, this, dst, src);
1816  }
1817 
1818 
1819 
1820  template <int dim,
1821  int fe_degree,
1822  int n_q_points_1d,
1823  int n_components,
1824  typename VectorType,
1825  typename VectorizedArrayType>
1826  void
1827  MassOperator<dim,
1828  fe_degree,
1829  n_q_points_1d,
1830  n_components,
1831  VectorType,
1832  VectorizedArrayType>::
1834  const MatrixFree<
1835  dim,
1836  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
1837  VectorizedArrayType> & data,
1838  VectorType & dst,
1839  const VectorType & src,
1840  const std::pair<unsigned int, unsigned int> &cell_range) const
1841  {
1842  using Number =
1843  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1844  FEEvaluation<dim,
1845  fe_degree,
1846  n_q_points_1d,
1847  n_components,
1848  Number,
1849  VectorizedArrayType>
1850  phi(data, this->selected_rows[0]);
1851  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
1852  {
1853  phi.reinit(cell);
1854  phi.read_dof_values(src);
1855  phi.evaluate(true, false, false);
1856  for (unsigned int q = 0; q < phi.n_q_points; ++q)
1857  phi.submit_value(phi.get_value(q), q);
1858  phi.integrate(true, false);
1859  phi.distribute_local_to_global(dst);
1860  }
1861  }
1862 
1863 
1864  //-----------------------------LaplaceOperator----------------------------------
1865 
1866  template <int dim,
1867  int fe_degree,
1868  int n_q_points_1d,
1869  int n_components,
1870  typename VectorType,
1871  typename VectorizedArrayType>
1872  LaplaceOperator<dim,
1873  fe_degree,
1874  n_q_points_1d,
1875  n_components,
1876  VectorType,
1877  VectorizedArrayType>::LaplaceOperator()
1878  : Base<dim, VectorType, VectorizedArrayType>()
1879  {}
1880 
1881 
1882 
1883  template <int dim,
1884  int fe_degree,
1885  int n_q_points_1d,
1886  int n_components,
1887  typename VectorType,
1888  typename VectorizedArrayType>
1889  void
1890  LaplaceOperator<dim,
1891  fe_degree,
1892  n_q_points_1d,
1893  n_components,
1894  VectorType,
1895  VectorizedArrayType>::clear()
1896  {
1898  scalar_coefficient.reset();
1899  }
1900 
1901 
1902 
1903  template <int dim,
1904  int fe_degree,
1905  int n_q_points_1d,
1906  int n_components,
1907  typename VectorType,
1908  typename VectorizedArrayType>
1909  void
1910  LaplaceOperator<dim,
1911  fe_degree,
1912  n_q_points_1d,
1913  n_components,
1914  VectorType,
1915  VectorizedArrayType>::
1917  const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient_)
1918  {
1919  scalar_coefficient = scalar_coefficient_;
1920  }
1921 
1922 
1923 
1924  template <int dim,
1925  int fe_degree,
1926  int n_q_points_1d,
1927  int n_components,
1928  typename VectorType,
1929  typename VectorizedArrayType>
1930  std::shared_ptr<Table<2, VectorizedArrayType>>
1931  LaplaceOperator<dim,
1932  fe_degree,
1933  n_q_points_1d,
1934  n_components,
1935  VectorType,
1936  VectorizedArrayType>::get_coefficient()
1937  {
1939  return scalar_coefficient;
1940  }
1941 
1942 
1943 
1944  template <int dim,
1945  int fe_degree,
1946  int n_q_points_1d,
1947  int n_components,
1948  typename VectorType,
1949  typename VectorizedArrayType>
1950  void
1951  LaplaceOperator<dim,
1952  fe_degree,
1953  n_q_points_1d,
1954  n_components,
1955  VectorType,
1956  VectorizedArrayType>::compute_diagonal()
1957  {
1958  using Number =
1959  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1961  ExcNotInitialized());
1962 
1964  this->diagonal_entries.reset(new DiagonalMatrix<VectorType>());
1965  VectorType &inverse_diagonal_vector =
1966  this->inverse_diagonal_entries->get_vector();
1967  VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1968  this->initialize_dof_vector(inverse_diagonal_vector);
1969  this->initialize_dof_vector(diagonal_vector);
1970 
1971  this->data->cell_loop(&LaplaceOperator::local_diagonal_cell,
1972  this,
1973  diagonal_vector,
1974  /*unused*/ diagonal_vector);
1975  this->set_constrained_entries_to_one(diagonal_vector);
1976 
1977  inverse_diagonal_vector = diagonal_vector;
1978 
1979  for (unsigned int i = 0; i < inverse_diagonal_vector.local_size(); ++i)
1980  if (std::abs(inverse_diagonal_vector.local_element(i)) >
1981  std::sqrt(std::numeric_limits<Number>::epsilon()))
1982  inverse_diagonal_vector.local_element(i) =
1983  1. / inverse_diagonal_vector.local_element(i);
1984  else
1985  inverse_diagonal_vector.local_element(i) = 1.;
1986 
1987  inverse_diagonal_vector.update_ghost_values();
1988  diagonal_vector.update_ghost_values();
1989  }
1990 
1991 
1992 
1993  template <int dim,
1994  int fe_degree,
1995  int n_q_points_1d,
1996  int n_components,
1997  typename VectorType,
1998  typename VectorizedArrayType>
1999  void
2000  LaplaceOperator<dim,
2001  fe_degree,
2002  n_q_points_1d,
2003  n_components,
2004  VectorType,
2005  VectorizedArrayType>::apply_add(VectorType & dst,
2006  const VectorType &src) const
2007  {
2009  &LaplaceOperator::local_apply_cell, this, dst, src);
2010  }
2011 
2012  namespace Implementation
2013  {
2014  template <typename VectorizedArrayType>
2015  bool
2016  non_negative(const VectorizedArrayType &n)
2017  {
2018  for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
2019  if (n[v] < 0.)
2020  return false;
2021 
2022  return true;
2023  }
2024  } // namespace Implementation
2025 
2026 
2027 
2028  template <int dim,
2029  int fe_degree,
2030  int n_q_points_1d,
2031  int n_components,
2032  typename VectorType,
2033  typename VectorizedArrayType>
2034  void
2035  LaplaceOperator<dim,
2036  fe_degree,
2037  n_q_points_1d,
2038  n_components,
2039  VectorType,
2040  VectorizedArrayType>::
2042  FEEvaluation<
2043  dim,
2044  fe_degree,
2045  n_q_points_1d,
2046  n_components,
2047  typename Base<dim, VectorType, VectorizedArrayType>::value_type> &phi,
2048  const unsigned int cell) const
2049  {
2050  phi.evaluate(false, true, false);
2051  if (scalar_coefficient.get())
2052  {
2053  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2054  {
2055  Assert(Implementation::non_negative((*scalar_coefficient)(cell, q)),
2056  ExcMessage("Coefficient must be non-negative"));
2057  phi.submit_gradient((*scalar_coefficient)(cell, q) *
2058  phi.get_gradient(q),
2059  q);
2060  }
2061  }
2062  else
2063  {
2064  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2065  {
2066  phi.submit_gradient(phi.get_gradient(q), q);
2067  }
2068  }
2069  phi.integrate(false, true);
2070  }
2071 
2072 
2073 
2074  template <int dim,
2075  int fe_degree,
2076  int n_q_points_1d,
2077  int n_components,
2078  typename VectorType,
2079  typename VectorizedArrayType>
2080  void
2081  LaplaceOperator<dim,
2082  fe_degree,
2083  n_q_points_1d,
2084  n_components,
2085  VectorType,
2086  VectorizedArrayType>::
2088  const MatrixFree<
2089  dim,
2090  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2091  VectorizedArrayType> & data,
2092  VectorType & dst,
2093  const VectorType & src,
2094  const std::pair<unsigned int, unsigned int> &cell_range) const
2095  {
2096  using Number =
2097  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2099  data, this->selected_rows[0]);
2100  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2101  {
2102  phi.reinit(cell);
2103  phi.read_dof_values(src);
2104  do_operation_on_cell(phi, cell);
2105  phi.distribute_local_to_global(dst);
2106  }
2107  }
2108 
2109 
2110  template <int dim,
2111  int fe_degree,
2112  int n_q_points_1d,
2113  int n_components,
2114  typename VectorType,
2115  typename VectorizedArrayType>
2116  void
2117  LaplaceOperator<dim,
2118  fe_degree,
2119  n_q_points_1d,
2120  n_components,
2121  VectorType,
2122  VectorizedArrayType>::
2124  const MatrixFree<
2125  dim,
2126  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2127  VectorizedArrayType> &data,
2128  VectorType & dst,
2129  const VectorType &,
2130  const std::pair<unsigned int, unsigned int> &cell_range) const
2131  {
2132  using Number =
2133  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2134 
2136  data, this->selected_rows[0]);
2137  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2138  {
2139  phi.reinit(cell);
2140  VectorizedArrayType local_diagonal_vector[phi.static_dofs_per_cell];
2141  for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2142  {
2143  for (unsigned int j = 0; j < phi.dofs_per_component; ++j)
2144  phi.begin_dof_values()[j] = VectorizedArrayType();
2145  phi.begin_dof_values()[i] = 1.;
2146  do_operation_on_cell(phi, cell);
2147  local_diagonal_vector[i] = phi.begin_dof_values()[i];
2148  }
2149  for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2150  for (unsigned int c = 0; c < phi.n_components; ++c)
2151  phi.begin_dof_values()[i + c * phi.dofs_per_component] =
2152  local_diagonal_vector[i];
2153  phi.distribute_local_to_global(dst);
2154  }
2155  }
2156 
2157 
2158 } // end of namespace MatrixFreeOperators
2159 
2160 
2161 DEAL_II_NAMESPACE_CLOSE
2162 
2163 #endif
typename OperatorType::value_type value_type
Definition: operators.h:543
VectorizedArrayType JxW(const unsigned int q_index) const
size_type m() 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:993
static const unsigned int invalid_unsigned_int
Definition: types.h:187
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal_inverse() const
Definition: operators.h:1607
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1571
virtual std::size_t memory_consumption() const
Definition: operators.h:1584
const FEEvaluationBase< dim, n_components, Number, false, VectorizedArrayType > & fe_eval
Definition: operators.h:709
void gauss_jordan()
size_type n() const
Definition: operators.h:1096
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:2123
void precondition_Jacobi(VectorType &dst, const VectorType &src, const value_type omega) const
Definition: operators.h:1642
std::shared_ptr< DiagonalMatrix< VectorType > > diagonal_entries
Definition: operators.h:439
void read_dof_values(const VectorType &src, const unsigned int first_index=0)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1641
std::vector< std::vector< std::pair< value_type, value_type > > > edge_constrained_values
Definition: operators.h:469
static ::ExceptionBase & ExcNotInitialized()
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > data
Definition: operators.h:433
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1570
#define AssertThrow(cond, exc)
Definition: exceptions.h:1523
constexpr SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
typename VectorType::value_type value_type
Definition: operators.h:195
typename OperatorType::size_type size_type
Definition: operators.h:548
void resize(const size_type size_in)
AlignedVector< Number > shape_values
Definition: shape_info.h:155
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > get_matrix_free() const
Definition: operators.h:1598
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1727
void vmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1331
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal() const
Definition: operators.h:1620
static ::ExceptionBase & ExcMessage(std::string arg1)
AlignedVector< VectorizedArrayType > inverse_shape
Definition: operators.h:714
void vmult_interface_up(VectorType &dst, const VectorType &src) const
Definition: operators.h:1529
size_type n() const
#define Assert(cond, exc)
Definition: exceptions.h:1411
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:2087
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1895
void postprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1441
constexpr unsigned int pow(const unsigned int base, const int iexp)
Definition: utilities.h:446
virtual void Tapply_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1631
std::shared_ptr< Table< 2, VectorizedArrayType > > scalar_coefficient
Definition: operators.h:935
const unsigned int dofs_per_component
std::vector< unsigned int > selected_rows
Definition: operators.h:451
void reinit(const unsigned int cell_batch_index)
virtual void clear()
Definition: operators.h:1110
void set_constrained_entries_to_one(VectorType &dst) const
Definition: operators.h:1299
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1135
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1687
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition: operators.h:1811
void Tvmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1342
static constexpr unsigned int n_components
virtual void apply_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:2005
size_type m() const
Definition: operators.h:1082
static constexpr unsigned int static_dofs_per_cell
typename VectorType::size_type size_type
Definition: operators.h:200
void mult_add(VectorType &dst, const VectorType &src, const bool transpose) const
Definition: operators.h:1421
std::shared_ptr< Table< 2, VectorizedArrayType > > get_coefficient()
Definition: operators.h:1936
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:2041
const VectorizedArrayType * begin_dof_values() const
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info() const
std::vector< unsigned int > selected_columns
Definition: operators.h:457
void vmult_interface_down(VectorType &dst, const VectorType &src) const
Definition: operators.h:1475
void distribute_local_to_global(VectorType &dst, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::n_array_elements > &mask=std::bitset< VectorizedArrayType::n_array_elements >().flip()) const
size_type size() const
static ::ExceptionBase & ExcNotImplemented()
value_type el(const unsigned int row, const unsigned int col) const
Definition: operators.h:1120
Definition: table.h:699
bool is_element(const size_type index) const
Definition: index_set.h:1739
void initialize(const OperatorType &operator_in)
Definition: operators.h:1677
SmartPointer< const OperatorType > mf_base_operator
Definition: operators.h:593
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1707
void apply(const AlignedVector< VectorizedArrayType > &inverse_coefficient, const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition: operators.h:1026
std::shared_ptr< DiagonalMatrix< VectorType > > inverse_diagonal_entries
Definition: operators.h:445
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1318
virtual void compute_diagonal() override
Definition: operators.h:1768
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:1833
std::vector< std::vector< unsigned int > > edge_constrained_indices
Definition: operators.h:463
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:1160
void preprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1391
void transform_from_q_points_to_basis(const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition: operators.h:1055
void set_coefficient(const std::shared_ptr< Table< 2, VectorizedArrayType >> &scalar_coefficient)
Definition: operators.h:1916
static ::ExceptionBase & ExcInternalError()
void adjust_ghost_range_if_necessary(const VectorType &vec, const bool is_row) const
Definition: operators.h:1353