Reference documentation for deal.II version 9.2.0
\(\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>
48  unsigned int>::type
49  n_blocks(const VectorType &vector)
50  {
51  return vector.n_blocks();
52  }
53 
54  template <typename VectorType>
56  unsigned int>::type
58  {
59  return 1;
60  }
61 
62  template <typename VectorType>
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>
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>
81  VectorType &>::type
82  subblock(VectorType &vector, unsigned int)
83  {
84  return vector;
85  }
86 
87  template <typename VectorType>
89  const VectorType &>::type
90  subblock(const VectorType &vector, unsigned int)
91  {
92  return vector;
93  }
94 
95  template <typename VectorType>
98  {
99  vector.collect_sizes();
100  }
101 
102  template <typename VectorType>
105  {}
106  } // namespace BlockHelper
107 
187  template <int dim,
189  typename VectorizedArrayType =
191  class Base : public Subscriptor
192  {
193  public:
197  using value_type = typename VectorType::value_type;
198 
203 
207  Base();
208 
212  virtual ~Base() override = default;
213 
218  virtual void
219  clear();
220 
237  void
238  initialize(std::shared_ptr<
240  const std::vector<unsigned int> &selected_row_blocks =
241  std::vector<unsigned int>(),
242  const std::vector<unsigned int> &selected_column_blocks =
243  std::vector<unsigned int>());
244 
258  void
259  initialize(std::shared_ptr<
261  const MGConstrainedDoFs & mg_constrained_dofs,
262  const unsigned int level,
263  const std::vector<unsigned int> &selected_row_blocks =
264  std::vector<unsigned int>());
265 
280  void
281  initialize(std::shared_ptr<
283  const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
284  const unsigned int level,
285  const std::vector<unsigned int> & selected_row_blocks =
286  std::vector<unsigned int>());
287 
291  size_type
292  m() const;
293 
297  size_type
298  n() const;
299 
303  void
304  vmult_interface_down(VectorType &dst, const VectorType &src) const;
305 
309  void
310  vmult_interface_up(VectorType &dst, const VectorType &src) const;
311 
315  void
316  vmult(VectorType &dst, const VectorType &src) const;
317 
321  void
322  Tvmult(VectorType &dst, const VectorType &src) const;
323 
327  void
328  vmult_add(VectorType &dst, const VectorType &src) const;
329 
333  void
334  Tvmult_add(VectorType &dst, const VectorType &src) const;
335 
340  value_type
341  el(const unsigned int row, const unsigned int col) const;
342 
347  virtual std::size_t
348  memory_consumption() const;
349 
353  void
354  initialize_dof_vector(VectorType &vec) const;
355 
363  virtual void
364  compute_diagonal() = 0;
365 
369  std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
370  get_matrix_free() const;
371 
375  const std::shared_ptr<DiagonalMatrix<VectorType>> &
376  get_matrix_diagonal_inverse() const;
377 
381  const std::shared_ptr<DiagonalMatrix<VectorType>> &
382  get_matrix_diagonal() const;
383 
384 
390  void
391  precondition_Jacobi(VectorType & dst,
392  const VectorType &src,
393  const value_type omega) const;
394 
395  protected:
400  void
401  preprocess_constraints(VectorType &dst, const VectorType &src) const;
402 
407  void
408  postprocess_constraints(VectorType &dst, const VectorType &src) const;
409 
414  void
415  set_constrained_entries_to_one(VectorType &dst) const;
416 
420  virtual void
421  apply_add(VectorType &dst, const VectorType &src) const = 0;
422 
428  virtual void
429  Tapply_add(VectorType &dst, const VectorType &src) const;
430 
434  std::shared_ptr<const MatrixFree<dim, value_type, VectorizedArrayType>>
436 
441  std::shared_ptr<DiagonalMatrix<VectorType>> diagonal_entries;
442 
447  std::shared_ptr<DiagonalMatrix<VectorType>> inverse_diagonal_entries;
448 
453  std::vector<unsigned int> selected_rows;
454 
459  std::vector<unsigned int> selected_columns;
460 
461  private:
465  std::vector<std::vector<unsigned int>> edge_constrained_indices;
466 
470  mutable std::vector<std::vector<std::pair<value_type, value_type>>>
472 
478 
483  void
484  mult_add(VectorType & dst,
485  const VectorType &src,
486  const bool transpose) const;
487 
495  void
496  adjust_ghost_range_if_necessary(const VectorType &vec,
497  const bool is_row) const;
498  };
499 
500 
501 
538  template <typename OperatorType>
540  {
541  public:
545  using value_type = typename OperatorType::value_type;
546 
551 
556 
560  void
561  clear();
562 
566  void
567  initialize(const OperatorType &operator_in);
568 
572  template <typename VectorType>
573  void
574  vmult(VectorType &dst, const VectorType &src) const;
575 
579  template <typename VectorType>
580  void
581  Tvmult(VectorType &dst, const VectorType &src) const;
582 
586  template <typename VectorType>
587  void
588  initialize_dof_vector(VectorType &vec) const;
589 
590 
591  private:
596  };
597 
598 
599 
619  template <int dim,
620  int fe_degree,
621  int n_components = 1,
622  typename Number = double,
623  typename VectorizedArrayType = VectorizedArray<Number>>
625  {
626  static_assert(
628  "Type of Number and of VectorizedArrayType do not match.");
629 
630  public:
636  const FEEvaluationBase<dim,
637  n_components,
638  Number,
639  false,
640  VectorizedArrayType> &fe_eval);
641 
650  void
651  apply(const AlignedVector<VectorizedArrayType> &inverse_coefficient,
652  const unsigned int n_actual_components,
653  const VectorizedArrayType * in_array,
654  VectorizedArrayType * out_array) const;
655 
667  void
668  apply(const VectorizedArrayType *in_array,
669  VectorizedArrayType * out_array) const;
670 
704  void
705  transform_from_q_points_to_basis(const unsigned int n_actual_components,
706  const VectorizedArrayType *in_array,
707  VectorizedArrayType *out_array) const;
708 
714  void
715  fill_inverse_JxW_values(
716  AlignedVector<VectorizedArrayType> &inverse_jxw) const;
717 
718  private:
722  const FEEvaluationBase<dim,
723  n_components,
724  Number,
725  false,
726  VectorizedArrayType> &fe_eval;
727  };
728 
729 
730 
740  template <int dim,
741  int fe_degree,
742  int n_q_points_1d = fe_degree + 1,
743  int n_components = 1,
744  typename VectorType = LinearAlgebra::distributed::Vector<double>,
745  typename VectorizedArrayType =
747  class MassOperator : public Base<dim, VectorType, VectorizedArrayType>
748  {
749  public:
753  using value_type =
755 
759  using size_type =
761 
765  MassOperator();
766 
771  virtual void
772  compute_diagonal() override;
773 
774  private:
780  virtual void
781  apply_add(VectorType &dst, const VectorType &src) const override;
782 
786  void
787  local_apply_cell(
789  VectorType & dst,
790  const VectorType & src,
791  const std::pair<unsigned int, unsigned int> & cell_range) const;
792  };
793 
794 
795 
808  template <int dim,
809  int fe_degree,
810  int n_q_points_1d = fe_degree + 1,
811  int n_components = 1,
812  typename VectorType = LinearAlgebra::distributed::Vector<double>,
813  typename VectorizedArrayType =
815  class LaplaceOperator : public Base<dim, VectorType, VectorizedArrayType>
816  {
817  public:
821  using value_type =
823 
827  using size_type =
829 
833  LaplaceOperator();
834 
841  virtual void
842  compute_diagonal() override;
843 
894  void
895  set_coefficient(
896  const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient);
897 
902  virtual void
903  clear() override;
904 
911  std::shared_ptr<Table<2, VectorizedArrayType>>
912  get_coefficient();
913 
914  private:
920  virtual void
921  apply_add(VectorType &dst, const VectorType &src) const override;
922 
926  void
927  local_apply_cell(
929  VectorType & dst,
930  const VectorType & src,
931  const std::pair<unsigned int, unsigned int> & cell_range) const;
932 
936  void
937  local_diagonal_cell(
939  VectorType & dst,
940  const VectorType &,
941  const std::pair<unsigned int, unsigned int> &cell_range) const;
942 
946  void
947  do_operation_on_cell(
949  & phi,
950  const unsigned int cell) const;
951 
955  std::shared_ptr<Table<2, VectorizedArrayType>> scalar_coefficient;
956  };
957 
958 
959 
960  // ------------------------------------ inline functions ---------------------
961 
962  template <int dim,
963  int fe_degree,
964  int n_components,
965  typename Number,
966  typename VectorizedArrayType>
967  inline CellwiseInverseMassMatrix<dim,
968  fe_degree,
969  n_components,
970  Number,
971  VectorizedArrayType>::
972  CellwiseInverseMassMatrix(
973  const FEEvaluationBase<dim,
974  n_components,
975  Number,
976  false,
977  VectorizedArrayType> &fe_eval)
978  : fe_eval(fe_eval)
979  {}
980 
981 
982 
983  template <int dim,
984  int fe_degree,
985  int n_components,
986  typename Number,
987  typename VectorizedArrayType>
988  inline void
990  fe_degree,
991  n_components,
992  Number,
993  VectorizedArrayType>::
995  AlignedVector<VectorizedArrayType> &inverse_jxw) const
996  {
997  constexpr unsigned int dofs_per_component_on_cell =
998  Utilities::pow(fe_degree + 1, dim);
999  Assert(inverse_jxw.size() > 0 &&
1000  inverse_jxw.size() % dofs_per_component_on_cell == 0,
1001  ExcMessage(
1002  "Expected diagonal to be a multiple of scalar dof per cells"));
1003 
1004  // temporarily reduce size of inverse_jxw to dofs_per_cell to get JxW values
1005  // from fe_eval (will not reallocate any memory)
1006  for (unsigned int q = 0; q < dofs_per_component_on_cell; ++q)
1007  inverse_jxw[q] = 1. / fe_eval.JxW(q);
1008  // copy values to rest of vector
1009  for (unsigned int q = dofs_per_component_on_cell; q < inverse_jxw.size();)
1010  for (unsigned int i = 0; i < dofs_per_component_on_cell; ++i, ++q)
1011  inverse_jxw[q] = inverse_jxw[i];
1012  }
1013 
1014 
1015 
1016  template <int dim,
1017  int fe_degree,
1018  int n_components,
1019  typename Number,
1020  typename VectorizedArrayType>
1021  inline void
1023  dim,
1024  fe_degree,
1025  n_components,
1026  Number,
1027  VectorizedArrayType>::apply(const VectorizedArrayType *in_array,
1028  VectorizedArrayType * out_array) const
1029  {
1031  dim,
1032  fe_degree,
1033  n_components,
1034  VectorizedArrayType>::apply(fe_eval, in_array, out_array);
1035  }
1036 
1037 
1038 
1039  template <int dim,
1040  int fe_degree,
1041  int n_components,
1042  typename Number,
1043  typename VectorizedArrayType>
1044  inline void
1046  fe_degree,
1047  n_components,
1048  Number,
1049  VectorizedArrayType>::
1050  apply(const AlignedVector<VectorizedArrayType> &inverse_coefficients,
1051  const unsigned int n_actual_components,
1052  const VectorizedArrayType * in_array,
1053  VectorizedArrayType * out_array) const
1054  {
1056  fe_degree,
1057  n_components,
1058  VectorizedArrayType>::
1059  apply(fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
1060  inverse_coefficients,
1061  n_actual_components,
1062  in_array,
1063  out_array);
1064  }
1065 
1066 
1067 
1068  template <int dim,
1069  int fe_degree,
1070  int n_components,
1071  typename Number,
1072  typename VectorizedArrayType>
1073  inline void
1075  fe_degree,
1076  n_components,
1077  Number,
1078  VectorizedArrayType>::
1079  transform_from_q_points_to_basis(const unsigned int n_actual_components,
1080  const VectorizedArrayType *in_array,
1081  VectorizedArrayType * out_array) const
1082  {
1084  fe_degree,
1085  n_components,
1086  VectorizedArrayType>::
1088  fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
1089  n_actual_components,
1090  in_array,
1091  out_array);
1092  }
1093 
1094 
1095 
1096  //----------------- Base operator -----------------------------
1097  template <int dim, typename VectorType, typename VectorizedArrayType>
1099  : Subscriptor()
1100  , have_interface_matrices(false)
1101  {}
1102 
1103 
1104 
1105  template <int dim, typename VectorType, typename VectorizedArrayType>
1108  {
1109  Assert(data.get() != nullptr, ExcNotInitialized());
1111  0;
1112  for (const unsigned int selected_row : selected_rows)
1113  total_size += data->get_vector_partitioner(selected_row)->size();
1114  return total_size;
1115  }
1116 
1117 
1118 
1119  template <int dim, typename VectorType, typename VectorizedArrayType>
1122  {
1123  Assert(data.get() != nullptr, ExcNotInitialized());
1125  0;
1126  for (const unsigned int selected_column : selected_columns)
1127  total_size += data->get_vector_partitioner(selected_column)->size();
1128  return total_size;
1129  }
1130 
1131 
1132 
1133  template <int dim, typename VectorType, typename VectorizedArrayType>
1134  void
1136  {
1137  data.reset();
1138  inverse_diagonal_entries.reset();
1139  }
1140 
1141 
1142 
1143  template <int dim, typename VectorType, typename VectorizedArrayType>
1146  const unsigned int col) const
1147  {
1148  (void)col;
1149  Assert(row == col, ExcNotImplemented());
1150  Assert(inverse_diagonal_entries.get() != nullptr &&
1151  inverse_diagonal_entries->m() > 0,
1152  ExcNotInitialized());
1153  return 1.0 / (*inverse_diagonal_entries)(row, row);
1154  }
1155 
1156 
1157 
1158  template <int dim, typename VectorType, typename VectorizedArrayType>
1159  void
1161  VectorType &vec) const
1162  {
1163  Assert(data.get() != nullptr, ExcNotInitialized());
1165  for (unsigned int i = 0; i < BlockHelper::n_blocks(vec); ++i)
1166  {
1167  const unsigned int index = selected_rows[i];
1168  if (!BlockHelper::subblock(vec, index)
1169  .partitioners_are_compatible(
1170  *data->get_dof_info(index).vector_partitioner))
1171  data->initialize_dof_vector(BlockHelper::subblock(vec, index), index);
1172 
1173  Assert(BlockHelper::subblock(vec, index)
1174  .partitioners_are_globally_compatible(
1175  *data->get_dof_info(index).vector_partitioner),
1176  ExcInternalError());
1177  }
1179  }
1180 
1181 
1182 
1183  template <int dim, typename VectorType, typename VectorizedArrayType>
1184  void
1187  data_,
1188  const std::vector<unsigned int> &given_row_selection,
1189  const std::vector<unsigned int> &given_column_selection)
1190  {
1191  data = data_;
1192 
1193  selected_rows.clear();
1194  selected_columns.clear();
1195  if (given_row_selection.empty())
1196  for (unsigned int i = 0; i < data_->n_components(); ++i)
1197  selected_rows.push_back(i);
1198  else
1199  {
1200  for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1201  {
1202  AssertIndexRange(given_row_selection[i], data_->n_components());
1203  for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1204  if (j != i)
1205  Assert(given_row_selection[j] != given_row_selection[i],
1206  ExcMessage("Given row indices must be unique"));
1207 
1208  selected_rows.push_back(given_row_selection[i]);
1209  }
1210  }
1211  if (given_column_selection.size() == 0)
1213  else
1214  {
1215  for (unsigned int i = 0; i < given_column_selection.size(); ++i)
1216  {
1217  AssertIndexRange(given_column_selection[i], data_->n_components());
1218  for (unsigned int j = 0; j < given_column_selection.size(); ++j)
1219  if (j != i)
1220  Assert(given_column_selection[j] != given_column_selection[i],
1221  ExcMessage("Given column indices must be unique"));
1222 
1223  selected_columns.push_back(given_column_selection[i]);
1224  }
1225  }
1226 
1227  edge_constrained_indices.clear();
1228  edge_constrained_indices.resize(selected_rows.size());
1229  edge_constrained_values.clear();
1230  edge_constrained_values.resize(selected_rows.size());
1231  have_interface_matrices = false;
1232  }
1233 
1234 
1235 
1236  template <int dim, typename VectorType, typename VectorizedArrayType>
1237  void
1240  data_,
1241  const MGConstrainedDoFs & mg_constrained_dofs,
1242  const unsigned int level,
1243  const std::vector<unsigned int> &given_row_selection)
1244  {
1245  std::vector<MGConstrainedDoFs> mg_constrained_dofs_vector(
1246  1, mg_constrained_dofs);
1247  initialize(data_, mg_constrained_dofs_vector, level, given_row_selection);
1248  }
1249 
1250 
1251 
1252  template <int dim, typename VectorType, typename VectorizedArrayType>
1253  void
1256  data_,
1257  const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
1258  const unsigned int level,
1259  const std::vector<unsigned int> & given_row_selection)
1260  {
1262  ExcMessage("level is not set"));
1263 
1264  selected_rows.clear();
1265  selected_columns.clear();
1266  if (given_row_selection.empty())
1267  for (unsigned int i = 0; i < data_->n_components(); ++i)
1268  selected_rows.push_back(i);
1269  else
1270  {
1271  for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1272  {
1273  AssertIndexRange(given_row_selection[i], data_->n_components());
1274  for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1275  if (j != i)
1276  Assert(given_row_selection[j] != given_row_selection[i],
1277  ExcMessage("Given row indices must be unique"));
1278 
1279  selected_rows.push_back(given_row_selection[i]);
1280  }
1281  }
1283 
1284  AssertDimension(mg_constrained_dofs.size(), selected_rows.size());
1285  edge_constrained_indices.clear();
1286  edge_constrained_indices.resize(selected_rows.size());
1287  edge_constrained_values.clear();
1288  edge_constrained_values.resize(selected_rows.size());
1289 
1290  data = data_;
1291 
1292  for (unsigned int j = 0; j < selected_rows.size(); ++j)
1293  {
1294  if (data_->n_macro_cells() > 0)
1295  {
1296  AssertDimension(level, data_->get_cell_iterator(0, 0, j)->level());
1297  }
1298 
1299  // setup edge_constrained indices
1300  std::vector<types::global_dof_index> interface_indices;
1301  mg_constrained_dofs[j]
1302  .get_refinement_edge_indices(level)
1303  .fill_index_vector(interface_indices);
1304  edge_constrained_indices[j].clear();
1305  edge_constrained_indices[j].reserve(interface_indices.size());
1306  edge_constrained_values[j].resize(interface_indices.size());
1307  const IndexSet &locally_owned =
1308  data->get_dof_handler(selected_rows[j]).locally_owned_mg_dofs(level);
1309  for (const auto interface_index : interface_indices)
1310  if (locally_owned.is_element(interface_index))
1311  edge_constrained_indices[j].push_back(
1312  locally_owned.index_within_set(interface_index));
1315  static_cast<unsigned int>(edge_constrained_indices[j].size()),
1316  data->get_vector_partitioner()->get_mpi_communicator()) > 0;
1317  }
1318  }
1319 
1320 
1321 
1322  template <int dim, typename VectorType, typename VectorizedArrayType>
1323  void
1325  VectorType &dst) const
1326  {
1327  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1328  {
1329  const std::vector<unsigned int> &constrained_dofs =
1330  data->get_constrained_dofs(selected_rows[j]);
1331  for (const auto constrained_dof : constrained_dofs)
1332  BlockHelper::subblock(dst, j).local_element(constrained_dof) = 1.;
1333  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1334  BlockHelper::subblock(dst, j).local_element(
1335  edge_constrained_indices[j][i]) = 1.;
1336  }
1337  }
1338 
1339 
1340 
1341  template <int dim, typename VectorType, typename VectorizedArrayType>
1342  void
1344  const VectorType &src) const
1345  {
1346  using Number =
1348  dst = Number(0.);
1349  vmult_add(dst, src);
1350  }
1351 
1352 
1353 
1354  template <int dim, typename VectorType, typename VectorizedArrayType>
1355  void
1357  VectorType & dst,
1358  const VectorType &src) const
1359  {
1360  mult_add(dst, src, false);
1361  }
1362 
1363 
1364 
1365  template <int dim, typename VectorType, typename VectorizedArrayType>
1366  void
1368  VectorType & dst,
1369  const VectorType &src) const
1370  {
1371  mult_add(dst, src, true);
1372  }
1373 
1374 
1375 
1376  template <int dim, typename VectorType, typename VectorizedArrayType>
1377  void
1379  const VectorType &src,
1380  const bool is_row) const
1381  {
1382  using Number =
1384  for (unsigned int i = 0; i < BlockHelper::n_blocks(src); ++i)
1385  {
1386  const unsigned int mf_component =
1387  is_row ? selected_rows[i] : selected_columns[i];
1388  // If both vectors use the same partitioner -> done
1389  if (BlockHelper::subblock(src, i).get_partitioner().get() ==
1390  data->get_dof_info(mf_component).vector_partitioner.get())
1391  continue;
1392 
1393  // If not, assert that the local ranges are the same and reset to the
1394  // current partitioner
1395  Assert(
1396  BlockHelper::subblock(src, i).get_partitioner()->local_size() ==
1397  data->get_dof_info(mf_component).vector_partitioner->local_size(),
1398  ExcMessage("The vector passed to the vmult() function does not have "
1399  "the correct size for compatibility with MatrixFree."));
1400 
1401  // copy the vector content to a temporary vector so that it does not get
1402  // lost
1404  BlockHelper::subblock(src, i));
1405  BlockHelper::subblock(const_cast<VectorType &>(src), i)
1406  .reinit(data->get_dof_info(mf_component).vector_partitioner);
1407  BlockHelper::subblock(const_cast<VectorType &>(src), i)
1408  .copy_locally_owned_data_from(copy_vec);
1409  }
1410  }
1411 
1412 
1413 
1414  template <int dim, typename VectorType, typename VectorizedArrayType>
1415  void
1417  VectorType & dst,
1418  const VectorType &src) const
1419  {
1420  using Number =
1422  adjust_ghost_range_if_necessary(src, false);
1424 
1425  // set zero Dirichlet values on the input vector (and remember the src and
1426  // dst values because we need to reset them at the end)
1427  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1428  {
1429  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1430  {
1431  edge_constrained_values[j][i] = std::pair<Number, Number>(
1432  BlockHelper::subblock(src, j).local_element(
1433  edge_constrained_indices[j][i]),
1434  BlockHelper::subblock(dst, j).local_element(
1435  edge_constrained_indices[j][i]));
1436  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1437  .local_element(edge_constrained_indices[j][i]) = 0.;
1438  }
1439  }
1440  }
1441 
1442 
1443 
1444  template <int dim, typename VectorType, typename VectorizedArrayType>
1445  void
1447  VectorType & dst,
1448  const VectorType &src,
1449  const bool transpose) const
1450  {
1451  AssertDimension(dst.size(), src.size());
1454  preprocess_constraints(dst, src);
1455  if (transpose)
1456  Tapply_add(dst, src);
1457  else
1458  apply_add(dst, src);
1459  postprocess_constraints(dst, src);
1460  }
1461 
1462 
1463 
1464  template <int dim, typename VectorType, typename VectorizedArrayType>
1465  void
1467  VectorType & dst,
1468  const VectorType &src) const
1469  {
1470  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1471  {
1472  const std::vector<unsigned int> &constrained_dofs =
1473  data->get_constrained_dofs(selected_rows[j]);
1474  for (const auto constrained_dof : constrained_dofs)
1475  BlockHelper::subblock(dst, j).local_element(constrained_dof) +=
1476  BlockHelper::subblock(src, j).local_element(constrained_dof);
1477  }
1478 
1479  // reset edge constrained values, multiply by unit matrix and add into
1480  // destination
1481  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1482  {
1483  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1484  {
1485  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1486  .local_element(edge_constrained_indices[j][i]) =
1487  edge_constrained_values[j][i].first;
1488  BlockHelper::subblock(dst, j).local_element(
1489  edge_constrained_indices[j][i]) =
1490  edge_constrained_values[j][i].second +
1491  edge_constrained_values[j][i].first;
1492  }
1493  }
1494  }
1495 
1496 
1497 
1498  template <int dim, typename VectorType, typename VectorizedArrayType>
1499  void
1501  VectorType & dst,
1502  const VectorType &src) const
1503  {
1504  using Number =
1506  AssertDimension(dst.size(), src.size());
1507  adjust_ghost_range_if_necessary(src, false);
1509 
1510  dst = Number(0.);
1511 
1513  return;
1514 
1515  // set zero Dirichlet values on the input vector (and remember the src and
1516  // dst values because we need to reset them at the end)
1517  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1518  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1519  {
1520  edge_constrained_values[j][i] = std::pair<Number, Number>(
1521  BlockHelper::subblock(src, j).local_element(
1522  edge_constrained_indices[j][i]),
1523  BlockHelper::subblock(dst, j).local_element(
1524  edge_constrained_indices[j][i]));
1525  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1526  .local_element(edge_constrained_indices[j][i]) = 0.;
1527  }
1528 
1529  apply_add(dst, src);
1530 
1531  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1532  {
1533  unsigned int c = 0;
1534  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1535  {
1536  for (; c < edge_constrained_indices[j][i]; ++c)
1537  BlockHelper::subblock(dst, j).local_element(c) = 0.;
1538  ++c;
1539 
1540  // reset the src values
1541  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1542  .local_element(edge_constrained_indices[j][i]) =
1543  edge_constrained_values[j][i].first;
1544  }
1545  for (; c < BlockHelper::subblock(dst, j).local_size(); ++c)
1546  BlockHelper::subblock(dst, j).local_element(c) = 0.;
1547  }
1548  }
1549 
1550 
1551 
1552  template <int dim, typename VectorType, typename VectorizedArrayType>
1553  void
1555  VectorType & dst,
1556  const VectorType &src) const
1557  {
1558  using Number =
1560  AssertDimension(dst.size(), src.size());
1561  adjust_ghost_range_if_necessary(src, false);
1563 
1564  dst = Number(0.);
1565 
1567  return;
1568 
1569  VectorType src_cpy(src);
1570  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1571  {
1572  unsigned int c = 0;
1573  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1574  {
1575  for (; c < edge_constrained_indices[j][i]; ++c)
1576  BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1577  ++c;
1578  }
1579  for (; c < BlockHelper::subblock(src_cpy, j).local_size(); ++c)
1580  BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1581  }
1582 
1583  apply_add(dst, src_cpy);
1584 
1585  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1586  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1587  BlockHelper::subblock(dst, j).local_element(
1588  edge_constrained_indices[j][i]) = 0.;
1589  }
1590 
1591 
1592 
1593  template <int dim, typename VectorType, typename VectorizedArrayType>
1594  void
1596  VectorType & dst,
1597  const VectorType &src) const
1598  {
1599  using Number =
1601  dst = Number(0.);
1602  Tvmult_add(dst, src);
1603  }
1604 
1605 
1606 
1607  template <int dim, typename VectorType, typename VectorizedArrayType>
1608  std::size_t
1610  {
1611  return inverse_diagonal_entries.get() != nullptr ?
1612  inverse_diagonal_entries->memory_consumption() :
1613  sizeof(*this);
1614  }
1615 
1616 
1617 
1618  template <int dim, typename VectorType, typename VectorizedArrayType>
1619  std::shared_ptr<const MatrixFree<
1620  dim,
1622  VectorizedArrayType>>
1624  {
1625  return data;
1626  }
1627 
1628 
1629 
1630  template <int dim, typename VectorType, typename VectorizedArrayType>
1631  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1633  const
1634  {
1635  Assert(inverse_diagonal_entries.get() != nullptr &&
1636  inverse_diagonal_entries->m() > 0,
1637  ExcNotInitialized());
1638  return inverse_diagonal_entries;
1639  }
1640 
1641 
1642 
1643  template <int dim, typename VectorType, typename VectorizedArrayType>
1644  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1646  {
1647  Assert(diagonal_entries.get() != nullptr && diagonal_entries->m() > 0,
1648  ExcNotInitialized());
1649  return diagonal_entries;
1650  }
1651 
1652 
1653 
1654  template <int dim, typename VectorType, typename VectorizedArrayType>
1655  void
1657  VectorType & dst,
1658  const VectorType &src) const
1659  {
1660  apply_add(dst, src);
1661  }
1662 
1663 
1664 
1665  template <int dim, typename VectorType, typename VectorizedArrayType>
1666  void
1668  VectorType & dst,
1669  const VectorType & src,
1670  const typename Base<dim, VectorType, VectorizedArrayType>::value_type omega)
1671  const
1672  {
1674  ExcNotInitialized());
1675  inverse_diagonal_entries->vmult(dst, src);
1676  dst *= omega;
1677  }
1678 
1679 
1680 
1681  //------------------------- MGInterfaceOperator ------------------------------
1682 
1683  template <typename OperatorType>
1685  : Subscriptor()
1686  , mf_base_operator(nullptr)
1687  {}
1688 
1689 
1690 
1691  template <typename OperatorType>
1692  void
1694  {
1695  mf_base_operator = nullptr;
1696  }
1697 
1698 
1699 
1700  template <typename OperatorType>
1701  void
1702  MGInterfaceOperator<OperatorType>::initialize(const OperatorType &operator_in)
1703  {
1704  mf_base_operator = &operator_in;
1705  }
1706 
1707 
1708 
1709  template <typename OperatorType>
1710  template <typename VectorType>
1711  void
1713  const VectorType &src) const
1714  {
1715 #ifndef DEAL_II_MSVC
1716  static_assert(
1718  "The vector type must be based on the same value type as this "
1719  "operator");
1720 #endif
1721 
1722  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1723 
1724  mf_base_operator->vmult_interface_down(dst, src);
1725  }
1726 
1727 
1728 
1729  template <typename OperatorType>
1730  template <typename VectorType>
1731  void
1733  const VectorType &src) const
1734  {
1735 #ifndef DEAL_II_MSVC
1736  static_assert(
1738  "The vector type must be based on the same value type as this "
1739  "operator");
1740 #endif
1741 
1742  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1743 
1744  mf_base_operator->vmult_interface_up(dst, src);
1745  }
1746 
1747 
1748 
1749  template <typename OperatorType>
1750  template <typename VectorType>
1751  void
1753  VectorType &vec) const
1754  {
1755  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1756 
1757  mf_base_operator->initialize_dof_vector(vec);
1758  }
1759 
1760 
1761 
1762  //-----------------------------MassOperator----------------------------------
1763 
1764  template <int dim,
1765  int fe_degree,
1766  int n_q_points_1d,
1767  int n_components,
1768  typename VectorType,
1769  typename VectorizedArrayType>
1770  MassOperator<dim,
1771  fe_degree,
1772  n_q_points_1d,
1773  n_components,
1774  VectorType,
1775  VectorizedArrayType>::MassOperator()
1776  : Base<dim, VectorType, VectorizedArrayType>()
1777  {}
1778 
1779 
1780 
1781  template <int dim,
1782  int fe_degree,
1783  int n_q_points_1d,
1784  int n_components,
1785  typename VectorType,
1786  typename VectorizedArrayType>
1787  void
1788  MassOperator<dim,
1789  fe_degree,
1790  n_q_points_1d,
1791  n_components,
1792  VectorType,
1793  VectorizedArrayType>::compute_diagonal()
1794  {
1795  using Number =
1796  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1798  ExcNotInitialized());
1799 
1800  this->inverse_diagonal_entries =
1801  std::make_shared<DiagonalMatrix<VectorType>>();
1802  this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1803  VectorType &inverse_diagonal_vector =
1804  this->inverse_diagonal_entries->get_vector();
1805  VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1806  this->initialize_dof_vector(inverse_diagonal_vector);
1807  this->initialize_dof_vector(diagonal_vector);
1808  inverse_diagonal_vector = Number(1.);
1809  apply_add(diagonal_vector, inverse_diagonal_vector);
1810 
1811  this->set_constrained_entries_to_one(diagonal_vector);
1812  inverse_diagonal_vector = diagonal_vector;
1813 
1814  const unsigned int local_size = inverse_diagonal_vector.local_size();
1815  for (unsigned int i = 0; i < local_size; ++i)
1816  inverse_diagonal_vector.local_element(i) =
1817  Number(1.) / inverse_diagonal_vector.local_element(i);
1818 
1819  inverse_diagonal_vector.update_ghost_values();
1820  diagonal_vector.update_ghost_values();
1821  }
1822 
1823 
1824 
1825  template <int dim,
1826  int fe_degree,
1827  int n_q_points_1d,
1828  int n_components,
1829  typename VectorType,
1830  typename VectorizedArrayType>
1831  void
1832  MassOperator<dim,
1833  fe_degree,
1834  n_q_points_1d,
1835  n_components,
1836  VectorType,
1837  VectorizedArrayType>::apply_add(VectorType & dst,
1838  const VectorType &src) const
1839  {
1841  &MassOperator::local_apply_cell, this, dst, src);
1842  }
1843 
1844 
1845 
1846  template <int dim,
1847  int fe_degree,
1848  int n_q_points_1d,
1849  int n_components,
1850  typename VectorType,
1851  typename VectorizedArrayType>
1852  void
1853  MassOperator<dim,
1854  fe_degree,
1855  n_q_points_1d,
1856  n_components,
1857  VectorType,
1858  VectorizedArrayType>::
1860  const MatrixFree<
1861  dim,
1862  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
1863  VectorizedArrayType> & data,
1864  VectorType & dst,
1865  const VectorType & src,
1866  const std::pair<unsigned int, unsigned int> &cell_range) const
1867  {
1868  using Number =
1869  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1870  FEEvaluation<dim,
1871  fe_degree,
1872  n_q_points_1d,
1873  n_components,
1874  Number,
1875  VectorizedArrayType>
1876  phi(data, this->selected_rows[0]);
1877  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
1878  {
1879  phi.reinit(cell);
1880  phi.read_dof_values(src);
1881  phi.evaluate(true, false, false);
1882  for (unsigned int q = 0; q < phi.n_q_points; ++q)
1883  phi.submit_value(phi.get_value(q), q);
1884  phi.integrate(true, false);
1885  phi.distribute_local_to_global(dst);
1886  }
1887  }
1888 
1889 
1890  //-----------------------------LaplaceOperator----------------------------------
1891 
1892  template <int dim,
1893  int fe_degree,
1894  int n_q_points_1d,
1895  int n_components,
1896  typename VectorType,
1897  typename VectorizedArrayType>
1898  LaplaceOperator<dim,
1899  fe_degree,
1900  n_q_points_1d,
1901  n_components,
1902  VectorType,
1903  VectorizedArrayType>::LaplaceOperator()
1904  : Base<dim, VectorType, VectorizedArrayType>()
1905  {}
1906 
1907 
1908 
1909  template <int dim,
1910  int fe_degree,
1911  int n_q_points_1d,
1912  int n_components,
1913  typename VectorType,
1914  typename VectorizedArrayType>
1915  void
1916  LaplaceOperator<dim,
1917  fe_degree,
1918  n_q_points_1d,
1919  n_components,
1920  VectorType,
1921  VectorizedArrayType>::clear()
1922  {
1924  scalar_coefficient.reset();
1925  }
1926 
1927 
1928 
1929  template <int dim,
1930  int fe_degree,
1931  int n_q_points_1d,
1932  int n_components,
1933  typename VectorType,
1934  typename VectorizedArrayType>
1935  void
1936  LaplaceOperator<dim,
1937  fe_degree,
1938  n_q_points_1d,
1939  n_components,
1940  VectorType,
1941  VectorizedArrayType>::
1943  const std::shared_ptr<Table<2, VectorizedArrayType>> &scalar_coefficient_)
1944  {
1945  scalar_coefficient = scalar_coefficient_;
1946  }
1947 
1948 
1949 
1950  template <int dim,
1951  int fe_degree,
1952  int n_q_points_1d,
1953  int n_components,
1954  typename VectorType,
1955  typename VectorizedArrayType>
1956  std::shared_ptr<Table<2, VectorizedArrayType>>
1957  LaplaceOperator<dim,
1958  fe_degree,
1959  n_q_points_1d,
1960  n_components,
1961  VectorType,
1962  VectorizedArrayType>::get_coefficient()
1963  {
1965  return scalar_coefficient;
1966  }
1967 
1968 
1969 
1970  template <int dim,
1971  int fe_degree,
1972  int n_q_points_1d,
1973  int n_components,
1974  typename VectorType,
1975  typename VectorizedArrayType>
1976  void
1977  LaplaceOperator<dim,
1978  fe_degree,
1979  n_q_points_1d,
1980  n_components,
1981  VectorType,
1982  VectorizedArrayType>::compute_diagonal()
1983  {
1984  using Number =
1985  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
1987  ExcNotInitialized());
1988 
1989  this->inverse_diagonal_entries =
1990  std::make_shared<DiagonalMatrix<VectorType>>();
1991  this->diagonal_entries = std::make_shared<DiagonalMatrix<VectorType>>();
1992  VectorType &inverse_diagonal_vector =
1993  this->inverse_diagonal_entries->get_vector();
1994  VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1995  this->initialize_dof_vector(inverse_diagonal_vector);
1996  this->initialize_dof_vector(diagonal_vector);
1997 
1998  this->data->cell_loop(&LaplaceOperator::local_diagonal_cell,
1999  this,
2000  diagonal_vector,
2001  /*unused*/ diagonal_vector);
2002  this->set_constrained_entries_to_one(diagonal_vector);
2003 
2004  inverse_diagonal_vector = diagonal_vector;
2005 
2006  for (unsigned int i = 0; i < inverse_diagonal_vector.local_size(); ++i)
2007  if (std::abs(inverse_diagonal_vector.local_element(i)) >
2009  inverse_diagonal_vector.local_element(i) =
2010  1. / inverse_diagonal_vector.local_element(i);
2011  else
2012  inverse_diagonal_vector.local_element(i) = 1.;
2013 
2014  inverse_diagonal_vector.update_ghost_values();
2015  diagonal_vector.update_ghost_values();
2016  }
2017 
2018 
2019 
2020  template <int dim,
2021  int fe_degree,
2022  int n_q_points_1d,
2023  int n_components,
2024  typename VectorType,
2025  typename VectorizedArrayType>
2026  void
2027  LaplaceOperator<dim,
2028  fe_degree,
2029  n_q_points_1d,
2030  n_components,
2031  VectorType,
2032  VectorizedArrayType>::apply_add(VectorType & dst,
2033  const VectorType &src) const
2034  {
2036  &LaplaceOperator::local_apply_cell, this, dst, src);
2037  }
2038 
2039  namespace Implementation
2040  {
2041  template <typename VectorizedArrayType>
2042  bool
2043  non_negative(const VectorizedArrayType &n)
2044  {
2045  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2046  if (n[v] < 0.)
2047  return false;
2048 
2049  return true;
2050  }
2051  } // namespace Implementation
2052 
2053 
2054 
2055  template <int dim,
2056  int fe_degree,
2057  int n_q_points_1d,
2058  int n_components,
2059  typename VectorType,
2060  typename VectorizedArrayType>
2061  void
2062  LaplaceOperator<dim,
2063  fe_degree,
2064  n_q_points_1d,
2065  n_components,
2066  VectorType,
2067  VectorizedArrayType>::
2069  FEEvaluation<
2070  dim,
2071  fe_degree,
2072  n_q_points_1d,
2073  n_components,
2074  typename Base<dim, VectorType, VectorizedArrayType>::value_type> &phi,
2075  const unsigned int cell) const
2076  {
2077  phi.evaluate(false, true, false);
2078  if (scalar_coefficient.get())
2079  {
2080  Assert(scalar_coefficient->size(1) == 1 ||
2081  scalar_coefficient->size(1) == phi.n_q_points,
2082  ExcMessage("The number of columns in the coefficient table must "
2083  "be either 1 or the number of quadrature points " +
2084  std::to_string(phi.n_q_points) +
2085  ", but the given value was " +
2086  std::to_string(scalar_coefficient->size(1))));
2087  if (scalar_coefficient->size(1) == phi.n_q_points)
2088  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2089  {
2091  (*scalar_coefficient)(cell, q)),
2092  ExcMessage("Coefficient must be non-negative"));
2093  phi.submit_gradient((*scalar_coefficient)(cell, q) *
2094  phi.get_gradient(q),
2095  q);
2096  }
2097  else
2098  {
2100  ExcMessage("Coefficient must be non-negative"));
2101  const VectorizedArrayType coefficient =
2102  (*scalar_coefficient)(cell, 0);
2103  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2104  phi.submit_gradient(coefficient * phi.get_gradient(q), q);
2105  }
2106  }
2107  else
2108  {
2109  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2110  {
2111  phi.submit_gradient(phi.get_gradient(q), q);
2112  }
2113  }
2114  phi.integrate(false, true);
2115  }
2116 
2117 
2118 
2119  template <int dim,
2120  int fe_degree,
2121  int n_q_points_1d,
2122  int n_components,
2123  typename VectorType,
2124  typename VectorizedArrayType>
2125  void
2126  LaplaceOperator<dim,
2127  fe_degree,
2128  n_q_points_1d,
2129  n_components,
2130  VectorType,
2131  VectorizedArrayType>::
2133  const MatrixFree<
2134  dim,
2135  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2136  VectorizedArrayType> & data,
2137  VectorType & dst,
2138  const VectorType & src,
2139  const std::pair<unsigned int, unsigned int> &cell_range) const
2140  {
2141  using Number =
2142  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2144  data, this->selected_rows[0]);
2145  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2146  {
2147  phi.reinit(cell);
2148  phi.read_dof_values(src);
2149  do_operation_on_cell(phi, cell);
2150  phi.distribute_local_to_global(dst);
2151  }
2152  }
2153 
2154 
2155  template <int dim,
2156  int fe_degree,
2157  int n_q_points_1d,
2158  int n_components,
2159  typename VectorType,
2160  typename VectorizedArrayType>
2161  void
2162  LaplaceOperator<dim,
2163  fe_degree,
2164  n_q_points_1d,
2165  n_components,
2166  VectorType,
2167  VectorizedArrayType>::
2169  const MatrixFree<
2170  dim,
2171  typename Base<dim, VectorType, VectorizedArrayType>::value_type,
2172  VectorizedArrayType> &data,
2173  VectorType & dst,
2174  const VectorType &,
2175  const std::pair<unsigned int, unsigned int> &cell_range) const
2176  {
2177  using Number =
2178  typename Base<dim, VectorType, VectorizedArrayType>::value_type;
2179 
2181  data, this->selected_rows[0]);
2182  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2183  {
2184  phi.reinit(cell);
2185  VectorizedArrayType local_diagonal_vector[phi.static_dofs_per_cell];
2186  for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2187  {
2188  for (unsigned int j = 0; j < phi.dofs_per_component; ++j)
2189  phi.begin_dof_values()[j] = VectorizedArrayType();
2190  phi.begin_dof_values()[i] = 1.;
2191  do_operation_on_cell(phi, cell);
2192  local_diagonal_vector[i] = phi.begin_dof_values()[i];
2193  }
2194  for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2195  for (unsigned int c = 0; c < phi.n_components; ++c)
2196  phi.begin_dof_values()[i + c * phi.dofs_per_component] =
2197  local_diagonal_vector[i];
2198  phi.distribute_local_to_global(dst);
2199  }
2200  }
2201 
2202 
2203 } // end of namespace MatrixFreeOperators
2204 
2205 
2207 
2208 #endif
typename OperatorType::value_type value_type
Definition: operators.h:545
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:994
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:1632
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
virtual std::size_t memory_consumption() const
Definition: operators.h:1609
const FEEvaluationBase< dim, n_components, Number, false, VectorizedArrayType > & fe_eval
Definition: operators.h:726
types::global_dof_index size_type
Definition: cuda_kernels.h:45
size_type n() const
Definition: operators.h:1121
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:2168
void precondition_Jacobi(VectorType &dst, const VectorType &src, const value_type omega) const
Definition: operators.h:1667
std::shared_ptr< DiagonalMatrix< VectorType > > diagonal_entries
Definition: operators.h:441
void read_dof_values(const VectorType &src, const unsigned int first_index=0)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
std::vector< std::vector< std::pair< value_type, value_type > > > edge_constrained_values
Definition: operators.h:471
static ::ExceptionBase & ExcNotInitialized()
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > data
Definition: operators.h:435
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1595
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
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:197
typename OperatorType::size_type size_type
Definition: operators.h:550
virtual void compute_diagonal() override
Definition: operators.h:1982
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:476
std::shared_ptr< const MatrixFree< dim, value_type, VectorizedArrayType > > get_matrix_free() const
Definition: operators.h:1623
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1752
void vmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1356
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal() const
Definition: operators.h:1645
static ::ExceptionBase & ExcMessage(std::string arg1)
void vmult_interface_up(VectorType &dst, const VectorType &src) const
Definition: operators.h:1554
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:1419
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:2132
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1922
void postprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1466
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:1656
std::shared_ptr< Table< 2, VectorizedArrayType > > scalar_coefficient
Definition: operators.h:955
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
auto apply(F &&fn, Tuple &&t) -> decltype(apply_impl(std::forward< F >(fn), std::forward< Tuple >(t), std_cxx14::make_index_sequence< std::tuple_size< typename std::remove_reference< Tuple >::type >::value >()))
Definition: tuple.h:40
unsigned int level
Definition: grid_out.cc:4355
const unsigned int dofs_per_component
std::vector< unsigned int > selected_rows
Definition: operators.h:453
std::string to_string(const T &t)
Definition: patterns.h:2360
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:1135
void set_constrained_entries_to_one(VectorType &dst) const
Definition: operators.h:1324
virtual void clear() override
Definition: operators.h:1921
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1160
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1712
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition: operators.h:2032
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition: operators.h:1837
void Tvmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1367
static constexpr unsigned int n_components
size_type m() const
Definition: operators.h:1107
static constexpr unsigned int static_dofs_per_cell
typename VectorType::size_type size_type
Definition: operators.h:202
void mult_add(VectorType &dst, const VectorType &src, const bool transpose) const
Definition: operators.h:1446
std::shared_ptr< Table< 2, VectorizedArrayType > > get_coefficient()
Definition: operators.h:1962
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:2068
const VectorizedArrayType * begin_dof_values() const
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
std::vector< unsigned int > selected_columns
Definition: operators.h:459
void vmult_interface_down(VectorType &dst, const VectorType &src) const
Definition: operators.h:1500
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:1145
Definition: table.h:699
std::vector< UnivariateShapeData< Number > > data
Definition: shape_info.h:355
bool is_element(const size_type index) const
Definition: index_set.h:1766
void initialize(const OperatorType &operator_in)
Definition: operators.h:1702
SmartPointer< const OperatorType > mf_base_operator
Definition: operators.h:595
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1732
void apply(const AlignedVector< VectorizedArrayType > &inverse_coefficient, const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition: operators.h:1050
std::shared_ptr< DiagonalMatrix< VectorType > > inverse_diagonal_entries
Definition: operators.h:447
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1343
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
bool non_negative(const VectorizedArrayType &n)
Definition: operators.h:2043
static const bool value
virtual void compute_diagonal() override
Definition: operators.h:1793
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:1859
std::vector< std::vector< unsigned int > > edge_constrained_indices
Definition: operators.h:465
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:1185
void preprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1416
void transform_from_q_points_to_basis(const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition: operators.h:1079
void set_coefficient(const std::shared_ptr< Table< 2, VectorizedArrayType >> &scalar_coefficient)
Definition: operators.h:1942
static ::ExceptionBase & ExcInternalError()
void adjust_ghost_range_if_necessary(const VectorType &vec, const bool is_row) const
Definition: operators.h:1378