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