Reference documentation for deal.II version Git cb0bd54b52 2019-09-21 16:31:22 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
operators.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2011 - 2019 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/exceptions.h>
22 #include <deal.II/base/subscriptor.h>
23 #include <deal.II/base/vectorization.h>
24 
25 #include <deal.II/lac/diagonal_matrix.h>
26 #include <deal.II/lac/la_parallel_vector.h>
27 
28 #include <deal.II/matrix_free/fe_evaluation.h>
29 #include <deal.II/matrix_free/matrix_free.h>
30 
31 #include <deal.II/multigrid/mg_constrained_dofs.h>
32 
33 
34 DEAL_II_NAMESPACE_OPEN
35 
36 
37 namespace MatrixFreeOperators
38 {
39  namespace BlockHelper
40  {
41  // workaroud for unifying non-block vector and block vector implementations
42  // a non-block vector has one block and the only subblock is the vector
43  // itself
44  template <typename VectorType>
45  typename std::enable_if<IsBlockVector<VectorType>::value,
46  unsigned int>::type
47  n_blocks(const VectorType &vector)
48  {
49  return vector.n_blocks();
50  }
51 
52  template <typename VectorType>
53  typename std::enable_if<!IsBlockVector<VectorType>::value,
54  unsigned int>::type
55  n_blocks(const VectorType &)
56  {
57  return 1;
58  }
59 
60  template <typename VectorType>
61  typename std::enable_if<IsBlockVector<VectorType>::value,
62  typename VectorType::BlockType &>::type
63  subblock(VectorType &vector, unsigned int block_no)
64  {
65  return vector.block(block_no);
66  }
67 
68  template <typename VectorType>
69  typename std::enable_if<IsBlockVector<VectorType>::value,
70  const typename VectorType::BlockType &>::type
71  subblock(const VectorType &vector, unsigned int block_no)
72  {
73  AssertIndexRange(block_no, vector.n_blocks());
74  return vector.block(block_no);
75  }
76 
77  template <typename VectorType>
78  typename std::enable_if<!IsBlockVector<VectorType>::value,
79  VectorType &>::type
80  subblock(VectorType &vector, unsigned int)
81  {
82  return vector;
83  }
84 
85  template <typename VectorType>
86  typename std::enable_if<!IsBlockVector<VectorType>::value,
87  const VectorType &>::type
88  subblock(const VectorType &vector, unsigned int)
89  {
90  return vector;
91  }
92 
93  template <typename VectorType>
94  typename std::enable_if<IsBlockVector<VectorType>::value, void>::type
95  collect_sizes(VectorType &vector)
96  {
97  vector.collect_sizes();
98  }
99 
100  template <typename VectorType>
101  typename std::enable_if<!IsBlockVector<VectorType>::value, void>::type
102  collect_sizes(const VectorType &)
103  {}
104  } // namespace BlockHelper
105 
185  template <int dim,
186  typename VectorType = LinearAlgebra::distributed::Vector<double>>
187  class Base : public Subscriptor
188  {
189  public:
193  using value_type = typename VectorType::value_type;
194 
198  using size_type = typename VectorType::size_type;
199 
203  Base();
204 
208  virtual ~Base() override = default;
209 
214  virtual void
215  clear();
216 
233  void
234  initialize(std::shared_ptr<const MatrixFree<dim, value_type>> data,
235  const std::vector<unsigned int> &selected_row_blocks =
236  std::vector<unsigned int>(),
237  const std::vector<unsigned int> &selected_column_blocks =
238  std::vector<unsigned int>());
239 
253  void
254  initialize(std::shared_ptr<const MatrixFree<dim, value_type>> data,
255  const MGConstrainedDoFs & mg_constrained_dofs,
256  const unsigned int level,
257  const std::vector<unsigned int> &selected_row_blocks =
258  std::vector<unsigned int>());
259 
274  void
275  initialize(std::shared_ptr<const MatrixFree<dim, value_type>> data_,
276  const std::vector<MGConstrainedDoFs> &mg_constrained_dofs,
277  const unsigned int level,
278  const std::vector<unsigned int> & selected_row_blocks =
279  std::vector<unsigned int>());
280 
284  size_type
285  m() const;
286 
290  size_type
291  n() const;
292 
296  void
297  vmult_interface_down(VectorType &dst, const VectorType &src) const;
298 
302  void
303  vmult_interface_up(VectorType &dst, const VectorType &src) const;
304 
308  void
309  vmult(VectorType &dst, const VectorType &src) const;
310 
314  void
315  Tvmult(VectorType &dst, const VectorType &src) const;
316 
320  void
321  vmult_add(VectorType &dst, const VectorType &src) const;
322 
326  void
327  Tvmult_add(VectorType &dst, const VectorType &src) const;
328 
333  value_type
334  el(const unsigned int row, const unsigned int col) const;
335 
340  virtual std::size_t
341  memory_consumption() const;
342 
346  void
347  initialize_dof_vector(VectorType &vec) const;
348 
356  virtual void
357  compute_diagonal() = 0;
358 
362  std::shared_ptr<const MatrixFree<dim, value_type>>
363  get_matrix_free() const;
364 
368  const std::shared_ptr<DiagonalMatrix<VectorType>> &
369  get_matrix_diagonal_inverse() const;
370 
374  const std::shared_ptr<DiagonalMatrix<VectorType>> &
375  get_matrix_diagonal() const;
376 
377 
383  void
384  precondition_Jacobi(VectorType & dst,
385  const VectorType &src,
386  const value_type omega) const;
387 
388  protected:
393  void
394  preprocess_constraints(VectorType &dst, const VectorType &src) const;
395 
400  void
401  postprocess_constraints(VectorType &dst, const VectorType &src) const;
402 
407  void
408  set_constrained_entries_to_one(VectorType &dst) const;
409 
413  virtual void
414  apply_add(VectorType &dst, const VectorType &src) const = 0;
415 
421  virtual void
422  Tapply_add(VectorType &dst, const VectorType &src) const;
423 
427  std::shared_ptr<const MatrixFree<dim, value_type>> data;
428 
433  std::shared_ptr<DiagonalMatrix<VectorType>> diagonal_entries;
434 
439  std::shared_ptr<DiagonalMatrix<VectorType>> inverse_diagonal_entries;
440 
445  std::vector<unsigned int> selected_rows;
446 
451  std::vector<unsigned int> selected_columns;
452 
453  private:
457  std::vector<std::vector<unsigned int>> edge_constrained_indices;
458 
462  mutable std::vector<std::vector<std::pair<value_type, value_type>>>
464 
470 
475  void
476  mult_add(VectorType & dst,
477  const VectorType &src,
478  const bool transpose) const;
479 
487  void
488  adjust_ghost_range_if_necessary(const VectorType &vec,
489  const bool is_row) const;
490  };
491 
492 
493 
530  template <typename OperatorType>
532  {
533  public:
537  using value_type = typename OperatorType::value_type;
538 
542  using size_type = typename OperatorType::size_type;
543 
548 
552  void
553  clear();
554 
558  void
559  initialize(const OperatorType &operator_in);
560 
564  template <typename VectorType>
565  void
566  vmult(VectorType &dst, const VectorType &src) const;
567 
571  template <typename VectorType>
572  void
573  Tvmult(VectorType &dst, const VectorType &src) const;
574 
578  template <typename VectorType>
579  void
580  initialize_dof_vector(VectorType &vec) const;
581 
582 
583  private:
588  };
589 
590 
591 
611  template <int dim,
612  int fe_degree,
613  int n_components = 1,
614  typename Number = double,
615  typename VectorizedArrayType = VectorizedArray<Number>>
617  {
618  static_assert(
619  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
620  "Type of Number and of VectorizedArrayType do not match.");
621 
622  public:
628  const FEEvaluationBase<dim,
629  n_components,
630  Number,
631  false,
632  VectorizedArrayType> &fe_eval);
633 
642  void
643  apply(const AlignedVector<VectorizedArrayType> &inverse_coefficient,
644  const unsigned int n_actual_components,
645  const VectorizedArrayType * in_array,
646  VectorizedArrayType * out_array) const;
647 
681  void
682  transform_from_q_points_to_basis(const unsigned int n_actual_components,
683  const VectorizedArrayType *in_array,
684  VectorizedArrayType *out_array) const;
685 
691  void
692  fill_inverse_JxW_values(
693  AlignedVector<VectorizedArrayType> &inverse_jxw) const;
694 
695  private:
699  const FEEvaluationBase<dim,
700  n_components,
701  Number,
702  false,
703  VectorizedArrayType> &fe_eval;
704 
709  };
710 
711 
712 
722  template <int dim,
723  int fe_degree,
724  int n_q_points_1d = fe_degree + 1,
725  int n_components = 1,
726  typename VectorType = LinearAlgebra::distributed::Vector<double>>
727  class MassOperator : public Base<dim, VectorType>
728  {
729  public:
734 
739 
743  MassOperator();
744 
749  virtual void
750  compute_diagonal() override;
751 
752  private:
758  virtual void
759  apply_add(VectorType &dst, const VectorType &src) const override;
760 
764  void
765  local_apply_cell(
766  const MatrixFree<dim, value_type> & data,
767  VectorType & dst,
768  const VectorType & src,
769  const std::pair<unsigned int, unsigned int> &cell_range) const;
770  };
771 
772 
773 
786  template <int dim,
787  int fe_degree,
788  int n_q_points_1d = fe_degree + 1,
789  int n_components = 1,
790  typename VectorType = LinearAlgebra::distributed::Vector<double>>
791  class LaplaceOperator : public Base<dim, VectorType>
792  {
793  public:
798 
803 
807  LaplaceOperator();
808 
815  virtual void
816  compute_diagonal();
817 
864  void
865  set_coefficient(const std::shared_ptr<Table<2, VectorizedArray<value_type>>>
866  &scalar_coefficient);
867 
868  virtual void
869  clear();
870 
877  std::shared_ptr<Table<2, VectorizedArray<value_type>>>
878  get_coefficient();
879 
880  private:
886  virtual void
887  apply_add(VectorType &dst, const VectorType &src) const;
888 
892  void
893  local_apply_cell(
894  const MatrixFree<dim, value_type> & data,
895  VectorType & dst,
896  const VectorType & src,
897  const std::pair<unsigned int, unsigned int> &cell_range) const;
898 
902  void
903  local_diagonal_cell(
904  const MatrixFree<dim, value_type> &data,
905  VectorType & dst,
906  const VectorType &,
907  const std::pair<unsigned int, unsigned int> &cell_range) const;
908 
912  void
913  do_operation_on_cell(
915  & phi,
916  const unsigned int cell) const;
917 
921  std::shared_ptr<Table<2, VectorizedArray<value_type>>> scalar_coefficient;
922  };
923 
924 
925 
926  // ------------------------------------ inline functions ---------------------
927 
928  template <int dim,
929  int fe_degree,
930  int n_components,
931  typename Number,
932  typename VectorizedArrayType>
933  inline CellwiseInverseMassMatrix<dim,
934  fe_degree,
935  n_components,
936  Number,
937  VectorizedArrayType>::
938  CellwiseInverseMassMatrix(
939  const FEEvaluationBase<dim,
940  n_components,
941  Number,
942  false,
943  VectorizedArrayType> &fe_eval)
944  : fe_eval(fe_eval)
945  {
946  FullMatrix<double> shapes_1d(fe_degree + 1, fe_degree + 1);
947  for (unsigned int i = 0, c = 0; i < shapes_1d.m(); ++i)
948  for (unsigned int j = 0; j < shapes_1d.n(); ++j, ++c)
949  shapes_1d(i, j) = fe_eval.get_shape_info().shape_values[c][0];
950  shapes_1d.gauss_jordan();
951  const unsigned int stride = (fe_degree + 2) / 2;
952  inverse_shape.resize(stride * (fe_degree + 1));
953  for (unsigned int i = 0; i < stride; ++i)
954  for (unsigned int q = 0; q < (fe_degree + 2) / 2; ++q)
955  {
956  inverse_shape[i * stride + q] =
957  0.5 * (shapes_1d(i, q) + shapes_1d(i, fe_degree - q));
958  inverse_shape[(fe_degree - i) * stride + q] =
959  0.5 * (shapes_1d(i, q) - shapes_1d(i, fe_degree - q));
960  }
961  if (fe_degree % 2 == 0)
962  for (unsigned int q = 0; q < (fe_degree + 2) / 2; ++q)
963  inverse_shape[fe_degree / 2 * stride + q] = shapes_1d(fe_degree / 2, q);
964  }
965 
966 
967 
968  template <int dim,
969  int fe_degree,
970  int n_components,
971  typename Number,
972  typename VectorizedArrayType>
973  inline void
975  fe_degree,
976  n_components,
977  Number,
978  VectorizedArrayType>::
980  AlignedVector<VectorizedArrayType> &inverse_jxw) const
981  {
982  constexpr unsigned int dofs_per_component_on_cell =
983  Utilities::pow(fe_degree + 1, dim);
984  Assert(inverse_jxw.size() > 0 &&
985  inverse_jxw.size() % dofs_per_component_on_cell == 0,
986  ExcMessage(
987  "Expected diagonal to be a multiple of scalar dof per cells"));
988 
989  // temporarily reduce size of inverse_jxw to dofs_per_cell to get JxW values
990  // from fe_eval (will not reallocate any memory)
991  for (unsigned int q = 0; q < dofs_per_component_on_cell; ++q)
992  inverse_jxw[q] = 1. / fe_eval.JxW(q);
993  // copy values to rest of vector
994  for (unsigned int q = dofs_per_component_on_cell; q < inverse_jxw.size();)
995  for (unsigned int i = 0; i < dofs_per_component_on_cell; ++i, ++q)
996  inverse_jxw[q] = inverse_jxw[i];
997  }
998 
999 
1000 
1001  template <int dim,
1002  int fe_degree,
1003  int n_components,
1004  typename Number,
1005  typename VectorizedArrayType>
1006  inline void
1008  fe_degree,
1009  n_components,
1010  Number,
1011  VectorizedArrayType>::
1012  apply(const AlignedVector<VectorizedArrayType> &inverse_coefficients,
1013  const unsigned int n_actual_components,
1014  const VectorizedArrayType * in_array,
1015  VectorizedArrayType * out_array) const
1016  {
1017  constexpr unsigned int dofs_per_component =
1018  Utilities::pow(fe_degree + 1, dim);
1019  Assert(inverse_coefficients.size() > 0 &&
1020  inverse_coefficients.size() % dofs_per_component == 0,
1021  ExcMessage(
1022  "Expected diagonal to be a multiple of scalar dof per cells"));
1023  if (inverse_coefficients.size() != dofs_per_component)
1024  AssertDimension(n_actual_components * dofs_per_component,
1025  inverse_coefficients.size());
1026 
1027  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
1028 
1030  dim,
1031  fe_degree + 1,
1032  fe_degree + 1,
1033  VectorizedArrayType>
1035 
1036  const unsigned int shift_coefficient =
1037  inverse_coefficients.size() > dofs_per_component ? dofs_per_component : 0;
1038  const VectorizedArrayType *inv_coefficient = inverse_coefficients.data();
1039  VectorizedArrayType temp_data_field[dofs_per_component];
1040  for (unsigned int d = 0; d < n_actual_components; ++d)
1041  {
1042  const VectorizedArrayType *in = in_array + d * dofs_per_component;
1043  VectorizedArrayType * out = out_array + d * dofs_per_component;
1044  // Need to select 'apply' method with hessian slot because values
1045  // assume symmetries that do not exist in the inverse shapes
1046  evaluator.template hessians<0, false, false>(in, temp_data_field);
1047  if (dim > 1)
1048  {
1049  evaluator.template hessians<1, false, false>(temp_data_field, out);
1050 
1051  if (dim == 3)
1052  {
1053  evaluator.template hessians<2, false, false>(out,
1054  temp_data_field);
1055  for (unsigned int q = 0; q < dofs_per_component; ++q)
1056  temp_data_field[q] *= inv_coefficient[q];
1057  evaluator.template hessians<2, true, false>(temp_data_field,
1058  out);
1059  }
1060  else if (dim == 2)
1061  for (unsigned int q = 0; q < dofs_per_component; ++q)
1062  out[q] *= inv_coefficient[q];
1063 
1064  evaluator.template hessians<1, true, false>(out, temp_data_field);
1065  }
1066  else
1067  {
1068  for (unsigned int q = 0; q < dofs_per_component; ++q)
1069  temp_data_field[q] *= inv_coefficient[q];
1070  }
1071  evaluator.template hessians<0, true, false>(temp_data_field, out);
1072 
1073  inv_coefficient += shift_coefficient;
1074  }
1075  }
1076 
1077 
1078 
1079  template <int dim,
1080  int fe_degree,
1081  int n_components,
1082  typename Number,
1083  typename VectorizedArrayType>
1084  inline void
1086  fe_degree,
1087  n_components,
1088  Number,
1089  VectorizedArrayType>::
1090  transform_from_q_points_to_basis(const unsigned int n_actual_components,
1091  const VectorizedArrayType *in_array,
1092  VectorizedArrayType * out_array) const
1093  {
1094  constexpr unsigned int dofs_per_cell = Utilities::pow(fe_degree + 1, dim);
1096  dim,
1097  fe_degree + 1,
1098  fe_degree + 1,
1099  VectorizedArrayType>
1102  inverse_shape);
1103 
1104  for (unsigned int d = 0; d < n_actual_components; ++d)
1105  {
1106  const VectorizedArrayType *in = in_array + d * dofs_per_cell;
1107  VectorizedArrayType * out = out_array + d * dofs_per_cell;
1108 
1109  if (dim == 3)
1110  {
1111  evaluator.template hessians<2, true, false>(in, out);
1112  evaluator.template hessians<1, true, false>(out, out);
1113  evaluator.template hessians<0, true, false>(out, out);
1114  }
1115  if (dim == 2)
1116  {
1117  evaluator.template hessians<1, true, false>(in, out);
1118  evaluator.template hessians<0, true, false>(out, out);
1119  }
1120  if (dim == 1)
1121  evaluator.template hessians<0, true, false>(in, out);
1122  }
1123  }
1124 
1125 
1126 
1127  //----------------- Base operator -----------------------------
1128  template <int dim, typename VectorType>
1130  : Subscriptor()
1131  , have_interface_matrices(false)
1132  {}
1133 
1134 
1135 
1136  template <int dim, typename VectorType>
1139  {
1140  Assert(data.get() != nullptr, ExcNotInitialized());
1141  typename Base<dim, VectorType>::size_type total_size = 0;
1142  for (unsigned int i = 0; i < selected_rows.size(); ++i)
1143  total_size += data->get_vector_partitioner(selected_rows[i])->size();
1144  return total_size;
1145  }
1146 
1147 
1148 
1149  template <int dim, typename VectorType>
1152  {
1153  Assert(data.get() != nullptr, ExcNotInitialized());
1154  typename Base<dim, VectorType>::size_type total_size = 0;
1155  for (unsigned int i = 0; i < selected_columns.size(); ++i)
1156  total_size += data->get_vector_partitioner(selected_columns[i])->size();
1157  return total_size;
1158  }
1159 
1160 
1161 
1162  template <int dim, typename VectorType>
1163  void
1165  {
1166  data.reset();
1167  inverse_diagonal_entries.reset();
1168  }
1169 
1170 
1171 
1172  template <int dim, typename VectorType>
1174  Base<dim, VectorType>::el(const unsigned int row,
1175  const unsigned int col) const
1176  {
1177  (void)col;
1178  Assert(row == col, ExcNotImplemented());
1179  Assert(inverse_diagonal_entries.get() != nullptr &&
1180  inverse_diagonal_entries->m() > 0,
1181  ExcNotInitialized());
1182  return 1.0 / (*inverse_diagonal_entries)(row, row);
1183  }
1184 
1185 
1186 
1187  template <int dim, typename VectorType>
1188  void
1190  {
1191  Assert(data.get() != nullptr, ExcNotInitialized());
1192  AssertDimension(BlockHelper::n_blocks(vec), selected_rows.size());
1193  for (unsigned int i = 0; i < BlockHelper::n_blocks(vec); ++i)
1194  {
1195  const unsigned int index = selected_rows[i];
1196  if (!BlockHelper::subblock(vec, index)
1197  .partitioners_are_compatible(
1198  *data->get_dof_info(index).vector_partitioner))
1199  data->initialize_dof_vector(BlockHelper::subblock(vec, index), index);
1200 
1201  Assert(BlockHelper::subblock(vec, index)
1202  .partitioners_are_globally_compatible(
1203  *data->get_dof_info(index).vector_partitioner),
1204  ExcInternalError());
1205  }
1206  BlockHelper::collect_sizes(vec);
1207  }
1208 
1209 
1210 
1211  template <int dim, typename VectorType>
1212  void
1214  std::shared_ptr<
1215  const MatrixFree<dim, typename Base<dim, VectorType>::value_type>> data_,
1216  const std::vector<unsigned int> &given_row_selection,
1217  const std::vector<unsigned int> &given_column_selection)
1218  {
1219  data = data_;
1220 
1221  selected_rows.clear();
1222  selected_columns.clear();
1223  if (given_row_selection.empty())
1224  for (unsigned int i = 0; i < data_->n_components(); ++i)
1225  selected_rows.push_back(i);
1226  else
1227  {
1228  for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1229  {
1230  AssertIndexRange(given_row_selection[i], data_->n_components());
1231  for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1232  if (j != i)
1233  Assert(given_row_selection[j] != given_row_selection[i],
1234  ExcMessage("Given row indices must be unique"));
1235 
1236  selected_rows.push_back(given_row_selection[i]);
1237  }
1238  }
1239  if (given_column_selection.size() == 0)
1241  else
1242  {
1243  for (unsigned int i = 0; i < given_column_selection.size(); ++i)
1244  {
1245  AssertIndexRange(given_column_selection[i], data_->n_components());
1246  for (unsigned int j = 0; j < given_column_selection.size(); ++j)
1247  if (j != i)
1248  Assert(given_column_selection[j] != given_column_selection[i],
1249  ExcMessage("Given column indices must be unique"));
1250 
1251  selected_columns.push_back(given_column_selection[i]);
1252  }
1253  }
1254 
1255  edge_constrained_indices.clear();
1256  edge_constrained_indices.resize(selected_rows.size());
1257  edge_constrained_values.clear();
1258  edge_constrained_values.resize(selected_rows.size());
1259  have_interface_matrices = false;
1260  }
1261 
1262 
1263 
1264  template <int dim, typename VectorType>
1265  void
1267  std::shared_ptr<
1268  const MatrixFree<dim, typename Base<dim, VectorType>::value_type>> data_,
1269  const MGConstrainedDoFs & mg_constrained_dofs,
1270  const unsigned int level,
1271  const std::vector<unsigned int> &given_row_selection)
1272  {
1273  std::vector<MGConstrainedDoFs> mg_constrained_dofs_vector(
1274  1, mg_constrained_dofs);
1275  initialize(data_, mg_constrained_dofs_vector, level, given_row_selection);
1276  }
1277 
1278 
1279 
1280  template <int dim, typename VectorType>
1281  void
1283  std::shared_ptr<const MatrixFree<dim, value_type>> data_,
1284  const std::vector<MGConstrainedDoFs> & mg_constrained_dofs,
1285  const unsigned int level,
1286  const std::vector<unsigned int> & given_row_selection)
1287  {
1289  ExcMessage("level is not set"));
1290 
1291  selected_rows.clear();
1292  selected_columns.clear();
1293  if (given_row_selection.empty())
1294  for (unsigned int i = 0; i < data_->n_components(); ++i)
1295  selected_rows.push_back(i);
1296  else
1297  {
1298  for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1299  {
1300  AssertIndexRange(given_row_selection[i], data_->n_components());
1301  for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1302  if (j != i)
1303  Assert(given_row_selection[j] != given_row_selection[i],
1304  ExcMessage("Given row indices must be unique"));
1305 
1306  selected_rows.push_back(given_row_selection[i]);
1307  }
1308  }
1310 
1311  AssertDimension(mg_constrained_dofs.size(), selected_rows.size());
1312  edge_constrained_indices.clear();
1313  edge_constrained_indices.resize(selected_rows.size());
1314  edge_constrained_values.clear();
1315  edge_constrained_values.resize(selected_rows.size());
1316 
1317  data = data_;
1318 
1319  for (unsigned int j = 0; j < selected_rows.size(); ++j)
1320  {
1321  if (data_->n_macro_cells() > 0)
1322  {
1323  AssertDimension(level, data_->get_cell_iterator(0, 0, j)->level());
1324  }
1325 
1326  // setup edge_constrained indices
1327  std::vector<types::global_dof_index> interface_indices;
1328  mg_constrained_dofs[j]
1329  .get_refinement_edge_indices(level)
1330  .fill_index_vector(interface_indices);
1331  edge_constrained_indices[j].clear();
1332  edge_constrained_indices[j].reserve(interface_indices.size());
1333  edge_constrained_values[j].resize(interface_indices.size());
1334  const IndexSet &locally_owned =
1335  data->get_dof_handler(selected_rows[j]).locally_owned_mg_dofs(level);
1336  for (const auto interface_index : interface_indices)
1337  if (locally_owned.is_element(interface_index))
1338  edge_constrained_indices[j].push_back(
1339  locally_owned.index_within_set(interface_index));
1342  static_cast<unsigned int>(edge_constrained_indices[j].size()),
1343  data->get_vector_partitioner()->get_mpi_communicator()) > 0;
1344  }
1345  }
1346 
1347 
1348 
1349  template <int dim, typename VectorType>
1350  void
1352  {
1353  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1354  {
1355  const std::vector<unsigned int> &constrained_dofs =
1356  data->get_constrained_dofs(selected_rows[j]);
1357  for (const auto constrained_dof : constrained_dofs)
1358  BlockHelper::subblock(dst, j).local_element(constrained_dof) = 1.;
1359  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1360  BlockHelper::subblock(dst, j).local_element(
1361  edge_constrained_indices[j][i]) = 1.;
1362  }
1363  }
1364 
1365 
1366 
1367  template <int dim, typename VectorType>
1368  void
1369  Base<dim, VectorType>::vmult(VectorType &dst, const VectorType &src) const
1370  {
1371  using Number = typename Base<dim, VectorType>::value_type;
1372  dst = Number(0.);
1373  vmult_add(dst, src);
1374  }
1375 
1376 
1377 
1378  template <int dim, typename VectorType>
1379  void
1380  Base<dim, VectorType>::vmult_add(VectorType &dst, const VectorType &src) const
1381  {
1382  mult_add(dst, src, false);
1383  }
1384 
1385 
1386 
1387  template <int dim, typename VectorType>
1388  void
1390  const VectorType &src) const
1391  {
1392  mult_add(dst, src, true);
1393  }
1394 
1395 
1396 
1397  template <int dim, typename VectorType>
1398  void
1400  const VectorType &src,
1401  const bool is_row) const
1402  {
1403  using Number = typename Base<dim, VectorType>::value_type;
1404  for (unsigned int i = 0; i < BlockHelper::n_blocks(src); ++i)
1405  {
1406  const unsigned int mf_component =
1407  is_row ? selected_rows[i] : selected_columns[i];
1408  // If both vectors use the same partitioner -> done
1409  if (BlockHelper::subblock(src, i).get_partitioner().get() ==
1410  data->get_dof_info(mf_component).vector_partitioner.get())
1411  continue;
1412 
1413  // If not, assert that the local ranges are the same and reset to the
1414  // current partitioner
1415  Assert(
1416  BlockHelper::subblock(src, i).get_partitioner()->local_size() ==
1417  data->get_dof_info(mf_component).vector_partitioner->local_size(),
1418  ExcMessage("The vector passed to the vmult() function does not have "
1419  "the correct size for compatibility with MatrixFree."));
1420 
1421  // copy the vector content to a temporary vector so that it does not get
1422  // lost
1424  BlockHelper::subblock(src, i));
1425  BlockHelper::subblock(const_cast<VectorType &>(src), i)
1426  .reinit(data->get_dof_info(mf_component).vector_partitioner);
1427  BlockHelper::subblock(const_cast<VectorType &>(src), i)
1428  .copy_locally_owned_data_from(copy_vec);
1429  }
1430  }
1431 
1432 
1433 
1434  template <int dim, typename VectorType>
1435  void
1437  const VectorType &src) const
1438  {
1439  using Number = typename Base<dim, VectorType>::value_type;
1440  adjust_ghost_range_if_necessary(src, false);
1442 
1443  // set zero Dirichlet values on the input vector (and remember the src and
1444  // dst values because we need to reset them at the end)
1445  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1446  {
1447  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1448  {
1449  edge_constrained_values[j][i] = std::pair<Number, Number>(
1450  BlockHelper::subblock(src, j).local_element(
1451  edge_constrained_indices[j][i]),
1452  BlockHelper::subblock(dst, j).local_element(
1453  edge_constrained_indices[j][i]));
1454  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1455  .local_element(edge_constrained_indices[j][i]) = 0.;
1456  }
1457  }
1458  }
1459 
1460 
1461 
1462  template <int dim, typename VectorType>
1463  void
1465  const VectorType &src,
1466  const bool transpose) const
1467  {
1468  AssertDimension(dst.size(), src.size());
1469  AssertDimension(BlockHelper::n_blocks(dst), BlockHelper::n_blocks(src));
1470  AssertDimension(BlockHelper::n_blocks(dst), selected_rows.size());
1471  preprocess_constraints(dst, src);
1472  if (transpose)
1473  Tapply_add(dst, src);
1474  else
1475  apply_add(dst, src);
1476  postprocess_constraints(dst, src);
1477  }
1478 
1479 
1480 
1481  template <int dim, typename VectorType>
1482  void
1484  const VectorType &src) const
1485  {
1486  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1487  {
1488  const std::vector<unsigned int> &constrained_dofs =
1489  data->get_constrained_dofs(selected_rows[j]);
1490  for (const auto constrained_dof : constrained_dofs)
1491  BlockHelper::subblock(dst, j).local_element(constrained_dof) +=
1492  BlockHelper::subblock(src, j).local_element(constrained_dof);
1493  }
1494 
1495  // reset edge constrained values, multiply by unit matrix and add into
1496  // destination
1497  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1498  {
1499  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1500  {
1501  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1502  .local_element(edge_constrained_indices[j][i]) =
1503  edge_constrained_values[j][i].first;
1504  BlockHelper::subblock(dst, j).local_element(
1505  edge_constrained_indices[j][i]) =
1506  edge_constrained_values[j][i].second +
1507  edge_constrained_values[j][i].first;
1508  }
1509  }
1510  }
1511 
1512 
1513 
1514  template <int dim, typename VectorType>
1515  void
1517  const VectorType &src) const
1518  {
1519  using Number = typename Base<dim, VectorType>::value_type;
1520  AssertDimension(dst.size(), src.size());
1521  adjust_ghost_range_if_necessary(src, false);
1523 
1524  dst = Number(0.);
1525 
1527  return;
1528 
1529  // set zero Dirichlet values on the input vector (and remember the src and
1530  // dst values because we need to reset them at the end)
1531  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1532  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1533  {
1534  edge_constrained_values[j][i] = std::pair<Number, Number>(
1535  BlockHelper::subblock(src, j).local_element(
1536  edge_constrained_indices[j][i]),
1537  BlockHelper::subblock(dst, j).local_element(
1538  edge_constrained_indices[j][i]));
1539  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1540  .local_element(edge_constrained_indices[j][i]) = 0.;
1541  }
1542 
1543  apply_add(dst, src);
1544 
1545  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1546  {
1547  unsigned int c = 0;
1548  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1549  {
1550  for (; c < edge_constrained_indices[j][i]; ++c)
1551  BlockHelper::subblock(dst, j).local_element(c) = 0.;
1552  ++c;
1553 
1554  // reset the src values
1555  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1556  .local_element(edge_constrained_indices[j][i]) =
1557  edge_constrained_values[j][i].first;
1558  }
1559  for (; c < BlockHelper::subblock(dst, j).local_size(); ++c)
1560  BlockHelper::subblock(dst, j).local_element(c) = 0.;
1561  }
1562  }
1563 
1564 
1565 
1566  template <int dim, typename VectorType>
1567  void
1569  const VectorType &src) const
1570  {
1571  using Number = typename Base<dim, VectorType>::value_type;
1572  AssertDimension(dst.size(), src.size());
1573  adjust_ghost_range_if_necessary(src, false);
1575 
1576  dst = Number(0.);
1577 
1579  return;
1580 
1581  VectorType src_cpy(src);
1582  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1583  {
1584  unsigned int c = 0;
1585  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1586  {
1587  for (; c < edge_constrained_indices[j][i]; ++c)
1588  BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1589  ++c;
1590  }
1591  for (; c < BlockHelper::subblock(src_cpy, j).local_size(); ++c)
1592  BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1593  }
1594 
1595  apply_add(dst, src_cpy);
1596 
1597  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1598  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1599  BlockHelper::subblock(dst, j).local_element(
1600  edge_constrained_indices[j][i]) = 0.;
1601  }
1602 
1603 
1604 
1605  template <int dim, typename VectorType>
1606  void
1607  Base<dim, VectorType>::Tvmult(VectorType &dst, const VectorType &src) const
1608  {
1609  using Number = typename Base<dim, VectorType>::value_type;
1610  dst = Number(0.);
1611  Tvmult_add(dst, src);
1612  }
1613 
1614 
1615 
1616  template <int dim, typename VectorType>
1617  std::size_t
1619  {
1620  return inverse_diagonal_entries.get() != nullptr ?
1621  inverse_diagonal_entries->memory_consumption() :
1622  sizeof(*this);
1623  }
1624 
1625 
1626 
1627  template <int dim, typename VectorType>
1628  std::shared_ptr<
1631  {
1632  return data;
1633  }
1634 
1635 
1636 
1637  template <int dim, typename VectorType>
1638  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1640  {
1641  Assert(inverse_diagonal_entries.get() != nullptr &&
1642  inverse_diagonal_entries->m() > 0,
1643  ExcNotInitialized());
1644  return inverse_diagonal_entries;
1645  }
1646 
1647 
1648 
1649  template <int dim, typename VectorType>
1650  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1652  {
1653  Assert(diagonal_entries.get() != nullptr && diagonal_entries->m() > 0,
1654  ExcNotInitialized());
1655  return diagonal_entries;
1656  }
1657 
1658 
1659 
1660  template <int dim, typename VectorType>
1661  void
1663  const VectorType &src) const
1664  {
1665  apply_add(dst, src);
1666  }
1667 
1668 
1669 
1670  template <int dim, typename VectorType>
1671  void
1673  VectorType & dst,
1674  const VectorType & src,
1675  const typename Base<dim, VectorType>::value_type omega) const
1676  {
1678  ExcNotInitialized());
1679  inverse_diagonal_entries->vmult(dst, src);
1680  dst *= omega;
1681  }
1682 
1683 
1684 
1685  //------------------------- MGInterfaceOperator ------------------------------
1686 
1687  template <typename OperatorType>
1689  : Subscriptor()
1690  , mf_base_operator(nullptr)
1691  {}
1692 
1693 
1694 
1695  template <typename OperatorType>
1696  void
1698  {
1699  mf_base_operator = nullptr;
1700  }
1701 
1702 
1703 
1704  template <typename OperatorType>
1705  void
1706  MGInterfaceOperator<OperatorType>::initialize(const OperatorType &operator_in)
1707  {
1708  mf_base_operator = &operator_in;
1709  }
1710 
1711 
1712 
1713  template <typename OperatorType>
1714  template <typename VectorType>
1715  void
1717  const VectorType &src) const
1718  {
1719 #ifndef DEAL_II_MSVC
1720  static_assert(
1721  std::is_same<typename VectorType::value_type, value_type>::value,
1722  "The vector type must be based on the same value type as this "
1723  "operator");
1724 #endif
1725 
1726  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1727 
1728  mf_base_operator->vmult_interface_down(dst, src);
1729  }
1730 
1731 
1732 
1733  template <typename OperatorType>
1734  template <typename VectorType>
1735  void
1737  const VectorType &src) const
1738  {
1739 #ifndef DEAL_II_MSVC
1740  static_assert(
1741  std::is_same<typename VectorType::value_type, value_type>::value,
1742  "The vector type must be based on the same value type as this "
1743  "operator");
1744 #endif
1745 
1746  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1747 
1748  mf_base_operator->vmult_interface_up(dst, src);
1749  }
1750 
1751 
1752 
1753  template <typename OperatorType>
1754  template <typename VectorType>
1755  void
1757  VectorType &vec) const
1758  {
1759  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1760 
1761  mf_base_operator->initialize_dof_vector(vec);
1762  }
1763 
1764 
1765 
1766  //-----------------------------MassOperator----------------------------------
1767 
1768  template <int dim,
1769  int fe_degree,
1770  int n_q_points_1d,
1771  int n_components,
1772  typename VectorType>
1775  : Base<dim, VectorType>()
1776  {}
1777 
1778 
1779 
1780  template <int dim,
1781  int fe_degree,
1782  int n_q_points_1d,
1783  int n_components,
1784  typename VectorType>
1785  void
1788  {
1789  using Number = typename Base<dim, VectorType>::value_type;
1790  Assert((Base<dim, VectorType>::data.get() != nullptr), ExcNotInitialized());
1791 
1793  this->diagonal_entries.reset(new DiagonalMatrix<VectorType>());
1794  VectorType &inverse_diagonal_vector =
1795  this->inverse_diagonal_entries->get_vector();
1796  VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1797  this->initialize_dof_vector(inverse_diagonal_vector);
1798  this->initialize_dof_vector(diagonal_vector);
1799  inverse_diagonal_vector = Number(1.);
1800  apply_add(diagonal_vector, inverse_diagonal_vector);
1801 
1802  this->set_constrained_entries_to_one(diagonal_vector);
1803  inverse_diagonal_vector = diagonal_vector;
1804 
1805  const unsigned int local_size = inverse_diagonal_vector.local_size();
1806  for (unsigned int i = 0; i < local_size; ++i)
1807  inverse_diagonal_vector.local_element(i) =
1808  Number(1.) / inverse_diagonal_vector.local_element(i);
1809 
1810  inverse_diagonal_vector.update_ghost_values();
1811  diagonal_vector.update_ghost_values();
1812  }
1813 
1814 
1815 
1816  template <int dim,
1817  int fe_degree,
1818  int n_q_points_1d,
1819  int n_components,
1820  typename VectorType>
1821  void
1823  apply_add(VectorType &dst, const VectorType &src) const
1824  {
1826  this,
1827  dst,
1828  src);
1829  }
1830 
1831 
1832 
1833  template <int dim,
1834  int fe_degree,
1835  int n_q_points_1d,
1836  int n_components,
1837  typename VectorType>
1838  void
1841  const MatrixFree<dim, typename Base<dim, VectorType>::value_type> &data,
1842  VectorType & dst,
1843  const VectorType & src,
1844  const std::pair<unsigned int, unsigned int> &cell_range) const
1845  {
1846  using Number = typename Base<dim, VectorType>::value_type;
1848  data, this->selected_rows[0]);
1849  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
1850  {
1851  phi.reinit(cell);
1852  phi.read_dof_values(src);
1853  phi.evaluate(true, false, false);
1854  for (unsigned int q = 0; q < phi.n_q_points; ++q)
1855  phi.submit_value(phi.get_value(q), q);
1856  phi.integrate(true, false);
1857  phi.distribute_local_to_global(dst);
1858  }
1859  }
1860 
1861 
1862  //-----------------------------LaplaceOperator----------------------------------
1863 
1864  template <int dim,
1865  int fe_degree,
1866  int n_q_points_1d,
1867  int n_components,
1868  typename VectorType>
1871  : Base<dim, VectorType>()
1872  {}
1873 
1874 
1875 
1876  template <int dim,
1877  int fe_degree,
1878  int n_q_points_1d,
1879  int n_components,
1880  typename VectorType>
1881  void
1884  {
1886  scalar_coefficient.reset();
1887  }
1888 
1889 
1890 
1891  template <int dim,
1892  int fe_degree,
1893  int n_q_points_1d,
1894  int n_components,
1895  typename VectorType>
1896  void
1899  const std::shared_ptr<
1901  &scalar_coefficient_)
1902  {
1903  scalar_coefficient = scalar_coefficient_;
1904  }
1905 
1906 
1907 
1908  template <int dim,
1909  int fe_degree,
1910  int n_q_points_1d,
1911  int n_components,
1912  typename VectorType>
1913  std::shared_ptr<
1914  Table<2,
1915  VectorizedArray<typename LaplaceOperator<dim,
1916  fe_degree,
1917  n_q_points_1d,
1918  n_components,
1919  VectorType>::value_type>>>
1922  {
1924  return scalar_coefficient;
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  void
1937  {
1938  using Number = typename Base<dim, VectorType>::value_type;
1939  Assert((Base<dim, VectorType>::data.get() != nullptr), ExcNotInitialized());
1940 
1942  this->diagonal_entries.reset(new DiagonalMatrix<VectorType>());
1943  VectorType &inverse_diagonal_vector =
1944  this->inverse_diagonal_entries->get_vector();
1945  VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1946  this->initialize_dof_vector(inverse_diagonal_vector);
1947  this->initialize_dof_vector(diagonal_vector);
1948 
1949  this->data->cell_loop(&LaplaceOperator::local_diagonal_cell,
1950  this,
1951  diagonal_vector,
1952  /*unused*/ diagonal_vector);
1953  this->set_constrained_entries_to_one(diagonal_vector);
1954 
1955  inverse_diagonal_vector = diagonal_vector;
1956 
1957  for (unsigned int i = 0; i < inverse_diagonal_vector.local_size(); ++i)
1958  if (std::abs(inverse_diagonal_vector.local_element(i)) >
1959  std::sqrt(std::numeric_limits<Number>::epsilon()))
1960  inverse_diagonal_vector.local_element(i) =
1961  1. / inverse_diagonal_vector.local_element(i);
1962  else
1963  inverse_diagonal_vector.local_element(i) = 1.;
1964 
1965  inverse_diagonal_vector.update_ghost_values();
1966  diagonal_vector.update_ghost_values();
1967  }
1968 
1969 
1970 
1971  template <int dim,
1972  int fe_degree,
1973  int n_q_points_1d,
1974  int n_components,
1975  typename VectorType>
1976  void
1978  apply_add(VectorType &dst, const VectorType &src) const
1979  {
1981  this,
1982  dst,
1983  src);
1984  }
1985 
1986  namespace Implementation
1987  {
1988  template <typename VectorizedArrayType>
1989  bool
1990  non_negative(const VectorizedArrayType &n)
1991  {
1992  for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
1993  if (n[v] < 0.)
1994  return false;
1995 
1996  return true;
1997  }
1998  } // namespace Implementation
1999 
2000 
2001 
2002  template <int dim,
2003  int fe_degree,
2004  int n_q_points_1d,
2005  int n_components,
2006  typename VectorType>
2007  void
2010  FEEvaluation<dim,
2011  fe_degree,
2012  n_q_points_1d,
2013  n_components,
2014  typename Base<dim, VectorType>::value_type> &phi,
2015  const unsigned int cell) const
2016  {
2017  phi.evaluate(false, true, false);
2018  if (scalar_coefficient.get())
2019  {
2020  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2021  {
2022  Assert(Implementation::non_negative((*scalar_coefficient)(cell, q)),
2023  ExcMessage("Coefficient must be non-negative"));
2024  phi.submit_gradient((*scalar_coefficient)(cell, q) *
2025  phi.get_gradient(q),
2026  q);
2027  }
2028  }
2029  else
2030  {
2031  for (unsigned int q = 0; q < phi.n_q_points; ++q)
2032  {
2033  phi.submit_gradient(phi.get_gradient(q), q);
2034  }
2035  }
2036  phi.integrate(false, true);
2037  }
2038 
2039 
2040 
2041  template <int dim,
2042  int fe_degree,
2043  int n_q_points_1d,
2044  int n_components,
2045  typename VectorType>
2046  void
2049  const MatrixFree<dim, typename Base<dim, VectorType>::value_type> &data,
2050  VectorType & dst,
2051  const VectorType & src,
2052  const std::pair<unsigned int, unsigned int> &cell_range) const
2053  {
2054  using Number = typename Base<dim, VectorType>::value_type;
2056  data, this->selected_rows[0]);
2057  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2058  {
2059  phi.reinit(cell);
2060  phi.read_dof_values(src);
2061  do_operation_on_cell(phi, cell);
2062  phi.distribute_local_to_global(dst);
2063  }
2064  }
2065 
2066 
2067  template <int dim,
2068  int fe_degree,
2069  int n_q_points_1d,
2070  int n_components,
2071  typename VectorType>
2072  void
2075  const MatrixFree<dim, typename Base<dim, VectorType>::value_type> &data,
2076  VectorType & dst,
2077  const VectorType &,
2078  const std::pair<unsigned int, unsigned int> &cell_range) const
2079  {
2080  using Number = typename Base<dim, VectorType>::value_type;
2081  using vectorized_value_type =
2083  vectorized_value_type;
2084 
2086  data, this->selected_rows[0]);
2087  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2088  {
2089  phi.reinit(cell);
2090  vectorized_value_type local_diagonal_vector[phi.static_dofs_per_cell];
2091  for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2092  {
2093  for (unsigned int j = 0; j < phi.dofs_per_component; ++j)
2094  phi.begin_dof_values()[j] = vectorized_value_type();
2095  phi.begin_dof_values()[i] = 1.;
2096  do_operation_on_cell(phi, cell);
2097  local_diagonal_vector[i] = phi.begin_dof_values()[i];
2098  }
2099  for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2100  for (unsigned int c = 0; c < phi.n_components; ++c)
2101  phi.begin_dof_values()[i + c * phi.dofs_per_component] =
2102  local_diagonal_vector[i];
2103  phi.distribute_local_to_global(dst);
2104  }
2105  }
2106 
2107 
2108 } // end of namespace MatrixFreeOperators
2109 
2110 
2111 DEAL_II_NAMESPACE_CLOSE
2112 
2113 #endif
typename OperatorType::value_type value_type
Definition: operators.h:537
VectorizedArrayType JxW(const unsigned int q_index) const
size_type m() const
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:2009
void fill_inverse_JxW_values(AlignedVector< VectorizedArrayType > &inverse_jxw) const
Definition: operators.h:979
static const unsigned int invalid_unsigned_int
Definition: types.h:187
void vmult_interface_down(VectorType &dst, const VectorType &src) const
Definition: operators.h:1516
std::vector< unsigned int > selected_rows
Definition: operators.h:445
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
void integrate(const bool integrate_values, const bool integrate_gradients)
virtual void apply_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1978
void local_apply_cell(const MatrixFree< dim, value_type > &data, VectorType &dst, const VectorType &src, const std::pair< unsigned int, unsigned int > &cell_range) const
Definition: operators.h:2048
const FEEvaluationBase< dim, n_components, Number, false, VectorizedArrayType > & fe_eval
Definition: operators.h:703
typename VectorType::size_type size_type
Definition: operators.h:198
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal_inverse() const
Definition: operators.h:1639
void gauss_jordan()
typename VectorType::value_type value_type
Definition: operators.h:193
void precondition_Jacobi(VectorType &dst, const VectorType &src, const value_type omega) const
Definition: operators.h:1672
void read_dof_values(const VectorType &src, const unsigned int first_index=0)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
pointer data()
std::vector< std::vector< std::pair< value_type, value_type > > > edge_constrained_values
Definition: operators.h:463
void mult_add(VectorType &dst, const VectorType &src, const bool transpose) const
Definition: operators.h:1464
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal() const
Definition: operators.h:1651
void evaluate(const bool evaluate_values, const bool evaluate_gradients, const bool evaluate_hessians=false)
static ::ExceptionBase & ExcNotInitialized()
void vmult_interface_up(VectorType &dst, const VectorType &src) const
Definition: operators.h:1568
std::shared_ptr< DiagonalMatrix< VectorType > > inverse_diagonal_entries
Definition: operators.h:439
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
void set_constrained_entries_to_one(VectorType &dst) const
Definition: operators.h:1351
std::shared_ptr< DiagonalMatrix< VectorType > > diagonal_entries
Definition: operators.h:433
constexpr SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
typename OperatorType::size_type size_type
Definition: operators.h:542
void resize(const size_type size_in)
AlignedVector< Number > shape_values
Definition: shape_info.h:155
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1756
std::vector< unsigned int > selected_columns
Definition: operators.h:451
value_type el(const unsigned int row, const unsigned int col) const
Definition: operators.h:1174
static ::ExceptionBase & ExcMessage(std::string arg1)
AlignedVector< VectorizedArrayType > inverse_shape
Definition: operators.h:708
size_type n() const
Definition: operators.h:1151
size_type n() const
virtual void Tapply_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1662
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1607
#define Assert(cond, exc)
Definition: exceptions.h:1407
const unsigned int n_q_points
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1888
void local_diagonal_cell(const MatrixFree< dim, value_type > &data, VectorType &dst, const VectorType &, const std::pair< unsigned int, unsigned int > &cell_range) const
Definition: operators.h:2074
constexpr unsigned int pow(const unsigned int base, const int iexp)
Definition: utilities.h:430
void set_coefficient(const std::shared_ptr< Table< 2, VectorizedArray< value_type >>> &scalar_coefficient)
Definition: operators.h:1898
const unsigned int dofs_per_component
void reinit(const unsigned int cell_batch_index)
std::vector< std::vector< unsigned int > > edge_constrained_indices
Definition: operators.h:457
void Tvmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1389
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1369
size_type m() const
Definition: operators.h:1138
virtual void apply_add(VectorType &dst, const VectorType &src) const =0
std::shared_ptr< Table< 2, VectorizedArray< value_type > > > scalar_coefficient
Definition: operators.h:921
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1716
static constexpr unsigned int n_components
static constexpr unsigned int static_dofs_per_cell
virtual void compute_diagonal() override
Definition: operators.h:1787
void preprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1436
const VectorizedArrayType * begin_dof_values() const
void local_apply_cell(const MatrixFree< dim, value_type > &data, VectorType &dst, const VectorType &src, const std::pair< unsigned int, unsigned int > &cell_range) const
Definition: operators.h:1840
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info() const
std::shared_ptr< const MatrixFree< dim, value_type > > get_matrix_free() const
Definition: operators.h:1630
std::shared_ptr< Table< 2, VectorizedArray< value_type > > > get_coefficient()
Definition: operators.h:1921
void distribute_local_to_global(VectorType &dst, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::n_array_elements > &mask=std::bitset< VectorizedArrayType::n_array_elements >().flip()) const
void submit_value(const value_type val_in, const unsigned int q_point)
size_type size() const
virtual std::size_t memory_consumption() const
Definition: operators.h:1618
std::shared_ptr< const MatrixFree< dim, value_type > > data
Definition: operators.h:427
virtual void apply_add(VectorType &dst, const VectorType &src) const override
Definition: operators.h:1823
void postprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1483
static ::ExceptionBase & ExcNotImplemented()
Definition: table.h:699
bool is_element(const size_type index) const
Definition: index_set.h:1732
void adjust_ghost_range_if_necessary(const VectorType &vec, const bool is_row) const
Definition: operators.h:1399
void vmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1380
void initialize(const OperatorType &operator_in)
Definition: operators.h:1706
virtual void clear()
Definition: operators.h:1164
SmartPointer< const OperatorType > mf_base_operator
Definition: operators.h:587
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1736
void apply(const AlignedVector< VectorizedArrayType > &inverse_coefficient, const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition: operators.h:1012
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1189
T max(const T &t, const MPI_Comm &mpi_communicator)
value_type get_value(const unsigned int q_point) const
void initialize(std::shared_ptr< const MatrixFree< dim, value_type >> 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 >())
void transform_from_q_points_to_basis(const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition: operators.h:1090
static ::ExceptionBase & ExcInternalError()