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