Reference documentation for deal.II version 9.1.1
\(\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>> &
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>
616  {
617  public:
624 
633  void
634  apply(const AlignedVector<VectorizedArray<Number>> &inverse_coefficient,
635  const unsigned int n_actual_components,
636  const VectorizedArray<Number> * in_array,
637  VectorizedArray<Number> * out_array) const;
638 
672  void
673  transform_from_q_points_to_basis(const unsigned int n_actual_components,
674  const VectorizedArray<Number> *in_array,
675  VectorizedArray<Number> *out_array) const;
676 
682  void
684  AlignedVector<VectorizedArray<Number>> &inverse_jxw) const;
685 
686  private:
691 
696  };
697 
698 
699 
709  template <int dim,
710  int fe_degree,
711  int n_q_points_1d = fe_degree + 1,
712  int n_components = 1,
713  typename VectorType = LinearAlgebra::distributed::Vector<double>>
714  class MassOperator : public Base<dim, VectorType>
715  {
716  public:
721 
726 
730  MassOperator();
731 
736  virtual void
737  compute_diagonal() override;
738 
739  private:
745  virtual void
746  apply_add(VectorType &dst, const VectorType &src) const override;
747 
751  void
754  VectorType & dst,
755  const VectorType & src,
756  const std::pair<unsigned int, unsigned int> &cell_range) const;
757  };
758 
759 
760 
773  template <int dim,
774  int fe_degree,
775  int n_q_points_1d = fe_degree + 1,
776  int n_components = 1,
777  typename VectorType = LinearAlgebra::distributed::Vector<double>>
778  class LaplaceOperator : public Base<dim, VectorType>
779  {
780  public:
785 
790 
794  LaplaceOperator();
795 
802  virtual void
804 
851  void
852  set_coefficient(const std::shared_ptr<Table<2, VectorizedArray<value_type>>>
854 
855  virtual void
856  clear();
857 
864  std::shared_ptr<Table<2, VectorizedArray<value_type>>>
865  get_coefficient();
866 
867  private:
873  virtual void
874  apply_add(VectorType &dst, const VectorType &src) const;
875 
879  void
882  VectorType & dst,
883  const VectorType & src,
884  const std::pair<unsigned int, unsigned int> &cell_range) const;
885 
889  void
892  VectorType & dst,
893  const VectorType &,
894  const std::pair<unsigned int, unsigned int> &cell_range) const;
895 
899  void
902  & phi,
903  const unsigned int cell) const;
904 
908  std::shared_ptr<Table<2, VectorizedArray<value_type>>> scalar_coefficient;
909  };
910 
911 
912 
913  // ------------------------------------ inline functions ---------------------
914 
915  template <int dim, int fe_degree, int n_components, typename Number>
919  : fe_eval(fe_eval)
920  {
921  FullMatrix<double> shapes_1d(fe_degree + 1, fe_degree + 1);
922  for (unsigned int i = 0, c = 0; i < shapes_1d.m(); ++i)
923  for (unsigned int j = 0; j < shapes_1d.n(); ++j, ++c)
924  shapes_1d(i, j) = fe_eval.get_shape_info().shape_values[c][0];
925  shapes_1d.gauss_jordan();
926  const unsigned int stride = (fe_degree + 2) / 2;
927  inverse_shape.resize(stride * (fe_degree + 1));
928  for (unsigned int i = 0; i < stride; ++i)
929  for (unsigned int q = 0; q < (fe_degree + 2) / 2; ++q)
930  {
931  inverse_shape[i * stride + q] =
932  0.5 * (shapes_1d(i, q) + shapes_1d(i, fe_degree - q));
933  inverse_shape[(fe_degree - i) * stride + q] =
934  0.5 * (shapes_1d(i, q) - shapes_1d(i, fe_degree - q));
935  }
936  if (fe_degree % 2 == 0)
937  for (unsigned int q = 0; q < (fe_degree + 2) / 2; ++q)
938  inverse_shape[fe_degree / 2 * stride + q] = shapes_1d(fe_degree / 2, q);
939  }
940 
941 
942 
943  template <int dim, int fe_degree, int n_components, typename Number>
944  inline void
947  AlignedVector<VectorizedArray<Number>> &inverse_jxw) const
948  {
949  constexpr unsigned int dofs_per_component_on_cell =
950  Utilities::pow(fe_degree + 1, dim);
951  Assert(inverse_jxw.size() > 0 &&
952  inverse_jxw.size() % dofs_per_component_on_cell == 0,
953  ExcMessage(
954  "Expected diagonal to be a multiple of scalar dof per cells"));
955 
956  // temporarily reduce size of inverse_jxw to dofs_per_cell to get JxW values
957  // from fe_eval (will not reallocate any memory)
958  for (unsigned int q = 0; q < dofs_per_component_on_cell; ++q)
959  inverse_jxw[q] = 1. / fe_eval.JxW(q);
960  // copy values to rest of vector
961  for (unsigned int q = dofs_per_component_on_cell; q < inverse_jxw.size();)
962  for (unsigned int i = 0; i < dofs_per_component_on_cell; ++i, ++q)
963  inverse_jxw[q] = inverse_jxw[i];
964  }
965 
966 
967 
968  template <int dim, int fe_degree, int n_components, typename Number>
969  inline void
971  const AlignedVector<VectorizedArray<Number>> &inverse_coefficients,
972  const unsigned int n_actual_components,
973  const VectorizedArray<Number> * in_array,
974  VectorizedArray<Number> * out_array) const
975  {
976  constexpr unsigned int dofs_per_component =
977  Utilities::pow(fe_degree + 1, dim);
978  Assert(inverse_coefficients.size() > 0 &&
979  inverse_coefficients.size() % dofs_per_component == 0,
980  ExcMessage(
981  "Expected diagonal to be a multiple of scalar dof per cells"));
982  if (inverse_coefficients.size() != dofs_per_component)
983  AssertDimension(n_actual_components * dofs_per_component,
984  inverse_coefficients.size());
985 
986  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
987 
989  dim,
990  fe_degree + 1,
991  fe_degree + 1,
993  evaluator(inverse_shape, inverse_shape, inverse_shape);
994 
995  const unsigned int shift_coefficient =
996  inverse_coefficients.size() > dofs_per_component ? dofs_per_component : 0;
997  const VectorizedArray<Number> *inv_coefficient =
998  inverse_coefficients.data();
999  VectorizedArray<Number> temp_data_field[dofs_per_component];
1000  for (unsigned int d = 0; d < n_actual_components; ++d)
1001  {
1002  const VectorizedArray<Number> *in = in_array + d * dofs_per_component;
1003  VectorizedArray<Number> * out = out_array + d * dofs_per_component;
1004  // Need to select 'apply' method with hessian slot because values
1005  // assume symmetries that do not exist in the inverse shapes
1006  evaluator.template hessians<0, false, false>(in, temp_data_field);
1007  if (dim > 1)
1008  {
1009  evaluator.template hessians<1, false, false>(temp_data_field, out);
1010 
1011  if (dim == 3)
1012  {
1013  evaluator.template hessians<2, false, false>(out,
1014  temp_data_field);
1015  for (unsigned int q = 0; q < dofs_per_component; ++q)
1016  temp_data_field[q] *= inv_coefficient[q];
1017  evaluator.template hessians<2, true, false>(temp_data_field,
1018  out);
1019  }
1020  else if (dim == 2)
1021  for (unsigned int q = 0; q < dofs_per_component; ++q)
1022  out[q] *= inv_coefficient[q];
1023 
1024  evaluator.template hessians<1, true, false>(out, temp_data_field);
1025  }
1026  else
1027  {
1028  for (unsigned int q = 0; q < dofs_per_component; ++q)
1029  temp_data_field[q] *= inv_coefficient[q];
1030  }
1031  evaluator.template hessians<0, true, false>(temp_data_field, out);
1032 
1033  inv_coefficient += shift_coefficient;
1034  }
1035  }
1036 
1037 
1038 
1039  template <int dim, int fe_degree, int n_components, typename Number>
1040  inline void
1042  transform_from_q_points_to_basis(const unsigned int n_actual_components,
1043  const VectorizedArray<Number> *in_array,
1044  VectorizedArray<Number> *out_array) const
1045  {
1046  constexpr unsigned int dofs_per_cell = Utilities::pow(fe_degree + 1, dim);
1048  dim,
1049  fe_degree + 1,
1050  fe_degree + 1,
1054  inverse_shape);
1055 
1056  for (unsigned int d = 0; d < n_actual_components; ++d)
1057  {
1058  const VectorizedArray<Number> *in = in_array + d * dofs_per_cell;
1059  VectorizedArray<Number> * out = out_array + d * dofs_per_cell;
1060 
1061  if (dim == 3)
1062  {
1063  evaluator.template hessians<2, true, false>(in, out);
1064  evaluator.template hessians<1, true, false>(out, out);
1065  evaluator.template hessians<0, true, false>(out, out);
1066  }
1067  if (dim == 2)
1068  {
1069  evaluator.template hessians<1, true, false>(in, out);
1070  evaluator.template hessians<0, true, false>(out, out);
1071  }
1072  if (dim == 1)
1073  evaluator.template hessians<0, true, false>(in, out);
1074  }
1075  }
1076 
1077 
1078 
1079  //----------------- Base operator -----------------------------
1080  template <int dim, typename VectorType>
1082  : Subscriptor()
1083  , have_interface_matrices(false)
1084  {}
1085 
1086 
1087 
1088  template <int dim, typename VectorType>
1091  {
1092  Assert(data.get() != nullptr, ExcNotInitialized());
1093  typename Base<dim, VectorType>::size_type total_size = 0;
1094  for (unsigned int i = 0; i < selected_rows.size(); ++i)
1095  total_size += data->get_vector_partitioner(selected_rows[i])->size();
1096  return total_size;
1097  }
1098 
1099 
1100 
1101  template <int dim, typename VectorType>
1104  {
1105  Assert(data.get() != nullptr, ExcNotInitialized());
1106  typename Base<dim, VectorType>::size_type total_size = 0;
1107  for (unsigned int i = 0; i < selected_columns.size(); ++i)
1108  total_size += data->get_vector_partitioner(selected_columns[i])->size();
1109  return total_size;
1110  }
1111 
1112 
1113 
1114  template <int dim, typename VectorType>
1115  void
1117  {
1118  data.reset();
1119  inverse_diagonal_entries.reset();
1120  }
1121 
1122 
1123 
1124  template <int dim, typename VectorType>
1126  Base<dim, VectorType>::el(const unsigned int row,
1127  const unsigned int col) const
1128  {
1129  (void)col;
1130  Assert(row == col, ExcNotImplemented());
1131  Assert(inverse_diagonal_entries.get() != nullptr &&
1132  inverse_diagonal_entries->m() > 0,
1133  ExcNotInitialized());
1134  return 1.0 / (*inverse_diagonal_entries)(row, row);
1135  }
1136 
1137 
1138 
1139  template <int dim, typename VectorType>
1140  void
1142  {
1143  Assert(data.get() != nullptr, ExcNotInitialized());
1144  AssertDimension(BlockHelper::n_blocks(vec), selected_rows.size());
1145  for (unsigned int i = 0; i < BlockHelper::n_blocks(vec); ++i)
1146  {
1147  const unsigned int index = selected_rows[i];
1148  if (!BlockHelper::subblock(vec, index)
1149  .partitioners_are_compatible(
1150  *data->get_dof_info(index).vector_partitioner))
1151  data->initialize_dof_vector(BlockHelper::subblock(vec, index), index);
1152 
1153  Assert(BlockHelper::subblock(vec, index)
1154  .partitioners_are_globally_compatible(
1155  *data->get_dof_info(index).vector_partitioner),
1156  ExcInternalError());
1157  }
1158  BlockHelper::collect_sizes(vec);
1159  }
1160 
1161 
1162 
1163  template <int dim, typename VectorType>
1164  void
1166  std::shared_ptr<
1167  const MatrixFree<dim, typename Base<dim, VectorType>::value_type>> data_,
1168  const std::vector<unsigned int> &given_row_selection,
1169  const std::vector<unsigned int> &given_column_selection)
1170  {
1171  data = data_;
1172 
1173  selected_rows.clear();
1174  selected_columns.clear();
1175  if (given_row_selection.empty())
1176  for (unsigned int i = 0; i < data_->n_components(); ++i)
1177  selected_rows.push_back(i);
1178  else
1179  {
1180  for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1181  {
1182  AssertIndexRange(given_row_selection[i], data_->n_components());
1183  for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1184  if (j != i)
1185  Assert(given_row_selection[j] != given_row_selection[i],
1186  ExcMessage("Given row indices must be unique"));
1187 
1188  selected_rows.push_back(given_row_selection[i]);
1189  }
1190  }
1191  if (given_column_selection.size() == 0)
1192  selected_columns = selected_rows;
1193  else
1194  {
1195  for (unsigned int i = 0; i < given_column_selection.size(); ++i)
1196  {
1197  AssertIndexRange(given_column_selection[i], data_->n_components());
1198  for (unsigned int j = 0; j < given_column_selection.size(); ++j)
1199  if (j != i)
1200  Assert(given_column_selection[j] != given_column_selection[i],
1201  ExcMessage("Given column indices must be unique"));
1202 
1203  selected_columns.push_back(given_column_selection[i]);
1204  }
1205  }
1206 
1207  edge_constrained_indices.clear();
1208  edge_constrained_indices.resize(selected_rows.size());
1209  edge_constrained_values.clear();
1210  edge_constrained_values.resize(selected_rows.size());
1211  have_interface_matrices = false;
1212  }
1213 
1214 
1215 
1216  template <int dim, typename VectorType>
1217  void
1219  std::shared_ptr<
1220  const MatrixFree<dim, typename Base<dim, VectorType>::value_type>> data_,
1221  const MGConstrainedDoFs & mg_constrained_dofs,
1222  const unsigned int level,
1223  const std::vector<unsigned int> &given_row_selection)
1224  {
1225  std::vector<MGConstrainedDoFs> mg_constrained_dofs_vector(
1226  1, mg_constrained_dofs);
1227  initialize(data_, mg_constrained_dofs_vector, level, given_row_selection);
1228  }
1229 
1230 
1231 
1232  template <int dim, typename VectorType>
1233  void
1235  std::shared_ptr<const MatrixFree<dim, value_type>> data_,
1236  const std::vector<MGConstrainedDoFs> & mg_constrained_dofs,
1237  const unsigned int level,
1238  const std::vector<unsigned int> & given_row_selection)
1239  {
1241  ExcMessage("level is not set"));
1242 
1243  selected_rows.clear();
1244  selected_columns.clear();
1245  if (given_row_selection.empty())
1246  for (unsigned int i = 0; i < data_->n_components(); ++i)
1247  selected_rows.push_back(i);
1248  else
1249  {
1250  for (unsigned int i = 0; i < given_row_selection.size(); ++i)
1251  {
1252  AssertIndexRange(given_row_selection[i], data_->n_components());
1253  for (unsigned int j = 0; j < given_row_selection.size(); ++j)
1254  if (j != i)
1255  Assert(given_row_selection[j] != given_row_selection[i],
1256  ExcMessage("Given row indices must be unique"));
1257 
1258  selected_rows.push_back(given_row_selection[i]);
1259  }
1260  }
1261  selected_columns = selected_rows;
1262 
1263  AssertDimension(mg_constrained_dofs.size(), selected_rows.size());
1264  edge_constrained_indices.clear();
1265  edge_constrained_indices.resize(selected_rows.size());
1266  edge_constrained_values.clear();
1267  edge_constrained_values.resize(selected_rows.size());
1268 
1269  data = data_;
1270 
1271  for (unsigned int j = 0; j < selected_rows.size(); ++j)
1272  {
1273  if (data_->n_macro_cells() > 0)
1274  {
1275  AssertDimension(level, data_->get_cell_iterator(0, 0, j)->level());
1276  }
1277 
1278  // setup edge_constrained indices
1279  std::vector<types::global_dof_index> interface_indices;
1280  mg_constrained_dofs[j]
1281  .get_refinement_edge_indices(level)
1282  .fill_index_vector(interface_indices);
1283  edge_constrained_indices[j].clear();
1284  edge_constrained_indices[j].reserve(interface_indices.size());
1285  edge_constrained_values[j].resize(interface_indices.size());
1286  const IndexSet &locally_owned =
1287  data->get_dof_handler(selected_rows[j]).locally_owned_mg_dofs(level);
1288  for (const auto interface_index : interface_indices)
1289  if (locally_owned.is_element(interface_index))
1290  edge_constrained_indices[j].push_back(
1291  locally_owned.index_within_set(interface_index));
1292  have_interface_matrices |=
1294  static_cast<unsigned int>(edge_constrained_indices[j].size()),
1295  data->get_vector_partitioner()->get_mpi_communicator()) > 0;
1296  }
1297  }
1298 
1299 
1300 
1301  template <int dim, typename VectorType>
1302  void
1304  {
1305  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1306  {
1307  const std::vector<unsigned int> &constrained_dofs =
1308  data->get_constrained_dofs(selected_rows[j]);
1309  for (const auto constrained_dof : constrained_dofs)
1310  BlockHelper::subblock(dst, j).local_element(constrained_dof) = 1.;
1311  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1312  BlockHelper::subblock(dst, j).local_element(
1313  edge_constrained_indices[j][i]) = 1.;
1314  }
1315  }
1316 
1317 
1318 
1319  template <int dim, typename VectorType>
1320  void
1321  Base<dim, VectorType>::vmult(VectorType &dst, const VectorType &src) const
1322  {
1323  using Number = typename Base<dim, VectorType>::value_type;
1324  dst = Number(0.);
1325  vmult_add(dst, src);
1326  }
1327 
1328 
1329 
1330  template <int dim, typename VectorType>
1331  void
1332  Base<dim, VectorType>::vmult_add(VectorType &dst, const VectorType &src) const
1333  {
1334  mult_add(dst, src, false);
1335  }
1336 
1337 
1338 
1339  template <int dim, typename VectorType>
1340  void
1342  const VectorType &src) const
1343  {
1344  mult_add(dst, src, true);
1345  }
1346 
1347 
1348 
1349  template <int dim, typename VectorType>
1350  void
1352  const VectorType &src,
1353  const bool is_row) const
1354  {
1355  using Number = typename Base<dim, VectorType>::value_type;
1356  for (unsigned int i = 0; i < BlockHelper::n_blocks(src); ++i)
1357  {
1358  const unsigned int mf_component =
1359  is_row ? selected_rows[i] : selected_columns[i];
1360  // If both vectors use the same partitioner -> done
1361  if (BlockHelper::subblock(src, i).get_partitioner().get() ==
1362  data->get_dof_info(mf_component).vector_partitioner.get())
1363  continue;
1364 
1365  // If not, assert that the local ranges are the same and reset to the
1366  // current partitioner
1367  Assert(
1368  BlockHelper::subblock(src, i).get_partitioner()->local_size() ==
1369  data->get_dof_info(mf_component).vector_partitioner->local_size(),
1370  ExcMessage("The vector passed to the vmult() function does not have "
1371  "the correct size for compatibility with MatrixFree."));
1372 
1373  // copy the vector content to a temporary vector so that it does not get
1374  // lost
1376  BlockHelper::subblock(src, i));
1377  BlockHelper::subblock(const_cast<VectorType &>(src), i)
1378  .reinit(data->get_dof_info(mf_component).vector_partitioner);
1379  BlockHelper::subblock(const_cast<VectorType &>(src), i)
1380  .copy_locally_owned_data_from(copy_vec);
1381  }
1382  }
1383 
1384 
1385 
1386  template <int dim, typename VectorType>
1387  void
1389  const VectorType &src) const
1390  {
1391  using Number = typename Base<dim, VectorType>::value_type;
1392  adjust_ghost_range_if_necessary(src, false);
1393  adjust_ghost_range_if_necessary(dst, true);
1394 
1395  // set zero Dirichlet values on the input vector (and remember the src and
1396  // dst values because we need to reset them at the end)
1397  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1398  {
1399  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1400  {
1401  edge_constrained_values[j][i] = std::pair<Number, Number>(
1402  BlockHelper::subblock(src, j).local_element(
1403  edge_constrained_indices[j][i]),
1404  BlockHelper::subblock(dst, j).local_element(
1405  edge_constrained_indices[j][i]));
1406  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1407  .local_element(edge_constrained_indices[j][i]) = 0.;
1408  }
1409  }
1410  }
1411 
1412 
1413 
1414  template <int dim, typename VectorType>
1415  void
1417  const VectorType &src,
1418  const bool transpose) const
1419  {
1420  AssertDimension(dst.size(), src.size());
1421  AssertDimension(BlockHelper::n_blocks(dst), BlockHelper::n_blocks(src));
1422  AssertDimension(BlockHelper::n_blocks(dst), selected_rows.size());
1423  preprocess_constraints(dst, src);
1424  if (transpose)
1425  Tapply_add(dst, src);
1426  else
1427  apply_add(dst, src);
1428  postprocess_constraints(dst, src);
1429  }
1430 
1431 
1432 
1433  template <int dim, typename VectorType>
1434  void
1436  const VectorType &src) const
1437  {
1438  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1439  {
1440  const std::vector<unsigned int> &constrained_dofs =
1441  data->get_constrained_dofs(selected_rows[j]);
1442  for (const auto constrained_dof : constrained_dofs)
1443  BlockHelper::subblock(dst, j).local_element(constrained_dof) +=
1444  BlockHelper::subblock(src, j).local_element(constrained_dof);
1445  }
1446 
1447  // reset edge constrained values, multiply by unit matrix and add into
1448  // destination
1449  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1450  {
1451  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1452  {
1453  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1454  .local_element(edge_constrained_indices[j][i]) =
1455  edge_constrained_values[j][i].first;
1456  BlockHelper::subblock(dst, j).local_element(
1457  edge_constrained_indices[j][i]) =
1458  edge_constrained_values[j][i].second +
1459  edge_constrained_values[j][i].first;
1460  }
1461  }
1462  }
1463 
1464 
1465 
1466  template <int dim, typename VectorType>
1467  void
1469  const VectorType &src) const
1470  {
1471  using Number = typename Base<dim, VectorType>::value_type;
1472  AssertDimension(dst.size(), src.size());
1473  adjust_ghost_range_if_necessary(src, false);
1474  adjust_ghost_range_if_necessary(dst, true);
1475 
1476  dst = Number(0.);
1477 
1478  if (!have_interface_matrices)
1479  return;
1480 
1481  // set zero Dirichlet values on the input vector (and remember the src and
1482  // dst values because we need to reset them at the end)
1483  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1484  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1485  {
1486  edge_constrained_values[j][i] = std::pair<Number, Number>(
1487  BlockHelper::subblock(src, j).local_element(
1488  edge_constrained_indices[j][i]),
1489  BlockHelper::subblock(dst, j).local_element(
1490  edge_constrained_indices[j][i]));
1491  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1492  .local_element(edge_constrained_indices[j][i]) = 0.;
1493  }
1494 
1495  apply_add(dst, src);
1496 
1497  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1498  {
1499  unsigned int c = 0;
1500  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1501  {
1502  for (; c < edge_constrained_indices[j][i]; ++c)
1503  BlockHelper::subblock(dst, j).local_element(c) = 0.;
1504  ++c;
1505 
1506  // reset the src values
1507  BlockHelper::subblock(const_cast<VectorType &>(src), j)
1508  .local_element(edge_constrained_indices[j][i]) =
1509  edge_constrained_values[j][i].first;
1510  }
1511  for (; c < BlockHelper::subblock(dst, j).local_size(); ++c)
1512  BlockHelper::subblock(dst, j).local_element(c) = 0.;
1513  }
1514  }
1515 
1516 
1517 
1518  template <int dim, typename VectorType>
1519  void
1521  const VectorType &src) const
1522  {
1523  using Number = typename Base<dim, VectorType>::value_type;
1524  AssertDimension(dst.size(), src.size());
1525  adjust_ghost_range_if_necessary(src, false);
1526  adjust_ghost_range_if_necessary(dst, true);
1527 
1528  dst = Number(0.);
1529 
1530  if (!have_interface_matrices)
1531  return;
1532 
1533  VectorType src_cpy(src);
1534  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1535  {
1536  unsigned int c = 0;
1537  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1538  {
1539  for (; c < edge_constrained_indices[j][i]; ++c)
1540  BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1541  ++c;
1542  }
1543  for (; c < BlockHelper::subblock(src_cpy, j).local_size(); ++c)
1544  BlockHelper::subblock(src_cpy, j).local_element(c) = 0.;
1545  }
1546 
1547  apply_add(dst, src_cpy);
1548 
1549  for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j)
1550  for (unsigned int i = 0; i < edge_constrained_indices[j].size(); ++i)
1551  BlockHelper::subblock(dst, j).local_element(
1552  edge_constrained_indices[j][i]) = 0.;
1553  }
1554 
1555 
1556 
1557  template <int dim, typename VectorType>
1558  void
1559  Base<dim, VectorType>::Tvmult(VectorType &dst, const VectorType &src) const
1560  {
1561  using Number = typename Base<dim, VectorType>::value_type;
1562  dst = Number(0.);
1563  Tvmult_add(dst, src);
1564  }
1565 
1566 
1567 
1568  template <int dim, typename VectorType>
1569  std::size_t
1571  {
1572  return inverse_diagonal_entries.get() != nullptr ?
1573  inverse_diagonal_entries->memory_consumption() :
1574  sizeof(*this);
1575  }
1576 
1577 
1578 
1579  template <int dim, typename VectorType>
1580  std::shared_ptr<
1583  {
1584  return data;
1585  }
1586 
1587 
1588 
1589  template <int dim, typename VectorType>
1590  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1592  {
1593  Assert(inverse_diagonal_entries.get() != nullptr &&
1594  inverse_diagonal_entries->m() > 0,
1595  ExcNotInitialized());
1596  return inverse_diagonal_entries;
1597  }
1598 
1599 
1600 
1601  template <int dim, typename VectorType>
1602  const std::shared_ptr<DiagonalMatrix<VectorType>> &
1604  {
1605  Assert(diagonal_entries.get() != nullptr && diagonal_entries->m() > 0,
1606  ExcNotInitialized());
1607  return diagonal_entries;
1608  }
1609 
1610 
1611 
1612  template <int dim, typename VectorType>
1613  void
1615  const VectorType &src) const
1616  {
1617  apply_add(dst, src);
1618  }
1619 
1620 
1621 
1622  template <int dim, typename VectorType>
1623  void
1625  VectorType & dst,
1626  const VectorType & src,
1627  const typename Base<dim, VectorType>::value_type omega) const
1628  {
1629  Assert(inverse_diagonal_entries.get() && inverse_diagonal_entries->m() > 0,
1630  ExcNotInitialized());
1631  inverse_diagonal_entries->vmult(dst, src);
1632  dst *= omega;
1633  }
1634 
1635 
1636 
1637  //------------------------- MGInterfaceOperator ------------------------------
1638 
1639  template <typename OperatorType>
1641  : Subscriptor()
1642  , mf_base_operator(nullptr)
1643  {}
1644 
1645 
1646 
1647  template <typename OperatorType>
1648  void
1650  {
1651  mf_base_operator = nullptr;
1652  }
1653 
1654 
1655 
1656  template <typename OperatorType>
1657  void
1658  MGInterfaceOperator<OperatorType>::initialize(const OperatorType &operator_in)
1659  {
1660  mf_base_operator = &operator_in;
1661  }
1662 
1663 
1664 
1665  template <typename OperatorType>
1666  template <typename VectorType>
1667  void
1669  const VectorType &src) const
1670  {
1671 #ifndef DEAL_II_MSVC
1672  static_assert(
1673  std::is_same<typename VectorType::value_type, value_type>::value,
1674  "The vector type must be based on the same value type as this"
1675  "operator");
1676 #endif
1677 
1678  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1679 
1680  mf_base_operator->vmult_interface_down(dst, src);
1681  }
1682 
1683 
1684 
1685  template <typename OperatorType>
1686  template <typename VectorType>
1687  void
1689  const VectorType &src) const
1690  {
1691 #ifndef DEAL_II_MSVC
1692  static_assert(
1693  std::is_same<typename VectorType::value_type, value_type>::value,
1694  "The vector type must be based on the same value type as this"
1695  "operator");
1696 #endif
1697 
1698  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1699 
1700  mf_base_operator->vmult_interface_up(dst, src);
1701  }
1702 
1703 
1704 
1705  template <typename OperatorType>
1706  template <typename VectorType>
1707  void
1709  VectorType &vec) const
1710  {
1711  Assert(mf_base_operator != nullptr, ExcNotInitialized());
1712 
1713  mf_base_operator->initialize_dof_vector(vec);
1714  }
1715 
1716 
1717 
1718  //-----------------------------MassOperator----------------------------------
1719 
1720  template <int dim,
1721  int fe_degree,
1722  int n_q_points_1d,
1723  int n_components,
1724  typename VectorType>
1727  : Base<dim, VectorType>()
1728  {}
1729 
1730 
1731 
1732  template <int dim,
1733  int fe_degree,
1734  int n_q_points_1d,
1735  int n_components,
1736  typename VectorType>
1737  void
1740  {
1741  using Number = typename Base<dim, VectorType>::value_type;
1742  Assert((Base<dim, VectorType>::data.get() != nullptr), ExcNotInitialized());
1743 
1744  this->inverse_diagonal_entries.reset(new DiagonalMatrix<VectorType>());
1745  this->diagonal_entries.reset(new DiagonalMatrix<VectorType>());
1746  VectorType &inverse_diagonal_vector =
1747  this->inverse_diagonal_entries->get_vector();
1748  VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1749  this->initialize_dof_vector(inverse_diagonal_vector);
1750  this->initialize_dof_vector(diagonal_vector);
1751  inverse_diagonal_vector = Number(1.);
1752  apply_add(diagonal_vector, inverse_diagonal_vector);
1753 
1754  this->set_constrained_entries_to_one(diagonal_vector);
1755  inverse_diagonal_vector = diagonal_vector;
1756 
1757  const unsigned int local_size = inverse_diagonal_vector.local_size();
1758  for (unsigned int i = 0; i < local_size; ++i)
1759  inverse_diagonal_vector.local_element(i) =
1760  Number(1.) / inverse_diagonal_vector.local_element(i);
1761 
1762  inverse_diagonal_vector.update_ghost_values();
1763  diagonal_vector.update_ghost_values();
1764  }
1765 
1766 
1767 
1768  template <int dim,
1769  int fe_degree,
1770  int n_q_points_1d,
1771  int n_components,
1772  typename VectorType>
1773  void
1775  apply_add(VectorType &dst, const VectorType &src) const
1776  {
1778  this,
1779  dst,
1780  src);
1781  }
1782 
1783 
1784 
1785  template <int dim,
1786  int fe_degree,
1787  int n_q_points_1d,
1788  int n_components,
1789  typename VectorType>
1790  void
1793  const MatrixFree<dim, typename Base<dim, VectorType>::value_type> &data,
1794  VectorType & dst,
1795  const VectorType & src,
1796  const std::pair<unsigned int, unsigned int> &cell_range) const
1797  {
1798  using Number = typename Base<dim, VectorType>::value_type;
1800  data, this->selected_rows[0]);
1801  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
1802  {
1803  phi.reinit(cell);
1804  phi.read_dof_values(src);
1805  phi.evaluate(true, false, false);
1806  for (unsigned int q = 0; q < phi.n_q_points; ++q)
1807  phi.submit_value(phi.get_value(q), q);
1808  phi.integrate(true, false);
1809  phi.distribute_local_to_global(dst);
1810  }
1811  }
1812 
1813 
1814  //-----------------------------LaplaceOperator----------------------------------
1815 
1816  template <int dim,
1817  int fe_degree,
1818  int n_q_points_1d,
1819  int n_components,
1820  typename VectorType>
1823  : Base<dim, VectorType>()
1824  {}
1825 
1826 
1827 
1828  template <int dim,
1829  int fe_degree,
1830  int n_q_points_1d,
1831  int n_components,
1832  typename VectorType>
1833  void
1836  {
1838  scalar_coefficient.reset();
1839  }
1840 
1841 
1842 
1843  template <int dim,
1844  int fe_degree,
1845  int n_q_points_1d,
1846  int n_components,
1847  typename VectorType>
1848  void
1851  const std::shared_ptr<
1853  &scalar_coefficient_)
1854  {
1855  scalar_coefficient = scalar_coefficient_;
1856  }
1857 
1858 
1859 
1860  template <int dim,
1861  int fe_degree,
1862  int n_q_points_1d,
1863  int n_components,
1864  typename VectorType>
1865  std::shared_ptr<
1866  Table<2,
1867  VectorizedArray<typename LaplaceOperator<dim,
1868  fe_degree,
1869  n_q_points_1d,
1870  n_components,
1871  VectorType>::value_type>>>
1874  {
1875  Assert(scalar_coefficient.get(), ExcNotInitialized());
1876  return scalar_coefficient;
1877  }
1878 
1879 
1880 
1881  template <int dim,
1882  int fe_degree,
1883  int n_q_points_1d,
1884  int n_components,
1885  typename VectorType>
1886  void
1889  {
1890  using Number = typename Base<dim, VectorType>::value_type;
1891  Assert((Base<dim, VectorType>::data.get() != nullptr), ExcNotInitialized());
1892 
1893  this->inverse_diagonal_entries.reset(new DiagonalMatrix<VectorType>());
1894  this->diagonal_entries.reset(new DiagonalMatrix<VectorType>());
1895  VectorType &inverse_diagonal_vector =
1896  this->inverse_diagonal_entries->get_vector();
1897  VectorType &diagonal_vector = this->diagonal_entries->get_vector();
1898  this->initialize_dof_vector(inverse_diagonal_vector);
1899  this->initialize_dof_vector(diagonal_vector);
1900 
1901  this->data->cell_loop(&LaplaceOperator::local_diagonal_cell,
1902  this,
1903  diagonal_vector,
1904  /*unused*/ diagonal_vector);
1905  this->set_constrained_entries_to_one(diagonal_vector);
1906 
1907  inverse_diagonal_vector = diagonal_vector;
1908 
1909  for (unsigned int i = 0; i < inverse_diagonal_vector.local_size(); ++i)
1910  if (std::abs(inverse_diagonal_vector.local_element(i)) >
1911  std::sqrt(std::numeric_limits<Number>::epsilon()))
1912  inverse_diagonal_vector.local_element(i) =
1913  1. / inverse_diagonal_vector.local_element(i);
1914  else
1915  inverse_diagonal_vector.local_element(i) = 1.;
1916 
1917  inverse_diagonal_vector.update_ghost_values();
1918  diagonal_vector.update_ghost_values();
1919  }
1920 
1921 
1922 
1923  template <int dim,
1924  int fe_degree,
1925  int n_q_points_1d,
1926  int n_components,
1927  typename VectorType>
1928  void
1930  apply_add(VectorType &dst, const VectorType &src) const
1931  {
1933  this,
1934  dst,
1935  src);
1936  }
1937 
1938  namespace Implementation
1939  {
1940  template <typename Number>
1941  bool
1942  non_negative(const VectorizedArray<Number> &n)
1943  {
1944  for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
1945  ++v)
1946  if (n[v] < 0.)
1947  return false;
1948 
1949  return true;
1950  }
1951  } // namespace Implementation
1952 
1953 
1954 
1955  template <int dim,
1956  int fe_degree,
1957  int n_q_points_1d,
1958  int n_components,
1959  typename VectorType>
1960  void
1963  FEEvaluation<dim,
1964  fe_degree,
1965  n_q_points_1d,
1966  n_components,
1967  typename Base<dim, VectorType>::value_type> &phi,
1968  const unsigned int cell) const
1969  {
1970  phi.evaluate(false, true, false);
1971  if (scalar_coefficient.get())
1972  {
1973  for (unsigned int q = 0; q < phi.n_q_points; ++q)
1974  {
1975  Assert(Implementation::non_negative((*scalar_coefficient)(cell, q)),
1976  ExcMessage("Coefficient must be non-negative"));
1977  phi.submit_gradient((*scalar_coefficient)(cell, q) *
1978  phi.get_gradient(q),
1979  q);
1980  }
1981  }
1982  else
1983  {
1984  for (unsigned int q = 0; q < phi.n_q_points; ++q)
1985  {
1986  phi.submit_gradient(phi.get_gradient(q), q);
1987  }
1988  }
1989  phi.integrate(false, true);
1990  }
1991 
1992 
1993 
1994  template <int dim,
1995  int fe_degree,
1996  int n_q_points_1d,
1997  int n_components,
1998  typename VectorType>
1999  void
2002  const MatrixFree<dim, typename Base<dim, VectorType>::value_type> &data,
2003  VectorType & dst,
2004  const VectorType & src,
2005  const std::pair<unsigned int, unsigned int> &cell_range) const
2006  {
2007  using Number = typename Base<dim, VectorType>::value_type;
2009  data, this->selected_rows[0]);
2010  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2011  {
2012  phi.reinit(cell);
2013  phi.read_dof_values(src);
2014  do_operation_on_cell(phi, cell);
2015  phi.distribute_local_to_global(dst);
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  void
2028  const MatrixFree<dim, typename Base<dim, VectorType>::value_type> &data,
2029  VectorType & dst,
2030  const VectorType &,
2031  const std::pair<unsigned int, unsigned int> &cell_range) const
2032  {
2033  using Number = typename Base<dim, VectorType>::value_type;
2035  data, this->selected_rows[0]);
2036  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2037  {
2038  phi.reinit(cell);
2039  VectorizedArray<Number> local_diagonal_vector[phi.static_dofs_per_cell];
2040  for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2041  {
2042  for (unsigned int j = 0; j < phi.dofs_per_component; ++j)
2044  phi.begin_dof_values()[i] = 1.;
2045  do_operation_on_cell(phi, cell);
2046  local_diagonal_vector[i] = phi.begin_dof_values()[i];
2047  }
2048  for (unsigned int i = 0; i < phi.dofs_per_component; ++i)
2049  for (unsigned int c = 0; c < phi.n_components; ++c)
2050  phi.begin_dof_values()[i + c * phi.dofs_per_component] =
2051  local_diagonal_vector[i];
2052  phi.distribute_local_to_global(dst);
2053  }
2054  }
2055 
2056 
2057 } // end of namespace MatrixFreeOperators
2058 
2059 
2060 DEAL_II_NAMESPACE_CLOSE
2061 
2062 #endif
typename OperatorType::value_type value_type
Definition: operators.h:537
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:1962
static const unsigned int invalid_unsigned_int
Definition: types.h:173
void vmult_interface_down(VectorType &dst, const VectorType &src) const
Definition: operators.h:1468
std::vector< unsigned int > selected_rows
Definition: operators.h:445
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
virtual void apply_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1930
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:2001
typename VectorType::size_type size_type
Definition: operators.h:198
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal_inverse() const
Definition: operators.h:1591
void gauss_jordan()
void evaluate(const bool evaluate_values, const bool evaluate_gradients, const bool evaluate_hessians=false)
void integrate(const bool integrate_values, const bool integrate_gradients)
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:1624
void fill_inverse_JxW_values(AlignedVector< VectorizedArray< Number >> &inverse_jxw) const
Definition: operators.h:946
static constexpr unsigned int static_dofs_per_cell
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
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:1416
const VectorizedArray< Number > * begin_dof_values() const
const std::shared_ptr< DiagonalMatrix< VectorType > > & get_matrix_diagonal() const
Definition: operators.h:1603
static ::ExceptionBase & ExcNotInitialized()
void vmult_interface_up(VectorType &dst, const VectorType &src) const
Definition: operators.h:1520
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:1303
std::shared_ptr< DiagonalMatrix< VectorType > > diagonal_entries
Definition: operators.h:433
void read_dof_values(const VectorType &src, const unsigned int first_index=0)
static constexpr unsigned int n_components
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
const unsigned int dofs_per_component
value_type get_value(const unsigned int q_point) const
void submit_value(const value_type val_in, const unsigned int q_point)
const FEEvaluationBase< dim, n_components, Number > & fe_eval
Definition: operators.h:690
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1708
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:1126
static ::ExceptionBase & ExcMessage(std::string arg1)
const unsigned int n_q_points
size_type n() const
Definition: operators.h:1103
size_type n() const
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArray< Number > > & get_shape_info() const
virtual void Tapply_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1614
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1559
#define Assert(cond, exc)
Definition: exceptions.h:1407
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1821
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:2027
constexpr unsigned int pow(const unsigned int base, const int iexp)
Definition: utilities.h:428
void set_coefficient(const std::shared_ptr< Table< 2, VectorizedArray< value_type >>> &scalar_coefficient)
Definition: operators.h:1850
virtual void compute_diagonal()=0
void apply(const AlignedVector< VectorizedArray< Number >> &inverse_coefficient, const unsigned int n_actual_components, const VectorizedArray< Number > *in_array, VectorizedArray< Number > *out_array) const
Definition: operators.h:970
std::vector< std::vector< unsigned int > > edge_constrained_indices
Definition: operators.h:457
AlignedVector< VectorizedArray< Number > > inverse_shape
Definition: operators.h:695
SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
void Tvmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1341
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1321
size_type m() const
Definition: operators.h:1090
void transform_from_q_points_to_basis(const unsigned int n_actual_components, const VectorizedArray< Number > *in_array, VectorizedArray< Number > *out_array) const
Definition: operators.h:1042
void reinit(const unsigned int cell_batch_index)
virtual void apply_add(VectorType &dst, const VectorType &src) const =0
std::shared_ptr< Table< 2, VectorizedArray< value_type > > > scalar_coefficient
Definition: operators.h:908
void vmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1668
void distribute_local_to_global(VectorType &dst, const unsigned int first_index=0, const std::bitset< VectorizedArray< Number >::n_array_elements > &mask=std::bitset< VectorizedArray< Number >::n_array_elements >().flip()) const
virtual void compute_diagonal() override
Definition: operators.h:1739
virtual ~Base() override=default
void preprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1388
CellwiseInverseMassMatrix(const FEEvaluationBase< dim, n_components, Number > &fe_eval)
Definition: operators.h:917
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:1792
std::shared_ptr< const MatrixFree< dim, value_type > > get_matrix_free() const
Definition: operators.h:1582
std::shared_ptr< Table< 2, VectorizedArray< value_type > > > get_coefficient()
Definition: operators.h:1873
virtual std::size_t memory_consumption() const
Definition: operators.h:1570
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:1775
void postprocess_constraints(VectorType &dst, const VectorType &src) const
Definition: operators.h:1435
static ::ExceptionBase & ExcNotImplemented()
Definition: table.h:37
bool is_element(const size_type index) const
Definition: index_set.h:1665
void adjust_ghost_range_if_necessary(const VectorType &vec, const bool is_row) const
Definition: operators.h:1351
void vmult_add(VectorType &dst, const VectorType &src) const
Definition: operators.h:1332
void initialize(const OperatorType &operator_in)
Definition: operators.h:1658
virtual void clear()
Definition: operators.h:1116
SmartPointer< const OperatorType > mf_base_operator
Definition: operators.h:587
void Tvmult(VectorType &dst, const VectorType &src) const
Definition: operators.h:1688
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
void initialize_dof_vector(VectorType &vec) const
Definition: operators.h:1141
T max(const T &t, const MPI_Comm &mpi_communicator)
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 >())
static ::ExceptionBase & ExcInternalError()