Reference documentation for deal.II version 9.4.0
\(\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 #include <deal.II/fe/mapping_q1.h>
34 
35 #include <deal.II/hp/dof_handler.h>
38 
43 
49 
50 #include <cstdlib>
51 #include <limits>
52 #include <list>
53 #include <memory>
54 
55 
57 
58 
59 
111 template <int dim,
112  typename Number = double,
113  typename VectorizedArrayType = VectorizedArray<Number>>
114 class MatrixFree : public Subscriptor
115 {
116  static_assert(
117  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
118  "Type of Number and of VectorizedArrayType do not match.");
119 
120 public:
125  using value_type = Number;
126  using vectorized_value_type = VectorizedArrayType;
127 
131  static constexpr unsigned int dimension = dim;
132 
180  {
185 
192  {
212  };
213 
219  const unsigned int tasks_block_size = 0,
225  const unsigned int mg_level = numbers::invalid_unsigned_int,
226  const bool store_plain_indices = true,
227  const bool initialize_indices = true,
228  const bool initialize_mapping = true,
229  const bool overlap_communication_computation = true,
230  const bool hold_all_faces_to_owned_cells = false,
231  const bool cell_vectorization_categories_strict = false,
232  const bool allow_ghosted_vectors_in_loops = true)
239  , mg_level(mg_level)
248  , communicator_sm(MPI_COMM_SELF)
249  {}
250 
263  , mg_level(other.mg_level)
275  {}
276 
281  operator=(const AdditionalData &other)
282  {
291  mg_level = other.mg_level;
303 
304  return *this;
305  }
306 
344 
354  unsigned int tasks_block_size;
355 
368 
388 
408 
436 
445  unsigned int mg_level;
446 
454 
465 
474 
488 
497 
514  std::vector<unsigned int> cell_vectorization_category;
515 
526 
538 
542  MPI_Comm communicator_sm;
543  };
544 
553 
558 
562  ~MatrixFree() override = default;
563 
576  template <typename QuadratureType, typename number2, typename MappingType>
577  void
578  reinit(const MappingType & mapping,
579  const DoFHandler<dim> & dof_handler,
580  const AffineConstraints<number2> &constraint,
581  const QuadratureType & quad,
582  const AdditionalData & additional_data = AdditionalData());
583 
590  template <typename QuadratureType, typename number2>
591  DEAL_II_DEPRECATED void
592  reinit(const DoFHandler<dim> & dof_handler,
593  const AffineConstraints<number2> &constraint,
594  const QuadratureType & quad,
595  const AdditionalData & additional_data = AdditionalData());
596 
618  template <typename QuadratureType, typename number2, typename MappingType>
619  void
620  reinit(const MappingType & mapping,
621  const std::vector<const DoFHandler<dim> *> & dof_handler,
622  const std::vector<const AffineConstraints<number2> *> &constraint,
623  const std::vector<QuadratureType> & quad,
624  const AdditionalData &additional_data = AdditionalData());
625 
631  template <typename QuadratureType,
632  typename number2,
633  typename DoFHandlerType,
634  typename MappingType>
635  DEAL_II_DEPRECATED void
636  reinit(const MappingType & mapping,
637  const std::vector<const DoFHandlerType *> & dof_handler,
638  const std::vector<const AffineConstraints<number2> *> &constraint,
639  const std::vector<QuadratureType> & quad,
640  const AdditionalData &additional_data = AdditionalData());
641 
648  template <typename QuadratureType, typename number2>
649  DEAL_II_DEPRECATED void
650  reinit(const std::vector<const DoFHandler<dim> *> & dof_handler,
651  const std::vector<const AffineConstraints<number2> *> &constraint,
652  const std::vector<QuadratureType> & quad,
653  const AdditionalData &additional_data = AdditionalData());
654 
660  template <typename QuadratureType, typename number2, typename DoFHandlerType>
661  DEAL_II_DEPRECATED void
662  reinit(const std::vector<const DoFHandlerType *> & dof_handler,
663  const std::vector<const AffineConstraints<number2> *> &constraint,
664  const std::vector<QuadratureType> & quad,
665  const AdditionalData &additional_data = AdditionalData());
666 
674  template <typename QuadratureType, typename number2, typename MappingType>
675  void
676  reinit(const MappingType & mapping,
677  const std::vector<const DoFHandler<dim> *> & dof_handler,
678  const std::vector<const AffineConstraints<number2> *> &constraint,
679  const QuadratureType & quad,
680  const AdditionalData &additional_data = AdditionalData());
681 
687  template <typename QuadratureType,
688  typename number2,
689  typename DoFHandlerType,
690  typename MappingType>
691  DEAL_II_DEPRECATED void
692  reinit(const MappingType & mapping,
693  const std::vector<const DoFHandlerType *> & dof_handler,
694  const std::vector<const AffineConstraints<number2> *> &constraint,
695  const QuadratureType & quad,
696  const AdditionalData &additional_data = AdditionalData());
697 
704  template <typename QuadratureType, typename number2>
705  DEAL_II_DEPRECATED void
706  reinit(const std::vector<const DoFHandler<dim> *> & dof_handler,
707  const std::vector<const AffineConstraints<number2> *> &constraint,
708  const QuadratureType & quad,
709  const AdditionalData &additional_data = AdditionalData());
710 
716  template <typename QuadratureType, typename number2, typename DoFHandlerType>
717  DEAL_II_DEPRECATED void
718  reinit(const std::vector<const DoFHandlerType *> & dof_handler,
719  const std::vector<const AffineConstraints<number2> *> &constraint,
720  const QuadratureType & quad,
721  const AdditionalData &additional_data = AdditionalData());
722 
728  void
730  const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free_base);
731 
741  void
742  update_mapping(const Mapping<dim> &mapping);
743 
747  void
748  update_mapping(const std::shared_ptr<hp::MappingCollection<dim>> &mapping);
749 
754  void
755  clear();
756 
758 
770  enum class DataAccessOnFaces
771  {
778  none,
779 
790  values,
791 
801 
812  gradients,
813 
823 
831  };
832 
877  template <typename OutVector, typename InVector>
878  void
879  cell_loop(const std::function<void(
881  OutVector &,
882  const InVector &,
883  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
884  OutVector & dst,
885  const InVector & src,
886  const bool zero_dst_vector = false) const;
887 
934  template <typename CLASS, typename OutVector, typename InVector>
935  void
936  cell_loop(void (CLASS::*cell_operation)(
937  const MatrixFree &,
938  OutVector &,
939  const InVector &,
940  const std::pair<unsigned int, unsigned int> &) const,
941  const CLASS * owning_class,
942  OutVector & dst,
943  const InVector &src,
944  const bool zero_dst_vector = false) const;
945 
949  template <typename CLASS, typename OutVector, typename InVector>
950  void
951  cell_loop(void (CLASS::*cell_operation)(
952  const MatrixFree &,
953  OutVector &,
954  const InVector &,
955  const std::pair<unsigned int, unsigned int> &),
956  CLASS * owning_class,
957  OutVector & dst,
958  const InVector &src,
959  const bool zero_dst_vector = false) const;
960 
1045  template <typename CLASS, typename OutVector, typename InVector>
1046  void
1047  cell_loop(void (CLASS::*cell_operation)(
1048  const MatrixFree &,
1049  OutVector &,
1050  const InVector &,
1051  const std::pair<unsigned int, unsigned int> &) const,
1052  const CLASS * owning_class,
1053  OutVector & dst,
1054  const InVector &src,
1055  const std::function<void(const unsigned int, const unsigned int)>
1056  &operation_before_loop,
1057  const std::function<void(const unsigned int, const unsigned int)>
1058  & operation_after_loop,
1059  const unsigned int dof_handler_index_pre_post = 0) const;
1060 
1064  template <typename CLASS, typename OutVector, typename InVector>
1065  void
1066  cell_loop(void (CLASS::*cell_operation)(
1067  const MatrixFree &,
1068  OutVector &,
1069  const InVector &,
1070  const std::pair<unsigned int, unsigned int> &),
1071  CLASS * owning_class,
1072  OutVector & dst,
1073  const InVector &src,
1074  const std::function<void(const unsigned int, const unsigned int)>
1075  &operation_before_loop,
1076  const std::function<void(const unsigned int, const unsigned int)>
1077  & operation_after_loop,
1078  const unsigned int dof_handler_index_pre_post = 0) const;
1079 
1084  template <typename OutVector, typename InVector>
1085  void
1086  cell_loop(const std::function<void(
1088  OutVector &,
1089  const InVector &,
1090  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1091  OutVector & dst,
1092  const InVector & src,
1093  const std::function<void(const unsigned int, const unsigned int)>
1094  &operation_before_loop,
1095  const std::function<void(const unsigned int, const unsigned int)>
1096  & operation_after_loop,
1097  const unsigned int dof_handler_index_pre_post = 0) const;
1098 
1174  template <typename OutVector, typename InVector>
1175  void
1176  loop(const std::function<
1178  OutVector &,
1179  const InVector &,
1180  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1181  const std::function<
1183  OutVector &,
1184  const InVector &,
1185  const std::pair<unsigned int, unsigned int> &)> &face_operation,
1186  const std::function<void(
1188  OutVector &,
1189  const InVector &,
1190  const std::pair<unsigned int, unsigned int> &)> &boundary_operation,
1191  OutVector & dst,
1192  const InVector & src,
1193  const bool zero_dst_vector = false,
1194  const DataAccessOnFaces dst_vector_face_access =
1196  const DataAccessOnFaces src_vector_face_access =
1198 
1284  template <typename CLASS, typename OutVector, typename InVector>
1285  void
1287  void (CLASS::*cell_operation)(const MatrixFree &,
1288  OutVector &,
1289  const InVector &,
1290  const std::pair<unsigned int, unsigned int> &)
1291  const,
1292  void (CLASS::*face_operation)(const MatrixFree &,
1293  OutVector &,
1294  const InVector &,
1295  const std::pair<unsigned int, unsigned int> &)
1296  const,
1297  void (CLASS::*boundary_operation)(
1298  const MatrixFree &,
1299  OutVector &,
1300  const InVector &,
1301  const std::pair<unsigned int, unsigned int> &) const,
1302  const CLASS * owning_class,
1303  OutVector & dst,
1304  const InVector & src,
1305  const bool zero_dst_vector = false,
1306  const DataAccessOnFaces dst_vector_face_access =
1308  const DataAccessOnFaces src_vector_face_access =
1310 
1314  template <typename CLASS, typename OutVector, typename InVector>
1315  void
1316  loop(void (CLASS::*cell_operation)(
1317  const MatrixFree &,
1318  OutVector &,
1319  const InVector &,
1320  const std::pair<unsigned int, unsigned int> &),
1321  void (CLASS::*face_operation)(
1322  const MatrixFree &,
1323  OutVector &,
1324  const InVector &,
1325  const std::pair<unsigned int, unsigned int> &),
1326  void (CLASS::*boundary_operation)(
1327  const MatrixFree &,
1328  OutVector &,
1329  const InVector &,
1330  const std::pair<unsigned int, unsigned int> &),
1331  CLASS * owning_class,
1332  OutVector & dst,
1333  const InVector & src,
1334  const bool zero_dst_vector = false,
1335  const DataAccessOnFaces dst_vector_face_access =
1337  const DataAccessOnFaces src_vector_face_access =
1339 
1404  template <typename CLASS, typename OutVector, typename InVector>
1405  void
1406  loop_cell_centric(void (CLASS::*cell_operation)(
1407  const MatrixFree &,
1408  OutVector &,
1409  const InVector &,
1410  const std::pair<unsigned int, unsigned int> &) const,
1411  const CLASS * owning_class,
1412  OutVector & dst,
1413  const InVector & src,
1414  const bool zero_dst_vector = false,
1415  const DataAccessOnFaces src_vector_face_access =
1417 
1421  template <typename CLASS, typename OutVector, typename InVector>
1422  void
1423  loop_cell_centric(void (CLASS::*cell_operation)(
1424  const MatrixFree &,
1425  OutVector &,
1426  const InVector &,
1427  const std::pair<unsigned int, unsigned int> &),
1428  CLASS * owning_class,
1429  OutVector & dst,
1430  const InVector & src,
1431  const bool zero_dst_vector = false,
1432  const DataAccessOnFaces src_vector_face_access =
1434 
1438  template <typename OutVector, typename InVector>
1439  void
1441  const std::function<void(const MatrixFree &,
1442  OutVector &,
1443  const InVector &,
1444  const std::pair<unsigned int, unsigned int> &)>
1445  & cell_operation,
1446  OutVector & dst,
1447  const InVector & src,
1448  const bool zero_dst_vector = false,
1449  const DataAccessOnFaces src_vector_face_access =
1451 
1459  std::pair<unsigned int, unsigned int>
1460  create_cell_subrange_hp(const std::pair<unsigned int, unsigned int> &range,
1461  const unsigned int fe_degree,
1462  const unsigned int dof_handler_index = 0) const;
1463 
1470  std::pair<unsigned int, unsigned int>
1472  const std::pair<unsigned int, unsigned int> &range,
1473  const unsigned int fe_index,
1474  const unsigned int dof_handler_index = 0) const;
1475 
1479  unsigned int
1481 
1485  unsigned int
1487  const std::pair<unsigned int, unsigned int> range) const;
1488 
1492  unsigned int
1493  get_face_active_fe_index(const std::pair<unsigned int, unsigned int> range,
1494  const bool is_interior_face = true) const;
1495 
1497 
1507  template <typename T>
1508  void
1510 
1516  template <typename T>
1517  void
1519 
1533  template <typename VectorType>
1534  void
1535  initialize_dof_vector(VectorType & vec,
1536  const unsigned int dof_handler_index = 0) const;
1537 
1555  template <typename Number2>
1556  void
1558  const unsigned int dof_handler_index = 0) const;
1559 
1570  const std::shared_ptr<const Utilities::MPI::Partitioner> &
1571  get_vector_partitioner(const unsigned int dof_handler_index = 0) const;
1572 
1576  const IndexSet &
1577  get_locally_owned_set(const unsigned int dof_handler_index = 0) const;
1578 
1582  const IndexSet &
1583  get_ghost_set(const unsigned int dof_handler_index = 0) const;
1584 
1594  const std::vector<unsigned int> &
1595  get_constrained_dofs(const unsigned int dof_handler_index = 0) const;
1596 
1607  void
1608  renumber_dofs(std::vector<types::global_dof_index> &renumbering,
1609  const unsigned int dof_handler_index = 0);
1610 
1612 
1620  template <int spacedim>
1621  static bool
1623 
1627  unsigned int
1628  n_components() const;
1629 
1634  unsigned int
1635  n_base_elements(const unsigned int dof_handler_index) const;
1636 
1644  unsigned int
1646 
1650  DEAL_II_DEPRECATED unsigned int
1651  n_macro_cells() const;
1652 
1662  unsigned int
1664 
1671  unsigned int
1673 
1683  unsigned int
1685 
1696  unsigned int
1698 
1703  unsigned int
1705 
1713  get_boundary_id(const unsigned int face_batch_index) const;
1714 
1719  std::array<types::boundary_id, VectorizedArrayType::size()>
1720  get_faces_by_cells_boundary_id(const unsigned int cell_batch_index,
1721  const unsigned int face_number) const;
1722 
1727  const DoFHandler<dim> &
1728  get_dof_handler(const unsigned int dof_handler_index = 0) const;
1729 
1740  template <typename DoFHandlerType>
1741  DEAL_II_DEPRECATED const DoFHandlerType &
1742  get_dof_handler(const unsigned int dof_handler_index = 0) const;
1743 
1757  get_cell_iterator(const unsigned int cell_batch_index,
1758  const unsigned int lane_index,
1759  const unsigned int dof_handler_index = 0) const;
1760 
1766  std::pair<int, int>
1767  get_cell_level_and_index(const unsigned int cell_batch_index,
1768  const unsigned int lane_index) const;
1769 
1782  std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>
1783  get_face_iterator(const unsigned int face_batch_index,
1784  const unsigned int lane_index,
1785  const bool interior = true,
1786  const unsigned int fe_component = 0) const;
1787 
1794  get_hp_cell_iterator(const unsigned int cell_batch_index,
1795  const unsigned int lane_index,
1796  const unsigned int dof_handler_index = 0) const;
1797 
1810  bool
1811  at_irregular_cell(const unsigned int cell_batch_index) const;
1812 
1816  DEAL_II_DEPRECATED unsigned int
1817  n_components_filled(const unsigned int cell_batch_number) const;
1818 
1828  unsigned int
1829  n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const;
1830 
1840  unsigned int
1841  n_active_entries_per_face_batch(const unsigned int face_batch_index) const;
1842 
1846  unsigned int
1847  get_dofs_per_cell(const unsigned int dof_handler_index = 0,
1848  const unsigned int hp_active_fe_index = 0) const;
1849 
1853  unsigned int
1854  get_n_q_points(const unsigned int quad_index = 0,
1855  const unsigned int hp_active_fe_index = 0) const;
1856 
1861  unsigned int
1862  get_dofs_per_face(const unsigned int dof_handler_index = 0,
1863  const unsigned int hp_active_fe_index = 0) const;
1864 
1869  unsigned int
1870  get_n_q_points_face(const unsigned int quad_index = 0,
1871  const unsigned int hp_active_fe_index = 0) const;
1872 
1876  const Quadrature<dim> &
1877  get_quadrature(const unsigned int quad_index = 0,
1878  const unsigned int hp_active_fe_index = 0) const;
1879 
1883  const Quadrature<dim - 1> &
1884  get_face_quadrature(const unsigned int quad_index = 0,
1885  const unsigned int hp_active_fe_index = 0) const;
1886 
1899  unsigned int
1901  const std::pair<unsigned int, unsigned int> cell_batch_range) const;
1902 
1907  std::pair<unsigned int, unsigned int>
1909  const std::pair<unsigned int, unsigned int> face_batch_range) const;
1910 
1922  unsigned int
1923  get_cell_category(const unsigned int cell_batch_index) const;
1924 
1929  std::pair<unsigned int, unsigned int>
1930  get_face_category(const unsigned int face_batch_index) const;
1931 
1935  bool
1937 
1942  bool
1944 
1949  unsigned int
1950  get_mg_level() const;
1951 
1956  std::size_t
1958 
1963  template <typename StreamType>
1964  void
1965  print_memory_consumption(StreamType &out) const;
1966 
1971  void
1972  print(std::ostream &out) const;
1973 
1975 
1986  get_task_info() const;
1987 
1988  /*
1989  * Return geometry-dependent information on the cells.
1990  */
1991  const internal::MatrixFreeFunctions::
1992  MappingInfo<dim, Number, VectorizedArrayType> &
1994 
1999  get_dof_info(const unsigned int dof_handler_index_component = 0) const;
2000 
2004  unsigned int
2006 
2011  const Number *
2012  constraint_pool_begin(const unsigned int pool_index) const;
2013 
2019  const Number *
2020  constraint_pool_end(const unsigned int pool_index) const;
2021 
2026  get_shape_info(const unsigned int dof_handler_index_component = 0,
2027  const unsigned int quad_index = 0,
2028  const unsigned int fe_base_element = 0,
2029  const unsigned int hp_active_fe_index = 0,
2030  const unsigned int hp_active_quad_index = 0) const;
2031 
2036  VectorizedArrayType::size()> &
2037  get_face_info(const unsigned int face_batch_index) const;
2038 
2039 
2045  const Table<3, unsigned int> &
2047 
2063 
2067  void
2069 
2081 
2085  void
2087  const AlignedVector<Number> *memory) const;
2088 
2090 
2091 private:
2096  template <typename number2, int q_dim>
2097  void
2099  const std::shared_ptr<hp::MappingCollection<dim>> & mapping,
2100  const std::vector<const DoFHandler<dim, dim> *> & dof_handlers,
2101  const std::vector<const AffineConstraints<number2> *> &constraint,
2102  const std::vector<IndexSet> & locally_owned_set,
2103  const std::vector<hp::QCollection<q_dim>> & quad,
2104  const AdditionalData & additional_data);
2105 
2112  template <typename number2>
2113  void
2115  const std::vector<const AffineConstraints<number2> *> &constraint,
2116  const std::vector<IndexSet> & locally_owned_set,
2117  const AdditionalData & additional_data);
2118 
2122  void
2124  const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
2125  const AdditionalData & additional_data);
2126 
2130  std::vector<SmartPointer<const DoFHandler<dim>>> dof_handlers;
2131 
2136  std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
2137 
2144  std::vector<Number> constraint_pool_data;
2145 
2150  std::vector<unsigned int> constraint_pool_row_index;
2151 
2158 
2164 
2171  std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
2172 
2173 
2181 
2188 
2193  internal::MatrixFreeFunctions::FaceInfo<VectorizedArrayType::size()>
2195 
2200 
2205 
2214  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2216 
2221  mutable std::list<std::pair<bool, AlignedVector<Number>>>
2223 
2227  unsigned int mg_level;
2228 };
2229 
2230 
2231 
2232 /*----------------------- Inline functions ----------------------------------*/
2233 
2234 #ifndef DOXYGEN
2235 
2236 
2237 
2238 template <int dim, typename Number, typename VectorizedArrayType>
2239 template <typename T>
2240 inline void
2242  AlignedVector<T> &vec) const
2243 {
2244  vec.resize(this->n_cell_batches() + this->n_ghost_cell_batches());
2245 }
2246 
2247 
2248 
2249 template <int dim, typename Number, typename VectorizedArrayType>
2250 template <typename T>
2251 inline void
2253  AlignedVector<T> &vec) const
2254 {
2255  vec.resize(this->n_inner_face_batches() + this->n_boundary_face_batches() +
2256  this->n_ghost_inner_face_batches());
2257 }
2258 
2259 
2260 
2261 template <int dim, typename Number, typename VectorizedArrayType>
2262 template <typename VectorType>
2263 inline void
2265  VectorType & vec,
2266  const unsigned int comp) const
2267 {
2268  static_assert(IsBlockVector<VectorType>::value == false,
2269  "This function is not supported for block vectors.");
2270 
2271  Assert(task_info.n_procs == 1,
2272  ExcMessage("This function can only be used in serial."));
2273 
2274  AssertIndexRange(comp, n_components());
2275  vec.reinit(dof_info[comp].vector_partitioner->size());
2276 }
2277 
2278 
2279 
2280 template <int dim, typename Number, typename VectorizedArrayType>
2281 template <typename Number2>
2282 inline void
2285  const unsigned int comp) const
2286 {
2287  AssertIndexRange(comp, n_components());
2288  vec.reinit(dof_info[comp].vector_partitioner, task_info.communicator_sm);
2289 }
2290 
2291 
2292 
2293 template <int dim, typename Number, typename VectorizedArrayType>
2294 inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2296  const unsigned int comp) const
2297 {
2298  AssertIndexRange(comp, n_components());
2299  return dof_info[comp].vector_partitioner;
2300 }
2301 
2302 
2303 
2304 template <int dim, typename Number, typename VectorizedArrayType>
2305 inline const std::vector<unsigned int> &
2307  const unsigned int comp) const
2308 {
2309  AssertIndexRange(comp, n_components());
2310  return dof_info[comp].constrained_dofs;
2311 }
2312 
2313 
2314 
2315 template <int dim, typename Number, typename VectorizedArrayType>
2316 inline unsigned int
2318 {
2319  AssertDimension(dof_handlers.size(), dof_info.size());
2320  return dof_handlers.size();
2321 }
2322 
2323 
2324 
2325 template <int dim, typename Number, typename VectorizedArrayType>
2326 inline unsigned int
2328  const unsigned int dof_no) const
2329 {
2330  AssertDimension(dof_handlers.size(), dof_info.size());
2331  AssertIndexRange(dof_no, dof_handlers.size());
2332  return dof_handlers[dof_no]->get_fe().n_base_elements();
2333 }
2334 
2335 
2336 
2337 template <int dim, typename Number, typename VectorizedArrayType>
2340 {
2341  return task_info;
2342 }
2343 
2344 
2345 
2346 template <int dim, typename Number, typename VectorizedArrayType>
2347 inline unsigned int
2349 {
2350  return *(task_info.cell_partition_data.end() - 2);
2351 }
2352 
2353 
2354 
2355 template <int dim, typename Number, typename VectorizedArrayType>
2356 inline unsigned int
2358 {
2359  return task_info.n_active_cells;
2360 }
2361 
2362 
2363 
2364 template <int dim, typename Number, typename VectorizedArrayType>
2365 inline unsigned int
2367 {
2368  return *(task_info.cell_partition_data.end() - 2);
2369 }
2370 
2371 
2372 
2373 template <int dim, typename Number, typename VectorizedArrayType>
2374 inline unsigned int
2376 {
2377  return *(task_info.cell_partition_data.end() - 1) -
2378  *(task_info.cell_partition_data.end() - 2);
2379 }
2380 
2381 
2382 
2383 template <int dim, typename Number, typename VectorizedArrayType>
2384 inline unsigned int
2386 {
2387  if (task_info.face_partition_data.size() == 0)
2388  return 0;
2389  return task_info.face_partition_data.back();
2390 }
2391 
2392 
2393 
2394 template <int dim, typename Number, typename VectorizedArrayType>
2395 inline unsigned int
2397 {
2398  if (task_info.face_partition_data.size() == 0)
2399  return 0;
2400  return task_info.boundary_partition_data.back() -
2402 }
2403 
2404 
2405 
2406 template <int dim, typename Number, typename VectorizedArrayType>
2407 inline unsigned int
2409 {
2410  if (task_info.face_partition_data.size() == 0)
2411  return 0;
2412  return face_info.faces.size() - task_info.boundary_partition_data.back();
2413 }
2414 
2415 
2416 
2417 template <int dim, typename Number, typename VectorizedArrayType>
2418 inline types::boundary_id
2420  const unsigned int face_batch_index) const
2421 {
2422  Assert(face_batch_index >= task_info.boundary_partition_data[0] &&
2423  face_batch_index < task_info.boundary_partition_data.back(),
2424  ExcIndexRange(face_batch_index,
2427  return types::boundary_id(face_info.faces[face_batch_index].exterior_face_no);
2428 }
2429 
2430 
2431 
2432 template <int dim, typename Number, typename VectorizedArrayType>
2433 inline std::array<types::boundary_id, VectorizedArrayType::size()>
2435  const unsigned int cell_batch_index,
2436  const unsigned int face_number) const
2437 {
2438  AssertIndexRange(cell_batch_index, n_cell_batches());
2441  ExcNotInitialized());
2442  std::array<types::boundary_id, VectorizedArrayType::size()> result;
2443  result.fill(numbers::invalid_boundary_id);
2444  for (unsigned int v = 0;
2445  v < n_active_entries_per_cell_batch(cell_batch_index);
2446  ++v)
2447  result[v] =
2448  face_info.cell_and_face_boundary_id(cell_batch_index, face_number, v);
2449  return result;
2450 }
2451 
2452 
2453 
2454 template <int dim, typename Number, typename VectorizedArrayType>
2455 inline const internal::MatrixFreeFunctions::
2456  MappingInfo<dim, Number, VectorizedArrayType> &
2458 {
2459  return mapping_info;
2460 }
2461 
2462 
2463 
2464 template <int dim, typename Number, typename VectorizedArrayType>
2467  const unsigned int dof_index) const
2468 {
2469  AssertIndexRange(dof_index, n_components());
2470  return dof_info[dof_index];
2471 }
2472 
2473 
2474 
2475 template <int dim, typename Number, typename VectorizedArrayType>
2476 inline unsigned int
2478 {
2479  return constraint_pool_row_index.size() - 1;
2480 }
2481 
2482 
2483 
2484 template <int dim, typename Number, typename VectorizedArrayType>
2485 inline const Number *
2487  const unsigned int row) const
2488 {
2490  return constraint_pool_data.empty() ?
2491  nullptr :
2493 }
2494 
2495 
2496 
2497 template <int dim, typename Number, typename VectorizedArrayType>
2498 inline const Number *
2500  const unsigned int row) const
2501 {
2503  return constraint_pool_data.empty() ?
2504  nullptr :
2506 }
2507 
2508 
2509 
2510 template <int dim, typename Number, typename VectorizedArrayType>
2511 inline std::pair<unsigned int, unsigned int>
2513  const std::pair<unsigned int, unsigned int> &range,
2514  const unsigned int degree,
2515  const unsigned int dof_handler_component) const
2516 {
2517  if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2518  {
2520  dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2522  dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2523  if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2524  return range;
2525  else
2526  return {range.second, range.second};
2527  }
2528 
2529  const unsigned int fe_index =
2530  dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2531  if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2532  return {range.second, range.second};
2533  else
2534  return create_cell_subrange_hp_by_index(range,
2535  fe_index,
2536  dof_handler_component);
2537 }
2538 
2539 
2540 
2541 template <int dim, typename Number, typename VectorizedArrayType>
2542 inline bool
2544  const unsigned int cell_batch_index) const
2545 {
2546  AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2547  return VectorizedArrayType::size() > 1 &&
2548  cell_level_index[(cell_batch_index + 1) * VectorizedArrayType::size() -
2549  1] == cell_level_index[(cell_batch_index + 1) *
2550  VectorizedArrayType::size() -
2551  2];
2552 }
2553 
2554 
2555 
2556 template <int dim, typename Number, typename VectorizedArrayType>
2557 unsigned int
2559 {
2560  return shape_info.size(2);
2561 }
2562 
2563 
2564 template <int dim, typename Number, typename VectorizedArrayType>
2565 unsigned int
2567  const std::pair<unsigned int, unsigned int> range) const
2568 {
2569  const auto &fe_indices = dof_info[0].cell_active_fe_index;
2570 
2571  if (fe_indices.empty() == true)
2572  return 0;
2573 
2574  const auto index = fe_indices[range.first];
2575 
2576  for (unsigned int i = range.first; i < range.second; ++i)
2577  AssertDimension(index, fe_indices[i]);
2578 
2579  return index;
2580 }
2581 
2582 
2583 
2584 template <int dim, typename Number, typename VectorizedArrayType>
2585 unsigned int
2587  const std::pair<unsigned int, unsigned int> range,
2588  const bool is_interior_face) const
2589 {
2590  const auto &fe_indices = dof_info[0].cell_active_fe_index;
2591 
2592  if (fe_indices.empty() == true)
2593  return 0;
2594 
2595  if (is_interior_face)
2596  {
2597  const unsigned int index =
2598  fe_indices[face_info.faces[range.first].cells_interior[0] /
2599  VectorizedArrayType::size()];
2600 
2601  for (unsigned int i = range.first; i < range.second; ++i)
2602  AssertDimension(index,
2603  fe_indices[face_info.faces[i].cells_interior[0] /
2604  VectorizedArrayType::size()]);
2605 
2606  return index;
2607  }
2608  else
2609  {
2610  const unsigned int index =
2611  fe_indices[face_info.faces[range.first].cells_exterior[0] /
2612  VectorizedArrayType::size()];
2613 
2614  for (unsigned int i = range.first; i < range.second; ++i)
2615  AssertDimension(index,
2616  fe_indices[face_info.faces[i].cells_exterior[0] /
2617  VectorizedArrayType::size()]);
2618 
2619  return index;
2620  }
2621 }
2622 
2623 
2624 
2625 template <int dim, typename Number, typename VectorizedArrayType>
2626 inline unsigned int
2628  const unsigned int cell_batch_index) const
2629 {
2630  return n_active_entries_per_cell_batch(cell_batch_index);
2631 }
2632 
2633 
2634 
2635 template <int dim, typename Number, typename VectorizedArrayType>
2636 inline unsigned int
2638  const unsigned int cell_batch_index) const
2639 {
2640  Assert(!dof_info.empty(), ExcNotInitialized());
2641  AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2642  const std::vector<unsigned char> &n_lanes_filled =
2643  dof_info[0].n_vectorization_lanes_filled
2645  AssertIndexRange(cell_batch_index, n_lanes_filled.size());
2646 
2647  return n_lanes_filled[cell_batch_index];
2648 }
2649 
2650 
2651 
2652 template <int dim, typename Number, typename VectorizedArrayType>
2653 inline unsigned int
2655  const unsigned int face_batch_index) const
2656 {
2657  AssertIndexRange(face_batch_index, face_info.faces.size());
2658  Assert(!dof_info.empty(), ExcNotInitialized());
2659  const std::vector<unsigned char> &n_lanes_filled =
2660  dof_info[0].n_vectorization_lanes_filled
2662  AssertIndexRange(face_batch_index, n_lanes_filled.size());
2663  return n_lanes_filled[face_batch_index];
2664 }
2665 
2666 
2667 
2668 template <int dim, typename Number, typename VectorizedArrayType>
2669 inline unsigned int
2671  const unsigned int dof_handler_index,
2672  const unsigned int active_fe_index) const
2673 {
2674  return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2675 }
2676 
2677 
2678 
2679 template <int dim, typename Number, typename VectorizedArrayType>
2680 inline unsigned int
2682  const unsigned int quad_index,
2683  const unsigned int active_fe_index) const
2684 {
2685  AssertIndexRange(quad_index, mapping_info.cell_data.size());
2686  return mapping_info.cell_data[quad_index]
2687  .descriptor[active_fe_index]
2688  .n_q_points;
2689 }
2690 
2691 
2692 
2693 template <int dim, typename Number, typename VectorizedArrayType>
2694 inline unsigned int
2696  const unsigned int dof_handler_index,
2697  const unsigned int active_fe_index) const
2698 {
2699  return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2700 }
2701 
2702 
2703 
2704 template <int dim, typename Number, typename VectorizedArrayType>
2705 inline unsigned int
2707  const unsigned int quad_index,
2708  const unsigned int active_fe_index) const
2709 {
2710  AssertIndexRange(quad_index, mapping_info.face_data.size());
2711  return mapping_info.face_data[quad_index]
2712  .descriptor[active_fe_index]
2713  .n_q_points;
2714 }
2715 
2716 
2717 
2718 template <int dim, typename Number, typename VectorizedArrayType>
2719 inline const IndexSet &
2721  const unsigned int dof_handler_index) const
2722 {
2723  return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2724 }
2725 
2726 
2727 
2728 template <int dim, typename Number, typename VectorizedArrayType>
2729 inline const IndexSet &
2731  const unsigned int dof_handler_index) const
2732 {
2733  return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2734 }
2735 
2736 
2737 
2738 template <int dim, typename Number, typename VectorizedArrayType>
2741  const unsigned int dof_handler_index,
2742  const unsigned int index_quad,
2743  const unsigned int index_fe,
2744  const unsigned int active_fe_index,
2745  const unsigned int active_quad_index) const
2746 {
2747  AssertIndexRange(dof_handler_index, dof_info.size());
2748  const unsigned int ind =
2749  dof_info[dof_handler_index].global_base_element_offset + index_fe;
2750  AssertIndexRange(ind, shape_info.size(0));
2751  AssertIndexRange(index_quad, shape_info.size(1));
2752  AssertIndexRange(active_fe_index, shape_info.size(2));
2753  AssertIndexRange(active_quad_index, shape_info.size(3));
2754  return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2755 }
2756 
2757 
2758 
2759 template <int dim, typename Number, typename VectorizedArrayType>
2761  VectorizedArrayType::size()> &
2763  const unsigned int face_batch_index) const
2764 {
2765  AssertIndexRange(face_batch_index, face_info.faces.size());
2766  return face_info.faces[face_batch_index];
2767 }
2768 
2769 
2770 
2771 template <int dim, typename Number, typename VectorizedArrayType>
2772 inline const Table<3, unsigned int> &
2774  const
2775 {
2777 }
2778 
2779 
2780 
2781 template <int dim, typename Number, typename VectorizedArrayType>
2782 inline const Quadrature<dim> &
2784  const unsigned int quad_index,
2785  const unsigned int active_fe_index) const
2786 {
2787  AssertIndexRange(quad_index, mapping_info.cell_data.size());
2788  return mapping_info.cell_data[quad_index]
2789  .descriptor[active_fe_index]
2790  .quadrature;
2791 }
2792 
2793 
2794 
2795 template <int dim, typename Number, typename VectorizedArrayType>
2796 inline const Quadrature<dim - 1> &
2798  const unsigned int quad_index,
2799  const unsigned int active_fe_index) const
2800 {
2801  AssertIndexRange(quad_index, mapping_info.face_data.size());
2802  return mapping_info.face_data[quad_index]
2803  .descriptor[active_fe_index]
2804  .quadrature;
2805 }
2806 
2807 
2808 
2809 template <int dim, typename Number, typename VectorizedArrayType>
2810 inline unsigned int
2812  const std::pair<unsigned int, unsigned int> range) const
2813 {
2814  auto result = get_cell_category(range.first);
2815 
2816  for (unsigned int i = range.first; i < range.second; ++i)
2817  result = std::max(result, get_cell_category(i));
2818 
2819  return result;
2820 }
2821 
2822 
2823 
2824 template <int dim, typename Number, typename VectorizedArrayType>
2825 inline std::pair<unsigned int, unsigned int>
2827  const std::pair<unsigned int, unsigned int> range) const
2828 {
2829  auto result = get_face_category(range.first);
2830 
2831  for (unsigned int i = range.first; i < range.second; ++i)
2832  {
2833  result.first = std::max(result.first, get_face_category(i).first);
2834  result.second = std::max(result.second, get_face_category(i).second);
2835  }
2836 
2837  return result;
2838 }
2839 
2840 
2841 
2842 template <int dim, typename Number, typename VectorizedArrayType>
2843 inline unsigned int
2845  const unsigned int cell_batch_index) const
2846 {
2847  AssertIndexRange(0, dof_info.size());
2848  AssertIndexRange(cell_batch_index, dof_info[0].cell_active_fe_index.size());
2849  if (dof_info[0].cell_active_fe_index.empty())
2850  return 0;
2851  else
2852  return dof_info[0].cell_active_fe_index[cell_batch_index];
2853 }
2854 
2855 
2856 
2857 template <int dim, typename Number, typename VectorizedArrayType>
2858 inline std::pair<unsigned int, unsigned int>
2860  const unsigned int face_batch_index) const
2861 {
2862  AssertIndexRange(face_batch_index, face_info.faces.size());
2863  if (dof_info[0].cell_active_fe_index.empty())
2864  return std::make_pair(0U, 0U);
2865 
2866  std::pair<unsigned int, unsigned int> result = std::make_pair(0U, 0U);
2867  for (unsigned int v = 0;
2868  v < VectorizedArrayType::size() &&
2869  face_info.faces[face_batch_index].cells_interior[v] !=
2871  ++v)
2872  result.first = std::max(
2873  result.first,
2874  dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
2875  .cells_interior[v] /
2876  VectorizedArrayType::size()]);
2877  if (face_info.faces[face_batch_index].cells_exterior[0] !=
2879  for (unsigned int v = 0;
2880  v < VectorizedArrayType::size() &&
2881  face_info.faces[face_batch_index].cells_exterior[v] !=
2883  ++v)
2884  result.second = std::max(
2885  result.second,
2886  dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
2887  .cells_exterior[v] /
2888  VectorizedArrayType::size()]);
2889  else
2890  result.second = numbers::invalid_unsigned_int;
2891  return result;
2892 }
2893 
2894 
2895 
2896 template <int dim, typename Number, typename VectorizedArrayType>
2897 inline bool
2899 {
2900  return indices_are_initialized;
2901 }
2902 
2903 
2904 
2905 template <int dim, typename Number, typename VectorizedArrayType>
2906 inline bool
2908 {
2909  return mapping_is_initialized;
2910 }
2911 
2912 
2913 template <int dim, typename Number, typename VectorizedArrayType>
2914 inline unsigned int
2916 {
2917  return mg_level;
2918 }
2919 
2920 
2921 
2922 template <int dim, typename Number, typename VectorizedArrayType>
2925 {
2926  using list_type =
2927  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
2928  list_type &data = scratch_pad.get();
2929  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
2930  if (it->first == false)
2931  {
2932  it->first = true;
2933  return &it->second;
2934  }
2935  data.emplace_front(true, AlignedVector<VectorizedArrayType>());
2936  return &data.front().second;
2937 }
2938 
2939 
2940 
2941 template <int dim, typename Number, typename VectorizedArrayType>
2942 void
2944  const AlignedVector<VectorizedArrayType> *scratch) const
2945 {
2946  using list_type =
2947  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
2948  list_type &data = scratch_pad.get();
2949  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
2950  if (&it->second == scratch)
2951  {
2952  Assert(it->first == true, ExcInternalError());
2953  it->first = false;
2954  return;
2955  }
2956  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2957 }
2958 
2959 
2960 
2961 template <int dim, typename Number, typename VectorizedArrayType>
2965 {
2966  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
2968  it != scratch_pad_non_threadsafe.end();
2969  ++it)
2970  if (it->first == false)
2971  {
2972  it->first = true;
2973  return &it->second;
2974  }
2975  scratch_pad_non_threadsafe.push_front(
2976  std::make_pair(true, AlignedVector<Number>()));
2977  return &scratch_pad_non_threadsafe.front().second;
2978 }
2979 
2980 
2981 
2982 template <int dim, typename Number, typename VectorizedArrayType>
2983 void
2986  const AlignedVector<Number> *scratch) const
2987 {
2988  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
2990  it != scratch_pad_non_threadsafe.end();
2991  ++it)
2992  if (&it->second == scratch)
2993  {
2994  Assert(it->first == true, ExcInternalError());
2995  it->first = false;
2996  return;
2997  }
2998  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2999 }
3000 
3001 
3002 
3003 // ------------------------------ reinit functions ---------------------------
3004 
3005 namespace internal
3006 {
3007  namespace MatrixFreeImplementation
3008  {
3009  template <int dim, int spacedim>
3010  inline std::vector<IndexSet>
3011  extract_locally_owned_index_sets(
3012  const std::vector<const ::DoFHandler<dim, spacedim> *> &dofh,
3013  const unsigned int level)
3014  {
3015  std::vector<IndexSet> locally_owned_set;
3016  locally_owned_set.reserve(dofh.size());
3017  for (unsigned int j = 0; j < dofh.size(); ++j)
3019  locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
3020  else
3021  locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
3022  return locally_owned_set;
3023  }
3024  } // namespace MatrixFreeImplementation
3025 } // namespace internal
3026 
3027 
3028 
3029 template <int dim, typename Number, typename VectorizedArrayType>
3030 template <typename QuadratureType, typename number2>
3031 void
3033  const DoFHandler<dim> & dof_handler,
3034  const AffineConstraints<number2> &constraints_in,
3035  const QuadratureType & quad,
3037  &additional_data)
3038 {
3039  std::vector<const DoFHandler<dim, dim> *> dof_handlers;
3040  std::vector<const AffineConstraints<number2> *> constraints;
3041  std::vector<QuadratureType> quads;
3042 
3043  dof_handlers.push_back(&dof_handler);
3044  constraints.push_back(&constraints_in);
3045  quads.push_back(quad);
3046 
3047  std::vector<IndexSet> locally_owned_sets =
3048  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3049  dof_handlers, additional_data.mg_level);
3050 
3051  std::vector<hp::QCollection<dim>> quad_hp;
3052  quad_hp.emplace_back(quad);
3053 
3054  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(
3056  dof_handlers,
3057  constraints,
3058  locally_owned_sets,
3059  quad_hp,
3060  additional_data);
3061 }
3062 
3063 
3064 
3065 template <int dim, typename Number, typename VectorizedArrayType>
3066 template <typename QuadratureType, typename number2, typename MappingType>
3067 void
3069  const MappingType & mapping,
3070  const DoFHandler<dim> & dof_handler,
3071  const AffineConstraints<number2> &constraints_in,
3072  const QuadratureType & quad,
3074  &additional_data)
3075 {
3076  std::vector<const DoFHandler<dim, dim> *> dof_handlers;
3077  std::vector<const AffineConstraints<number2> *> constraints;
3078 
3079  dof_handlers.push_back(&dof_handler);
3080  constraints.push_back(&constraints_in);
3081 
3082  std::vector<IndexSet> locally_owned_sets =
3083  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3084  dof_handlers, additional_data.mg_level);
3085 
3086  std::vector<hp::QCollection<dim>> quad_hp;
3087  quad_hp.emplace_back(quad);
3088 
3089  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3090  dof_handlers,
3091  constraints,
3092  locally_owned_sets,
3093  quad_hp,
3094  additional_data);
3095 }
3096 
3097 
3098 
3099 template <int dim, typename Number, typename VectorizedArrayType>
3100 template <typename QuadratureType, typename number2>
3101 void
3103  const std::vector<const DoFHandler<dim> *> & dof_handler,
3104  const std::vector<const AffineConstraints<number2> *> &constraint,
3105  const std::vector<QuadratureType> & quad,
3107  &additional_data)
3108 {
3109  std::vector<IndexSet> locally_owned_set =
3110  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3111  dof_handler, additional_data.mg_level);
3112  std::vector<hp::QCollection<dim>> quad_hp;
3113  for (unsigned int q = 0; q < quad.size(); ++q)
3114  quad_hp.emplace_back(quad[q]);
3115 
3116  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(
3118  dof_handler,
3119  constraint,
3120  locally_owned_set,
3121  quad_hp,
3122  additional_data);
3123 }
3124 
3125 
3126 
3127 template <int dim, typename Number, typename VectorizedArrayType>
3128 template <typename QuadratureType,
3129  typename number2,
3130  typename DoFHandlerType,
3131  typename MappingType>
3132 void
3134  const MappingType & mapping,
3135  const std::vector<const DoFHandlerType *> & dof_handler,
3136  const std::vector<const AffineConstraints<number2> *> &constraint,
3137  const std::vector<QuadratureType> & quad,
3138  const AdditionalData & additional_data)
3139 {
3140  static_assert(dim == DoFHandlerType::dimension,
3141  "Dimension dim not equal to DoFHandlerType::dimension.");
3142 
3143  std::vector<const DoFHandler<dim> *> dof_handlers;
3144 
3145  for (const auto dh : dof_handler)
3146  dof_handlers.push_back(dh);
3147 
3148  this->reinit(mapping, dof_handlers, constraint, quad, additional_data);
3149 }
3150 
3151 
3152 
3153 template <int dim, typename Number, typename VectorizedArrayType>
3154 template <typename QuadratureType, typename number2>
3155 void
3157  const std::vector<const DoFHandler<dim> *> & dof_handler,
3158  const std::vector<const AffineConstraints<number2> *> &constraint,
3159  const QuadratureType & quad,
3161  &additional_data)
3162 {
3163  std::vector<IndexSet> locally_owned_set =
3164  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3165  dof_handler, additional_data.mg_level);
3166  std::vector<hp::QCollection<dim>> quad_hp;
3167  quad_hp.emplace_back(quad);
3168 
3169  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(
3171  dof_handler,
3172  constraint,
3173  locally_owned_set,
3174  quad_hp,
3175  additional_data);
3176 }
3177 
3178 
3179 
3180 template <int dim, typename Number, typename VectorizedArrayType>
3181 template <typename QuadratureType, typename number2, typename DoFHandlerType>
3182 void
3184  const std::vector<const DoFHandlerType *> & dof_handler,
3185  const std::vector<const AffineConstraints<number2> *> &constraint,
3186  const std::vector<QuadratureType> & quad,
3187  const AdditionalData & additional_data)
3188 {
3189  static_assert(dim == DoFHandlerType::dimension,
3190  "Dimension dim not equal to DoFHandlerType::dimension.");
3191 
3192  std::vector<const DoFHandler<dim> *> dof_handlers;
3193 
3194  for (const auto dh : dof_handler)
3195  dof_handlers.push_back(dh);
3196 
3197  this->reinit(dof_handlers, constraint, quad, additional_data);
3198 }
3199 
3200 
3201 
3202 template <int dim, typename Number, typename VectorizedArrayType>
3203 template <typename QuadratureType, typename number2, typename MappingType>
3204 void
3206  const MappingType & mapping,
3207  const std::vector<const DoFHandler<dim> *> & dof_handler,
3208  const std::vector<const AffineConstraints<number2> *> &constraint,
3209  const QuadratureType & quad,
3211  &additional_data)
3212 {
3213  std::vector<IndexSet> locally_owned_set =
3214  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3215  dof_handler, additional_data.mg_level);
3216  std::vector<hp::QCollection<dim>> quad_hp;
3217  quad_hp.emplace_back(quad);
3218 
3219  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3220  dof_handler,
3221  constraint,
3222  locally_owned_set,
3223  quad_hp,
3224  additional_data);
3225 }
3226 
3227 
3228 
3229 template <int dim, typename Number, typename VectorizedArrayType>
3230 template <typename QuadratureType, typename number2, typename MappingType>
3231 void
3233  const MappingType & mapping,
3234  const std::vector<const DoFHandler<dim> *> & dof_handler,
3235  const std::vector<const AffineConstraints<number2> *> &constraint,
3236  const std::vector<QuadratureType> & quad,
3238  &additional_data)
3239 {
3240  std::vector<IndexSet> locally_owned_set =
3241  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3242  dof_handler, additional_data.mg_level);
3243  std::vector<hp::QCollection<dim>> quad_hp;
3244  for (unsigned int q = 0; q < quad.size(); ++q)
3245  quad_hp.emplace_back(quad[q]);
3246 
3247  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3248  dof_handler,
3249  constraint,
3250  locally_owned_set,
3251  quad_hp,
3252  additional_data);
3253 }
3254 
3255 
3256 
3257 template <int dim, typename Number, typename VectorizedArrayType>
3258 template <typename QuadratureType,
3259  typename number2,
3260  typename DoFHandlerType,
3261  typename MappingType>
3262 void
3264  const MappingType & mapping,
3265  const std::vector<const DoFHandlerType *> & dof_handler,
3266  const std::vector<const AffineConstraints<number2> *> &constraint,
3267  const QuadratureType & quad,
3268  const AdditionalData & additional_data)
3269 {
3270  static_assert(dim == DoFHandlerType::dimension,
3271  "Dimension dim not equal to DoFHandlerType::dimension.");
3272 
3273  std::vector<const DoFHandler<dim> *> dof_handlers;
3274 
3275  for (const auto dh : dof_handler)
3276  dof_handlers.push_back(dh);
3277 
3278  this->reinit(mapping, dof_handlers, constraint, quad, additional_data);
3279 }
3280 
3281 
3282 
3283 template <int dim, typename Number, typename VectorizedArrayType>
3284 template <typename QuadratureType, typename number2, typename DoFHandlerType>
3285 void
3287  const std::vector<const DoFHandlerType *> & dof_handler,
3288  const std::vector<const AffineConstraints<number2> *> &constraint,
3289  const QuadratureType & quad,
3290  const AdditionalData & additional_data)
3291 {
3292  static_assert(dim == DoFHandlerType::dimension,
3293  "Dimension dim not equal to DoFHandlerType::dimension.");
3294 
3295  std::vector<const DoFHandler<dim> *> dof_handlers;
3296 
3297  for (const auto dh : dof_handler)
3298  dof_handlers.push_back(dh);
3299 
3300  this->reinit(dof_handlers, constraint, quad, additional_data);
3301 }
3302 
3303 
3304 
3305 // ------------------------------ implementation of loops --------------------
3306 
3307 // internal helper functions that define how to call MPI data exchange
3308 // functions: for generic vectors, do nothing at all. For distributed vectors,
3309 // call update_ghost_values_start function and so on. If we have collections
3310 // of vectors, just do the individual functions of the components. In order to
3311 // keep ghost values consistent (whether we are in read or write mode), we
3312 // also reset the values at the end. the whole situation is a bit complicated
3313 // by the fact that we need to treat block vectors differently, which use some
3314 // additional helper functions to select the blocks and template magic.
3315 namespace internal
3316 {
3320  template <int dim, typename Number, typename VectorizedArrayType>
3321  struct VectorDataExchange
3322  {
3323  // A shift for the MPI messages to reduce the risk for accidental
3324  // interaction with other open communications that a user program might
3325  // set up (parallel vectors support unfinished communication). We let
3326  // the other vectors use the first 20 assigned numbers and start the
3327  // matrix-free communication.
3328  static constexpr unsigned int channel_shift = 20;
3329 
3330 
3331 
3336  VectorDataExchange(
3337  const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3338  const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3339  DataAccessOnFaces vector_face_access,
3340  const unsigned int n_components)
3341  : matrix_free(matrix_free)
3342  , vector_face_access(
3343  matrix_free.get_task_info().face_partition_data.empty() ?
3344  ::MatrixFree<dim, Number, VectorizedArrayType>::
3345  DataAccessOnFaces::unspecified :
3346  vector_face_access)
3347  , ghosts_were_set(false)
3348 # ifdef DEAL_II_WITH_MPI
3349  , tmp_data(n_components)
3350  , requests(n_components)
3351 # endif
3352  {
3353  (void)n_components;
3354  if (this->vector_face_access !=
3356  DataAccessOnFaces::unspecified)
3357  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3359  matrix_free.get_dof_info(c).vector_exchanger_face_variants.size(),
3360  5);
3361  }
3362 
3363 
3364 
3368  ~VectorDataExchange() // NOLINT
3369  {
3370 # ifdef DEAL_II_WITH_MPI
3371  for (unsigned int i = 0; i < tmp_data.size(); ++i)
3372  if (tmp_data[i] != nullptr)
3373  matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3374 # endif
3375  }
3376 
3377 
3378 
3383  template <typename VectorType>
3384  unsigned int
3385  find_vector_in_mf(const VectorType &vec,
3386  const bool check_global_compatibility = true) const
3387  {
3388  // case 1: vector was set up with MatrixFree::initialize_dof_vector()
3389  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3390  if (vec.get_partitioner().get() ==
3391  matrix_free.get_dof_info(c).vector_partitioner.get())
3392  return c;
3393 
3394  // case 2: user provided own partitioner (compatibility mode)
3395  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3396  if (check_global_compatibility ?
3397  vec.get_partitioner()->is_globally_compatible(
3398  *matrix_free.get_dof_info(c).vector_partitioner) :
3399  vec.get_partitioner()->is_compatible(
3400  *matrix_free.get_dof_info(c).vector_partitioner))
3401  return c;
3402 
3403  Assert(false,
3404  ExcNotImplemented("Could not find partitioner that fits vector"));
3405 
3407  }
3408 
3409 
3410 
3416  get_partitioner(const unsigned int mf_component) const
3417  {
3418  AssertDimension(matrix_free.get_dof_info(mf_component)
3419  .vector_exchanger_face_variants.size(),
3420  5);
3421  if (vector_face_access ==
3424  return *matrix_free.get_dof_info(mf_component)
3425  .vector_exchanger_face_variants[0];
3426  else if (vector_face_access ==
3429  return *matrix_free.get_dof_info(mf_component)
3430  .vector_exchanger_face_variants[1];
3431  else if (vector_face_access ==
3434  return *matrix_free.get_dof_info(mf_component)
3435  .vector_exchanger_face_variants[2];
3436  else if (vector_face_access ==
3438  DataAccessOnFaces::values_all_faces)
3439  return *matrix_free.get_dof_info(mf_component)
3440  .vector_exchanger_face_variants[3];
3441  else if (vector_face_access ==
3443  DataAccessOnFaces::gradients_all_faces)
3444  return *matrix_free.get_dof_info(mf_component)
3445  .vector_exchanger_face_variants[4];
3446  else
3447  return *matrix_free.get_dof_info(mf_component).vector_exchanger.get();
3448  }
3449 
3450 
3451 
3455  template <typename VectorType,
3456  typename std::enable_if<is_not_parallel_vector<VectorType>,
3457  VectorType>::type * = nullptr>
3458  void
3459  update_ghost_values_start(const unsigned int /*component_in_block_vector*/,
3460  const VectorType & /*vec*/)
3461  {}
3462 
3463 
3468  template <
3469  typename VectorType,
3470  typename std::enable_if<!has_update_ghost_values_start<VectorType> &&
3471  !is_not_parallel_vector<VectorType>,
3472  VectorType>::type * = nullptr>
3473  void
3474  update_ghost_values_start(const unsigned int component_in_block_vector,
3475  const VectorType & vec)
3476  {
3477  (void)component_in_block_vector;
3478  bool ghosts_set = vec.has_ghost_elements();
3479 
3480  Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3481  ghosts_set == false,
3482  ExcNotImplemented());
3483 
3484  if (ghosts_set)
3485  ghosts_were_set = true;
3486 
3487  vec.update_ghost_values();
3488  }
3489 
3490 
3491 
3497  template <
3498  typename VectorType,
3499  typename std::enable_if<has_update_ghost_values_start<VectorType> &&
3500  !has_exchange_on_subset<VectorType>,
3501  VectorType>::type * = nullptr>
3502  void
3503  update_ghost_values_start(const unsigned int component_in_block_vector,
3504  const VectorType & vec)
3505  {
3506  (void)component_in_block_vector;
3507  bool ghosts_set = vec.has_ghost_elements();
3508 
3509  Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3510  ghosts_set == false,
3511  ExcNotImplemented());
3512 
3513  if (ghosts_set)
3514  ghosts_were_set = true;
3515 
3516  vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3517  }
3518 
3519 
3520 
3527  template <
3528  typename VectorType,
3529  typename std::enable_if<has_update_ghost_values_start<VectorType> &&
3530  has_exchange_on_subset<VectorType>,
3531  VectorType>::type * = nullptr>
3532  void
3533  update_ghost_values_start(const unsigned int component_in_block_vector,
3534  const VectorType & vec)
3535  {
3536  static_assert(
3537  std::is_same<Number, typename VectorType::value_type>::value,
3538  "Type mismatch between VectorType and VectorDataExchange");
3539  (void)component_in_block_vector;
3540  bool ghosts_set = vec.has_ghost_elements();
3541 
3542  Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3543  ghosts_set == false,
3544  ExcNotImplemented());
3545 
3546  if (ghosts_set)
3547  ghosts_were_set = true;
3548 
3549  if (vec.size() != 0)
3550  {
3551 # ifdef DEAL_II_WITH_MPI
3552  const unsigned int mf_component = find_vector_in_mf(vec);
3553 
3554  const auto &part = get_partitioner(mf_component);
3555 
3556  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3557  part.n_import_sm_procs() == 0)
3558  return;
3559 
3560  tmp_data[component_in_block_vector] =
3561  matrix_free.acquire_scratch_data_non_threadsafe();
3562  tmp_data[component_in_block_vector]->resize_fast(
3563  part.n_import_indices());
3564  AssertDimension(requests.size(), tmp_data.size());
3565 
3566  part.export_to_ghosted_array_start(
3567  component_in_block_vector * 2 + channel_shift,
3568  ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3569  vec.shared_vector_data(),
3570  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3571  part.locally_owned_size(),
3572  matrix_free.get_dof_info(mf_component)
3573  .vector_partitioner->n_ghost_indices()),
3574  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3575  part.n_import_indices()),
3576  this->requests[component_in_block_vector]);
3577 # endif
3578  }
3579  }
3580 
3581 
3582 
3587  template <
3588  typename VectorType,
3589  typename std::enable_if<!has_update_ghost_values_start<VectorType>,
3590  VectorType>::type * = nullptr>
3591  void
3592  update_ghost_values_finish(const unsigned int /*component_in_block_vector*/,
3593  const VectorType & /*vec*/)
3594  {}
3595 
3596 
3597 
3603  template <
3604  typename VectorType,
3605  typename std::enable_if<has_update_ghost_values_start<VectorType> &&
3606  !has_exchange_on_subset<VectorType>,
3607  VectorType>::type * = nullptr>
3608  void
3609  update_ghost_values_finish(const unsigned int component_in_block_vector,
3610  const VectorType & vec)
3611  {
3612  (void)component_in_block_vector;
3613  vec.update_ghost_values_finish();
3614  }
3615 
3616 
3617 
3624  template <
3625  typename VectorType,
3626  typename std::enable_if<has_update_ghost_values_start<VectorType> &&
3627  has_exchange_on_subset<VectorType>,
3628  VectorType>::type * = nullptr>
3629  void
3630  update_ghost_values_finish(const unsigned int component_in_block_vector,
3631  const VectorType & vec)
3632  {
3633  static_assert(
3634  std::is_same<Number, typename VectorType::value_type>::value,
3635  "Type mismatch between VectorType and VectorDataExchange");
3636  (void)component_in_block_vector;
3637 
3638  if (vec.size() != 0)
3639  {
3640 # ifdef DEAL_II_WITH_MPI
3641  AssertIndexRange(component_in_block_vector, tmp_data.size());
3642  AssertDimension(requests.size(), tmp_data.size());
3643 
3644  const unsigned int mf_component = find_vector_in_mf(vec);
3645 
3646  const auto &part = get_partitioner(mf_component);
3647 
3648  if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3649  part.n_import_sm_procs() != 0)
3650  {
3651  part.export_to_ghosted_array_finish(
3652  ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3653  vec.shared_vector_data(),
3654  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3655  part.locally_owned_size(),
3656  matrix_free.get_dof_info(mf_component)
3657  .vector_partitioner->n_ghost_indices()),
3658  this->requests[component_in_block_vector]);
3659 
3660  matrix_free.release_scratch_data_non_threadsafe(
3661  tmp_data[component_in_block_vector]);
3662  tmp_data[component_in_block_vector] = nullptr;
3663  }
3664 # endif
3665  }
3666  // let vector know that ghosts are being updated and we can read from
3667  // them
3668  vec.set_ghost_state(true);
3669  }
3670 
3671 
3672 
3676  template <typename VectorType,
3677  typename std::enable_if<is_not_parallel_vector<VectorType>,
3678  VectorType>::type * = nullptr>
3679  void
3680  compress_start(const unsigned int /*component_in_block_vector*/,
3681  VectorType & /*vec*/)
3682  {}
3683 
3684 
3685 
3690  template <typename VectorType,
3691  typename std::enable_if<!has_compress_start<VectorType> &&
3692  !is_not_parallel_vector<VectorType>,
3693  VectorType>::type * = nullptr>
3694  void
3695  compress_start(const unsigned int component_in_block_vector,
3696  VectorType & vec)
3697  {
3698  (void)component_in_block_vector;
3699  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3700  vec.compress(::VectorOperation::add);
3701  }
3702 
3703 
3704 
3710  template <typename VectorType,
3711  typename std::enable_if<has_compress_start<VectorType> &&
3712  !has_exchange_on_subset<VectorType>,
3713  VectorType>::type * = nullptr>
3714  void
3715  compress_start(const unsigned int component_in_block_vector,
3716  VectorType & vec)
3717  {
3718  (void)component_in_block_vector;
3719  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3720  vec.compress_start(component_in_block_vector + channel_shift);
3721  }
3722 
3723 
3724 
3731  template <typename VectorType,
3732  typename std::enable_if<has_compress_start<VectorType> &&
3733  has_exchange_on_subset<VectorType>,
3734  VectorType>::type * = nullptr>
3735  void
3736  compress_start(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  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3744 
3745  if (vec.size() != 0)
3746  {
3747 # ifdef DEAL_II_WITH_MPI
3748  const unsigned int mf_component = find_vector_in_mf(vec);
3749 
3750  const auto &part = get_partitioner(mf_component);
3751 
3752  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3753  part.n_import_sm_procs() == 0)
3754  return;
3755 
3756  tmp_data[component_in_block_vector] =
3757  matrix_free.acquire_scratch_data_non_threadsafe();
3758  tmp_data[component_in_block_vector]->resize_fast(
3759  part.n_import_indices());
3760  AssertDimension(requests.size(), tmp_data.size());
3761 
3762  part.import_from_ghosted_array_start(
3764  component_in_block_vector * 2 + channel_shift,
3765  ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3766  vec.shared_vector_data(),
3767  ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3768  matrix_free.get_dof_info(mf_component)
3769  .vector_partitioner->n_ghost_indices()),
3770  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3771  part.n_import_indices()),
3772  this->requests[component_in_block_vector]);
3773 # endif
3774  }
3775  }
3776 
3777 
3778 
3783  template <typename VectorType,
3784  typename std::enable_if<!has_compress_start<VectorType>,
3785  VectorType>::type * = nullptr>
3786  void
3787  compress_finish(const unsigned int /*component_in_block_vector*/,
3788  VectorType & /*vec*/)
3789  {}
3790 
3791 
3792 
3798  template <typename VectorType,
3799  typename std::enable_if<has_compress_start<VectorType> &&
3800  !has_exchange_on_subset<VectorType>,
3801  VectorType>::type * = nullptr>
3802  void
3803  compress_finish(const unsigned int component_in_block_vector,
3804  VectorType & vec)
3805  {
3806  (void)component_in_block_vector;
3807  vec.compress_finish(::VectorOperation::add);
3808  }
3809 
3810 
3811 
3818  template <typename VectorType,
3819  typename std::enable_if<has_compress_start<VectorType> &&
3820  has_exchange_on_subset<VectorType>,
3821  VectorType>::type * = nullptr>
3822  void
3823  compress_finish(const unsigned int component_in_block_vector,
3824  VectorType & vec)
3825  {
3826  static_assert(
3827  std::is_same<Number, typename VectorType::value_type>::value,
3828  "Type mismatch between VectorType and VectorDataExchange");
3829  (void)component_in_block_vector;
3830  if (vec.size() != 0)
3831  {
3832 # ifdef DEAL_II_WITH_MPI
3833  AssertIndexRange(component_in_block_vector, tmp_data.size());
3834  AssertDimension(requests.size(), tmp_data.size());
3835 
3836  const unsigned int mf_component = find_vector_in_mf(vec);
3837 
3838  const auto &part = get_partitioner(mf_component);
3839 
3840  if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3841  part.n_import_sm_procs() != 0)
3842  {
3843  part.import_from_ghosted_array_finish(
3845  ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3846  vec.shared_vector_data(),
3847  ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3848  matrix_free.get_dof_info(mf_component)
3849  .vector_partitioner->n_ghost_indices()),
3851  tmp_data[component_in_block_vector]->begin(),
3852  part.n_import_indices()),
3853  this->requests[component_in_block_vector]);
3854 
3855  matrix_free.release_scratch_data_non_threadsafe(
3856  tmp_data[component_in_block_vector]);
3857  tmp_data[component_in_block_vector] = nullptr;
3858  }
3859 
3861  {
3862  const int ierr =
3863  MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3864  AssertThrowMPI(ierr);
3865  }
3866 # endif
3867  }
3868  }
3869 
3870 
3871 
3875  template <typename VectorType,
3876  typename std::enable_if<is_not_parallel_vector<VectorType>,
3877  VectorType>::type * = nullptr>
3878  void
3879  reset_ghost_values(const VectorType & /*vec*/) const
3880  {}
3881 
3882 
3883 
3888  template <typename VectorType,
3889  typename std::enable_if<!has_exchange_on_subset<VectorType> &&
3890  !is_not_parallel_vector<VectorType>,
3891  VectorType>::type * = nullptr>
3892  void
3893  reset_ghost_values(const VectorType &vec) const
3894  {
3895  if (ghosts_were_set == true)
3896  return;
3897 
3898  vec.zero_out_ghost_values();
3899  }
3900 
3901 
3902 
3908  template <typename VectorType,
3909  typename std::enable_if<has_exchange_on_subset<VectorType>,
3910  VectorType>::type * = nullptr>
3911  void
3912  reset_ghost_values(const VectorType &vec) const
3913  {
3914  static_assert(
3915  std::is_same<Number, typename VectorType::value_type>::value,
3916  "Type mismatch between VectorType and VectorDataExchange");
3917  if (ghosts_were_set == true)
3918  return;
3919 
3920  if (vec.size() != 0)
3921  {
3922 # ifdef DEAL_II_WITH_MPI
3923  AssertDimension(requests.size(), tmp_data.size());
3924 
3925  const unsigned int mf_component = find_vector_in_mf(vec);
3926 
3927  const auto &part = get_partitioner(mf_component);
3928 
3929  if (part.n_ghost_indices() > 0)
3930  {
3931  part.reset_ghost_values(ArrayView<Number>(
3933  .begin() +
3934  part.locally_owned_size(),
3935  matrix_free.get_dof_info(mf_component)
3936  .vector_partitioner->n_ghost_indices()));
3937  }
3938 
3939 # endif
3940  }
3941  // let vector know that it's not ghosted anymore
3942  vec.set_ghost_state(false);
3943  }
3944 
3945 
3946 
3952  template <typename VectorType,
3953  typename std::enable_if<has_exchange_on_subset<VectorType>,
3954  VectorType>::type * = nullptr>
3955  void
3956  zero_vector_region(const unsigned int range_index, VectorType &vec) const
3957  {
3958  static_assert(
3959  std::is_same<Number, typename VectorType::value_type>::value,
3960  "Type mismatch between VectorType and VectorDataExchange");
3961  if (range_index == numbers::invalid_unsigned_int)
3962  vec = Number();
3963  else
3964  {
3965  const unsigned int mf_component = find_vector_in_mf(vec, false);
3966  const internal::MatrixFreeFunctions::DoFInfo &dof_info =
3967  matrix_free.get_dof_info(mf_component);
3968  Assert(dof_info.vector_zero_range_list_index.empty() == false,
3969  ExcNotInitialized());
3970 
3971  Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3972  ExcInternalError());
3973  AssertIndexRange(range_index,
3974  dof_info.vector_zero_range_list_index.size() - 1);
3975  for (unsigned int id =
3976  dof_info.vector_zero_range_list_index[range_index];
3977  id != dof_info.vector_zero_range_list_index[range_index + 1];
3978  ++id)
3979  std::memset(vec.begin() + dof_info.vector_zero_range_list[id].first,
3980  0,
3981  (dof_info.vector_zero_range_list[id].second -
3982  dof_info.vector_zero_range_list[id].first) *
3983  sizeof(Number));
3984  }
3985  }
3986 
3987 
3988 
3994  template <typename VectorType,
3995  typename std::enable_if<!has_exchange_on_subset<VectorType>,
3996  VectorType>::type * = nullptr,
3997  typename VectorType::value_type * = nullptr>
3998  void
3999  zero_vector_region(const unsigned int range_index, VectorType &vec) const
4000  {
4001  if (range_index == numbers::invalid_unsigned_int || range_index == 0)
4002  vec = typename VectorType::value_type();
4003  }
4004 
4005 
4006 
4011  void
4012  zero_vector_region(const unsigned int, ...) const
4013  {
4014  Assert(false,
4015  ExcNotImplemented("Zeroing is only implemented for vector types "
4016  "which provide VectorType::value_type"));
4017  }
4018 
4019 
4020 
4021  const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
4022  const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
4023  DataAccessOnFaces vector_face_access;
4024  bool ghosts_were_set;
4025 # ifdef DEAL_II_WITH_MPI
4026  std::vector<AlignedVector<Number> *> tmp_data;
4027  std::vector<std::vector<MPI_Request>> requests;
4028 # endif
4029  }; // VectorDataExchange
4030 
4031  template <typename VectorStruct>
4032  unsigned int
4033  n_components(const VectorStruct &vec);
4034 
4035  template <typename VectorStruct>
4036  unsigned int
4037  n_components_block(const VectorStruct &vec,
4038  std::integral_constant<bool, true>)
4039  {
4040  unsigned int components = 0;
4041  for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
4042  components += n_components(vec.block(bl));
4043  return components;
4044  }
4045 
4046  template <typename VectorStruct>
4047  unsigned int
4048  n_components_block(const VectorStruct &, std::integral_constant<bool, false>)
4049  {
4050  return 1;
4051  }
4052 
4053  template <typename VectorStruct>
4054  unsigned int
4055  n_components(const VectorStruct &vec)
4056  {
4057  return n_components_block(
4058  vec, std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
4059  }
4060 
4061  template <typename VectorStruct>
4062  inline unsigned int
4063  n_components(const std::vector<VectorStruct> &vec)
4064  {
4065  unsigned int components = 0;
4066  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4067  components += n_components_block(
4068  vec[comp],
4069  std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
4070  return components;
4071  }
4072 
4073  template <typename VectorStruct>
4074  inline unsigned int
4075  n_components(const std::vector<VectorStruct *> &vec)
4076  {
4077  unsigned int components = 0;
4078  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4079  components += n_components_block(
4080  *vec[comp],
4081  std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
4082  return components;
4083  }
4084 
4085 
4086 
4087  // A helper function to identify block vectors with many components where we
4088  // should not try to overlap computations and communication because there
4089  // would be too many outstanding communication requests.
4090 
4091  // default value for vectors that do not have communication_block_size
4092  template <typename VectorStruct,
4093  typename std::enable_if<!has_communication_block_size<VectorStruct>,
4094  VectorStruct>::type * = nullptr>
4095  constexpr unsigned int
4096  get_communication_block_size(const VectorStruct &)
4097  {
4099  }
4100 
4101 
4102 
4103  template <typename VectorStruct,
4104  typename std::enable_if<has_communication_block_size<VectorStruct>,
4105  VectorStruct>::type * = nullptr>
4106  constexpr unsigned int
4107  get_communication_block_size(const VectorStruct &)
4108  {
4109  return VectorStruct::communication_block_size;
4110  }
4111 
4112 
4113 
4114  // --------------------------------------------------------------------------
4115  // below we have wrappers to distinguish between block and non-block vectors.
4116  // --------------------------------------------------------------------------
4117 
4118  //
4119  // update_ghost_values_start
4120  //
4121 
4122  // update_ghost_values for block vectors
4123  template <int dim,
4124  typename VectorStruct,
4125  typename Number,
4126  typename VectorizedArrayType,
4127  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4128  VectorStruct>::type * = nullptr>
4129  void
4130  update_ghost_values_start(
4131  const VectorStruct & vec,
4132  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4133  const unsigned int channel = 0)
4134  {
4135  if (get_communication_block_size(vec) < vec.n_blocks())
4136  {
4137  // don't forget to set ghosts_were_set, that otherwise happens
4138  // inside VectorDataExchange::update_ghost_values_start()
4139  exchanger.ghosts_were_set = vec.has_ghost_elements();
4140  vec.update_ghost_values();
4141  }
4142  else
4143  {
4144  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4145  update_ghost_values_start(vec.block(i), exchanger, channel + i);
4146  }
4147  }
4148 
4149 
4150 
4151  // update_ghost_values for non-block vectors
4152  template <int dim,
4153  typename VectorStruct,
4154  typename Number,
4155  typename VectorizedArrayType,
4156  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4157  VectorStruct>::type * = nullptr>
4158  void
4159  update_ghost_values_start(
4160  const VectorStruct & vec,
4161  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4162  const unsigned int channel = 0)
4163  {
4164  exchanger.update_ghost_values_start(channel, vec);
4165  }
4166 
4167 
4168 
4169  // update_ghost_values_start() for vector of vectors
4170  template <int dim,
4171  typename VectorStruct,
4172  typename Number,
4173  typename VectorizedArrayType>
4174  inline void
4175  update_ghost_values_start(
4176  const std::vector<VectorStruct> & vec,
4177  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4178  {
4179  unsigned int component_index = 0;
4180  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4181  {
4182  update_ghost_values_start(vec[comp], exchanger, component_index);
4183  component_index += n_components(vec[comp]);
4184  }
4185  }
4186 
4187 
4188 
4189  // update_ghost_values_start() for vector of pointers to vectors
4190  template <int dim,
4191  typename VectorStruct,
4192  typename Number,
4193  typename VectorizedArrayType>
4194  inline void
4195  update_ghost_values_start(
4196  const std::vector<VectorStruct *> & vec,
4197  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4198  {
4199  unsigned int component_index = 0;
4200  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4201  {
4202  update_ghost_values_start(*vec[comp], exchanger, component_index);
4203  component_index += n_components(*vec[comp]);
4204  }
4205  }
4206 
4207 
4208 
4209  //
4210  // update_ghost_values_finish
4211  //
4212 
4213  // for block vectors
4214  template <int dim,
4215  typename VectorStruct,
4216  typename Number,
4217  typename VectorizedArrayType,
4218  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4219  VectorStruct>::type * = nullptr>
4220  void
4221  update_ghost_values_finish(
4222  const VectorStruct & vec,
4223  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4224  const unsigned int channel = 0)
4225  {
4226  if (get_communication_block_size(vec) < vec.n_blocks())
4227  {
4228  // do nothing, everything has already been completed in the _start()
4229  // call
4230  }
4231  else
4232  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4233  update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4234  }
4235 
4236 
4237 
4238  // for non-block vectors
4239  template <int dim,
4240  typename VectorStruct,
4241  typename Number,
4242  typename VectorizedArrayType,
4243  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4244  VectorStruct>::type * = nullptr>
4245  void
4246  update_ghost_values_finish(
4247  const VectorStruct & vec,
4248  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4249  const unsigned int channel = 0)
4250  {
4251  exchanger.update_ghost_values_finish(channel, vec);
4252  }
4253 
4254 
4255 
4256  // for vector of vectors
4257  template <int dim,
4258  typename VectorStruct,
4259  typename Number,
4260  typename VectorizedArrayType>
4261  inline void
4262  update_ghost_values_finish(
4263  const std::vector<VectorStruct> & vec,
4264  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4265  {
4266  unsigned int component_index = 0;
4267  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4268  {
4269  update_ghost_values_finish(vec[comp], exchanger, component_index);
4270  component_index += n_components(vec[comp]);
4271  }
4272  }
4273 
4274 
4275 
4276  // for vector of pointers to vectors
4277  template <int dim,
4278  typename VectorStruct,
4279  typename Number,
4280  typename VectorizedArrayType>
4281  inline void
4282  update_ghost_values_finish(
4283  const 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  update_ghost_values_finish(*vec[comp], exchanger, component_index);
4290  component_index += n_components(*vec[comp]);
4291  }
4292  }
4293 
4294 
4295 
4296  //
4297  // compress_start
4298  //
4299 
4300  // for block vectors
4301  template <int dim,
4302  typename VectorStruct,
4303  typename Number,
4304  typename VectorizedArrayType,
4305  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4306  VectorStruct>::type * = nullptr>
4307  inline void
4308  compress_start(
4309  VectorStruct & vec,
4310  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4311  const unsigned int channel = 0)
4312  {
4313  if (get_communication_block_size(vec) < vec.n_blocks())
4314  vec.compress(::VectorOperation::add);
4315  else
4316  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4317  compress_start(vec.block(i), exchanger, channel + i);
4318  }
4319 
4320 
4321 
4322  // for non-block vectors
4323  template <int dim,
4324  typename VectorStruct,
4325  typename Number,
4326  typename VectorizedArrayType,
4327  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4328  VectorStruct>::type * = nullptr>
4329  inline void
4330  compress_start(
4331  VectorStruct & vec,
4332  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4333  const unsigned int channel = 0)
4334  {
4335  exchanger.compress_start(channel, vec);
4336  }
4337 
4338 
4339 
4340  // for std::vector of vectors
4341  template <int dim,
4342  typename VectorStruct,
4343  typename Number,
4344  typename VectorizedArrayType>
4345  inline void
4346  compress_start(
4347  std::vector<VectorStruct> & vec,
4348  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4349  {
4350  unsigned int component_index = 0;
4351  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4352  {
4353  compress_start(vec[comp], exchanger, component_index);
4354  component_index += n_components(vec[comp]);
4355  }
4356  }
4357 
4358 
4359 
4360  // for std::vector of pointer to vectors
4361  template <int dim,
4362  typename VectorStruct,
4363  typename Number,
4364  typename VectorizedArrayType>
4365  inline void
4366  compress_start(
4367  std::vector<VectorStruct *> & vec,
4368  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4369  {
4370  unsigned int component_index = 0;
4371  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4372  {
4373  compress_start(*vec[comp], exchanger, component_index);
4374  component_index += n_components(*vec[comp]);
4375  }
4376  }
4377 
4378 
4379 
4380  //
4381  // compress_finish
4382  //
4383 
4384  // for block vectors
4385  template <int dim,
4386  typename VectorStruct,
4387  typename Number,
4388  typename VectorizedArrayType,
4389  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4390  VectorStruct>::type * = nullptr>
4391  inline void
4392  compress_finish(
4393  VectorStruct & vec,
4394  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4395  const unsigned int channel = 0)
4396  {
4397  if (get_communication_block_size(vec) < vec.n_blocks())
4398  {
4399  // do nothing, everything has already been completed in the _start()
4400  // call
4401  }
4402  else
4403  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4404  compress_finish(vec.block(i), exchanger, channel + i);
4405  }
4406 
4407 
4408 
4409  // for non-block vectors
4410  template <int dim,
4411  typename VectorStruct,
4412  typename Number,
4413  typename VectorizedArrayType,
4414  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4415  VectorStruct>::type * = nullptr>
4416  inline void
4417  compress_finish(
4418  VectorStruct & vec,
4419  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4420  const unsigned int channel = 0)
4421  {
4422  exchanger.compress_finish(channel, vec);
4423  }
4424 
4425 
4426 
4427  // for std::vector of vectors
4428  template <int dim,
4429  typename VectorStruct,
4430  typename Number,
4431  typename VectorizedArrayType>
4432  inline void
4433  compress_finish(
4434  std::vector<VectorStruct> & vec,
4435  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4436  {
4437  unsigned int component_index = 0;
4438  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4439  {
4440  compress_finish(vec[comp], exchanger, component_index);
4441  component_index += n_components(vec[comp]);
4442  }
4443  }
4444 
4445 
4446 
4447  // for std::vector of pointer to vectors
4448  template <int dim,
4449  typename VectorStruct,
4450  typename Number,
4451  typename VectorizedArrayType>
4452  inline void
4453  compress_finish(
4454  std::vector<VectorStruct *> & vec,
4455  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4456  {
4457  unsigned int component_index = 0;
4458  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4459  {
4460  compress_finish(*vec[comp], exchanger, component_index);
4461  component_index += n_components(*vec[comp]);
4462  }
4463  }
4464 
4465 
4466 
4467  //
4468  // reset_ghost_values:
4469  //
4470  // if the input vector did not have ghosts imported, clear them here again
4471  // in order to avoid subsequent operations e.g. in linear solvers to work
4472  // with ghosts all the time
4473  //
4474 
4475  // for block vectors
4476  template <int dim,
4477  typename VectorStruct,
4478  typename Number,
4479  typename VectorizedArrayType,
4480  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4481  VectorStruct>::type * = nullptr>
4482  inline void
4483  reset_ghost_values(
4484  const VectorStruct & vec,
4485  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4486  {
4487  // return immediately if there is nothing to do.
4488  if (exchanger.ghosts_were_set == true)
4489  return;
4490 
4491  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4492  reset_ghost_values(vec.block(i), exchanger);
4493  }
4494 
4495 
4496 
4497  // for non-block vectors
4498  template <int dim,
4499  typename VectorStruct,
4500  typename Number,
4501  typename VectorizedArrayType,
4502  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4503  VectorStruct>::type * = nullptr>
4504  inline void
4505  reset_ghost_values(
4506  const VectorStruct & vec,
4507  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4508  {
4509  exchanger.reset_ghost_values(vec);
4510  }
4511 
4512 
4513 
4514  // for std::vector of vectors
4515  template <int dim,
4516  typename VectorStruct,
4517  typename Number,
4518  typename VectorizedArrayType>
4519  inline void
4520  reset_ghost_values(
4521  const std::vector<VectorStruct> & vec,
4522  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4523  {
4524  // return immediately if there is nothing to do.
4525  if (exchanger.ghosts_were_set == true)
4526  return;
4527 
4528  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4529  reset_ghost_values(vec[comp], exchanger);
4530  }
4531 
4532 
4533 
4534  // for std::vector of pointer to vectors
4535  template <int dim,
4536  typename VectorStruct,
4537  typename Number,
4538  typename VectorizedArrayType>
4539  inline void
4540  reset_ghost_values(
4541  const std::vector<VectorStruct *> & vec,
4542  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4543  {
4544  // return immediately if there is nothing to do.
4545  if (exchanger.ghosts_were_set == true)
4546  return;
4547 
4548  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4549  reset_ghost_values(*vec[comp], exchanger);
4550  }
4551 
4552 
4553 
4554  //
4555  // zero_vector_region
4556  //
4557 
4558  // for block vectors
4559  template <int dim,
4560  typename VectorStruct,
4561  typename Number,
4562  typename VectorizedArrayType,
4563  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4564  VectorStruct>::type * = nullptr>
4565  inline void
4566  zero_vector_region(
4567  const unsigned int range_index,
4568  VectorStruct & vec,
4569  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4570  {
4571  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4572  exchanger.zero_vector_region(range_index, vec.block(i));
4573  }
4574 
4575 
4576 
4577  // for non-block vectors
4578  template <int dim,
4579  typename VectorStruct,
4580  typename Number,
4581  typename VectorizedArrayType,
4582  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4583  VectorStruct>::type * = nullptr>
4584  inline void
4585  zero_vector_region(
4586  const unsigned int range_index,
4587  VectorStruct & vec,
4588  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4589  {
4590  exchanger.zero_vector_region(range_index, vec);
4591  }
4592 
4593 
4594 
4595  // for std::vector of vectors
4596  template <int dim,
4597  typename VectorStruct,
4598  typename Number,
4599  typename VectorizedArrayType>
4600  inline void
4601  zero_vector_region(
4602  const unsigned int range_index,
4603  std::vector<VectorStruct> & vec,
4604  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4605  {
4606  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4607  zero_vector_region(range_index, vec[comp], exchanger);
4608  }
4609 
4610 
4611 
4612  // for std::vector of pointers to vectors
4613  template <int dim,
4614  typename VectorStruct,
4615  typename Number,
4616  typename VectorizedArrayType>
4617  inline void
4618  zero_vector_region(
4619  const unsigned int range_index,
4620  std::vector<VectorStruct *> & vec,
4621  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4622  {
4623  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4624  zero_vector_region(range_index, *vec[comp], exchanger);
4625  }
4626 
4627 
4628 
4629  // Apply a unit matrix operation to constrained DoFs: Default cases where we
4630  // cannot detect a LinearAlgebra::distributed::Vector, we do not do
4631  // anything, else we apply the constraints as a unit operation
4632  template <typename VectorStruct1, typename VectorStruct2>
4633  inline void
4634  apply_operation_to_constrained_dofs(const std::vector<unsigned int> &,
4635  const VectorStruct1 &,
4636  VectorStruct2 &)
4637  {}
4638 
4639  template <typename Number>
4640  inline void
4641  apply_operation_to_constrained_dofs(
4642  const std::vector<unsigned int> & constrained_dofs,
4645  {
4646  for (const unsigned int i : constrained_dofs)
4647  dst.local_element(i) = src.local_element(i);
4648  }
4649 
4650 
4651  namespace MatrixFreeFunctions
4652  {
4653  // struct to select between a const interface and a non-const interface
4654  // for MFWorker
4655  template <typename, typename, typename, typename, bool>
4656  struct InterfaceSelector
4657  {};
4658 
4659  // Version for constant functions
4660  template <typename MF,
4661  typename InVector,
4662  typename OutVector,
4663  typename Container>
4664  struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4665  {
4666  using function_type = void (Container::*)(
4667  const MF &,
4668  OutVector &,
4669  const InVector &,
4670  const std::pair<unsigned int, unsigned int> &) const;
4671  };
4672 
4673  // Version for non-constant functions
4674  template <typename MF,
4675  typename InVector,
4676  typename OutVector,
4677  typename Container>
4678  struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4679  {
4680  using function_type =
4681  void (Container::*)(const MF &,
4682  OutVector &,
4683  const InVector &,
4684  const std::pair<unsigned int, unsigned int> &);
4685  };
4686  } // namespace MatrixFreeFunctions
4687 
4688 
4689 
4690  // A implementation class for the worker object that runs the various
4691  // operations we want to perform during the matrix-free loop
4692  template <typename MF,
4693  typename InVector,
4694  typename OutVector,
4695  typename Container,
4696  bool is_constant>
4697  class MFWorker : public MFWorkerInterface
4698  {
4699  public:
4700  // An alias to make the arguments further down more readable
4701  using function_type = typename MatrixFreeFunctions::
4702  InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4703  function_type;
4704 
4705  // constructor, binds all the arguments to this class
4706  MFWorker(const MF & matrix_free,
4707  const InVector & src,
4708  OutVector & dst,
4709  const bool zero_dst_vector_setting,
4710  const Container & container,
4711  function_type cell_function,
4712  function_type face_function,
4713  function_type boundary_function,
4714  const typename MF::DataAccessOnFaces src_vector_face_access =
4716  const typename MF::DataAccessOnFaces dst_vector_face_access =
4718  const std::function<void(const unsigned int, const unsigned int)>
4719  &operation_before_loop = {},
4720  const std::function<void(const unsigned int, const unsigned int)>
4721  & operation_after_loop = {},
4722  const unsigned int dof_handler_index_pre_post = 0)
4723  : matrix_free(matrix_free)
4724  , container(const_cast<Container &>(container))
4725  , cell_function(cell_function)
4726  , face_function(face_function)
4727  , boundary_function(boundary_function)
4728  , src(src)
4729  , dst(dst)
4730  , src_data_exchanger(matrix_free,
4731  src_vector_face_access,
4732  n_components(src))
4733  , dst_data_exchanger(matrix_free,
4734  dst_vector_face_access,
4735  n_components(dst))
4736  , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4737  , zero_dst_vector_setting(zero_dst_vector_setting &&
4738  !src_and_dst_are_same)
4739  , operation_before_loop(operation_before_loop)
4740  , operation_after_loop(operation_after_loop)
4741  , dof_handler_index_pre_post(dof_handler_index_pre_post)
4742  {}
4743 
4744  // Runs the cell work. If no function is given, nothing is done
4745  virtual void
4746  cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4747  {
4748  if (cell_function != nullptr && cell_range.second > cell_range.first)
4749  for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4750  {
4751  const auto cell_subrange =
4752  matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4753 
4754  if (cell_subrange.second <= cell_subrange.first)
4755  continue;
4756 
4757  (container.*
4758  cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4759  }
4760  }
4761 
4762  virtual void
4763  cell(const unsigned int range_index) override
4764  {
4765  process_range(cell_function,
4766  matrix_free.get_task_info().cell_partition_data_hp_ptr,
4767  matrix_free.get_task_info().cell_partition_data_hp,
4768  range_index);
4769  }
4770 
4771  virtual void
4772  face(const unsigned int range_index) override
4773  {
4774  process_range(face_function,
4775  matrix_free.get_task_info().face_partition_data_hp_ptr,
4776  matrix_free.get_task_info().face_partition_data_hp,
4777  range_index);
4778  }
4779 
4780  virtual void
4781  boundary(const unsigned int range_index) override
4782  {
4783  process_range(boundary_function,
4784  matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4785  matrix_free.get_task_info().boundary_partition_data_hp,
4786  range_index);
4787  }
4788 
4789  private:
4790  void
4791  process_range(const function_type & fu,
4792  const std::vector<unsigned int> &ptr,
4793  const std::vector<unsigned int> &data,
4794  const unsigned int range_index)
4795  {
4796  if (fu == nullptr)
4797  return;
4798 
4799  AssertIndexRange(range_index + 1, ptr.size());
4800  for (unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4801  {
4802  AssertIndexRange(2 * i + 1, data.size());
4803  (container.*fu)(matrix_free,
4804  this->dst,
4805  this->src,
4806  std::make_pair(data[2 * i], data[2 * i + 1]));
4807  }
4808  }
4809 
4810  public:
4811  // Starts the communication for the update ghost values operation. We
4812  // cannot call this update if ghost and destination are the same because
4813  // that would introduce spurious entries in the destination (there is also
4814  // the problem that reading from a vector that we also write to is usually
4815  // not intended in case there is overlap, but this is up to the
4816  // application code to decide and we cannot catch this case here).
4817  virtual void
4818  vector_update_ghosts_start() override
4819  {
4820  if (!src_and_dst_are_same)
4821  internal::update_ghost_values_start(src, src_data_exchanger);
4822  }
4823 
4824  // Finishes the communication for the update ghost values operation
4825  virtual void
4826  vector_update_ghosts_finish() override
4827  {
4828  if (!src_and_dst_are_same)
4829  internal::update_ghost_values_finish(src, src_data_exchanger);
4830  }
4831 
4832  // Starts the communication for the vector compress operation
4833  virtual void
4834  vector_compress_start() override
4835  {
4836  internal::compress_start(dst, dst_data_exchanger);
4837  }
4838 
4839  // Finishes the communication for the vector compress operation
4840  virtual void
4841  vector_compress_finish() override
4842  {
4843  internal::compress_finish(dst, dst_data_exchanger);
4844  if (!src_and_dst_are_same)
4845  internal::reset_ghost_values(src, src_data_exchanger);
4846  }
4847 
4848  // Zeros the given input vector
4849  virtual void
4850  zero_dst_vector_range(const unsigned int range_index) override
4851  {
4852  if (zero_dst_vector_setting)
4853  internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4854  }
4855 
4856  virtual void
4857  cell_loop_pre_range(const unsigned int range_index) override
4858  {
4859  if (operation_before_loop)
4860  {
4861  const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4862  matrix_free.get_dof_info(dof_handler_index_pre_post);
4863  if (range_index == numbers::invalid_unsigned_int)
4864  {
4865  // Case with threaded loop -> currently no overlap implemented
4867  0U,
4868  dof_info.vector_partitioner->locally_owned_size(),
4869  operation_before_loop,
4871  }
4872  else
4873  {
4874  AssertIndexRange(range_index,
4875  dof_info.cell_loop_pre_list_index.size() - 1);
4876  for (unsigned int id =
4877  dof_info.cell_loop_pre_list_index[range_index];
4878  id != dof_info.cell_loop_pre_list_index[range_index + 1];
4879  ++id)
4880  operation_before_loop(dof_info.cell_loop_pre_list[id].first,
4881  dof_info.cell_loop_pre_list[id].second);
4882  }
4883  }
4884  }
4885 
4886  virtual void
4887  cell_loop_post_range(const unsigned int range_index) override
4888  {
4889  if (operation_after_loop)
4890  {
4891  // Run unit matrix operation on constrained dofs if we are at the
4892  // last range
4893  const std::vector<unsigned int> &partition_row_index =
4894  matrix_free.get_task_info().partition_row_index;
4895  if (range_index ==
4896  partition_row_index[partition_row_index.size() - 2] - 1)
4897  apply_operation_to_constrained_dofs(
4898  matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
4899  src,
4900  dst);
4901 
4902  const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4903  matrix_free.get_dof_info(dof_handler_index_pre_post);
4904  if (range_index == numbers::invalid_unsigned_int)
4905  {
4906  // Case with threaded loop -> currently no overlap implemented
4908  0U,
4909  dof_info.vector_partitioner->locally_owned_size(),
4910  operation_after_loop,
4912  }
4913  else
4914  {
4915  AssertIndexRange(range_index,
4916  dof_info.cell_loop_post_list_index.size() - 1);
4917  for (unsigned int id =
4918  dof_info.cell_loop_post_list_index[range_index];
4919  id != dof_info.cell_loop_post_list_index[range_index + 1];
4920  ++id)
4921  operation_after_loop(dof_info.cell_loop_post_list[id].first,
4922  dof_info.cell_loop_post_list[id].second);
4923  }
4924  }
4925  }
4926 
4927  private:
4928  const MF & matrix_free;
4929  Container & container;
4930  function_type cell_function;
4931  function_type face_function;
4932  function_type boundary_function;
4933 
4934  const InVector &src;
4935  OutVector & dst;
4936  VectorDataExchange<MF::dimension,
4937  typename MF::value_type,
4938  typename MF::vectorized_value_type>
4939  src_data_exchanger;
4940  VectorDataExchange<MF::dimension,
4941  typename MF::value_type,
4942  typename MF::vectorized_value_type>
4943  dst_data_exchanger;
4944  const bool src_and_dst_are_same;
4945  const bool zero_dst_vector_setting;
4946  const std::function<void(const unsigned int, const unsigned int)>
4947  operation_before_loop;
4948  const std::function<void(const unsigned int, const unsigned int)>
4949  operation_after_loop;
4950  const unsigned int dof_handler_index_pre_post;
4951  };
4952 
4953 
4954 
4959  template <class MF, typename InVector, typename OutVector>
4960  struct MFClassWrapper
4961  {
4962  using function_type =
4963  std::function<void(const MF &,
4964  OutVector &,
4965  const InVector &,
4966  const std::pair<unsigned int, unsigned int> &)>;
4967 
4968  MFClassWrapper(const function_type cell,
4969  const function_type face,
4970  const function_type boundary)
4971  : cell(cell)
4972  , face(face)
4973  , boundary(boundary)
4974  {}
4975 
4976  void
4977  cell_integrator(const MF & mf,
4978  OutVector & dst,
4979  const InVector & src,
4980  const std::pair<unsigned int, unsigned int> &range) const
4981  {
4982  if (cell)
4983  cell(mf, dst, src, range);
4984  }
4985 
4986  void
4987  face_integrator(const MF & mf,
4988  OutVector & dst,
4989  const InVector & src,
4990  const std::pair<unsigned int, unsigned int> &range) const
4991  {
4992  if (face)
4993  face(mf, dst, src, range);
4994  }
4995 
4996  void
4997  boundary_integrator(
4998  const MF & mf,
4999  OutVector & dst,
5000  const InVector & src,
5001  const std::pair<unsigned int, unsigned int> &range) const
5002  {
5003  if (boundary)
5004  boundary(mf, dst, src, range);
5005  }
5006 
5007  const function_type cell;
5008  const function_type face;
5009  const function_type boundary;
5010  };
5011 
5012 } // end of namespace internal
5013 
5014 
5015 
5016 template <int dim, typename Number, typename VectorizedArrayType>
5017 template <typename OutVector, typename InVector>
5018 inline void
5020  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5021  OutVector &,
5022  const InVector &,
5023  const std::pair<unsigned int, unsigned int> &)>
5024  & cell_operation,
5025  OutVector & dst,
5026  const InVector &src,
5027  const bool zero_dst_vector) const
5028 {
5029  using Wrapper =
5030  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5031  InVector,
5032  OutVector>;
5033  Wrapper wrap(cell_operation, nullptr, nullptr);
5034  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5035  InVector,
5036  OutVector,
5037  Wrapper,
5038  true>
5039  worker(*this,
5040  src,
5041  dst,
5042  zero_dst_vector,
5043  wrap,
5044  &Wrapper::cell_integrator,
5045  &Wrapper::face_integrator,
5046  &Wrapper::boundary_integrator);
5047 
5048  task_info.loop(worker);
5049 }
5050 
5051 
5052 
5053 template <int dim, typename Number, typename VectorizedArrayType>
5054 template <typename OutVector, typename InVector>
5055 inline void
5057  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5058  OutVector &,
5059  const InVector &,
5060  const std::pair<unsigned int, unsigned int> &)>
5061  & cell_operation,
5062  OutVector & dst,
5063  const InVector &src,
5064  const std::function<void(const unsigned int, const unsigned int)>
5065  &operation_before_loop,
5066  const std::function<void(const unsigned int, const unsigned int)>
5067  & operation_after_loop,
5068  const unsigned int dof_handler_index_pre_post) const
5069 {
5070  using Wrapper =
5071  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5072  InVector,
5073  OutVector>;
5074  Wrapper wrap(cell_operation, nullptr, nullptr);
5075  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5076  InVector,
5077  OutVector,
5078  Wrapper,
5079  true>
5080  worker(*this,
5081  src,
5082  dst,
5083  false,
5084  wrap,
5085  &Wrapper::cell_integrator,
5086  &Wrapper::face_integrator,
5087  &Wrapper::boundary_integrator,
5090  operation_before_loop,
5091  operation_after_loop,
5092  dof_handler_index_pre_post);
5093 
5094  task_info.loop(worker);
5095 }
5096 
5097 
5098 
5099 template <int dim, typename Number, typename VectorizedArrayType>
5100 template <typename OutVector, typename InVector>
5101 inline void
5103  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5104  OutVector &,
5105  const InVector &,
5106  const std::pair<unsigned int, unsigned int> &)>
5107  &cell_operation,
5108  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5109  OutVector &,
5110  const InVector &,
5111  const std::pair<unsigned int, unsigned int> &)>
5112  &face_operation,
5113  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5114  OutVector &,
5115  const InVector &,
5116  const std::pair<unsigned int, unsigned int> &)>
5117  & boundary_operation,
5118  OutVector & dst,
5119  const InVector & src,
5120  const bool zero_dst_vector,
5121  const DataAccessOnFaces dst_vector_face_access,
5122  const DataAccessOnFaces src_vector_face_access) const
5123 {
5124  using Wrapper =
5125  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5126  InVector,
5127  OutVector>;
5128  Wrapper wrap(cell_operation, face_operation, boundary_operation);
5129  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5130  InVector,
5131  OutVector,
5132  Wrapper,
5133  true>
5134  worker(*this,
5135  src,
5136  dst,
5137  zero_dst_vector,
5138  wrap,
5139  &Wrapper::cell_integrator,
5140  &Wrapper::face_integrator,
5141  &Wrapper::boundary_integrator,
5142  src_vector_face_access,
5143  dst_vector_face_access);
5144 
5145  task_info.loop(worker);
5146 }
5147 
5148 
5149 
5150 template <int dim, typename Number, typename VectorizedArrayType>
5151 template <typename CLASS, typename OutVector, typename InVector>
5152 inline void
5154  void (CLASS::*function_pointer)(
5156  OutVector &,
5157  const InVector &,
5158  const std::pair<unsigned int, unsigned int> &) const,
5159  const CLASS * owning_class,
5160  OutVector & dst,
5161  const InVector &src,
5162  const bool zero_dst_vector) const
5163 {
5164  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5165  InVector,
5166  OutVector,
5167  CLASS,
5168  true>
5169  worker(*this,
5170  src,
5171  dst,
5172  zero_dst_vector,
5173  *owning_class,
5174  function_pointer,
5175  nullptr,
5176  nullptr);
5177  task_info.loop(worker);
5178 }
5179 
5180 
5181 
5182 template <int dim, typename Number, typename VectorizedArrayType>
5183 template <typename CLASS, typename OutVector, typename InVector>
5184 inline void
5186  void (CLASS::*function_pointer)(
5188  OutVector &,
5189  const InVector &,
5190  const std::pair<unsigned int, unsigned int> &) const,
5191  const CLASS * owning_class,
5192  OutVector & dst,
5193  const InVector &src,
5194  const std::function<void(const unsigned int, const unsigned int)>
5195  &operation_before_loop,
5196  const std::function<void(const unsigned int, const unsigned int)>
5197  & operation_after_loop,
5198  const unsigned int dof_handler_index_pre_post) const
5199 {
5200  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5201  InVector,
5202  OutVector,
5203  CLASS,
5204  true>
5205  worker(*this,
5206  src,
5207  dst,
5208  false,
5209  *owning_class,
5210  function_pointer,
5211  nullptr,
5212  nullptr,
5215  operation_before_loop,
5216  operation_after_loop,
5217  dof_handler_index_pre_post);
5218  task_info.loop(worker);
5219 }
5220 
5221 
5222 
5223 template <int dim, typename Number, typename VectorizedArrayType>
5224 template <typename CLASS, typename OutVector, typename InVector>
5225 inline void
5227  void (CLASS::*cell_operation)(
5229  OutVector &,
5230  const InVector &,
5231  const std::pair<unsigned int, unsigned int> &) const,
5232  void (CLASS::*face_operation)(
5234  OutVector &,
5235  const InVector &,
5236  const std::pair<unsigned int, unsigned int> &) const,
5237  void (CLASS::*boundary_operation)(
5239  OutVector &,
5240  const InVector &,
5241  const std::pair<unsigned int, unsigned int> &) const,
5242  const CLASS * owning_class,
5243  OutVector & dst,
5244  const InVector & src,
5245  const bool zero_dst_vector,
5246  const DataAccessOnFaces dst_vector_face_access,
5247  const DataAccessOnFaces src_vector_face_access) const
5248 {
5249  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5250  InVector,
5251  OutVector,
5252  CLASS,
5253  true>
5254  worker(*this,
5255  src,
5256  dst,
5257  zero_dst_vector,
5258  *owning_class,
5259  cell_operation,
5260  face_operation,
5261  boundary_operation,
5262  src_vector_face_access,
5263  dst_vector_face_access);
5264  task_info.loop(worker);
5265 }
5266 
5267 
5268 
5269 template <int dim, typename Number, typename VectorizedArrayType>
5270 template <typename CLASS, typename OutVector, typename InVector>
5271 inline void
5273  void (CLASS::*function_pointer)(
5275  OutVector &,
5276  const InVector &,
5277  const std::pair<unsigned int, unsigned int> &),
5278  CLASS * owning_class,
5279  OutVector & dst,
5280  const InVector &src,
5281  const bool zero_dst_vector) const
5282 {
5283  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5284  InVector,
5285  OutVector,
5286  CLASS,
5287  false>
5288  worker(*this,
5289  src,
5290  dst,
5291  zero_dst_vector,
5292  *owning_class,
5293  function_pointer,
5294  nullptr,
5295  nullptr);
5296  task_info.loop(worker);
5297 }
5298 
5299 
5300 
5301 template <int dim, typename Number, typename VectorizedArrayType>
5302 template <typename CLASS, typename OutVector, typename InVector>
5303 inline void
5305  void (CLASS::*function_pointer)(
5307  OutVector &,
5308  const InVector &,
5309  const std::pair<unsigned int, unsigned int> &),
5310  CLASS * owning_class,
5311  OutVector & dst,
5312  const InVector &src,
5313  const std::function<void(const unsigned int, const unsigned int)>
5314  &operation_before_loop,
5315  const std::function<void(const unsigned int, const unsigned int)>
5316  & operation_after_loop,
5317  const unsigned int dof_handler_index_pre_post) const
5318 {
5319  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5320  InVector,
5321  OutVector,
5322  CLASS,
5323  false>
5324  worker(*this,
5325  src,
5326  dst,
5327  false,
5328  *owning_class,
5329  function_pointer,
5330  nullptr,
5331  nullptr,
5334  operation_before_loop,
5335  operation_after_loop,
5336  dof_handler_index_pre_post);
5337  task_info.loop(worker);
5338 }
5339 
5340 
5341 
5342 template <int dim, typename Number, typename VectorizedArrayType>
5343 template <typename CLASS, typename OutVector, typename InVector>
5344 inline void
5346  void (CLASS::*cell_operation)(
5348  OutVector &,
5349  const InVector &,
5350  const std::pair<unsigned int, unsigned int> &),
5351  void (CLASS::*face_operation)(
5353  OutVector &,
5354  const InVector &,
5355  const std::pair<unsigned int, unsigned int> &),
5356  void (CLASS::*boundary_operation)(
5358  OutVector &,
5359  const InVector &,
5360  const std::pair<unsigned int, unsigned int> &),
5361  CLASS * owning_class,
5362  OutVector & dst,
5363  const InVector & src,
5364  const bool zero_dst_vector,
5365  const DataAccessOnFaces dst_vector_face_access,
5366  const DataAccessOnFaces src_vector_face_access) const
5367 {
5368  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5369  InVector,
5370  OutVector,
5371  CLASS,
5372  false>
5373  worker(*this,
5374  src,
5375  dst,
5376  zero_dst_vector,
5377  *owning_class,
5378  cell_operation,
5379  face_operation,
5380  boundary_operation,
5381  src_vector_face_access,
5382  dst_vector_face_access);
5383  task_info.loop(worker);
5384 }
5385 
5386 
5387 
5388 template <int dim, typename Number, typename VectorizedArrayType>
5389 template <typename CLASS, typename OutVector, typename InVector>
5390 inline void
5392  void (CLASS::*function_pointer)(
5394  OutVector &,
5395  const InVector &,
5396  const std::pair<unsigned int, unsigned int> &) const,
5397  const CLASS * owning_class,
5398  OutVector & dst,
5399  const InVector & src,
5400  const bool zero_dst_vector,
5401  const DataAccessOnFaces src_vector_face_access) const
5402 {
5403  auto src_vector_face_access_temp = src_vector_face_access;
5404  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5405  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5406  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5407  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5408 
5409  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5410  InVector,
5411  OutVector,
5412  CLASS,
5413  true>
5414  worker(*this,
5415  src,
5416  dst,
5417  zero_dst_vector,
5418  *owning_class,
5419  function_pointer,
5420  nullptr,
5421  nullptr,
5422  src_vector_face_access_temp,
5424  task_info.loop(worker);
5425 }
5426 
5427 
5428 
5429 template <int dim, typename Number, typename VectorizedArrayType>
5430 template <typename CLASS, typename OutVector, typename InVector>
5431 inline void
5433  void (CLASS::*function_pointer)(
5435  OutVector &,
5436  const InVector &,
5437  const std::pair<unsigned int, unsigned int> &),
5438  CLASS * owning_class,
5439  OutVector & dst,
5440  const InVector & src,
5441  const bool zero_dst_vector,
5442  const DataAccessOnFaces src_vector_face_access) const
5443 {
5444  auto src_vector_face_access_temp = src_vector_face_access;
5445  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5446  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5447  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5448  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5449 
5450  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5451  InVector,
5452  OutVector,
5453  CLASS,
5454  false>
5455  worker(*this,
5456  src,
5457  dst,
5458  zero_dst_vector,
5459  *owning_class,
5460  function_pointer,
5461  nullptr,
5462  nullptr,
5463  src_vector_face_access_temp,
5465  task_info.loop(worker);
5466 }
5467 
5468 
5469 
5470 template <int dim, typename Number, typename VectorizedArrayType>
5471 template <typename OutVector, typename InVector>
5472 inline void
5474  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5475  OutVector &,
5476  const InVector &,
5477  const std::pair<unsigned int, unsigned int> &)>
5478  & cell_operation,
5479  OutVector & dst,
5480  const InVector & src,
5481  const bool zero_dst_vector,
5482  const DataAccessOnFaces src_vector_face_access) const
5483 {
5484  auto src_vector_face_access_temp = src_vector_face_access;
5485  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5486  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5487  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5488  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5489 
5490  using Wrapper =
5491  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5492  InVector,
5493  OutVector>;
5494  Wrapper wrap(cell_operation, nullptr, nullptr);
5495 
5496  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5497  InVector,
5498  OutVector,
5499  Wrapper,
5500  true>
5501  worker(*this,
5502  src,
5503  dst,
5504  zero_dst_vector,
5505  wrap,
5506  &Wrapper::cell_integrator,
5507  &Wrapper::face_integrator,
5508  &Wrapper::boundary_integrator,
5509  src_vector_face_access_temp,
5511  task_info.loop(worker);
5512 }
5513 
5514 
5515 #endif // ifndef DOXYGEN
5516 
5517 
5518 
5520 
5521 #endif
void resize(const size_type new_size)
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:2157
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:2163
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
void reinit(const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
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
DoFHandler< dim >::active_cell_iterator get_hp_cell_iterator(const unsigned int cell_batch_index, const unsigned int lane_index, const unsigned int dof_handler_index=0) const
void update_mapping(const Mapping< dim > &mapping)
unsigned int get_mg_level() const
unsigned int n_components_filled(const unsigned int cell_batch_number) const
void print_memory_consumption(StreamType &out) const
void update_mapping(const std::shared_ptr< hp::MappingCollection< dim >> &mapping)
void reinit(const std::vector< const DoFHandlerType * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< QuadratureType > &quad, const AdditionalData &additional_data=AdditionalData())
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:2187
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:2204
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:2227
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
const DoFHandlerType & get_dof_handler(const unsigned int dof_handler_index=0) const
std::vector< std::pair< unsigned int, unsigned int > > cell_level_index
Definition: matrix_free.h:2171
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
void reinit(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< Number > constraint_pool_data
Definition: matrix_free.h:2144
VectorizedArrayType vectorized_value_type
Definition: matrix_free.h:126
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
void reinit(const MappingType &mapping, const std::vector< const DoFHandlerType * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
Threads::ThreadLocalStorage< std::list< std::pair< bool, AlignedVector< VectorizedArrayType > > > > scratch_pad
Definition: matrix_free.h:2215
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_index) const
Number value_type
Definition: matrix_free.h:125
std::size_t memory_consumption() 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:2222
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:2130
void reinit(const std::vector< const DoFHandlerType * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
void initialize_indices(const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const AdditionalData &additional_data)
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())
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
unsigned int n_macro_cells() const
void reinit(const MappingType &mapping, const std::vector< const DoFHandlerType * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< QuadratureType > &quad, const AdditionalData &additional_data=AdditionalData())
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:2136
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:2199
internal::MatrixFreeFunctions::FaceInfo< VectorizedArrayType::size()> face_info
Definition: matrix_free.h:2194
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:2150
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner(const unsigned int dof_handler_index=0) const
void reinit(const DoFHandler< dim > &dof_handler, const AffineConstraints< number2 > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
unsigned int cell_level_index_end_local
Definition: matrix_free.h:2180
unsigned int n_base_elements(const unsigned int dof_handler_index) 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:131
A class that provides a separate storage location on each thread that accesses the object.
#define DEAL_II_DEPRECATED
Definition: config.h:164
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define DEAL_II_WITH_MPI
Definition: config.h:54
UpdateFlags
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_default
No update.
Point< 2 > second
Definition: grid_out.cc:4604
Point< 2 > first
Definition: grid_out.cc:4603
unsigned int level
Definition: grid_out.cc:4606
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1790
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
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:1583
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:466
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:438
Number local_element(const size_type local_index) const
void reinit(const size_type size, const bool omit_zeroing_entries=false)
static const char U
VectorType::value_type * begin(VectorType &V)
bool job_supports_mpi()
Definition: mpi.cc:1014
unsigned int minimum_parallel_grain_size
Definition: parallel.cc:34
const types::boundary_id invalid_boundary_id
Definition: types.h:244
static const unsigned int invalid_unsigned_int
Definition: types.h:201
void apply_to_subranges(const tbb::blocked_range< RangeType > &range, const Function &f)
Definition: parallel.h:302
unsigned int boundary_id
Definition: types.h:129
TasksParallelScheme tasks_parallel_scheme
Definition: matrix_free.h:343
UpdateFlags mapping_update_flags_inner_faces
Definition: matrix_free.h:407
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:217
AdditionalData & operator=(const AdditionalData &other)
Definition: matrix_free.h:281
std::vector< unsigned int > cell_vectorization_category
Definition: matrix_free.h:514
UpdateFlags mapping_update_flags_boundary_faces
Definition: matrix_free.h:387
UpdateFlags mapping_update_flags
Definition: matrix_free.h:367
UpdateFlags mapping_update_flags_faces_by_cells
Definition: matrix_free.h:435
AdditionalData(const AdditionalData &other)
Definition: matrix_free.h:254
unsigned int tasks_block_size
Definition: matrix_free.h:354
static bool equal(const T *p1, const T *p2)
std::vector< unsigned int > cell_loop_pre_list_index
Definition: dof_info.h:777
std::vector< unsigned int > cell_loop_post_list_index
Definition: dof_info.h:790
std::vector< std::pair< unsigned int, unsigned int > > vector_zero_range_list
Definition: dof_info.h:770
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition: dof_info.h:621
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_pre_list
Definition: dof_info.h:783
std::vector< unsigned int > vector_zero_range_list_index
Definition: dof_info.h:765
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_post_list
Definition: dof_info.h:796
::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:521
void loop(MFWorkerInterface &worker) const
Definition: task_info.cc:350
std::vector< unsigned int > cell_partition_data
Definition: task_info.h:477
std::vector< unsigned int > face_partition_data
Definition: task_info.h:499