Reference documentation for deal.II version Git 0943bc0020 2021-10-22 11:23:14 -0400
\(\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 
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
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
99  {
100  vector.collect_sizes();
101  }
102 
103  template <typename VectorType>
104  typename std::enable_if<!IsBlockVector<VectorType>::value, void>::type
106  {}
107  } // namespace BlockHelper
108 
186  template <int dim,
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
347  memory_consumption() const;
348 
352  void
353  initialize_dof_vector(VectorType &vec) const;
354 
367  virtual void
368  compute_diagonal() = 0;
369 
373  std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
374  get_matrix_free() const;
375 
379  const std::shared_ptr<DiagonalMatrix<VectorType>> &
380  get_matrix_diagonal_inverse() const;
381 
385  const std::shared_ptr<DiagonalMatrix<VectorType>> &
386  get_matrix_diagonal() const;
387 
388 
394  void
395  precondition_Jacobi(VectorType & dst,
396  const VectorType &src,
397  const value_type omega) const;
398 
399  protected:
404  void
405  preprocess_constraints(VectorType &dst, const VectorType &src) const;
406 
411  void
412  postprocess_constraints(VectorType &dst, const VectorType &src) const;
413 
418  void
419  set_constrained_entries_to_one(VectorType &dst) const;
420 
424  virtual void
425  apply_add(VectorType &dst, const VectorType &src) const = 0;
426 
432  virtual void
433  Tapply_add(VectorType &dst, const VectorType &src) const;
434 
438  std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
440 
445  std::shared_ptr<DiagonalMatrix<VectorType>> diagonal_entries;
446 
451  std::shared_ptr<DiagonalMatrix<VectorType>> inverse_diagonal_entries;
452 
457  std::vector<unsigned int> selected_rows;
458 
463  std::vector<unsigned int> selected_columns;
464 
465  private:
469  std::vector<std::vector<unsigned int>> edge_constrained_indices;
470 
474  mutable std::vector<std::vector<std::pair<value_type, value_type>>>
476 
482 
487  void
488  mult_add(VectorType & dst,
489  const VectorType &src,
490  const bool transpose) const;
491 
499  void
500  adjust_ghost_range_if_necessary(const VectorType &vec,
501  const bool is_row) const;
502  };
503 
504 
505 
540  template <typename OperatorType>
542  {
543  public:
547  using value_type = typename OperatorType::value_type;
548 
553 
558 
562  void
563  clear();
564 
568  void
569  initialize(const OperatorType &operator_in);
570 
574  template <typename VectorType>
575  void
576  vmult(VectorType &dst, const VectorType &src) const;
577 
581  template <typename VectorType>
582  void
583  Tvmult(VectorType &dst, const VectorType &src) const;
584 
588  template <typename VectorType>
589  void
590  initialize_dof_vector(VectorType &vec) const;
591 
592 
593  private:
598  };
599 
600 
601 
619  template <int dim,
620  int fe_degree,
621  int n_components = 1,
622  typename Number = double,
623  typename VectorizedArrayType = VectorizedArray<Number>>
625  {
626  static_assert(
627  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
628  "Type of Number and of VectorizedArrayType do not match.");
629 
630  public:
636  const FEEvaluationBase<dim,
637  n_components,
638  Number,
639  false,
640  VectorizedArrayType> &fe_eval);
641 
650  void
651  apply(const AlignedVector<VectorizedArrayType> &inverse_coefficient,
652  const unsigned int n_actual_components,
653  const VectorizedArrayType * in_array,
654  VectorizedArrayType * out_array) const;
655 
667  void
668  apply(const VectorizedArrayType *in_array,
669  VectorizedArrayType * out_array) const;
670 
704  void
705  transform_from_q_points_to_basis(const unsigned int n_actual_components,
706  const VectorizedArrayType *in_array,
707  VectorizedArrayType *out_array) const;
708 
714  void
715  fill_inverse_JxW_values(
716  AlignedVector<VectorizedArrayType> &inverse_jxw) const;
717 
718  private:
722  const FEEvaluationBase<dim,
723  n_components,
724  Number,
725  false,
726  VectorizedArrayType> &fe_eval;
727  };
728 
729 
730 
738  template <int dim,
739  int fe_degree,
740  int n_q_points_1d = fe_degree + 1,
741  int n_components = 1,
742  typename VectorType = LinearAlgebra::distributed::Vector<double>,
743  typename VectorizedArrayType =
745  class MassOperator : public Base<dim, VectorType, VectorizedArrayType>
746  {
747  public:
751  using value_type =
753 
757  using size_type =
759 
763  MassOperator();
764 
768  virtual void
769  compute_diagonal() override;
770 
786  void
787  compute_lumped_diagonal();
788 
792  const std::shared_ptr<DiagonalMatrix<VectorType>> &
793  get_matrix_lumped_diagonal() const;
794 
798  const std::shared_ptr<DiagonalMatrix<VectorType>> &
799  get_matrix_lumped_diagonal_inverse() const;
800 
801  private:
807  virtual void
808  apply_add(VectorType &dst, const VectorType &src) const override;
809 
813  void
814  local_apply_cell(
816  VectorType & dst,
817  const VectorType & src,
818  const std::pair<unsigned int, unsigned int> & cell_range) const;
819 
824  std::shared_ptr<DiagonalMatrix<VectorType>> lumped_diagonal_entries;
825 
830  std::shared_ptr<DiagonalMatrix<VectorType>> inverse_lumped_diagonal_entries;
831  };
832 
833 
834 
845  template <int dim,
846  int fe_degree,
847  int n_q_points_1d = fe_degree + 1,
848  int n_components = 1,
849  typename VectorType = LinearAlgebra::distributed::Vector<double>,
850  typename VectorizedArrayType =
852  class LaplaceOperator : public Base<dim, VectorType, VectorizedArrayType>
853  {
854  public:
858  using value_type =
860 
864  using size_type =
866 
870  LaplaceOperator();
871 
878  virtual void
879  compute_diagonal() override;
880 
931  void
932  set_coefficient(
933  const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient);
934 
939  virtual void
940  clear() override;
941 
948  std::shared_ptr<Table<2, VectorizedArrayType>>
949  get_coefficient();
950 
951  private:
957  virtual void
958  apply_add(VectorType &dst, const VectorType &src) const override;
959 
963  void
964  local_apply_cell(
966  VectorType & dst,
967  const VectorType & src,
968  const std::pair<unsigned int, unsigned int> & cell_range) const;
969 
973  void
974  local_diagonal_cell(
976  VectorType & dst,
977  const VectorType &,
978  const std::pair<unsigned int, unsigned int> &cell_range) const;
979 
983  void
984  do_operation_on_cell(
986  & phi,
987  const unsigned int cell) const;
988 
992  std::shared_ptr<Table<2, VectorizedArrayType>> scalar_coefficient;
993  };
994 
995 
996 
997  // ------------------------------------ inline functions ---------------------
998 
999  template <int dim,
1000  int fe_degree,
1001  int n_components,
1002  typename Number,
1003  typename VectorizedArrayType>
1004  inline CellwiseInverseMassMatrix<dim,
1005  fe_degree,
1006  n_components,
1007  Number,
1008  VectorizedArrayType>::
1009  CellwiseInverseMassMatrix(
1010  const FEEvaluationBase<dim,
1011  n_components,
1012  Number,
1013  false,
1014  VectorizedArrayType> &fe_eval)
1015  : fe_eval(fe_eval)
1016  {}
1017 
1018 
1019 
1020  template <int dim,
1021  int fe_degree,
1022  int n_components,
1023  typename Number,
1024  typename VectorizedArrayType>
1025  inline void
1027  fe_degree,
1028  n_components,
1029  Number,
1030  VectorizedArrayType>::
1032  AlignedVector<VectorizedArrayType> &inverse_jxw) const
1033  {
1034  constexpr unsigned int dofs_per_component_on_cell =
1035  Utilities::pow(fe_degree + 1, dim);
1036  Assert(inverse_jxw.size() > 0 &&
1037  inverse_jxw.size() % dofs_per_component_on_cell == 0,
1038  ExcMessage(
1039  "Expected diagonal to be a multiple of scalar dof per cells"));
1040 
1041  // compute values for the first component
1042  for (unsigned int q = 0; q < dofs_per_component_on_cell; ++q)
1043  inverse_jxw[q] = 1. / fe_eval.JxW(q);
1044  // copy values to rest of vector
1045  for (unsigned int q = dofs_per_component_on_cell; q < inverse_jxw.size();)
1046  for (unsigned int i = 0; i < dofs_per_component_on_cell; ++i, ++q)
1047  inverse_jxw[q] = inverse_jxw[i];
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  dim,
1060  fe_degree,
1061  n_components,
1062  Number,
1063  VectorizedArrayType>::apply(const VectorizedArrayType *in_array,
1064  VectorizedArrayType * out_array) const
1065  {
1067  template run<fe_degree>(n_components, fe_eval, in_array, out_array);
1068  }
1069 
1070 
1071 
1072  template <int dim,
1073  int fe_degree,
1074  int n_components,
1075  typename Number,
1076  typename VectorizedArrayType>
1077  inline void
1079  fe_degree,
1080  n_components,
1081  Number,
1082  VectorizedArrayType>::
1083  apply(const AlignedVector<VectorizedArrayType> &inverse_coefficients,
1084  const unsigned int n_actual_components,
1085  const VectorizedArrayType * in_array,
1086  VectorizedArrayType * out_array) const
1087  {
1089  template run<fe_degree>(
1090  n_actual_components,
1091  fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
1092  inverse_coefficients,
1093  in_array,
1094  out_array);
1095  }
1096 
1097 
1098 
1099  template <int dim,
1100  int fe_degree,
1101  int n_components,
1102  typename Number,
1103  typename VectorizedArrayType>
1104  inline void
1106  fe_degree,
1107  n_components,
1108  Number,
1109  VectorizedArrayType>::
1110  transform_from_q_points_to_basis(const unsigned int n_actual_components,
1111  const VectorizedArrayType *in_array,
1112  VectorizedArrayType * out_array) const
1113  {
1114  const auto n_q_points_1d = fe_eval.get_shape_info().data[0].n_q_points_1d;
1115 
1116  if (fe_degree > -1 && (fe_degree + 1 == n_q_points_1d))
1118  dim,
1119  VectorizedArrayType>::template run<fe_degree,
1120  fe_degree + 1>(n_actual_components,
1121  fe_eval,
1122  in_array,
1123  out_array);
1124  else
1126  transform_from_q_points_to_basis(n_actual_components,
1127  fe_degree,
1128  n_q_points_1d,
1129  fe_eval,
1130  in_array,
1131  out_array);
1132  }
1133 
1134 
1135 
1136  //----------------- Base operator -----------------------------
1137  template <int dim, typename VectorType, typename VectorizedArrayType>
1139  : Subscriptor()
1140  , have_interface_matrices(false)
1141  {}
1142 
1143 
1144 
1145  template <int dim, typename VectorType, typename VectorizedArrayType>
1148  {
1149  Assert(data.get() != nullptr, ExcNotInitialized());
1151  0;
1152  for (const unsigned int selected_row : selected_rows)
1153  total_size += data->get_vector_partitioner(selected_row)->size();
1154  return total_size;
1155  }
1156 
1157 
1158 
1159  template <int dim, typename VectorType, typename VectorizedArrayType>
1162  {
1163  Assert(data.get() != nullptr, ExcNotInitialized());
1165  0;
1166  for (const unsigned int selected_column : selected_columns)
1167  total_size += data->get_vector_partitioner(selected_column)->size();
1168  return total_size;
1169  }
1170 
1171 
1172 
1173  template <int dim, typename VectorType, typename VectorizedArrayType>
1174  void
1176  {
1177  data.reset();
1178  inverse_diagonal_entries.reset();
1179  }
1180 
1181 
1182 
1183  template <int dim, typename VectorType, typename VectorizedArrayType>
1186  const unsigned int col) const
1187  {
1188  (void)col;
1189  Assert(row == col, ExcNotImplemented());
1190  Assert(inverse_diagonal_entries.get() != nullptr &&
1191  inverse_diagonal_entries->m() > 0,
1192  ExcNotInitialized());
1193  return 1.0 / (*inverse_diagonal_entries)(row, row);
1194  }
1195 
1196 
1197 
1198  template <int dim, typename VectorType, typename VectorizedArrayType>
1199  void
1201  VectorType &vec) const
1202  {
1203  Assert(data.get() != nullptr, ExcNotInitialized());
1205  for (unsigned int i = 0; i < BlockHelper::n_blocks(vec); ++i)
1206  {
1207  const unsigned int index = selected_rows[i];
1208  if (!BlockHelper::subblock(vec, index)
1209  .partitioners_are_compatible(
1210  *data->get_dof_info(index).vector_partitioner))
1211  data->initialize_dof_vector(BlockHelper::subblock(vec, index), index);
1212 
1213  Assert(BlockHelper::subblock(vec, index)
1214  .partitioners_are_globally_compatible(
1215  *data->get_dof_info(index).vector_partitioner),
1216  ExcInternalError());
1217  }
1219  }
1220 
1221 
1222 
1223  template <int dim, typename VectorType, typename VectorizedArrayType>
1224  void
1227  data_,
1228  const std::vector<unsigned int> &given_row_selection,
1229  const std::vector<unsigned int> &given_column_selection)
1230  {
1231  data = data_;
1232 
1233  selected_rows.clear();
1234  selected_columns.clear();
1235  if (given_row_selection.empty())
1236  for (unsigned int i = 0; i < data_->n_components(); ++i)
1237  selected_rows.push_back(i);
1238  else
1239  {
1240  for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1241  {
1242  AssertIndexRange(given_row_selection[i], data_->n_components());
1243  for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1244  if (j != i)
1245  Assert(given_row_selection[j] != given_row_selection[i],
1246  ExcMessage("Given row indices must be unique"));
1247 
1248  selected_rows.push_back(given_row_selection[i]);
1249  }
1250  }
1251  if (given_column_selection.size() == 0)
1253  else
1254  {
1255  for (unsigned int i = 0; i < given_column_selection.size(); ++i)
1256  {
1257  AssertIndexRange(given_column_selection[i], data_->n_components());
1258  for (unsigned int j = 0; j < given_column_selection.size(); ++j)
1259  if (j != i)
1260  Assert(given_column_selection[j] != given_column_selection[i],
1261  ExcMessage("Given column indices must be unique"));
1262 
1263  selected_columns.push_back(given_column_selection[i]);
1264  }
1265  }
1266 
1267  edge_constrained_indices.clear();
1268  edge_constrained_indices.resize(selected_rows.size());
1269  edge_constrained_values.clear();
1270  edge_constrained_values.resize(selected_rows.size());
1271  have_interface_matrices = false;
1272  }
1273 
1274 
1275 
1276  template <int dim, typename VectorType, typename VectorizedArrayType>
1277  void
1280  data_,
1281  const MGConstrainedDoFs & mg_constrained_dofs,
1282  const unsigned int level,
1283  const std::vector<unsigned int> &given_row_selection)
1284  {
1285  std::vector<MGConstrainedDoFs> mg_constrained_dofs_vector(
1286  1, mg_constrained_dofs);
1287  initialize(data_, mg_constrained_dofs_vector, level, given_row_selection);
1288  }
1289 
1290 
1291 
1292  template <int dim, typename VectorType, typename VectorizedArrayType>
1293  void
1296  data_,
1297  const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
1298  const unsigned int level,
1299  const std::vector<unsigned int> & given_row_selection)
1300  {
1302  ExcMessage("level is not set"));
1303 
1304  selected_rows.clear();
1305  selected_columns.clear();
1306  if (given_row_selection.empty())
1307  for (unsigned int i = 0; i < data_->n_components(); ++i)
1308  selected_rows.push_back(i);
1309  else
1310  {
1311  for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1312  {
1313  AssertIndexRange(given_row_selection[i], data_->n_components());
1314  for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1315  if (j != i)
1316  Assert(given_row_selection[j] != given_row_selection[i],
1317  ExcMessage("Given row indices must be unique"));
1318 
1319  selected_rows.push_back(given_row_selection[i]);
1320  }
1321  }
1323 
1324  AssertDimension(mg_constrained_dofs.size(), selected_rows.size());
1325  edge_constrained_indices.clear();
1326  edge_constrained_indices.resize(selected_rows.size());
1327  edge_constrained_values.clear();
1328  edge_constrained_values.resize(selected_rows.size());
1329 
1330  data = data_;
1331 
1332  for (unsigned int j = 0; j < selected_rows.size(); ++j)
1333  {
1334  if (data_->n_cell_batches() > 0)
1335  {
1336  AssertDimension(level, data_->get_cell_iterator(0, 0, j)->level());
1337  }
1338 
1339  // setup edge_constrained indices
1340  std::vector<types::global_dof_index> interface_indices;
1341  mg_constrained_dofs[j]
1342  .get_refinement_edge_indices(level)
1343  .fill_index_vector(interface_indices);
1344  edge_constrained_indices[j].clear();
1345  edge_constrained_indices[j].reserve(interface_indices.size());
1346  edge_constrained_values[j].resize(interface_indices.size());
1347  const IndexSet &locally_owned =
1348  data->get_dof_handler(selected_rows[j]).locally_owned_mg_dofs(level);
1349  for (const auto interface_index : interface_indices)
1350  if (locally_owned.is_element(interface_index))
1351  edge_constrained_indices[j].push_back(
1352  locally_owned.index_within_set(interface_index));
1355  static_cast<unsigned int>(edge_constrained_indices[j].size()),
1356  data->get_vector_partitioner()->get_mpi_communicator()) > 0;
1357  }
1358  }
1359 
1360 
1361 
1362  template <int dim, typename VectorType, typename VectorizedArrayType>
1363  void
1365  VectorType &dst) const
1366  {
1367  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1368  {
1369  const std::vector<unsigned int> &constrained_dofs =
1370  data->get_constrained_dofs(selected_rows[j]);
1371  for (const auto constrained_dof : constrained_dofs)
1372  BlockHelper::subblock(dst, j).local_element(constrained_dof) = 1.;
1373  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1374  BlockHelper::subblock(dst, j).local_element(
1375  edge_constrained_indices[j][i]) = 1.;
1376  }
1377  }
1378 
1379 
1380 
1381  template <int dim, typename VectorType, typename VectorizedArrayType>
1382  void
1384  const VectorType &src) const
1385  {
1386  using Number =
1388  dst = Number(0.);
1389  vmult_add(dst, src);
1390  }
1391 
1392 
1393 
1394  template <int dim, typename VectorType, typename VectorizedArrayType>
1395  void
1397  VectorType & dst,
1398  const VectorType &src) const
1399  {
1400  mult_add(dst, src, false);
1401  }
1402 
1403 
1404 
1405  template <int dim, typename VectorType, typename VectorizedArrayType>
1406  void
1408  VectorType & dst,
1409  const VectorType &src) const
1410  {
1411  mult_add(dst, src, true);
1412  }
1413 
1414 
1415 
1416  template <int dim, typename VectorType, typename VectorizedArrayType>
1417  void
1419  const VectorType &src,
1420  const bool is_row) const
1421  {
1422  using Number =
1424  for (unsigned int i = 0; i < BlockHelper::n_blocks(src); ++i)
1425  {
1426  const unsigned int mf_component =
1427  is_row ? selected_rows[i] : selected_columns[i];
1428  // If both vectors use the same partitioner -> done
1429  if (BlockHelper::subblock(src, i).get_partitioner().get() ==
1430  data->get_dof_info(mf_component).vector_partitioner.get())
1431  continue;
1432 
1433  // If not, assert that the local ranges are the same and reset to the
1434  // current partitioner
1436  .get_partitioner()
1437  ->locally_owned_size() ==
1438  data->get_dof_info(mf_component)
1439  .vector_partitioner->locally_owned_size(),
1440  ExcMessage(
1441  "The vector passed to the vmult() function does not have "
1442  "the correct size for compatibility with MatrixFree."));
1443 
1444  // copy the vector content to a temporary vector so that it does not get
1445  // lost
1447  BlockHelper::subblock(src, i));
1448  this->data->initialize_dof_vector(
1449  BlockHelper::subblock(const_cast<VectorType &>(src), i),
1450  mf_component);
1451  BlockHelper::subblock(const_cast<VectorType &>(src), i)
1452  .copy_locally_owned_data_from(copy_vec);
1453  }
1454  }
1455 
1456 
1457 
1458  template <int dim, typename VectorType, typename VectorizedArrayType>
1459  void
1461  VectorType & dst,
1462  const VectorType &src) const
1463  {
1464  using Number =
1466  adjust_ghost_range_if_necessary(src, false);
1468 
1469  // set zero Dirichlet values on the input vector (and remember the src and
1470  // dst values because we need to reset them at the end)
1471  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1472  {
1473  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1474  {
1475  edge_constrained_values[j][i] = std::pair<Number, Number>(
1476  BlockHelper::subblock(src, j).local_element(
1477  edge_constrained_indices[j][i]),
1478  BlockHelper::subblock(dst, j).local_element(
1479  edge_constrained_indices[j][i]));
1480  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1481  .local_element(edge_constrained_indices[j][i]) = 0.;
1482  }
1483  }
1484  }
1485 
1486 
1487 
1488  template <int dim, typename VectorType, typename VectorizedArrayType>
1489  void
1491  VectorType & dst,
1492  const VectorType &src,
1493  const bool transpose) const
1494  {
1495  AssertDimension(dst.size(), src.size());
1498  preprocess_constraints(dst, src);
1499  if (transpose)
1500  Tapply_add(dst, src);
1501  else
1502  apply_add(dst, src);
1503  postprocess_constraints(dst, src);
1504  }
1505 
1506 
1507 
1508  template <int dim, typename VectorType, typename VectorizedArrayType>
1509  void
1511  VectorType & dst,
1512  const VectorType &src) const
1513  {
1514  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1515  {
1516  const std::vector<unsigned int> &constrained_dofs =
1517  data->get_constrained_dofs(selected_rows[j]);
1518  for (const auto constrained_dof : constrained_dofs)
1519  BlockHelper::subblock(dst, j).local_element(constrained_dof) +=
1520  BlockHelper::subblock(src, j).local_element(constrained_dof);
1521  }
1522 
1523  // reset edge constrained values, multiply by unit matrix and add into
1524  // destination
1525  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1526  {
1527  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1528  {
1529  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1530  .local_element(edge_constrained_indices[j][i]) =
1531  edge_constrained_values[j][i].first;
1532  BlockHelper::subblock(dst, j).local_element(
1533  edge_constrained_indices[j][i]) =
1534  edge_constrained_values[j][i].second +
1535  edge_constrained_values[j][i].first;
1536  }
1537  }
1538  }
1539 
1540 
1541 
1542  template <int dim, typename VectorType, typename VectorizedArrayType>
1543  void
1545  VectorType & dst,
1546  const VectorType &src) const
1547  {
1548  using Number =
1550  AssertDimension(dst.size(), src.size());
1551  adjust_ghost_range_if_necessary(src, false);
1553 
1554  dst = Number(0.);
1555 
1557  return;
1558 
1559  // set zero Dirichlet values on the input vector (and remember the src and
1560  // dst values because we need to reset them at the end)
1561  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1562  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1563  {
1564  edge_constrained_values[j][i] = std::pair<Number, Number>(
1565  BlockHelper::subblock(src, j).local_element(
1566  edge_constrained_indices[j][i]),
1567  BlockHelper::subblock(dst, j).local_element(
1568  edge_constrained_indices[j][i]));
1569  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1570  .local_element(edge_constrained_indices[j][i]) = 0.;
1571  }
1572 
1573  apply_add(dst, src);
1574 
1575  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1576  {
1577  unsigned int c = 0;
1578  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1579  {
1580  for (; c < edge_constrained_indices[j][i]; ++c)
1581  BlockHelper::subblock(dst, j).local_element(c) = 0.;
1582  ++c;
1583 
1584  // reset the src values
1585  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1586  .local_element(edge_constrained_indices[j][i]) =
1587  edge_constrained_values[j][i].first;
1588  }
1589  for (; c < BlockHelper::subblock(dst, j).locally_owned_size(); ++c)
1590  BlockHelper::subblock(dst, j).local_element(c) = 0.;
1591  }
1592  }
1593 
1594 
1595 
1596  template <int dim, typename VectorType, typename VectorizedArrayType>
1597  void
1599  VectorType & dst,
1600  const VectorType &src) const
1601  {
1602  using Number =
1604  AssertDimension(dst.size(), src.size());
1605  adjust_ghost_range_if_necessary(src, false);
1607 
1608  dst = Number(0.);
1609 
1611  return;
1612 
1613  VectorType src_cpy(src);
1614  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1615  {
1616  unsigned int c = 0;
1617  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1618  {
1619  for (; c < edge_constrained_indices[j][i]; ++c)
1620  BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1621  ++c;
1622  }
1623  for (; c < BlockHelper::subblock(src_cpy, j).locally_owned_size(); ++c)
1624  BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1625  }
1626 
1627  apply_add(dst, src_cpy);
1628 
1629  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1630  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1631  BlockHelper::subblock(dst, j).local_element(
1632  edge_constrained_indices[j][i]) = 0.;
1633  }
1634 
1635 
1636 
1637  template <int dim, typename VectorType, typename VectorizedArrayType>
1638  void
1640  VectorType & dst,
1641  const VectorType &src) const
1642  {
1643  using Number =
1645  dst = Number(0.);
1646  Tvmult_add(dst, src);
1647  }
1648 
1649 
1650 
1651  template <int dim, typename VectorType, typename VectorizedArrayType>
1652  std::size_t
1654  {
1655  return inverse_diagonal_entries.get() != nullptr ?
1656  inverse_diagonal_entries->memory_consumption() :
1657  sizeof(*this);
1658  }
1659 
1660 
1661 
1662  template <int dim, typename VectorType, typename VectorizedArrayType>
1663  std::shared_ptr<const MatrixFree<
1664  dim,
1666  VectorizedArrayType>>
1668  {
1669  return data;
1670  }
1671 
1672 
1673 
1674  template <int dim, typename VectorType, typename VectorizedArrayType>
1675  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1677  const
1678  {
1679  Assert(inverse_diagonal_entries.get() != nullptr &&
1680  inverse_diagonal_entries->m() > 0,
1681  ExcNotInitialized());
1682  return inverse_diagonal_entries;
1683  }
1684 
1685 
1686 
1687  template <int dim, typename VectorType, typename VectorizedArrayType>
1688  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1690  {
1691  Assert(diagonal_entries.get() != nullptr && diagonal_entries->m() > 0,
1692  ExcNotInitialized());
1693  return diagonal_entries;
1694  }
1695 
1696 
1697 
1698  template <int dim, typename VectorType, typename VectorizedArrayType>
1699  void
1701  VectorType & dst,
1702  const VectorType &src) const
1703  {
1704  apply_add(dst, src);
1705  }
1706 
1707 
1708 
1709  template <int dim, typename VectorType, typename VectorizedArrayType>
1710  void
1712  VectorType & dst,
1713  const VectorType & src,
1714  const typename Base<dim, VectorType, VectorizedArrayType>::value_type omega)
1715  const
1716  {
1718  ExcNotInitialized());
1719  inverse_diagonal_entries->vmult(dst, src);
1720  dst *= omega;
1721  }
1722 
1723 
1724 
1725  //------------------------- MGInterfaceOperator ------------------------------
1726 
1727  template <typename OperatorType>
1729  : Subscriptor()
1730  , mf_base_operator(nullptr)
1731  {}
1732 
1733 
1734 
1735  template <typename OperatorType>
1736  void
1738  {
1739  mf_base_operator = nullptr;
1740  }
1741 
1742 
1743 
1744  template <typename OperatorType>
1745  void
1746  MGInterfaceOperator<OperatorType>::initialize(const OperatorType &operator_in)
1747  {
1748  mf_base_operator = &operator_in;
1749  }
1750 
1751 
1752 
1753  template <typename OperatorType>
1754  template <typename VectorType>
1755  void
1757  const VectorType &src) const
1758  {
1759 #ifndef DEAL_II_MSVC
1760  static_assert(
1761  std::is_same<typename VectorType::value_type, value_type>::value,
1762  "The vector type must be based on the same value type as this "
1763  "operator");
1764 #endif
1765 
1766  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1767 
1768  mf_base_operator->vmult_interface_down(dst, src);
1769  }
1770 
1771 
1772 
1773  template <typename OperatorType>
1774  template <typename VectorType>
1775  void
1777  const VectorType &src) const
1778  {
1779 #ifndef DEAL_II_MSVC
1780  static_assert(
1781  std::is_same<typename VectorType::value_type, value_type>::value,
1782  "The vector type must be based on the same value type as this "
1783  "operator");
1784 #endif
1785 
1786  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1787 
1788  mf_base_operator->vmult_interface_up(dst, src);
1789  }
1790 
1791 
1792 
1793  template <typename OperatorType>
1794  template <typename VectorType>
1795  void
1797  VectorType &vec) const
1798  {
1799  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1800 
1801  mf_base_operator->initialize_dof_vector(vec);
1802  }
1803 
1804 
1805 
1806  //-----------------------------MassOperator----------------------------------
1807 
1808  template <int dim,
1809  int fe_degree,
1810  int n_q_points_1d,
1811  int n_components,
1812  typename VectorType,
1813  typename VectorizedArrayType>
1814  MassOperator<dim,
1815  fe_degree,
1816  n_q_points_1d,
1817  n_components,
1818  VectorType,
1819  VectorizedArrayType>::MassOperator()
1820  : Base<dim, VectorType, VectorizedArrayType>()
1821  {}
1822 
1823 
1824 
1825  template <int dim,
1826  int fe_degree,
1827  int n_q_points_1d,
1828  int n_components,
1829  typename VectorType,
1830  typename VectorizedArrayType>
1831  void
1832  MassOperator<dim,
1833  fe_degree,
1834  n_q_points_1d,
1835  n_components,
1836  VectorType,
1837  VectorizedArrayType>::compute_diagonal()
1838  {
1839  using Number =
1840  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1842  ExcNotInitialized());
1843  Assert(this->selected_rows == this->selected_columns,
1844  ExcMessage("This function is only implemented for square (not "
1845  "rectangular) operators."));
1846 
1847  this->inverse_diagonal_entries =
1848  std::make_shared<DiagonalMatrix<VectorType>>();
1849  this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1850  VectorType &inverse_diagonal_vector =
1851  this->inverse_diagonal_entries->get_vector();
1852  VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1853  this->initialize_dof_vector(inverse_diagonal_vector);
1854  this->initialize_dof_vector(diagonal_vector);
1855 
1856  // Set up the action of the mass matrix in a way that's compatible with
1857  // MatrixFreeTools::compute_diagonal:
1858  auto diagonal_evaluation = [](auto &integrator) {
1859  integrator.evaluate(EvaluationFlags::values);
1860  for (unsigned int q = 0; q < integrator.n_q_points; ++q)
1861  integrator.submit_value(integrator.get_value(q), q);
1862  integrator.integrate(EvaluationFlags::values);
1863  };
1864 
1865  std::function<void(
1866  FEEvaluation<
1867  dim,
1868  fe_degree,
1869  n_q_points_1d,
1870  n_components,
1871  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
1872  VectorizedArrayType> &)>
1873  diagonal_evaluation_f(diagonal_evaluation);
1874 
1875  Assert(this->selected_rows.size() > 0, ExcInternalError());
1876  for (unsigned int block_n = 0; block_n < this->selected_rows.size();
1877  ++block_n)
1879  BlockHelper::subblock(diagonal_vector,
1880  block_n),
1881  diagonal_evaluation_f,
1882  this->selected_rows[block_n]);
1883 
1884  // Constrained entries will create zeros on the main diagonal, which we
1885  // don't want
1886  this->set_constrained_entries_to_one(diagonal_vector);
1887 
1888  inverse_diagonal_vector = diagonal_vector;
1889 
1890  for (unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
1891  ++i)
1892  {
1893  Assert(diagonal_vector.local_element(i) > Number(0),
1894  ExcInternalError());
1895  inverse_diagonal_vector.local_element(i) =
1896  1. / inverse_diagonal_vector.local_element(i);
1897  }
1898 
1899  inverse_diagonal_vector.update_ghost_values();
1900  diagonal_vector.update_ghost_values();
1901  }
1902 
1903 
1904 
1905  template <int dim,
1906  int fe_degree,
1907  int n_q_points_1d,
1908  int n_components,
1909  typename VectorType,
1910  typename VectorizedArrayType>
1911  void
1912  MassOperator<dim,
1913  fe_degree,
1914  n_q_points_1d,
1915  n_components,
1916  VectorType,
1917  VectorizedArrayType>::compute_lumped_diagonal()
1918  {
1919  using Number =
1920  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1922  ExcNotInitialized());
1923  Assert(this->selected_rows == this->selected_columns,
1924  ExcMessage("This function is only implemented for square (not "
1925  "rectangular) operators."));
1926 
1928  std::make_shared<DiagonalMatrix<VectorType>>();
1929  lumped_diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1930  VectorType &inverse_lumped_diagonal_vector =
1931  inverse_lumped_diagonal_entries->get_vector();
1932  VectorType &lumped_diagonal_vector = lumped_diagonal_entries->get_vector();
1933  this->initialize_dof_vector(inverse_lumped_diagonal_vector);
1934  this->initialize_dof_vector(lumped_diagonal_vector);
1935 
1936  // Re-use the inverse_lumped_diagonal_vector as the vector of 1s
1937  inverse_lumped_diagonal_vector = Number(1.);
1938  apply_add(lumped_diagonal_vector, inverse_lumped_diagonal_vector);
1939  this->set_constrained_entries_to_one(lumped_diagonal_vector);
1940 
1941  const size_type locally_owned_size =
1942  inverse_lumped_diagonal_vector.locally_owned_size();
1943  // A caller may request a lumped diagonal matrix when it doesn't make sense
1944  // (e.g., an element with negative-mean basis functions). Avoid division by
1945  // zero so we don't cause a floating point exception but permit negative
1946  // entries here.
1947  for (size_type i = 0; i < locally_owned_size; ++i)
1948  {
1949  if (inverse_lumped_diagonal_vector.local_element(i) == Number(0.))
1950  inverse_lumped_diagonal_vector.local_element(i) = Number(1.);
1951  else
1952  inverse_lumped_diagonal_vector.local_element(i) =
1953  Number(1.) / lumped_diagonal_vector.local_element(i);
1954  }
1955 
1956  inverse_lumped_diagonal_vector.update_ghost_values();
1957  lumped_diagonal_vector.update_ghost_values();
1958  }
1959 
1960 
1961 
1962  template <int dim,
1963  int fe_degree,
1964  int n_q_points_1d,
1965  int n_components,
1966  typename VectorType,
1967  typename VectorizedArrayType>
1968  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1969  MassOperator<dim,
1970  fe_degree,
1971  n_q_points_1d,
1972  n_components,
1973  VectorType,
1974  VectorizedArrayType>::get_matrix_lumped_diagonal_inverse() const
1975  {
1976  Assert(inverse_lumped_diagonal_entries.get() != nullptr &&
1978  ExcNotInitialized());
1980  }
1981 
1982 
1983 
1984  template <int dim,
1985  int fe_degree,
1986  int n_q_points_1d,
1987  int n_components,
1988  typename VectorType,
1989  typename VectorizedArrayType>
1990  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1991  MassOperator<dim,
1992  fe_degree,
1993  n_q_points_1d,
1994  n_components,
1995  VectorType,
1996  VectorizedArrayType>::get_matrix_lumped_diagonal() const
1997  {
1998  Assert(lumped_diagonal_entries.get() != nullptr &&
1999  lumped_diagonal_entries->m() > 0,
2000  ExcNotInitialized());
2001  return lumped_diagonal_entries;
2002  }
2003 
2004 
2005 
2006  template <int dim,
2007  int fe_degree,
2008  int n_q_points_1d,
2009  int n_components,
2010  typename VectorType,
2011  typename VectorizedArrayType>
2012  void
2013  MassOperator<dim,
2014  fe_degree,
2015  n_q_points_1d,
2016  n_components,
2017  VectorType,
2018  VectorizedArrayType>::apply_add(VectorType & dst,
2019  const VectorType &src) const
2020  {
2022  &MassOperator::local_apply_cell, this, dst, src);
2023  }
2024 
2025 
2026 
2027  template <int dim,
2028  int fe_degree,
2029  int n_q_points_1d,
2030  int n_components,
2031  typename VectorType,
2032  typename VectorizedArrayType>
2033  void
2034  MassOperator<dim,
2035  fe_degree,
2036  n_q_points_1d,
2037  n_components,
2038  VectorType,
2039  VectorizedArrayType>::
2041  const MatrixFree<
2042  dim,
2043  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2044  VectorizedArrayType> & data,
2045  VectorType & dst,
2046  const VectorType & src,
2047  const std::pair<unsigned int, unsigned int> &cell_range) const
2048  {
2049  using Number =
2050  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2051  FEEvaluation<dim,
2052  fe_degree,
2053  n_q_points_1d,
2054  n_components,
2055  Number,
2056  VectorizedArrayType>
2057  phi(data, this->selected_rows[0]);
2058  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2059  {
2060  phi.reinit(cell);
2061  phi.read_dof_values(src);
2062  phi.evaluate(EvaluationFlags::values);
2063  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2064  phi.submit_value(phi.get_value(q), q);
2065  phi.integrate(EvaluationFlags::values);
2066  phi.distribute_local_to_global(dst);
2067  }
2068  }
2069 
2070 
2071  //-----------------------------LaplaceOperator----------------------------------
2072 
2073  template <int dim,
2074  int fe_degree,
2075  int n_q_points_1d,
2076  int n_components,
2077  typename VectorType,
2078  typename VectorizedArrayType>
2079  LaplaceOperator<dim,
2080  fe_degree,
2081  n_q_points_1d,
2082  n_components,
2083  VectorType,
2084  VectorizedArrayType>::LaplaceOperator()
2085  : Base<dim, VectorType, VectorizedArrayType>()
2086  {}
2087 
2088 
2089 
2090  template <int dim,
2091  int fe_degree,
2092  int n_q_points_1d,
2093  int n_components,
2094  typename VectorType,
2095  typename VectorizedArrayType>
2096  void
2097  LaplaceOperator<dim,
2098  fe_degree,
2099  n_q_points_1d,
2100  n_components,
2101  VectorType,
2102  VectorizedArrayType>::clear()
2103  {
2105  scalar_coefficient.reset();
2106  }
2107 
2108 
2109 
2110  template <int dim,
2111  int fe_degree,
2112  int n_q_points_1d,
2113  int n_components,
2114  typename VectorType,
2115  typename VectorizedArrayType>
2116  void
2117  LaplaceOperator<dim,
2118  fe_degree,
2119  n_q_points_1d,
2120  n_components,
2121  VectorType,
2122  VectorizedArrayType>::
2124  const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient_)
2125  {
2126  scalar_coefficient = scalar_coefficient_;
2127  }
2128 
2129 
2130 
2131  template <int dim,
2132  int fe_degree,
2133  int n_q_points_1d,
2134  int n_components,
2135  typename VectorType,
2136  typename VectorizedArrayType>
2137  std::shared_ptr<Table<2, VectorizedArrayType>>
2138  LaplaceOperator<dim,
2139  fe_degree,
2140  n_q_points_1d,
2141  n_components,
2142  VectorType,
2143  VectorizedArrayType>::get_coefficient()
2144  {
2146  return scalar_coefficient;
2147  }
2148 
2149 
2150 
2151  template <int dim,
2152  int fe_degree,
2153  int n_q_points_1d,
2154  int n_components,
2155  typename VectorType,
2156  typename VectorizedArrayType>
2157  void
2158  LaplaceOperator<dim,
2159  fe_degree,
2160  n_q_points_1d,
2161  n_components,
2162  VectorType,
2163  VectorizedArrayType>::compute_diagonal()
2164  {
2165  using Number =
2166  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2168  ExcNotInitialized());
2169 
2170  this->inverse_diagonal_entries =
2171  std::make_shared<DiagonalMatrix<VectorType>>();
2172  this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
2173  VectorType &inverse_diagonal_vector =
2174  this->inverse_diagonal_entries->get_vector();
2175  VectorType &diagonal_vector = this->diagonal_entries->get_vector();
2176  this->initialize_dof_vector(inverse_diagonal_vector);
2177  this->initialize_dof_vector(diagonal_vector);
2178 
2179  this->data->cell_loop(&LaplaceOperator::local_diagonal_cell,
2180  this,
2181  diagonal_vector,
2182  /*unused*/ diagonal_vector);
2183  this->set_constrained_entries_to_one(diagonal_vector);
2184 
2185  inverse_diagonal_vector = diagonal_vector;
2186 
2187  for (unsigned int i = 0; i < inverse_diagonal_vector.locally_owned_size();
2188  ++i)
2189  if (std::abs(inverse_diagonal_vector.local_element(i)) >
2191  inverse_diagonal_vector.local_element(i) =
2192  1. / inverse_diagonal_vector.local_element(i);
2193  else
2194  inverse_diagonal_vector.local_element(i) = 1.;
2195 
2196  inverse_diagonal_vector.update_ghost_values();
2197  diagonal_vector.update_ghost_values();
2198  }
2199 
2200 
2201 
2202  template <int dim,
2203  int fe_degree,
2204  int n_q_points_1d,
2205  int n_components,
2206  typename VectorType,
2207  typename VectorizedArrayType>
2208  void
2209  LaplaceOperator<dim,
2210  fe_degree,
2211  n_q_points_1d,
2212  n_components,
2213  VectorType,
2214  VectorizedArrayType>::apply_add(VectorType & dst,
2215  const VectorType &src) const
2216  {
2218  &LaplaceOperator::local_apply_cell, this, dst, src);
2219  }
2220 
2221  namespace Implementation
2222  {
2223  template <typename VectorizedArrayType>
2224  bool
2225  non_negative(const VectorizedArrayType &n)
2226  {
2227  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2228  if (n[v] < 0.)
2229  return false;
2230 
2231  return true;
2232  }
2233  } // namespace Implementation
2234 
2235 
2236 
2237  template <int dim,
2238  int fe_degree,
2239  int n_q_points_1d,
2240  int n_components,
2241  typename VectorType,
2242  typename VectorizedArrayType>
2243  void
2244  LaplaceOperator<dim,
2245  fe_degree,
2246  n_q_points_1d,
2247  n_components,
2248  VectorType,
2249  VectorizedArrayType>::
2251  FEEvaluation<
2252  dim,
2253  fe_degree,
2254  n_q_points_1d,
2255  n_components,
2256  typename Base<dim, VectorType, VectorizedArrayType>::value_type> &phi,
2257  const unsigned int cell) const
2258  {
2259  phi.evaluate(EvaluationFlags::gradients);
2260  if (scalar_coefficient.get())
2261  {
2262  Assert(scalar_coefficient->size(1) == 1 ||
2263  scalar_coefficient->size(1) == phi.n_q_points,
2264  ExcMessage("The number of columns in the coefficient table must "
2265  "be either 1 or the number of quadrature points " +
2266  std::to_string(phi.n_q_points) +
2267  ", but the given value was " +
2268  std::to_string(scalar_coefficient->size(1))));
2269  if (scalar_coefficient->size(1) == phi.n_q_points)
2270  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2271  {
2273  (*scalar_coefficient)(cell, q)),
2274  ExcMessage("Coefficient must be non-negative"));
2275  phi.submit_gradient((*scalar_coefficient)(cell, q) *
2276  phi.get_gradient(q),
2277  q);
2278  }
2279  else
2280  {
2282  ExcMessage("Coefficient must be non-negative"));
2283  const VectorizedArrayType coefficient =
2284  (*scalar_coefficient)(cell, 0);
2285  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2286  phi.submit_gradient(coefficient * phi.get_gradient(q), q);
2287  }
2288  }
2289  else
2290  {
2291  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2292  {
2293  phi.submit_gradient(phi.get_gradient(q), q);
2294  }
2295  }
2296  phi.integrate(EvaluationFlags::gradients);
2297  }
2298 
2299 
2300 
2301  template <int dim,
2302  int fe_degree,
2303  int n_q_points_1d,
2304  int n_components,
2305  typename VectorType,
2306  typename VectorizedArrayType>
2307  void
2308  LaplaceOperator<dim,
2309  fe_degree,
2310  n_q_points_1d,
2311  n_components,
2312  VectorType,
2313  VectorizedArrayType>::
2315  const MatrixFree<
2316  dim,
2317  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2318  VectorizedArrayType> & data,
2319  VectorType & dst,
2320  const VectorType & src,
2321  const std::pair<unsigned int, unsigned int> &cell_range) const
2322  {
2323  using Number =
2324  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2326  data, this->selected_rows[0]);
2327  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2328  {
2329  phi.reinit(cell);
2330  phi.read_dof_values(src);
2331  do_operation_on_cell(phi, cell);
2332  phi.distribute_local_to_global(dst);
2333  }
2334  }
2335 
2336 
2337  template <int dim,
2338  int fe_degree,
2339  int n_q_points_1d,
2340  int n_components,
2341  typename VectorType,
2342  typename VectorizedArrayType>
2343  void
2344  LaplaceOperator<dim,
2345  fe_degree,
2346  n_q_points_1d,
2347  n_components,
2348  VectorType,
2349  VectorizedArrayType>::
2351  const MatrixFree<
2352  dim,
2353  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2354  VectorizedArrayType> &data,
2355  VectorType & dst,
2356  const VectorType &,
2357  const std::pair<unsigned int, unsigned int> &cell_range) const
2358  {
2359  using Number =
2360  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2361 
2363  data, this->selected_rows[0]);
2364  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2365  {
2366  phi.reinit(cell);
2367  VectorizedArrayType local_diagonal_vector[phi.static_dofs_per_cell];
2368  for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2369  {
2370  for (unsigned int j = 0; j < phi.dofs_per_component; ++j)
2371  phi.begin_dof_values()[j] = VectorizedArrayType();
2372  phi.begin_dof_values()[i] = 1.;
2373  do_operation_on_cell(phi, cell);
2374  local_diagonal_vector[i] = phi.begin_dof_values()[i];
2375  }
2376  for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2377  for (unsigned int c = 0; c < phi.n_components; ++c)
2378  phi.begin_dof_values()[i + c * phi.dofs_per_component] =
2379  local_diagonal_vector[i];
2380  phi.distribute_local_to_global(dst);
2381  }
2382  }
2383 
2384 
2385 } // end of namespace MatrixFreeOperators
2386 
2387 
2389 
2390 #endif
typename OperatorType::value_type value_type
Definition: operators.h:547
virtual void apply_add(VectorType &dst, const VectorType &src) const =0
void fill_inverse_JxW_values(AlignedVector< VectorizedArrayType > &inverse_jxw) const
Definition: operators.h:1031
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:1676
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1655
virtual std::size_t memory_consumption() const
Definition: operators.h:1653
const FEEvaluationBase< dim, n_components, Number, false, VectorizedArrayType > & fe_eval
Definition: operators.h:726
types::global_dof_index size_type
Definition: cuda_kernels.h:45
size_type n() const
Definition: operators.h:1161
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:2350
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_lumped_diagonal_inverse() const
Definition: operators.h:1974
void precondition_Jacobi(VectorType &dst, const VectorType &src, const value_type omega) const
Definition: operators.h:1711
std::shared_ptr< DiagonalMatrix< VectorType > > diagonal_entries
Definition: operators.h:445
void read_dof_values(const VectorType &src, const unsigned int first_index=0)
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_lumped_diagonal() const
Definition: operators.h:1996
#define AssertIndexRange(index, range)
Definition: exceptions.h:1720
std::vector< std::vector< std::pair< value_type, value_type > > > edge_constrained_values
Definition: operators.h:475
static ::ExceptionBase & ExcNotInitialized()
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > data
Definition: operators.h:439
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1639
#define AssertThrow(cond, exc)
Definition: exceptions.h:1571
std::enable_if< IsBlockVector< VectorType >::value, void >::type collect_sizes(VectorType &vector)
Definition: operators.h:98
typename VectorType::value_type value_type
Definition: operators.h:196
typename OperatorType::size_type size_type
Definition: operators.h:552
virtual void compute_diagonal() override
Definition: operators.h:2163
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:1667
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1796
void vmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1396
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:1689
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:1598
std::enable_if< IsBlockVector< VectorType >::value, typename VectorType::BlockType & >::type subblock(VectorType &vector, unsigned int block_no)
Definition: operators.h:66
#define Assert(cond, exc)
Definition: exceptions.h:1461
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:2314
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:1926
void postprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1510
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:1700
std::shared_ptr< Table< 2, VectorizedArrayType > > scalar_coefficient
Definition: operators.h:992
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
unsigned int level
Definition: grid_out.cc:4589
const unsigned int dofs_per_component
std::vector< unsigned int > selected_rows
Definition: operators.h:457
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:50
virtual void clear()
Definition: operators.h:1175
void set_constrained_entries_to_one(VectorType &dst) const
Definition: operators.h:1364
virtual void clear() override
Definition: operators.h:2102
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1200
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1756
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition: operators.h:2214
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition: operators.h:2018
void Tvmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1407
static constexpr unsigned int n_components
size_type m() const
Definition: operators.h:1147
static constexpr unsigned int static_dofs_per_cell
typename VectorType::size_type size_type
Definition: operators.h:201
void mult_add(VectorType &dst, const VectorType &src, const bool transpose) const
Definition: operators.h:1490
std::shared_ptr< Table< 2, VectorizedArrayType > > get_coefficient()
Definition: operators.h:2143
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:2250
std::shared_ptr< DiagonalMatrix< VectorType > > lumped_diagonal_entries
Definition: operators.h:824
const VectorizedArrayType * begin_dof_values() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
std::vector< unsigned int > selected_columns
Definition: operators.h:463
void vmult_interface_down(VectorType &dst, const VectorType &src) const
Definition: operators.h:1544
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:1185
std::vector< UnivariateShapeData< Number > > data
Definition: shape_info.h:392
bool is_element(const size_type index) const
Definition: index_set.h:1770
void initialize(const OperatorType &operator_in)
Definition: operators.h:1746
SmartPointer< const OperatorType > mf_base_operator
Definition: operators.h:597
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1776
void apply(const AlignedVector< VectorizedArrayType > &inverse_coefficient, const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition: operators.h:1083
std::shared_ptr< DiagonalMatrix< VectorType > > inverse_diagonal_entries
Definition: operators.h:451
std::shared_ptr< DiagonalMatrix< VectorType > > inverse_lumped_diagonal_entries
Definition: operators.h:830
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1383
bool non_negative(const VectorizedArrayType &n)
Definition: operators.h:2225
virtual void compute_diagonal() override
Definition: operators.h:1837
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:2040
std::vector< std::vector< unsigned int > > edge_constrained_indices
Definition: operators.h:469
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:1225
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:1460
void transform_from_q_points_to_basis(const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition: operators.h:1110
void set_coefficient(const std::shared_ptr< Table< 2, VectorizedArrayType >> &scalar_coefficient)
Definition: operators.h:2123
static ::ExceptionBase & ExcInternalError()
void adjust_ghost_range_if_necessary(const VectorType &vec, const bool is_row) const
Definition: operators.h:1418