Reference documentation for deal.II version Git abc1da1e07 2020-07-06 12:09:19 -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 - 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  // temporarily reduce size of inverse_jxw to dofs_per_cell to get JxW values
995  // from fe_eval (will not reallocate any memory)
996  for (unsigned int q = 0; q < dofs_per_component_on_cell; ++q)
997  inverse_jxw[q] = 1. / fe_eval.JxW(q);
998  // copy values to rest of vector
999  for (unsigned int q = dofs_per_component_on_cell; q < inverse_jxw.size();)
1000  for (unsigned int i = 0; i < dofs_per_component_on_cell; ++i, ++q)
1001  inverse_jxw[q] = inverse_jxw[i];
1002  }
1003 
1004 
1005 
1006  template <int dim,
1007  int fe_degree,
1008  int n_components,
1009  typename Number,
1010  typename VectorizedArrayType>
1011  inline void
1013  dim,
1014  fe_degree,
1015  n_components,
1016  Number,
1017  VectorizedArrayType>::apply(const VectorizedArrayType *in_array,
1018  VectorizedArrayType * out_array) const
1019  {
1021  dim,
1022  fe_degree,
1023  n_components,
1024  VectorizedArrayType>::apply(fe_eval, in_array, out_array);
1025  }
1026 
1027 
1028 
1029  template <int dim,
1030  int fe_degree,
1031  int n_components,
1032  typename Number,
1033  typename VectorizedArrayType>
1034  inline void
1036  fe_degree,
1037  n_components,
1038  Number,
1039  VectorizedArrayType>::
1040  apply(const AlignedVector<VectorizedArrayType> &inverse_coefficients,
1041  const unsigned int n_actual_components,
1042  const VectorizedArrayType * in_array,
1043  VectorizedArrayType * out_array) const
1044  {
1046  fe_degree,
1047  n_components,
1048  VectorizedArrayType>::
1049  apply(fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
1050  inverse_coefficients,
1051  n_actual_components,
1052  in_array,
1053  out_array);
1054  }
1055 
1056 
1057 
1058  template <int dim,
1059  int fe_degree,
1060  int n_components,
1061  typename Number,
1062  typename VectorizedArrayType>
1063  inline void
1065  fe_degree,
1066  n_components,
1067  Number,
1068  VectorizedArrayType>::
1069  transform_from_q_points_to_basis(const unsigned int n_actual_components,
1070  const VectorizedArrayType *in_array,
1071  VectorizedArrayType * out_array) const
1072  {
1074  fe_degree,
1075  n_components,
1076  VectorizedArrayType>::
1078  fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
1079  n_actual_components,
1080  in_array,
1081  out_array);
1082  }
1083 
1084 
1085 
1086  //----------------- Base operator -----------------------------
1087  template <int dim, typename VectorType, typename VectorizedArrayType>
1089  : Subscriptor()
1090  , have_interface_matrices(false)
1091  {}
1092 
1093 
1094 
1095  template <int dim, typename VectorType, typename VectorizedArrayType>
1098  {
1099  Assert(data.get() != nullptr, ExcNotInitialized());
1101  0;
1102  for (const unsigned int selected_row : selected_rows)
1103  total_size += data->get_vector_partitioner(selected_row)->size();
1104  return total_size;
1105  }
1106 
1107 
1108 
1109  template <int dim, typename VectorType, typename VectorizedArrayType>
1112  {
1113  Assert(data.get() != nullptr, ExcNotInitialized());
1115  0;
1116  for (const unsigned int selected_column : selected_columns)
1117  total_size += data->get_vector_partitioner(selected_column)->size();
1118  return total_size;
1119  }
1120 
1121 
1122 
1123  template <int dim, typename VectorType, typename VectorizedArrayType>
1124  void
1126  {
1127  data.reset();
1128  inverse_diagonal_entries.reset();
1129  }
1130 
1131 
1132 
1133  template <int dim, typename VectorType, typename VectorizedArrayType>
1136  const unsigned int col) const
1137  {
1138  (void)col;
1139  Assert(row == col, ExcNotImplemented());
1140  Assert(inverse_diagonal_entries.get() != nullptr &&
1141  inverse_diagonal_entries->m() > 0,
1142  ExcNotInitialized());
1143  return 1.0 / (*inverse_diagonal_entries)(row, row);
1144  }
1145 
1146 
1147 
1148  template <int dim, typename VectorType, typename VectorizedArrayType>
1149  void
1151  VectorType &vec) const
1152  {
1153  Assert(data.get() != nullptr, ExcNotInitialized());
1155  for (unsigned int i = 0; i < BlockHelper::n_blocks(vec); ++i)
1156  {
1157  const unsigned int index = selected_rows[i];
1158  if (!BlockHelper::subblock(vec, index)
1159  .partitioners_are_compatible(
1160  *data->get_dof_info(index).vector_partitioner))
1161  data->initialize_dof_vector(BlockHelper::subblock(vec, index), index);
1162 
1163  Assert(BlockHelper::subblock(vec, index)
1164  .partitioners_are_globally_compatible(
1165  *data->get_dof_info(index).vector_partitioner),
1166  ExcInternalError());
1167  }
1169  }
1170 
1171 
1172 
1173  template <int dim, typename VectorType, typename VectorizedArrayType>
1174  void
1177  data_,
1178  const std::vector<unsigned int> &given_row_selection,
1179  const std::vector<unsigned int> &given_column_selection)
1180  {
1181  data = data_;
1182 
1183  selected_rows.clear();
1184  selected_columns.clear();
1185  if (given_row_selection.empty())
1186  for (unsigned int i = 0; i < data_->n_components(); ++i)
1187  selected_rows.push_back(i);
1188  else
1189  {
1190  for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1191  {
1192  AssertIndexRange(given_row_selection[i], data_->n_components());
1193  for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1194  if (j != i)
1195  Assert(given_row_selection[j] != given_row_selection[i],
1196  ExcMessage("Given row indices must be unique"));
1197 
1198  selected_rows.push_back(given_row_selection[i]);
1199  }
1200  }
1201  if (given_column_selection.size() == 0)
1203  else
1204  {
1205  for (unsigned int i = 0; i < given_column_selection.size(); ++i)
1206  {
1207  AssertIndexRange(given_column_selection[i], data_->n_components());
1208  for (unsigned int j = 0; j < given_column_selection.size(); ++j)
1209  if (j != i)
1210  Assert(given_column_selection[j] != given_column_selection[i],
1211  ExcMessage("Given column indices must be unique"));
1212 
1213  selected_columns.push_back(given_column_selection[i]);
1214  }
1215  }
1216 
1217  edge_constrained_indices.clear();
1218  edge_constrained_indices.resize(selected_rows.size());
1219  edge_constrained_values.clear();
1220  edge_constrained_values.resize(selected_rows.size());
1221  have_interface_matrices = false;
1222  }
1223 
1224 
1225 
1226  template <int dim, typename VectorType, typename VectorizedArrayType>
1227  void
1230  data_,
1231  const MGConstrainedDoFs & mg_constrained_dofs,
1232  const unsigned int level,
1233  const std::vector<unsigned int> &given_row_selection)
1234  {
1235  std::vector<MGConstrainedDoFs> mg_constrained_dofs_vector(
1236  1, mg_constrained_dofs);
1237  initialize(data_, mg_constrained_dofs_vector, level, given_row_selection);
1238  }
1239 
1240 
1241 
1242  template <int dim, typename VectorType, typename VectorizedArrayType>
1243  void
1246  data_,
1247  const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
1248  const unsigned int level,
1249  const std::vector<unsigned int> & given_row_selection)
1250  {
1252  ExcMessage("level is not set"));
1253 
1254  selected_rows.clear();
1255  selected_columns.clear();
1256  if (given_row_selection.empty())
1257  for (unsigned int i = 0; i < data_->n_components(); ++i)
1258  selected_rows.push_back(i);
1259  else
1260  {
1261  for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1262  {
1263  AssertIndexRange(given_row_selection[i], data_->n_components());
1264  for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1265  if (j != i)
1266  Assert(given_row_selection[j] != given_row_selection[i],
1267  ExcMessage("Given row indices must be unique"));
1268 
1269  selected_rows.push_back(given_row_selection[i]);
1270  }
1271  }
1273 
1274  AssertDimension(mg_constrained_dofs.size(), selected_rows.size());
1275  edge_constrained_indices.clear();
1276  edge_constrained_indices.resize(selected_rows.size());
1277  edge_constrained_values.clear();
1278  edge_constrained_values.resize(selected_rows.size());
1279 
1280  data = data_;
1281 
1282  for (unsigned int j = 0; j < selected_rows.size(); ++j)
1283  {
1284  if (data_->n_macro_cells() > 0)
1285  {
1286  AssertDimension(level, data_->get_cell_iterator(0, 0, j)->level());
1287  }
1288 
1289  // setup edge_constrained indices
1290  std::vector<types::global_dof_index> interface_indices;
1291  mg_constrained_dofs[j]
1292  .get_refinement_edge_indices(level)
1293  .fill_index_vector(interface_indices);
1294  edge_constrained_indices[j].clear();
1295  edge_constrained_indices[j].reserve(interface_indices.size());
1296  edge_constrained_values[j].resize(interface_indices.size());
1297  const IndexSet &locally_owned =
1298  data->get_dof_handler(selected_rows[j]).locally_owned_mg_dofs(level);
1299  for (const auto interface_index : interface_indices)
1300  if (locally_owned.is_element(interface_index))
1301  edge_constrained_indices[j].push_back(
1302  locally_owned.index_within_set(interface_index));
1305  static_cast<unsigned int>(edge_constrained_indices[j].size()),
1306  data->get_vector_partitioner()->get_mpi_communicator()) > 0;
1307  }
1308  }
1309 
1310 
1311 
1312  template <int dim, typename VectorType, typename VectorizedArrayType>
1313  void
1315  VectorType &dst) const
1316  {
1317  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1318  {
1319  const std::vector<unsigned int> &constrained_dofs =
1320  data->get_constrained_dofs(selected_rows[j]);
1321  for (const auto constrained_dof : constrained_dofs)
1322  BlockHelper::subblock(dst, j).local_element(constrained_dof) = 1.;
1323  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1324  BlockHelper::subblock(dst, j).local_element(
1325  edge_constrained_indices[j][i]) = 1.;
1326  }
1327  }
1328 
1329 
1330 
1331  template <int dim, typename VectorType, typename VectorizedArrayType>
1332  void
1334  const VectorType &src) const
1335  {
1336  using Number =
1338  dst = Number(0.);
1339  vmult_add(dst, src);
1340  }
1341 
1342 
1343 
1344  template <int dim, typename VectorType, typename VectorizedArrayType>
1345  void
1347  VectorType & dst,
1348  const VectorType &src) const
1349  {
1350  mult_add(dst, src, false);
1351  }
1352 
1353 
1354 
1355  template <int dim, typename VectorType, typename VectorizedArrayType>
1356  void
1358  VectorType & dst,
1359  const VectorType &src) const
1360  {
1361  mult_add(dst, src, true);
1362  }
1363 
1364 
1365 
1366  template <int dim, typename VectorType, typename VectorizedArrayType>
1367  void
1369  const VectorType &src,
1370  const bool is_row) const
1371  {
1372  using Number =
1374  for (unsigned int i = 0; i < BlockHelper::n_blocks(src); ++i)
1375  {
1376  const unsigned int mf_component =
1377  is_row ? selected_rows[i] : selected_columns[i];
1378  // If both vectors use the same partitioner -> done
1379  if (BlockHelper::subblock(src, i).get_partitioner().get() ==
1380  data->get_dof_info(mf_component).vector_partitioner.get())
1381  continue;
1382 
1383  // If not, assert that the local ranges are the same and reset to the
1384  // current partitioner
1385  Assert(
1386  BlockHelper::subblock(src, i).get_partitioner()->local_size() ==
1387  data->get_dof_info(mf_component).vector_partitioner->local_size(),
1388  ExcMessage("The vector passed to the vmult() function does not have "
1389  "the correct size for compatibility with MatrixFree."));
1390 
1391  // copy the vector content to a temporary vector so that it does not get
1392  // lost
1394  BlockHelper::subblock(src, i));
1395  BlockHelper::subblock(const_cast<VectorType &>(src), i)
1396  .reinit(data->get_dof_info(mf_component).vector_partitioner);
1397  BlockHelper::subblock(const_cast<VectorType &>(src), i)
1398  .copy_locally_owned_data_from(copy_vec);
1399  }
1400  }
1401 
1402 
1403 
1404  template <int dim, typename VectorType, typename VectorizedArrayType>
1405  void
1407  VectorType & dst,
1408  const VectorType &src) const
1409  {
1410  using Number =
1412  adjust_ghost_range_if_necessary(src, false);
1414 
1415  // set zero Dirichlet values on the input vector (and remember the src and
1416  // dst values because we need to reset them at the end)
1417  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1418  {
1419  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1420  {
1421  edge_constrained_values[j][i] = std::pair<Number, Number>(
1422  BlockHelper::subblock(src, j).local_element(
1423  edge_constrained_indices[j][i]),
1424  BlockHelper::subblock(dst, j).local_element(
1425  edge_constrained_indices[j][i]));
1426  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1427  .local_element(edge_constrained_indices[j][i]) = 0.;
1428  }
1429  }
1430  }
1431 
1432 
1433 
1434  template <int dim, typename VectorType, typename VectorizedArrayType>
1435  void
1437  VectorType & dst,
1438  const VectorType &src,
1439  const bool transpose) const
1440  {
1441  AssertDimension(dst.size(), src.size());
1444  preprocess_constraints(dst, src);
1445  if (transpose)
1446  Tapply_add(dst, src);
1447  else
1448  apply_add(dst, src);
1449  postprocess_constraints(dst, src);
1450  }
1451 
1452 
1453 
1454  template <int dim, typename VectorType, typename VectorizedArrayType>
1455  void
1457  VectorType & dst,
1458  const VectorType &src) const
1459  {
1460  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1461  {
1462  const std::vector<unsigned int> &constrained_dofs =
1463  data->get_constrained_dofs(selected_rows[j]);
1464  for (const auto constrained_dof : constrained_dofs)
1465  BlockHelper::subblock(dst, j).local_element(constrained_dof) +=
1466  BlockHelper::subblock(src, j).local_element(constrained_dof);
1467  }
1468 
1469  // reset edge constrained values, multiply by unit matrix and add into
1470  // destination
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  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1476  .local_element(edge_constrained_indices[j][i]) =
1477  edge_constrained_values[j][i].first;
1478  BlockHelper::subblock(dst, j).local_element(
1479  edge_constrained_indices[j][i]) =
1480  edge_constrained_values[j][i].second +
1481  edge_constrained_values[j][i].first;
1482  }
1483  }
1484  }
1485 
1486 
1487 
1488  template <int dim, typename VectorType, typename VectorizedArrayType>
1489  void
1491  VectorType & dst,
1492  const VectorType &src) const
1493  {
1494  using Number =
1496  AssertDimension(dst.size(), src.size());
1497  adjust_ghost_range_if_necessary(src, false);
1499 
1500  dst = Number(0.);
1501 
1503  return;
1504 
1505  // set zero Dirichlet values on the input vector (and remember the src and
1506  // dst values because we need to reset them at the end)
1507  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1508  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1509  {
1510  edge_constrained_values[j][i] = std::pair<Number, Number>(
1511  BlockHelper::subblock(src, j).local_element(
1512  edge_constrained_indices[j][i]),
1513  BlockHelper::subblock(dst, j).local_element(
1514  edge_constrained_indices[j][i]));
1515  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1516  .local_element(edge_constrained_indices[j][i]) = 0.;
1517  }
1518 
1519  apply_add(dst, src);
1520 
1521  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1522  {
1523  unsigned int c = 0;
1524  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1525  {
1526  for (; c < edge_constrained_indices[j][i]; ++c)
1527  BlockHelper::subblock(dst, j).local_element(c) = 0.;
1528  ++c;
1529 
1530  // reset the src values
1531  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1532  .local_element(edge_constrained_indices[j][i]) =
1533  edge_constrained_values[j][i].first;
1534  }
1535  for (; c < BlockHelper::subblock(dst, j).local_size(); ++c)
1536  BlockHelper::subblock(dst, j).local_element(c) = 0.;
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  VectorType src_cpy(src);
1560  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1561  {
1562  unsigned int c = 0;
1563  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1564  {
1565  for (; c < edge_constrained_indices[j][i]; ++c)
1566  BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1567  ++c;
1568  }
1569  for (; c < BlockHelper::subblock(src_cpy, j).local_size(); ++c)
1570  BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1571  }
1572 
1573  apply_add(dst, src_cpy);
1574 
1575  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1576  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1577  BlockHelper::subblock(dst, j).local_element(
1578  edge_constrained_indices[j][i]) = 0.;
1579  }
1580 
1581 
1582 
1583  template <int dim, typename VectorType, typename VectorizedArrayType>
1584  void
1586  VectorType & dst,
1587  const VectorType &src) const
1588  {
1589  using Number =
1591  dst = Number(0.);
1592  Tvmult_add(dst, src);
1593  }
1594 
1595 
1596 
1597  template <int dim, typename VectorType, typename VectorizedArrayType>
1598  std::size_t
1600  {
1601  return inverse_diagonal_entries.get() != nullptr ?
1602  inverse_diagonal_entries->memory_consumption() :
1603  sizeof(*this);
1604  }
1605 
1606 
1607 
1608  template <int dim, typename VectorType, typename VectorizedArrayType>
1609  std::shared_ptr<const MatrixFree<
1610  dim,
1612  VectorizedArrayType>>
1614  {
1615  return data;
1616  }
1617 
1618 
1619 
1620  template <int dim, typename VectorType, typename VectorizedArrayType>
1621  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1623  const
1624  {
1625  Assert(inverse_diagonal_entries.get() != nullptr &&
1626  inverse_diagonal_entries->m() > 0,
1627  ExcNotInitialized());
1628  return inverse_diagonal_entries;
1629  }
1630 
1631 
1632 
1633  template <int dim, typename VectorType, typename VectorizedArrayType>
1634  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1636  {
1637  Assert(diagonal_entries.get() != nullptr && diagonal_entries->m() > 0,
1638  ExcNotInitialized());
1639  return diagonal_entries;
1640  }
1641 
1642 
1643 
1644  template <int dim, typename VectorType, typename VectorizedArrayType>
1645  void
1647  VectorType & dst,
1648  const VectorType &src) const
1649  {
1650  apply_add(dst, src);
1651  }
1652 
1653 
1654 
1655  template <int dim, typename VectorType, typename VectorizedArrayType>
1656  void
1658  VectorType & dst,
1659  const VectorType & src,
1660  const typename Base<dim, VectorType, VectorizedArrayType>::value_type omega)
1661  const
1662  {
1664  ExcNotInitialized());
1665  inverse_diagonal_entries->vmult(dst, src);
1666  dst *= omega;
1667  }
1668 
1669 
1670 
1671  //------------------------- MGInterfaceOperator ------------------------------
1672 
1673  template <typename OperatorType>
1675  : Subscriptor()
1676  , mf_base_operator(nullptr)
1677  {}
1678 
1679 
1680 
1681  template <typename OperatorType>
1682  void
1684  {
1685  mf_base_operator = nullptr;
1686  }
1687 
1688 
1689 
1690  template <typename OperatorType>
1691  void
1692  MGInterfaceOperator<OperatorType>::initialize(const OperatorType &operator_in)
1693  {
1694  mf_base_operator = &operator_in;
1695  }
1696 
1697 
1698 
1699  template <typename OperatorType>
1700  template <typename VectorType>
1701  void
1703  const VectorType &src) const
1704  {
1705 #ifndef DEAL_II_MSVC
1706  static_assert(
1707  std::is_same<typename VectorType::value_type, value_type>::value,
1708  "The vector type must be based on the same value type as this "
1709  "operator");
1710 #endif
1711 
1712  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1713 
1714  mf_base_operator->vmult_interface_down(dst, src);
1715  }
1716 
1717 
1718 
1719  template <typename OperatorType>
1720  template <typename VectorType>
1721  void
1723  const VectorType &src) const
1724  {
1725 #ifndef DEAL_II_MSVC
1726  static_assert(
1727  std::is_same<typename VectorType::value_type, value_type>::value,
1728  "The vector type must be based on the same value type as this "
1729  "operator");
1730 #endif
1731 
1732  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1733 
1734  mf_base_operator->vmult_interface_up(dst, src);
1735  }
1736 
1737 
1738 
1739  template <typename OperatorType>
1740  template <typename VectorType>
1741  void
1743  VectorType &vec) const
1744  {
1745  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1746 
1747  mf_base_operator->initialize_dof_vector(vec);
1748  }
1749 
1750 
1751 
1752  //-----------------------------MassOperator----------------------------------
1753 
1754  template <int dim,
1755  int fe_degree,
1756  int n_q_points_1d,
1757  int n_components,
1758  typename VectorType,
1759  typename VectorizedArrayType>
1760  MassOperator<dim,
1761  fe_degree,
1762  n_q_points_1d,
1763  n_components,
1764  VectorType,
1765  VectorizedArrayType>::MassOperator()
1766  : Base<dim, VectorType, VectorizedArrayType>()
1767  {}
1768 
1769 
1770 
1771  template <int dim,
1772  int fe_degree,
1773  int n_q_points_1d,
1774  int n_components,
1775  typename VectorType,
1776  typename VectorizedArrayType>
1777  void
1778  MassOperator<dim,
1779  fe_degree,
1780  n_q_points_1d,
1781  n_components,
1782  VectorType,
1783  VectorizedArrayType>::compute_diagonal()
1784  {
1785  using Number =
1786  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1788  ExcNotInitialized());
1789 
1790  this->inverse_diagonal_entries =
1791  std::make_shared<DiagonalMatrix<VectorType>>();
1792  this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1793  VectorType &inverse_diagonal_vector =
1794  this->inverse_diagonal_entries->get_vector();
1795  VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1796  this->initialize_dof_vector(inverse_diagonal_vector);
1797  this->initialize_dof_vector(diagonal_vector);
1798  inverse_diagonal_vector = Number(1.);
1799  apply_add(diagonal_vector, inverse_diagonal_vector);
1800 
1801  this->set_constrained_entries_to_one(diagonal_vector);
1802  inverse_diagonal_vector = diagonal_vector;
1803 
1804  const unsigned int local_size = inverse_diagonal_vector.local_size();
1805  for (unsigned int i = 0; i < local_size; ++i)
1806  inverse_diagonal_vector.local_element(i) =
1807  Number(1.) / inverse_diagonal_vector.local_element(i);
1808 
1809  inverse_diagonal_vector.update_ghost_values();
1810  diagonal_vector.update_ghost_values();
1811  }
1812 
1813 
1814 
1815  template <int dim,
1816  int fe_degree,
1817  int n_q_points_1d,
1818  int n_components,
1819  typename VectorType,
1820  typename VectorizedArrayType>
1821  void
1822  MassOperator<dim,
1823  fe_degree,
1824  n_q_points_1d,
1825  n_components,
1826  VectorType,
1827  VectorizedArrayType>::apply_add(VectorType & dst,
1828  const VectorType &src) const
1829  {
1831  &MassOperator::local_apply_cell, this, dst, src);
1832  }
1833 
1834 
1835 
1836  template <int dim,
1837  int fe_degree,
1838  int n_q_points_1d,
1839  int n_components,
1840  typename VectorType,
1841  typename VectorizedArrayType>
1842  void
1843  MassOperator<dim,
1844  fe_degree,
1845  n_q_points_1d,
1846  n_components,
1847  VectorType,
1848  VectorizedArrayType>::
1850  const MatrixFree<
1851  dim,
1852  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
1853  VectorizedArrayType> & data,
1854  VectorType & dst,
1855  const VectorType & src,
1856  const std::pair<unsigned int, unsigned int> &cell_range) const
1857  {
1858  using Number =
1859  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1860  FEEvaluation<dim,
1861  fe_degree,
1862  n_q_points_1d,
1863  n_components,
1864  Number,
1865  VectorizedArrayType>
1866  phi(data, this->selected_rows[0]);
1867  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
1868  {
1869  phi.reinit(cell);
1870  phi.read_dof_values(src);
1871  phi.evaluate(EvaluationFlags::values);
1872  for (unsigned int q = 0; q < phi.n_q_points; ++q)
1873  phi.submit_value(phi.get_value(q), q);
1874  phi.integrate(EvaluationFlags::values);
1875  phi.distribute_local_to_global(dst);
1876  }
1877  }
1878 
1879 
1880  //-----------------------------LaplaceOperator----------------------------------
1881 
1882  template <int dim,
1883  int fe_degree,
1884  int n_q_points_1d,
1885  int n_components,
1886  typename VectorType,
1887  typename VectorizedArrayType>
1888  LaplaceOperator<dim,
1889  fe_degree,
1890  n_q_points_1d,
1891  n_components,
1892  VectorType,
1893  VectorizedArrayType>::LaplaceOperator()
1894  : Base<dim, VectorType, VectorizedArrayType>()
1895  {}
1896 
1897 
1898 
1899  template <int dim,
1900  int fe_degree,
1901  int n_q_points_1d,
1902  int n_components,
1903  typename VectorType,
1904  typename VectorizedArrayType>
1905  void
1906  LaplaceOperator<dim,
1907  fe_degree,
1908  n_q_points_1d,
1909  n_components,
1910  VectorType,
1911  VectorizedArrayType>::clear()
1912  {
1914  scalar_coefficient.reset();
1915  }
1916 
1917 
1918 
1919  template <int dim,
1920  int fe_degree,
1921  int n_q_points_1d,
1922  int n_components,
1923  typename VectorType,
1924  typename VectorizedArrayType>
1925  void
1926  LaplaceOperator<dim,
1927  fe_degree,
1928  n_q_points_1d,
1929  n_components,
1930  VectorType,
1931  VectorizedArrayType>::
1933  const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient_)
1934  {
1935  scalar_coefficient = scalar_coefficient_;
1936  }
1937 
1938 
1939 
1940  template <int dim,
1941  int fe_degree,
1942  int n_q_points_1d,
1943  int n_components,
1944  typename VectorType,
1945  typename VectorizedArrayType>
1946  std::shared_ptr<Table<2, VectorizedArrayType>>
1947  LaplaceOperator<dim,
1948  fe_degree,
1949  n_q_points_1d,
1950  n_components,
1951  VectorType,
1952  VectorizedArrayType>::get_coefficient()
1953  {
1955  return scalar_coefficient;
1956  }
1957 
1958 
1959 
1960  template <int dim,
1961  int fe_degree,
1962  int n_q_points_1d,
1963  int n_components,
1964  typename VectorType,
1965  typename VectorizedArrayType>
1966  void
1967  LaplaceOperator<dim,
1968  fe_degree,
1969  n_q_points_1d,
1970  n_components,
1971  VectorType,
1972  VectorizedArrayType>::compute_diagonal()
1973  {
1974  using Number =
1975  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1977  ExcNotInitialized());
1978 
1979  this->inverse_diagonal_entries =
1980  std::make_shared<DiagonalMatrix<VectorType>>();
1981  this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1982  VectorType &inverse_diagonal_vector =
1983  this->inverse_diagonal_entries->get_vector();
1984  VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1985  this->initialize_dof_vector(inverse_diagonal_vector);
1986  this->initialize_dof_vector(diagonal_vector);
1987 
1988  this->data->cell_loop(&LaplaceOperator::local_diagonal_cell,
1989  this,
1990  diagonal_vector,
1991  /*unused*/ diagonal_vector);
1992  this->set_constrained_entries_to_one(diagonal_vector);
1993 
1994  inverse_diagonal_vector = diagonal_vector;
1995 
1996  for (unsigned int i = 0; i < inverse_diagonal_vector.local_size(); ++i)
1997  if (std::abs(inverse_diagonal_vector.local_element(i)) >
1999  inverse_diagonal_vector.local_element(i) =
2000  1. / inverse_diagonal_vector.local_element(i);
2001  else
2002  inverse_diagonal_vector.local_element(i) = 1.;
2003 
2004  inverse_diagonal_vector.update_ghost_values();
2005  diagonal_vector.update_ghost_values();
2006  }
2007 
2008 
2009 
2010  template <int dim,
2011  int fe_degree,
2012  int n_q_points_1d,
2013  int n_components,
2014  typename VectorType,
2015  typename VectorizedArrayType>
2016  void
2017  LaplaceOperator<dim,
2018  fe_degree,
2019  n_q_points_1d,
2020  n_components,
2021  VectorType,
2022  VectorizedArrayType>::apply_add(VectorType & dst,
2023  const VectorType &src) const
2024  {
2026  &LaplaceOperator::local_apply_cell, this, dst, src);
2027  }
2028 
2029  namespace Implementation
2030  {
2031  template <typename VectorizedArrayType>
2032  bool
2033  non_negative(const VectorizedArrayType &n)
2034  {
2035  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2036  if (n[v] < 0.)
2037  return false;
2038 
2039  return true;
2040  }
2041  } // namespace Implementation
2042 
2043 
2044 
2045  template <int dim,
2046  int fe_degree,
2047  int n_q_points_1d,
2048  int n_components,
2049  typename VectorType,
2050  typename VectorizedArrayType>
2051  void
2052  LaplaceOperator<dim,
2053  fe_degree,
2054  n_q_points_1d,
2055  n_components,
2056  VectorType,
2057  VectorizedArrayType>::
2059  FEEvaluation<
2060  dim,
2061  fe_degree,
2062  n_q_points_1d,
2063  n_components,
2064  typename Base<dim, VectorType, VectorizedArrayType>::value_type> &phi,
2065  const unsigned int cell) const
2066  {
2067  phi.evaluate(EvaluationFlags::gradients);
2068  if (scalar_coefficient.get())
2069  {
2070  Assert(scalar_coefficient->size(1) == 1 ||
2071  scalar_coefficient->size(1) == phi.n_q_points,
2072  ExcMessage("The number of columns in the coefficient table must "
2073  "be either 1 or the number of quadrature points " +
2074  std::to_string(phi.n_q_points) +
2075  ", but the given value was " +
2076  std::to_string(scalar_coefficient->size(1))));
2077  if (scalar_coefficient->size(1) == phi.n_q_points)
2078  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2079  {
2081  (*scalar_coefficient)(cell, q)),
2082  ExcMessage("Coefficient must be non-negative"));
2083  phi.submit_gradient((*scalar_coefficient)(cell, q) *
2084  phi.get_gradient(q),
2085  q);
2086  }
2087  else
2088  {
2090  ExcMessage("Coefficient must be non-negative"));
2091  const VectorizedArrayType coefficient =
2092  (*scalar_coefficient)(cell, 0);
2093  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2094  phi.submit_gradient(coefficient * phi.get_gradient(q), q);
2095  }
2096  }
2097  else
2098  {
2099  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2100  {
2101  phi.submit_gradient(phi.get_gradient(q), q);
2102  }
2103  }
2104  phi.integrate(EvaluationFlags::gradients);
2105  }
2106 
2107 
2108 
2109  template <int dim,
2110  int fe_degree,
2111  int n_q_points_1d,
2112  int n_components,
2113  typename VectorType,
2114  typename VectorizedArrayType>
2115  void
2116  LaplaceOperator<dim,
2117  fe_degree,
2118  n_q_points_1d,
2119  n_components,
2120  VectorType,
2121  VectorizedArrayType>::
2123  const MatrixFree<
2124  dim,
2125  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2126  VectorizedArrayType> & data,
2127  VectorType & dst,
2128  const VectorType & src,
2129  const std::pair<unsigned int, unsigned int> &cell_range) const
2130  {
2131  using Number =
2132  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2134  data, this->selected_rows[0]);
2135  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2136  {
2137  phi.reinit(cell);
2138  phi.read_dof_values(src);
2139  do_operation_on_cell(phi, cell);
2140  phi.distribute_local_to_global(dst);
2141  }
2142  }
2143 
2144 
2145  template <int dim,
2146  int fe_degree,
2147  int n_q_points_1d,
2148  int n_components,
2149  typename VectorType,
2150  typename VectorizedArrayType>
2151  void
2152  LaplaceOperator<dim,
2153  fe_degree,
2154  n_q_points_1d,
2155  n_components,
2156  VectorType,
2157  VectorizedArrayType>::
2159  const MatrixFree<
2160  dim,
2161  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2162  VectorizedArrayType> &data,
2163  VectorType & dst,
2164  const VectorType &,
2165  const std::pair<unsigned int, unsigned int> &cell_range) const
2166  {
2167  using Number =
2168  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2169 
2171  data, this->selected_rows[0]);
2172  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2173  {
2174  phi.reinit(cell);
2175  VectorizedArrayType local_diagonal_vector[phi.static_dofs_per_cell];
2176  for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2177  {
2178  for (unsigned int j = 0; j < phi.dofs_per_component; ++j)
2179  phi.begin_dof_values()[j] = VectorizedArrayType();
2180  phi.begin_dof_values()[i] = 1.;
2181  do_operation_on_cell(phi, cell);
2182  local_diagonal_vector[i] = phi.begin_dof_values()[i];
2183  }
2184  for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2185  for (unsigned int c = 0; c < phi.n_components; ++c)
2186  phi.begin_dof_values()[i + c * phi.dofs_per_component] =
2187  local_diagonal_vector[i];
2188  phi.distribute_local_to_global(dst);
2189  }
2190  }
2191 
2192 
2193 } // end of namespace MatrixFreeOperators
2194 
2195 
2197 
2198 #endif
typename OperatorType::value_type value_type
Definition: operators.h:541
VectorizedArrayType JxW(const unsigned int q_index) const
virtual void apply_add(VectorType &dst, const VectorType &src) const =0
void fill_inverse_JxW_values(AlignedVector< VectorizedArrayType > &inverse_jxw) const
Definition: operators.h:984
static const unsigned int invalid_unsigned_int
Definition: types.h:191
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal_inverse() const
Definition: operators.h:1622
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1560
virtual std::size_t memory_consumption() const
Definition: operators.h:1599
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:1111
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:2158
void precondition_Jacobi(VectorType &dst, const VectorType &src, const value_type omega) const
Definition: operators.h:1657
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:1628
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:1585
#define AssertThrow(cond, exc)
Definition: exceptions.h:1513
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:1972
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:1613
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1742
void vmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1346
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal() const
Definition: operators.h:1635
static ::ExceptionBase & ExcMessage(std::string arg1)
void vmult_interface_up(VectorType &dst, const VectorType &src) const
Definition: operators.h:1544
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:1403
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:2122
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:1919
void postprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1456
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:1646
std::shared_ptr< Table< 2, VectorizedArrayType > > scalar_coefficient
Definition: operators.h:945
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
unsigned int level
Definition: grid_out.cc:4355
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:2341
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:1125
void set_constrained_entries_to_one(VectorType &dst) const
Definition: operators.h:1314
virtual void clear() override
Definition: operators.h:1911
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1150
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1702
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition: operators.h:2022
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition: operators.h:1827
void Tvmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1357
static constexpr unsigned int n_components
size_type m() const
Definition: operators.h:1097
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:1436
std::shared_ptr< Table< 2, VectorizedArrayType > > get_coefficient()
Definition: operators.h:1952
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:2058
const VectorizedArrayType * begin_dof_values() const
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
std::vector< unsigned int > selected_columns
Definition: operators.h:457
void vmult_interface_down(VectorType &dst, const VectorType &src) const
Definition: operators.h:1490
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:1135
Definition: table.h:687
std::vector< UnivariateShapeData< Number > > data
Definition: shape_info.h:351
bool is_element(const size_type index) const
Definition: index_set.h:1763
void initialize(const OperatorType &operator_in)
Definition: operators.h:1692
SmartPointer< const OperatorType > mf_base_operator
Definition: operators.h:591
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1722
void apply(const AlignedVector< VectorizedArrayType > &inverse_coefficient, const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition: operators.h:1040
std::shared_ptr< DiagonalMatrix< VectorType > > inverse_diagonal_entries
Definition: operators.h:445
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1333
bool non_negative(const VectorizedArrayType &n)
Definition: operators.h:2033
virtual void compute_diagonal() override
Definition: operators.h:1783
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:1849
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:1175
void preprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1406
void transform_from_q_points_to_basis(const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition: operators.h:1069
void set_coefficient(const std::shared_ptr< Table< 2, VectorizedArrayType >> &scalar_coefficient)
Definition: operators.h:1932
static ::ExceptionBase & ExcInternalError()
void adjust_ghost_range_if_necessary(const VectorType &vec, const bool is_row) const
Definition: operators.h:1368