Reference documentation for deal.II version Git 6d836ad036 2021-06-16 20:57:45 +0200
\(\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 - 2021 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii_matrix_free_operators_h
18 #define dealii_matrix_free_operators_h
19 
20 
21 #include <deal.II/base/config.h>
22 
26 
29 
32 
34 
35 
37 
38 
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
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
98  {
99  vector.collect_sizes();
100  }
101 
102  template <typename VectorType>
103  typename std::enable_if<!IsBlockVector<VectorType>::value, void>::type
105  {}
106  } // namespace BlockHelper
107 
185  template <int dim,
187  typename VectorizedArrayType =
189  class Base : public Subscriptor
190  {
191  public:
195  using value_type = typename VectorType::value_type;
196 
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 
534  template <typename OperatorType>
536  {
537  public:
541  using value_type = typename OperatorType::value_type;
542 
547 
552 
556  void
557  clear();
558 
562  void
563  initialize(const OperatorType &operator_in);
564 
568  template <typename VectorType>
569  void
570  vmult(VectorType &dst, const VectorType &src) const;
571 
575  template <typename VectorType>
576  void
577  Tvmult(VectorType &dst, const VectorType &src) const;
578 
582  template <typename VectorType>
583  void
584  initialize_dof_vector(VectorType &vec) const;
585 
586 
587  private:
592  };
593 
594 
595 
613  template <int dim,
614  int fe_degree,
615  int n_components = 1,
616  typename Number = double,
617  typename VectorizedArrayType = VectorizedArray<Number>>
619  {
620  static_assert(
621  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
622  "Type of Number and of VectorizedArrayType do not match.");
623 
624  public:
630  const FEEvaluationBase<dim,
631  n_components,
632  Number,
633  false,
634  VectorizedArrayType> &fe_eval);
635 
644  void
645  apply(const AlignedVector<VectorizedArrayType> &inverse_coefficient,
646  const unsigned int n_actual_components,
647  const VectorizedArrayType * in_array,
648  VectorizedArrayType * out_array) const;
649 
661  void
662  apply(const VectorizedArrayType *in_array,
663  VectorizedArrayType * out_array) const;
664 
698  void
699  transform_from_q_points_to_basis(const unsigned int n_actual_components,
700  const VectorizedArrayType *in_array,
701  VectorizedArrayType *out_array) const;
702 
708  void
709  fill_inverse_JxW_values(
710  AlignedVector<VectorizedArrayType> &inverse_jxw) const;
711 
712  private:
716  const FEEvaluationBase<dim,
717  n_components,
718  Number,
719  false,
720  VectorizedArrayType> &fe_eval;
721  };
722 
723 
724 
732  template <int dim,
733  int fe_degree,
734  int n_q_points_1d = fe_degree + 1,
735  int n_components = 1,
736  typename VectorType = LinearAlgebra::distributed::Vector<double>,
737  typename VectorizedArrayType =
739  class MassOperator : public Base<dim, VectorType, VectorizedArrayType>
740  {
741  public:
745  using value_type =
747 
751  using size_type =
753 
757  MassOperator();
758 
763  virtual void
764  compute_diagonal() override;
765 
766  private:
772  virtual void
773  apply_add(VectorType &dst, const VectorType &src) const override;
774 
778  void
779  local_apply_cell(
781  VectorType & dst,
782  const VectorType & src,
783  const std::pair<unsigned int, unsigned int> & cell_range) const;
784  };
785 
786 
787 
798  template <int dim,
799  int fe_degree,
800  int n_q_points_1d = fe_degree + 1,
801  int n_components = 1,
802  typename VectorType = LinearAlgebra::distributed::Vector<double>,
803  typename VectorizedArrayType =
805  class LaplaceOperator : public Base<dim, VectorType, VectorizedArrayType>
806  {
807  public:
811  using value_type =
813 
817  using size_type =
819 
823  LaplaceOperator();
824 
831  virtual void
832  compute_diagonal() override;
833 
884  void
885  set_coefficient(
886  const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient);
887 
892  virtual void
893  clear() override;
894 
901  std::shared_ptr<Table<2, VectorizedArrayType>>
902  get_coefficient();
903 
904  private:
910  virtual void
911  apply_add(VectorType &dst, const VectorType &src) const override;
912 
916  void
917  local_apply_cell(
919  VectorType & dst,
920  const VectorType & src,
921  const std::pair<unsigned int, unsigned int> & cell_range) const;
922 
926  void
927  local_diagonal_cell(
929  VectorType & dst,
930  const VectorType &,
931  const std::pair<unsigned int, unsigned int> &cell_range) const;
932 
936  void
937  do_operation_on_cell(
939  & phi,
940  const unsigned int cell) const;
941 
945  std::shared_ptr<Table<2, VectorizedArrayType>> scalar_coefficient;
946  };
947 
948 
949 
950  // ------------------------------------ inline functions ---------------------
951 
952  template <int dim,
953  int fe_degree,
954  int n_components,
955  typename Number,
956  typename VectorizedArrayType>
957  inline CellwiseInverseMassMatrix<dim,
958  fe_degree,
959  n_components,
960  Number,
961  VectorizedArrayType>::
962  CellwiseInverseMassMatrix(
963  const FEEvaluationBase<dim,
964  n_components,
965  Number,
966  false,
967  VectorizedArrayType> &fe_eval)
968  : fe_eval(fe_eval)
969  {}
970 
971 
972 
973  template <int dim,
974  int fe_degree,
975  int n_components,
976  typename Number,
977  typename VectorizedArrayType>
978  inline void
980  fe_degree,
981  n_components,
982  Number,
983  VectorizedArrayType>::
985  AlignedVector<VectorizedArrayType> &inverse_jxw) const
986  {
987  constexpr unsigned int dofs_per_component_on_cell =
988  Utilities::pow(fe_degree + 1, dim);
989  Assert(inverse_jxw.size() > 0 &&
990  inverse_jxw.size() % dofs_per_component_on_cell == 0,
991  ExcMessage(
992  "Expected diagonal to be a multiple of scalar dof per cells"));
993 
994  // compute values for the first component
995  for (unsigned int q = 0; q < dofs_per_component_on_cell; ++q)
996  inverse_jxw[q] = 1. / fe_eval.JxW(q);
997  // copy values to rest of vector
998  for (unsigned int q = dofs_per_component_on_cell; q < inverse_jxw.size();)
999  for (unsigned int i = 0; i < dofs_per_component_on_cell; ++i, ++q)
1000  inverse_jxw[q] = inverse_jxw[i];
1001  }
1002 
1003 
1004 
1005  template <int dim,
1006  int fe_degree,
1007  int n_components,
1008  typename Number,
1009  typename VectorizedArrayType>
1010  inline void
1012  dim,
1013  fe_degree,
1014  n_components,
1015  Number,
1016  VectorizedArrayType>::apply(const VectorizedArrayType *in_array,
1017  VectorizedArrayType * out_array) const
1018  {
1020  template run<fe_degree>(n_components, fe_eval, in_array, out_array);
1021  }
1022 
1023 
1024 
1025  template <int dim,
1026  int fe_degree,
1027  int n_components,
1028  typename Number,
1029  typename VectorizedArrayType>
1030  inline void
1032  fe_degree,
1033  n_components,
1034  Number,
1035  VectorizedArrayType>::
1036  apply(const AlignedVector<VectorizedArrayType> &inverse_coefficients,
1037  const unsigned int n_actual_components,
1038  const VectorizedArrayType * in_array,
1039  VectorizedArrayType * out_array) const
1040  {
1042  template run<fe_degree>(
1043  n_actual_components,
1044  fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
1045  inverse_coefficients,
1046  in_array,
1047  out_array);
1048  }
1049 
1050 
1051 
1052  template <int dim,
1053  int fe_degree,
1054  int n_components,
1055  typename Number,
1056  typename VectorizedArrayType>
1057  inline void
1059  fe_degree,
1060  n_components,
1061  Number,
1062  VectorizedArrayType>::
1063  transform_from_q_points_to_basis(const unsigned int n_actual_components,
1064  const VectorizedArrayType *in_array,
1065  VectorizedArrayType * out_array) const
1066  {
1067  const auto n_q_points_1d = fe_eval.get_shape_info().data[0].n_q_points_1d;
1068 
1069  if (fe_degree > -1 && (fe_degree + 1 == n_q_points_1d))
1071  dim,
1072  VectorizedArrayType>::template run<fe_degree,
1073  fe_degree + 1>(n_actual_components,
1074  fe_eval,
1075  in_array,
1076  out_array);
1077  else
1079  transform_from_q_points_to_basis(n_actual_components,
1080  fe_degree,
1081  n_q_points_1d,
1082  fe_eval,
1083  in_array,
1084  out_array);
1085  }
1086 
1087 
1088 
1089  //----------------- Base operator -----------------------------
1090  template <int dim, typename VectorType, typename VectorizedArrayType>
1092  : Subscriptor()
1093  , have_interface_matrices(false)
1094  {}
1095 
1096 
1097 
1098  template <int dim, typename VectorType, typename VectorizedArrayType>
1101  {
1102  Assert(data.get() != nullptr, ExcNotInitialized());
1104  0;
1105  for (const unsigned int selected_row : selected_rows)
1106  total_size += data->get_vector_partitioner(selected_row)->size();
1107  return total_size;
1108  }
1109 
1110 
1111 
1112  template <int dim, typename VectorType, typename VectorizedArrayType>
1115  {
1116  Assert(data.get() != nullptr, ExcNotInitialized());
1118  0;
1119  for (const unsigned int selected_column : selected_columns)
1120  total_size += data->get_vector_partitioner(selected_column)->size();
1121  return total_size;
1122  }
1123 
1124 
1125 
1126  template <int dim, typename VectorType, typename VectorizedArrayType>
1127  void
1129  {
1130  data.reset();
1131  inverse_diagonal_entries.reset();
1132  }
1133 
1134 
1135 
1136  template <int dim, typename VectorType, typename VectorizedArrayType>
1139  const unsigned int col) const
1140  {
1141  (void)col;
1142  Assert(row == col, ExcNotImplemented());
1143  Assert(inverse_diagonal_entries.get() != nullptr &&
1144  inverse_diagonal_entries->m() > 0,
1145  ExcNotInitialized());
1146  return 1.0 / (*inverse_diagonal_entries)(row, row);
1147  }
1148 
1149 
1150 
1151  template <int dim, typename VectorType, typename VectorizedArrayType>
1152  void
1154  VectorType &vec) const
1155  {
1156  Assert(data.get() != nullptr, ExcNotInitialized());
1158  for (unsigned int i = 0; i < BlockHelper::n_blocks(vec); ++i)
1159  {
1160  const unsigned int index = selected_rows[i];
1161  if (!BlockHelper::subblock(vec, index)
1162  .partitioners_are_compatible(
1163  *data->get_dof_info(index).vector_partitioner))
1164  data->initialize_dof_vector(BlockHelper::subblock(vec, index), index);
1165 
1166  Assert(BlockHelper::subblock(vec, index)
1167  .partitioners_are_globally_compatible(
1168  *data->get_dof_info(index).vector_partitioner),
1169  ExcInternalError());
1170  }
1172  }
1173 
1174 
1175 
1176  template <int dim, typename VectorType, typename VectorizedArrayType>
1177  void
1180  data_,
1181  const std::vector<unsigned int> &given_row_selection,
1182  const std::vector<unsigned int> &given_column_selection)
1183  {
1184  data = data_;
1185 
1186  selected_rows.clear();
1187  selected_columns.clear();
1188  if (given_row_selection.empty())
1189  for (unsigned int i = 0; i < data_->n_components(); ++i)
1190  selected_rows.push_back(i);
1191  else
1192  {
1193  for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1194  {
1195  AssertIndexRange(given_row_selection[i], data_->n_components());
1196  for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1197  if (j != i)
1198  Assert(given_row_selection[j] != given_row_selection[i],
1199  ExcMessage("Given row indices must be unique"));
1200 
1201  selected_rows.push_back(given_row_selection[i]);
1202  }
1203  }
1204  if (given_column_selection.size() == 0)
1206  else
1207  {
1208  for (unsigned int i = 0; i < given_column_selection.size(); ++i)
1209  {
1210  AssertIndexRange(given_column_selection[i], data_->n_components());
1211  for (unsigned int j = 0; j < given_column_selection.size(); ++j)
1212  if (j != i)
1213  Assert(given_column_selection[j] != given_column_selection[i],
1214  ExcMessage("Given column indices must be unique"));
1215 
1216  selected_columns.push_back(given_column_selection[i]);
1217  }
1218  }
1219 
1220  edge_constrained_indices.clear();
1221  edge_constrained_indices.resize(selected_rows.size());
1222  edge_constrained_values.clear();
1223  edge_constrained_values.resize(selected_rows.size());
1224  have_interface_matrices = false;
1225  }
1226 
1227 
1228 
1229  template <int dim, typename VectorType, typename VectorizedArrayType>
1230  void
1233  data_,
1234  const MGConstrainedDoFs & mg_constrained_dofs,
1235  const unsigned int level,
1236  const std::vector<unsigned int> &given_row_selection)
1237  {
1238  std::vector<MGConstrainedDoFs> mg_constrained_dofs_vector(
1239  1, mg_constrained_dofs);
1240  initialize(data_, mg_constrained_dofs_vector, level, given_row_selection);
1241  }
1242 
1243 
1244 
1245  template <int dim, typename VectorType, typename VectorizedArrayType>
1246  void
1249  data_,
1250  const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
1251  const unsigned int level,
1252  const std::vector<unsigned int> & given_row_selection)
1253  {
1255  ExcMessage("level is not set"));
1256 
1257  selected_rows.clear();
1258  selected_columns.clear();
1259  if (given_row_selection.empty())
1260  for (unsigned int i = 0; i < data_->n_components(); ++i)
1261  selected_rows.push_back(i);
1262  else
1263  {
1264  for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1265  {
1266  AssertIndexRange(given_row_selection[i], data_->n_components());
1267  for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1268  if (j != i)
1269  Assert(given_row_selection[j] != given_row_selection[i],
1270  ExcMessage("Given row indices must be unique"));
1271 
1272  selected_rows.push_back(given_row_selection[i]);
1273  }
1274  }
1276 
1277  AssertDimension(mg_constrained_dofs.size(), selected_rows.size());
1278  edge_constrained_indices.clear();
1279  edge_constrained_indices.resize(selected_rows.size());
1280  edge_constrained_values.clear();
1281  edge_constrained_values.resize(selected_rows.size());
1282 
1283  data = data_;
1284 
1285  for (unsigned int j = 0; j < selected_rows.size(); ++j)
1286  {
1287  if (data_->n_cell_batches() > 0)
1288  {
1289  AssertDimension(level, data_->get_cell_iterator(0, 0, j)->level());
1290  }
1291 
1292  // setup edge_constrained indices
1293  std::vector<types::global_dof_index> interface_indices;
1294  mg_constrained_dofs[j]
1295  .get_refinement_edge_indices(level)
1296  .fill_index_vector(interface_indices);
1297  edge_constrained_indices[j].clear();
1298  edge_constrained_indices[j].reserve(interface_indices.size());
1299  edge_constrained_values[j].resize(interface_indices.size());
1300  const IndexSet &locally_owned =
1301  data->get_dof_handler(selected_rows[j]).locally_owned_mg_dofs(level);
1302  for (const auto interface_index : interface_indices)
1303  if (locally_owned.is_element(interface_index))
1304  edge_constrained_indices[j].push_back(
1305  locally_owned.index_within_set(interface_index));
1308  static_cast<unsigned int>(edge_constrained_indices[j].size()),
1309  data->get_vector_partitioner()->get_mpi_communicator()) > 0;
1310  }
1311  }
1312 
1313 
1314 
1315  template <int dim, typename VectorType, typename VectorizedArrayType>
1316  void
1318  VectorType &dst) const
1319  {
1320  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1321  {
1322  const std::vector<unsigned int> &constrained_dofs =
1323  data->get_constrained_dofs(selected_rows[j]);
1324  for (const auto constrained_dof : constrained_dofs)
1325  BlockHelper::subblock(dst, j).local_element(constrained_dof) = 1.;
1326  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1327  BlockHelper::subblock(dst, j).local_element(
1328  edge_constrained_indices[j][i]) = 1.;
1329  }
1330  }
1331 
1332 
1333 
1334  template <int dim, typename VectorType, typename VectorizedArrayType>
1335  void
1337  const VectorType &src) const
1338  {
1339  using Number =
1341  dst = Number(0.);
1342  vmult_add(dst, src);
1343  }
1344 
1345 
1346 
1347  template <int dim, typename VectorType, typename VectorizedArrayType>
1348  void
1350  VectorType & dst,
1351  const VectorType &src) const
1352  {
1353  mult_add(dst, src, false);
1354  }
1355 
1356 
1357 
1358  template <int dim, typename VectorType, typename VectorizedArrayType>
1359  void
1361  VectorType & dst,
1362  const VectorType &src) const
1363  {
1364  mult_add(dst, src, true);
1365  }
1366 
1367 
1368 
1369  template <int dim, typename VectorType, typename VectorizedArrayType>
1370  void
1372  const VectorType &src,
1373  const bool is_row) const
1374  {
1375  using Number =
1377  for (unsigned int i = 0; i < BlockHelper::n_blocks(src); ++i)
1378  {
1379  const unsigned int mf_component =
1380  is_row ? selected_rows[i] : selected_columns[i];
1381  // If both vectors use the same partitioner -> done
1382  if (BlockHelper::subblock(src, i).get_partitioner().get() ==
1383  data->get_dof_info(mf_component).vector_partitioner.get())
1384  continue;
1385 
1386  // If not, assert that the local ranges are the same and reset to the
1387  // current partitioner
1389  .get_partitioner()
1390  ->locally_owned_size() ==
1391  data->get_dof_info(mf_component)
1392  .vector_partitioner->locally_owned_size(),
1393  ExcMessage(
1394  "The vector passed to the vmult() function does not have "
1395  "the correct size for compatibility with MatrixFree."));
1396 
1397  // copy the vector content to a temporary vector so that it does not get
1398  // lost
1400  BlockHelper::subblock(src, i));
1401  this->data->initialize_dof_vector(
1402  BlockHelper::subblock(const_cast<VectorType &>(src), i),
1403  mf_component);
1404  BlockHelper::subblock(const_cast<VectorType &>(src), i)
1405  .copy_locally_owned_data_from(copy_vec);
1406  }
1407  }
1408 
1409 
1410 
1411  template <int dim, typename VectorType, typename VectorizedArrayType>
1412  void
1414  VectorType & dst,
1415  const VectorType &src) const
1416  {
1417  using Number =
1419  adjust_ghost_range_if_necessary(src, false);
1421 
1422  // set zero Dirichlet values on the input vector (and remember the src and
1423  // dst values because we need to reset them at the end)
1424  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1425  {
1426  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1427  {
1428  edge_constrained_values[j][i] = std::pair<Number, Number>(
1429  BlockHelper::subblock(src, j).local_element(
1430  edge_constrained_indices[j][i]),
1431  BlockHelper::subblock(dst, j).local_element(
1432  edge_constrained_indices[j][i]));
1433  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1434  .local_element(edge_constrained_indices[j][i]) = 0.;
1435  }
1436  }
1437  }
1438 
1439 
1440 
1441  template <int dim, typename VectorType, typename VectorizedArrayType>
1442  void
1444  VectorType & dst,
1445  const VectorType &src,
1446  const bool transpose) const
1447  {
1448  AssertDimension(dst.size(), src.size());
1451  preprocess_constraints(dst, src);
1452  if (transpose)
1453  Tapply_add(dst, src);
1454  else
1455  apply_add(dst, src);
1456  postprocess_constraints(dst, src);
1457  }
1458 
1459 
1460 
1461  template <int dim, typename VectorType, typename VectorizedArrayType>
1462  void
1464  VectorType & dst,
1465  const VectorType &src) const
1466  {
1467  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1468  {
1469  const std::vector<unsigned int> &constrained_dofs =
1470  data->get_constrained_dofs(selected_rows[j]);
1471  for (const auto constrained_dof : constrained_dofs)
1472  BlockHelper::subblock(dst, j).local_element(constrained_dof) +=
1473  BlockHelper::subblock(src, j).local_element(constrained_dof);
1474  }
1475 
1476  // reset edge constrained values, multiply by unit matrix and add into
1477  // destination
1478  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1479  {
1480  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1481  {
1482  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1483  .local_element(edge_constrained_indices[j][i]) =
1484  edge_constrained_values[j][i].first;
1485  BlockHelper::subblock(dst, j).local_element(
1486  edge_constrained_indices[j][i]) =
1487  edge_constrained_values[j][i].second +
1488  edge_constrained_values[j][i].first;
1489  }
1490  }
1491  }
1492 
1493 
1494 
1495  template <int dim, typename VectorType, typename VectorizedArrayType>
1496  void
1498  VectorType & dst,
1499  const VectorType &src) const
1500  {
1501  using Number =
1503  AssertDimension(dst.size(), src.size());
1504  adjust_ghost_range_if_necessary(src, false);
1506 
1507  dst = Number(0.);
1508 
1510  return;
1511 
1512  // set zero Dirichlet values on the input vector (and remember the src and
1513  // dst values because we need to reset them at the end)
1514  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1515  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1516  {
1517  edge_constrained_values[j][i] = std::pair<Number, Number>(
1518  BlockHelper::subblock(src, j).local_element(
1519  edge_constrained_indices[j][i]),
1520  BlockHelper::subblock(dst, j).local_element(
1521  edge_constrained_indices[j][i]));
1522  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1523  .local_element(edge_constrained_indices[j][i]) = 0.;
1524  }
1525 
1526  apply_add(dst, src);
1527 
1528  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1529  {
1530  unsigned int c = 0;
1531  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1532  {
1533  for (; c < edge_constrained_indices[j][i]; ++c)
1534  BlockHelper::subblock(dst, j).local_element(c) = 0.;
1535  ++c;
1536 
1537  // reset the src values
1538  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1539  .local_element(edge_constrained_indices[j][i]) =
1540  edge_constrained_values[j][i].first;
1541  }
1542  for (; c < BlockHelper::subblock(dst, j).locally_owned_size(); ++c)
1543  BlockHelper::subblock(dst, j).local_element(c) = 0.;
1544  }
1545  }
1546 
1547 
1548 
1549  template <int dim, typename VectorType, typename VectorizedArrayType>
1550  void
1552  VectorType & dst,
1553  const VectorType &src) const
1554  {
1555  using Number =
1557  AssertDimension(dst.size(), src.size());
1558  adjust_ghost_range_if_necessary(src, false);
1560 
1561  dst = Number(0.);
1562 
1564  return;
1565 
1566  VectorType src_cpy(src);
1567  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1568  {
1569  unsigned int c = 0;
1570  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1571  {
1572  for (; c < edge_constrained_indices[j][i]; ++c)
1573  BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1574  ++c;
1575  }
1576  for (; c < BlockHelper::subblock(src_cpy, j).locally_owned_size(); ++c)
1577  BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1578  }
1579 
1580  apply_add(dst, src_cpy);
1581 
1582  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1583  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1584  BlockHelper::subblock(dst, j).local_element(
1585  edge_constrained_indices[j][i]) = 0.;
1586  }
1587 
1588 
1589 
1590  template <int dim, typename VectorType, typename VectorizedArrayType>
1591  void
1593  VectorType & dst,
1594  const VectorType &src) const
1595  {
1596  using Number =
1598  dst = Number(0.);
1599  Tvmult_add(dst, src);
1600  }
1601 
1602 
1603 
1604  template <int dim, typename VectorType, typename VectorizedArrayType>
1605  std::size_t
1607  {
1608  return inverse_diagonal_entries.get() != nullptr ?
1609  inverse_diagonal_entries->memory_consumption() :
1610  sizeof(*this);
1611  }
1612 
1613 
1614 
1615  template <int dim, typename VectorType, typename VectorizedArrayType>
1616  std::shared_ptr<const MatrixFree<
1617  dim,
1619  VectorizedArrayType>>
1621  {
1622  return data;
1623  }
1624 
1625 
1626 
1627  template <int dim, typename VectorType, typename VectorizedArrayType>
1628  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1630  const
1631  {
1632  Assert(inverse_diagonal_entries.get() != nullptr &&
1633  inverse_diagonal_entries->m() > 0,
1634  ExcNotInitialized());
1635  return inverse_diagonal_entries;
1636  }
1637 
1638 
1639 
1640  template <int dim, typename VectorType, typename VectorizedArrayType>
1641  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1643  {
1644  Assert(diagonal_entries.get() != nullptr && diagonal_entries->m() > 0,
1645  ExcNotInitialized());
1646  return diagonal_entries;
1647  }
1648 
1649 
1650 
1651  template <int dim, typename VectorType, typename VectorizedArrayType>
1652  void
1654  VectorType & dst,
1655  const VectorType &src) const
1656  {
1657  apply_add(dst, src);
1658  }
1659 
1660 
1661 
1662  template <int dim, typename VectorType, typename VectorizedArrayType>
1663  void
1665  VectorType & dst,
1666  const VectorType & src,
1667  const typename Base<dim, VectorType, VectorizedArrayType>::value_type omega)
1668  const
1669  {
1671  ExcNotInitialized());
1672  inverse_diagonal_entries->vmult(dst, src);
1673  dst *= omega;
1674  }
1675 
1676 
1677 
1678  //------------------------- MGInterfaceOperator ------------------------------
1679 
1680  template <typename OperatorType>
1682  : Subscriptor()
1683  , mf_base_operator(nullptr)
1684  {}
1685 
1686 
1687 
1688  template <typename OperatorType>
1689  void
1691  {
1692  mf_base_operator = nullptr;
1693  }
1694 
1695 
1696 
1697  template <typename OperatorType>
1698  void
1699  MGInterfaceOperator<OperatorType>::initialize(const OperatorType &operator_in)
1700  {
1701  mf_base_operator = &operator_in;
1702  }
1703 
1704 
1705 
1706  template <typename OperatorType>
1707  template <typename VectorType>
1708  void
1710  const VectorType &src) const
1711  {
1712 #ifndef DEAL_II_MSVC
1713  static_assert(
1714  std::is_same<typename VectorType::value_type, value_type>::value,
1715  "The vector type must be based on the same value type as this "
1716  "operator");
1717 #endif
1718 
1719  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1720 
1721  mf_base_operator->vmult_interface_down(dst, src);
1722  }
1723 
1724 
1725 
1726  template <typename OperatorType>
1727  template <typename VectorType>
1728  void
1730  const VectorType &src) const
1731  {
1732 #ifndef DEAL_II_MSVC
1733  static_assert(
1734  std::is_same<typename VectorType::value_type, value_type>::value,
1735  "The vector type must be based on the same value type as this "
1736  "operator");
1737 #endif
1738 
1739  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1740 
1741  mf_base_operator->vmult_interface_up(dst, src);
1742  }
1743 
1744 
1745 
1746  template <typename OperatorType>
1747  template <typename VectorType>
1748  void
1750  VectorType &vec) const
1751  {
1752  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1753 
1754  mf_base_operator->initialize_dof_vector(vec);
1755  }
1756 
1757 
1758 
1759  //-----------------------------MassOperator----------------------------------
1760 
1761  template <int dim,
1762  int fe_degree,
1763  int n_q_points_1d,
1764  int n_components,
1765  typename VectorType,
1766  typename VectorizedArrayType>
1767  MassOperator<dim,
1768  fe_degree,
1769  n_q_points_1d,
1770  n_components,
1771  VectorType,
1772  VectorizedArrayType>::MassOperator()
1773  : Base<dim, VectorType, VectorizedArrayType>()
1774  {}
1775 
1776 
1777 
1778  template <int dim,
1779  int fe_degree,
1780  int n_q_points_1d,
1781  int n_components,
1782  typename VectorType,
1783  typename VectorizedArrayType>
1784  void
1785  MassOperator<dim,
1786  fe_degree,
1787  n_q_points_1d,
1788  n_components,
1789  VectorType,
1790  VectorizedArrayType>::compute_diagonal()
1791  {
1792  using Number =
1793  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1795  ExcNotInitialized());
1796 
1797  this->inverse_diagonal_entries =
1798  std::make_shared<DiagonalMatrix<VectorType>>();
1799  this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1800  VectorType &inverse_diagonal_vector =
1801  this->inverse_diagonal_entries->get_vector();
1802  VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1803  this->initialize_dof_vector(inverse_diagonal_vector);
1804  this->initialize_dof_vector(diagonal_vector);
1805  inverse_diagonal_vector = Number(1.);
1806  apply_add(diagonal_vector, inverse_diagonal_vector);
1807 
1808  this->set_constrained_entries_to_one(diagonal_vector);
1809  inverse_diagonal_vector = diagonal_vector;
1810 
1811  const unsigned int locally_owned_size =
1812  inverse_diagonal_vector.locally_owned_size();
1813  for (unsigned int i = 0; i < locally_owned_size; ++i)
1814  inverse_diagonal_vector.local_element(i) =
1815  Number(1.) / inverse_diagonal_vector.local_element(i);
1816 
1817  inverse_diagonal_vector.update_ghost_values();
1818  diagonal_vector.update_ghost_values();
1819  }
1820 
1821 
1822 
1823  template <int dim,
1824  int fe_degree,
1825  int n_q_points_1d,
1826  int n_components,
1827  typename VectorType,
1828  typename VectorizedArrayType>
1829  void
1830  MassOperator<dim,
1831  fe_degree,
1832  n_q_points_1d,
1833  n_components,
1834  VectorType,
1835  VectorizedArrayType>::apply_add(VectorType & dst,
1836  const VectorType &src) const
1837  {
1839  &MassOperator::local_apply_cell, this, dst, src);
1840  }
1841 
1842 
1843 
1844  template <int dim,
1845  int fe_degree,
1846  int n_q_points_1d,
1847  int n_components,
1848  typename VectorType,
1849  typename VectorizedArrayType>
1850  void
1851  MassOperator<dim,
1852  fe_degree,
1853  n_q_points_1d,
1854  n_components,
1855  VectorType,
1856  VectorizedArrayType>::
1858  const MatrixFree<
1859  dim,
1860  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
1861  VectorizedArrayType> & data,
1862  VectorType & dst,
1863  const VectorType & src,
1864  const std::pair<unsigned int, unsigned int> &cell_range) const
1865  {
1866  using Number =
1867  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1868  FEEvaluation<dim,
1869  fe_degree,
1870  n_q_points_1d,
1871  n_components,
1872  Number,
1873  VectorizedArrayType>
1874  phi(data, this->selected_rows[0]);
1875  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
1876  {
1877  phi.reinit(cell);
1878  phi.read_dof_values(src);
1879  phi.evaluate(EvaluationFlags::values);
1880  for (unsigned int q = 0; q < phi.n_q_points; ++q)
1881  phi.submit_value(phi.get_value(q), q);
1882  phi.integrate(EvaluationFlags::values);
1883  phi.distribute_local_to_global(dst);
1884  }
1885  }
1886 
1887 
1888  //-----------------------------LaplaceOperator----------------------------------
1889 
1890  template <int dim,
1891  int fe_degree,
1892  int n_q_points_1d,
1893  int n_components,
1894  typename VectorType,
1895  typename VectorizedArrayType>
1896  LaplaceOperator<dim,
1897  fe_degree,
1898  n_q_points_1d,
1899  n_components,
1900  VectorType,
1901  VectorizedArrayType>::LaplaceOperator()
1902  : Base<dim, VectorType, VectorizedArrayType>()
1903  {}
1904 
1905 
1906 
1907  template <int dim,
1908  int fe_degree,
1909  int n_q_points_1d,
1910  int n_components,
1911  typename VectorType,
1912  typename VectorizedArrayType>
1913  void
1914  LaplaceOperator<dim,
1915  fe_degree,
1916  n_q_points_1d,
1917  n_components,
1918  VectorType,
1919  VectorizedArrayType>::clear()
1920  {
1922  scalar_coefficient.reset();
1923  }
1924 
1925 
1926 
1927  template <int dim,
1928  int fe_degree,
1929  int n_q_points_1d,
1930  int n_components,
1931  typename VectorType,
1932  typename VectorizedArrayType>
1933  void
1934  LaplaceOperator<dim,
1935  fe_degree,
1936  n_q_points_1d,
1937  n_components,
1938  VectorType,
1939  VectorizedArrayType>::
1941  const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient_)
1942  {
1943  scalar_coefficient = scalar_coefficient_;
1944  }
1945 
1946 
1947 
1948  template <int dim,
1949  int fe_degree,
1950  int n_q_points_1d,
1951  int n_components,
1952  typename VectorType,
1953  typename VectorizedArrayType>
1954  std::shared_ptr<Table<2, VectorizedArrayType>>
1955  LaplaceOperator<dim,
1956  fe_degree,
1957  n_q_points_1d,
1958  n_components,
1959  VectorType,
1960  VectorizedArrayType>::get_coefficient()
1961  {
1963  return scalar_coefficient;
1964  }
1965 
1966 
1967 
1968  template <int dim,
1969  int fe_degree,
1970  int n_q_points_1d,
1971  int n_components,
1972  typename VectorType,
1973  typename VectorizedArrayType>
1974  void
1975  LaplaceOperator<dim,
1976  fe_degree,
1977  n_q_points_1d,
1978  n_components,
1979  VectorType,
1980  VectorizedArrayType>::compute_diagonal()
1981  {
1982  using Number =
1983  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1985  ExcNotInitialized());
1986 
1987  this->inverse_diagonal_entries =
1988  std::make_shared<DiagonalMatrix<VectorType>>();
1989  this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1990  VectorType &inverse_diagonal_vector =
1991  this->inverse_diagonal_entries->get_vector();
1992  VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1993  this->initialize_dof_vector(inverse_diagonal_vector);
1994  this->initialize_dof_vector(diagonal_vector);
1995 
1996  this->data->cell_loop(&LaplaceOperator::local_diagonal_cell,
1997  this,
1998  diagonal_vector,
1999  /*unused*/ diagonal_vector);
2000  this->set_constrained_entries_to_one(diagonal_vector);
2001 
2002  inverse_diagonal_vector = diagonal_vector;
2003 
2004  for (unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
2005  ++i)
2006  if (std::abs(inverse_diagonal_vector.local_element(i)) >
2008  inverse_diagonal_vector.local_element(i) =
2009  1. / inverse_diagonal_vector.local_element(i);
2010  else
2011  inverse_diagonal_vector.local_element(i) = 1.;
2012 
2013  inverse_diagonal_vector.update_ghost_values();
2014  diagonal_vector.update_ghost_values();
2015  }
2016 
2017 
2018 
2019  template <int dim,
2020  int fe_degree,
2021  int n_q_points_1d,
2022  int n_components,
2023  typename VectorType,
2024  typename VectorizedArrayType>
2025  void
2026  LaplaceOperator<dim,
2027  fe_degree,
2028  n_q_points_1d,
2029  n_components,
2030  VectorType,
2031  VectorizedArrayType>::apply_add(VectorType & dst,
2032  const VectorType &src) const
2033  {
2035  &LaplaceOperator::local_apply_cell, this, dst, src);
2036  }
2037 
2038  namespace Implementation
2039  {
2040  template <typename VectorizedArrayType>
2041  bool
2042  non_negative(const VectorizedArrayType &n)
2043  {
2044  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2045  if (n[v] < 0.)
2046  return false;
2047 
2048  return true;
2049  }
2050  } // namespace Implementation
2051 
2052 
2053 
2054  template <int dim,
2055  int fe_degree,
2056  int n_q_points_1d,
2057  int n_components,
2058  typename VectorType,
2059  typename VectorizedArrayType>
2060  void
2061  LaplaceOperator<dim,
2062  fe_degree,
2063  n_q_points_1d,
2064  n_components,
2065  VectorType,
2066  VectorizedArrayType>::
2068  FEEvaluation<
2069  dim,
2070  fe_degree,
2071  n_q_points_1d,
2072  n_components,
2073  typename Base<dim, VectorType, VectorizedArrayType>::value_type> &phi,
2074  const unsigned int cell) const
2075  {
2076  phi.evaluate(EvaluationFlags::gradients);
2077  if (scalar_coefficient.get())
2078  {
2079  Assert(scalar_coefficient->size(1) == 1 ||
2080  scalar_coefficient->size(1) == phi.n_q_points,
2081  ExcMessage("The number of columns in the coefficient table must "
2082  "be either 1 or the number of quadrature points " +
2083  std::to_string(phi.n_q_points) +
2084  ", but the given value was " +
2085  std::to_string(scalar_coefficient->size(1))));
2086  if (scalar_coefficient->size(1) == phi.n_q_points)
2087  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2088  {
2090  (*scalar_coefficient)(cell, q)),
2091  ExcMessage("Coefficient must be non-negative"));
2092  phi.submit_gradient((*scalar_coefficient)(cell, q) *
2093  phi.get_gradient(q),
2094  q);
2095  }
2096  else
2097  {
2099  ExcMessage("Coefficient must be non-negative"));
2100  const VectorizedArrayType coefficient =
2101  (*scalar_coefficient)(cell, 0);
2102  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2103  phi.submit_gradient(coefficient * phi.get_gradient(q), q);
2104  }
2105  }
2106  else
2107  {
2108  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2109  {
2110  phi.submit_gradient(phi.get_gradient(q), q);
2111  }
2112  }
2113  phi.integrate(EvaluationFlags::gradients);
2114  }
2115 
2116 
2117 
2118  template <int dim,
2119  int fe_degree,
2120  int n_q_points_1d,
2121  int n_components,
2122  typename VectorType,
2123  typename VectorizedArrayType>
2124  void
2125  LaplaceOperator<dim,
2126  fe_degree,
2127  n_q_points_1d,
2128  n_components,
2129  VectorType,
2130  VectorizedArrayType>::
2132  const MatrixFree<
2133  dim,
2134  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2135  VectorizedArrayType> & data,
2136  VectorType & dst,
2137  const VectorType & src,
2138  const std::pair<unsigned int, unsigned int> &cell_range) const
2139  {
2140  using Number =
2141  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2143  data, this->selected_rows[0]);
2144  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2145  {
2146  phi.reinit(cell);
2147  phi.read_dof_values(src);
2148  do_operation_on_cell(phi, cell);
2149  phi.distribute_local_to_global(dst);
2150  }
2151  }
2152 
2153 
2154  template <int dim,
2155  int fe_degree,
2156  int n_q_points_1d,
2157  int n_components,
2158  typename VectorType,
2159  typename VectorizedArrayType>
2160  void
2161  LaplaceOperator<dim,
2162  fe_degree,
2163  n_q_points_1d,
2164  n_components,
2165  VectorType,
2166  VectorizedArrayType>::
2168  const MatrixFree<
2169  dim,
2170  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2171  VectorizedArrayType> &data,
2172  VectorType & dst,
2173  const VectorType &,
2174  const std::pair<unsigned int, unsigned int> &cell_range) const
2175  {
2176  using Number =
2177  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2178 
2180  data, this->selected_rows[0]);
2181  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2182  {
2183  phi.reinit(cell);
2184  VectorizedArrayType local_diagonal_vector[phi.static_dofs_per_cell];
2185  for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2186  {
2187  for (unsigned int j = 0; j < phi.dofs_per_component; ++j)
2188  phi.begin_dof_values()[j] = VectorizedArrayType();
2189  phi.begin_dof_values()[i] = 1.;
2190  do_operation_on_cell(phi, cell);
2191  local_diagonal_vector[i] = phi.begin_dof_values()[i];
2192  }
2193  for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2194  for (unsigned int c = 0; c < phi.n_components; ++c)
2195  phi.begin_dof_values()[i + c * phi.dofs_per_component] =
2196  local_diagonal_vector[i];
2197  phi.distribute_local_to_global(dst);
2198  }
2199  }
2200 
2201 
2202 } // end of namespace MatrixFreeOperators
2203 
2204 
2206 
2207 #endif
typename OperatorType::value_type value_type
Definition: operators.h:541
virtual void apply_add(VectorType &dst, const VectorType &src) const =0
void fill_inverse_JxW_values(AlignedVector< VectorizedArrayType > &inverse_jxw) const
Definition: operators.h:984
static const unsigned int invalid_unsigned_int
Definition: types.h:196
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal_inverse() const
Definition: operators.h:1629
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1656
virtual std::size_t memory_consumption() const
Definition: operators.h:1606
const FEEvaluationBase< dim, n_components, Number, false, VectorizedArrayType > & fe_eval
Definition: operators.h:720
types::global_dof_index size_type
Definition: cuda_kernels.h:45
size_type n() const
Definition: operators.h:1114
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:2167
void precondition_Jacobi(VectorType &dst, const VectorType &src, const value_type omega) const
Definition: operators.h:1664
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:1721
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:1592
#define AssertThrow(cond, exc)
Definition: exceptions.h:1588
std::enable_if< IsBlockVector< VectorType >::value, void >::type collect_sizes(VectorType &vector)
Definition: operators.h:97
typename VectorType::value_type value_type
Definition: operators.h:195
typename OperatorType::size_type size_type
Definition: operators.h:546
virtual void compute_diagonal() override
Definition: operators.h:1980
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:461
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > get_matrix_free() const
Definition: operators.h:1620
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1749
void vmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1349
static void transform_from_q_points_to_basis(const unsigned int n_components, const unsigned int fe_degree, const unsigned int n_q_points_1d, const FEEvaluationBaseData< dim, Number, false, VectorizedArrayType > &fe_eval, const VectorizedArrayType *in_array, VectorizedArrayType *out_array)
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal() const
Definition: operators.h:1642
static ::ExceptionBase & ExcMessage(std::string arg1)
VectorizedArrayType JxW(const unsigned int q_point) const
void vmult_interface_up(VectorType &dst, const VectorType &src) const
Definition: operators.h:1551
std::enable_if< IsBlockVector< VectorType >::value, typename VectorType::BlockType & >::type subblock(VectorType &vector, unsigned int block_no)
Definition: operators.h:65
#define Assert(cond, exc)
Definition: exceptions.h:1478
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:2131
auto apply(F &&fn, Tuple &&t) -> decltype(apply_impl(std::forward< F >(fn), std::forward< Tuple >(t), std::make_index_sequence< std::tuple_size< typename std::remove_reference< Tuple >::type >::value >()))
Definition: tuple.h:36
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1921
void postprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1463
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
virtual void Tapply_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1653
std::shared_ptr< Table< 2, VectorizedArrayType > > scalar_coefficient
Definition: operators.h:945
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:399
unsigned int level
Definition: grid_out.cc:4590
const unsigned int dofs_per_component
std::vector< unsigned int > selected_rows
Definition: operators.h:451
std::string to_string(const T &t)
Definition: patterns.h:2329
void reinit(const unsigned int cell_batch_index)
std::enable_if< IsBlockVector< VectorType >::value, unsigned int >::type n_blocks(const VectorType &vector)
Definition: operators.h:49
virtual void clear()
Definition: operators.h:1128
void set_constrained_entries_to_one(VectorType &dst) const
Definition: operators.h:1317
virtual void clear() override
Definition: operators.h:1919
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1153
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1709
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition: operators.h:2031
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition: operators.h:1835
void Tvmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1360
static constexpr unsigned int n_components
size_type m() const
Definition: operators.h:1100
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:1443
std::shared_ptr< Table< 2, VectorizedArrayType > > get_coefficient()
Definition: operators.h:1960
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:2067
const VectorizedArrayType * begin_dof_values() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:398
std::vector< unsigned int > selected_columns
Definition: operators.h:457
void vmult_interface_down(VectorType &dst, const VectorType &src) const
Definition: operators.h:1497
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info() const
void distribute_local_to_global(VectorType &dst, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip()) const
size_type size() const
static ::ExceptionBase & ExcNotImplemented()
value_type el(const unsigned int row, const unsigned int col) const
Definition: operators.h:1138
std::vector< UnivariateShapeData< Number > > data
Definition: shape_info.h:386
bool is_element(const size_type index) const
Definition: index_set.h:1765
void initialize(const OperatorType &operator_in)
Definition: operators.h:1699
SmartPointer< const OperatorType > mf_base_operator
Definition: operators.h:591
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1729
void apply(const AlignedVector< VectorizedArrayType > &inverse_coefficient, const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition: operators.h:1036
std::shared_ptr< DiagonalMatrix< VectorType > > inverse_diagonal_entries
Definition: operators.h:445
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1336
bool non_negative(const VectorizedArrayType &n)
Definition: operators.h:2042
virtual void compute_diagonal() override
Definition: operators.h:1790
void compute_diagonal(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, LinearAlgebra::distributed::Vector< Number > &diagonal_global, const std::function< void(FEEvaluation< dim, fe_degree, n_q_points_1d, n_components, Number, VectorizedArrayType > &)> &local_vmult, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
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:1857
std::vector< std::vector< unsigned int > > edge_constrained_indices
Definition: operators.h:463
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
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:1178
void run(const Iterator &begin, const typename identity< Iterator >::type &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
Definition: work_stream.h:691
void preprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1413
void transform_from_q_points_to_basis(const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition: operators.h:1063
void set_coefficient(const std::shared_ptr< Table< 2, VectorizedArrayType >> &scalar_coefficient)
Definition: operators.h:1940
static ::ExceptionBase & ExcInternalError()
void adjust_ghost_range_if_necessary(const VectorType &vec, const bool is_row) const
Definition: operators.h:1371