Reference documentation for deal.II version GIT 05ffa62ef0 2022-11-26 15:30:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
matrix_free.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2011 - 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 
178  {
183 
190  {
210  };
211 
217  const unsigned int tasks_block_size = 0,
223  const unsigned int mg_level = numbers::invalid_unsigned_int,
224  const bool store_plain_indices = true,
225  const bool initialize_indices = true,
226  const bool initialize_mapping = true,
227  const bool overlap_communication_computation = true,
228  const bool hold_all_faces_to_owned_cells = false,
229  const bool cell_vectorization_categories_strict = false,
230  const bool allow_ghosted_vectors_in_loops = true)
237  , mg_level(mg_level)
246  , communicator_sm(MPI_COMM_SELF)
247  {}
248 
261  , mg_level(other.mg_level)
273  {}
274 
279  operator=(const AdditionalData &other)
280  {
289  mg_level = other.mg_level;
301 
302  return *this;
303  }
304 
342 
352  unsigned int tasks_block_size;
353 
366 
386 
406 
434 
443  unsigned int mg_level;
444 
452 
463 
472 
486 
495 
517  std::vector<unsigned int> cell_vectorization_category;
518 
529 
541 
545  MPI_Comm communicator_sm;
546  };
547 
556 
561 
565  ~MatrixFree() override = default;
566 
579  template <typename QuadratureType, typename number2, typename MappingType>
580  void
581  reinit(const MappingType & mapping,
582  const DoFHandler<dim> & dof_handler,
583  const AffineConstraints<number2> &constraint,
584  const QuadratureType & quad,
585  const AdditionalData & additional_data = AdditionalData());
586 
608  template <typename QuadratureType, typename number2, typename MappingType>
609  void
610  reinit(const MappingType & mapping,
611  const std::vector<const DoFHandler<dim> *> & dof_handler,
612  const std::vector<const AffineConstraints<number2> *> &constraint,
613  const std::vector<QuadratureType> & quad,
614  const AdditionalData &additional_data = AdditionalData());
615 
623  template <typename QuadratureType, typename number2, typename MappingType>
624  void
625  reinit(const MappingType & mapping,
626  const std::vector<const DoFHandler<dim> *> & dof_handler,
627  const std::vector<const AffineConstraints<number2> *> &constraint,
628  const QuadratureType & quad,
629  const AdditionalData &additional_data = AdditionalData());
630 
636  void
638  const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free_base);
639 
649  void
650  update_mapping(const Mapping<dim> &mapping);
651 
655  void
656  update_mapping(const std::shared_ptr<hp::MappingCollection<dim>> &mapping);
657 
662  void
663  clear();
664 
678  enum class DataAccessOnFaces
679  {
686  none,
687 
698  values,
699 
709 
720  gradients,
721 
731 
739  };
740 
785  template <typename OutVector, typename InVector>
786  void
787  cell_loop(const std::function<void(
789  OutVector &,
790  const InVector &,
791  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
792  OutVector & dst,
793  const InVector & src,
794  const bool zero_dst_vector = false) const;
795 
842  template <typename CLASS, typename OutVector, typename InVector>
843  void
844  cell_loop(void (CLASS::*cell_operation)(
845  const MatrixFree &,
846  OutVector &,
847  const InVector &,
848  const std::pair<unsigned int, unsigned int> &) const,
849  const CLASS * owning_class,
850  OutVector & dst,
851  const InVector &src,
852  const bool zero_dst_vector = false) const;
853 
857  template <typename CLASS, typename OutVector, typename InVector>
858  void
859  cell_loop(void (CLASS::*cell_operation)(
860  const MatrixFree &,
861  OutVector &,
862  const InVector &,
863  const std::pair<unsigned int, unsigned int> &),
864  CLASS * owning_class,
865  OutVector & dst,
866  const InVector &src,
867  const bool zero_dst_vector = false) const;
868 
953  template <typename CLASS, typename OutVector, typename InVector>
954  void
955  cell_loop(void (CLASS::*cell_operation)(
956  const MatrixFree &,
957  OutVector &,
958  const InVector &,
959  const std::pair<unsigned int, unsigned int> &) const,
960  const CLASS * owning_class,
961  OutVector & dst,
962  const InVector &src,
963  const std::function<void(const unsigned int, const unsigned int)>
964  &operation_before_loop,
965  const std::function<void(const unsigned int, const unsigned int)>
966  & operation_after_loop,
967  const unsigned int dof_handler_index_pre_post = 0) const;
968 
972  template <typename CLASS, typename OutVector, typename InVector>
973  void
974  cell_loop(void (CLASS::*cell_operation)(
975  const MatrixFree &,
976  OutVector &,
977  const InVector &,
978  const std::pair<unsigned int, unsigned int> &),
979  CLASS * owning_class,
980  OutVector & dst,
981  const InVector &src,
982  const std::function<void(const unsigned int, const unsigned int)>
983  &operation_before_loop,
984  const std::function<void(const unsigned int, const unsigned int)>
985  & operation_after_loop,
986  const unsigned int dof_handler_index_pre_post = 0) const;
987 
992  template <typename OutVector, typename InVector>
993  void
994  cell_loop(const std::function<void(
996  OutVector &,
997  const InVector &,
998  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
999  OutVector & dst,
1000  const InVector & src,
1001  const std::function<void(const unsigned int, const unsigned int)>
1002  &operation_before_loop,
1003  const std::function<void(const unsigned int, const unsigned int)>
1004  & operation_after_loop,
1005  const unsigned int dof_handler_index_pre_post = 0) const;
1006 
1082  template <typename OutVector, typename InVector>
1083  void
1084  loop(const std::function<
1086  OutVector &,
1087  const InVector &,
1088  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1089  const std::function<
1091  OutVector &,
1092  const InVector &,
1093  const std::pair<unsigned int, unsigned int> &)> &face_operation,
1094  const std::function<void(
1096  OutVector &,
1097  const InVector &,
1098  const std::pair<unsigned int, unsigned int> &)> &boundary_operation,
1099  OutVector & dst,
1100  const InVector & src,
1101  const bool zero_dst_vector = false,
1102  const DataAccessOnFaces dst_vector_face_access =
1104  const DataAccessOnFaces src_vector_face_access =
1106 
1192  template <typename CLASS, typename OutVector, typename InVector>
1193  void
1195  void (CLASS::*cell_operation)(const MatrixFree &,
1196  OutVector &,
1197  const InVector &,
1198  const std::pair<unsigned int, unsigned int> &)
1199  const,
1200  void (CLASS::*face_operation)(const MatrixFree &,
1201  OutVector &,
1202  const InVector &,
1203  const std::pair<unsigned int, unsigned int> &)
1204  const,
1205  void (CLASS::*boundary_operation)(
1206  const MatrixFree &,
1207  OutVector &,
1208  const InVector &,
1209  const std::pair<unsigned int, unsigned int> &) const,
1210  const CLASS * owning_class,
1211  OutVector & dst,
1212  const InVector & src,
1213  const bool zero_dst_vector = false,
1214  const DataAccessOnFaces dst_vector_face_access =
1216  const DataAccessOnFaces src_vector_face_access =
1218 
1222  template <typename CLASS, typename OutVector, typename InVector>
1223  void
1224  loop(void (CLASS::*cell_operation)(
1225  const MatrixFree &,
1226  OutVector &,
1227  const InVector &,
1228  const std::pair<unsigned int, unsigned int> &),
1229  void (CLASS::*face_operation)(
1230  const MatrixFree &,
1231  OutVector &,
1232  const InVector &,
1233  const std::pair<unsigned int, unsigned int> &),
1234  void (CLASS::*boundary_operation)(
1235  const MatrixFree &,
1236  OutVector &,
1237  const InVector &,
1238  const std::pair<unsigned int, unsigned int> &),
1239  CLASS * owning_class,
1240  OutVector & dst,
1241  const InVector & src,
1242  const bool zero_dst_vector = false,
1243  const DataAccessOnFaces dst_vector_face_access =
1245  const DataAccessOnFaces src_vector_face_access =
1247 
1367  template <typename CLASS, typename OutVector, typename InVector>
1368  void
1370  void (CLASS::*cell_operation)(const MatrixFree &,
1371  OutVector &,
1372  const InVector &,
1373  const std::pair<unsigned int, unsigned int> &)
1374  const,
1375  void (CLASS::*face_operation)(const MatrixFree &,
1376  OutVector &,
1377  const InVector &,
1378  const std::pair<unsigned int, unsigned int> &)
1379  const,
1380  void (CLASS::*boundary_operation)(
1381  const MatrixFree &,
1382  OutVector &,
1383  const InVector &,
1384  const std::pair<unsigned int, unsigned int> &) const,
1385  const CLASS * owning_class,
1386  OutVector & dst,
1387  const InVector &src,
1388  const std::function<void(const unsigned int, const unsigned int)>
1389  &operation_before_loop,
1390  const std::function<void(const unsigned int, const unsigned int)>
1391  & operation_after_loop,
1392  const unsigned int dof_handler_index_pre_post = 0,
1393  const DataAccessOnFaces dst_vector_face_access =
1395  const DataAccessOnFaces src_vector_face_access =
1397 
1401  template <typename CLASS, typename OutVector, typename InVector>
1402  void
1403  loop(void (CLASS::*cell_operation)(
1404  const MatrixFree &,
1405  OutVector &,
1406  const InVector &,
1407  const std::pair<unsigned int, unsigned int> &),
1408  void (CLASS::*face_operation)(
1409  const MatrixFree &,
1410  OutVector &,
1411  const InVector &,
1412  const std::pair<unsigned int, unsigned int> &),
1413  void (CLASS::*boundary_operation)(
1414  const MatrixFree &,
1415  OutVector &,
1416  const InVector &,
1417  const std::pair<unsigned int, unsigned int> &),
1418  const CLASS * owning_class,
1419  OutVector & dst,
1420  const InVector &src,
1421  const std::function<void(const unsigned int, const unsigned int)>
1422  &operation_before_loop,
1423  const std::function<void(const unsigned int, const unsigned int)>
1424  & operation_after_loop,
1425  const unsigned int dof_handler_index_pre_post = 0,
1426  const DataAccessOnFaces dst_vector_face_access =
1428  const DataAccessOnFaces src_vector_face_access =
1430 
1436  template <typename OutVector, typename InVector>
1437  void
1438  loop(const std::function<
1440  OutVector &,
1441  const InVector &,
1442  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1443  const std::function<
1445  OutVector &,
1446  const InVector &,
1447  const std::pair<unsigned int, unsigned int> &)> &face_operation,
1448  const std::function<void(
1450  OutVector &,
1451  const InVector &,
1452  const std::pair<unsigned int, unsigned int> &)> &boundary_operation,
1453  OutVector & dst,
1454  const InVector & src,
1455  const std::function<void(const unsigned int, const unsigned int)>
1456  &operation_before_loop,
1457  const std::function<void(const unsigned int, const unsigned int)>
1458  & operation_after_loop,
1459  const unsigned int dof_handler_index_pre_post = 0,
1460  const DataAccessOnFaces dst_vector_face_access =
1462  const DataAccessOnFaces src_vector_face_access =
1464 
1529  template <typename CLASS, typename OutVector, typename InVector>
1530  void
1531  loop_cell_centric(void (CLASS::*cell_operation)(
1532  const MatrixFree &,
1533  OutVector &,
1534  const InVector &,
1535  const std::pair<unsigned int, unsigned int> &) const,
1536  const CLASS * owning_class,
1537  OutVector & dst,
1538  const InVector & src,
1539  const bool zero_dst_vector = false,
1540  const DataAccessOnFaces src_vector_face_access =
1542 
1546  template <typename CLASS, typename OutVector, typename InVector>
1547  void
1548  loop_cell_centric(void (CLASS::*cell_operation)(
1549  const MatrixFree &,
1550  OutVector &,
1551  const InVector &,
1552  const std::pair<unsigned int, unsigned int> &),
1553  CLASS * owning_class,
1554  OutVector & dst,
1555  const InVector & src,
1556  const bool zero_dst_vector = false,
1557  const DataAccessOnFaces src_vector_face_access =
1559 
1563  template <typename OutVector, typename InVector>
1564  void
1566  const std::function<void(const MatrixFree &,
1567  OutVector &,
1568  const InVector &,
1569  const std::pair<unsigned int, unsigned int> &)>
1570  & cell_operation,
1571  OutVector & dst,
1572  const InVector & src,
1573  const bool zero_dst_vector = false,
1574  const DataAccessOnFaces src_vector_face_access =
1576 
1584  std::pair<unsigned int, unsigned int>
1585  create_cell_subrange_hp(const std::pair<unsigned int, unsigned int> &range,
1586  const unsigned int fe_degree,
1587  const unsigned int dof_handler_index = 0) const;
1588 
1595  std::pair<unsigned int, unsigned int>
1597  const std::pair<unsigned int, unsigned int> &range,
1598  const unsigned int fe_index,
1599  const unsigned int dof_handler_index = 0) const;
1600 
1604  unsigned int
1606 
1610  unsigned int
1612  const std::pair<unsigned int, unsigned int> range) const;
1613 
1617  unsigned int
1618  get_face_active_fe_index(const std::pair<unsigned int, unsigned int> range,
1619  const bool is_interior_face = true) const;
1620 
1632  template <typename T>
1633  void
1635 
1641  template <typename T>
1642  void
1644 
1658  template <typename VectorType>
1659  void
1660  initialize_dof_vector(VectorType & vec,
1661  const unsigned int dof_handler_index = 0) const;
1662 
1680  template <typename Number2>
1681  void
1683  const unsigned int dof_handler_index = 0) const;
1684 
1695  const std::shared_ptr<const Utilities::MPI::Partitioner> &
1696  get_vector_partitioner(const unsigned int dof_handler_index = 0) const;
1697 
1701  const IndexSet &
1702  get_locally_owned_set(const unsigned int dof_handler_index = 0) const;
1703 
1707  const IndexSet &
1708  get_ghost_set(const unsigned int dof_handler_index = 0) const;
1709 
1719  const std::vector<unsigned int> &
1720  get_constrained_dofs(const unsigned int dof_handler_index = 0) const;
1721 
1732  void
1733  renumber_dofs(std::vector<types::global_dof_index> &renumbering,
1734  const unsigned int dof_handler_index = 0);
1735 
1745  template <int spacedim>
1746  static bool
1748 
1752  unsigned int
1753  n_components() const;
1754 
1759  unsigned int
1760  n_base_elements(const unsigned int dof_handler_index) const;
1761 
1769  unsigned int
1771 
1781  unsigned int
1783 
1790  unsigned int
1792 
1802  unsigned int
1804 
1815  unsigned int
1817 
1822  unsigned int
1824 
1836  get_boundary_id(const unsigned int face_batch_index) const;
1837 
1842  std::array<types::boundary_id, VectorizedArrayType::size()>
1843  get_faces_by_cells_boundary_id(const unsigned int cell_batch_index,
1844  const unsigned int face_number) const;
1845 
1850  const DoFHandler<dim> &
1851  get_dof_handler(const unsigned int dof_handler_index = 0) const;
1852 
1860  get_affine_constraints(const unsigned int dof_handler_index = 0) const;
1861 
1875  get_cell_iterator(const unsigned int cell_batch_index,
1876  const unsigned int lane_index,
1877  const unsigned int dof_handler_index = 0) const;
1878 
1884  std::pair<int, int>
1885  get_cell_level_and_index(const unsigned int cell_batch_index,
1886  const unsigned int lane_index) const;
1887 
1894  unsigned int
1896  const typename Triangulation<dim>::cell_iterator &cell) const;
1897 
1910  std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>
1911  get_face_iterator(const unsigned int face_batch_index,
1912  const unsigned int lane_index,
1913  const bool interior = true,
1914  const unsigned int fe_component = 0) const;
1915 
1928  bool
1929  at_irregular_cell(const unsigned int cell_batch_index) const;
1930 
1940  unsigned int
1941  n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const;
1942 
1952  unsigned int
1953  n_active_entries_per_face_batch(const unsigned int face_batch_index) const;
1954 
1958  unsigned int
1959  get_dofs_per_cell(const unsigned int dof_handler_index = 0,
1960  const unsigned int hp_active_fe_index = 0) const;
1961 
1965  unsigned int
1966  get_n_q_points(const unsigned int quad_index = 0,
1967  const unsigned int hp_active_fe_index = 0) const;
1968 
1973  unsigned int
1974  get_dofs_per_face(const unsigned int dof_handler_index = 0,
1975  const unsigned int hp_active_fe_index = 0) const;
1976 
1981  unsigned int
1982  get_n_q_points_face(const unsigned int quad_index = 0,
1983  const unsigned int hp_active_fe_index = 0) const;
1984 
1988  const Quadrature<dim> &
1989  get_quadrature(const unsigned int quad_index = 0,
1990  const unsigned int hp_active_fe_index = 0) const;
1991 
1995  const Quadrature<dim - 1> &
1996  get_face_quadrature(const unsigned int quad_index = 0,
1997  const unsigned int hp_active_fe_index = 0) const;
1998 
2011  unsigned int
2013  const std::pair<unsigned int, unsigned int> cell_batch_range) const;
2014 
2019  std::pair<unsigned int, unsigned int>
2021  const std::pair<unsigned int, unsigned int> face_batch_range) const;
2022 
2034  unsigned int
2035  get_cell_category(const unsigned int cell_batch_index) const;
2036 
2041  std::pair<unsigned int, unsigned int>
2042  get_face_category(const unsigned int face_batch_index) const;
2043 
2047  bool
2049 
2054  bool
2056 
2061  unsigned int
2062  get_mg_level() const;
2063 
2068  std::size_t
2070 
2075  template <typename StreamType>
2076  void
2077  print_memory_consumption(StreamType &out) const;
2078 
2083  void
2084  print(std::ostream &out) const;
2085 
2098  get_task_info() const;
2099 
2100  /*
2101  * Return geometry-dependent information on the cells.
2102  */
2103  const internal::MatrixFreeFunctions::
2104  MappingInfo<dim, Number, VectorizedArrayType> &
2106 
2111  get_dof_info(const unsigned int dof_handler_index_component = 0) const;
2112 
2116  unsigned int
2118 
2123  const Number *
2124  constraint_pool_begin(const unsigned int pool_index) const;
2125 
2131  const Number *
2132  constraint_pool_end(const unsigned int pool_index) const;
2133 
2138  get_shape_info(const unsigned int dof_handler_index_component = 0,
2139  const unsigned int quad_index = 0,
2140  const unsigned int fe_base_element = 0,
2141  const unsigned int hp_active_fe_index = 0,
2142  const unsigned int hp_active_quad_index = 0) const;
2143 
2148  VectorizedArrayType::size()> &
2149  get_face_info(const unsigned int face_batch_index) const;
2150 
2151 
2157  const Table<3, unsigned int> &
2159 
2175 
2179  void
2181 
2193 
2197  void
2199  const AlignedVector<Number> *memory) const;
2200 
2203 private:
2208  template <typename number2, int q_dim>
2209  void
2211  const std::shared_ptr<hp::MappingCollection<dim>> & mapping,
2212  const std::vector<const DoFHandler<dim, dim> *> & dof_handlers,
2213  const std::vector<const AffineConstraints<number2> *> &constraint,
2214  const std::vector<IndexSet> & locally_owned_set,
2215  const std::vector<hp::QCollection<q_dim>> & quad,
2216  const AdditionalData & additional_data);
2217 
2224  template <typename number2>
2225  void
2227  const std::vector<const AffineConstraints<number2> *> &constraint,
2228  const std::vector<IndexSet> & locally_owned_set,
2229  const AdditionalData & additional_data);
2230 
2234  void
2236  const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
2237  const AdditionalData & additional_data);
2238 
2242  std::vector<SmartPointer<const DoFHandler<dim>>> dof_handlers;
2243 
2250  std::vector<SmartPointer<const AffineConstraints<Number>>> affine_constraints;
2251 
2256  std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
2257 
2264  std::vector<Number> constraint_pool_data;
2265 
2270  std::vector<unsigned int> constraint_pool_row_index;
2271 
2278 
2284 
2291  std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
2292 
2297  std::vector<unsigned int> mf_cell_indices;
2298 
2306 
2313 
2318  internal::MatrixFreeFunctions::FaceInfo<VectorizedArrayType::size()>
2320 
2325 
2330 
2339  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2341 
2346  mutable std::list<std::pair<bool, AlignedVector<Number>>>
2348 
2352  unsigned int mg_level;
2353 };
2354 
2355 
2356 
2357 /*----------------------- Inline functions ----------------------------------*/
2358 
2359 #ifndef DOXYGEN
2360 
2361 
2362 
2363 template <int dim, typename Number, typename VectorizedArrayType>
2364 template <typename T>
2365 inline void
2367  AlignedVector<T> &vec) const
2368 {
2369  vec.resize(this->n_cell_batches() + this->n_ghost_cell_batches());
2370 }
2371 
2372 
2373 
2374 template <int dim, typename Number, typename VectorizedArrayType>
2375 template <typename T>
2376 inline void
2378  AlignedVector<T> &vec) const
2379 {
2380  vec.resize(this->n_inner_face_batches() + this->n_boundary_face_batches() +
2381  this->n_ghost_inner_face_batches());
2382 }
2383 
2384 
2385 
2386 template <int dim, typename Number, typename VectorizedArrayType>
2387 template <typename VectorType>
2388 inline void
2390  VectorType & vec,
2391  const unsigned int comp) const
2392 {
2393  static_assert(IsBlockVector<VectorType>::value == false,
2394  "This function is not supported for block vectors.");
2395 
2396  Assert(task_info.n_procs == 1,
2397  ExcMessage("This function can only be used in serial."));
2398 
2399  AssertIndexRange(comp, n_components());
2400  vec.reinit(dof_info[comp].vector_partitioner->size());
2401 }
2402 
2403 
2404 
2405 template <int dim, typename Number, typename VectorizedArrayType>
2406 template <typename Number2>
2407 inline void
2410  const unsigned int comp) const
2411 {
2412  AssertIndexRange(comp, n_components());
2413  vec.reinit(dof_info[comp].vector_partitioner, task_info.communicator_sm);
2414 }
2415 
2416 
2417 
2418 template <int dim, typename Number, typename VectorizedArrayType>
2419 inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2421  const unsigned int comp) const
2422 {
2423  AssertIndexRange(comp, n_components());
2424  return dof_info[comp].vector_partitioner;
2425 }
2426 
2427 
2428 
2429 template <int dim, typename Number, typename VectorizedArrayType>
2430 inline const std::vector<unsigned int> &
2432  const unsigned int comp) const
2433 {
2434  AssertIndexRange(comp, n_components());
2435  return dof_info[comp].constrained_dofs;
2436 }
2437 
2438 
2439 
2440 template <int dim, typename Number, typename VectorizedArrayType>
2441 inline unsigned int
2443 {
2444  AssertDimension(dof_handlers.size(), dof_info.size());
2445  return dof_handlers.size();
2446 }
2447 
2448 
2449 
2450 template <int dim, typename Number, typename VectorizedArrayType>
2451 inline unsigned int
2453  const unsigned int dof_no) const
2454 {
2455  AssertDimension(dof_handlers.size(), dof_info.size());
2456  AssertIndexRange(dof_no, dof_handlers.size());
2457  return dof_handlers[dof_no]->get_fe().n_base_elements();
2458 }
2459 
2460 
2461 
2462 template <int dim, typename Number, typename VectorizedArrayType>
2465 {
2466  return task_info;
2467 }
2468 
2469 
2470 
2471 template <int dim, typename Number, typename VectorizedArrayType>
2472 inline unsigned int
2474 {
2475  return task_info.n_active_cells;
2476 }
2477 
2478 
2479 
2480 template <int dim, typename Number, typename VectorizedArrayType>
2481 inline unsigned int
2483 {
2484  return *(task_info.cell_partition_data.end() - 2);
2485 }
2486 
2487 
2488 
2489 template <int dim, typename Number, typename VectorizedArrayType>
2490 inline unsigned int
2492 {
2493  return *(task_info.cell_partition_data.end() - 1) -
2494  *(task_info.cell_partition_data.end() - 2);
2495 }
2496 
2497 
2498 
2499 template <int dim, typename Number, typename VectorizedArrayType>
2500 inline unsigned int
2502 {
2503  if (task_info.face_partition_data.size() == 0)
2504  return 0;
2505  return task_info.face_partition_data.back();
2506 }
2507 
2508 
2509 
2510 template <int dim, typename Number, typename VectorizedArrayType>
2511 inline unsigned int
2513 {
2514  if (task_info.face_partition_data.size() == 0)
2515  return 0;
2516  return task_info.boundary_partition_data.back() -
2518 }
2519 
2520 
2521 
2522 template <int dim, typename Number, typename VectorizedArrayType>
2523 inline unsigned int
2525 {
2526  if (task_info.face_partition_data.size() == 0)
2527  return 0;
2528  return face_info.faces.size() - task_info.boundary_partition_data.back();
2529 }
2530 
2531 
2532 
2533 template <int dim, typename Number, typename VectorizedArrayType>
2534 inline types::boundary_id
2536  const unsigned int face_batch_index) const
2537 {
2538  Assert(face_batch_index >= task_info.boundary_partition_data[0] &&
2539  face_batch_index < task_info.boundary_partition_data.back(),
2540  ExcIndexRange(face_batch_index,
2543  return types::boundary_id(face_info.faces[face_batch_index].exterior_face_no);
2544 }
2545 
2546 
2547 
2548 template <int dim, typename Number, typename VectorizedArrayType>
2549 inline std::array<types::boundary_id, VectorizedArrayType::size()>
2551  const unsigned int cell_batch_index,
2552  const unsigned int face_number) const
2553 {
2554  AssertIndexRange(cell_batch_index, n_cell_batches());
2557  ExcNotInitialized());
2558  std::array<types::boundary_id, VectorizedArrayType::size()> result;
2559  result.fill(numbers::invalid_boundary_id);
2560  for (unsigned int v = 0;
2561  v < n_active_entries_per_cell_batch(cell_batch_index);
2562  ++v)
2563  result[v] =
2564  face_info.cell_and_face_boundary_id(cell_batch_index, face_number, v);
2565  return result;
2566 }
2567 
2568 
2569 
2570 template <int dim, typename Number, typename VectorizedArrayType>
2571 inline const internal::MatrixFreeFunctions::
2572  MappingInfo<dim, Number, VectorizedArrayType> &
2574 {
2575  return mapping_info;
2576 }
2577 
2578 
2579 
2580 template <int dim, typename Number, typename VectorizedArrayType>
2583  const unsigned int dof_index) const
2584 {
2585  AssertIndexRange(dof_index, n_components());
2586  return dof_info[dof_index];
2587 }
2588 
2589 
2590 
2591 template <int dim, typename Number, typename VectorizedArrayType>
2592 inline unsigned int
2594 {
2595  return constraint_pool_row_index.size() - 1;
2596 }
2597 
2598 
2599 
2600 template <int dim, typename Number, typename VectorizedArrayType>
2601 inline const Number *
2603  const unsigned int row) const
2604 {
2606  return constraint_pool_data.empty() ?
2607  nullptr :
2609 }
2610 
2611 
2612 
2613 template <int dim, typename Number, typename VectorizedArrayType>
2614 inline const Number *
2616  const unsigned int row) const
2617 {
2619  return constraint_pool_data.empty() ?
2620  nullptr :
2622 }
2623 
2624 
2625 
2626 template <int dim, typename Number, typename VectorizedArrayType>
2627 inline std::pair<unsigned int, unsigned int>
2629  const std::pair<unsigned int, unsigned int> &range,
2630  const unsigned int degree,
2631  const unsigned int dof_handler_component) const
2632 {
2633  if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2634  {
2636  dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2638  dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2639  if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2640  return range;
2641  else
2642  return {range.second, range.second};
2643  }
2644 
2645  const unsigned int fe_index =
2646  dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2647  if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2648  return {range.second, range.second};
2649  else
2650  return create_cell_subrange_hp_by_index(range,
2651  fe_index,
2652  dof_handler_component);
2653 }
2654 
2655 
2656 
2657 template <int dim, typename Number, typename VectorizedArrayType>
2658 inline bool
2660  const unsigned int cell_batch_index) const
2661 {
2662  AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2663  return VectorizedArrayType::size() > 1 &&
2664  cell_level_index[(cell_batch_index + 1) * VectorizedArrayType::size() -
2665  1] == cell_level_index[(cell_batch_index + 1) *
2666  VectorizedArrayType::size() -
2667  2];
2668 }
2669 
2670 
2671 
2672 template <int dim, typename Number, typename VectorizedArrayType>
2673 unsigned int
2675 {
2676  return shape_info.size(2);
2677 }
2678 
2679 
2680 template <int dim, typename Number, typename VectorizedArrayType>
2681 unsigned int
2683  const std::pair<unsigned int, unsigned int> range) const
2684 {
2685  const auto &fe_indices = dof_info[0].cell_active_fe_index;
2686 
2687  if (fe_indices.empty() == true)
2688  return 0;
2689 
2690  const auto index = fe_indices[range.first];
2691 
2692  for (unsigned int i = range.first; i < range.second; ++i)
2693  AssertDimension(index, fe_indices[i]);
2694 
2695  return index;
2696 }
2697 
2698 
2699 
2700 template <int dim, typename Number, typename VectorizedArrayType>
2701 unsigned int
2703  const std::pair<unsigned int, unsigned int> range,
2704  const bool is_interior_face) const
2705 {
2706  const auto &fe_indices = dof_info[0].cell_active_fe_index;
2707 
2708  if (fe_indices.empty() == true)
2709  return 0;
2710 
2711  if (is_interior_face)
2712  {
2713  const unsigned int index =
2714  fe_indices[face_info.faces[range.first].cells_interior[0] /
2715  VectorizedArrayType::size()];
2716 
2717  for (unsigned int i = range.first; i < range.second; ++i)
2718  AssertDimension(index,
2719  fe_indices[face_info.faces[i].cells_interior[0] /
2720  VectorizedArrayType::size()]);
2721 
2722  return index;
2723  }
2724  else
2725  {
2726  const unsigned int index =
2727  fe_indices[face_info.faces[range.first].cells_exterior[0] /
2728  VectorizedArrayType::size()];
2729 
2730  for (unsigned int i = range.first; i < range.second; ++i)
2731  AssertDimension(index,
2732  fe_indices[face_info.faces[i].cells_exterior[0] /
2733  VectorizedArrayType::size()]);
2734 
2735  return index;
2736  }
2737 }
2738 
2739 
2740 
2741 template <int dim, typename Number, typename VectorizedArrayType>
2742 inline unsigned int
2744  const unsigned int cell_batch_index) const
2745 {
2746  Assert(!dof_info.empty(), ExcNotInitialized());
2747  AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2748  const std::vector<unsigned char> &n_lanes_filled =
2749  dof_info[0].n_vectorization_lanes_filled
2751  AssertIndexRange(cell_batch_index, n_lanes_filled.size());
2752 
2753  return n_lanes_filled[cell_batch_index];
2754 }
2755 
2756 
2757 
2758 template <int dim, typename Number, typename VectorizedArrayType>
2759 inline unsigned int
2761  const unsigned int face_batch_index) const
2762 {
2763  AssertIndexRange(face_batch_index, face_info.faces.size());
2764  Assert(!dof_info.empty(), ExcNotInitialized());
2765  const std::vector<unsigned char> &n_lanes_filled =
2766  dof_info[0].n_vectorization_lanes_filled
2768  AssertIndexRange(face_batch_index, n_lanes_filled.size());
2769  return n_lanes_filled[face_batch_index];
2770 }
2771 
2772 
2773 
2774 template <int dim, typename Number, typename VectorizedArrayType>
2775 inline unsigned int
2777  const unsigned int dof_handler_index,
2778  const unsigned int active_fe_index) const
2779 {
2780  return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2781 }
2782 
2783 
2784 
2785 template <int dim, typename Number, typename VectorizedArrayType>
2786 inline unsigned int
2788  const unsigned int quad_index,
2789  const unsigned int active_fe_index) const
2790 {
2791  AssertIndexRange(quad_index, mapping_info.cell_data.size());
2792  return mapping_info.cell_data[quad_index]
2793  .descriptor[active_fe_index]
2794  .n_q_points;
2795 }
2796 
2797 
2798 
2799 template <int dim, typename Number, typename VectorizedArrayType>
2800 inline unsigned int
2802  const unsigned int dof_handler_index,
2803  const unsigned int active_fe_index) const
2804 {
2805  return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2806 }
2807 
2808 
2809 
2810 template <int dim, typename Number, typename VectorizedArrayType>
2811 inline unsigned int
2813  const unsigned int quad_index,
2814  const unsigned int active_fe_index) const
2815 {
2816  AssertIndexRange(quad_index, mapping_info.face_data.size());
2817  return mapping_info.face_data[quad_index]
2818  .descriptor[active_fe_index]
2819  .n_q_points;
2820 }
2821 
2822 
2823 
2824 template <int dim, typename Number, typename VectorizedArrayType>
2825 inline const IndexSet &
2827  const unsigned int dof_handler_index) const
2828 {
2829  return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2830 }
2831 
2832 
2833 
2834 template <int dim, typename Number, typename VectorizedArrayType>
2835 inline const IndexSet &
2837  const unsigned int dof_handler_index) const
2838 {
2839  return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2840 }
2841 
2842 
2843 
2844 template <int dim, typename Number, typename VectorizedArrayType>
2847  const unsigned int dof_handler_index,
2848  const unsigned int index_quad,
2849  const unsigned int index_fe,
2850  const unsigned int active_fe_index,
2851  const unsigned int active_quad_index) const
2852 {
2853  AssertIndexRange(dof_handler_index, dof_info.size());
2854  const unsigned int ind =
2855  dof_info[dof_handler_index].global_base_element_offset + index_fe;
2856  AssertIndexRange(ind, shape_info.size(0));
2857  AssertIndexRange(index_quad, shape_info.size(1));
2858  AssertIndexRange(active_fe_index, shape_info.size(2));
2859  AssertIndexRange(active_quad_index, shape_info.size(3));
2860  return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2861 }
2862 
2863 
2864 
2865 template <int dim, typename Number, typename VectorizedArrayType>
2867  VectorizedArrayType::size()> &
2869  const unsigned int face_batch_index) const
2870 {
2871  AssertIndexRange(face_batch_index, face_info.faces.size());
2872  return face_info.faces[face_batch_index];
2873 }
2874 
2875 
2876 
2877 template <int dim, typename Number, typename VectorizedArrayType>
2878 inline const Table<3, unsigned int> &
2880  const
2881 {
2883 }
2884 
2885 
2886 
2887 template <int dim, typename Number, typename VectorizedArrayType>
2888 inline const Quadrature<dim> &
2890  const unsigned int quad_index,
2891  const unsigned int active_fe_index) const
2892 {
2893  AssertIndexRange(quad_index, mapping_info.cell_data.size());
2894  return mapping_info.cell_data[quad_index]
2895  .descriptor[active_fe_index]
2896  .quadrature;
2897 }
2898 
2899 
2900 
2901 template <int dim, typename Number, typename VectorizedArrayType>
2902 inline const Quadrature<dim - 1> &
2904  const unsigned int quad_index,
2905  const unsigned int active_fe_index) const
2906 {
2907  AssertIndexRange(quad_index, mapping_info.face_data.size());
2908  return mapping_info.face_data[quad_index]
2909  .descriptor[active_fe_index]
2910  .quadrature;
2911 }
2912 
2913 
2914 
2915 template <int dim, typename Number, typename VectorizedArrayType>
2916 inline unsigned int
2918  const std::pair<unsigned int, unsigned int> range) const
2919 {
2920  auto result = get_cell_category(range.first);
2921 
2922  for (unsigned int i = range.first; i < range.second; ++i)
2923  result = std::max(result, get_cell_category(i));
2924 
2925  return result;
2926 }
2927 
2928 
2929 
2930 template <int dim, typename Number, typename VectorizedArrayType>
2931 inline std::pair<unsigned int, unsigned int>
2933  const std::pair<unsigned int, unsigned int> range) const
2934 {
2935  auto result = get_face_category(range.first);
2936 
2937  for (unsigned int i = range.first; i < range.second; ++i)
2938  {
2939  result.first = std::max(result.first, get_face_category(i).first);
2940  result.second = std::max(result.second, get_face_category(i).second);
2941  }
2942 
2943  return result;
2944 }
2945 
2946 
2947 
2948 template <int dim, typename Number, typename VectorizedArrayType>
2949 inline unsigned int
2951  const unsigned int cell_batch_index) const
2952 {
2953  AssertIndexRange(0, dof_info.size());
2954  AssertIndexRange(cell_batch_index, dof_info[0].cell_active_fe_index.size());
2955  if (dof_info[0].cell_active_fe_index.empty())
2956  return 0;
2957  else
2958  return dof_info[0].cell_active_fe_index[cell_batch_index];
2959 }
2960 
2961 
2962 
2963 template <int dim, typename Number, typename VectorizedArrayType>
2964 inline std::pair<unsigned int, unsigned int>
2966  const unsigned int face_batch_index) const
2967 {
2968  AssertIndexRange(face_batch_index, face_info.faces.size());
2969  if (dof_info[0].cell_active_fe_index.empty())
2970  return std::make_pair(0U, 0U);
2971 
2972  std::pair<unsigned int, unsigned int> result = std::make_pair(0U, 0U);
2973  for (unsigned int v = 0;
2974  v < VectorizedArrayType::size() &&
2975  face_info.faces[face_batch_index].cells_interior[v] !=
2977  ++v)
2978  result.first = std::max(
2979  result.first,
2980  dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
2981  .cells_interior[v] /
2982  VectorizedArrayType::size()]);
2983  if (face_info.faces[face_batch_index].cells_exterior[0] !=
2985  for (unsigned int v = 0;
2986  v < VectorizedArrayType::size() &&
2987  face_info.faces[face_batch_index].cells_exterior[v] !=
2989  ++v)
2990  result.second = std::max(
2991  result.second,
2992  dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
2993  .cells_exterior[v] /
2994  VectorizedArrayType::size()]);
2995  else
2996  result.second = numbers::invalid_unsigned_int;
2997  return result;
2998 }
2999 
3000 
3001 
3002 template <int dim, typename Number, typename VectorizedArrayType>
3003 inline bool
3005 {
3006  return indices_are_initialized;
3007 }
3008 
3009 
3010 
3011 template <int dim, typename Number, typename VectorizedArrayType>
3012 inline bool
3014 {
3015  return mapping_is_initialized;
3016 }
3017 
3018 
3019 template <int dim, typename Number, typename VectorizedArrayType>
3020 inline unsigned int
3022 {
3023  return mg_level;
3024 }
3025 
3026 
3027 
3028 template <int dim, typename Number, typename VectorizedArrayType>
3031 {
3032  using list_type =
3033  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3034  list_type &data = scratch_pad.get();
3035  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
3036  if (it->first == false)
3037  {
3038  it->first = true;
3039  return &it->second;
3040  }
3041  data.emplace_front(true, AlignedVector<VectorizedArrayType>());
3042  return &data.front().second;
3043 }
3044 
3045 
3046 
3047 template <int dim, typename Number, typename VectorizedArrayType>
3048 void
3050  const AlignedVector<VectorizedArrayType> *scratch) const
3051 {
3052  using list_type =
3053  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3054  list_type &data = scratch_pad.get();
3055  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
3056  if (&it->second == scratch)
3057  {
3058  Assert(it->first == true, ExcInternalError());
3059  it->first = false;
3060  return;
3061  }
3062  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
3063 }
3064 
3065 
3066 
3067 template <int dim, typename Number, typename VectorizedArrayType>
3071 {
3072  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
3074  it != scratch_pad_non_threadsafe.end();
3075  ++it)
3076  if (it->first == false)
3077  {
3078  it->first = true;
3079  return &it->second;
3080  }
3081  scratch_pad_non_threadsafe.push_front(
3082  std::make_pair(true, AlignedVector<Number>()));
3083  return &scratch_pad_non_threadsafe.front().second;
3084 }
3085 
3086 
3087 
3088 template <int dim, typename Number, typename VectorizedArrayType>
3089 void
3092  const AlignedVector<Number> *scratch) const
3093 {
3094  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
3096  it != scratch_pad_non_threadsafe.end();
3097  ++it)
3098  if (&it->second == scratch)
3099  {
3100  Assert(it->first == true, ExcInternalError());
3101  it->first = false;
3102  return;
3103  }
3104  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
3105 }
3106 
3107 
3108 
3109 // ------------------------------ reinit functions ---------------------------
3110 
3111 namespace internal
3112 {
3113  namespace MatrixFreeImplementation
3114  {
3115  template <int dim, int spacedim>
3116  inline std::vector<IndexSet>
3117  extract_locally_owned_index_sets(
3118  const std::vector<const ::DoFHandler<dim, spacedim> *> &dofh,
3119  const unsigned int level)
3120  {
3121  std::vector<IndexSet> locally_owned_set;
3122  locally_owned_set.reserve(dofh.size());
3123  for (unsigned int j = 0; j < dofh.size(); ++j)
3125  locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
3126  else
3127  locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
3128  return locally_owned_set;
3129  }
3130  } // namespace MatrixFreeImplementation
3131 } // namespace internal
3132 
3133 
3134 
3135 template <int dim, typename Number, typename VectorizedArrayType>
3136 template <typename QuadratureType, typename number2, typename MappingType>
3137 void
3139  const MappingType & mapping,
3140  const DoFHandler<dim> & dof_handler,
3141  const AffineConstraints<number2> &constraints_in,
3142  const QuadratureType & quad,
3144  &additional_data)
3145 {
3146  std::vector<const DoFHandler<dim, dim> *> dof_handlers;
3147  std::vector<const AffineConstraints<number2> *> constraints;
3148 
3149  dof_handlers.push_back(&dof_handler);
3150  constraints.push_back(&constraints_in);
3151 
3152  std::vector<IndexSet> locally_owned_sets =
3153  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3154  dof_handlers, additional_data.mg_level);
3155 
3156  std::vector<hp::QCollection<dim>> quad_hp;
3157  quad_hp.emplace_back(quad);
3158 
3159  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3160  dof_handlers,
3161  constraints,
3162  locally_owned_sets,
3163  quad_hp,
3164  additional_data);
3165 }
3166 
3167 
3168 
3169 template <int dim, typename Number, typename VectorizedArrayType>
3170 template <typename QuadratureType, typename number2, typename MappingType>
3171 void
3173  const MappingType & mapping,
3174  const std::vector<const DoFHandler<dim> *> & dof_handler,
3175  const std::vector<const AffineConstraints<number2> *> &constraint,
3176  const QuadratureType & quad,
3178  &additional_data)
3179 {
3180  std::vector<IndexSet> locally_owned_set =
3181  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3182  dof_handler, additional_data.mg_level);
3183  std::vector<hp::QCollection<dim>> quad_hp;
3184  quad_hp.emplace_back(quad);
3185 
3186  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3187  dof_handler,
3188  constraint,
3189  locally_owned_set,
3190  quad_hp,
3191  additional_data);
3192 }
3193 
3194 
3195 
3196 template <int dim, typename Number, typename VectorizedArrayType>
3197 template <typename QuadratureType, typename number2, typename MappingType>
3198 void
3200  const MappingType & mapping,
3201  const std::vector<const DoFHandler<dim> *> & dof_handler,
3202  const std::vector<const AffineConstraints<number2> *> &constraint,
3203  const std::vector<QuadratureType> & quad,
3205  &additional_data)
3206 {
3207  std::vector<IndexSet> locally_owned_set =
3208  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3209  dof_handler, additional_data.mg_level);
3210  std::vector<hp::QCollection<dim>> quad_hp;
3211  for (unsigned int q = 0; q < quad.size(); ++q)
3212  quad_hp.emplace_back(quad[q]);
3213 
3214  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3215  dof_handler,
3216  constraint,
3217  locally_owned_set,
3218  quad_hp,
3219  additional_data);
3220 }
3221 
3222 
3223 
3224 // ------------------------------ implementation of loops --------------------
3225 
3226 // internal helper functions that define how to call MPI data exchange
3227 // functions: for generic vectors, do nothing at all. For distributed vectors,
3228 // call update_ghost_values_start function and so on. If we have collections
3229 // of vectors, just do the individual functions of the components. In order to
3230 // keep ghost values consistent (whether we are in read or write mode), we
3231 // also reset the values at the end. the whole situation is a bit complicated
3232 // by the fact that we need to treat block vectors differently, which use some
3233 // additional helper functions to select the blocks and template magic.
3234 namespace internal
3235 {
3239  template <int dim, typename Number, typename VectorizedArrayType>
3240  struct VectorDataExchange
3241  {
3242  // A shift for the MPI messages to reduce the risk for accidental
3243  // interaction with other open communications that a user program might
3244  // set up (parallel vectors support unfinished communication). We let
3245  // the other vectors use the first 20 assigned numbers and start the
3246  // matrix-free communication.
3247  static constexpr unsigned int channel_shift = 20;
3248 
3249 
3250 
3255  VectorDataExchange(
3256  const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3257  const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3258  DataAccessOnFaces vector_face_access,
3259  const unsigned int n_components)
3260  : matrix_free(matrix_free)
3261  , vector_face_access(
3262  matrix_free.get_task_info().face_partition_data.empty() ?
3263  ::MatrixFree<dim, Number, VectorizedArrayType>::
3264  DataAccessOnFaces::unspecified :
3265  vector_face_access)
3266  , ghosts_were_set(false)
3267 # ifdef DEAL_II_WITH_MPI
3268  , tmp_data(n_components)
3269  , requests(n_components)
3270 # endif
3271  {
3272  (void)n_components;
3273  if (this->vector_face_access !=
3275  DataAccessOnFaces::unspecified)
3276  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3278  matrix_free.get_dof_info(c).vector_exchanger_face_variants.size(),
3279  5);
3280  }
3281 
3282 
3283 
3287  ~VectorDataExchange() // NOLINT
3288  {
3289 # ifdef DEAL_II_WITH_MPI
3290  for (unsigned int i = 0; i < tmp_data.size(); ++i)
3291  if (tmp_data[i] != nullptr)
3292  matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3293 # endif
3294  }
3295 
3296 
3297 
3302  template <typename VectorType>
3303  unsigned int
3304  find_vector_in_mf(const VectorType &vec,
3305  const bool check_global_compatibility = true) const
3306  {
3307  // case 1: vector was set up with MatrixFree::initialize_dof_vector()
3308  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3309  if (vec.get_partitioner().get() ==
3310  matrix_free.get_dof_info(c).vector_partitioner.get())
3311  return c;
3312 
3313  // case 2: user provided own partitioner (compatibility mode)
3314  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3315  if (check_global_compatibility ?
3316  vec.get_partitioner()->is_globally_compatible(
3317  *matrix_free.get_dof_info(c).vector_partitioner) :
3318  vec.get_partitioner()->is_compatible(
3319  *matrix_free.get_dof_info(c).vector_partitioner))
3320  return c;
3321 
3322  Assert(false,
3323  ExcNotImplemented("Could not find partitioner that fits vector"));
3324 
3326  }
3327 
3328 
3329 
3335  get_partitioner(const unsigned int mf_component) const
3336  {
3337  AssertDimension(matrix_free.get_dof_info(mf_component)
3338  .vector_exchanger_face_variants.size(),
3339  5);
3340  if (vector_face_access ==
3343  return *matrix_free.get_dof_info(mf_component)
3344  .vector_exchanger_face_variants[0];
3345  else if (vector_face_access ==
3348  return *matrix_free.get_dof_info(mf_component)
3349  .vector_exchanger_face_variants[1];
3350  else if (vector_face_access ==
3353  return *matrix_free.get_dof_info(mf_component)
3354  .vector_exchanger_face_variants[2];
3355  else if (vector_face_access ==
3357  DataAccessOnFaces::values_all_faces)
3358  return *matrix_free.get_dof_info(mf_component)
3359  .vector_exchanger_face_variants[3];
3360  else if (vector_face_access ==
3362  DataAccessOnFaces::gradients_all_faces)
3363  return *matrix_free.get_dof_info(mf_component)
3364  .vector_exchanger_face_variants[4];
3365  else
3366  return *matrix_free.get_dof_info(mf_component).vector_exchanger.get();
3367  }
3368 
3369 
3370 
3374  template <typename VectorType,
3375  std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3376  * = nullptr>
3377  void
3378  update_ghost_values_start(const unsigned int /*component_in_block_vector*/,
3379  const VectorType & /*vec*/)
3380  {}
3381 
3382 
3387  template <typename VectorType,
3388  std::enable_if_t<!has_update_ghost_values_start<VectorType> &&
3389  !is_not_parallel_vector<VectorType>,
3390  VectorType> * = nullptr>
3391  void
3392  update_ghost_values_start(const unsigned int component_in_block_vector,
3393  const VectorType & vec)
3394  {
3395  (void)component_in_block_vector;
3396  bool ghosts_set = vec.has_ghost_elements();
3397 
3398  Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3399  ghosts_set == false,
3400  ExcNotImplemented());
3401 
3402  if (ghosts_set)
3403  ghosts_were_set = true;
3404 
3405  vec.update_ghost_values();
3406  }
3407 
3408 
3409 
3415  template <typename VectorType,
3416  std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3417  !has_exchange_on_subset<VectorType>,
3418  VectorType> * = nullptr>
3419  void
3420  update_ghost_values_start(const unsigned int component_in_block_vector,
3421  const VectorType & vec)
3422  {
3423  (void)component_in_block_vector;
3424  bool ghosts_set = vec.has_ghost_elements();
3425 
3426  Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3427  ghosts_set == false,
3428  ExcNotImplemented());
3429 
3430  if (ghosts_set)
3431  ghosts_were_set = true;
3432 
3433  vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3434  }
3435 
3436 
3437 
3444  template <typename VectorType,
3445  std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3446  has_exchange_on_subset<VectorType>,
3447  VectorType> * = nullptr>
3448  void
3449  update_ghost_values_start(const unsigned int component_in_block_vector,
3450  const VectorType & vec)
3451  {
3452  static_assert(
3453  std::is_same<Number, typename VectorType::value_type>::value,
3454  "Type mismatch between VectorType and VectorDataExchange");
3455  (void)component_in_block_vector;
3456  bool ghosts_set = vec.has_ghost_elements();
3457 
3458  Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3459  ghosts_set == false,
3460  ExcNotImplemented());
3461 
3462  if (ghosts_set)
3463  ghosts_were_set = true;
3464 
3465  if (vec.size() != 0)
3466  {
3467 # ifdef DEAL_II_WITH_MPI
3468  const unsigned int mf_component = find_vector_in_mf(vec);
3469 
3470  const auto &part = get_partitioner(mf_component);
3471 
3472  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3473  part.n_import_sm_procs() == 0)
3474  return;
3475 
3476  tmp_data[component_in_block_vector] =
3477  matrix_free.acquire_scratch_data_non_threadsafe();
3478  tmp_data[component_in_block_vector]->resize_fast(
3479  part.n_import_indices());
3480  AssertDimension(requests.size(), tmp_data.size());
3481 
3482  part.export_to_ghosted_array_start(
3483  component_in_block_vector * 2 + channel_shift,
3484  ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3485  vec.shared_vector_data(),
3486  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3487  part.locally_owned_size(),
3488  matrix_free.get_dof_info(mf_component)
3489  .vector_partitioner->n_ghost_indices()),
3490  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3491  part.n_import_indices()),
3492  this->requests[component_in_block_vector]);
3493 # endif
3494  }
3495  }
3496 
3497 
3498 
3503  template <typename VectorType,
3504  std::enable_if_t<!has_update_ghost_values_start<VectorType>,
3505  VectorType> * = nullptr>
3506  void
3507  update_ghost_values_finish(const unsigned int /*component_in_block_vector*/,
3508  const VectorType & /*vec*/)
3509  {}
3510 
3511 
3512 
3518  template <typename VectorType,
3519  std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3520  !has_exchange_on_subset<VectorType>,
3521  VectorType> * = nullptr>
3522  void
3523  update_ghost_values_finish(const unsigned int component_in_block_vector,
3524  const VectorType & vec)
3525  {
3526  (void)component_in_block_vector;
3527  vec.update_ghost_values_finish();
3528  }
3529 
3530 
3531 
3538  template <typename VectorType,
3539  std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3540  has_exchange_on_subset<VectorType>,
3541  VectorType> * = nullptr>
3542  void
3543  update_ghost_values_finish(const unsigned int component_in_block_vector,
3544  const VectorType & vec)
3545  {
3546  static_assert(
3547  std::is_same<Number, typename VectorType::value_type>::value,
3548  "Type mismatch between VectorType and VectorDataExchange");
3549  (void)component_in_block_vector;
3550 
3551  if (vec.size() != 0)
3552  {
3553 # ifdef DEAL_II_WITH_MPI
3554  AssertIndexRange(component_in_block_vector, tmp_data.size());
3555  AssertDimension(requests.size(), tmp_data.size());
3556 
3557  const unsigned int mf_component = find_vector_in_mf(vec);
3558 
3559  const auto &part = get_partitioner(mf_component);
3560 
3561  if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3562  part.n_import_sm_procs() != 0)
3563  {
3564  part.export_to_ghosted_array_finish(
3565  ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3566  vec.shared_vector_data(),
3567  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3568  part.locally_owned_size(),
3569  matrix_free.get_dof_info(mf_component)
3570  .vector_partitioner->n_ghost_indices()),
3571  this->requests[component_in_block_vector]);
3572 
3573  matrix_free.release_scratch_data_non_threadsafe(
3574  tmp_data[component_in_block_vector]);
3575  tmp_data[component_in_block_vector] = nullptr;
3576  }
3577 # endif
3578  }
3579  // let vector know that ghosts are being updated and we can read from
3580  // them
3581  vec.set_ghost_state(true);
3582  }
3583 
3584 
3585 
3589  template <typename VectorType,
3590  std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3591  * = nullptr>
3592  void
3593  compress_start(const unsigned int /*component_in_block_vector*/,
3594  VectorType & /*vec*/)
3595  {}
3596 
3597 
3598 
3603  template <typename VectorType,
3604  std::enable_if_t<!has_compress_start<VectorType> &&
3605  !is_not_parallel_vector<VectorType>,
3606  VectorType> * = nullptr>
3607  void
3608  compress_start(const unsigned int component_in_block_vector,
3609  VectorType & vec)
3610  {
3611  (void)component_in_block_vector;
3612  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3613  vec.compress(::VectorOperation::add);
3614  }
3615 
3616 
3617 
3623  template <typename VectorType,
3624  std::enable_if_t<has_compress_start<VectorType> &&
3625  !has_exchange_on_subset<VectorType>,
3626  VectorType> * = nullptr>
3627  void
3628  compress_start(const unsigned int component_in_block_vector,
3629  VectorType & vec)
3630  {
3631  (void)component_in_block_vector;
3632  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3633  vec.compress_start(component_in_block_vector + channel_shift);
3634  }
3635 
3636 
3637 
3644  template <typename VectorType,
3645  std::enable_if_t<has_compress_start<VectorType> &&
3646  has_exchange_on_subset<VectorType>,
3647  VectorType> * = nullptr>
3648  void
3649  compress_start(const unsigned int component_in_block_vector,
3650  VectorType & vec)
3651  {
3652  static_assert(
3653  std::is_same<Number, typename VectorType::value_type>::value,
3654  "Type mismatch between VectorType and VectorDataExchange");
3655  (void)component_in_block_vector;
3656  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3657 
3658  if (vec.size() != 0)
3659  {
3660 # ifdef DEAL_II_WITH_MPI
3661  const unsigned int mf_component = find_vector_in_mf(vec);
3662 
3663  const auto &part = get_partitioner(mf_component);
3664 
3665  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3666  part.n_import_sm_procs() == 0)
3667  return;
3668 
3669  tmp_data[component_in_block_vector] =
3670  matrix_free.acquire_scratch_data_non_threadsafe();
3671  tmp_data[component_in_block_vector]->resize_fast(
3672  part.n_import_indices());
3673  AssertDimension(requests.size(), tmp_data.size());
3674 
3675  part.import_from_ghosted_array_start(
3677  component_in_block_vector * 2 + channel_shift,
3678  ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3679  vec.shared_vector_data(),
3680  ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3681  matrix_free.get_dof_info(mf_component)
3682  .vector_partitioner->n_ghost_indices()),
3683  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3684  part.n_import_indices()),
3685  this->requests[component_in_block_vector]);
3686 # endif
3687  }
3688  }
3689 
3690 
3691 
3696  template <
3697  typename VectorType,
3698  std::enable_if_t<!has_compress_start<VectorType>, VectorType> * = nullptr>
3699  void
3700  compress_finish(const unsigned int /*component_in_block_vector*/,
3701  VectorType & /*vec*/)
3702  {}
3703 
3704 
3705 
3711  template <typename VectorType,
3712  std::enable_if_t<has_compress_start<VectorType> &&
3713  !has_exchange_on_subset<VectorType>,
3714  VectorType> * = nullptr>
3715  void
3716  compress_finish(const unsigned int component_in_block_vector,
3717  VectorType & vec)
3718  {
3719  (void)component_in_block_vector;
3720  vec.compress_finish(::VectorOperation::add);
3721  }
3722 
3723 
3724 
3731  template <typename VectorType,
3732  std::enable_if_t<has_compress_start<VectorType> &&
3733  has_exchange_on_subset<VectorType>,
3734  VectorType> * = nullptr>
3735  void
3736  compress_finish(const unsigned int component_in_block_vector,
3737  VectorType & vec)
3738  {
3739  static_assert(
3740  std::is_same<Number, typename VectorType::value_type>::value,
3741  "Type mismatch between VectorType and VectorDataExchange");
3742  (void)component_in_block_vector;
3743  if (vec.size() != 0)
3744  {
3745 # ifdef DEAL_II_WITH_MPI
3746  AssertIndexRange(component_in_block_vector, tmp_data.size());
3747  AssertDimension(requests.size(), tmp_data.size());
3748 
3749  const unsigned int mf_component = find_vector_in_mf(vec);
3750 
3751  const auto &part = get_partitioner(mf_component);
3752 
3753  if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3754  part.n_import_sm_procs() != 0)
3755  {
3756  part.import_from_ghosted_array_finish(
3758  ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3759  vec.shared_vector_data(),
3760  ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3761  matrix_free.get_dof_info(mf_component)
3762  .vector_partitioner->n_ghost_indices()),
3764  tmp_data[component_in_block_vector]->begin(),
3765  part.n_import_indices()),
3766  this->requests[component_in_block_vector]);
3767 
3768  matrix_free.release_scratch_data_non_threadsafe(
3769  tmp_data[component_in_block_vector]);
3770  tmp_data[component_in_block_vector] = nullptr;
3771  }
3772 
3774  {
3775  const int ierr =
3776  MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3777  AssertThrowMPI(ierr);
3778  }
3779 # endif
3780  }
3781  }
3782 
3783 
3784 
3788  template <typename VectorType,
3789  std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3790  * = nullptr>
3791  void
3792  reset_ghost_values(const VectorType & /*vec*/) const
3793  {}
3794 
3795 
3796 
3801  template <typename VectorType,
3802  std::enable_if_t<!has_exchange_on_subset<VectorType> &&
3803  !is_not_parallel_vector<VectorType>,
3804  VectorType> * = nullptr>
3805  void
3806  reset_ghost_values(const VectorType &vec) const
3807  {
3808  if (ghosts_were_set == true)
3809  return;
3810 
3811  vec.zero_out_ghost_values();
3812  }
3813 
3814 
3815 
3821  template <typename VectorType,
3822  std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3823  * = nullptr>
3824  void
3825  reset_ghost_values(const VectorType &vec) const
3826  {
3827  static_assert(
3828  std::is_same<Number, typename VectorType::value_type>::value,
3829  "Type mismatch between VectorType and VectorDataExchange");
3830  if (ghosts_were_set == true)
3831  return;
3832 
3833  if (vec.size() != 0)
3834  {
3835 # ifdef DEAL_II_WITH_MPI
3836  AssertDimension(requests.size(), tmp_data.size());
3837 
3838  const unsigned int mf_component = find_vector_in_mf(vec);
3839 
3840  const auto &part = get_partitioner(mf_component);
3841 
3842  if (part.n_ghost_indices() > 0)
3843  {
3844  part.reset_ghost_values(ArrayView<Number>(
3846  .begin() +
3847  part.locally_owned_size(),
3848  matrix_free.get_dof_info(mf_component)
3849  .vector_partitioner->n_ghost_indices()));
3850  }
3851 
3852 # endif
3853  }
3854  // let vector know that it's not ghosted anymore
3855  vec.set_ghost_state(false);
3856  }
3857 
3858 
3859 
3865  template <typename VectorType,
3866  std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3867  * = nullptr>
3868  void
3869  zero_vector_region(const unsigned int range_index, VectorType &vec) const
3870  {
3871  static_assert(
3872  std::is_same<Number, typename VectorType::value_type>::value,
3873  "Type mismatch between VectorType and VectorDataExchange");
3874  if (range_index == numbers::invalid_unsigned_int)
3875  vec = Number();
3876  else
3877  {
3878  const unsigned int mf_component = find_vector_in_mf(vec, false);
3879  const internal::MatrixFreeFunctions::DoFInfo &dof_info =
3880  matrix_free.get_dof_info(mf_component);
3881  Assert(dof_info.vector_zero_range_list_index.empty() == false,
3882  ExcNotInitialized());
3883 
3884  Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3885  ExcInternalError());
3886  AssertIndexRange(range_index,
3887  dof_info.vector_zero_range_list_index.size() - 1);
3888  for (unsigned int id =
3889  dof_info.vector_zero_range_list_index[range_index];
3890  id != dof_info.vector_zero_range_list_index[range_index + 1];
3891  ++id)
3892  std::memset(vec.begin() + dof_info.vector_zero_range_list[id].first,
3893  0,
3894  (dof_info.vector_zero_range_list[id].second -
3895  dof_info.vector_zero_range_list[id].first) *
3896  sizeof(Number));
3897  }
3898  }
3899 
3900 
3901 
3907  template <typename VectorType,
3908  std::enable_if_t<!has_exchange_on_subset<VectorType>, VectorType>
3909  * = nullptr,
3910  typename VectorType::value_type * = nullptr>
3911  void
3912  zero_vector_region(const unsigned int range_index, VectorType &vec) const
3913  {
3914  if (range_index == numbers::invalid_unsigned_int || range_index == 0)
3915  vec = typename VectorType::value_type();
3916  }
3917 
3918 
3919 
3924  void
3925  zero_vector_region(const unsigned int, ...) const
3926  {
3927  Assert(false,
3928  ExcNotImplemented("Zeroing is only implemented for vector types "
3929  "which provide VectorType::value_type"));
3930  }
3931 
3932 
3933 
3934  const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
3935  const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3936  DataAccessOnFaces vector_face_access;
3937  bool ghosts_were_set;
3938 # ifdef DEAL_II_WITH_MPI
3939  std::vector<AlignedVector<Number> *> tmp_data;
3940  std::vector<std::vector<MPI_Request>> requests;
3941 # endif
3942  }; // VectorDataExchange
3943 
3944  template <typename VectorStruct>
3945  unsigned int
3946  n_components(const VectorStruct &vec);
3947 
3948  template <typename VectorStruct>
3949  unsigned int
3950  n_components_block(const VectorStruct &vec,
3951  std::integral_constant<bool, true>)
3952  {
3953  unsigned int components = 0;
3954  for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
3955  components += n_components(vec.block(bl));
3956  return components;
3957  }
3958 
3959  template <typename VectorStruct>
3960  unsigned int
3961  n_components_block(const VectorStruct &, std::integral_constant<bool, false>)
3962  {
3963  return 1;
3964  }
3965 
3966  template <typename VectorStruct>
3967  unsigned int
3968  n_components(const VectorStruct &vec)
3969  {
3970  return n_components_block(
3971  vec, std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3972  }
3973 
3974  template <typename VectorStruct>
3975  inline unsigned int
3976  n_components(const std::vector<VectorStruct> &vec)
3977  {
3978  unsigned int components = 0;
3979  for (unsigned int comp = 0; comp < vec.size(); ++comp)
3980  components += n_components_block(
3981  vec[comp],
3982  std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3983  return components;
3984  }
3985 
3986  template <typename VectorStruct>
3987  inline unsigned int
3988  n_components(const std::vector<VectorStruct *> &vec)
3989  {
3990  unsigned int components = 0;
3991  for (unsigned int comp = 0; comp < vec.size(); ++comp)
3992  components += n_components_block(
3993  *vec[comp],
3994  std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3995  return components;
3996  }
3997 
3998 
3999 
4000  // A helper function to identify block vectors with many components where we
4001  // should not try to overlap computations and communication because there
4002  // would be too many outstanding communication requests.
4003 
4004  // default value for vectors that do not have communication_block_size
4005  template <typename VectorStruct,
4006  std::enable_if_t<!has_communication_block_size<VectorStruct>,
4007  VectorStruct> * = nullptr>
4008  constexpr unsigned int
4009  get_communication_block_size(const VectorStruct &)
4010  {
4012  }
4013 
4014 
4015 
4016  template <typename VectorStruct,
4017  std::enable_if_t<has_communication_block_size<VectorStruct>,
4018  VectorStruct> * = nullptr>
4019  constexpr unsigned int
4020  get_communication_block_size(const VectorStruct &)
4021  {
4022  return VectorStruct::communication_block_size;
4023  }
4024 
4025 
4026 
4027  template <typename VectorType,
4028  std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType> * =
4029  nullptr>
4030  bool
4031  has_ghost_elements(const VectorType &vec)
4032  {
4033  (void)vec;
4034  return false;
4035  }
4036 
4037 
4038 
4039  template <typename VectorType,
4040  std::enable_if_t<!is_not_parallel_vector<VectorType>, VectorType>
4041  * = nullptr>
4042  bool
4043  has_ghost_elements(const VectorType &vec)
4044  {
4045  return vec.has_ghost_elements();
4046  }
4047 
4048 
4049 
4050  // --------------------------------------------------------------------------
4051  // below we have wrappers to distinguish between block and non-block vectors.
4052  // --------------------------------------------------------------------------
4053 
4054  //
4055  // update_ghost_values_start
4056  //
4057 
4058  // update_ghost_values for block vectors
4059  template <int dim,
4060  typename VectorStruct,
4061  typename Number,
4062  typename VectorizedArrayType,
4063  std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4064  * = nullptr>
4065  void
4066  update_ghost_values_start(
4067  const VectorStruct & vec,
4068  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4069  const unsigned int channel = 0)
4070  {
4071  if (get_communication_block_size(vec) < vec.n_blocks())
4072  {
4073  // don't forget to set ghosts_were_set, that otherwise happens
4074  // inside VectorDataExchange::update_ghost_values_start()
4075  exchanger.ghosts_were_set = vec.has_ghost_elements();
4076  vec.update_ghost_values();
4077  }
4078  else
4079  {
4080  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4081  update_ghost_values_start(vec.block(i), exchanger, channel + i);
4082  }
4083  }
4084 
4085 
4086 
4087  // update_ghost_values for non-block vectors
4088  template <int dim,
4089  typename VectorStruct,
4090  typename Number,
4091  typename VectorizedArrayType,
4092  std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4093  * = nullptr>
4094  void
4095  update_ghost_values_start(
4096  const VectorStruct & vec,
4097  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4098  const unsigned int channel = 0)
4099  {
4100  exchanger.update_ghost_values_start(channel, vec);
4101  }
4102 
4103 
4104 
4105  // update_ghost_values_start() for vector of vectors
4106  template <int dim,
4107  typename VectorStruct,
4108  typename Number,
4109  typename VectorizedArrayType>
4110  inline void
4111  update_ghost_values_start(
4112  const std::vector<VectorStruct> & vec,
4113  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4114  {
4115  unsigned int component_index = 0;
4116  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4117  {
4118  update_ghost_values_start(vec[comp], exchanger, component_index);
4119  component_index += n_components(vec[comp]);
4120  }
4121  }
4122 
4123 
4124 
4125  // update_ghost_values_start() for vector of pointers to vectors
4126  template <int dim,
4127  typename VectorStruct,
4128  typename Number,
4129  typename VectorizedArrayType>
4130  inline void
4131  update_ghost_values_start(
4132  const std::vector<VectorStruct *> & vec,
4133  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4134  {
4135  unsigned int component_index = 0;
4136  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4137  {
4138  update_ghost_values_start(*vec[comp], exchanger, component_index);
4139  component_index += n_components(*vec[comp]);
4140  }
4141  }
4142 
4143 
4144 
4145  //
4146  // update_ghost_values_finish
4147  //
4148 
4149  // for block vectors
4150  template <int dim,
4151  typename VectorStruct,
4152  typename Number,
4153  typename VectorizedArrayType,
4154  std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4155  * = nullptr>
4156  void
4157  update_ghost_values_finish(
4158  const VectorStruct & vec,
4159  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4160  const unsigned int channel = 0)
4161  {
4162  if (get_communication_block_size(vec) < vec.n_blocks())
4163  {
4164  // do nothing, everything has already been completed in the _start()
4165  // call
4166  }
4167  else
4168  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4169  update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4170  }
4171 
4172 
4173 
4174  // for non-block vectors
4175  template <int dim,
4176  typename VectorStruct,
4177  typename Number,
4178  typename VectorizedArrayType,
4179  std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4180  * = nullptr>
4181  void
4182  update_ghost_values_finish(
4183  const VectorStruct & vec,
4184  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4185  const unsigned int channel = 0)
4186  {
4187  exchanger.update_ghost_values_finish(channel, vec);
4188  }
4189 
4190 
4191 
4192  // for vector of vectors
4193  template <int dim,
4194  typename VectorStruct,
4195  typename Number,
4196  typename VectorizedArrayType>
4197  inline void
4198  update_ghost_values_finish(
4199  const std::vector<VectorStruct> & vec,
4200  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4201  {
4202  unsigned int component_index = 0;
4203  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4204  {
4205  update_ghost_values_finish(vec[comp], exchanger, component_index);
4206  component_index += n_components(vec[comp]);
4207  }
4208  }
4209 
4210 
4211 
4212  // for vector of pointers to vectors
4213  template <int dim,
4214  typename VectorStruct,
4215  typename Number,
4216  typename VectorizedArrayType>
4217  inline void
4218  update_ghost_values_finish(
4219  const std::vector<VectorStruct *> & vec,
4220  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4221  {
4222  unsigned int component_index = 0;
4223  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4224  {
4225  update_ghost_values_finish(*vec[comp], exchanger, component_index);
4226  component_index += n_components(*vec[comp]);
4227  }
4228  }
4229 
4230 
4231 
4232  //
4233  // compress_start
4234  //
4235 
4236  // for block vectors
4237  template <int dim,
4238  typename VectorStruct,
4239  typename Number,
4240  typename VectorizedArrayType,
4241  std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4242  * = nullptr>
4243  inline void
4244  compress_start(
4245  VectorStruct & vec,
4246  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4247  const unsigned int channel = 0)
4248  {
4249  if (get_communication_block_size(vec) < vec.n_blocks())
4250  vec.compress(::VectorOperation::add);
4251  else
4252  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4253  compress_start(vec.block(i), exchanger, channel + i);
4254  }
4255 
4256 
4257 
4258  // for non-block vectors
4259  template <int dim,
4260  typename VectorStruct,
4261  typename Number,
4262  typename VectorizedArrayType,
4263  std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4264  * = nullptr>
4265  inline void
4266  compress_start(
4267  VectorStruct & vec,
4268  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4269  const unsigned int channel = 0)
4270  {
4271  exchanger.compress_start(channel, vec);
4272  }
4273 
4274 
4275 
4276  // for std::vector of vectors
4277  template <int dim,
4278  typename VectorStruct,
4279  typename Number,
4280  typename VectorizedArrayType>
4281  inline void
4282  compress_start(
4283  std::vector<VectorStruct> & vec,
4284  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4285  {
4286  unsigned int component_index = 0;
4287  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4288  {
4289  compress_start(vec[comp], exchanger, component_index);
4290  component_index += n_components(vec[comp]);
4291  }
4292  }
4293 
4294 
4295 
4296  // for std::vector of pointer to vectors
4297  template <int dim,
4298  typename VectorStruct,
4299  typename Number,
4300  typename VectorizedArrayType>
4301  inline void
4302  compress_start(
4303  std::vector<VectorStruct *> & vec,
4304  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4305  {
4306  unsigned int component_index = 0;
4307  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4308  {
4309  compress_start(*vec[comp], exchanger, component_index);
4310  component_index += n_components(*vec[comp]);
4311  }
4312  }
4313 
4314 
4315 
4316  //
4317  // compress_finish
4318  //
4319 
4320  // for block vectors
4321  template <int dim,
4322  typename VectorStruct,
4323  typename Number,
4324  typename VectorizedArrayType,
4325  std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4326  * = nullptr>
4327  inline void
4328  compress_finish(
4329  VectorStruct & vec,
4330  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4331  const unsigned int channel = 0)
4332  {
4333  if (get_communication_block_size(vec) < vec.n_blocks())
4334  {
4335  // do nothing, everything has already been completed in the _start()
4336  // call
4337  }
4338  else
4339  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4340  compress_finish(vec.block(i), exchanger, channel + i);
4341  }
4342 
4343 
4344 
4345  // for non-block vectors
4346  template <int dim,
4347  typename VectorStruct,
4348  typename Number,
4349  typename VectorizedArrayType,
4350  std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4351  * = nullptr>
4352  inline void
4353  compress_finish(
4354  VectorStruct & vec,
4355  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4356  const unsigned int channel = 0)
4357  {
4358  exchanger.compress_finish(channel, vec);
4359  }
4360 
4361 
4362 
4363  // for std::vector of vectors
4364  template <int dim,
4365  typename VectorStruct,
4366  typename Number,
4367  typename VectorizedArrayType>
4368  inline void
4369  compress_finish(
4370  std::vector<VectorStruct> & vec,
4371  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4372  {
4373  unsigned int component_index = 0;
4374  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4375  {
4376  compress_finish(vec[comp], exchanger, component_index);
4377  component_index += n_components(vec[comp]);
4378  }
4379  }
4380 
4381 
4382 
4383  // for std::vector of pointer to vectors
4384  template <int dim,
4385  typename VectorStruct,
4386  typename Number,
4387  typename VectorizedArrayType>
4388  inline void
4389  compress_finish(
4390  std::vector<VectorStruct *> & vec,
4391  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4392  {
4393  unsigned int component_index = 0;
4394  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4395  {
4396  compress_finish(*vec[comp], exchanger, component_index);
4397  component_index += n_components(*vec[comp]);
4398  }
4399  }
4400 
4401 
4402 
4403  //
4404  // reset_ghost_values:
4405  //
4406  // if the input vector did not have ghosts imported, clear them here again
4407  // in order to avoid subsequent operations e.g. in linear solvers to work
4408  // with ghosts all the time
4409  //
4410 
4411  // for block vectors
4412  template <int dim,
4413  typename VectorStruct,
4414  typename Number,
4415  typename VectorizedArrayType,
4416  std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4417  * = nullptr>
4418  inline void
4419  reset_ghost_values(
4420  const VectorStruct & vec,
4421  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4422  {
4423  // return immediately if there is nothing to do.
4424  if (exchanger.ghosts_were_set == true)
4425  return;
4426 
4427  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4428  reset_ghost_values(vec.block(i), exchanger);
4429  }
4430 
4431 
4432 
4433  // for non-block vectors
4434  template <int dim,
4435  typename VectorStruct,
4436  typename Number,
4437  typename VectorizedArrayType,
4438  std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4439  * = nullptr>
4440  inline void
4441  reset_ghost_values(
4442  const VectorStruct & vec,
4443  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4444  {
4445  exchanger.reset_ghost_values(vec);
4446  }
4447 
4448 
4449 
4450  // for std::vector of vectors
4451  template <int dim,
4452  typename VectorStruct,
4453  typename Number,
4454  typename VectorizedArrayType>
4455  inline void
4456  reset_ghost_values(
4457  const std::vector<VectorStruct> & vec,
4458  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4459  {
4460  // return immediately if there is nothing to do.
4461  if (exchanger.ghosts_were_set == true)
4462  return;
4463 
4464  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4465  reset_ghost_values(vec[comp], exchanger);
4466  }
4467 
4468 
4469 
4470  // for std::vector of pointer to vectors
4471  template <int dim,
4472  typename VectorStruct,
4473  typename Number,
4474  typename VectorizedArrayType>
4475  inline void
4476  reset_ghost_values(
4477  const std::vector<VectorStruct *> & vec,
4478  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4479  {
4480  // return immediately if there is nothing to do.
4481  if (exchanger.ghosts_were_set == true)
4482  return;
4483 
4484  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4485  reset_ghost_values(*vec[comp], exchanger);
4486  }
4487 
4488 
4489 
4490  //
4491  // zero_vector_region
4492  //
4493 
4494  // for block vectors
4495  template <int dim,
4496  typename VectorStruct,
4497  typename Number,
4498  typename VectorizedArrayType,
4499  std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4500  * = nullptr>
4501  inline void
4502  zero_vector_region(
4503  const unsigned int range_index,
4504  VectorStruct & vec,
4505  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4506  {
4507  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4508  exchanger.zero_vector_region(range_index, vec.block(i));
4509  }
4510 
4511 
4512 
4513  // for non-block vectors
4514  template <int dim,
4515  typename VectorStruct,
4516  typename Number,
4517  typename VectorizedArrayType,
4518  std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4519  * = nullptr>
4520  inline void
4521  zero_vector_region(
4522  const unsigned int range_index,
4523  VectorStruct & vec,
4524  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4525  {
4526  exchanger.zero_vector_region(range_index, vec);
4527  }
4528 
4529 
4530 
4531  // for std::vector of vectors
4532  template <int dim,
4533  typename VectorStruct,
4534  typename Number,
4535  typename VectorizedArrayType>
4536  inline void
4537  zero_vector_region(
4538  const unsigned int range_index,
4539  std::vector<VectorStruct> & vec,
4540  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4541  {
4542  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4543  zero_vector_region(range_index, vec[comp], exchanger);
4544  }
4545 
4546 
4547 
4548  // for std::vector of pointers to vectors
4549  template <int dim,
4550  typename VectorStruct,
4551  typename Number,
4552  typename VectorizedArrayType>
4553  inline void
4554  zero_vector_region(
4555  const unsigned int range_index,
4556  std::vector<VectorStruct *> & vec,
4557  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4558  {
4559  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4560  zero_vector_region(range_index, *vec[comp], exchanger);
4561  }
4562 
4563 
4564 
4565  // Apply a unit matrix operation to constrained DoFs: Default cases where we
4566  // cannot detect a LinearAlgebra::distributed::Vector, we do not do
4567  // anything, else we apply the constraints as a unit operation
4568  template <typename VectorStruct1, typename VectorStruct2>
4569  inline void
4570  apply_operation_to_constrained_dofs(const std::vector<unsigned int> &,
4571  const VectorStruct1 &,
4572  VectorStruct2 &)
4573  {}
4574 
4575  template <typename Number>
4576  inline void
4577  apply_operation_to_constrained_dofs(
4578  const std::vector<unsigned int> & constrained_dofs,
4581  {
4582  for (const unsigned int i : constrained_dofs)
4583  dst.local_element(i) = src.local_element(i);
4584  }
4585 
4586 
4587  namespace MatrixFreeFunctions
4588  {
4589  // struct to select between a const interface and a non-const interface
4590  // for MFWorker
4591  template <typename, typename, typename, typename, bool>
4592  struct InterfaceSelector
4593  {};
4594 
4595  // Version for constant functions
4596  template <typename MF,
4597  typename InVector,
4598  typename OutVector,
4599  typename Container>
4600  struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4601  {
4602  using function_type = void (Container::*)(
4603  const MF &,
4604  OutVector &,
4605  const InVector &,
4606  const std::pair<unsigned int, unsigned int> &) const;
4607  };
4608 
4609  // Version for non-constant functions
4610  template <typename MF,
4611  typename InVector,
4612  typename OutVector,
4613  typename Container>
4614  struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4615  {
4616  using function_type =
4617  void (Container::*)(const MF &,
4618  OutVector &,
4619  const InVector &,
4620  const std::pair<unsigned int, unsigned int> &);
4621  };
4622  } // namespace MatrixFreeFunctions
4623 
4624 
4625 
4626  // A implementation class for the worker object that runs the various
4627  // operations we want to perform during the matrix-free loop
4628  template <typename MF,
4629  typename InVector,
4630  typename OutVector,
4631  typename Container,
4632  bool is_constant>
4633  class MFWorker : public MFWorkerInterface
4634  {
4635  public:
4636  // An alias to make the arguments further down more readable
4637  using function_type = typename MatrixFreeFunctions::
4638  InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4639  function_type;
4640 
4641  // constructor, binds all the arguments to this class
4642  MFWorker(const MF & matrix_free,
4643  const InVector & src,
4644  OutVector & dst,
4645  const bool zero_dst_vector_setting,
4646  const Container & container,
4647  function_type cell_function,
4648  function_type face_function,
4649  function_type boundary_function,
4650  const typename MF::DataAccessOnFaces src_vector_face_access =
4652  const typename MF::DataAccessOnFaces dst_vector_face_access =
4654  const std::function<void(const unsigned int, const unsigned int)>
4655  &operation_before_loop = {},
4656  const std::function<void(const unsigned int, const unsigned int)>
4657  & operation_after_loop = {},
4658  const unsigned int dof_handler_index_pre_post = 0)
4659  : matrix_free(matrix_free)
4660  , container(const_cast<Container &>(container))
4661  , cell_function(cell_function)
4662  , face_function(face_function)
4663  , boundary_function(boundary_function)
4664  , src(src)
4665  , dst(dst)
4666  , src_data_exchanger(matrix_free,
4667  src_vector_face_access,
4668  n_components(src))
4669  , dst_data_exchanger(matrix_free,
4670  dst_vector_face_access,
4671  n_components(dst))
4672  , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4673  , zero_dst_vector_setting(zero_dst_vector_setting &&
4674  !src_and_dst_are_same)
4675  , operation_before_loop(operation_before_loop)
4676  , operation_after_loop(operation_after_loop)
4677  , dof_handler_index_pre_post(dof_handler_index_pre_post)
4678  {
4679  Assert(!has_ghost_elements(dst),
4680  ExcMessage("The destination vector passed to the matrix-free "
4681  "loop is ghosted. This is not allowed."));
4682  }
4683 
4684  // Runs the cell work. If no function is given, nothing is done
4685  virtual void
4686  cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4687  {
4688  if (cell_function != nullptr && cell_range.second > cell_range.first)
4689  for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4690  {
4691  const auto cell_subrange =
4692  matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4693 
4694  if (cell_subrange.second <= cell_subrange.first)
4695  continue;
4696 
4697  (container.*
4698  cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4699  }
4700  }
4701 
4702  virtual void
4703  cell(const unsigned int range_index) override
4704  {
4705  process_range(cell_function,
4706  matrix_free.get_task_info().cell_partition_data_hp_ptr,
4707  matrix_free.get_task_info().cell_partition_data_hp,
4708  range_index);
4709  }
4710 
4711  virtual void
4712  face(const unsigned int range_index) override
4713  {
4714  process_range(face_function,
4715  matrix_free.get_task_info().face_partition_data_hp_ptr,
4716  matrix_free.get_task_info().face_partition_data_hp,
4717  range_index);
4718  }
4719 
4720  virtual void
4721  boundary(const unsigned int range_index) override
4722  {
4723  process_range(boundary_function,
4724  matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4725  matrix_free.get_task_info().boundary_partition_data_hp,
4726  range_index);
4727  }
4728 
4729  private:
4730  void
4731  process_range(const function_type & fu,
4732  const std::vector<unsigned int> &ptr,
4733  const std::vector<unsigned int> &data,
4734  const unsigned int range_index)
4735  {
4736  if (fu == nullptr)
4737  return;
4738 
4739  AssertIndexRange(range_index + 1, ptr.size());
4740  for (unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4741  {
4742  AssertIndexRange(2 * i + 1, data.size());
4743  (container.*fu)(matrix_free,
4744  this->dst,
4745  this->src,
4746  std::make_pair(data[2 * i], data[2 * i + 1]));
4747  }
4748  }
4749 
4750  public:
4751  // Starts the communication for the update ghost values operation. We
4752  // cannot call this update if ghost and destination are the same because
4753  // that would introduce spurious entries in the destination (there is also
4754  // the problem that reading from a vector that we also write to is usually
4755  // not intended in case there is overlap, but this is up to the
4756  // application code to decide and we cannot catch this case here).
4757  virtual void
4758  vector_update_ghosts_start() override
4759  {
4760  if (!src_and_dst_are_same)
4761  internal::update_ghost_values_start(src, src_data_exchanger);
4762  }
4763 
4764  // Finishes the communication for the update ghost values operation
4765  virtual void
4766  vector_update_ghosts_finish() override
4767  {
4768  if (!src_and_dst_are_same)
4769  internal::update_ghost_values_finish(src, src_data_exchanger);
4770  }
4771 
4772  // Starts the communication for the vector compress operation
4773  virtual void
4774  vector_compress_start() override
4775  {
4776  internal::compress_start(dst, dst_data_exchanger);
4777  }
4778 
4779  // Finishes the communication for the vector compress operation
4780  virtual void
4781  vector_compress_finish() override
4782  {
4783  internal::compress_finish(dst, dst_data_exchanger);
4784  if (!src_and_dst_are_same)
4785  internal::reset_ghost_values(src, src_data_exchanger);
4786  }
4787 
4788  // Zeros the given input vector
4789  virtual void
4790  zero_dst_vector_range(const unsigned int range_index) override
4791  {
4792  if (zero_dst_vector_setting)
4793  internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4794  }
4795 
4796  virtual void
4797  cell_loop_pre_range(const unsigned int range_index) override
4798  {
4799  if (operation_before_loop)
4800  {
4801  const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4802  matrix_free.get_dof_info(dof_handler_index_pre_post);
4803  if (range_index == numbers::invalid_unsigned_int)
4804  {
4805  // Case with threaded loop -> currently no overlap implemented
4807  0U,
4808  dof_info.vector_partitioner->locally_owned_size(),
4809  operation_before_loop,
4811  }
4812  else
4813  {
4814  AssertIndexRange(range_index,
4815  dof_info.cell_loop_pre_list_index.size() - 1);
4816  for (unsigned int id =
4817  dof_info.cell_loop_pre_list_index[range_index];
4818  id != dof_info.cell_loop_pre_list_index[range_index + 1];
4819  ++id)
4820  operation_before_loop(dof_info.cell_loop_pre_list[id].first,
4821  dof_info.cell_loop_pre_list[id].second);
4822  }
4823  }
4824  }
4825 
4826  virtual void
4827  cell_loop_post_range(const unsigned int range_index) override
4828  {
4829  if (operation_after_loop)
4830  {
4831  // Run unit matrix operation on constrained dofs if we are at the
4832  // last range
4833  const std::vector<unsigned int> &partition_row_index =
4834  matrix_free.get_task_info().partition_row_index;
4835  if (range_index ==
4836  partition_row_index[partition_row_index.size() - 2] - 1)
4837  apply_operation_to_constrained_dofs(
4838  matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
4839  src,
4840  dst);
4841 
4842  const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4843  matrix_free.get_dof_info(dof_handler_index_pre_post);
4844  if (range_index == numbers::invalid_unsigned_int)
4845  {
4846  // Case with threaded loop -> currently no overlap implemented
4848  0U,
4849  dof_info.vector_partitioner->locally_owned_size(),
4850  operation_after_loop,
4852  }
4853  else
4854  {
4855  AssertIndexRange(range_index,
4856  dof_info.cell_loop_post_list_index.size() - 1);
4857  for (unsigned int id =
4858  dof_info.cell_loop_post_list_index[range_index];
4859  id != dof_info.cell_loop_post_list_index[range_index + 1];
4860  ++id)
4861  operation_after_loop(dof_info.cell_loop_post_list[id].first,
4862  dof_info.cell_loop_post_list[id].second);
4863  }
4864  }
4865  }
4866 
4867  private:
4868  const MF & matrix_free;
4869  Container & container;
4870  function_type cell_function;
4871  function_type face_function;
4872  function_type boundary_function;
4873 
4874  const InVector &src;
4875  OutVector & dst;
4876  VectorDataExchange<MF::dimension,
4877  typename MF::value_type,
4878  typename MF::vectorized_value_type>
4879  src_data_exchanger;
4880  VectorDataExchange<MF::dimension,
4881  typename MF::value_type,
4882  typename MF::vectorized_value_type>
4883  dst_data_exchanger;
4884  const bool src_and_dst_are_same;
4885  const bool zero_dst_vector_setting;
4886  const std::function<void(const unsigned int, const unsigned int)>
4887  operation_before_loop;
4888  const std::function<void(const unsigned int, const unsigned int)>
4889  operation_after_loop;
4890  const unsigned int dof_handler_index_pre_post;
4891  };
4892 
4893 
4894 
4899  template <class MF, typename InVector, typename OutVector>
4900  struct MFClassWrapper
4901  {
4902  using function_type =
4903  std::function<void(const MF &,
4904  OutVector &,
4905  const InVector &,
4906  const std::pair<unsigned int, unsigned int> &)>;
4907 
4908  MFClassWrapper(const function_type cell,
4909  const function_type face,
4910  const function_type boundary)
4911  : cell(cell)
4912  , face(face)
4913  , boundary(boundary)
4914  {}
4915 
4916  void
4917  cell_integrator(const MF & mf,
4918  OutVector & dst,
4919  const InVector & src,
4920  const std::pair<unsigned int, unsigned int> &range) const
4921  {
4922  if (cell)
4923  cell(mf, dst, src, range);
4924  }
4925 
4926  void
4927  face_integrator(const MF & mf,
4928  OutVector & dst,
4929  const InVector & src,
4930  const std::pair<unsigned int, unsigned int> &range) const
4931  {
4932  if (face)
4933  face(mf, dst, src, range);
4934  }
4935 
4936  void
4937  boundary_integrator(
4938  const MF & mf,
4939  OutVector & dst,
4940  const InVector & src,
4941  const std::pair<unsigned int, unsigned int> &range) const
4942  {
4943  if (boundary)
4944  boundary(mf, dst, src, range);
4945  }
4946 
4947  const function_type cell;
4948  const function_type face;
4949  const function_type boundary;
4950  };
4951 
4952 } // end of namespace internal
4953 
4954 
4955 
4956 template <int dim, typename Number, typename VectorizedArrayType>
4957 template <typename OutVector, typename InVector>
4958 inline void
4960  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4961  OutVector &,
4962  const InVector &,
4963  const std::pair<unsigned int, unsigned int> &)>
4964  & cell_operation,
4965  OutVector & dst,
4966  const InVector &src,
4967  const bool zero_dst_vector) const
4968 {
4969  using Wrapper =
4970  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4971  InVector,
4972  OutVector>;
4973  Wrapper wrap(cell_operation, nullptr, nullptr);
4974  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4975  InVector,
4976  OutVector,
4977  Wrapper,
4978  true>
4979  worker(*this,
4980  src,
4981  dst,
4982  zero_dst_vector,
4983  wrap,
4984  &Wrapper::cell_integrator,
4985  &Wrapper::face_integrator,
4986  &Wrapper::boundary_integrator);
4987 
4988  task_info.loop(worker);
4989 }
4990 
4991 
4992 
4993 template <int dim, typename Number, typename VectorizedArrayType>
4994 template <typename OutVector, typename InVector>
4995 inline void
4997  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4998  OutVector &,
4999  const InVector &,
5000  const std::pair<unsigned int, unsigned int> &)>
5001  & cell_operation,
5002  OutVector & dst,
5003  const InVector &src,
5004  const std::function<void(const unsigned int, const unsigned int)>
5005  &operation_before_loop,
5006  const std::function<void(const unsigned int, const unsigned int)>
5007  & operation_after_loop,
5008  const unsigned int dof_handler_index_pre_post) const
5009 {
5010  using Wrapper =
5011  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5012  InVector,
5013  OutVector>;
5014  Wrapper wrap(cell_operation, nullptr, nullptr);
5015  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5016  InVector,
5017  OutVector,
5018  Wrapper,
5019  true>
5020  worker(*this,
5021  src,
5022  dst,
5023  false,
5024  wrap,
5025  &Wrapper::cell_integrator,
5026  &Wrapper::face_integrator,
5027  &Wrapper::boundary_integrator,
5030  operation_before_loop,
5031  operation_after_loop,
5032  dof_handler_index_pre_post);
5033 
5034  task_info.loop(worker);
5035 }
5036 
5037 
5038 
5039 template <int dim, typename Number, typename VectorizedArrayType>
5040 template <typename OutVector, typename InVector>
5041 inline void
5043  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5044  OutVector &,
5045  const InVector &,
5046  const std::pair<unsigned int, unsigned int> &)>
5047  &cell_operation,
5048  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5049  OutVector &,
5050  const InVector &,
5051  const std::pair<unsigned int, unsigned int> &)>
5052  &face_operation,
5053  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5054  OutVector &,
5055  const InVector &,
5056  const std::pair<unsigned int, unsigned int> &)>
5057  & boundary_operation,
5058  OutVector & dst,
5059  const InVector & src,
5060  const bool zero_dst_vector,
5061  const DataAccessOnFaces dst_vector_face_access,
5062  const DataAccessOnFaces src_vector_face_access) const
5063 {
5064  using Wrapper =
5065  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5066  InVector,
5067  OutVector>;
5068  Wrapper wrap(cell_operation, face_operation, boundary_operation);
5069  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5070  InVector,
5071  OutVector,
5072  Wrapper,
5073  true>
5074  worker(*this,
5075  src,
5076  dst,
5077  zero_dst_vector,
5078  wrap,
5079  &Wrapper::cell_integrator,
5080  &Wrapper::face_integrator,
5081  &Wrapper::boundary_integrator,
5082  src_vector_face_access,
5083  dst_vector_face_access);
5084 
5085  task_info.loop(worker);
5086 }
5087 
5088 
5089 
5090 template <int dim, typename Number, typename VectorizedArrayType>
5091 template <typename CLASS, typename OutVector, typename InVector>
5092 inline void
5094  void (CLASS::*function_pointer)(
5096  OutVector &,
5097  const InVector &,
5098  const std::pair<unsigned int, unsigned int> &) const,
5099  const CLASS * owning_class,
5100  OutVector & dst,
5101  const InVector &src,
5102  const bool zero_dst_vector) const
5103 {
5104  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5105  InVector,
5106  OutVector,
5107  CLASS,
5108  true>
5109  worker(*this,
5110  src,
5111  dst,
5112  zero_dst_vector,
5113  *owning_class,
5114  function_pointer,
5115  nullptr,
5116  nullptr);
5117  task_info.loop(worker);
5118 }
5119 
5120 
5121 
5122 template <int dim, typename Number, typename VectorizedArrayType>
5123 template <typename CLASS, typename OutVector, typename InVector>
5124 inline void
5126  void (CLASS::*function_pointer)(
5128  OutVector &,
5129  const InVector &,
5130  const std::pair<unsigned int, unsigned int> &) const,
5131  const CLASS * owning_class,
5132  OutVector & dst,
5133  const InVector &src,
5134  const std::function<void(const unsigned int, const unsigned int)>
5135  &operation_before_loop,
5136  const std::function<void(const unsigned int, const unsigned int)>
5137  & operation_after_loop,
5138  const unsigned int dof_handler_index_pre_post) const
5139 {
5140  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5141  InVector,
5142  OutVector,
5143  CLASS,
5144  true>
5145  worker(*this,
5146  src,
5147  dst,
5148  false,
5149  *owning_class,
5150  function_pointer,
5151  nullptr,
5152  nullptr,
5155  operation_before_loop,
5156  operation_after_loop,
5157  dof_handler_index_pre_post);
5158  task_info.loop(worker);
5159 }
5160 
5161 
5162 
5163 template <int dim, typename Number, typename VectorizedArrayType>
5164 template <typename CLASS, typename OutVector, typename InVector>
5165 inline void
5167  void (CLASS::*cell_operation)(
5169  OutVector &,
5170  const InVector &,
5171  const std::pair<unsigned int, unsigned int> &) const,
5172  void (CLASS::*face_operation)(
5174  OutVector &,
5175  const InVector &,
5176  const std::pair<unsigned int, unsigned int> &) const,
5177  void (CLASS::*boundary_operation)(
5179  OutVector &,
5180  const InVector &,
5181  const std::pair<unsigned int, unsigned int> &) const,
5182  const CLASS * owning_class,
5183  OutVector & dst,
5184  const InVector & src,
5185  const bool zero_dst_vector,
5186  const DataAccessOnFaces dst_vector_face_access,
5187  const DataAccessOnFaces src_vector_face_access) const
5188 {
5189  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5190  InVector,
5191  OutVector,
5192  CLASS,
5193  true>
5194  worker(*this,
5195  src,
5196  dst,
5197  zero_dst_vector,
5198  *owning_class,
5199  cell_operation,
5200  face_operation,
5201  boundary_operation,
5202  src_vector_face_access,
5203  dst_vector_face_access);
5204  task_info.loop(worker);
5205 }
5206 
5207 
5208 
5209 template <int dim, typename Number, typename VectorizedArrayType>
5210 template <typename CLASS, typename OutVector, typename InVector>
5211 inline void
5213  void (CLASS::*function_pointer)(
5215  OutVector &,
5216  const InVector &,
5217  const std::pair<unsigned int, unsigned int> &),
5218  CLASS * owning_class,
5219  OutVector & dst,
5220  const InVector &src,
5221  const bool zero_dst_vector) const
5222 {
5223  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5224  InVector,
5225  OutVector,
5226  CLASS,
5227  false>
5228  worker(*this,
5229  src,
5230  dst,
5231  zero_dst_vector,
5232  *owning_class,
5233  function_pointer,
5234  nullptr,
5235  nullptr);
5236  task_info.loop(worker);
5237 }
5238 
5239 
5240 
5241 template <int dim, typename Number, typename VectorizedArrayType>
5242 template <typename CLASS, typename OutVector, typename InVector>
5243 inline void
5245  void (CLASS::*function_pointer)(
5247  OutVector &,
5248  const InVector &,
5249  const std::pair<unsigned int, unsigned int> &),
5250  CLASS * owning_class,
5251  OutVector & dst,
5252  const InVector &src,
5253  const std::function<void(const unsigned int, const unsigned int)>
5254  &operation_before_loop,
5255  const std::function<void(const unsigned int, const unsigned int)>
5256  & operation_after_loop,
5257  const unsigned int dof_handler_index_pre_post) const
5258 {
5259  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5260  InVector,
5261  OutVector,
5262  CLASS,
5263  false>
5264  worker(*this,
5265  src,
5266  dst,
5267  false,
5268  *owning_class,
5269  function_pointer,
5270  nullptr,
5271  nullptr,
5274  operation_before_loop,
5275  operation_after_loop,
5276  dof_handler_index_pre_post);
5277  task_info.loop(worker);
5278 }
5279 
5280 
5281 
5282 template <int dim, typename Number, typename VectorizedArrayType>
5283 template <typename CLASS, typename OutVector, typename InVector>
5284 inline void
5286  void (CLASS::*cell_operation)(
5288  OutVector &,
5289  const InVector &,
5290  const std::pair<unsigned int, unsigned int> &),
5291  void (CLASS::*face_operation)(
5293  OutVector &,
5294  const InVector &,
5295  const std::pair<unsigned int, unsigned int> &),
5296  void (CLASS::*boundary_operation)(
5298  OutVector &,
5299  const InVector &,
5300  const std::pair<unsigned int, unsigned int> &),
5301  CLASS * owning_class,
5302  OutVector & dst,
5303  const InVector & src,
5304  const bool zero_dst_vector,
5305  const DataAccessOnFaces dst_vector_face_access,
5306  const DataAccessOnFaces src_vector_face_access) const
5307 {
5308  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5309  InVector,
5310  OutVector,
5311  CLASS,
5312  false>
5313  worker(*this,
5314  src,
5315  dst,
5316  zero_dst_vector,
5317  *owning_class,
5318  cell_operation,
5319  face_operation,
5320  boundary_operation,
5321  src_vector_face_access,
5322  dst_vector_face_access);
5323  task_info.loop(worker);
5324 }
5325 
5326 
5327 
5328 template <int dim, typename Number, typename VectorizedArrayType>
5329 template <typename OutVector, typename InVector>
5330 inline void
5332  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5333  OutVector &,
5334  const InVector &,
5335  const std::pair<unsigned int, unsigned int> &)>
5336  &cell_operation,
5337  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5338  OutVector &,
5339  const InVector &,
5340  const std::pair<unsigned int, unsigned int> &)>
5341  &face_operation,
5342  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5343  OutVector &,
5344  const InVector &,
5345  const std::pair<unsigned int, unsigned int> &)>
5346  & boundary_operation,
5347  OutVector & dst,
5348  const InVector &src,
5349  const std::function<void(const unsigned int, const unsigned int)>
5350  &operation_before_loop,
5351  const std::function<void(const unsigned int, const unsigned int)>
5352  & operation_after_loop,
5353  const unsigned int dof_handler_index_pre_post,
5354  const DataAccessOnFaces dst_vector_face_access,
5355  const DataAccessOnFaces src_vector_face_access) const
5356 {
5357  using Wrapper =
5358  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5359  InVector,
5360  OutVector>;
5361  Wrapper wrap(cell_operation, face_operation, boundary_operation);
5362  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5363  InVector,
5364  OutVector,
5365  Wrapper,
5366  true>
5367  worker(*this,
5368  src,
5369  dst,
5370  false,
5371  wrap,
5372  &Wrapper::cell_integrator,
5373  &Wrapper::face_integrator,
5374  &Wrapper::boundary_integrator,
5375  src_vector_face_access,
5376  dst_vector_face_access,
5377  operation_before_loop,
5378  operation_after_loop,
5379  dof_handler_index_pre_post);
5380 
5381  task_info.loop(worker);
5382 }
5383 
5384 
5385 
5386 template <int dim, typename Number, typename VectorizedArrayType>
5387 template <typename CLASS, typename OutVector, typename InVector>
5388 inline void
5390  void (CLASS::*cell_operation)(const MatrixFree &,
5391  OutVector &,
5392  const InVector &,
5393  const std::pair<unsigned int, unsigned int> &)
5394  const,
5395  void (CLASS::*face_operation)(const MatrixFree &,
5396  OutVector &,
5397  const InVector &,
5398  const std::pair<unsigned int, unsigned int> &)
5399  const,
5400  void (CLASS::*boundary_operation)(
5401  const MatrixFree &,
5402  OutVector &,
5403  const InVector &,
5404  const std::pair<unsigned int, unsigned int> &) const,
5405  const CLASS * owning_class,
5406  OutVector & dst,
5407  const InVector &src,
5408  const std::function<void(const unsigned int, const unsigned int)>
5409  &operation_before_loop,
5410  const std::function<void(const unsigned int, const unsigned int)>
5411  & operation_after_loop,
5412  const unsigned int dof_handler_index_pre_post,
5413  const DataAccessOnFaces dst_vector_face_access,
5414  const DataAccessOnFaces src_vector_face_access) const
5415 {
5416  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5417  InVector,
5418  OutVector,
5419  CLASS,
5420  true>
5421  worker(*this,
5422  src,
5423  dst,
5424  false,
5425  *owning_class,
5426  cell_operation,
5427  face_operation,
5428  boundary_operation,
5429  src_vector_face_access,
5430  dst_vector_face_access,
5431  operation_before_loop,
5432  operation_after_loop,
5433  dof_handler_index_pre_post);
5434  task_info.loop(worker);
5435 }
5436 
5437 
5438 
5439 template <int dim, typename Number, typename VectorizedArrayType>
5440 template <typename CLASS, typename OutVector, typename InVector>
5441 inline void
5443  void (CLASS::*cell_operation)(const MatrixFree &,
5444  OutVector &,
5445  const InVector &,
5446  const std::pair<unsigned int, unsigned int> &),
5447  void (CLASS::*face_operation)(const MatrixFree &,
5448  OutVector &,
5449  const InVector &,
5450  const std::pair<unsigned int, unsigned int> &),
5451  void (CLASS::*boundary_operation)(
5452  const MatrixFree &,
5453  OutVector &,
5454  const InVector &,
5455  const std::pair<unsigned int, unsigned int> &),
5456  const CLASS * owning_class,
5457  OutVector & dst,
5458  const InVector &src,
5459  const std::function<void(const unsigned int, const unsigned int)>
5460  &operation_before_loop,
5461  const std::function<void(const unsigned int, const unsigned int)>
5462  & operation_after_loop,
5463  const unsigned int dof_handler_index_pre_post,
5464  const DataAccessOnFaces dst_vector_face_access,
5465  const DataAccessOnFaces src_vector_face_access) const
5466 {
5467  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5468  InVector,
5469  OutVector,
5470  CLASS,
5471  false>
5472  worker(*this,
5473  src,
5474  dst,
5475  false,
5476  *owning_class,
5477  cell_operation,
5478  face_operation,
5479  boundary_operation,
5480  src_vector_face_access,
5481  dst_vector_face_access,
5482  operation_before_loop,
5483  operation_after_loop,
5484  dof_handler_index_pre_post);
5485  task_info.loop(worker);
5486 }
5487 
5488 
5489 
5490 template <int dim, typename Number, typename VectorizedArrayType>
5491 template <typename CLASS, typename OutVector, typename InVector>
5492 inline void
5494  void (CLASS::*function_pointer)(
5496  OutVector &,
5497  const InVector &,
5498  const std::pair<unsigned int, unsigned int> &) const,
5499  const CLASS * owning_class,
5500  OutVector & dst,
5501  const InVector & src,
5502  const bool zero_dst_vector,
5503  const DataAccessOnFaces src_vector_face_access) const
5504 {
5505  auto src_vector_face_access_temp = src_vector_face_access;
5506  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5507  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5508  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5509  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5510 
5511  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5512  InVector,
5513  OutVector,
5514  CLASS,
5515  true>
5516  worker(*this,
5517  src,
5518  dst,
5519  zero_dst_vector,
5520  *owning_class,
5521  function_pointer,
5522  nullptr,
5523  nullptr,
5524  src_vector_face_access_temp,
5526  task_info.loop(worker);
5527 }
5528 
5529 
5530 
5531 template <int dim, typename Number, typename VectorizedArrayType>
5532 template <typename CLASS, typename OutVector, typename InVector>
5533 inline void
5535  void (CLASS::*function_pointer)(
5537  OutVector &,
5538  const InVector &,
5539  const std::pair<unsigned int, unsigned int> &),
5540  CLASS * owning_class,
5541  OutVector & dst,
5542  const InVector & src,
5543  const bool zero_dst_vector,
5544  const DataAccessOnFaces src_vector_face_access) const
5545 {
5546  auto src_vector_face_access_temp = src_vector_face_access;
5547  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5548  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5549  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5550  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5551 
5552  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5553  InVector,
5554  OutVector,
5555  CLASS,
5556  false>
5557  worker(*this,
5558  src,
5559  dst,
5560  zero_dst_vector,
5561  *owning_class,
5562  function_pointer,
5563  nullptr,
5564  nullptr,
5565  src_vector_face_access_temp,
5567  task_info.loop(worker);
5568 }
5569 
5570 
5571 
5572 template <int dim, typename Number, typename VectorizedArrayType>
5573 template <typename OutVector, typename InVector>
5574 inline void
5576  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5577  OutVector &,
5578  const InVector &,
5579  const std::pair<unsigned int, unsigned int> &)>
5580  & cell_operation,
5581  OutVector & dst,
5582  const InVector & src,
5583  const bool zero_dst_vector,
5584  const DataAccessOnFaces src_vector_face_access) const
5585 {
5586  auto src_vector_face_access_temp = src_vector_face_access;
5587  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5588  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5589  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5590  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5591 
5592  using Wrapper =
5593  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5594  InVector,
5595  OutVector>;
5596  Wrapper wrap(cell_operation, nullptr, nullptr);
5597 
5598  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5599  InVector,
5600  OutVector,
5601  Wrapper,
5602  true>
5603  worker(*this,
5604  src,
5605  dst,
5606  zero_dst_vector,
5607  wrap,
5608  &Wrapper::cell_integrator,
5609  &Wrapper::face_integrator,
5610  &Wrapper::boundary_integrator,
5611  src_vector_face_access_temp,
5613  task_info.loop(worker);
5614 }
5615 
5616 
5617 #endif // ifndef DOXYGEN
5618 
5619 
5620 
5622 
5623 #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:311
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:2277
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:2283
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:2312
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:2329
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:2352
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:2291
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:2264
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:2340
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:2347
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:2242
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:2250
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:2297
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:2256
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:2324
internal::MatrixFreeFunctions::FaceInfo< VectorizedArrayType::size()> face_info
Definition: matrix_free.h:2319
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:2270
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:2305
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:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
Point< 2 > second
Definition: grid_out.cc:4606
Point< 2 > first
Definition: grid_out.cc:4605
unsigned int level
Definition: grid_out.cc:4608
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1695
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1818
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
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:1611
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:465
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:1028
unsigned int minimum_parallel_grain_size
Definition: parallel.cc:34
const types::boundary_id invalid_boundary_id
Definition: types.h:254
static const unsigned int invalid_unsigned_int
Definition: types.h:206
void apply_to_subranges(const tbb::blocked_range< RangeType > &range, const Function &f)
Definition: parallel.h:302
unsigned short int fe_index
Definition: types.h:59
unsigned int boundary_id
Definition: types.h:134
TasksParallelScheme tasks_parallel_scheme
Definition: matrix_free.h:341
UpdateFlags mapping_update_flags_inner_faces
Definition: matrix_free.h:405
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:215
AdditionalData & operator=(const AdditionalData &other)
Definition: matrix_free.h:279
std::vector< unsigned int > cell_vectorization_category
Definition: matrix_free.h:517
UpdateFlags mapping_update_flags_boundary_faces
Definition: matrix_free.h:385
UpdateFlags mapping_update_flags
Definition: matrix_free.h:365
UpdateFlags mapping_update_flags_faces_by_cells
Definition: matrix_free.h:433
AdditionalData(const AdditionalData &other)
Definition: matrix_free.h:252
unsigned int tasks_block_size
Definition: matrix_free.h:352
static bool equal(const T *p1, const T *p2)
std::vector< unsigned int > cell_loop_pre_list_index
Definition: dof_info.h:752
std::vector< unsigned int > cell_loop_post_list_index
Definition: dof_info.h:765
std::vector< std::pair< unsigned int, unsigned int > > vector_zero_range_list
Definition: dof_info.h:745
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition: dof_info.h:596
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_pre_list
Definition: dof_info.h:758
std::vector< unsigned int > vector_zero_range_list_index
Definition: dof_info.h:740
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_post_list
Definition: dof_info.h:771
::Table< 3, unsigned int > cell_and_face_to_plain_faces
Definition: face_info.h:172
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:178
std::vector< unsigned int > boundary_partition_data
Definition: task_info.h:519
void loop(MFWorkerInterface &worker) const
Definition: task_info.cc:350
std::vector< unsigned int > cell_partition_data
Definition: task_info.h:475
std::vector< unsigned int > face_partition_data
Definition: task_info.h:497