Reference documentation for deal.II version GIT 5dcc62ab57 2022-07-04 21:05:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
operators.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2011 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii_matrix_free_operators_h
18 #define dealii_matrix_free_operators_h
19 
20 
21 #include <deal.II/base/config.h>
22 
26 
29 
33 
35 
36 
38 
39 
41 {
42  namespace BlockHelper
43  {
44  // workaround for unifying non-block vector and block vector implementations
45  // a non-block vector has one block and the only subblock is the vector
46  // itself
47  template <typename VectorType>
48  typename std::enable_if<IsBlockVector<VectorType>::value,
49  unsigned int>::type
50  n_blocks(const VectorType &vector)
51  {
52  return vector.n_blocks();
53  }
54 
55  template <typename VectorType>
56  typename std::enable_if<!IsBlockVector<VectorType>::value,
57  unsigned int>::type
58  n_blocks(const VectorType &)
59  {
60  return 1;
61  }
62 
63  template <typename VectorType>
64  typename std::enable_if<IsBlockVector<VectorType>::value,
65  typename VectorType::BlockType &>::type
66  subblock(VectorType &vector, unsigned int block_no)
67  {
68  return vector.block(block_no);
69  }
70 
71  template <typename VectorType>
72  typename std::enable_if<IsBlockVector<VectorType>::value,
73  const typename VectorType::BlockType &>::type
74  subblock(const VectorType &vector, unsigned int block_no)
75  {
76  AssertIndexRange(block_no, vector.n_blocks());
77  return vector.block(block_no);
78  }
79 
80  template <typename VectorType>
81  typename std::enable_if<!IsBlockVector<VectorType>::value,
82  VectorType &>::type
83  subblock(VectorType &vector, unsigned int)
84  {
85  return vector;
86  }
87 
88  template <typename VectorType>
89  typename std::enable_if<!IsBlockVector<VectorType>::value,
90  const VectorType &>::type
91  subblock(const VectorType &vector, unsigned int)
92  {
93  return vector;
94  }
95 
96  template <typename VectorType>
97  typename std::enable_if<IsBlockVector<VectorType>::value, void>::type
98  collect_sizes(VectorType &vector)
99  {
100  vector.collect_sizes();
101  }
102 
103  template <typename VectorType>
104  typename std::enable_if<!IsBlockVector<VectorType>::value, void>::type
105  collect_sizes(const VectorType &)
106  {}
107  } // namespace BlockHelper
108 
186  template <int dim,
187  typename VectorType = LinearAlgebra::distributed::Vector<double>,
188  typename VectorizedArrayType =
190  class Base : public Subscriptor
191  {
192  public:
196  using value_type = typename VectorType::value_type;
197 
202 
206  Base();
207 
211  virtual ~Base() override = default;
212 
217  virtual void
218  clear();
219 
236  void
237  initialize(std::shared_ptr<
239  const std::vector<unsigned int> &selected_row_blocks =
240  std::vector<unsigned int>(),
241  const std::vector<unsigned int> &selected_column_blocks =
242  std::vector<unsigned int>());
243 
257  void
258  initialize(std::shared_ptr<
260  const MGConstrainedDoFs & mg_constrained_dofs,
261  const unsigned int level,
262  const std::vector<unsigned int> &selected_row_blocks =
263  std::vector<unsigned int>());
264 
279  void
280  initialize(std::shared_ptr<
282  const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
283  const unsigned int level,
284  const std::vector<unsigned int> & selected_row_blocks =
285  std::vector<unsigned int>());
286 
290  size_type
291  m() const;
292 
296  size_type
297  n() const;
298 
302  void
303  vmult_interface_down(VectorType &dst, const VectorType &src) const;
304 
308  void
309  vmult_interface_up(VectorType &dst, const VectorType &src) const;
310 
314  void
315  vmult(VectorType &dst, const VectorType &src) const;
316 
320  void
321  Tvmult(VectorType &dst, const VectorType &src) const;
322 
326  void
327  vmult_add(VectorType &dst, const VectorType &src) const;
328 
332  void
333  Tvmult_add(VectorType &dst, const VectorType &src) const;
334 
339  value_type
340  el(const unsigned int row, const unsigned int col) const;
341 
346  virtual std::size_t
348 
352  void
353  initialize_dof_vector(VectorType &vec) const;
354 
367  virtual void
369 
373  std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
375 
379  const std::shared_ptr<DiagonalMatrix<VectorType>> &
381 
385  const std::shared_ptr<DiagonalMatrix<VectorType>> &
387 
393  void
394  precondition_Jacobi(VectorType & dst,
395  const VectorType &src,
396  const value_type omega) const;
397 
398  protected:
403  void
404  preprocess_constraints(VectorType &dst, const VectorType &src) const;
405 
410  void
411  postprocess_constraints(VectorType &dst, const VectorType &src) const;
412 
417  void
418  set_constrained_entries_to_one(VectorType &dst) const;
419 
423  virtual void
424  apply_add(VectorType &dst, const VectorType &src) const = 0;
425 
431  virtual void
432  Tapply_add(VectorType &dst, const VectorType &src) const;
433 
437  std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
439 
444  std::shared_ptr<DiagonalMatrix<VectorType>> diagonal_entries;
445 
450  std::shared_ptr<DiagonalMatrix<VectorType>> inverse_diagonal_entries;
451 
456  std::vector<unsigned int> selected_rows;
457 
462  std::vector<unsigned int> selected_columns;
463 
464  private:
468  std::vector<std::vector<unsigned int>> edge_constrained_indices;
469 
473  mutable std::vector<std::vector<std::pair<value_type, value_type>>>
475 
481 
486  void
487  mult_add(VectorType & dst,
488  const VectorType &src,
489  const bool transpose) const;
490 
498  void
499  adjust_ghost_range_if_necessary(const VectorType &vec,
500  const bool is_row) const;
501  };
502 
503 
504 
539  template <typename OperatorType>
541  {
542  public:
546  using value_type = typename OperatorType::value_type;
547 
552 
557 
561  void
562  clear();
563 
567  void
568  initialize(const OperatorType &operator_in);
569 
573  template <typename VectorType>
574  void
575  vmult(VectorType &dst, const VectorType &src) const;
576 
580  template <typename VectorType>
581  void
582  Tvmult(VectorType &dst, const VectorType &src) const;
583 
587  template <typename VectorType>
588  void
589  initialize_dof_vector(VectorType &vec) const;
590 
591 
592  private:
597  };
598 
599 
600 
618  template <int dim,
619  int fe_degree,
620  int n_components = 1,
621  typename Number = double,
622  typename VectorizedArrayType = VectorizedArray<Number>>
624  {
625  static_assert(
626  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
627  "Type of Number and of VectorizedArrayType do not match.");
628 
629  public:
635  const FEEvaluationBase<dim,
636  n_components,
637  Number,
638  false,
639  VectorizedArrayType> &fe_eval);
640 
649  void
650  apply(const AlignedVector<VectorizedArrayType> &inverse_coefficient,
651  const unsigned int n_actual_components,
652  const VectorizedArrayType * in_array,
653  VectorizedArrayType * out_array) const;
654 
666  void
667  apply(const VectorizedArrayType *in_array,
668  VectorizedArrayType * out_array) const;
669 
703  void
704  transform_from_q_points_to_basis(const unsigned int n_actual_components,
705  const VectorizedArrayType *in_array,
706  VectorizedArrayType *out_array) const;
707 
713  void
715  AlignedVector<VectorizedArrayType> &inverse_jxw) const;
716 
717  private:
721  const FEEvaluationBase<dim,
722  n_components,
723  Number,
724  false,
725  VectorizedArrayType> &fe_eval;
726  };
727 
728 
729 
737  template <int dim,
738  int fe_degree,
739  int n_q_points_1d = fe_degree + 1,
740  int n_components = 1,
741  typename VectorType = LinearAlgebra::distributed::Vector<double>,
742  typename VectorizedArrayType =
744  class MassOperator : public Base<dim, VectorType, VectorizedArrayType>
745  {
746  public:
750  using value_type =
752 
756  using size_type =
758 
762  MassOperator();
763 
767  virtual void
768  compute_diagonal() override;
769 
785  void
787 
791  const std::shared_ptr<DiagonalMatrix<VectorType>> &
793 
797  const std::shared_ptr<DiagonalMatrix<VectorType>> &
799 
800  private:
806  virtual void
807  apply_add(VectorType &dst, const VectorType &src) const override;
808 
812  void
815  VectorType & dst,
816  const VectorType & src,
817  const std::pair<unsigned int, unsigned int> & cell_range) const;
818 
823  std::shared_ptr<DiagonalMatrix<VectorType>> lumped_diagonal_entries;
824 
829  std::shared_ptr<DiagonalMatrix<VectorType>> inverse_lumped_diagonal_entries;
830  };
831 
832 
833 
844  template <int dim,
845  int fe_degree,
846  int n_q_points_1d = fe_degree + 1,
847  int n_components = 1,
848  typename VectorType = LinearAlgebra::distributed::Vector<double>,
849  typename VectorizedArrayType =
851  class LaplaceOperator : public Base<dim, VectorType, VectorizedArrayType>
852  {
853  public:
857  using value_type =
859 
863  using size_type =
865 
869  LaplaceOperator();
870 
877  virtual void
878  compute_diagonal() override;
879 
930  void
932  const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient);
933 
938  virtual void
939  clear() override;
940 
947  std::shared_ptr<Table<2, VectorizedArrayType>>
948  get_coefficient();
949 
950  private:
956  virtual void
957  apply_add(VectorType &dst, const VectorType &src) const override;
958 
962  void
965  VectorType & dst,
966  const VectorType & src,
967  const std::pair<unsigned int, unsigned int> & cell_range) const;
968 
972  void
975  VectorType & dst,
976  const VectorType &,
977  const std::pair<unsigned int, unsigned int> &cell_range) const;
978 
982  void
985  & phi,
986  const unsigned int cell) const;
987 
991  std::shared_ptr<Table<2, VectorizedArrayType>> scalar_coefficient;
992  };
993 
994 
995 
996  // ------------------------------------ inline functions ---------------------
997 
998  template <int dim,
999  int fe_degree,
1000  int n_components,
1001  typename Number,
1002  typename VectorizedArrayType>
1003  inline CellwiseInverseMassMatrix<dim,
1004  fe_degree,
1005  n_components,
1006  Number,
1007  VectorizedArrayType>::
1008  CellwiseInverseMassMatrix(
1009  const FEEvaluationBase<dim,
1010  n_components,
1011  Number,
1012  false,
1013  VectorizedArrayType> &fe_eval)
1014  : fe_eval(fe_eval)
1015  {}
1016 
1017 
1018 
1019  template <int dim,
1020  int fe_degree,
1021  int n_components,
1022  typename Number,
1023  typename VectorizedArrayType>
1024  inline void
1026  fe_degree,
1027  n_components,
1028  Number,
1029  VectorizedArrayType>::
1030  fill_inverse_JxW_values(
1031  AlignedVector<VectorizedArrayType> &inverse_jxw) const
1032  {
1033  constexpr unsigned int dofs_per_component_on_cell =
1034  Utilities::pow(fe_degree + 1, dim);
1035  Assert(inverse_jxw.size() > 0 &&
1036  inverse_jxw.size() % dofs_per_component_on_cell == 0,
1037  ExcMessage(
1038  "Expected diagonal to be a multiple of scalar dof per cells"));
1039 
1040  // compute values for the first component
1041  for (unsigned int q = 0; q < dofs_per_component_on_cell; ++q)
1042  inverse_jxw[q] = 1. / fe_eval.JxW(q);
1043  // copy values to rest of vector
1044  for (unsigned int q = dofs_per_component_on_cell; q < inverse_jxw.size();)
1045  for (unsigned int i = 0; i < dofs_per_component_on_cell; ++i, ++q)
1046  inverse_jxw[q] = inverse_jxw[i];
1047  }
1048 
1049 
1050 
1051  template <int dim,
1052  int fe_degree,
1053  int n_components,
1054  typename Number,
1055  typename VectorizedArrayType>
1056  inline void
1058  dim,
1059  fe_degree,
1060  n_components,
1061  Number,
1062  VectorizedArrayType>::apply(const VectorizedArrayType *in_array,
1063  VectorizedArrayType * out_array) const
1064  {
1065  if (fe_degree > -1)
1067  template run<fe_degree>(n_components, fe_eval, in_array, out_array);
1068  else
1070  n_components, fe_eval, in_array, out_array);
1071  }
1072 
1073 
1074 
1075  template <int dim,
1076  int fe_degree,
1077  int n_components,
1078  typename Number,
1079  typename VectorizedArrayType>
1080  inline void
1082  fe_degree,
1083  n_components,
1084  Number,
1085  VectorizedArrayType>::
1086  apply(const AlignedVector<VectorizedArrayType> &inverse_coefficients,
1087  const unsigned int n_actual_components,
1088  const VectorizedArrayType * in_array,
1089  VectorizedArrayType * out_array) const
1090  {
1091  const unsigned int given_degree =
1092  fe_eval.get_shape_info().data[0].fe_degree;
1093  if (fe_degree > -1)
1095  VectorizedArrayType>::
1096  template run<fe_degree>(
1097  n_actual_components,
1098  fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
1099  inverse_coefficients,
1100  in_array,
1101  out_array);
1102  else
1104  n_actual_components,
1105  given_degree,
1106  fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
1107  inverse_coefficients,
1108  in_array,
1109  out_array);
1110  }
1111 
1112 
1113 
1114  template <int dim,
1115  int fe_degree,
1116  int n_components,
1117  typename Number,
1118  typename VectorizedArrayType>
1119  inline void
1121  fe_degree,
1122  n_components,
1123  Number,
1124  VectorizedArrayType>::
1125  transform_from_q_points_to_basis(const unsigned int n_actual_components,
1126  const VectorizedArrayType *in_array,
1127  VectorizedArrayType * out_array) const
1128  {
1129  const auto n_q_points_1d = fe_eval.get_shape_info().data[0].n_q_points_1d;
1130 
1131  if (fe_degree > -1 && (fe_degree + 1 == n_q_points_1d))
1133  dim,
1134  VectorizedArrayType>::template run<fe_degree,
1135  fe_degree + 1>(n_actual_components,
1136  fe_eval,
1137  in_array,
1138  out_array);
1139  else
1141  transform_from_q_points_to_basis(n_actual_components,
1142  fe_eval,
1143  in_array,
1144  out_array);
1145  }
1146 
1147 
1148 
1149  //----------------- Base operator -----------------------------
1150  template <int dim, typename VectorType, typename VectorizedArrayType>
1152  : Subscriptor()
1153  , have_interface_matrices(false)
1154  {}
1155 
1156 
1157 
1158  template <int dim, typename VectorType, typename VectorizedArrayType>
1161  {
1162  Assert(data.get() != nullptr, ExcNotInitialized());
1164  0;
1165  for (const unsigned int selected_row : selected_rows)
1166  total_size += data->get_vector_partitioner(selected_row)->size();
1167  return total_size;
1168  }
1169 
1170 
1171 
1172  template <int dim, typename VectorType, typename VectorizedArrayType>
1175  {
1176  Assert(data.get() != nullptr, ExcNotInitialized());
1178  0;
1179  for (const unsigned int selected_column : selected_columns)
1180  total_size += data->get_vector_partitioner(selected_column)->size();
1181  return total_size;
1182  }
1183 
1184 
1185 
1186  template <int dim, typename VectorType, typename VectorizedArrayType>
1187  void
1189  {
1190  data.reset();
1191  inverse_diagonal_entries.reset();
1192  }
1193 
1194 
1195 
1196  template <int dim, typename VectorType, typename VectorizedArrayType>
1199  const unsigned int col) const
1200  {
1201  (void)col;
1202  Assert(row == col, ExcNotImplemented());
1203  Assert(inverse_diagonal_entries.get() != nullptr &&
1204  inverse_diagonal_entries->m() > 0,
1205  ExcNotInitialized());
1206  return 1.0 / (*inverse_diagonal_entries)(row, row);
1207  }
1208 
1209 
1210 
1211  template <int dim, typename VectorType, typename VectorizedArrayType>
1212  void
1214  VectorType &vec) const
1215  {
1216  Assert(data.get() != nullptr, ExcNotInitialized());
1217  AssertDimension(BlockHelper::n_blocks(vec), selected_rows.size());
1218  for (unsigned int i = 0; i < BlockHelper::n_blocks(vec); ++i)
1219  {
1220  const unsigned int index = selected_rows[i];
1221  if (!BlockHelper::subblock(vec, index)
1222  .partitioners_are_compatible(
1223  *data->get_dof_info(index).vector_partitioner))
1224  data->initialize_dof_vector(BlockHelper::subblock(vec, index), index);
1225 
1226  Assert(BlockHelper::subblock(vec, index)
1227  .partitioners_are_globally_compatible(
1228  *data->get_dof_info(index).vector_partitioner),
1229  ExcInternalError());
1230  }
1232  }
1233 
1234 
1235 
1236  template <int dim, typename VectorType, typename VectorizedArrayType>
1237  void
1240  data_,
1241  const std::vector<unsigned int> &given_row_selection,
1242  const std::vector<unsigned int> &given_column_selection)
1243  {
1244  data = data_;
1245 
1246  selected_rows.clear();
1247  selected_columns.clear();
1248  if (given_row_selection.empty())
1249  for (unsigned int i = 0; i < data_->n_components(); ++i)
1250  selected_rows.push_back(i);
1251  else
1252  {
1253  for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1254  {
1255  AssertIndexRange(given_row_selection[i], data_->n_components());
1256  for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1257  if (j != i)
1258  Assert(given_row_selection[j] != given_row_selection[i],
1259  ExcMessage("Given row indices must be unique"));
1260 
1261  selected_rows.push_back(given_row_selection[i]);
1262  }
1263  }
1264  if (given_column_selection.size() == 0)
1265  selected_columns = selected_rows;
1266  else
1267  {
1268  for (unsigned int i = 0; i < given_column_selection.size(); ++i)
1269  {
1270  AssertIndexRange(given_column_selection[i], data_->n_components());
1271  for (unsigned int j = 0; j < given_column_selection.size(); ++j)
1272  if (j != i)
1273  Assert(given_column_selection[j] != given_column_selection[i],
1274  ExcMessage("Given column indices must be unique"));
1275 
1276  selected_columns.push_back(given_column_selection[i]);
1277  }
1278  }
1279 
1280  edge_constrained_indices.clear();
1281  edge_constrained_indices.resize(selected_rows.size());
1282  edge_constrained_values.clear();
1283  edge_constrained_values.resize(selected_rows.size());
1284  have_interface_matrices = false;
1285  }
1286 
1287 
1288 
1289  template <int dim, typename VectorType, typename VectorizedArrayType>
1290  void
1293  data_,
1294  const MGConstrainedDoFs & mg_constrained_dofs,
1295  const unsigned int level,
1296  const std::vector<unsigned int> &given_row_selection)
1297  {
1298  std::vector<MGConstrainedDoFs> mg_constrained_dofs_vector(
1299  1, mg_constrained_dofs);
1300  initialize(data_, mg_constrained_dofs_vector, level, given_row_selection);
1301  }
1302 
1303 
1304 
1305  template <int dim, typename VectorType, typename VectorizedArrayType>
1306  void
1309  data_,
1310  const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
1311  const unsigned int level,
1312  const std::vector<unsigned int> & given_row_selection)
1313  {
1315  ExcMessage("level is not set"));
1316 
1317  selected_rows.clear();
1318  selected_columns.clear();
1319  if (given_row_selection.empty())
1320  for (unsigned int i = 0; i < data_->n_components(); ++i)
1321  selected_rows.push_back(i);
1322  else
1323  {
1324  for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1325  {
1326  AssertIndexRange(given_row_selection[i], data_->n_components());
1327  for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1328  if (j != i)
1329  Assert(given_row_selection[j] != given_row_selection[i],
1330  ExcMessage("Given row indices must be unique"));
1331 
1332  selected_rows.push_back(given_row_selection[i]);
1333  }
1334  }
1335  selected_columns = selected_rows;
1336 
1337  AssertDimension(mg_constrained_dofs.size(), selected_rows.size());
1338  edge_constrained_indices.clear();
1339  edge_constrained_indices.resize(selected_rows.size());
1340  edge_constrained_values.clear();
1341  edge_constrained_values.resize(selected_rows.size());
1342 
1343  data = data_;
1344 
1345  for (unsigned int j = 0; j < selected_rows.size(); ++j)
1346  {
1347  if (data_->n_cell_batches() > 0)
1348  {
1349  AssertDimension(level, data_->get_cell_iterator(0, 0, j)->level());
1350  }
1351 
1352  // setup edge_constrained indices
1353  std::vector<types::global_dof_index> interface_indices;
1354  mg_constrained_dofs[j]
1355  .get_refinement_edge_indices(level)
1356  .fill_index_vector(interface_indices);
1357  edge_constrained_indices[j].clear();
1358  edge_constrained_indices[j].reserve(interface_indices.size());
1359  edge_constrained_values[j].resize(interface_indices.size());
1360  const IndexSet &locally_owned =
1361  data->get_dof_handler(selected_rows[j]).locally_owned_mg_dofs(level);
1362  for (const auto interface_index : interface_indices)
1363  if (locally_owned.is_element(interface_index))
1364  edge_constrained_indices[j].push_back(
1365  locally_owned.index_within_set(interface_index));
1366  have_interface_matrices |=
1368  static_cast<unsigned int>(edge_constrained_indices[j].size()),
1369  data->get_vector_partitioner()->get_mpi_communicator()) > 0;
1370  }
1371  }
1372 
1373 
1374 
1375  template <int dim, typename VectorType, typename VectorizedArrayType>
1376  void
1378  VectorType &dst) const
1379  {
1380  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1381  {
1382  const std::vector<unsigned int> &constrained_dofs =
1383  data->get_constrained_dofs(selected_rows[j]);
1384  for (const auto constrained_dof : constrained_dofs)
1385  BlockHelper::subblock(dst, j).local_element(constrained_dof) = 1.;
1386  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1387  BlockHelper::subblock(dst, j).local_element(
1388  edge_constrained_indices[j][i]) = 1.;
1389  }
1390  }
1391 
1392 
1393 
1394  template <int dim, typename VectorType, typename VectorizedArrayType>
1395  void
1397  const VectorType &src) const
1398  {
1399  using Number =
1401  dst = Number(0.);
1402  vmult_add(dst, src);
1403  }
1404 
1405 
1406 
1407  template <int dim, typename VectorType, typename VectorizedArrayType>
1408  void
1410  VectorType & dst,
1411  const VectorType &src) const
1412  {
1413  mult_add(dst, src, false);
1414  }
1415 
1416 
1417 
1418  template <int dim, typename VectorType, typename VectorizedArrayType>
1419  void
1421  VectorType & dst,
1422  const VectorType &src) const
1423  {
1424  mult_add(dst, src, true);
1425  }
1426 
1427 
1428 
1429  template <int dim, typename VectorType, typename VectorizedArrayType>
1430  void
1432  const VectorType &src,
1433  const bool is_row) const
1434  {
1435  using Number =
1437  for (unsigned int i = 0; i < BlockHelper::n_blocks(src); ++i)
1438  {
1439  const unsigned int mf_component =
1440  is_row ? selected_rows[i] : selected_columns[i];
1441  // If both vectors use the same partitioner -> done
1442  if (BlockHelper::subblock(src, i).get_partitioner().get() ==
1443  data->get_dof_info(mf_component).vector_partitioner.get())
1444  continue;
1445 
1446  // If not, assert that the local ranges are the same and reset to the
1447  // current partitioner
1449  .get_partitioner()
1450  ->locally_owned_size() ==
1451  data->get_dof_info(mf_component)
1452  .vector_partitioner->locally_owned_size(),
1453  ExcMessage(
1454  "The vector passed to the vmult() function does not have "
1455  "the correct size for compatibility with MatrixFree."));
1456 
1457  // copy the vector content to a temporary vector so that it does not get
1458  // lost
1460  BlockHelper::subblock(src, i));
1461  this->data->initialize_dof_vector(
1462  BlockHelper::subblock(const_cast<VectorType &>(src), i),
1463  mf_component);
1464  BlockHelper::subblock(const_cast<VectorType &>(src), i)
1465  .copy_locally_owned_data_from(copy_vec);
1466  }
1467  }
1468 
1469 
1470 
1471  template <int dim, typename VectorType, typename VectorizedArrayType>
1472  void
1474  VectorType & dst,
1475  const VectorType &src) const
1476  {
1477  using Number =
1479  adjust_ghost_range_if_necessary(src, false);
1480  adjust_ghost_range_if_necessary(dst, true);
1481 
1482  // set zero Dirichlet values on the input vector (and remember the src and
1483  // dst values because we need to reset them at the end)
1484  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1485  {
1486  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1487  {
1488  edge_constrained_values[j][i] = std::pair<Number, Number>(
1489  BlockHelper::subblock(src, j).local_element(
1490  edge_constrained_indices[j][i]),
1491  BlockHelper::subblock(dst, j).local_element(
1492  edge_constrained_indices[j][i]));
1493  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1494  .local_element(edge_constrained_indices[j][i]) = 0.;
1495  }
1496  }
1497  }
1498 
1499 
1500 
1501  template <int dim, typename VectorType, typename VectorizedArrayType>
1502  void
1504  VectorType & dst,
1505  const VectorType &src,
1506  const bool transpose) const
1507  {
1508  AssertDimension(dst.size(), src.size());
1510  AssertDimension(BlockHelper::n_blocks(dst), selected_rows.size());
1511  preprocess_constraints(dst, src);
1512  if (transpose)
1513  Tapply_add(dst, src);
1514  else
1515  apply_add(dst, src);
1516  postprocess_constraints(dst, src);
1517  }
1518 
1519 
1520 
1521  template <int dim, typename VectorType, typename VectorizedArrayType>
1522  void
1524  VectorType & dst,
1525  const VectorType &src) const
1526  {
1527  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1528  {
1529  const std::vector<unsigned int> &constrained_dofs =
1530  data->get_constrained_dofs(selected_rows[j]);
1531  for (const auto constrained_dof : constrained_dofs)
1532  BlockHelper::subblock(dst, j).local_element(constrained_dof) +=
1533  BlockHelper::subblock(src, j).local_element(constrained_dof);
1534  }
1535 
1536  // reset edge constrained values, multiply by unit matrix and add into
1537  // destination
1538  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1539  {
1540  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1541  {
1542  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1543  .local_element(edge_constrained_indices[j][i]) =
1544  edge_constrained_values[j][i].first;
1545  BlockHelper::subblock(dst, j).local_element(
1546  edge_constrained_indices[j][i]) =
1547  edge_constrained_values[j][i].second +
1548  edge_constrained_values[j][i].first;
1549  }
1550  }
1551  }
1552 
1553 
1554 
1555  template <int dim, typename VectorType, typename VectorizedArrayType>
1556  void
1558  VectorType & dst,
1559  const VectorType &src) const
1560  {
1561  using Number =
1563  AssertDimension(dst.size(), src.size());
1564  adjust_ghost_range_if_necessary(src, false);
1565  adjust_ghost_range_if_necessary(dst, true);
1566 
1567  dst = Number(0.);
1568 
1569  if (!have_interface_matrices)
1570  return;
1571 
1572  // set zero Dirichlet values on the input vector (and remember the src and
1573  // dst values because we need to reset them at the end)
1574  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1575  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1576  {
1577  edge_constrained_values[j][i] = std::pair<Number, Number>(
1578  BlockHelper::subblock(src, j).local_element(
1579  edge_constrained_indices[j][i]),
1580  BlockHelper::subblock(dst, j).local_element(
1581  edge_constrained_indices[j][i]));
1582  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1583  .local_element(edge_constrained_indices[j][i]) = 0.;
1584  }
1585 
1586  apply_add(dst, src);
1587 
1588  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1589  {
1590  unsigned int c = 0;
1591  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1592  {
1593  for (; c < edge_constrained_indices[j][i]; ++c)
1594  BlockHelper::subblock(dst, j).local_element(c) = 0.;
1595  ++c;
1596 
1597  // reset the src values
1598  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1599  .local_element(edge_constrained_indices[j][i]) =
1600  edge_constrained_values[j][i].first;
1601  }
1602  for (; c < BlockHelper::subblock(dst, j).locally_owned_size(); ++c)
1603  BlockHelper::subblock(dst, j).local_element(c) = 0.;
1604  }
1605  }
1606 
1607 
1608 
1609  template <int dim, typename VectorType, typename VectorizedArrayType>
1610  void
1612  VectorType & dst,
1613  const VectorType &src) const
1614  {
1615  using Number =
1617  AssertDimension(dst.size(), src.size());
1618  adjust_ghost_range_if_necessary(src, false);
1619  adjust_ghost_range_if_necessary(dst, true);
1620 
1621  dst = Number(0.);
1622 
1623  if (!have_interface_matrices)
1624  return;
1625 
1626  VectorType src_cpy(src);
1627  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1628  {
1629  unsigned int c = 0;
1630  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1631  {
1632  for (; c < edge_constrained_indices[j][i]; ++c)
1633  BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1634  ++c;
1635  }
1636  for (; c < BlockHelper::subblock(src_cpy, j).locally_owned_size(); ++c)
1637  BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1638  }
1639 
1640  apply_add(dst, src_cpy);
1641 
1642  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1643  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1644  BlockHelper::subblock(dst, j).local_element(
1645  edge_constrained_indices[j][i]) = 0.;
1646  }
1647 
1648 
1649 
1650  template <int dim, typename VectorType, typename VectorizedArrayType>
1651  void
1653  VectorType & dst,
1654  const VectorType &src) const
1655  {
1656  using Number =
1658  dst = Number(0.);
1659  Tvmult_add(dst, src);
1660  }
1661 
1662 
1663 
1664  template <int dim, typename VectorType, typename VectorizedArrayType>
1665  std::size_t
1667  {
1668  return inverse_diagonal_entries.get() != nullptr ?
1669  inverse_diagonal_entries->memory_consumption() :
1670  sizeof(*this);
1671  }
1672 
1673 
1674 
1675  template <int dim, typename VectorType, typename VectorizedArrayType>
1676  std::shared_ptr<const MatrixFree<
1677  dim,
1679  VectorizedArrayType>>
1681  {
1682  return data;
1683  }
1684 
1685 
1686 
1687  template <int dim, typename VectorType, typename VectorizedArrayType>
1688  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1690  const
1691  {
1692  Assert(inverse_diagonal_entries.get() != nullptr &&
1693  inverse_diagonal_entries->m() > 0,
1694  ExcNotInitialized());
1695  return inverse_diagonal_entries;
1696  }
1697 
1698 
1699 
1700  template <int dim, typename VectorType, typename VectorizedArrayType>
1701  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1703  {
1704  Assert(diagonal_entries.get() != nullptr && diagonal_entries->m() > 0,
1705  ExcNotInitialized());
1706  return diagonal_entries;
1707  }
1708 
1709 
1710 
1711  template <int dim, typename VectorType, typename VectorizedArrayType>
1712  void
1714  VectorType & dst,
1715  const VectorType &src) const
1716  {
1717  apply_add(dst, src);
1718  }
1719 
1720 
1721 
1722  template <int dim, typename VectorType, typename VectorizedArrayType>
1723  void
1725  VectorType & dst,
1726  const VectorType & src,
1728  const
1729  {
1730  Assert(inverse_diagonal_entries.get() && inverse_diagonal_entries->m() > 0,
1731  ExcNotInitialized());
1732  inverse_diagonal_entries->vmult(dst, src);
1733  dst *= omega;
1734  }
1735 
1736 
1737 
1738  //------------------------- MGInterfaceOperator ------------------------------
1739 
1740  template <typename OperatorType>
1742  : Subscriptor()
1743  , mf_base_operator(nullptr)
1744  {}
1745 
1746 
1747 
1748  template <typename OperatorType>
1749  void
1751  {
1752  mf_base_operator = nullptr;
1753  }
1754 
1755 
1756 
1757  template <typename OperatorType>
1758  void
1759  MGInterfaceOperator<OperatorType>::initialize(const OperatorType &operator_in)
1760  {
1761  mf_base_operator = &operator_in;
1762  }
1763 
1764 
1765 
1766  template <typename OperatorType>
1767  template <typename VectorType>
1768  void
1770  const VectorType &src) const
1771  {
1772 #ifndef DEAL_II_MSVC
1773  static_assert(
1774  std::is_same<typename VectorType::value_type, value_type>::value,
1775  "The vector type must be based on the same value type as this "
1776  "operator");
1777 #endif
1778 
1779  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1780 
1781  mf_base_operator->vmult_interface_down(dst, src);
1782  }
1783 
1784 
1785 
1786  template <typename OperatorType>
1787  template <typename VectorType>
1788  void
1790  const VectorType &src) const
1791  {
1792 #ifndef DEAL_II_MSVC
1793  static_assert(
1794  std::is_same<typename VectorType::value_type, value_type>::value,
1795  "The vector type must be based on the same value type as this "
1796  "operator");
1797 #endif
1798 
1799  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1800 
1801  mf_base_operator->vmult_interface_up(dst, src);
1802  }
1803 
1804 
1805 
1806  template <typename OperatorType>
1807  template <typename VectorType>
1808  void
1810  VectorType &vec) const
1811  {
1812  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1813 
1814  mf_base_operator->initialize_dof_vector(vec);
1815  }
1816 
1817 
1818 
1819  //-----------------------------MassOperator----------------------------------
1820 
1821  template <int dim,
1822  int fe_degree,
1823  int n_q_points_1d,
1824  int n_components,
1825  typename VectorType,
1826  typename VectorizedArrayType>
1827  MassOperator<dim,
1828  fe_degree,
1829  n_q_points_1d,
1830  n_components,
1831  VectorType,
1832  VectorizedArrayType>::MassOperator()
1833  : Base<dim, VectorType, VectorizedArrayType>()
1834  {}
1835 
1836 
1837 
1838  template <int dim,
1839  int fe_degree,
1840  int n_q_points_1d,
1841  int n_components,
1842  typename VectorType,
1843  typename VectorizedArrayType>
1844  void
1845  MassOperator<dim,
1846  fe_degree,
1847  n_q_points_1d,
1848  n_components,
1849  VectorType,
1850  VectorizedArrayType>::compute_diagonal()
1851  {
1852  using Number =
1855  ExcNotInitialized());
1856  Assert(this->selected_rows == this->selected_columns,
1857  ExcMessage("This function is only implemented for square (not "
1858  "rectangular) operators."));
1859 
1860  this->inverse_diagonal_entries =
1861  std::make_shared<DiagonalMatrix<VectorType>>();
1862  this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1863  VectorType &inverse_diagonal_vector =
1864  this->inverse_diagonal_entries->get_vector();
1865  VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1866  this->initialize_dof_vector(inverse_diagonal_vector);
1867  this->initialize_dof_vector(diagonal_vector);
1868 
1869  // Set up the action of the mass matrix in a way that's compatible with
1870  // MatrixFreeTools::compute_diagonal:
1871  auto diagonal_evaluation = [](auto &integrator) {
1872  integrator.evaluate(EvaluationFlags::values);
1873  for (unsigned int q = 0; q < integrator.n_q_points; ++q)
1874  integrator.submit_value(integrator.get_value(q), q);
1875  integrator.integrate(EvaluationFlags::values);
1876  };
1877 
1878  std::function<void(
1879  FEEvaluation<
1880  dim,
1881  fe_degree,
1882  n_q_points_1d,
1883  n_components,
1885  VectorizedArrayType> &)>
1886  diagonal_evaluation_f(diagonal_evaluation);
1887 
1888  Assert(this->selected_rows.size() > 0, ExcInternalError());
1889  for (unsigned int block_n = 0; block_n < this->selected_rows.size();
1890  ++block_n)
1892  BlockHelper::subblock(diagonal_vector,
1893  block_n),
1894  diagonal_evaluation_f,
1895  this->selected_rows[block_n]);
1896 
1897  // Constrained entries will create zeros on the main diagonal, which we
1898  // don't want
1899  this->set_constrained_entries_to_one(diagonal_vector);
1900 
1901  inverse_diagonal_vector = diagonal_vector;
1902 
1903  for (unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
1904  ++i)
1905  {
1906  Assert(diagonal_vector.local_element(i) > Number(0),
1907  ExcInternalError());
1908  inverse_diagonal_vector.local_element(i) =
1909  1. / inverse_diagonal_vector.local_element(i);
1910  }
1911 
1912  // We never need ghost values so don't update them
1913  }
1914 
1915 
1916 
1917  template <int dim,
1918  int fe_degree,
1919  int n_q_points_1d,
1920  int n_components,
1921  typename VectorType,
1922  typename VectorizedArrayType>
1923  void
1924  MassOperator<dim,
1925  fe_degree,
1926  n_q_points_1d,
1927  n_components,
1928  VectorType,
1929  VectorizedArrayType>::compute_lumped_diagonal()
1930  {
1931  using Number =
1934  ExcNotInitialized());
1935  Assert(this->selected_rows == this->selected_columns,
1936  ExcMessage("This function is only implemented for square (not "
1937  "rectangular) operators."));
1938 
1939  inverse_lumped_diagonal_entries =
1940  std::make_shared<DiagonalMatrix<VectorType>>();
1941  lumped_diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1942  VectorType &inverse_lumped_diagonal_vector =
1943  inverse_lumped_diagonal_entries->get_vector();
1944  VectorType &lumped_diagonal_vector = lumped_diagonal_entries->get_vector();
1945  this->initialize_dof_vector(inverse_lumped_diagonal_vector);
1946  this->initialize_dof_vector(lumped_diagonal_vector);
1947 
1948  // Re-use the inverse_lumped_diagonal_vector as the vector of 1s
1949  inverse_lumped_diagonal_vector = Number(1.);
1950  apply_add(lumped_diagonal_vector, inverse_lumped_diagonal_vector);
1951  this->set_constrained_entries_to_one(lumped_diagonal_vector);
1952 
1953  const size_type locally_owned_size =
1954  inverse_lumped_diagonal_vector.locally_owned_size();
1955  // A caller may request a lumped diagonal matrix when it doesn't make sense
1956  // (e.g., an element with negative-mean basis functions). Avoid division by
1957  // zero so we don't cause a floating point exception but permit negative
1958  // entries here.
1959  for (size_type i = 0; i < locally_owned_size; ++i)
1960  {
1961  if (lumped_diagonal_vector.local_element(i) == Number(0.))
1962  inverse_lumped_diagonal_vector.local_element(i) = Number(1.);
1963  else
1964  inverse_lumped_diagonal_vector.local_element(i) =
1965  Number(1.) / lumped_diagonal_vector.local_element(i);
1966  }
1967 
1968  inverse_lumped_diagonal_vector.update_ghost_values();
1969  lumped_diagonal_vector.update_ghost_values();
1970  }
1971 
1972 
1973 
1974  template <int dim,
1975  int fe_degree,
1976  int n_q_points_1d,
1977  int n_components,
1978  typename VectorType,
1979  typename VectorizedArrayType>
1980  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1981  MassOperator<dim,
1982  fe_degree,
1983  n_q_points_1d,
1984  n_components,
1985  VectorType,
1986  VectorizedArrayType>::get_matrix_lumped_diagonal_inverse() const
1987  {
1988  Assert(inverse_lumped_diagonal_entries.get() != nullptr &&
1989  inverse_lumped_diagonal_entries->m() > 0,
1990  ExcNotInitialized());
1991  return inverse_lumped_diagonal_entries;
1992  }
1993 
1994 
1995 
1996  template <int dim,
1997  int fe_degree,
1998  int n_q_points_1d,
1999  int n_components,
2000  typename VectorType,
2001  typename VectorizedArrayType>
2002  const std::shared_ptr<DiagonalMatrix<VectorType>> &
2003  MassOperator<dim,
2004  fe_degree,
2005  n_q_points_1d,
2006  n_components,
2007  VectorType,
2008  VectorizedArrayType>::get_matrix_lumped_diagonal() const
2009  {
2010  Assert(lumped_diagonal_entries.get() != nullptr &&
2011  lumped_diagonal_entries->m() > 0,
2012  ExcNotInitialized());
2013  return lumped_diagonal_entries;
2014  }
2015 
2016 
2017 
2018  template <int dim,
2019  int fe_degree,
2020  int n_q_points_1d,
2021  int n_components,
2022  typename VectorType,
2023  typename VectorizedArrayType>
2024  void
2025  MassOperator<dim,
2026  fe_degree,
2027  n_q_points_1d,
2028  n_components,
2029  VectorType,
2030  VectorizedArrayType>::apply_add(VectorType & dst,
2031  const VectorType &src) const
2032  {
2034  &MassOperator::local_apply_cell, this, dst, src);
2035  }
2036 
2037 
2038 
2039  template <int dim,
2040  int fe_degree,
2041  int n_q_points_1d,
2042  int n_components,
2043  typename VectorType,
2044  typename VectorizedArrayType>
2045  void
2046  MassOperator<dim,
2047  fe_degree,
2048  n_q_points_1d,
2049  n_components,
2050  VectorType,
2051  VectorizedArrayType>::
2052  local_apply_cell(
2053  const MatrixFree<
2054  dim,
2056  VectorizedArrayType> & data,
2057  VectorType & dst,
2058  const VectorType & src,
2059  const std::pair<unsigned int, unsigned int> &cell_range) const
2060  {
2061  using Number =
2063  FEEvaluation<dim,
2064  fe_degree,
2065  n_q_points_1d,
2066  n_components,
2067  Number,
2068  VectorizedArrayType>
2069  phi(data, this->selected_rows[0]);
2070  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2071  {
2072  phi.reinit(cell);
2073  phi.read_dof_values(src);
2074  phi.evaluate(EvaluationFlags::values);
2075  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2076  phi.submit_value(phi.get_value(q), q);
2077  phi.integrate(EvaluationFlags::values);
2078  phi.distribute_local_to_global(dst);
2079  }
2080  }
2081 
2082 
2083  //-----------------------------LaplaceOperator----------------------------------
2084 
2085  template <int dim,
2086  int fe_degree,
2087  int n_q_points_1d,
2088  int n_components,
2089  typename VectorType,
2090  typename VectorizedArrayType>
2091  LaplaceOperator<dim,
2092  fe_degree,
2093  n_q_points_1d,
2094  n_components,
2095  VectorType,
2096  VectorizedArrayType>::LaplaceOperator()
2097  : Base<dim, VectorType, VectorizedArrayType>()
2098  {}
2099 
2100 
2101 
2102  template <int dim,
2103  int fe_degree,
2104  int n_q_points_1d,
2105  int n_components,
2106  typename VectorType,
2107  typename VectorizedArrayType>
2108  void
2109  LaplaceOperator<dim,
2110  fe_degree,
2111  n_q_points_1d,
2112  n_components,
2113  VectorType,
2114  VectorizedArrayType>::clear()
2115  {
2117  scalar_coefficient.reset();
2118  }
2119 
2120 
2121 
2122  template <int dim,
2123  int fe_degree,
2124  int n_q_points_1d,
2125  int n_components,
2126  typename VectorType,
2127  typename VectorizedArrayType>
2128  void
2129  LaplaceOperator<dim,
2130  fe_degree,
2131  n_q_points_1d,
2132  n_components,
2133  VectorType,
2134  VectorizedArrayType>::
2135  set_coefficient(
2136  const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient_)
2137  {
2138  scalar_coefficient = scalar_coefficient_;
2139  }
2140 
2141 
2142 
2143  template <int dim,
2144  int fe_degree,
2145  int n_q_points_1d,
2146  int n_components,
2147  typename VectorType,
2148  typename VectorizedArrayType>
2149  std::shared_ptr<Table<2, VectorizedArrayType>>
2150  LaplaceOperator<dim,
2151  fe_degree,
2152  n_q_points_1d,
2153  n_components,
2154  VectorType,
2155  VectorizedArrayType>::get_coefficient()
2156  {
2157  Assert(scalar_coefficient.get(), ExcNotInitialized());
2158  return scalar_coefficient;
2159  }
2160 
2161 
2162 
2163  template <int dim,
2164  int fe_degree,
2165  int n_q_points_1d,
2166  int n_components,
2167  typename VectorType,
2168  typename VectorizedArrayType>
2169  void
2170  LaplaceOperator<dim,
2171  fe_degree,
2172  n_q_points_1d,
2173  n_components,
2174  VectorType,
2175  VectorizedArrayType>::compute_diagonal()
2176  {
2177  using Number =
2180  ExcNotInitialized());
2181 
2182  this->inverse_diagonal_entries =
2183  std::make_shared<DiagonalMatrix<VectorType>>();
2184  this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
2185  VectorType &inverse_diagonal_vector =
2186  this->inverse_diagonal_entries->get_vector();
2187  VectorType &diagonal_vector = this->diagonal_entries->get_vector();
2188  this->initialize_dof_vector(inverse_diagonal_vector);
2189  this->initialize_dof_vector(diagonal_vector);
2190 
2191  this->data->cell_loop(&LaplaceOperator::local_diagonal_cell,
2192  this,
2193  diagonal_vector,
2194  /*unused*/ diagonal_vector);
2195  this->set_constrained_entries_to_one(diagonal_vector);
2196 
2197  inverse_diagonal_vector = diagonal_vector;
2198 
2199  for (unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
2200  ++i)
2201  if (std::abs(inverse_diagonal_vector.local_element(i)) >
2203  inverse_diagonal_vector.local_element(i) =
2204  1. / inverse_diagonal_vector.local_element(i);
2205  else
2206  inverse_diagonal_vector.local_element(i) = 1.;
2207 
2208  // We never need ghost values so don't update them
2209  }
2210 
2211 
2212 
2213  template <int dim,
2214  int fe_degree,
2215  int n_q_points_1d,
2216  int n_components,
2217  typename VectorType,
2218  typename VectorizedArrayType>
2219  void
2220  LaplaceOperator<dim,
2221  fe_degree,
2222  n_q_points_1d,
2223  n_components,
2224  VectorType,
2225  VectorizedArrayType>::apply_add(VectorType & dst,
2226  const VectorType &src) const
2227  {
2229  &LaplaceOperator::local_apply_cell, this, dst, src);
2230  }
2231 
2232  namespace Implementation
2233  {
2234  template <typename VectorizedArrayType>
2235  bool
2236  non_negative(const VectorizedArrayType &n)
2237  {
2238  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2239  if (n[v] < 0.)
2240  return false;
2241 
2242  return true;
2243  }
2244  } // namespace Implementation
2245 
2246 
2247 
2248  template <int dim,
2249  int fe_degree,
2250  int n_q_points_1d,
2251  int n_components,
2252  typename VectorType,
2253  typename VectorizedArrayType>
2254  void
2255  LaplaceOperator<dim,
2256  fe_degree,
2257  n_q_points_1d,
2258  n_components,
2259  VectorType,
2260  VectorizedArrayType>::
2261  do_operation_on_cell(
2262  FEEvaluation<
2263  dim,
2264  fe_degree,
2265  n_q_points_1d,
2266  n_components,
2268  const unsigned int cell) const
2269  {
2270  phi.evaluate(EvaluationFlags::gradients);
2271  if (scalar_coefficient.get())
2272  {
2273  Assert(scalar_coefficient->size(1) == 1 ||
2274  scalar_coefficient->size(1) == phi.n_q_points,
2275  ExcMessage("The number of columns in the coefficient table must "
2276  "be either 1 or the number of quadrature points " +
2277  std::to_string(phi.n_q_points) +
2278  ", but the given value was " +
2279  std::to_string(scalar_coefficient->size(1))));
2280  if (scalar_coefficient->size(1) == phi.n_q_points)
2281  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2282  {
2284  (*scalar_coefficient)(cell, q)),
2285  ExcMessage("Coefficient must be non-negative"));
2286  phi.submit_gradient((*scalar_coefficient)(cell, q) *
2287  phi.get_gradient(q),
2288  q);
2289  }
2290  else
2291  {
2292  Assert(Implementation::non_negative((*scalar_coefficient)(cell, 0)),
2293  ExcMessage("Coefficient must be non-negative"));
2294  const VectorizedArrayType coefficient =
2295  (*scalar_coefficient)(cell, 0);
2296  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2297  phi.submit_gradient(coefficient * phi.get_gradient(q), q);
2298  }
2299  }
2300  else
2301  {
2302  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2303  {
2304  phi.submit_gradient(phi.get_gradient(q), q);
2305  }
2306  }
2307  phi.integrate(EvaluationFlags::gradients);
2308  }
2309 
2310 
2311 
2312  template <int dim,
2313  int fe_degree,
2314  int n_q_points_1d,
2315  int n_components,
2316  typename VectorType,
2317  typename VectorizedArrayType>
2318  void
2319  LaplaceOperator<dim,
2320  fe_degree,
2321  n_q_points_1d,
2322  n_components,
2323  VectorType,
2324  VectorizedArrayType>::
2325  local_apply_cell(
2326  const MatrixFree<
2327  dim,
2329  VectorizedArrayType> & data,
2330  VectorType & dst,
2331  const VectorType & src,
2332  const std::pair<unsigned int, unsigned int> &cell_range) const
2333  {
2334  using Number =
2337  data, this->selected_rows[0]);
2338  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2339  {
2340  phi.reinit(cell);
2341  phi.read_dof_values(src);
2342  do_operation_on_cell(phi, cell);
2343  phi.distribute_local_to_global(dst);
2344  }
2345  }
2346 
2347 
2348  template <int dim,
2349  int fe_degree,
2350  int n_q_points_1d,
2351  int n_components,
2352  typename VectorType,
2353  typename VectorizedArrayType>
2354  void
2355  LaplaceOperator<dim,
2356  fe_degree,
2357  n_q_points_1d,
2358  n_components,
2359  VectorType,
2360  VectorizedArrayType>::
2361  local_diagonal_cell(
2362  const MatrixFree<
2363  dim,
2365  VectorizedArrayType> &data,
2366  VectorType & dst,
2367  const VectorType &,
2368  const std::pair<unsigned int, unsigned int> &cell_range) const
2369  {
2370  using Number =
2372 
2374  data, this->selected_rows[0]);
2375  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2376  {
2377  phi.reinit(cell);
2378  VectorizedArrayType local_diagonal_vector[phi.static_dofs_per_cell];
2379  for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2380  {
2381  for (unsigned int j = 0; j < phi.dofs_per_component; ++j)
2382  phi.begin_dof_values()[j] = VectorizedArrayType();
2383  phi.begin_dof_values()[i] = 1.;
2384  do_operation_on_cell(phi, cell);
2385  local_diagonal_vector[i] = phi.begin_dof_values()[i];
2386  }
2387  for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2388  for (unsigned int c = 0; c < phi.n_components; ++c)
2389  phi.begin_dof_values()[i + c * phi.dofs_per_component] =
2390  local_diagonal_vector[i];
2391  phi.distribute_local_to_global(dst);
2392  }
2393  }
2394 
2395 
2396 } // end of namespace MatrixFreeOperators
2397 
2398 
2400 
2401 #endif
size_type size() 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
void read_dof_values(const VectorType &src, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip())
const Number * begin_dof_values() const
const unsigned int dofs_per_component
void reinit(const unsigned int cell_batch_index)
static constexpr unsigned int static_dofs_per_cell
static constexpr unsigned int n_components
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1850
bool is_element(const size_type index) const
Definition: index_set.h:1743
virtual ~Base() override=default
std::vector< unsigned int > selected_rows
Definition: operators.h:456
void Tvmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1420
void set_constrained_entries_to_one(VectorType &dst) const
Definition: operators.h:1377
virtual void compute_diagonal()=0
void vmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1409
void vmult_interface_down(VectorType &dst, const VectorType &src) const
Definition: operators.h:1557
size_type n() const
Definition: operators.h:1174
void preprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1473
void mult_add(VectorType &dst, const VectorType &src, const bool transpose) const
Definition: operators.h:1503
std::vector< std::vector< unsigned int > > edge_constrained_indices
Definition: operators.h:468
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal() const
Definition: operators.h:1702
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1652
std::vector< std::vector< std::pair< value_type, value_type > > > edge_constrained_values
Definition: operators.h:474
size_type m() const
Definition: operators.h:1160
std::vector< unsigned int > selected_columns
Definition: operators.h:462
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > get_matrix_free() const
Definition: operators.h:1680
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal_inverse() const
Definition: operators.h:1689
virtual void clear()
Definition: operators.h:1188
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1213
std::shared_ptr< DiagonalMatrix< VectorType > > diagonal_entries
Definition: operators.h:444
virtual void Tapply_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1713
void initialize(std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType >> data, const MGConstrainedDoFs &mg_constrained_dofs, const unsigned int level, const std::vector< unsigned int > &selected_row_blocks=std::vector< unsigned int >())
Definition: operators.h:1291
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > data
Definition: operators.h:438
void adjust_ghost_range_if_necessary(const VectorType &vec, const bool is_row) const
Definition: operators.h:1431
virtual std::size_t memory_consumption() const
Definition: operators.h:1666
void initialize(std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType >> data_, const std::vector< MGConstrainedDoFs > &mg_constrained_dofs, const unsigned int level, const std::vector< unsigned int > &selected_row_blocks=std::vector< unsigned int >())
Definition: operators.h:1307
std::shared_ptr< DiagonalMatrix< VectorType > > inverse_diagonal_entries
Definition: operators.h:450
value_type el(const unsigned int row, const unsigned int col) const
Definition: operators.h:1198
virtual void apply_add(VectorType &dst, const VectorType &src) const =0
typename VectorType::size_type size_type
Definition: operators.h:201
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:1238
void precondition_Jacobi(VectorType &dst, const VectorType &src, const value_type omega) const
Definition: operators.h:1724
typename VectorType::value_type value_type
Definition: operators.h:196
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1396
void vmult_interface_up(VectorType &dst, const VectorType &src) const
Definition: operators.h:1611
void postprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1523
const FEEvaluationBase< dim, n_components, Number, false, VectorizedArrayType > & fe_eval
Definition: operators.h:725
void fill_inverse_JxW_values(AlignedVector< VectorizedArrayType > &inverse_jxw) const
Definition: operators.h:1030
CellwiseInverseMassMatrix(const FEEvaluationBase< dim, n_components, Number, false, VectorizedArrayType > &fe_eval)
Definition: operators.h:1008
void transform_from_q_points_to_basis(const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition: operators.h:1125
void apply(const AlignedVector< VectorizedArrayType > &inverse_coefficient, const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition: operators.h:1086
std::shared_ptr< Table< 2, VectorizedArrayType > > get_coefficient()
Definition: operators.h:2155
typename Base< dim, VectorType, VectorizedArrayType >::value_type value_type
Definition: operators.h:858
virtual void compute_diagonal() override
Definition: operators.h:2175
std::shared_ptr< Table< 2, VectorizedArrayType > > scalar_coefficient
Definition: operators.h:991
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:2361
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:2261
void set_coefficient(const std::shared_ptr< Table< 2, VectorizedArrayType >> &scalar_coefficient)
Definition: operators.h:2135
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:2325
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition: operators.h:2225
typename Base< dim, VectorType, VectorizedArrayType >::size_type size_type
Definition: operators.h:864
virtual void clear() override
Definition: operators.h:2114
typename OperatorType::value_type value_type
Definition: operators.h:546
SmartPointer< const OperatorType > mf_base_operator
Definition: operators.h:596
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1769
void initialize(const OperatorType &operator_in)
Definition: operators.h:1759
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1789
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1809
typename OperatorType::size_type size_type
Definition: operators.h:551
std::shared_ptr< DiagonalMatrix< VectorType > > lumped_diagonal_entries
Definition: operators.h:823
typename Base< dim, VectorType, VectorizedArrayType >::size_type size_type
Definition: operators.h:757
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition: operators.h:2030
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_lumped_diagonal() const
Definition: operators.h:2008
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_lumped_diagonal_inverse() const
Definition: operators.h:1986
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:2052
virtual void compute_diagonal() override
Definition: operators.h:1850
typename Base< dim, VectorType, VectorizedArrayType >::value_type value_type
Definition: operators.h:751
std::shared_ptr< DiagonalMatrix< VectorType > > inverse_lumped_diagonal_entries
Definition: operators.h:829
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
unsigned int level
Definition: grid_out.cc:4607
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1473
std::string to_string(const T &t)
Definition: patterns.h:2403
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
types::global_dof_index size_type
Definition: cuda_kernels.h:45
std::enable_if< IsBlockVector< VectorType >::value, unsigned int >::type n_blocks(const VectorType &vector)
Definition: operators.h:50
std::enable_if< IsBlockVector< VectorType >::value, void >::type collect_sizes(VectorType &vector)
Definition: operators.h:98
std::enable_if< IsBlockVector< VectorType >::value, typename VectorType::BlockType & >::type subblock(VectorType &vector, unsigned int block_no)
Definition: operators.h:66
bool non_negative(const VectorizedArrayType &n)
Definition: operators.h:2236
void compute_diagonal(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, VectorType &diagonal_global, const std::function< void(FEEvaluation< dim, fe_degree, n_q_points_1d, n_components, Number, VectorizedArrayType > &)> &local_vmult, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
T max(const T &t, const MPI_Comm &mpi_communicator)
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:462
void run(const Iterator &begin, const typename identity< Iterator >::type &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
Definition: work_stream.h:474
static const unsigned int invalid_unsigned_int
Definition: types.h:201
auto apply(F &&fn, Tuple &&t) -> decltype(apply_impl(std::forward< F >(fn), std::forward< Tuple >(t), std::make_index_sequence< std::tuple_size< typename std::remove_reference< Tuple >::type >::value >()))
Definition: tuple.h:36
static void apply(const unsigned int n_components, const FEEvaluationData< dim, Number, false > &fe_eval, const Number *in_array, Number *out_array)
static void transform_from_q_points_to_basis(const unsigned int n_components, const FEEvaluationData< dim, Number, false > &fe_eval, const Number *in_array, Number *out_array)