Reference documentation for deal.II version GIT 7b2de2f2f9 2023-09-24 11:00:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
matrix_free.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2011 - 2023 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_h
18 #define dealii_matrix_free_h
19 
20 #include <deal.II/base/config.h>
21 
28 
30 
31 #include <deal.II/fe/fe.h>
32 #include <deal.II/fe/mapping.h>
33 
36 
41 
48 
49 #include <cstdlib>
50 #include <limits>
51 #include <list>
52 #include <memory>
53 
54 
56 
57 
58 
110 template <int dim,
111  typename Number = double,
112  typename VectorizedArrayType = VectorizedArray<Number>>
113 class MatrixFree : public Subscriptor
114 {
115  static_assert(
116  std::is_same_v<Number, typename VectorizedArrayType::value_type>,
117  "Type of Number and of VectorizedArrayType do not match.");
118 
119 public:
124  using value_type = Number;
125  using vectorized_value_type = VectorizedArrayType;
126 
130  static constexpr unsigned int dimension = dim;
131 
185  {
190 
197  {
217  };
218 
224  const unsigned int tasks_block_size = 0,
230  const unsigned int mg_level = numbers::invalid_unsigned_int,
231  const bool store_plain_indices = true,
232  const bool initialize_indices = true,
233  const bool initialize_mapping = true,
234  const bool overlap_communication_computation = true,
235  const bool hold_all_faces_to_owned_cells = false,
236  const bool cell_vectorization_categories_strict = false,
237  const bool allow_ghosted_vectors_in_loops = true)
244  , mg_level(mg_level)
253  , communicator_sm(MPI_COMM_SELF)
254  {}
255 
268  , mg_level(other.mg_level)
280  {}
281 
286  operator=(const AdditionalData &other)
287  {
296  mg_level = other.mg_level;
308 
309  return *this;
310  }
311 
349 
359  unsigned int tasks_block_size;
360 
375 
395 
415 
443 
452  unsigned int mg_level;
453 
461 
472 
481 
495 
504 
526  std::vector<unsigned int> cell_vectorization_category;
527 
538 
550 
554  MPI_Comm communicator_sm;
555  };
556 
565 
570 
574  ~MatrixFree() override = default;
575 
588  template <typename QuadratureType, typename number2, typename MappingType>
589  void
590  reinit(const MappingType &mapping,
591  const DoFHandler<dim> &dof_handler,
592  const AffineConstraints<number2> &constraint,
593  const QuadratureType &quad,
594  const AdditionalData &additional_data = AdditionalData());
595 
617  template <typename QuadratureType, typename number2, typename MappingType>
618  void
619  reinit(const MappingType &mapping,
620  const std::vector<const DoFHandler<dim> *> &dof_handler,
621  const std::vector<const AffineConstraints<number2> *> &constraint,
622  const std::vector<QuadratureType> &quad,
623  const AdditionalData &additional_data = AdditionalData());
624 
632  template <typename QuadratureType, typename number2, typename MappingType>
633  void
634  reinit(const MappingType &mapping,
635  const std::vector<const DoFHandler<dim> *> &dof_handler,
636  const std::vector<const AffineConstraints<number2> *> &constraint,
637  const QuadratureType &quad,
638  const AdditionalData &additional_data = AdditionalData());
639 
645  void
647  const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free_base);
648 
658  void
659  update_mapping(const Mapping<dim> &mapping);
660 
664  void
665  update_mapping(const std::shared_ptr<hp::MappingCollection<dim>> &mapping);
666 
671  void
672  clear();
673 
687  enum class DataAccessOnFaces
688  {
695  none,
696 
707  values,
708 
718 
729  gradients,
730 
740 
748  };
749 
794  template <typename OutVector, typename InVector>
795  void
796  cell_loop(const std::function<void(
798  OutVector &,
799  const InVector &,
800  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
801  OutVector &dst,
802  const InVector &src,
803  const bool zero_dst_vector = false) const;
804 
851  template <typename CLASS, typename OutVector, typename InVector>
852  void
853  cell_loop(void (CLASS::*cell_operation)(
854  const MatrixFree &,
855  OutVector &,
856  const InVector &,
857  const std::pair<unsigned int, unsigned int> &) const,
858  const CLASS *owning_class,
859  OutVector &dst,
860  const InVector &src,
861  const bool zero_dst_vector = false) const;
862 
866  template <typename CLASS, typename OutVector, typename InVector>
867  void
868  cell_loop(void (CLASS::*cell_operation)(
869  const MatrixFree &,
870  OutVector &,
871  const InVector &,
872  const std::pair<unsigned int, unsigned int> &),
873  CLASS *owning_class,
874  OutVector &dst,
875  const InVector &src,
876  const bool zero_dst_vector = false) const;
877 
962  template <typename CLASS, typename OutVector, typename InVector>
963  void
964  cell_loop(void (CLASS::*cell_operation)(
965  const MatrixFree &,
966  OutVector &,
967  const InVector &,
968  const std::pair<unsigned int, unsigned int> &) const,
969  const CLASS *owning_class,
970  OutVector &dst,
971  const InVector &src,
972  const std::function<void(const unsigned int, const unsigned int)>
973  &operation_before_loop,
974  const std::function<void(const unsigned int, const unsigned int)>
975  &operation_after_loop,
976  const unsigned int dof_handler_index_pre_post = 0) const;
977 
981  template <typename CLASS, typename OutVector, typename InVector>
982  void
983  cell_loop(void (CLASS::*cell_operation)(
984  const MatrixFree &,
985  OutVector &,
986  const InVector &,
987  const std::pair<unsigned int, unsigned int> &),
988  CLASS *owning_class,
989  OutVector &dst,
990  const InVector &src,
991  const std::function<void(const unsigned int, const unsigned int)>
992  &operation_before_loop,
993  const std::function<void(const unsigned int, const unsigned int)>
994  &operation_after_loop,
995  const unsigned int dof_handler_index_pre_post = 0) const;
996 
1001  template <typename OutVector, typename InVector>
1002  void
1003  cell_loop(const std::function<void(
1005  OutVector &,
1006  const InVector &,
1007  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1008  OutVector &dst,
1009  const InVector &src,
1010  const std::function<void(const unsigned int, const unsigned int)>
1011  &operation_before_loop,
1012  const std::function<void(const unsigned int, const unsigned int)>
1013  &operation_after_loop,
1014  const unsigned int dof_handler_index_pre_post = 0) const;
1015 
1091  template <typename OutVector, typename InVector>
1092  void
1093  loop(const std::function<
1095  OutVector &,
1096  const InVector &,
1097  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1098  const std::function<
1100  OutVector &,
1101  const InVector &,
1102  const std::pair<unsigned int, unsigned int> &)> &face_operation,
1103  const std::function<void(
1105  OutVector &,
1106  const InVector &,
1107  const std::pair<unsigned int, unsigned int> &)> &boundary_operation,
1108  OutVector &dst,
1109  const InVector &src,
1110  const bool zero_dst_vector = false,
1111  const DataAccessOnFaces dst_vector_face_access =
1113  const DataAccessOnFaces src_vector_face_access =
1115 
1201  template <typename CLASS, typename OutVector, typename InVector>
1202  void
1204  void (CLASS::*cell_operation)(const MatrixFree &,
1205  OutVector &,
1206  const InVector &,
1207  const std::pair<unsigned int, unsigned int> &)
1208  const,
1209  void (CLASS::*face_operation)(const MatrixFree &,
1210  OutVector &,
1211  const InVector &,
1212  const std::pair<unsigned int, unsigned int> &)
1213  const,
1214  void (CLASS::*boundary_operation)(
1215  const MatrixFree &,
1216  OutVector &,
1217  const InVector &,
1218  const std::pair<unsigned int, unsigned int> &) const,
1219  const CLASS *owning_class,
1220  OutVector &dst,
1221  const InVector &src,
1222  const bool zero_dst_vector = false,
1223  const DataAccessOnFaces dst_vector_face_access =
1225  const DataAccessOnFaces src_vector_face_access =
1227 
1231  template <typename CLASS, typename OutVector, typename InVector>
1232  void
1233  loop(void (CLASS::*cell_operation)(
1234  const MatrixFree &,
1235  OutVector &,
1236  const InVector &,
1237  const std::pair<unsigned int, unsigned int> &),
1238  void (CLASS::*face_operation)(
1239  const MatrixFree &,
1240  OutVector &,
1241  const InVector &,
1242  const std::pair<unsigned int, unsigned int> &),
1243  void (CLASS::*boundary_operation)(
1244  const MatrixFree &,
1245  OutVector &,
1246  const InVector &,
1247  const std::pair<unsigned int, unsigned int> &),
1248  CLASS *owning_class,
1249  OutVector &dst,
1250  const InVector &src,
1251  const bool zero_dst_vector = false,
1252  const DataAccessOnFaces dst_vector_face_access =
1254  const DataAccessOnFaces src_vector_face_access =
1256 
1376  template <typename CLASS, typename OutVector, typename InVector>
1377  void
1379  void (CLASS::*cell_operation)(const MatrixFree &,
1380  OutVector &,
1381  const InVector &,
1382  const std::pair<unsigned int, unsigned int> &)
1383  const,
1384  void (CLASS::*face_operation)(const MatrixFree &,
1385  OutVector &,
1386  const InVector &,
1387  const std::pair<unsigned int, unsigned int> &)
1388  const,
1389  void (CLASS::*boundary_operation)(
1390  const MatrixFree &,
1391  OutVector &,
1392  const InVector &,
1393  const std::pair<unsigned int, unsigned int> &) const,
1394  const CLASS *owning_class,
1395  OutVector &dst,
1396  const InVector &src,
1397  const std::function<void(const unsigned int, const unsigned int)>
1398  &operation_before_loop,
1399  const std::function<void(const unsigned int, const unsigned int)>
1400  &operation_after_loop,
1401  const unsigned int dof_handler_index_pre_post = 0,
1402  const DataAccessOnFaces dst_vector_face_access =
1404  const DataAccessOnFaces src_vector_face_access =
1406 
1410  template <typename CLASS, typename OutVector, typename InVector>
1411  void
1412  loop(void (CLASS::*cell_operation)(
1413  const MatrixFree &,
1414  OutVector &,
1415  const InVector &,
1416  const std::pair<unsigned int, unsigned int> &),
1417  void (CLASS::*face_operation)(
1418  const MatrixFree &,
1419  OutVector &,
1420  const InVector &,
1421  const std::pair<unsigned int, unsigned int> &),
1422  void (CLASS::*boundary_operation)(
1423  const MatrixFree &,
1424  OutVector &,
1425  const InVector &,
1426  const std::pair<unsigned int, unsigned int> &),
1427  const CLASS *owning_class,
1428  OutVector &dst,
1429  const InVector &src,
1430  const std::function<void(const unsigned int, const unsigned int)>
1431  &operation_before_loop,
1432  const std::function<void(const unsigned int, const unsigned int)>
1433  &operation_after_loop,
1434  const unsigned int dof_handler_index_pre_post = 0,
1435  const DataAccessOnFaces dst_vector_face_access =
1437  const DataAccessOnFaces src_vector_face_access =
1439 
1445  template <typename OutVector, typename InVector>
1446  void
1447  loop(const std::function<
1449  OutVector &,
1450  const InVector &,
1451  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1452  const std::function<
1454  OutVector &,
1455  const InVector &,
1456  const std::pair<unsigned int, unsigned int> &)> &face_operation,
1457  const std::function<void(
1459  OutVector &,
1460  const InVector &,
1461  const std::pair<unsigned int, unsigned int> &)> &boundary_operation,
1462  OutVector &dst,
1463  const InVector &src,
1464  const std::function<void(const unsigned int, const unsigned int)>
1465  &operation_before_loop,
1466  const std::function<void(const unsigned int, const unsigned int)>
1467  &operation_after_loop,
1468  const unsigned int dof_handler_index_pre_post = 0,
1469  const DataAccessOnFaces dst_vector_face_access =
1471  const DataAccessOnFaces src_vector_face_access =
1473 
1538  template <typename CLASS, typename OutVector, typename InVector>
1539  void
1540  loop_cell_centric(void (CLASS::*cell_operation)(
1541  const MatrixFree &,
1542  OutVector &,
1543  const InVector &,
1544  const std::pair<unsigned int, unsigned int> &) const,
1545  const CLASS *owning_class,
1546  OutVector &dst,
1547  const InVector &src,
1548  const bool zero_dst_vector = false,
1549  const DataAccessOnFaces src_vector_face_access =
1551 
1555  template <typename CLASS, typename OutVector, typename InVector>
1556  void
1557  loop_cell_centric(void (CLASS::*cell_operation)(
1558  const MatrixFree &,
1559  OutVector &,
1560  const InVector &,
1561  const std::pair<unsigned int, unsigned int> &),
1562  CLASS *owning_class,
1563  OutVector &dst,
1564  const InVector &src,
1565  const bool zero_dst_vector = false,
1566  const DataAccessOnFaces src_vector_face_access =
1568 
1572  template <typename OutVector, typename InVector>
1573  void
1575  const std::function<void(const MatrixFree &,
1576  OutVector &,
1577  const InVector &,
1578  const std::pair<unsigned int, unsigned int> &)>
1579  &cell_operation,
1580  OutVector &dst,
1581  const InVector &src,
1582  const bool zero_dst_vector = false,
1583  const DataAccessOnFaces src_vector_face_access =
1585 
1593  std::pair<unsigned int, unsigned int>
1594  create_cell_subrange_hp(const std::pair<unsigned int, unsigned int> &range,
1595  const unsigned int fe_degree,
1596  const unsigned int dof_handler_index = 0) const;
1597 
1604  std::pair<unsigned int, unsigned int>
1606  const std::pair<unsigned int, unsigned int> &range,
1607  const unsigned int fe_index,
1608  const unsigned int dof_handler_index = 0) const;
1609 
1613  unsigned int
1615 
1619  unsigned int
1621  const std::pair<unsigned int, unsigned int> range) const;
1622 
1626  unsigned int
1627  get_face_active_fe_index(const std::pair<unsigned int, unsigned int> range,
1628  const bool is_interior_face = true) const;
1629 
1641  template <typename T>
1642  void
1644 
1650  template <typename T>
1651  void
1653 
1667  template <typename VectorType>
1668  void
1669  initialize_dof_vector(VectorType &vec,
1670  const unsigned int dof_handler_index = 0) const;
1671 
1689  template <typename Number2>
1690  void
1692  const unsigned int dof_handler_index = 0) const;
1693 
1704  const std::shared_ptr<const Utilities::MPI::Partitioner> &
1705  get_vector_partitioner(const unsigned int dof_handler_index = 0) const;
1706 
1710  const IndexSet &
1711  get_locally_owned_set(const unsigned int dof_handler_index = 0) const;
1712 
1716  const IndexSet &
1717  get_ghost_set(const unsigned int dof_handler_index = 0) const;
1718 
1728  const std::vector<unsigned int> &
1729  get_constrained_dofs(const unsigned int dof_handler_index = 0) const;
1730 
1741  void
1742  renumber_dofs(std::vector<types::global_dof_index> &renumbering,
1743  const unsigned int dof_handler_index = 0);
1744 
1754  template <int spacedim>
1755  static bool
1757 
1761  unsigned int
1762  n_components() const;
1763 
1768  unsigned int
1769  n_base_elements(const unsigned int dof_handler_index) const;
1770 
1778  unsigned int
1780 
1790  unsigned int
1792 
1799  unsigned int
1801 
1811  unsigned int
1813 
1824  unsigned int
1826 
1831  unsigned int
1833 
1845  get_boundary_id(const unsigned int face_batch_index) const;
1846 
1851  std::array<types::boundary_id, VectorizedArrayType::size()>
1852  get_faces_by_cells_boundary_id(const unsigned int cell_batch_index,
1853  const unsigned int face_number) const;
1854 
1859  const DoFHandler<dim> &
1860  get_dof_handler(const unsigned int dof_handler_index = 0) const;
1861 
1869  get_affine_constraints(const unsigned int dof_handler_index = 0) const;
1870 
1884  get_cell_iterator(const unsigned int cell_batch_index,
1885  const unsigned int lane_index,
1886  const unsigned int dof_handler_index = 0) const;
1887 
1893  std::pair<int, int>
1894  get_cell_level_and_index(const unsigned int cell_batch_index,
1895  const unsigned int lane_index) const;
1896 
1903  unsigned int
1905  const typename Triangulation<dim>::cell_iterator &cell) const;
1906 
1919  std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>
1920  get_face_iterator(const unsigned int face_batch_index,
1921  const unsigned int lane_index,
1922  const bool interior = true,
1923  const unsigned int fe_component = 0) const;
1924 
1937  bool
1938  at_irregular_cell(const unsigned int cell_batch_index) const;
1939 
1949  unsigned int
1950  n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const;
1951 
1961  unsigned int
1962  n_active_entries_per_face_batch(const unsigned int face_batch_index) const;
1963 
1967  unsigned int
1968  get_dofs_per_cell(const unsigned int dof_handler_index = 0,
1969  const unsigned int hp_active_fe_index = 0) const;
1970 
1974  unsigned int
1975  get_n_q_points(const unsigned int quad_index = 0,
1976  const unsigned int hp_active_fe_index = 0) const;
1977 
1982  unsigned int
1983  get_dofs_per_face(const unsigned int dof_handler_index = 0,
1984  const unsigned int hp_active_fe_index = 0) const;
1985 
1990  unsigned int
1991  get_n_q_points_face(const unsigned int quad_index = 0,
1992  const unsigned int hp_active_fe_index = 0) const;
1993 
1997  const Quadrature<dim> &
1998  get_quadrature(const unsigned int quad_index = 0,
1999  const unsigned int hp_active_fe_index = 0) const;
2000 
2004  const Quadrature<dim - 1> &
2005  get_face_quadrature(const unsigned int quad_index = 0,
2006  const unsigned int hp_active_fe_index = 0) const;
2007 
2020  unsigned int
2022  const std::pair<unsigned int, unsigned int> cell_batch_range) const;
2023 
2028  std::pair<unsigned int, unsigned int>
2030  const std::pair<unsigned int, unsigned int> face_batch_range) const;
2031 
2043  unsigned int
2044  get_cell_category(const unsigned int cell_batch_index) const;
2045 
2050  std::pair<unsigned int, unsigned int>
2051  get_face_category(const unsigned int face_batch_index) const;
2052 
2056  bool
2058 
2063  bool
2065 
2070  unsigned int
2071  get_mg_level() const;
2072 
2077  std::size_t
2079 
2084  template <typename StreamType>
2085  void
2086  print_memory_consumption(StreamType &out) const;
2087 
2092  void
2093  print(std::ostream &out) const;
2094 
2107  get_task_info() const;
2108 
2109  /*
2110  * Return geometry-dependent information on the cells.
2111  */
2112  const internal::MatrixFreeFunctions::
2113  MappingInfo<dim, Number, VectorizedArrayType> &
2115 
2120  get_dof_info(const unsigned int dof_handler_index_component = 0) const;
2121 
2125  unsigned int
2127 
2132  const Number *
2133  constraint_pool_begin(const unsigned int pool_index) const;
2134 
2140  const Number *
2141  constraint_pool_end(const unsigned int pool_index) const;
2142 
2147  get_shape_info(const unsigned int dof_handler_index_component = 0,
2148  const unsigned int quad_index = 0,
2149  const unsigned int fe_base_element = 0,
2150  const unsigned int hp_active_fe_index = 0,
2151  const unsigned int hp_active_quad_index = 0) const;
2152 
2157  VectorizedArrayType::size()> &
2158  get_face_info(const unsigned int face_batch_index) const;
2159 
2160 
2166  const Table<3, unsigned int> &
2168 
2184 
2188  void
2190 
2202 
2206  void
2208  const AlignedVector<Number> *memory) const;
2209 
2212 private:
2217  template <typename number2, int q_dim>
2218  void
2220  const std::shared_ptr<hp::MappingCollection<dim>> &mapping,
2221  const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
2222  const std::vector<const AffineConstraints<number2> *> &constraint,
2223  const std::vector<IndexSet> &locally_owned_set,
2224  const std::vector<hp::QCollection<q_dim>> &quad,
2225  const AdditionalData &additional_data);
2226 
2233  template <typename number2>
2234  void
2236  const std::vector<const AffineConstraints<number2> *> &constraint,
2237  const std::vector<IndexSet> &locally_owned_set,
2238  const AdditionalData &additional_data);
2239 
2243  void
2245  const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
2246  const AdditionalData &additional_data);
2247 
2251  std::vector<SmartPointer<const DoFHandler<dim>>> dof_handlers;
2252 
2259  std::vector<SmartPointer<const AffineConstraints<Number>>> affine_constraints;
2260 
2265  std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
2266 
2273  std::vector<Number> constraint_pool_data;
2274 
2279  std::vector<unsigned int> constraint_pool_row_index;
2280 
2287 
2293 
2300  std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
2301 
2306  std::vector<unsigned int> mf_cell_indices;
2307 
2315 
2322 
2327  internal::MatrixFreeFunctions::FaceInfo<VectorizedArrayType::size()>
2329 
2334 
2339 
2348  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2350 
2355  mutable std::list<std::pair<bool, AlignedVector<Number>>>
2357 
2361  unsigned int mg_level;
2362 };
2363 
2364 
2365 
2366 /*----------------------- Inline functions ----------------------------------*/
2367 
2368 #ifndef DOXYGEN
2369 
2370 
2371 
2372 template <int dim, typename Number, typename VectorizedArrayType>
2373 template <typename T>
2374 inline void
2376  AlignedVector<T> &vec) const
2377 {
2378  vec.resize(this->n_cell_batches() + this->n_ghost_cell_batches());
2379 }
2380 
2381 
2382 
2383 template <int dim, typename Number, typename VectorizedArrayType>
2384 template <typename T>
2385 inline void
2387  AlignedVector<T> &vec) const
2388 {
2389  vec.resize(this->n_inner_face_batches() + this->n_boundary_face_batches() +
2390  this->n_ghost_inner_face_batches());
2391 }
2392 
2393 
2394 
2395 template <int dim, typename Number, typename VectorizedArrayType>
2396 template <typename VectorType>
2397 inline void
2399  VectorType &vec,
2400  const unsigned int comp) const
2401 {
2402  static_assert(IsBlockVector<VectorType>::value == false,
2403  "This function is not supported for block vectors.");
2404 
2405  Assert(task_info.n_procs == 1,
2406  ExcMessage("This function can only be used in serial."));
2407 
2408  AssertIndexRange(comp, n_components());
2409  vec.reinit(dof_info[comp].vector_partitioner->size());
2410 }
2411 
2412 
2413 
2414 template <int dim, typename Number, typename VectorizedArrayType>
2415 template <typename Number2>
2416 inline void
2419  const unsigned int comp) const
2420 {
2421  AssertIndexRange(comp, n_components());
2422  vec.reinit(dof_info[comp].vector_partitioner, task_info.communicator_sm);
2423 }
2424 
2425 
2426 
2427 template <int dim, typename Number, typename VectorizedArrayType>
2428 inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2430  const unsigned int comp) const
2431 {
2432  AssertIndexRange(comp, n_components());
2433  return dof_info[comp].vector_partitioner;
2434 }
2435 
2436 
2437 
2438 template <int dim, typename Number, typename VectorizedArrayType>
2439 inline const std::vector<unsigned int> &
2441  const unsigned int comp) const
2442 {
2443  AssertIndexRange(comp, n_components());
2444  return dof_info[comp].constrained_dofs;
2445 }
2446 
2447 
2448 
2449 template <int dim, typename Number, typename VectorizedArrayType>
2450 inline unsigned int
2452 {
2453  AssertDimension(dof_handlers.size(), dof_info.size());
2454  return dof_handlers.size();
2455 }
2456 
2457 
2458 
2459 template <int dim, typename Number, typename VectorizedArrayType>
2460 inline unsigned int
2462  const unsigned int dof_no) const
2463 {
2464  AssertDimension(dof_handlers.size(), dof_info.size());
2465  AssertIndexRange(dof_no, dof_handlers.size());
2466  return dof_handlers[dof_no]->get_fe().n_base_elements();
2467 }
2468 
2469 
2470 
2471 template <int dim, typename Number, typename VectorizedArrayType>
2474 {
2475  return task_info;
2476 }
2477 
2478 
2479 
2480 template <int dim, typename Number, typename VectorizedArrayType>
2481 inline unsigned int
2483 {
2484  return task_info.n_active_cells;
2485 }
2486 
2487 
2488 
2489 template <int dim, typename Number, typename VectorizedArrayType>
2490 inline unsigned int
2492 {
2493  return *(task_info.cell_partition_data.end() - 2);
2494 }
2495 
2496 
2497 
2498 template <int dim, typename Number, typename VectorizedArrayType>
2499 inline unsigned int
2501 {
2502  return *(task_info.cell_partition_data.end() - 1) -
2503  *(task_info.cell_partition_data.end() - 2);
2504 }
2505 
2506 
2507 
2508 template <int dim, typename Number, typename VectorizedArrayType>
2509 inline unsigned int
2511 {
2512  if (task_info.face_partition_data.empty())
2513  return 0;
2514  return task_info.face_partition_data.back();
2515 }
2516 
2517 
2518 
2519 template <int dim, typename Number, typename VectorizedArrayType>
2520 inline unsigned int
2522 {
2523  if (task_info.face_partition_data.empty())
2524  return 0;
2525  return task_info.boundary_partition_data.back() -
2527 }
2528 
2529 
2530 
2531 template <int dim, typename Number, typename VectorizedArrayType>
2532 inline unsigned int
2534 {
2535  if (task_info.face_partition_data.empty())
2536  return 0;
2537  return face_info.faces.size() - task_info.boundary_partition_data.back();
2538 }
2539 
2540 
2541 
2542 template <int dim, typename Number, typename VectorizedArrayType>
2543 inline types::boundary_id
2545  const unsigned int face_batch_index) const
2546 {
2547  Assert(face_batch_index >= task_info.boundary_partition_data[0] &&
2548  face_batch_index < task_info.boundary_partition_data.back(),
2549  ExcIndexRange(face_batch_index,
2552  return types::boundary_id(face_info.faces[face_batch_index].exterior_face_no);
2553 }
2554 
2555 
2556 
2557 template <int dim, typename Number, typename VectorizedArrayType>
2558 inline std::array<types::boundary_id, VectorizedArrayType::size()>
2560  const unsigned int cell_batch_index,
2561  const unsigned int face_number) const
2562 {
2563  AssertIndexRange(cell_batch_index, n_cell_batches());
2566  ExcNotInitialized());
2567  std::array<types::boundary_id, VectorizedArrayType::size()> result;
2568  result.fill(numbers::invalid_boundary_id);
2569  for (unsigned int v = 0;
2570  v < n_active_entries_per_cell_batch(cell_batch_index);
2571  ++v)
2572  result[v] =
2573  face_info.cell_and_face_boundary_id(cell_batch_index, face_number, v);
2574  return result;
2575 }
2576 
2577 
2578 
2579 template <int dim, typename Number, typename VectorizedArrayType>
2580 inline const internal::MatrixFreeFunctions::
2581  MappingInfo<dim, Number, VectorizedArrayType> &
2583 {
2584  return mapping_info;
2585 }
2586 
2587 
2588 
2589 template <int dim, typename Number, typename VectorizedArrayType>
2592  const unsigned int dof_index) const
2593 {
2594  AssertIndexRange(dof_index, n_components());
2595  return dof_info[dof_index];
2596 }
2597 
2598 
2599 
2600 template <int dim, typename Number, typename VectorizedArrayType>
2601 inline unsigned int
2603 {
2604  return constraint_pool_row_index.size() - 1;
2605 }
2606 
2607 
2608 
2609 template <int dim, typename Number, typename VectorizedArrayType>
2610 inline const Number *
2612  const unsigned int row) const
2613 {
2615  return constraint_pool_data.empty() ?
2616  nullptr :
2618 }
2619 
2620 
2621 
2622 template <int dim, typename Number, typename VectorizedArrayType>
2623 inline const Number *
2625  const unsigned int row) const
2626 {
2628  return constraint_pool_data.empty() ?
2629  nullptr :
2631 }
2632 
2633 
2634 
2635 template <int dim, typename Number, typename VectorizedArrayType>
2636 inline std::pair<unsigned int, unsigned int>
2638  const std::pair<unsigned int, unsigned int> &range,
2639  const unsigned int degree,
2640  const unsigned int dof_handler_component) const
2641 {
2642  if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2643  {
2645  dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2647  dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2648  if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2649  return range;
2650  else
2651  return {range.second, range.second};
2652  }
2653 
2654  const unsigned int fe_index =
2655  dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2656  if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2657  return {range.second, range.second};
2658  else
2659  return create_cell_subrange_hp_by_index(range,
2660  fe_index,
2661  dof_handler_component);
2662 }
2663 
2664 
2665 
2666 template <int dim, typename Number, typename VectorizedArrayType>
2667 inline bool
2669  const unsigned int cell_batch_index) const
2670 {
2671  AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2672  return VectorizedArrayType::size() > 1 &&
2673  cell_level_index[(cell_batch_index + 1) * VectorizedArrayType::size() -
2674  1] == cell_level_index[(cell_batch_index + 1) *
2675  VectorizedArrayType::size() -
2676  2];
2677 }
2678 
2679 
2680 
2681 template <int dim, typename Number, typename VectorizedArrayType>
2682 unsigned int
2684 {
2685  return shape_info.size(2);
2686 }
2687 
2688 
2689 template <int dim, typename Number, typename VectorizedArrayType>
2690 unsigned int
2692  const std::pair<unsigned int, unsigned int> range) const
2693 {
2694  const auto &fe_indices = dof_info[0].cell_active_fe_index;
2695 
2696  if (fe_indices.empty() == true)
2697  return 0;
2698 
2699  const auto index = fe_indices[range.first];
2700 
2701  for (unsigned int i = range.first; i < range.second; ++i)
2702  AssertDimension(index, fe_indices[i]);
2703 
2704  return index;
2705 }
2706 
2707 
2708 
2709 template <int dim, typename Number, typename VectorizedArrayType>
2710 unsigned int
2712  const std::pair<unsigned int, unsigned int> range,
2713  const bool is_interior_face) const
2714 {
2715  const auto &fe_indices = dof_info[0].cell_active_fe_index;
2716 
2717  if (fe_indices.empty() == true)
2718  return 0;
2719 
2720  if (is_interior_face)
2721  {
2722  const unsigned int index =
2723  fe_indices[face_info.faces[range.first].cells_interior[0] /
2724  VectorizedArrayType::size()];
2725 
2726  for (unsigned int i = range.first; i < range.second; ++i)
2727  AssertDimension(index,
2728  fe_indices[face_info.faces[i].cells_interior[0] /
2729  VectorizedArrayType::size()]);
2730 
2731  return index;
2732  }
2733  else
2734  {
2735  const unsigned int index =
2736  fe_indices[face_info.faces[range.first].cells_exterior[0] /
2737  VectorizedArrayType::size()];
2738 
2739  for (unsigned int i = range.first; i < range.second; ++i)
2740  AssertDimension(index,
2741  fe_indices[face_info.faces[i].cells_exterior[0] /
2742  VectorizedArrayType::size()]);
2743 
2744  return index;
2745  }
2746 }
2747 
2748 
2749 
2750 template <int dim, typename Number, typename VectorizedArrayType>
2751 inline unsigned int
2753  const unsigned int cell_batch_index) const
2754 {
2755  Assert(!dof_info.empty(), ExcNotInitialized());
2756  AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2757  const std::vector<unsigned char> &n_lanes_filled =
2758  dof_info[0].n_vectorization_lanes_filled
2760  AssertIndexRange(cell_batch_index, n_lanes_filled.size());
2761 
2762  return n_lanes_filled[cell_batch_index];
2763 }
2764 
2765 
2766 
2767 template <int dim, typename Number, typename VectorizedArrayType>
2768 inline unsigned int
2770  const unsigned int face_batch_index) const
2771 {
2772  AssertIndexRange(face_batch_index, face_info.faces.size());
2773  Assert(!dof_info.empty(), ExcNotInitialized());
2774  const std::vector<unsigned char> &n_lanes_filled =
2775  dof_info[0].n_vectorization_lanes_filled
2777  AssertIndexRange(face_batch_index, n_lanes_filled.size());
2778  return n_lanes_filled[face_batch_index];
2779 }
2780 
2781 
2782 
2783 template <int dim, typename Number, typename VectorizedArrayType>
2784 inline unsigned int
2786  const unsigned int dof_handler_index,
2787  const unsigned int active_fe_index) const
2788 {
2789  return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2790 }
2791 
2792 
2793 
2794 template <int dim, typename Number, typename VectorizedArrayType>
2795 inline unsigned int
2797  const unsigned int quad_index,
2798  const unsigned int active_fe_index) const
2799 {
2800  AssertIndexRange(quad_index, mapping_info.cell_data.size());
2801  return mapping_info.cell_data[quad_index]
2802  .descriptor[active_fe_index]
2803  .n_q_points;
2804 }
2805 
2806 
2807 
2808 template <int dim, typename Number, typename VectorizedArrayType>
2809 inline unsigned int
2811  const unsigned int dof_handler_index,
2812  const unsigned int active_fe_index) const
2813 {
2814  return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2815 }
2816 
2817 
2818 
2819 template <int dim, typename Number, typename VectorizedArrayType>
2820 inline unsigned int
2822  const unsigned int quad_index,
2823  const unsigned int active_fe_index) const
2824 {
2825  AssertIndexRange(quad_index, mapping_info.face_data.size());
2826  return mapping_info.face_data[quad_index]
2827  .descriptor[active_fe_index]
2828  .n_q_points;
2829 }
2830 
2831 
2832 
2833 template <int dim, typename Number, typename VectorizedArrayType>
2834 inline const IndexSet &
2836  const unsigned int dof_handler_index) const
2837 {
2838  return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2839 }
2840 
2841 
2842 
2843 template <int dim, typename Number, typename VectorizedArrayType>
2844 inline const IndexSet &
2846  const unsigned int dof_handler_index) const
2847 {
2848  return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2849 }
2850 
2851 
2852 
2853 template <int dim, typename Number, typename VectorizedArrayType>
2856  const unsigned int dof_handler_index,
2857  const unsigned int index_quad,
2858  const unsigned int index_fe,
2859  const unsigned int active_fe_index,
2860  const unsigned int active_quad_index) const
2861 {
2862  AssertIndexRange(dof_handler_index, dof_info.size());
2863  const unsigned int ind =
2864  dof_info[dof_handler_index].global_base_element_offset + index_fe;
2865  AssertIndexRange(ind, shape_info.size(0));
2866  AssertIndexRange(index_quad, shape_info.size(1));
2867  AssertIndexRange(active_fe_index, shape_info.size(2));
2868  AssertIndexRange(active_quad_index, shape_info.size(3));
2869  return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2870 }
2871 
2872 
2873 
2874 template <int dim, typename Number, typename VectorizedArrayType>
2876  VectorizedArrayType::size()> &
2878  const unsigned int face_batch_index) const
2879 {
2880  AssertIndexRange(face_batch_index, face_info.faces.size());
2881  return face_info.faces[face_batch_index];
2882 }
2883 
2884 
2885 
2886 template <int dim, typename Number, typename VectorizedArrayType>
2887 inline const Table<3, unsigned int> &
2889  const
2890 {
2892 }
2893 
2894 
2895 
2896 template <int dim, typename Number, typename VectorizedArrayType>
2897 inline const Quadrature<dim> &
2899  const unsigned int quad_index,
2900  const unsigned int active_fe_index) const
2901 {
2902  AssertIndexRange(quad_index, mapping_info.cell_data.size());
2903  return mapping_info.cell_data[quad_index]
2904  .descriptor[active_fe_index]
2905  .quadrature;
2906 }
2907 
2908 
2909 
2910 template <int dim, typename Number, typename VectorizedArrayType>
2911 inline const Quadrature<dim - 1> &
2913  const unsigned int quad_index,
2914  const unsigned int active_fe_index) const
2915 {
2916  AssertIndexRange(quad_index, mapping_info.face_data.size());
2917  return mapping_info.face_data[quad_index]
2918  .descriptor[active_fe_index]
2919  .quadrature;
2920 }
2921 
2922 
2923 
2924 template <int dim, typename Number, typename VectorizedArrayType>
2925 inline unsigned int
2927  const std::pair<unsigned int, unsigned int> range) const
2928 {
2929  auto result = get_cell_category(range.first);
2930 
2931  for (unsigned int i = range.first; i < range.second; ++i)
2932  result = std::max(result, get_cell_category(i));
2933 
2934  return result;
2935 }
2936 
2937 
2938 
2939 template <int dim, typename Number, typename VectorizedArrayType>
2940 inline std::pair<unsigned int, unsigned int>
2942  const std::pair<unsigned int, unsigned int> range) const
2943 {
2944  auto result = get_face_category(range.first);
2945 
2946  for (unsigned int i = range.first; i < range.second; ++i)
2947  {
2948  result.first = std::max(result.first, get_face_category(i).first);
2949  result.second = std::max(result.second, get_face_category(i).second);
2950  }
2951 
2952  return result;
2953 }
2954 
2955 
2956 
2957 template <int dim, typename Number, typename VectorizedArrayType>
2958 inline unsigned int
2960  const unsigned int cell_batch_index) const
2961 {
2962  AssertIndexRange(0, dof_info.size());
2963  AssertIndexRange(cell_batch_index, dof_info[0].cell_active_fe_index.size());
2964  if (dof_info[0].cell_active_fe_index.empty())
2965  return 0;
2966  else
2967  return dof_info[0].cell_active_fe_index[cell_batch_index];
2968 }
2969 
2970 
2971 
2972 template <int dim, typename Number, typename VectorizedArrayType>
2973 inline std::pair<unsigned int, unsigned int>
2975  const unsigned int face_batch_index) const
2976 {
2977  AssertIndexRange(face_batch_index, face_info.faces.size());
2978  if (dof_info[0].cell_active_fe_index.empty())
2979  return std::make_pair(0U, 0U);
2980 
2981  std::pair<unsigned int, unsigned int> result = std::make_pair(0U, 0U);
2982  for (unsigned int v = 0;
2983  v < VectorizedArrayType::size() &&
2984  face_info.faces[face_batch_index].cells_interior[v] !=
2986  ++v)
2987  result.first = std::max(
2988  result.first,
2989  dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
2990  .cells_interior[v] /
2991  VectorizedArrayType::size()]);
2992  if (face_info.faces[face_batch_index].cells_exterior[0] !=
2994  for (unsigned int v = 0;
2995  v < VectorizedArrayType::size() &&
2996  face_info.faces[face_batch_index].cells_exterior[v] !=
2998  ++v)
2999  result.second = std::max(
3000  result.second,
3001  dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
3002  .cells_exterior[v] /
3003  VectorizedArrayType::size()]);
3004  else
3005  result.second = numbers::invalid_unsigned_int;
3006  return result;
3007 }
3008 
3009 
3010 
3011 template <int dim, typename Number, typename VectorizedArrayType>
3012 inline bool
3014 {
3015  return indices_are_initialized;
3016 }
3017 
3018 
3019 
3020 template <int dim, typename Number, typename VectorizedArrayType>
3021 inline bool
3023 {
3024  return mapping_is_initialized;
3025 }
3026 
3027 
3028 template <int dim, typename Number, typename VectorizedArrayType>
3029 inline unsigned int
3031 {
3032  return mg_level;
3033 }
3034 
3035 
3036 
3037 template <int dim, typename Number, typename VectorizedArrayType>
3040 {
3041  using list_type =
3042  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3043  list_type &data = scratch_pad.get();
3044  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
3045  if (it->first == false)
3046  {
3047  it->first = true;
3048  return &it->second;
3049  }
3050  data.emplace_front(true, AlignedVector<VectorizedArrayType>());
3051  return &data.front().second;
3052 }
3053 
3054 
3055 
3056 template <int dim, typename Number, typename VectorizedArrayType>
3057 void
3059  const AlignedVector<VectorizedArrayType> *scratch) const
3060 {
3061  using list_type =
3062  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3063  list_type &data = scratch_pad.get();
3064  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
3065  if (&it->second == scratch)
3066  {
3067  Assert(it->first == true, ExcInternalError());
3068  it->first = false;
3069  return;
3070  }
3071  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
3072 }
3073 
3074 
3075 
3076 template <int dim, typename Number, typename VectorizedArrayType>
3080 {
3081  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
3083  it != scratch_pad_non_threadsafe.end();
3084  ++it)
3085  if (it->first == false)
3086  {
3087  it->first = true;
3088  return &it->second;
3089  }
3090  scratch_pad_non_threadsafe.push_front(
3091  std::make_pair(true, AlignedVector<Number>()));
3092  return &scratch_pad_non_threadsafe.front().second;
3093 }
3094 
3095 
3096 
3097 template <int dim, typename Number, typename VectorizedArrayType>
3098 void
3101  const AlignedVector<Number> *scratch) const
3102 {
3103  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
3105  it != scratch_pad_non_threadsafe.end();
3106  ++it)
3107  if (&it->second == scratch)
3108  {
3109  Assert(it->first == true, ExcInternalError());
3110  it->first = false;
3111  return;
3112  }
3113  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
3114 }
3115 
3116 
3117 
3118 // ------------------------------ reinit functions ---------------------------
3119 
3120 namespace internal
3121 {
3122  namespace MatrixFreeImplementation
3123  {
3124  template <int dim, int spacedim>
3125  inline std::vector<IndexSet>
3126  extract_locally_owned_index_sets(
3127  const std::vector<const DoFHandler<dim, spacedim> *> &dofh,
3128  const unsigned int level)
3129  {
3130  std::vector<IndexSet> locally_owned_set;
3131  locally_owned_set.reserve(dofh.size());
3132  for (unsigned int j = 0; j < dofh.size(); ++j)
3134  locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
3135  else
3136  locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
3137  return locally_owned_set;
3138  }
3139  } // namespace MatrixFreeImplementation
3140 } // namespace internal
3141 
3142 
3143 
3144 template <int dim, typename Number, typename VectorizedArrayType>
3145 template <typename QuadratureType, typename number2, typename MappingType>
3146 void
3148  const MappingType &mapping,
3149  const DoFHandler<dim> &dof_handler,
3150  const AffineConstraints<number2> &constraints_in,
3151  const QuadratureType &quad,
3153  &additional_data)
3154 {
3155  std::vector<const DoFHandler<dim, dim> *> dof_handlers;
3156  std::vector<const AffineConstraints<number2> *> constraints;
3157 
3158  dof_handlers.push_back(&dof_handler);
3159  constraints.push_back(&constraints_in);
3160 
3161  std::vector<IndexSet> locally_owned_sets =
3162  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3163  dof_handlers, additional_data.mg_level);
3164 
3165  std::vector<hp::QCollection<dim>> quad_hp;
3166  quad_hp.emplace_back(quad);
3167 
3168  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3169  dof_handlers,
3170  constraints,
3171  locally_owned_sets,
3172  quad_hp,
3173  additional_data);
3174 }
3175 
3176 
3177 
3178 template <int dim, typename Number, typename VectorizedArrayType>
3179 template <typename QuadratureType, typename number2, typename MappingType>
3180 void
3182  const MappingType &mapping,
3183  const std::vector<const DoFHandler<dim> *> &dof_handler,
3184  const std::vector<const AffineConstraints<number2> *> &constraint,
3185  const QuadratureType &quad,
3187  &additional_data)
3188 {
3189  std::vector<IndexSet> locally_owned_set =
3190  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3191  dof_handler, additional_data.mg_level);
3192  std::vector<hp::QCollection<dim>> quad_hp;
3193  quad_hp.emplace_back(quad);
3194 
3195  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3196  dof_handler,
3197  constraint,
3198  locally_owned_set,
3199  quad_hp,
3200  additional_data);
3201 }
3202 
3203 
3204 
3205 template <int dim, typename Number, typename VectorizedArrayType>
3206 template <typename QuadratureType, typename number2, typename MappingType>
3207 void
3209  const MappingType &mapping,
3210  const std::vector<const DoFHandler<dim> *> &dof_handler,
3211  const std::vector<const AffineConstraints<number2> *> &constraint,
3212  const std::vector<QuadratureType> &quad,
3214  &additional_data)
3215 {
3216  std::vector<IndexSet> locally_owned_set =
3217  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3218  dof_handler, additional_data.mg_level);
3219  std::vector<hp::QCollection<dim>> quad_hp;
3220  for (unsigned int q = 0; q < quad.size(); ++q)
3221  quad_hp.emplace_back(quad[q]);
3222 
3223  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3224  dof_handler,
3225  constraint,
3226  locally_owned_set,
3227  quad_hp,
3228  additional_data);
3229 }
3230 
3231 
3232 
3233 // ------------------------------ implementation of loops --------------------
3234 
3235 // internal helper functions that define how to call MPI data exchange
3236 // functions: for generic vectors, do nothing at all. For distributed vectors,
3237 // call update_ghost_values_start function and so on. If we have collections
3238 // of vectors, just do the individual functions of the components. In order to
3239 // keep ghost values consistent (whether we are in read or write mode), we
3240 // also reset the values at the end. the whole situation is a bit complicated
3241 // by the fact that we need to treat block vectors differently, which use some
3242 // additional helper functions to select the blocks and template magic.
3243 namespace internal
3244 {
3248  template <int dim, typename Number, typename VectorizedArrayType>
3249  struct VectorDataExchange
3250  {
3251  // A shift for the MPI messages to reduce the risk for accidental
3252  // interaction with other open communications that a user program might
3253  // set up (parallel vectors support unfinished communication). We let
3254  // the other vectors use the first 20 assigned numbers and start the
3255  // matrix-free communication.
3256  static constexpr unsigned int channel_shift = 20;
3257 
3258 
3259 
3264  VectorDataExchange(
3265  const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3266  const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3267  DataAccessOnFaces vector_face_access,
3268  const unsigned int n_components)
3269  : matrix_free(matrix_free)
3270  , vector_face_access(
3271  matrix_free.get_task_info().face_partition_data.empty() ?
3272  ::MatrixFree<dim, Number, VectorizedArrayType>::
3273  DataAccessOnFaces::unspecified :
3274  vector_face_access)
3275  , ghosts_were_set(false)
3276 # ifdef DEAL_II_WITH_MPI
3277  , tmp_data(n_components)
3278  , requests(n_components)
3279 # endif
3280  {
3281  (void)n_components;
3282  if (this->vector_face_access !=
3284  DataAccessOnFaces::unspecified)
3285  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3287  matrix_free.get_dof_info(c).vector_exchanger_face_variants.size(),
3288  5);
3289  }
3290 
3291 
3292 
3296  ~VectorDataExchange() // NOLINT
3297  {
3298 # ifdef DEAL_II_WITH_MPI
3299  for (unsigned int i = 0; i < tmp_data.size(); ++i)
3300  if (tmp_data[i] != nullptr)
3301  matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3302 # endif
3303  }
3304 
3305 
3306 
3311  template <typename VectorType>
3312  unsigned int
3313  find_vector_in_mf(const VectorType &vec,
3314  const bool check_global_compatibility = true) const
3315  {
3316  // case 1: vector was set up with MatrixFree::initialize_dof_vector()
3317  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3318  if (vec.get_partitioner().get() ==
3319  matrix_free.get_dof_info(c).vector_partitioner.get())
3320  return c;
3321 
3322  // case 2: user provided own partitioner (compatibility mode)
3323  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3324  if (check_global_compatibility ?
3325  vec.get_partitioner()->is_globally_compatible(
3326  *matrix_free.get_dof_info(c).vector_partitioner) :
3327  vec.get_partitioner()->is_compatible(
3328  *matrix_free.get_dof_info(c).vector_partitioner))
3329  return c;
3330 
3331  Assert(false,
3332  ExcNotImplemented("Could not find partitioner that fits vector"));
3333 
3335  }
3336 
3337 
3338 
3344  get_partitioner(const unsigned int mf_component) const
3345  {
3346  AssertDimension(matrix_free.get_dof_info(mf_component)
3347  .vector_exchanger_face_variants.size(),
3348  5);
3349  if (vector_face_access ==
3352  return *matrix_free.get_dof_info(mf_component)
3353  .vector_exchanger_face_variants[0];
3354  else if (vector_face_access ==
3357  return *matrix_free.get_dof_info(mf_component)
3358  .vector_exchanger_face_variants[1];
3359  else if (vector_face_access ==
3362  return *matrix_free.get_dof_info(mf_component)
3363  .vector_exchanger_face_variants[2];
3364  else if (vector_face_access ==
3366  DataAccessOnFaces::values_all_faces)
3367  return *matrix_free.get_dof_info(mf_component)
3368  .vector_exchanger_face_variants[3];
3369  else if (vector_face_access ==
3371  DataAccessOnFaces::gradients_all_faces)
3372  return *matrix_free.get_dof_info(mf_component)
3373  .vector_exchanger_face_variants[4];
3374  else
3375  return *matrix_free.get_dof_info(mf_component).vector_exchanger.get();
3376  }
3377 
3378 
3379 
3383  template <typename VectorType,
3384  std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3385  * = nullptr>
3386  void
3387  update_ghost_values_start(const unsigned int /*component_in_block_vector*/,
3388  const VectorType & /*vec*/)
3389  {}
3390 
3391 
3396  template <typename VectorType,
3397  std::enable_if_t<!has_update_ghost_values_start<VectorType> &&
3398  !is_not_parallel_vector<VectorType>,
3399  VectorType> * = nullptr>
3400  void
3401  update_ghost_values_start(const unsigned int component_in_block_vector,
3402  const VectorType &vec)
3403  {
3404  (void)component_in_block_vector;
3405  const bool ghosts_set = vec.has_ghost_elements();
3406 
3407  Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3408  ghosts_set == false,
3409  ExcNotImplemented());
3410 
3411  if (ghosts_set)
3412  ghosts_were_set = true;
3413 
3414  vec.update_ghost_values();
3415  }
3416 
3417 
3418 
3424  template <typename VectorType,
3425  std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3426  !has_exchange_on_subset<VectorType>,
3427  VectorType> * = nullptr>
3428  void
3429  update_ghost_values_start(const unsigned int component_in_block_vector,
3430  const VectorType &vec)
3431  {
3432  (void)component_in_block_vector;
3433  const bool ghosts_set = vec.has_ghost_elements();
3434 
3435  Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3436  ghosts_set == false,
3437  ExcNotImplemented());
3438 
3439  if (ghosts_set)
3440  ghosts_were_set = true;
3441 
3442  vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3443  }
3444 
3445 
3446 
3453  template <typename VectorType,
3454  std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3455  has_exchange_on_subset<VectorType>,
3456  VectorType> * = nullptr>
3457  void
3458  update_ghost_values_start(const unsigned int component_in_block_vector,
3459  const VectorType &vec)
3460  {
3461  static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3462  "Type mismatch between VectorType and VectorDataExchange");
3463  (void)component_in_block_vector;
3464  const bool ghosts_set = vec.has_ghost_elements();
3465 
3466  Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3467  ghosts_set == false,
3468  ExcNotImplemented());
3469 
3470  if (ghosts_set)
3471  ghosts_were_set = true;
3472 
3473  if (vec.size() != 0)
3474  {
3475 # ifdef DEAL_II_WITH_MPI
3476  const unsigned int mf_component = find_vector_in_mf(vec);
3477 
3478  const auto &part = get_partitioner(mf_component);
3479 
3480  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3481  part.n_import_sm_procs() == 0)
3482  return;
3483 
3484  tmp_data[component_in_block_vector] =
3485  matrix_free.acquire_scratch_data_non_threadsafe();
3486  tmp_data[component_in_block_vector]->resize_fast(
3487  part.n_import_indices());
3488  AssertDimension(requests.size(), tmp_data.size());
3489 
3490  part.export_to_ghosted_array_start(
3491  component_in_block_vector * 2 + channel_shift,
3492  ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3493  vec.shared_vector_data(),
3494  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3495  part.locally_owned_size(),
3496  matrix_free.get_dof_info(mf_component)
3497  .vector_partitioner->n_ghost_indices()),
3498  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3499  part.n_import_indices()),
3500  this->requests[component_in_block_vector]);
3501 # endif
3502  }
3503  }
3504 
3505 
3506 
3511  template <typename VectorType,
3512  std::enable_if_t<!has_update_ghost_values_start<VectorType>,
3513  VectorType> * = nullptr>
3514  void
3515  update_ghost_values_finish(const unsigned int /*component_in_block_vector*/,
3516  const VectorType & /*vec*/)
3517  {}
3518 
3519 
3520 
3526  template <typename VectorType,
3527  std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3528  !has_exchange_on_subset<VectorType>,
3529  VectorType> * = nullptr>
3530  void
3531  update_ghost_values_finish(const unsigned int component_in_block_vector,
3532  const VectorType &vec)
3533  {
3534  (void)component_in_block_vector;
3535  vec.update_ghost_values_finish();
3536  }
3537 
3538 
3539 
3546  template <typename VectorType,
3547  std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3548  has_exchange_on_subset<VectorType>,
3549  VectorType> * = nullptr>
3550  void
3551  update_ghost_values_finish(const unsigned int component_in_block_vector,
3552  const VectorType &vec)
3553  {
3554  static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3555  "Type mismatch between VectorType and VectorDataExchange");
3556  (void)component_in_block_vector;
3557 
3558  if (vec.size() != 0)
3559  {
3560 # ifdef DEAL_II_WITH_MPI
3561  AssertIndexRange(component_in_block_vector, tmp_data.size());
3562  AssertDimension(requests.size(), tmp_data.size());
3563 
3564  const unsigned int mf_component = find_vector_in_mf(vec);
3565 
3566  const auto &part = get_partitioner(mf_component);
3567 
3568  if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3569  part.n_import_sm_procs() != 0)
3570  {
3571  part.export_to_ghosted_array_finish(
3572  ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3573  vec.shared_vector_data(),
3574  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3575  part.locally_owned_size(),
3576  matrix_free.get_dof_info(mf_component)
3577  .vector_partitioner->n_ghost_indices()),
3578  this->requests[component_in_block_vector]);
3579 
3580  matrix_free.release_scratch_data_non_threadsafe(
3581  tmp_data[component_in_block_vector]);
3582  tmp_data[component_in_block_vector] = nullptr;
3583  }
3584 # endif
3585  }
3586  // let vector know that ghosts are being updated and we can read from
3587  // them
3588  vec.set_ghost_state(true);
3589  }
3590 
3591 
3592 
3596  template <typename VectorType,
3597  std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3598  * = nullptr>
3599  void
3600  compress_start(const unsigned int /*component_in_block_vector*/,
3601  VectorType & /*vec*/)
3602  {}
3603 
3604 
3605 
3610  template <typename VectorType,
3611  std::enable_if_t<!has_compress_start<VectorType> &&
3612  !is_not_parallel_vector<VectorType>,
3613  VectorType> * = nullptr>
3614  void
3615  compress_start(const unsigned int component_in_block_vector,
3616  VectorType &vec)
3617  {
3618  (void)component_in_block_vector;
3619  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3620  vec.compress(VectorOperation::add);
3621  }
3622 
3623 
3624 
3630  template <typename VectorType,
3631  std::enable_if_t<has_compress_start<VectorType> &&
3632  !has_exchange_on_subset<VectorType>,
3633  VectorType> * = nullptr>
3634  void
3635  compress_start(const unsigned int component_in_block_vector,
3636  VectorType &vec)
3637  {
3638  (void)component_in_block_vector;
3639  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3640  vec.compress_start(component_in_block_vector + channel_shift);
3641  }
3642 
3643 
3644 
3651  template <typename VectorType,
3652  std::enable_if_t<has_compress_start<VectorType> &&
3653  has_exchange_on_subset<VectorType>,
3654  VectorType> * = nullptr>
3655  void
3656  compress_start(const unsigned int component_in_block_vector,
3657  VectorType &vec)
3658  {
3659  static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3660  "Type mismatch between VectorType and VectorDataExchange");
3661  (void)component_in_block_vector;
3662  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3663 
3664  if (vec.size() != 0)
3665  {
3666 # ifdef DEAL_II_WITH_MPI
3667  const unsigned int mf_component = find_vector_in_mf(vec);
3668 
3669  const auto &part = get_partitioner(mf_component);
3670 
3671  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3672  part.n_import_sm_procs() == 0)
3673  return;
3674 
3675  tmp_data[component_in_block_vector] =
3676  matrix_free.acquire_scratch_data_non_threadsafe();
3677  tmp_data[component_in_block_vector]->resize_fast(
3678  part.n_import_indices());
3679  AssertDimension(requests.size(), tmp_data.size());
3680 
3681  part.import_from_ghosted_array_start(
3683  component_in_block_vector * 2 + channel_shift,
3684  ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3685  vec.shared_vector_data(),
3686  ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3687  matrix_free.get_dof_info(mf_component)
3688  .vector_partitioner->n_ghost_indices()),
3689  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3690  part.n_import_indices()),
3691  this->requests[component_in_block_vector]);
3692 # endif
3693  }
3694  }
3695 
3696 
3697 
3702  template <
3703  typename VectorType,
3704  std::enable_if_t<!has_compress_start<VectorType>, VectorType> * = nullptr>
3705  void
3706  compress_finish(const unsigned int /*component_in_block_vector*/,
3707  VectorType & /*vec*/)
3708  {}
3709 
3710 
3711 
3717  template <typename VectorType,
3718  std::enable_if_t<has_compress_start<VectorType> &&
3719  !has_exchange_on_subset<VectorType>,
3720  VectorType> * = nullptr>
3721  void
3722  compress_finish(const unsigned int component_in_block_vector,
3723  VectorType &vec)
3724  {
3725  (void)component_in_block_vector;
3726  vec.compress_finish(VectorOperation::add);
3727  }
3728 
3729 
3730 
3737  template <typename VectorType,
3738  std::enable_if_t<has_compress_start<VectorType> &&
3739  has_exchange_on_subset<VectorType>,
3740  VectorType> * = nullptr>
3741  void
3742  compress_finish(const unsigned int component_in_block_vector,
3743  VectorType &vec)
3744  {
3745  static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3746  "Type mismatch between VectorType and VectorDataExchange");
3747  (void)component_in_block_vector;
3748  if (vec.size() != 0)
3749  {
3750 # ifdef DEAL_II_WITH_MPI
3751  AssertIndexRange(component_in_block_vector, tmp_data.size());
3752  AssertDimension(requests.size(), tmp_data.size());
3753 
3754  const unsigned int mf_component = find_vector_in_mf(vec);
3755 
3756  const auto &part = get_partitioner(mf_component);
3757 
3758  if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3759  part.n_import_sm_procs() != 0)
3760  {
3761  part.import_from_ghosted_array_finish(
3763  ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3764  vec.shared_vector_data(),
3765  ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3766  matrix_free.get_dof_info(mf_component)
3767  .vector_partitioner->n_ghost_indices()),
3769  tmp_data[component_in_block_vector]->begin(),
3770  part.n_import_indices()),
3771  this->requests[component_in_block_vector]);
3772 
3773  matrix_free.release_scratch_data_non_threadsafe(
3774  tmp_data[component_in_block_vector]);
3775  tmp_data[component_in_block_vector] = nullptr;
3776  }
3777 
3779  {
3780  const int ierr =
3781  MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3782  AssertThrowMPI(ierr);
3783  }
3784 # endif
3785  }
3786  }
3787 
3788 
3789 
3793  template <typename VectorType,
3794  std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3795  * = nullptr>
3796  void
3797  reset_ghost_values(const VectorType & /*vec*/) const
3798  {}
3799 
3800 
3801 
3806  template <typename VectorType,
3807  std::enable_if_t<!has_exchange_on_subset<VectorType> &&
3808  !is_not_parallel_vector<VectorType>,
3809  VectorType> * = nullptr>
3810  void
3811  reset_ghost_values(const VectorType &vec) const
3812  {
3813  if (ghosts_were_set == true)
3814  return;
3815 
3816  vec.zero_out_ghost_values();
3817  }
3818 
3819 
3820 
3826  template <typename VectorType,
3827  std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3828  * = nullptr>
3829  void
3830  reset_ghost_values(const VectorType &vec) const
3831  {
3832  static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3833  "Type mismatch between VectorType and VectorDataExchange");
3834  if (ghosts_were_set == true)
3835  return;
3836 
3837  if (vec.size() != 0)
3838  {
3839 # ifdef DEAL_II_WITH_MPI
3840  AssertDimension(requests.size(), tmp_data.size());
3841 
3842  const unsigned int mf_component = find_vector_in_mf(vec);
3843 
3844  const auto &part = get_partitioner(mf_component);
3845 
3846  if (part.n_ghost_indices() > 0)
3847  {
3848  part.reset_ghost_values(ArrayView<Number>(
3850  .begin() +
3851  part.locally_owned_size(),
3852  matrix_free.get_dof_info(mf_component)
3853  .vector_partitioner->n_ghost_indices()));
3854  }
3855 
3856 # endif
3857  }
3858  // let vector know that it's not ghosted anymore
3859  vec.set_ghost_state(false);
3860  }
3861 
3862 
3863 
3869  template <typename VectorType,
3870  std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3871  * = nullptr>
3872  void
3873  zero_vector_region(const unsigned int range_index, VectorType &vec) const
3874  {
3875  static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3876  "Type mismatch between VectorType and VectorDataExchange");
3877  if (range_index == numbers::invalid_unsigned_int)
3878  vec = Number();
3879  else
3880  {
3881  const unsigned int mf_component = find_vector_in_mf(vec, false);
3882  const internal::MatrixFreeFunctions::DoFInfo &dof_info =
3883  matrix_free.get_dof_info(mf_component);
3884  Assert(dof_info.vector_zero_range_list_index.empty() == false,
3885  ExcNotInitialized());
3886 
3887  Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3888  ExcInternalError());
3889  AssertIndexRange(range_index,
3890  dof_info.vector_zero_range_list_index.size() - 1);
3891  for (unsigned int id =
3892  dof_info.vector_zero_range_list_index[range_index];
3893  id != dof_info.vector_zero_range_list_index[range_index + 1];
3894  ++id)
3895  std::memset(vec.begin() + dof_info.vector_zero_range_list[id].first,
3896  0,
3897  (dof_info.vector_zero_range_list[id].second -
3898  dof_info.vector_zero_range_list[id].first) *
3899  sizeof(Number));
3900  }
3901  }
3902 
3903 
3904 
3910  template <typename VectorType,
3911  std::enable_if_t<!has_exchange_on_subset<VectorType>, VectorType>
3912  * = nullptr,
3913  typename VectorType::value_type * = nullptr>
3914  void
3915  zero_vector_region(const unsigned int range_index, VectorType &vec) const
3916  {
3917  if (range_index == numbers::invalid_unsigned_int || range_index == 0)
3918  vec = typename VectorType::value_type();
3919  }
3920 
3921 
3922 
3927  void
3928  zero_vector_region(const unsigned int, ...) const
3929  {
3930  Assert(false,
3931  ExcNotImplemented("Zeroing is only implemented for vector types "
3932  "which provide VectorType::value_type"));
3933  }
3934 
3935 
3936 
3937  const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
3938  const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3939  DataAccessOnFaces vector_face_access;
3940  bool ghosts_were_set;
3941 # ifdef DEAL_II_WITH_MPI
3942  std::vector<AlignedVector<Number> *> tmp_data;
3943  std::vector<std::vector<MPI_Request>> requests;
3944 # endif
3945  }; // VectorDataExchange
3946 
3947  template <typename VectorStruct>
3948  unsigned int
3949  n_components(const VectorStruct &vec);
3950 
3951  template <typename VectorStruct>
3952  unsigned int
3953  n_components_block(const VectorStruct &vec,
3954  std::integral_constant<bool, true>)
3955  {
3956  unsigned int components = 0;
3957  for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
3958  components += n_components(vec.block(bl));
3959  return components;
3960  }
3961 
3962  template <typename VectorStruct>
3963  unsigned int
3964  n_components_block(const VectorStruct &, std::integral_constant<bool, false>)
3965  {
3966  return 1;
3967  }
3968 
3969  template <typename VectorStruct>
3970  unsigned int
3971  n_components(const VectorStruct &vec)
3972  {
3973  return n_components_block(
3974  vec, std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3975  }
3976 
3977  template <typename VectorStruct>
3978  inline unsigned int
3979  n_components(const std::vector<VectorStruct> &vec)
3980  {
3981  unsigned int components = 0;
3982  for (unsigned int comp = 0; comp < vec.size(); ++comp)
3983  components += n_components_block(
3984  vec[comp],
3985  std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3986  return components;
3987  }
3988 
3989  template <typename VectorStruct>
3990  inline unsigned int
3991  n_components(const std::vector<VectorStruct *> &vec)
3992  {
3993  unsigned int components = 0;
3994  for (unsigned int comp = 0; comp < vec.size(); ++comp)
3995  components += n_components_block(
3996  *vec[comp],
3997  std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3998  return components;
3999  }
4000 
4001 
4002 
4003  // A helper function to identify block vectors with many components where we
4004  // should not try to overlap computations and communication because there
4005  // would be too many outstanding communication requests.
4006 
4007  // default value for vectors that do not have communication_block_size
4008  template <typename VectorStruct,
4009  std::enable_if_t<!has_communication_block_size<VectorStruct>,
4010  VectorStruct> * = nullptr>
4011  constexpr unsigned int
4012  get_communication_block_size(const VectorStruct &)
4013  {
4015  }
4016 
4017 
4018 
4019  template <typename VectorStruct,
4020  std::enable_if_t<has_communication_block_size<VectorStruct>,
4021  VectorStruct> * = nullptr>
4022  constexpr unsigned int
4023  get_communication_block_size(const VectorStruct &)
4024  {
4025  return VectorStruct::communication_block_size;
4026  }
4027 
4028 
4029 
4030  template <typename VectorType,
4031  std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType> * =
4032  nullptr>
4033  bool
4034  has_ghost_elements(const VectorType &vec)
4035  {
4036  (void)vec;
4037  return false;
4038  }
4039 
4040 
4041 
4042  template <typename VectorType,
4043  std::enable_if_t<!is_not_parallel_vector<VectorType>, VectorType>
4044  * = nullptr>
4045  bool
4046  has_ghost_elements(const VectorType &vec)
4047  {
4048  return vec.has_ghost_elements();
4049  }
4050 
4051 
4052 
4053  // --------------------------------------------------------------------------
4054  // below we have wrappers to distinguish between block and non-block vectors.
4055  // --------------------------------------------------------------------------
4056 
4057  //
4058  // update_ghost_values_start
4059  //
4060 
4061  // update_ghost_values for block vectors
4062  template <int dim,
4063  typename VectorStruct,
4064  typename Number,
4065  typename VectorizedArrayType,
4066  std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4067  * = nullptr>
4068  void
4069  update_ghost_values_start(
4070  const VectorStruct &vec,
4071  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4072  const unsigned int channel = 0)
4073  {
4074  if (get_communication_block_size(vec) < vec.n_blocks())
4075  {
4076  const bool ghosts_set = vec.has_ghost_elements();
4077 
4078  Assert(exchanger.matrix_free.get_task_info()
4079  .allow_ghosted_vectors_in_loops ||
4080  ghosts_set == false,
4081  ExcNotImplemented());
4082 
4083  if (ghosts_set)
4084  exchanger.ghosts_were_set = true;
4085 
4086  vec.update_ghost_values();
4087  }
4088  else
4089  {
4090  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4091  update_ghost_values_start(vec.block(i), exchanger, channel + i);
4092  }
4093  }
4094 
4095 
4096 
4097  // update_ghost_values for non-block vectors
4098  template <int dim,
4099  typename VectorStruct,
4100  typename Number,
4101  typename VectorizedArrayType,
4102  std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4103  * = nullptr>
4104  void
4105  update_ghost_values_start(
4106  const VectorStruct &vec,
4107  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4108  const unsigned int channel = 0)
4109  {
4110  exchanger.update_ghost_values_start(channel, vec);
4111  }
4112 
4113 
4114 
4115  // update_ghost_values_start() for vector of vectors
4116  template <int dim,
4117  typename VectorStruct,
4118  typename Number,
4119  typename VectorizedArrayType>
4120  inline void
4121  update_ghost_values_start(
4122  const std::vector<VectorStruct> &vec,
4123  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4124  {
4125  unsigned int component_index = 0;
4126  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4127  {
4128  update_ghost_values_start(vec[comp], exchanger, component_index);
4129  component_index += n_components(vec[comp]);
4130  }
4131  }
4132 
4133 
4134 
4135  // update_ghost_values_start() for vector of pointers to vectors
4136  template <int dim,
4137  typename VectorStruct,
4138  typename Number,
4139  typename VectorizedArrayType>
4140  inline void
4141  update_ghost_values_start(
4142  const std::vector<VectorStruct *> &vec,
4143  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4144  {
4145  unsigned int component_index = 0;
4146  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4147  {
4148  update_ghost_values_start(*vec[comp], exchanger, component_index);
4149  component_index += n_components(*vec[comp]);
4150  }
4151  }
4152 
4153 
4154 
4155  //
4156  // update_ghost_values_finish
4157  //
4158 
4159  // for block vectors
4160  template <int dim,
4161  typename VectorStruct,
4162  typename Number,
4163  typename VectorizedArrayType,
4164  std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4165  * = nullptr>
4166  void
4167  update_ghost_values_finish(
4168  const VectorStruct &vec,
4169  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4170  const unsigned int channel = 0)
4171  {
4172  if (get_communication_block_size(vec) < vec.n_blocks())
4173  {
4174  // do nothing, everything has already been completed in the _start()
4175  // call
4176  }
4177  else
4178  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4179  update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4180  }
4181 
4182 
4183 
4184  // for non-block vectors
4185  template <int dim,
4186  typename VectorStruct,
4187  typename Number,
4188  typename VectorizedArrayType,
4189  std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4190  * = nullptr>
4191  void
4192  update_ghost_values_finish(
4193  const VectorStruct &vec,
4194  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4195  const unsigned int channel = 0)
4196  {
4197  exchanger.update_ghost_values_finish(channel, vec);
4198  }
4199 
4200 
4201 
4202  // for vector of vectors
4203  template <int dim,
4204  typename VectorStruct,
4205  typename Number,
4206  typename VectorizedArrayType>
4207  inline void
4208  update_ghost_values_finish(
4209  const std::vector<VectorStruct> &vec,
4210  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4211  {
4212  unsigned int component_index = 0;
4213  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4214  {
4215  update_ghost_values_finish(vec[comp], exchanger, component_index);
4216  component_index += n_components(vec[comp]);
4217  }
4218  }
4219 
4220 
4221 
4222  // for vector of pointers to vectors
4223  template <int dim,
4224  typename VectorStruct,
4225  typename Number,
4226  typename VectorizedArrayType>
4227  inline void
4228  update_ghost_values_finish(
4229  const std::vector<VectorStruct *> &vec,
4230  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4231  {
4232  unsigned int component_index = 0;
4233  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4234  {
4235  update_ghost_values_finish(*vec[comp], exchanger, component_index);
4236  component_index += n_components(*vec[comp]);
4237  }
4238  }
4239 
4240 
4241 
4242  //
4243  // compress_start
4244  //
4245 
4246  // for block vectors
4247  template <int dim,
4248  typename VectorStruct,
4249  typename Number,
4250  typename VectorizedArrayType,
4251  std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4252  * = nullptr>
4253  inline void
4254  compress_start(
4255  VectorStruct &vec,
4256  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4257  const unsigned int channel = 0)
4258  {
4259  if (get_communication_block_size(vec) < vec.n_blocks())
4260  vec.compress(VectorOperation::add);
4261  else
4262  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4263  compress_start(vec.block(i), exchanger, channel + i);
4264  }
4265 
4266 
4267 
4268  // for non-block vectors
4269  template <int dim,
4270  typename VectorStruct,
4271  typename Number,
4272  typename VectorizedArrayType,
4273  std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4274  * = nullptr>
4275  inline void
4276  compress_start(
4277  VectorStruct &vec,
4278  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4279  const unsigned int channel = 0)
4280  {
4281  exchanger.compress_start(channel, vec);
4282  }
4283 
4284 
4285 
4286  // for std::vector of vectors
4287  template <int dim,
4288  typename VectorStruct,
4289  typename Number,
4290  typename VectorizedArrayType>
4291  inline void
4292  compress_start(
4293  std::vector<VectorStruct> &vec,
4294  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4295  {
4296  unsigned int component_index = 0;
4297  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4298  {
4299  compress_start(vec[comp], exchanger, component_index);
4300  component_index += n_components(vec[comp]);
4301  }
4302  }
4303 
4304 
4305 
4306  // for std::vector of pointer to vectors
4307  template <int dim,
4308  typename VectorStruct,
4309  typename Number,
4310  typename VectorizedArrayType>
4311  inline void
4312  compress_start(
4313  std::vector<VectorStruct *> &vec,
4314  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4315  {
4316  unsigned int component_index = 0;
4317  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4318  {
4319  compress_start(*vec[comp], exchanger, component_index);
4320  component_index += n_components(*vec[comp]);
4321  }
4322  }
4323 
4324 
4325 
4326  //
4327  // compress_finish
4328  //
4329 
4330  // for block vectors
4331  template <int dim,
4332  typename VectorStruct,
4333  typename Number,
4334  typename VectorizedArrayType,
4335  std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4336  * = nullptr>
4337  inline void
4338  compress_finish(
4339  VectorStruct &vec,
4340  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4341  const unsigned int channel = 0)
4342  {
4343  if (get_communication_block_size(vec) < vec.n_blocks())
4344  {
4345  // do nothing, everything has already been completed in the _start()
4346  // call
4347  }
4348  else
4349  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4350  compress_finish(vec.block(i), exchanger, channel + i);
4351  }
4352 
4353 
4354 
4355  // for non-block vectors
4356  template <int dim,
4357  typename VectorStruct,
4358  typename Number,
4359  typename VectorizedArrayType,
4360  std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4361  * = nullptr>
4362  inline void
4363  compress_finish(
4364  VectorStruct &vec,
4365  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4366  const unsigned int channel = 0)
4367  {
4368  exchanger.compress_finish(channel, vec);
4369  }
4370 
4371 
4372 
4373  // for std::vector of vectors
4374  template <int dim,
4375  typename VectorStruct,
4376  typename Number,
4377  typename VectorizedArrayType>
4378  inline void
4379  compress_finish(
4380  std::vector<VectorStruct> &vec,
4381  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4382  {
4383  unsigned int component_index = 0;
4384  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4385  {
4386  compress_finish(vec[comp], exchanger, component_index);
4387  component_index += n_components(vec[comp]);
4388  }
4389  }
4390 
4391 
4392 
4393  // for std::vector of pointer to vectors
4394  template <int dim,
4395  typename VectorStruct,
4396  typename Number,
4397  typename VectorizedArrayType>
4398  inline void
4399  compress_finish(
4400  std::vector<VectorStruct *> &vec,
4401  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4402  {
4403  unsigned int component_index = 0;
4404  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4405  {
4406  compress_finish(*vec[comp], exchanger, component_index);
4407  component_index += n_components(*vec[comp]);
4408  }
4409  }
4410 
4411 
4412 
4413  //
4414  // reset_ghost_values:
4415  //
4416  // if the input vector did not have ghosts imported, clear them here again
4417  // in order to avoid subsequent operations e.g. in linear solvers to work
4418  // with ghosts all the time
4419  //
4420 
4421  // for block vectors
4422  template <int dim,
4423  typename VectorStruct,
4424  typename Number,
4425  typename VectorizedArrayType,
4426  std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4427  * = nullptr>
4428  inline void
4429  reset_ghost_values(
4430  const VectorStruct &vec,
4431  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4432  {
4433  // return immediately if there is nothing to do.
4434  if (exchanger.ghosts_were_set == true)
4435  return;
4436 
4437  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4438  reset_ghost_values(vec.block(i), exchanger);
4439  }
4440 
4441 
4442 
4443  // for non-block vectors
4444  template <int dim,
4445  typename VectorStruct,
4446  typename Number,
4447  typename VectorizedArrayType,
4448  std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4449  * = nullptr>
4450  inline void
4451  reset_ghost_values(
4452  const VectorStruct &vec,
4453  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4454  {
4455  exchanger.reset_ghost_values(vec);
4456  }
4457 
4458 
4459 
4460  // for std::vector of vectors
4461  template <int dim,
4462  typename VectorStruct,
4463  typename Number,
4464  typename VectorizedArrayType>
4465  inline void
4466  reset_ghost_values(
4467  const std::vector<VectorStruct> &vec,
4468  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4469  {
4470  // return immediately if there is nothing to do.
4471  if (exchanger.ghosts_were_set == true)
4472  return;
4473 
4474  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4475  reset_ghost_values(vec[comp], exchanger);
4476  }
4477 
4478 
4479 
4480  // for std::vector of pointer to vectors
4481  template <int dim,
4482  typename VectorStruct,
4483  typename Number,
4484  typename VectorizedArrayType>
4485  inline void
4486  reset_ghost_values(
4487  const std::vector<VectorStruct *> &vec,
4488  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4489  {
4490  // return immediately if there is nothing to do.
4491  if (exchanger.ghosts_were_set == true)
4492  return;
4493 
4494  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4495  reset_ghost_values(*vec[comp], exchanger);
4496  }
4497 
4498 
4499 
4500  //
4501  // zero_vector_region
4502  //
4503 
4504  // for block vectors
4505  template <int dim,
4506  typename VectorStruct,
4507  typename Number,
4508  typename VectorizedArrayType,
4509  std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4510  * = nullptr>
4511  inline void
4512  zero_vector_region(
4513  const unsigned int range_index,
4514  VectorStruct &vec,
4515  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4516  {
4517  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4518  exchanger.zero_vector_region(range_index, vec.block(i));
4519  }
4520 
4521 
4522 
4523  // for non-block vectors
4524  template <int dim,
4525  typename VectorStruct,
4526  typename Number,
4527  typename VectorizedArrayType,
4528  std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4529  * = nullptr>
4530  inline void
4531  zero_vector_region(
4532  const unsigned int range_index,
4533  VectorStruct &vec,
4534  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4535  {
4536  exchanger.zero_vector_region(range_index, vec);
4537  }
4538 
4539 
4540 
4541  // for std::vector of vectors
4542  template <int dim,
4543  typename VectorStruct,
4544  typename Number,
4545  typename VectorizedArrayType>
4546  inline void
4547  zero_vector_region(
4548  const unsigned int range_index,
4549  std::vector<VectorStruct> &vec,
4550  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4551  {
4552  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4553  zero_vector_region(range_index, vec[comp], exchanger);
4554  }
4555 
4556 
4557 
4558  // for std::vector of pointers to vectors
4559  template <int dim,
4560  typename VectorStruct,
4561  typename Number,
4562  typename VectorizedArrayType>
4563  inline void
4564  zero_vector_region(
4565  const unsigned int range_index,
4566  std::vector<VectorStruct *> &vec,
4567  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4568  {
4569  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4570  zero_vector_region(range_index, *vec[comp], exchanger);
4571  }
4572 
4573 
4574 
4575  // Apply a unit matrix operation to constrained DoFs: Default cases where we
4576  // cannot detect a LinearAlgebra::distributed::Vector, we do not do
4577  // anything, else we apply the constraints as a unit operation
4578  template <typename VectorStruct1, typename VectorStruct2>
4579  inline void
4580  apply_operation_to_constrained_dofs(const std::vector<unsigned int> &,
4581  const VectorStruct1 &,
4582  VectorStruct2 &)
4583  {}
4584 
4585  template <typename Number>
4586  inline void
4587  apply_operation_to_constrained_dofs(
4588  const std::vector<unsigned int> &constrained_dofs,
4591  {
4592  for (const unsigned int i : constrained_dofs)
4593  dst.local_element(i) = src.local_element(i);
4594  }
4595 
4596 
4597  namespace MatrixFreeFunctions
4598  {
4599  // struct to select between a const interface and a non-const interface
4600  // for MFWorker
4601  template <typename, typename, typename, typename, bool>
4602  struct InterfaceSelector
4603  {};
4604 
4605  // Version for constant functions
4606  template <typename MF,
4607  typename InVector,
4608  typename OutVector,
4609  typename Container>
4610  struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4611  {
4612  using function_type = void (Container::*)(
4613  const MF &,
4614  OutVector &,
4615  const InVector &,
4616  const std::pair<unsigned int, unsigned int> &) const;
4617  };
4618 
4619  // Version for non-constant functions
4620  template <typename MF,
4621  typename InVector,
4622  typename OutVector,
4623  typename Container>
4624  struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4625  {
4626  using function_type =
4627  void (Container::*)(const MF &,
4628  OutVector &,
4629  const InVector &,
4630  const std::pair<unsigned int, unsigned int> &);
4631  };
4632  } // namespace MatrixFreeFunctions
4633 
4634 
4635 
4636  // A implementation class for the worker object that runs the various
4637  // operations we want to perform during the matrix-free loop
4638  template <typename MF,
4639  typename InVector,
4640  typename OutVector,
4641  typename Container,
4642  bool is_constant>
4643  class MFWorker : public MFWorkerInterface
4644  {
4645  public:
4646  // An alias to make the arguments further down more readable
4647  using function_type = typename MatrixFreeFunctions::
4648  InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4649  function_type;
4650 
4651  // constructor, binds all the arguments to this class
4652  MFWorker(const MF &matrix_free,
4653  const InVector &src,
4654  OutVector &dst,
4655  const bool zero_dst_vector_setting,
4656  const Container &container,
4657  function_type cell_function,
4658  function_type face_function,
4659  function_type boundary_function,
4660  const typename MF::DataAccessOnFaces src_vector_face_access =
4662  const typename MF::DataAccessOnFaces dst_vector_face_access =
4664  const std::function<void(const unsigned int, const unsigned int)>
4665  &operation_before_loop = {},
4666  const std::function<void(const unsigned int, const unsigned int)>
4667  &operation_after_loop = {},
4668  const unsigned int dof_handler_index_pre_post = 0)
4669  : matrix_free(matrix_free)
4670  , container(const_cast<Container &>(container))
4671  , cell_function(cell_function)
4672  , face_function(face_function)
4673  , boundary_function(boundary_function)
4674  , src(src)
4675  , dst(dst)
4676  , src_data_exchanger(matrix_free,
4677  src_vector_face_access,
4678  n_components(src))
4679  , dst_data_exchanger(matrix_free,
4680  dst_vector_face_access,
4681  n_components(dst))
4682  , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4683  , zero_dst_vector_setting(zero_dst_vector_setting &&
4684  !src_and_dst_are_same)
4685  , operation_before_loop(operation_before_loop)
4686  , operation_after_loop(operation_after_loop)
4687  , dof_handler_index_pre_post(dof_handler_index_pre_post)
4688  {
4689  Assert(!has_ghost_elements(dst),
4690  ExcMessage("The destination vector passed to the matrix-free "
4691  "loop is ghosted. This is not allowed."));
4692  }
4693 
4694  // Runs the cell work. If no function is given, nothing is done
4695  virtual void
4696  cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4697  {
4698  if (cell_function != nullptr && cell_range.second > cell_range.first)
4699  for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4700  {
4701  const auto cell_subrange =
4702  matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4703 
4704  if (cell_subrange.second <= cell_subrange.first)
4705  continue;
4706 
4707  (container.*
4708  cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4709  }
4710  }
4711 
4712  virtual void
4713  cell(const unsigned int range_index) override
4714  {
4715  process_range(cell_function,
4716  matrix_free.get_task_info().cell_partition_data_hp_ptr,
4717  matrix_free.get_task_info().cell_partition_data_hp,
4718  range_index);
4719  }
4720 
4721  virtual void
4722  face(const unsigned int range_index) override
4723  {
4724  process_range(face_function,
4725  matrix_free.get_task_info().face_partition_data_hp_ptr,
4726  matrix_free.get_task_info().face_partition_data_hp,
4727  range_index);
4728  }
4729 
4730  virtual void
4731  boundary(const unsigned int range_index) override
4732  {
4733  process_range(boundary_function,
4734  matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4735  matrix_free.get_task_info().boundary_partition_data_hp,
4736  range_index);
4737  }
4738 
4739  private:
4740  void
4741  process_range(const function_type &fu,
4742  const std::vector<unsigned int> &ptr,
4743  const std::vector<unsigned int> &data,
4744  const unsigned int range_index)
4745  {
4746  if (fu == nullptr)
4747  return;
4748 
4749  AssertIndexRange(range_index + 1, ptr.size());
4750  for (unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4751  {
4752  AssertIndexRange(2 * i + 1, data.size());
4753  (container.*fu)(matrix_free,
4754  this->dst,
4755  this->src,
4756  std::make_pair(data[2 * i], data[2 * i + 1]));
4757  }
4758  }
4759 
4760  public:
4761  // Starts the communication for the update ghost values operation. We
4762  // cannot call this update if ghost and destination are the same because
4763  // that would introduce spurious entries in the destination (there is also
4764  // the problem that reading from a vector that we also write to is usually
4765  // not intended in case there is overlap, but this is up to the
4766  // application code to decide and we cannot catch this case here).
4767  virtual void
4768  vector_update_ghosts_start() override
4769  {
4770  if (!src_and_dst_are_same)
4771  internal::update_ghost_values_start(src, src_data_exchanger);
4772  }
4773 
4774  // Finishes the communication for the update ghost values operation
4775  virtual void
4776  vector_update_ghosts_finish() override
4777  {
4778  if (!src_and_dst_are_same)
4779  internal::update_ghost_values_finish(src, src_data_exchanger);
4780  }
4781 
4782  // Starts the communication for the vector compress operation
4783  virtual void
4784  vector_compress_start() override
4785  {
4786  internal::compress_start(dst, dst_data_exchanger);
4787  }
4788 
4789  // Finishes the communication for the vector compress operation
4790  virtual void
4791  vector_compress_finish() override
4792  {
4793  internal::compress_finish(dst, dst_data_exchanger);
4794  if (!src_and_dst_are_same)
4795  internal::reset_ghost_values(src, src_data_exchanger);
4796  }
4797 
4798  // Zeros the given input vector
4799  virtual void
4800  zero_dst_vector_range(const unsigned int range_index) override
4801  {
4802  if (zero_dst_vector_setting)
4803  internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4804  }
4805 
4806  virtual void
4807  cell_loop_pre_range(const unsigned int range_index) override
4808  {
4809  if (operation_before_loop)
4810  {
4811  const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4812  matrix_free.get_dof_info(dof_handler_index_pre_post);
4813  if (range_index == numbers::invalid_unsigned_int)
4814  {
4815  // Case with threaded loop -> currently no overlap implemented
4817  0U,
4818  dof_info.vector_partitioner->locally_owned_size(),
4819  operation_before_loop,
4821  }
4822  else
4823  {
4824  AssertIndexRange(range_index,
4825  dof_info.cell_loop_pre_list_index.size() - 1);
4826  for (unsigned int id =
4827  dof_info.cell_loop_pre_list_index[range_index];
4828  id != dof_info.cell_loop_pre_list_index[range_index + 1];
4829  ++id)
4830  operation_before_loop(dof_info.cell_loop_pre_list[id].first,
4831  dof_info.cell_loop_pre_list[id].second);
4832  }
4833  }
4834  }
4835 
4836  virtual void
4837  cell_loop_post_range(const unsigned int range_index) override
4838  {
4839  if (operation_after_loop)
4840  {
4841  // Run unit matrix operation on constrained dofs if we are at the
4842  // last range
4843  const std::vector<unsigned int> &partition_row_index =
4844  matrix_free.get_task_info().partition_row_index;
4845  if (range_index ==
4846  partition_row_index[partition_row_index.size() - 2] - 1)
4847  apply_operation_to_constrained_dofs(
4848  matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
4849  src,
4850  dst);
4851 
4852  const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4853  matrix_free.get_dof_info(dof_handler_index_pre_post);
4854  if (range_index == numbers::invalid_unsigned_int)
4855  {
4856  // Case with threaded loop -> currently no overlap implemented
4858  0U,
4859  dof_info.vector_partitioner->locally_owned_size(),
4860  operation_after_loop,
4862  }
4863  else
4864  {
4865  AssertIndexRange(range_index,
4866  dof_info.cell_loop_post_list_index.size() - 1);
4867  for (unsigned int id =
4868  dof_info.cell_loop_post_list_index[range_index];
4869  id != dof_info.cell_loop_post_list_index[range_index + 1];
4870  ++id)
4871  operation_after_loop(dof_info.cell_loop_post_list[id].first,
4872  dof_info.cell_loop_post_list[id].second);
4873  }
4874  }
4875  }
4876 
4877  private:
4878  const MF &matrix_free;
4879  Container &container;
4880  function_type cell_function;
4881  function_type face_function;
4882  function_type boundary_function;
4883 
4884  const InVector &src;
4885  OutVector &dst;
4886  VectorDataExchange<MF::dimension,
4887  typename MF::value_type,
4888  typename MF::vectorized_value_type>
4889  src_data_exchanger;
4890  VectorDataExchange<MF::dimension,
4891  typename MF::value_type,
4892  typename MF::vectorized_value_type>
4893  dst_data_exchanger;
4894  const bool src_and_dst_are_same;
4895  const bool zero_dst_vector_setting;
4896  const std::function<void(const unsigned int, const unsigned int)>
4897  operation_before_loop;
4898  const std::function<void(const unsigned int, const unsigned int)>
4899  operation_after_loop;
4900  const unsigned int dof_handler_index_pre_post;
4901  };
4902 
4903 
4904 
4909  template <class MF, typename InVector, typename OutVector>
4910  struct MFClassWrapper
4911  {
4912  using function_type =
4913  std::function<void(const MF &,
4914  OutVector &,
4915  const InVector &,
4916  const std::pair<unsigned int, unsigned int> &)>;
4917 
4918  MFClassWrapper(const function_type cell,
4919  const function_type face,
4920  const function_type boundary)
4921  : cell(cell)
4922  , face(face)
4923  , boundary(boundary)
4924  {}
4925 
4926  void
4927  cell_integrator(const MF &mf,
4928  OutVector &dst,
4929  const InVector &src,
4930  const std::pair<unsigned int, unsigned int> &range) const
4931  {
4932  if (cell)
4933  cell(mf, dst, src, range);
4934  }
4935 
4936  void
4937  face_integrator(const MF &mf,
4938  OutVector &dst,
4939  const InVector &src,
4940  const std::pair<unsigned int, unsigned int> &range) const
4941  {
4942  if (face)
4943  face(mf, dst, src, range);
4944  }
4945 
4946  void
4947  boundary_integrator(
4948  const MF &mf,
4949  OutVector &dst,
4950  const InVector &src,
4951  const std::pair<unsigned int, unsigned int> &range) const
4952  {
4953  if (boundary)
4954  boundary(mf, dst, src, range);
4955  }
4956 
4957  const function_type cell;
4958  const function_type face;
4959  const function_type boundary;
4960  };
4961 
4962 } // end of namespace internal
4963 
4964 
4965 
4966 template <int dim, typename Number, typename VectorizedArrayType>
4967 template <typename OutVector, typename InVector>
4968 inline void
4970  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4971  OutVector &,
4972  const InVector &,
4973  const std::pair<unsigned int, unsigned int> &)>
4974  &cell_operation,
4975  OutVector &dst,
4976  const InVector &src,
4977  const bool zero_dst_vector) const
4978 {
4979  using Wrapper =
4980  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4981  InVector,
4982  OutVector>;
4983  Wrapper wrap(cell_operation, nullptr, nullptr);
4984  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4985  InVector,
4986  OutVector,
4987  Wrapper,
4988  true>
4989  worker(*this,
4990  src,
4991  dst,
4992  zero_dst_vector,
4993  wrap,
4994  &Wrapper::cell_integrator,
4995  &Wrapper::face_integrator,
4996  &Wrapper::boundary_integrator);
4997 
4998  task_info.loop(worker);
4999 }
5000 
5001 
5002 
5003 template <int dim, typename Number, typename VectorizedArrayType>
5004 template <typename OutVector, typename InVector>
5005 inline void
5007  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5008  OutVector &,
5009  const InVector &,
5010  const std::pair<unsigned int, unsigned int> &)>
5011  &cell_operation,
5012  OutVector &dst,
5013  const InVector &src,
5014  const std::function<void(const unsigned int, const unsigned int)>
5015  &operation_before_loop,
5016  const std::function<void(const unsigned int, const unsigned int)>
5017  &operation_after_loop,
5018  const unsigned int dof_handler_index_pre_post) const
5019 {
5020  using Wrapper =
5021  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5022  InVector,
5023  OutVector>;
5024  Wrapper wrap(cell_operation, nullptr, nullptr);
5025  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5026  InVector,
5027  OutVector,
5028  Wrapper,
5029  true>
5030  worker(*this,
5031  src,
5032  dst,
5033  false,
5034  wrap,
5035  &Wrapper::cell_integrator,
5036  &Wrapper::face_integrator,
5037  &Wrapper::boundary_integrator,
5040  operation_before_loop,
5041  operation_after_loop,
5042  dof_handler_index_pre_post);
5043 
5044  task_info.loop(worker);
5045 }
5046 
5047 
5048 
5049 template <int dim, typename Number, typename VectorizedArrayType>
5050 template <typename OutVector, typename InVector>
5051 inline void
5053  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5054  OutVector &,
5055  const InVector &,
5056  const std::pair<unsigned int, unsigned int> &)>
5057  &cell_operation,
5058  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5059  OutVector &,
5060  const InVector &,
5061  const std::pair<unsigned int, unsigned int> &)>
5062  &face_operation,
5063  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5064  OutVector &,
5065  const InVector &,
5066  const std::pair<unsigned int, unsigned int> &)>
5067  &boundary_operation,
5068  OutVector &dst,
5069  const InVector &src,
5070  const bool zero_dst_vector,
5071  const DataAccessOnFaces dst_vector_face_access,
5072  const DataAccessOnFaces src_vector_face_access) const
5073 {
5074  using Wrapper =
5075  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5076  InVector,
5077  OutVector>;
5078  Wrapper wrap(cell_operation, face_operation, boundary_operation);
5079  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5080  InVector,
5081  OutVector,
5082  Wrapper,
5083  true>
5084  worker(*this,
5085  src,
5086  dst,
5087  zero_dst_vector,
5088  wrap,
5089  &Wrapper::cell_integrator,
5090  &Wrapper::face_integrator,
5091  &Wrapper::boundary_integrator,
5092  src_vector_face_access,
5093  dst_vector_face_access);
5094 
5095  task_info.loop(worker);
5096 }
5097 
5098 
5099 
5100 template <int dim, typename Number, typename VectorizedArrayType>
5101 template <typename CLASS, typename OutVector, typename InVector>
5102 inline void
5104  void (CLASS::*function_pointer)(
5106  OutVector &,
5107  const InVector &,
5108  const std::pair<unsigned int, unsigned int> &) const,
5109  const CLASS *owning_class,
5110  OutVector &dst,
5111  const InVector &src,
5112  const bool zero_dst_vector) const
5113 {
5114  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5115  InVector,
5116  OutVector,
5117  CLASS,
5118  true>
5119  worker(*this,
5120  src,
5121  dst,
5122  zero_dst_vector,
5123  *owning_class,
5124  function_pointer,
5125  nullptr,
5126  nullptr);
5127  task_info.loop(worker);
5128 }
5129 
5130 
5131 
5132 template <int dim, typename Number, typename VectorizedArrayType>
5133 template <typename CLASS, typename OutVector, typename InVector>
5134 inline void
5136  void (CLASS::*function_pointer)(
5138  OutVector &,
5139  const InVector &,
5140  const std::pair<unsigned int, unsigned int> &) const,
5141  const CLASS *owning_class,
5142  OutVector &dst,
5143  const InVector &src,
5144  const std::function<void(const unsigned int, const unsigned int)>
5145  &operation_before_loop,
5146  const std::function<void(const unsigned int, const unsigned int)>
5147  &operation_after_loop,
5148  const unsigned int dof_handler_index_pre_post) const
5149 {
5150  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5151  InVector,
5152  OutVector,
5153  CLASS,
5154  true>
5155  worker(*this,
5156  src,
5157  dst,
5158  false,
5159  *owning_class,
5160  function_pointer,
5161  nullptr,
5162  nullptr,
5165  operation_before_loop,
5166  operation_after_loop,
5167  dof_handler_index_pre_post);
5168  task_info.loop(worker);
5169 }
5170 
5171 
5172 
5173 template <int dim, typename Number, typename VectorizedArrayType>
5174 template <typename CLASS, typename OutVector, typename InVector>
5175 inline void
5177  void (CLASS::*cell_operation)(
5179  OutVector &,
5180  const InVector &,
5181  const std::pair<unsigned int, unsigned int> &) const,
5182  void (CLASS::*face_operation)(
5184  OutVector &,
5185  const InVector &,
5186  const std::pair<unsigned int, unsigned int> &) const,
5187  void (CLASS::*boundary_operation)(
5189  OutVector &,
5190  const InVector &,
5191  const std::pair<unsigned int, unsigned int> &) const,
5192  const CLASS *owning_class,
5193  OutVector &dst,
5194  const InVector &src,
5195  const bool zero_dst_vector,
5196  const DataAccessOnFaces dst_vector_face_access,
5197  const DataAccessOnFaces src_vector_face_access) const
5198 {
5199  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5200  InVector,
5201  OutVector,
5202  CLASS,
5203  true>
5204  worker(*this,
5205  src,
5206  dst,
5207  zero_dst_vector,
5208  *owning_class,
5209  cell_operation,
5210  face_operation,
5211  boundary_operation,
5212  src_vector_face_access,
5213  dst_vector_face_access);
5214  task_info.loop(worker);
5215 }
5216 
5217 
5218 
5219 template <int dim, typename Number, typename VectorizedArrayType>
5220 template <typename CLASS, typename OutVector, typename InVector>
5221 inline void
5223  void (CLASS::*function_pointer)(
5225  OutVector &,
5226  const InVector &,
5227  const std::pair<unsigned int, unsigned int> &),
5228  CLASS *owning_class,
5229  OutVector &dst,
5230  const InVector &src,
5231  const bool zero_dst_vector) const
5232 {
5233  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5234  InVector,
5235  OutVector,
5236  CLASS,
5237  false>
5238  worker(*this,
5239  src,
5240  dst,
5241  zero_dst_vector,
5242  *owning_class,
5243  function_pointer,
5244  nullptr,
5245  nullptr);
5246  task_info.loop(worker);
5247 }
5248 
5249 
5250 
5251 template <int dim, typename Number, typename VectorizedArrayType>
5252 template <typename CLASS, typename OutVector, typename InVector>
5253 inline void
5255  void (CLASS::*function_pointer)(
5257  OutVector &,
5258  const InVector &,
5259  const std::pair<unsigned int, unsigned int> &),
5260  CLASS *owning_class,
5261  OutVector &dst,
5262  const InVector &src,
5263  const std::function<void(const unsigned int, const unsigned int)>
5264  &operation_before_loop,
5265  const std::function<void(const unsigned int, const unsigned int)>
5266  &operation_after_loop,
5267  const unsigned int dof_handler_index_pre_post) const
5268 {
5269  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5270  InVector,
5271  OutVector,
5272  CLASS,
5273  false>
5274  worker(*this,
5275  src,
5276  dst,
5277  false,
5278  *owning_class,
5279  function_pointer,
5280  nullptr,
5281  nullptr,
5284  operation_before_loop,
5285  operation_after_loop,
5286  dof_handler_index_pre_post);
5287  task_info.loop(worker);
5288 }
5289 
5290 
5291 
5292 template <int dim, typename Number, typename VectorizedArrayType>
5293 template <typename CLASS, typename OutVector, typename InVector>
5294 inline void
5296  void (CLASS::*cell_operation)(
5298  OutVector &,
5299  const InVector &,
5300  const std::pair<unsigned int, unsigned int> &),
5301  void (CLASS::*face_operation)(
5303  OutVector &,
5304  const InVector &,
5305  const std::pair<unsigned int, unsigned int> &),
5306  void (CLASS::*boundary_operation)(
5308  OutVector &,
5309  const InVector &,
5310  const std::pair<unsigned int, unsigned int> &),
5311  CLASS *owning_class,
5312  OutVector &dst,
5313  const InVector &src,
5314  const bool zero_dst_vector,
5315  const DataAccessOnFaces dst_vector_face_access,
5316  const DataAccessOnFaces src_vector_face_access) const
5317 {
5318  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5319  InVector,
5320  OutVector,
5321  CLASS,
5322  false>
5323  worker(*this,
5324  src,
5325  dst,
5326  zero_dst_vector,
5327  *owning_class,
5328  cell_operation,
5329  face_operation,
5330  boundary_operation,
5331  src_vector_face_access,
5332  dst_vector_face_access);
5333  task_info.loop(worker);
5334 }
5335 
5336 
5337 
5338 template <int dim, typename Number, typename VectorizedArrayType>
5339 template <typename OutVector, typename InVector>
5340 inline void
5342  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5343  OutVector &,
5344  const InVector &,
5345  const std::pair<unsigned int, unsigned int> &)>
5346  &cell_operation,
5347  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5348  OutVector &,
5349  const InVector &,
5350  const std::pair<unsigned int, unsigned int> &)>
5351  &face_operation,
5352  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5353  OutVector &,
5354  const InVector &,
5355  const std::pair<unsigned int, unsigned int> &)>
5356  &boundary_operation,
5357  OutVector &dst,
5358  const InVector &src,
5359  const std::function<void(const unsigned int, const unsigned int)>
5360  &operation_before_loop,
5361  const std::function<void(const unsigned int, const unsigned int)>
5362  &operation_after_loop,
5363  const unsigned int dof_handler_index_pre_post,
5364  const DataAccessOnFaces dst_vector_face_access,
5365  const DataAccessOnFaces src_vector_face_access) const
5366 {
5367  using Wrapper =
5368  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5369  InVector,
5370  OutVector>;
5371  Wrapper wrap(cell_operation, face_operation, boundary_operation);
5372  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5373  InVector,
5374  OutVector,
5375  Wrapper,
5376  true>
5377  worker(*this,
5378  src,
5379  dst,
5380  false,
5381  wrap,
5382  &Wrapper::cell_integrator,
5383  &Wrapper::face_integrator,
5384  &Wrapper::boundary_integrator,
5385  src_vector_face_access,
5386  dst_vector_face_access,
5387  operation_before_loop,
5388  operation_after_loop,
5389  dof_handler_index_pre_post);
5390 
5391  task_info.loop(worker);
5392 }
5393 
5394 
5395 
5396 template <int dim, typename Number, typename VectorizedArrayType>
5397 template <typename CLASS, typename OutVector, typename InVector>
5398 inline void
5400  void (CLASS::*cell_operation)(const MatrixFree &,
5401  OutVector &,
5402  const InVector &,
5403  const std::pair<unsigned int, unsigned int> &)
5404  const,
5405  void (CLASS::*face_operation)(const MatrixFree &,
5406  OutVector &,
5407  const InVector &,
5408  const std::pair<unsigned int, unsigned int> &)
5409  const,
5410  void (CLASS::*boundary_operation)(
5411  const MatrixFree &,
5412  OutVector &,
5413  const InVector &,
5414  const std::pair<unsigned int, unsigned int> &) const,
5415  const CLASS *owning_class,
5416  OutVector &dst,
5417  const InVector &src,
5418  const std::function<void(const unsigned int, const unsigned int)>
5419  &operation_before_loop,
5420  const std::function<void(const unsigned int, const unsigned int)>
5421  &operation_after_loop,
5422  const unsigned int dof_handler_index_pre_post,
5423  const DataAccessOnFaces dst_vector_face_access,
5424  const DataAccessOnFaces src_vector_face_access) const
5425 {
5426  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5427  InVector,
5428  OutVector,
5429  CLASS,
5430  true>
5431  worker(*this,
5432  src,
5433  dst,
5434  false,
5435  *owning_class,
5436  cell_operation,
5437  face_operation,
5438  boundary_operation,
5439  src_vector_face_access,
5440  dst_vector_face_access,
5441  operation_before_loop,
5442  operation_after_loop,
5443  dof_handler_index_pre_post);
5444  task_info.loop(worker);
5445 }
5446 
5447 
5448 
5449 template <int dim, typename Number, typename VectorizedArrayType>
5450 template <typename CLASS, typename OutVector, typename InVector>
5451 inline void
5453  void (CLASS::*cell_operation)(const MatrixFree &,
5454  OutVector &,
5455  const InVector &,
5456  const std::pair<unsigned int, unsigned int> &),
5457  void (CLASS::*face_operation)(const MatrixFree &,
5458  OutVector &,
5459  const InVector &,
5460  const std::pair<unsigned int, unsigned int> &),
5461  void (CLASS::*boundary_operation)(
5462  const MatrixFree &,
5463  OutVector &,
5464  const InVector &,
5465  const std::pair<unsigned int, unsigned int> &),
5466  const CLASS *owning_class,
5467  OutVector &dst,
5468  const InVector &src,
5469  const std::function<void(const unsigned int, const unsigned int)>
5470  &operation_before_loop,
5471  const std::function<void(const unsigned int, const unsigned int)>
5472  &operation_after_loop,
5473  const unsigned int dof_handler_index_pre_post,
5474  const DataAccessOnFaces dst_vector_face_access,
5475  const DataAccessOnFaces src_vector_face_access) const
5476 {
5477  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5478  InVector,
5479  OutVector,
5480  CLASS,
5481  false>
5482  worker(*this,
5483  src,
5484  dst,
5485  false,
5486  *owning_class,
5487  cell_operation,
5488  face_operation,
5489  boundary_operation,
5490  src_vector_face_access,
5491  dst_vector_face_access,
5492  operation_before_loop,
5493  operation_after_loop,
5494  dof_handler_index_pre_post);
5495  task_info.loop(worker);
5496 }
5497 
5498 
5499 
5500 template <int dim, typename Number, typename VectorizedArrayType>
5501 template <typename CLASS, typename OutVector, typename InVector>
5502 inline void
5504  void (CLASS::*function_pointer)(
5506  OutVector &,
5507  const InVector &,
5508  const std::pair<unsigned int, unsigned int> &) const,
5509  const CLASS *owning_class,
5510  OutVector &dst,
5511  const InVector &src,
5512  const bool zero_dst_vector,
5513  const DataAccessOnFaces src_vector_face_access) const
5514 {
5515  auto src_vector_face_access_temp = src_vector_face_access;
5516  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5517  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5518  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5519  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5520 
5521  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5522  InVector,
5523  OutVector,
5524  CLASS,
5525  true>
5526  worker(*this,
5527  src,
5528  dst,
5529  zero_dst_vector,
5530  *owning_class,
5531  function_pointer,
5532  nullptr,
5533  nullptr,
5534  src_vector_face_access_temp,
5536  task_info.loop(worker);
5537 }
5538 
5539 
5540 
5541 template <int dim, typename Number, typename VectorizedArrayType>
5542 template <typename CLASS, typename OutVector, typename InVector>
5543 inline void
5545  void (CLASS::*function_pointer)(
5547  OutVector &,
5548  const InVector &,
5549  const std::pair<unsigned int, unsigned int> &),
5550  CLASS *owning_class,
5551  OutVector &dst,
5552  const InVector &src,
5553  const bool zero_dst_vector,
5554  const DataAccessOnFaces src_vector_face_access) const
5555 {
5556  auto src_vector_face_access_temp = src_vector_face_access;
5557  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5558  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5559  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5560  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5561 
5562  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5563  InVector,
5564  OutVector,
5565  CLASS,
5566  false>
5567  worker(*this,
5568  src,
5569  dst,
5570  zero_dst_vector,
5571  *owning_class,
5572  function_pointer,
5573  nullptr,
5574  nullptr,
5575  src_vector_face_access_temp,
5577  task_info.loop(worker);
5578 }
5579 
5580 
5581 
5582 template <int dim, typename Number, typename VectorizedArrayType>
5583 template <typename OutVector, typename InVector>
5584 inline void
5586  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5587  OutVector &,
5588  const InVector &,
5589  const std::pair<unsigned int, unsigned int> &)>
5590  &cell_operation,
5591  OutVector &dst,
5592  const InVector &src,
5593  const bool zero_dst_vector,
5594  const DataAccessOnFaces src_vector_face_access) const
5595 {
5596  auto src_vector_face_access_temp = src_vector_face_access;
5597  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5598  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5599  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5600  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5601 
5602  using Wrapper =
5603  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5604  InVector,
5605  OutVector>;
5606  Wrapper wrap(cell_operation, nullptr, nullptr);
5607 
5608  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5609  InVector,
5610  OutVector,
5611  Wrapper,
5612  true>
5613  worker(*this,
5614  src,
5615  dst,
5616  zero_dst_vector,
5617  wrap,
5618  &Wrapper::cell_integrator,
5619  &Wrapper::face_integrator,
5620  &Wrapper::boundary_integrator,
5621  src_vector_face_access_temp,
5623  task_info.loop(worker);
5624 }
5625 
5626 
5627 #endif // ifndef DOXYGEN
5628 
5629 
5630 
5632 
5633 #endif
void resize(const size_type new_size)
Number local_element(const size_type local_index) const
void reinit(const size_type size, const bool omit_zeroing_entries=false)
Abstract base class for mapping classes.
Definition: mapping.h:317
unsigned int n_active_fe_indices() const
void loop_cell_centric(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > mapping_info
Definition: matrix_free.h:2286
unsigned int get_cell_active_fe_index(const std::pair< unsigned int, unsigned int > range) const
AlignedVector< Number > * acquire_scratch_data_non_threadsafe() const
unsigned int n_ghost_cell_batches() const
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void initialize_dof_handlers(const std::vector< const DoFHandler< dim, dim > * > &dof_handlers, const AdditionalData &additional_data)
types::boundary_id get_boundary_id(const unsigned int face_batch_index) const
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > > shape_info
Definition: matrix_free.h:2292
const Quadrature< dim > & get_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
std::pair< unsigned int, unsigned int > create_cell_subrange_hp(const std::pair< unsigned int, unsigned int > &range, const unsigned int fe_degree, const unsigned int dof_handler_index=0) const
const IndexSet & get_ghost_set(const unsigned int dof_handler_index=0) const
void print(std::ostream &out) const
const Quadrature< dim - 1 > & get_face_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int n_inner_face_batches() const
void internal_reinit(const std::shared_ptr< hp::MappingCollection< dim >> &mapping, const std::vector< const DoFHandler< dim, dim > * > &dof_handlers, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const std::vector< hp::QCollection< q_dim >> &quad, const AdditionalData &additional_data)
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
void update_mapping(const Mapping< dim > &mapping)
unsigned int get_mg_level() const
void print_memory_consumption(StreamType &out) const
void update_mapping(const std::shared_ptr< hp::MappingCollection< dim >> &mapping)
bool mapping_initialized() const
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
const Number * constraint_pool_begin(const unsigned int pool_index) const
void clear()
unsigned int get_n_q_points_face(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
std::pair< unsigned int, unsigned int > create_cell_subrange_hp_by_index(const std::pair< unsigned int, unsigned int > &range, const unsigned int fe_index, const unsigned int dof_handler_index=0) const
bool at_irregular_cell(const unsigned int cell_batch_index) const
~MatrixFree() override=default
void copy_from(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free_base)
void initialize_face_data_vector(AlignedVector< T > &vec) const
unsigned int get_dofs_per_cell(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int n_constraint_pool_entries() const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
std::pair< int, int > get_cell_level_and_index(const unsigned int cell_batch_index, const unsigned int lane_index) const
std::pair< unsigned int, unsigned int > get_face_range_category(const std::pair< unsigned int, unsigned int > face_batch_range) const
const std::vector< unsigned int > & get_constrained_dofs(const unsigned int dof_handler_index=0) const
std::pair< unsigned int, unsigned int > get_face_category(const unsigned int face_batch_index) const
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
unsigned int get_dofs_per_face(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
const Number * constraint_pool_end(const unsigned int pool_index) const
internal::MatrixFreeFunctions::TaskInfo task_info
Definition: matrix_free.h:2321
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*boundary_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
std::pair< typename DoFHandler< dim >::cell_iterator, unsigned int > get_face_iterator(const unsigned int face_batch_index, const unsigned int lane_index, const bool interior=true, const unsigned int fe_component=0) const
bool mapping_is_initialized
Definition: matrix_free.h:2338
void loop_cell_centric(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void initialize_dof_vector(LinearAlgebra::distributed::Vector< Number2 > &vec, const unsigned int dof_handler_index=0) const
MatrixFree(const MatrixFree< dim, Number, VectorizedArrayType > &other)
unsigned int mg_level
Definition: matrix_free.h:2361
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
std::vector< std::pair< unsigned int, unsigned int > > cell_level_index
Definition: matrix_free.h:2300
std::array< types::boundary_id, VectorizedArrayType::size()> get_faces_by_cells_boundary_id(const unsigned int cell_batch_index, const unsigned int face_number) const
std::vector< Number > constraint_pool_data
Definition: matrix_free.h:2273
VectorizedArrayType vectorized_value_type
Definition: matrix_free.h:125
unsigned int n_cell_batches() const
bool indices_initialized() const
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
unsigned int get_cell_category(const unsigned int cell_batch_index) const
Threads::ThreadLocalStorage< std::list< std::pair< bool, AlignedVector< VectorizedArrayType > > > > scratch_pad
Definition: matrix_free.h:2349
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_index) const
unsigned int get_matrix_free_cell_index(const typename Triangulation< dim >::cell_iterator &cell) const
Number value_type
Definition: matrix_free.h:124
std::size_t memory_consumption() const
const AffineConstraints< Number > & get_affine_constraints(const unsigned int dof_handler_index=0) const
unsigned int n_boundary_face_batches() const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
std::list< std::pair< bool, AlignedVector< Number > > > scratch_pad_non_threadsafe
Definition: matrix_free.h:2356
DoFHandler< dim >::cell_iterator get_cell_iterator(const unsigned int cell_batch_index, const unsigned int lane_index, const unsigned int dof_handler_index=0) const
void loop_cell_centric(const std::function< void(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
std::vector< SmartPointer< const DoFHandler< dim > > > dof_handlers
Definition: matrix_free.h:2251
void initialize_indices(const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const AdditionalData &additional_data)
std::vector< SmartPointer< const AffineConstraints< Number > > > affine_constraints
Definition: matrix_free.h:2259
unsigned int n_components() const
unsigned int n_physical_cells() const
void initialize_dof_vector(VectorType &vec, const unsigned int dof_handler_index=0) const
unsigned int get_n_q_points(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
void reinit(const MappingType &mapping, const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
std::vector< unsigned int > mf_cell_indices
Definition: matrix_free.h:2306
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*boundary_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
std::vector< internal::MatrixFreeFunctions::DoFInfo > dof_info
Definition: matrix_free.h:2265
void initialize_cell_data_vector(AlignedVector< T > &vec) const
unsigned int n_ghost_inner_face_batches() const
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
void release_scratch_data_non_threadsafe(const AlignedVector< Number > *memory) const
void renumber_dofs(std::vector< types::global_dof_index > &renumbering, const unsigned int dof_handler_index=0)
unsigned int get_face_active_fe_index(const std::pair< unsigned int, unsigned int > range, const bool is_interior_face=true) const
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
bool indices_are_initialized
Definition: matrix_free.h:2333
internal::MatrixFreeFunctions::FaceInfo< VectorizedArrayType::size()> face_info
Definition: matrix_free.h:2328
void reinit(const MappingType &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< number2 > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
unsigned int get_cell_range_category(const std::pair< unsigned int, unsigned int > cell_batch_range) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*boundary_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void reinit(const MappingType &mapping, const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< QuadratureType > &quad, const AdditionalData &additional_data=AdditionalData())
std::vector< unsigned int > constraint_pool_row_index
Definition: matrix_free.h:2279
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner(const unsigned int dof_handler_index=0) const
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_operation, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
unsigned int cell_level_index_end_local
Definition: matrix_free.h:2314
unsigned int n_base_elements(const unsigned int dof_handler_index) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*boundary_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
const IndexSet & get_locally_owned_set(const unsigned int dof_handler_index=0) const
static bool is_supported(const FiniteElement< dim, spacedim > &fe)
static constexpr unsigned int dimension
Definition: matrix_free.h:130
A class that provides a separate storage location on each thread that accesses the object.
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
Point< 2 > second
Definition: grid_out.cc:4615
Point< 2 > first
Definition: grid_out.cc:4614
unsigned int level
Definition: grid_out.cc:4617
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1916
#define AssertIndexRange(index, range)
Definition: exceptions.h:1857
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:469
UpdateFlags
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_default
No update.
static const char U
VectorType::value_type * begin(VectorType &V)
bool job_supports_mpi()
Definition: mpi.cc:1072
unsigned int minimum_parallel_grain_size
Definition: parallel.cc:34
const types::boundary_id invalid_boundary_id
Definition: types.h:261
static const unsigned int invalid_unsigned_int
Definition: types.h:213
void apply_to_subranges(const tbb::blocked_range< Iterator > &range, const Function &f)
Definition: parallel.h:371
unsigned short int fe_index
Definition: types.h:60
unsigned int boundary_id
Definition: types.h:141
TasksParallelScheme tasks_parallel_scheme
Definition: matrix_free.h:348
UpdateFlags mapping_update_flags_inner_faces
Definition: matrix_free.h:414
AdditionalData(const TasksParallelScheme tasks_parallel_scheme=partition_partition, const unsigned int tasks_block_size=0, const UpdateFlags mapping_update_flags=update_gradients|update_JxW_values, const UpdateFlags mapping_update_flags_boundary_faces=update_default, const UpdateFlags mapping_update_flags_inner_faces=update_default, const UpdateFlags mapping_update_flags_faces_by_cells=update_default, const unsigned int mg_level=numbers::invalid_unsigned_int, const bool store_plain_indices=true, const bool initialize_indices=true, const bool initialize_mapping=true, const bool overlap_communication_computation=true, const bool hold_all_faces_to_owned_cells=false, const bool cell_vectorization_categories_strict=false, const bool allow_ghosted_vectors_in_loops=true)
Definition: matrix_free.h:222
AdditionalData & operator=(const AdditionalData &other)
Definition: matrix_free.h:286
std::vector< unsigned int > cell_vectorization_category
Definition: matrix_free.h:526
UpdateFlags mapping_update_flags_boundary_faces
Definition: matrix_free.h:394
UpdateFlags mapping_update_flags
Definition: matrix_free.h:374
UpdateFlags mapping_update_flags_faces_by_cells
Definition: matrix_free.h:442
AdditionalData(const AdditionalData &other)
Definition: matrix_free.h:259
unsigned int tasks_block_size
Definition: matrix_free.h:359
static bool equal(const T *p1, const T *p2)
std::vector< unsigned int > cell_loop_pre_list_index
Definition: dof_info.h:749
std::vector< unsigned int > cell_loop_post_list_index
Definition: dof_info.h:762
std::vector< std::pair< unsigned int, unsigned int > > vector_zero_range_list
Definition: dof_info.h:742
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition: dof_info.h:593
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_pre_list
Definition: dof_info.h:755
std::vector< unsigned int > vector_zero_range_list_index
Definition: dof_info.h:737
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_post_list
Definition: dof_info.h:768
::Table< 3, unsigned int > cell_and_face_to_plain_faces
Definition: face_info.h:173
std::vector< FaceToCellTopology< vectorization_width > > faces
Definition: face_info.h:165
::Table< 3, types::boundary_id > cell_and_face_boundary_id
Definition: face_info.h:179
std::vector< unsigned int > boundary_partition_data
Definition: task_info.h:519
void loop(MFWorkerInterface &worker) const
Definition: task_info.cc:351
std::vector< unsigned int > cell_partition_data
Definition: task_info.h:475
std::vector< unsigned int > face_partition_data
Definition: task_info.h:497