Reference documentation for deal.II version 9.3.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 - 2021 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 
36 
37 #include <deal.II/hp/dof_handler.h>
40 
45 
51 
52 #include <cstdlib>
53 #include <limits>
54 #include <list>
55 #include <memory>
56 
57 
59 
60 
61 
113 template <int dim,
114  typename Number = double,
115  typename VectorizedArrayType = VectorizedArray<Number>>
116 class MatrixFree : public Subscriptor
117 {
118  static_assert(
119  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
120  "Type of Number and of VectorizedArrayType do not match.");
121 
122 public:
127  using value_type = Number;
128  using vectorized_value_type = VectorizedArrayType;
129 
133  static const unsigned int dimension = dim;
134 
182  {
189  {
209  };
210 
211  // remove with level_mg_handler
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)
238  , mg_level(mg_level)
239  , level_mg_handler(this->mg_level)
247  , communicator_sm(MPI_COMM_SELF)
248  {}
249 
262  , mg_level(other.mg_level)
263  , level_mg_handler(this->mg_level)
274  {}
275 
276  // remove with level_mg_handler
278 
283  operator=(const AdditionalData &other)
284  {
293  mg_level = other.mg_level;
304 
305  return *this;
306  }
307 
345 
355  unsigned int tasks_block_size;
356 
369 
390 
411 
439 
448  unsigned int mg_level;
449 
456 
464 
475 
484 
498 
507 
524  std::vector<unsigned int> cell_vectorization_category;
525 
536 
541  };
542 
550  MatrixFree();
551 
556 
560  ~MatrixFree() override = default;
561 
574  template <typename QuadratureType, typename number2, typename MappingType>
575  void
576  reinit(const MappingType & mapping,
577  const DoFHandler<dim> & dof_handler,
578  const AffineConstraints<number2> &constraint,
579  const QuadratureType & quad,
580  const AdditionalData & additional_data = AdditionalData());
581 
588  template <typename QuadratureType, typename number2>
589  DEAL_II_DEPRECATED void
590  reinit(const DoFHandler<dim> & dof_handler,
591  const AffineConstraints<number2> &constraint,
592  const QuadratureType & quad,
593  const AdditionalData & additional_data = AdditionalData());
594 
616  template <typename QuadratureType, typename number2, typename MappingType>
617  void
618  reinit(const MappingType & mapping,
619  const std::vector<const DoFHandler<dim> *> & dof_handler,
620  const std::vector<const AffineConstraints<number2> *> &constraint,
621  const std::vector<QuadratureType> & quad,
622  const AdditionalData &additional_data = AdditionalData());
623 
629  template <typename QuadratureType,
630  typename number2,
631  typename DoFHandlerType,
632  typename MappingType>
633  DEAL_II_DEPRECATED void
634  reinit(const MappingType & mapping,
635  const std::vector<const DoFHandlerType *> & dof_handler,
636  const std::vector<const AffineConstraints<number2> *> &constraint,
637  const std::vector<QuadratureType> & quad,
638  const AdditionalData &additional_data = AdditionalData());
639 
646  template <typename QuadratureType, typename number2>
647  DEAL_II_DEPRECATED void
648  reinit(const std::vector<const DoFHandler<dim> *> & dof_handler,
649  const std::vector<const AffineConstraints<number2> *> &constraint,
650  const std::vector<QuadratureType> & quad,
651  const AdditionalData &additional_data = AdditionalData());
652 
658  template <typename QuadratureType, typename number2, typename DoFHandlerType>
659  DEAL_II_DEPRECATED void
660  reinit(const std::vector<const DoFHandlerType *> & dof_handler,
661  const std::vector<const AffineConstraints<number2> *> &constraint,
662  const std::vector<QuadratureType> & quad,
663  const AdditionalData &additional_data = AdditionalData());
664 
672  template <typename QuadratureType, typename number2, typename MappingType>
673  void
674  reinit(const MappingType & mapping,
675  const std::vector<const DoFHandler<dim> *> & dof_handler,
676  const std::vector<const AffineConstraints<number2> *> &constraint,
677  const QuadratureType & quad,
678  const AdditionalData &additional_data = AdditionalData());
679 
685  template <typename QuadratureType,
686  typename number2,
687  typename DoFHandlerType,
688  typename MappingType>
689  DEAL_II_DEPRECATED void
690  reinit(const MappingType & mapping,
691  const std::vector<const DoFHandlerType *> & dof_handler,
692  const std::vector<const AffineConstraints<number2> *> &constraint,
693  const QuadratureType & quad,
694  const AdditionalData &additional_data = AdditionalData());
695 
702  template <typename QuadratureType, typename number2>
703  DEAL_II_DEPRECATED void
704  reinit(const std::vector<const DoFHandler<dim> *> & dof_handler,
705  const std::vector<const AffineConstraints<number2> *> &constraint,
706  const QuadratureType & quad,
707  const AdditionalData &additional_data = AdditionalData());
708 
714  template <typename QuadratureType, typename number2, typename DoFHandlerType>
715  DEAL_II_DEPRECATED void
716  reinit(const std::vector<const DoFHandlerType *> & dof_handler,
717  const std::vector<const AffineConstraints<number2> *> &constraint,
718  const QuadratureType & quad,
719  const AdditionalData &additional_data = AdditionalData());
720 
726  void
727  copy_from(
728  const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free_base);
729 
739  void
740  update_mapping(const Mapping<dim> &mapping);
741 
745  void
746  update_mapping(const std::shared_ptr<hp::MappingCollection<dim>> &mapping);
747 
752  void
753  clear();
754 
756 
768  enum class DataAccessOnFaces
769  {
776  none,
777 
788  values,
789 
798  values_all_faces,
799 
810  gradients,
811 
820  gradients_all_faces,
821 
828  unspecified
829  };
830 
875  template <typename OutVector, typename InVector>
876  void
877  cell_loop(const std::function<void(
879  OutVector &,
880  const InVector &,
881  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
882  OutVector & dst,
883  const InVector & src,
884  const bool zero_dst_vector = false) const;
885 
932  template <typename CLASS, typename OutVector, typename InVector>
933  void
934  cell_loop(void (CLASS::*cell_operation)(
935  const MatrixFree &,
936  OutVector &,
937  const InVector &,
938  const std::pair<unsigned int, unsigned int> &) const,
939  const CLASS * owning_class,
940  OutVector & dst,
941  const InVector &src,
942  const bool zero_dst_vector = false) const;
943 
947  template <typename CLASS, typename OutVector, typename InVector>
948  void
949  cell_loop(void (CLASS::*cell_operation)(
950  const MatrixFree &,
951  OutVector &,
952  const InVector &,
953  const std::pair<unsigned int, unsigned int> &),
954  CLASS * owning_class,
955  OutVector & dst,
956  const InVector &src,
957  const bool zero_dst_vector = false) const;
958 
1043  template <typename CLASS, typename OutVector, typename InVector>
1044  void
1045  cell_loop(void (CLASS::*cell_operation)(
1046  const MatrixFree &,
1047  OutVector &,
1048  const InVector &,
1049  const std::pair<unsigned int, unsigned int> &) const,
1050  const CLASS * owning_class,
1051  OutVector & dst,
1052  const InVector &src,
1053  const std::function<void(const unsigned int, const unsigned int)>
1054  &operation_before_loop,
1055  const std::function<void(const unsigned int, const unsigned int)>
1056  & operation_after_loop,
1057  const unsigned int dof_handler_index_pre_post = 0) const;
1058 
1062  template <typename CLASS, typename OutVector, typename InVector>
1063  void
1064  cell_loop(void (CLASS::*cell_operation)(
1065  const MatrixFree &,
1066  OutVector &,
1067  const InVector &,
1068  const std::pair<unsigned int, unsigned int> &),
1069  CLASS * owning_class,
1070  OutVector & dst,
1071  const InVector &src,
1072  const std::function<void(const unsigned int, const unsigned int)>
1073  &operation_before_loop,
1074  const std::function<void(const unsigned int, const unsigned int)>
1075  & operation_after_loop,
1076  const unsigned int dof_handler_index_pre_post = 0) const;
1077 
1082  template <typename OutVector, typename InVector>
1083  void
1084  cell_loop(const std::function<void(
1086  OutVector &,
1087  const InVector &,
1088  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1089  OutVector & dst,
1090  const InVector & src,
1091  const std::function<void(const unsigned int, const unsigned int)>
1092  &operation_before_loop,
1093  const std::function<void(const unsigned int, const unsigned int)>
1094  & operation_after_loop,
1095  const unsigned int dof_handler_index_pre_post = 0) const;
1096 
1172  template <typename OutVector, typename InVector>
1173  void
1174  loop(const std::function<
1176  OutVector &,
1177  const InVector &,
1178  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1179  const std::function<
1181  OutVector &,
1182  const InVector &,
1183  const std::pair<unsigned int, unsigned int> &)> &face_operation,
1184  const std::function<void(
1186  OutVector &,
1187  const InVector &,
1188  const std::pair<unsigned int, unsigned int> &)> &boundary_operation,
1189  OutVector & dst,
1190  const InVector & src,
1191  const bool zero_dst_vector = false,
1192  const DataAccessOnFaces dst_vector_face_access =
1194  const DataAccessOnFaces src_vector_face_access =
1196 
1282  template <typename CLASS, typename OutVector, typename InVector>
1283  void
1284  loop(
1285  void (CLASS::*cell_operation)(const MatrixFree &,
1286  OutVector &,
1287  const InVector &,
1288  const std::pair<unsigned int, unsigned int> &)
1289  const,
1290  void (CLASS::*face_operation)(const MatrixFree &,
1291  OutVector &,
1292  const InVector &,
1293  const std::pair<unsigned int, unsigned int> &)
1294  const,
1295  void (CLASS::*boundary_operation)(
1296  const MatrixFree &,
1297  OutVector &,
1298  const InVector &,
1299  const std::pair<unsigned int, unsigned int> &) const,
1300  const CLASS * owning_class,
1301  OutVector & dst,
1302  const InVector & src,
1303  const bool zero_dst_vector = false,
1304  const DataAccessOnFaces dst_vector_face_access =
1306  const DataAccessOnFaces src_vector_face_access =
1308 
1312  template <typename CLASS, typename OutVector, typename InVector>
1313  void
1314  loop(void (CLASS::*cell_operation)(
1315  const MatrixFree &,
1316  OutVector &,
1317  const InVector &,
1318  const std::pair<unsigned int, unsigned int> &),
1319  void (CLASS::*face_operation)(
1320  const MatrixFree &,
1321  OutVector &,
1322  const InVector &,
1323  const std::pair<unsigned int, unsigned int> &),
1324  void (CLASS::*boundary_operation)(
1325  const MatrixFree &,
1326  OutVector &,
1327  const InVector &,
1328  const std::pair<unsigned int, unsigned int> &),
1329  CLASS * owning_class,
1330  OutVector & dst,
1331  const InVector & src,
1332  const bool zero_dst_vector = false,
1333  const DataAccessOnFaces dst_vector_face_access =
1335  const DataAccessOnFaces src_vector_face_access =
1337 
1402  template <typename CLASS, typename OutVector, typename InVector>
1403  void
1404  loop_cell_centric(void (CLASS::*cell_operation)(
1405  const MatrixFree &,
1406  OutVector &,
1407  const InVector &,
1408  const std::pair<unsigned int, unsigned int> &) const,
1409  const CLASS * owning_class,
1410  OutVector & dst,
1411  const InVector & src,
1412  const bool zero_dst_vector = false,
1413  const DataAccessOnFaces src_vector_face_access =
1415 
1419  template <typename CLASS, typename OutVector, typename InVector>
1420  void
1421  loop_cell_centric(void (CLASS::*cell_operation)(
1422  const MatrixFree &,
1423  OutVector &,
1424  const InVector &,
1425  const std::pair<unsigned int, unsigned int> &),
1426  CLASS * owning_class,
1427  OutVector & dst,
1428  const InVector & src,
1429  const bool zero_dst_vector = false,
1430  const DataAccessOnFaces src_vector_face_access =
1432 
1436  template <typename OutVector, typename InVector>
1437  void
1439  const std::function<void(const MatrixFree &,
1440  OutVector &,
1441  const InVector &,
1442  const std::pair<unsigned int, unsigned int> &)>
1443  & cell_operation,
1444  OutVector & dst,
1445  const InVector & src,
1446  const bool zero_dst_vector = false,
1447  const DataAccessOnFaces src_vector_face_access =
1449 
1457  std::pair<unsigned int, unsigned int>
1458  create_cell_subrange_hp(const std::pair<unsigned int, unsigned int> &range,
1459  const unsigned int fe_degree,
1460  const unsigned int dof_handler_index = 0) const;
1461 
1468  std::pair<unsigned int, unsigned int>
1470  const std::pair<unsigned int, unsigned int> &range,
1471  const unsigned int fe_index,
1472  const unsigned int dof_handler_index = 0) const;
1473 
1477  unsigned int
1478  n_active_fe_indices() const;
1479 
1483  unsigned int
1485  const std::pair<unsigned int, unsigned int> range) const;
1486 
1490  unsigned int
1491  get_face_active_fe_index(const std::pair<unsigned int, unsigned int> range,
1492  const bool is_interior_face = true) const;
1493 
1495 
1520  template <typename VectorType>
1521  void
1523  const unsigned int dof_handler_index = 0) const;
1524 
1545  template <typename Number2>
1546  void
1548  const unsigned int dof_handler_index = 0) const;
1549 
1560  const std::shared_ptr<const Utilities::MPI::Partitioner> &
1561  get_vector_partitioner(const unsigned int dof_handler_index = 0) const;
1562 
1566  const IndexSet &
1567  get_locally_owned_set(const unsigned int dof_handler_index = 0) const;
1568 
1572  const IndexSet &
1573  get_ghost_set(const unsigned int dof_handler_index = 0) const;
1574 
1584  const std::vector<unsigned int> &
1585  get_constrained_dofs(const unsigned int dof_handler_index = 0) const;
1586 
1597  void
1598  renumber_dofs(std::vector<types::global_dof_index> &renumbering,
1599  const unsigned int dof_handler_index = 0);
1600 
1602 
1610  template <int spacedim>
1611  static bool
1613 
1617  unsigned int
1618  n_components() const;
1619 
1624  unsigned int
1625  n_base_elements(const unsigned int dof_handler_index) const;
1626 
1634  unsigned int
1635  n_physical_cells() const;
1636 
1640  DEAL_II_DEPRECATED unsigned int
1641  n_macro_cells() const;
1642 
1652  unsigned int
1653  n_cell_batches() const;
1654 
1661  unsigned int
1662  n_ghost_cell_batches() const;
1663 
1671  unsigned int
1672  n_inner_face_batches() const;
1673 
1682  unsigned int
1683  n_boundary_face_batches() const;
1684 
1689  unsigned int
1691 
1699  get_boundary_id(const unsigned int macro_face) const;
1700 
1705  std::array<types::boundary_id, VectorizedArrayType::size()>
1706  get_faces_by_cells_boundary_id(const unsigned int cell_batch_index,
1707  const unsigned int face_number) const;
1708 
1713  const DoFHandler<dim> &
1714  get_dof_handler(const unsigned int dof_handler_index = 0) const;
1715 
1726  template <typename DoFHandlerType>
1727  DEAL_II_DEPRECATED const DoFHandlerType &
1728  get_dof_handler(const unsigned int dof_handler_index = 0) const;
1729 
1743  get_cell_iterator(const unsigned int cell_batch_index,
1744  const unsigned int lane_index,
1745  const unsigned int dof_handler_index = 0) const;
1746 
1752  std::pair<int, int>
1753  get_cell_level_and_index(const unsigned int cell_batch_index,
1754  const unsigned int lane_index) const;
1755 
1768  std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>
1769  get_face_iterator(const unsigned int face_batch_index,
1770  const unsigned int lane_index,
1771  const bool interior = true,
1772  const unsigned int fe_component = 0) const;
1773 
1780  get_hp_cell_iterator(const unsigned int cell_batch_index,
1781  const unsigned int lane_index,
1782  const unsigned int dof_handler_index = 0) const;
1783 
1796  bool
1797  at_irregular_cell(const unsigned int cell_batch_index) const;
1798 
1802  DEAL_II_DEPRECATED unsigned int
1803  n_components_filled(const unsigned int cell_batch_number) const;
1804 
1814  unsigned int
1815  n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const;
1816 
1826  unsigned int
1827  n_active_entries_per_face_batch(const unsigned int face_batch_index) const;
1828 
1832  unsigned int
1833  get_dofs_per_cell(const unsigned int dof_handler_index = 0,
1834  const unsigned int hp_active_fe_index = 0) const;
1835 
1839  unsigned int
1840  get_n_q_points(const unsigned int quad_index = 0,
1841  const unsigned int hp_active_fe_index = 0) const;
1842 
1847  unsigned int
1848  get_dofs_per_face(const unsigned int dof_handler_index = 0,
1849  const unsigned int hp_active_fe_index = 0) const;
1850 
1855  unsigned int
1856  get_n_q_points_face(const unsigned int quad_index = 0,
1857  const unsigned int hp_active_fe_index = 0) const;
1858 
1862  const Quadrature<dim> &
1863  get_quadrature(const unsigned int quad_index = 0,
1864  const unsigned int hp_active_fe_index = 0) const;
1865 
1869  const Quadrature<dim - 1> &
1870  get_face_quadrature(const unsigned int quad_index = 0,
1871  const unsigned int hp_active_fe_index = 0) const;
1872 
1879  unsigned int
1880  get_cell_category(const unsigned int cell_batch_index) const;
1881 
1886  std::pair<unsigned int, unsigned int>
1887  get_face_category(const unsigned int macro_face) const;
1888 
1892  bool
1893  indices_initialized() const;
1894 
1899  bool
1900  mapping_initialized() const;
1901 
1906  unsigned int
1907  get_mg_level() const;
1908 
1913  std::size_t
1914  memory_consumption() const;
1915 
1920  template <typename StreamType>
1921  void
1922  print_memory_consumption(StreamType &out) const;
1923 
1928  void
1929  print(std::ostream &out) const;
1930 
1932 
1943  get_task_info() const;
1944 
1945  /*
1946  * Return geometry-dependent information on the cells.
1947  */
1948  const internal::MatrixFreeFunctions::
1949  MappingInfo<dim, Number, VectorizedArrayType> &
1950  get_mapping_info() const;
1951 
1956  get_dof_info(const unsigned int dof_handler_index_component = 0) const;
1957 
1961  unsigned int
1962  n_constraint_pool_entries() const;
1963 
1968  const Number *
1969  constraint_pool_begin(const unsigned int pool_index) const;
1970 
1976  const Number *
1977  constraint_pool_end(const unsigned int pool_index) const;
1978 
1983  get_shape_info(const unsigned int dof_handler_index_component = 0,
1984  const unsigned int quad_index = 0,
1985  const unsigned int fe_base_element = 0,
1986  const unsigned int hp_active_fe_index = 0,
1987  const unsigned int hp_active_quad_index = 0) const;
1988 
1993  VectorizedArrayType::size()> &
1994  get_face_info(const unsigned int face_batch_index) const;
1995 
1996 
2002  const Table<3, unsigned int> &
2004 
2019  acquire_scratch_data() const;
2020 
2024  void
2026 
2038 
2042  void
2044  const AlignedVector<Number> *memory) const;
2045 
2047 
2048 private:
2053  template <typename number2, int q_dim>
2054  void
2056  const std::shared_ptr<hp::MappingCollection<dim>> & mapping,
2057  const std::vector<const DoFHandler<dim, dim> *> & dof_handlers,
2058  const std::vector<const AffineConstraints<number2> *> &constraint,
2059  const std::vector<IndexSet> & locally_owned_set,
2060  const std::vector<hp::QCollection<q_dim>> & quad,
2061  const AdditionalData & additional_data);
2062 
2069  template <typename number2>
2070  void
2072  const std::vector<const AffineConstraints<number2> *> &constraint,
2073  const std::vector<IndexSet> & locally_owned_set,
2074  const AdditionalData & additional_data);
2075 
2079  void
2081  const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
2082  const AdditionalData & additional_data);
2083 
2087  std::vector<SmartPointer<const DoFHandler<dim>>> dof_handlers;
2088 
2093  std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
2094 
2101  std::vector<Number> constraint_pool_data;
2102 
2107  std::vector<unsigned int> constraint_pool_row_index;
2108 
2115 
2121 
2128  std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
2129 
2130 
2138 
2145 
2152 
2157 
2162 
2171  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2173 
2178  mutable std::list<std::pair<bool, AlignedVector<Number>>>
2180 
2184  unsigned int mg_level;
2185 };
2186 
2187 
2188 
2189 /*----------------------- Inline functions ----------------------------------*/
2190 
2191 #ifndef DOXYGEN
2192 
2193 
2194 
2195 template <int dim, typename Number, typename VectorizedArrayType>
2196 template <typename VectorType>
2197 inline void
2199  VectorType & vec,
2200  const unsigned int comp) const
2201 {
2202  AssertIndexRange(comp, n_components());
2203  vec.reinit(dof_info[comp].vector_partitioner->size());
2204 }
2205 
2206 
2207 
2208 template <int dim, typename Number, typename VectorizedArrayType>
2209 template <typename Number2>
2210 inline void
2213  const unsigned int comp) const
2214 {
2215  AssertIndexRange(comp, n_components());
2216  vec.reinit(dof_info[comp].vector_partitioner, task_info.communicator_sm);
2217 }
2218 
2219 
2220 
2221 template <int dim, typename Number, typename VectorizedArrayType>
2222 inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2224  const unsigned int comp) const
2225 {
2226  AssertIndexRange(comp, n_components());
2227  return dof_info[comp].vector_partitioner;
2228 }
2229 
2230 
2231 
2232 template <int dim, typename Number, typename VectorizedArrayType>
2233 inline const std::vector<unsigned int> &
2235  const unsigned int comp) const
2236 {
2237  AssertIndexRange(comp, n_components());
2238  return dof_info[comp].constrained_dofs;
2239 }
2240 
2241 
2242 
2243 template <int dim, typename Number, typename VectorizedArrayType>
2244 inline unsigned int
2246 {
2247  AssertDimension(dof_handlers.size(), dof_info.size());
2248  return dof_handlers.size();
2249 }
2250 
2251 
2252 
2253 template <int dim, typename Number, typename VectorizedArrayType>
2254 inline unsigned int
2256  const unsigned int dof_no) const
2257 {
2258  AssertDimension(dof_handlers.size(), dof_info.size());
2259  AssertIndexRange(dof_no, dof_handlers.size());
2260  return dof_handlers[dof_no]->get_fe().n_base_elements();
2261 }
2262 
2263 
2264 
2265 template <int dim, typename Number, typename VectorizedArrayType>
2268 {
2269  return task_info;
2270 }
2271 
2272 
2273 
2274 template <int dim, typename Number, typename VectorizedArrayType>
2275 inline unsigned int
2277 {
2278  return *(task_info.cell_partition_data.end() - 2);
2279 }
2280 
2281 
2282 
2283 template <int dim, typename Number, typename VectorizedArrayType>
2284 inline unsigned int
2286 {
2287  return task_info.n_active_cells;
2288 }
2289 
2290 
2291 
2292 template <int dim, typename Number, typename VectorizedArrayType>
2293 inline unsigned int
2295 {
2296  return *(task_info.cell_partition_data.end() - 2);
2297 }
2298 
2299 
2300 
2301 template <int dim, typename Number, typename VectorizedArrayType>
2302 inline unsigned int
2304 {
2305  return *(task_info.cell_partition_data.end() - 1) -
2306  *(task_info.cell_partition_data.end() - 2);
2307 }
2308 
2309 
2310 
2311 template <int dim, typename Number, typename VectorizedArrayType>
2312 inline unsigned int
2314 {
2315  if (task_info.face_partition_data.size() == 0)
2316  return 0;
2317  return task_info.face_partition_data.back();
2318 }
2319 
2320 
2321 
2322 template <int dim, typename Number, typename VectorizedArrayType>
2323 inline unsigned int
2325 {
2326  if (task_info.face_partition_data.size() == 0)
2327  return 0;
2328  return task_info.boundary_partition_data.back() -
2330 }
2331 
2332 
2333 
2334 template <int dim, typename Number, typename VectorizedArrayType>
2335 inline unsigned int
2337 {
2338  if (task_info.face_partition_data.size() == 0)
2339  return 0;
2340  return face_info.faces.size() - task_info.boundary_partition_data.back();
2341 }
2342 
2343 
2344 
2345 template <int dim, typename Number, typename VectorizedArrayType>
2346 inline types::boundary_id
2348  const unsigned int macro_face) const
2349 {
2350  Assert(macro_face >= task_info.boundary_partition_data[0] &&
2351  macro_face < task_info.boundary_partition_data.back(),
2352  ExcIndexRange(macro_face,
2355  return types::boundary_id(face_info.faces[macro_face].exterior_face_no);
2356 }
2357 
2358 
2359 
2360 template <int dim, typename Number, typename VectorizedArrayType>
2361 inline std::array<types::boundary_id, VectorizedArrayType::size()>
2363  const unsigned int cell_batch_index,
2364  const unsigned int face_number) const
2365 {
2366  AssertIndexRange(cell_batch_index, n_cell_batches());
2369  ExcNotInitialized());
2370  std::array<types::boundary_id, VectorizedArrayType::size()> result;
2371  result.fill(numbers::invalid_boundary_id);
2372  for (unsigned int v = 0;
2373  v < n_active_entries_per_cell_batch(cell_batch_index);
2374  ++v)
2375  result[v] =
2376  face_info.cell_and_face_boundary_id(cell_batch_index, face_number, v);
2377  return result;
2378 }
2379 
2380 
2381 
2382 template <int dim, typename Number, typename VectorizedArrayType>
2383 inline const internal::MatrixFreeFunctions::
2384  MappingInfo<dim, Number, VectorizedArrayType> &
2386 {
2387  return mapping_info;
2388 }
2389 
2390 
2391 
2392 template <int dim, typename Number, typename VectorizedArrayType>
2395  const unsigned int dof_index) const
2396 {
2397  AssertIndexRange(dof_index, n_components());
2398  return dof_info[dof_index];
2399 }
2400 
2401 
2402 
2403 template <int dim, typename Number, typename VectorizedArrayType>
2404 inline unsigned int
2406 {
2407  return constraint_pool_row_index.size() - 1;
2408 }
2409 
2410 
2411 
2412 template <int dim, typename Number, typename VectorizedArrayType>
2413 inline const Number *
2415  const unsigned int row) const
2416 {
2418  return constraint_pool_data.empty() ?
2419  nullptr :
2421 }
2422 
2423 
2424 
2425 template <int dim, typename Number, typename VectorizedArrayType>
2426 inline const Number *
2428  const unsigned int row) const
2429 {
2431  return constraint_pool_data.empty() ?
2432  nullptr :
2434 }
2435 
2436 
2437 
2438 template <int dim, typename Number, typename VectorizedArrayType>
2439 inline std::pair<unsigned int, unsigned int>
2441  const std::pair<unsigned int, unsigned int> &range,
2442  const unsigned int degree,
2443  const unsigned int dof_handler_component) const
2444 {
2445  if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2446  {
2448  dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2450  dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2451  if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2452  return range;
2453  else
2454  return {range.second, range.second};
2455  }
2456 
2457  const unsigned int fe_index =
2458  dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2459  if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2460  return {range.second, range.second};
2461  else
2462  return create_cell_subrange_hp_by_index(range,
2463  fe_index,
2464  dof_handler_component);
2465 }
2466 
2467 
2468 
2469 template <int dim, typename Number, typename VectorizedArrayType>
2470 inline bool
2472  const unsigned int cell_batch_index) const
2473 {
2474  AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2475  return VectorizedArrayType::size() > 1 &&
2476  cell_level_index[(cell_batch_index + 1) * VectorizedArrayType::size() -
2477  1] == cell_level_index[(cell_batch_index + 1) *
2478  VectorizedArrayType::size() -
2479  2];
2480 }
2481 
2482 
2483 
2484 template <int dim, typename Number, typename VectorizedArrayType>
2485 unsigned int
2487 {
2488  return shape_info.size(2);
2489 }
2490 
2491 
2492 template <int dim, typename Number, typename VectorizedArrayType>
2493 unsigned int
2495  const std::pair<unsigned int, unsigned int> range) const
2496 {
2497  const auto &fe_indices = dof_info[0].cell_active_fe_index;
2498 
2499  if (fe_indices.empty() == true)
2500  return 0;
2501 
2502  const auto index = fe_indices[range.first];
2503 
2504  for (unsigned int i = range.first; i < range.second; ++i)
2505  AssertDimension(index, fe_indices[i]);
2506 
2507  return index;
2508 }
2509 
2510 
2511 
2512 template <int dim, typename Number, typename VectorizedArrayType>
2513 unsigned int
2515  const std::pair<unsigned int, unsigned int> range,
2516  const bool is_interior_face) const
2517 {
2518  const auto &fe_indices = dof_info[0].cell_active_fe_index;
2519 
2520  if (fe_indices.empty() == true)
2521  return 0;
2522 
2523  if (is_interior_face)
2524  {
2525  const unsigned int index =
2526  fe_indices[face_info.faces[range.first].cells_interior[0] /
2527  VectorizedArrayType::size()];
2528 
2529  for (unsigned int i = range.first; i < range.second; ++i)
2530  AssertDimension(index,
2531  fe_indices[face_info.faces[i].cells_interior[0] /
2532  VectorizedArrayType::size()]);
2533 
2534  return index;
2535  }
2536  else
2537  {
2538  const unsigned int index =
2539  fe_indices[face_info.faces[range.first].cells_exterior[0] /
2540  VectorizedArrayType::size()];
2541 
2542  for (unsigned int i = range.first; i < range.second; ++i)
2543  AssertDimension(index,
2544  fe_indices[face_info.faces[i].cells_exterior[0] /
2545  VectorizedArrayType::size()]);
2546 
2547  return index;
2548  }
2549 }
2550 
2551 
2552 
2553 template <int dim, typename Number, typename VectorizedArrayType>
2554 inline unsigned int
2556  const unsigned int cell_batch_index) const
2557 {
2558  return n_active_entries_per_cell_batch(cell_batch_index);
2559 }
2560 
2561 
2562 
2563 template <int dim, typename Number, typename VectorizedArrayType>
2564 inline unsigned int
2566  const unsigned int cell_batch_index) const
2567 {
2568  AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2569  unsigned int n_lanes = VectorizedArrayType::size();
2570  while (n_lanes > 1 &&
2571  cell_level_index[cell_batch_index * VectorizedArrayType::size() +
2572  n_lanes - 1] ==
2573  cell_level_index[cell_batch_index * VectorizedArrayType::size() +
2574  n_lanes - 2])
2575  --n_lanes;
2576  AssertIndexRange(n_lanes - 1, VectorizedArrayType::size());
2577  return n_lanes;
2578 }
2579 
2580 
2581 
2582 template <int dim, typename Number, typename VectorizedArrayType>
2583 inline unsigned int
2585  const unsigned int face_batch_index) const
2586 {
2587  AssertIndexRange(face_batch_index, face_info.faces.size());
2588  unsigned int n_lanes = VectorizedArrayType::size();
2589  while (n_lanes > 1 &&
2590  face_info.faces[face_batch_index].cells_interior[n_lanes - 1] ==
2592  --n_lanes;
2593  AssertIndexRange(n_lanes - 1, VectorizedArrayType::size());
2594  return n_lanes;
2595 }
2596 
2597 
2598 
2599 template <int dim, typename Number, typename VectorizedArrayType>
2600 inline unsigned int
2602  const unsigned int dof_handler_index,
2603  const unsigned int active_fe_index) const
2604 {
2605  return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2606 }
2607 
2608 
2609 
2610 template <int dim, typename Number, typename VectorizedArrayType>
2611 inline unsigned int
2613  const unsigned int quad_index,
2614  const unsigned int active_fe_index) const
2615 {
2616  AssertIndexRange(quad_index, mapping_info.cell_data.size());
2617  return mapping_info.cell_data[quad_index]
2618  .descriptor[active_fe_index]
2619  .n_q_points;
2620 }
2621 
2622 
2623 
2624 template <int dim, typename Number, typename VectorizedArrayType>
2625 inline unsigned int
2627  const unsigned int dof_handler_index,
2628  const unsigned int active_fe_index) const
2629 {
2630  return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2631 }
2632 
2633 
2634 
2635 template <int dim, typename Number, typename VectorizedArrayType>
2636 inline unsigned int
2638  const unsigned int quad_index,
2639  const unsigned int active_fe_index) const
2640 {
2641  AssertIndexRange(quad_index, mapping_info.face_data.size());
2642  return mapping_info.face_data[quad_index]
2643  .descriptor[active_fe_index]
2644  .n_q_points;
2645 }
2646 
2647 
2648 
2649 template <int dim, typename Number, typename VectorizedArrayType>
2650 inline const IndexSet &
2652  const unsigned int dof_handler_index) const
2653 {
2654  return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2655 }
2656 
2657 
2658 
2659 template <int dim, typename Number, typename VectorizedArrayType>
2660 inline const IndexSet &
2662  const unsigned int dof_handler_index) const
2663 {
2664  return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2665 }
2666 
2667 
2668 
2669 template <int dim, typename Number, typename VectorizedArrayType>
2672  const unsigned int dof_handler_index,
2673  const unsigned int index_quad,
2674  const unsigned int index_fe,
2675  const unsigned int active_fe_index,
2676  const unsigned int active_quad_index) const
2677 {
2678  AssertIndexRange(dof_handler_index, dof_info.size());
2679  const unsigned int ind =
2680  dof_info[dof_handler_index].global_base_element_offset + index_fe;
2681  AssertIndexRange(ind, shape_info.size(0));
2682  AssertIndexRange(index_quad, shape_info.size(1));
2683  AssertIndexRange(active_fe_index, shape_info.size(2));
2684  AssertIndexRange(active_quad_index, shape_info.size(3));
2685  return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2686 }
2687 
2688 
2689 
2690 template <int dim, typename Number, typename VectorizedArrayType>
2692  VectorizedArrayType::size()> &
2694  const unsigned int macro_face) const
2695 {
2696  AssertIndexRange(macro_face, face_info.faces.size());
2697  return face_info.faces[macro_face];
2698 }
2699 
2700 
2701 
2702 template <int dim, typename Number, typename VectorizedArrayType>
2703 inline const Table<3, unsigned int> &
2705  const
2706 {
2708 }
2709 
2710 
2711 
2712 template <int dim, typename Number, typename VectorizedArrayType>
2713 inline const Quadrature<dim> &
2715  const unsigned int quad_index,
2716  const unsigned int active_fe_index) const
2717 {
2718  AssertIndexRange(quad_index, mapping_info.cell_data.size());
2719  return mapping_info.cell_data[quad_index]
2720  .descriptor[active_fe_index]
2721  .quadrature;
2722 }
2723 
2724 
2725 
2726 template <int dim, typename Number, typename VectorizedArrayType>
2727 inline const Quadrature<dim - 1> &
2729  const unsigned int quad_index,
2730  const unsigned int active_fe_index) const
2731 {
2732  AssertIndexRange(quad_index, mapping_info.face_data.size());
2733  return mapping_info.face_data[quad_index]
2734  .descriptor[active_fe_index]
2735  .quadrature;
2736 }
2737 
2738 
2739 
2740 template <int dim, typename Number, typename VectorizedArrayType>
2741 inline unsigned int
2743  const unsigned int cell_batch_index) const
2744 {
2745  AssertIndexRange(0, dof_info.size());
2746  AssertIndexRange(cell_batch_index, dof_info[0].cell_active_fe_index.size());
2747  if (dof_info[0].cell_active_fe_index.empty())
2748  return 0;
2749  else
2750  return dof_info[0].cell_active_fe_index[cell_batch_index];
2751 }
2752 
2753 
2754 
2755 template <int dim, typename Number, typename VectorizedArrayType>
2756 inline std::pair<unsigned int, unsigned int>
2758  const unsigned int macro_face) const
2759 {
2760  AssertIndexRange(macro_face, face_info.faces.size());
2761  if (dof_info[0].cell_active_fe_index.empty())
2762  return std::make_pair(0U, 0U);
2763 
2764  std::pair<unsigned int, unsigned int> result;
2765  for (unsigned int v = 0; v < VectorizedArrayType::size() &&
2766  face_info.faces[macro_face].cells_interior[v] !=
2768  ++v)
2769  result.first = std::max(
2770  result.first,
2771  dof_info[0]
2772  .cell_active_fe_index[face_info.faces[macro_face].cells_interior[v]]);
2773  if (face_info.faces[macro_face].cells_exterior[0] !=
2775  for (unsigned int v = 0; v < VectorizedArrayType::size() &&
2776  face_info.faces[macro_face].cells_exterior[v] !=
2778  ++v)
2779  result.second = std::max(
2780  result.first,
2781  dof_info[0]
2782  .cell_active_fe_index[face_info.faces[macro_face].cells_exterior[v]]);
2783  else
2784  result.second = numbers::invalid_unsigned_int;
2785  return result;
2786 }
2787 
2788 
2789 
2790 template <int dim, typename Number, typename VectorizedArrayType>
2791 inline bool
2793 {
2794  return indices_are_initialized;
2795 }
2796 
2797 
2798 
2799 template <int dim, typename Number, typename VectorizedArrayType>
2800 inline bool
2802 {
2803  return mapping_is_initialized;
2804 }
2805 
2806 
2807 template <int dim, typename Number, typename VectorizedArrayType>
2808 inline unsigned int
2810 {
2811  return mg_level;
2812 }
2813 
2814 
2815 
2816 template <int dim, typename Number, typename VectorizedArrayType>
2819 {
2820  using list_type =
2821  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
2822  list_type &data = scratch_pad.get();
2823  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
2824  if (it->first == false)
2825  {
2826  it->first = true;
2827  return &it->second;
2828  }
2829  data.emplace_front(true, AlignedVector<VectorizedArrayType>());
2830  return &data.front().second;
2831 }
2832 
2833 
2834 
2835 template <int dim, typename Number, typename VectorizedArrayType>
2836 void
2838  const AlignedVector<VectorizedArrayType> *scratch) const
2839 {
2840  using list_type =
2841  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
2842  list_type &data = scratch_pad.get();
2843  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
2844  if (&it->second == scratch)
2845  {
2846  Assert(it->first == true, ExcInternalError());
2847  it->first = false;
2848  return;
2849  }
2850  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2851 }
2852 
2853 
2854 
2855 template <int dim, typename Number, typename VectorizedArrayType>
2859 {
2860  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
2862  it != scratch_pad_non_threadsafe.end();
2863  ++it)
2864  if (it->first == false)
2865  {
2866  it->first = true;
2867  return &it->second;
2868  }
2869  scratch_pad_non_threadsafe.push_front(
2870  std::make_pair(true, AlignedVector<Number>()));
2871  return &scratch_pad_non_threadsafe.front().second;
2872 }
2873 
2874 
2875 
2876 template <int dim, typename Number, typename VectorizedArrayType>
2877 void
2880  const AlignedVector<Number> *scratch) const
2881 {
2882  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
2884  it != scratch_pad_non_threadsafe.end();
2885  ++it)
2886  if (&it->second == scratch)
2887  {
2888  Assert(it->first == true, ExcInternalError());
2889  it->first = false;
2890  return;
2891  }
2892  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2893 }
2894 
2895 
2896 
2897 // ------------------------------ reinit functions ---------------------------
2898 
2899 namespace internal
2900 {
2901  namespace MatrixFreeImplementation
2902  {
2903  template <int dim, int spacedim>
2904  inline std::vector<IndexSet>
2905  extract_locally_owned_index_sets(
2906  const std::vector<const ::DoFHandler<dim, spacedim> *> &dofh,
2907  const unsigned int level)
2908  {
2909  std::vector<IndexSet> locally_owned_set;
2910  locally_owned_set.reserve(dofh.size());
2911  for (unsigned int j = 0; j < dofh.size(); j++)
2912  if (level == numbers::invalid_unsigned_int)
2913  locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
2914  else
2915  locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
2916  return locally_owned_set;
2917  }
2918  } // namespace MatrixFreeImplementation
2919 } // namespace internal
2920 
2921 
2922 
2923 template <int dim, typename Number, typename VectorizedArrayType>
2924 template <typename QuadratureType, typename number2>
2925 void
2927  const DoFHandler<dim> & dof_handler,
2928  const AffineConstraints<number2> &constraints_in,
2929  const QuadratureType & quad,
2931  &additional_data)
2932 {
2933  std::vector<const DoFHandler<dim, dim> *> dof_handlers;
2934  std::vector<const AffineConstraints<number2> *> constraints;
2935  std::vector<QuadratureType> quads;
2936 
2937  dof_handlers.push_back(&dof_handler);
2938  constraints.push_back(&constraints_in);
2939  quads.push_back(quad);
2940 
2941  std::vector<IndexSet> locally_owned_sets =
2942  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2943  dof_handlers, additional_data.mg_level);
2944 
2945  std::vector<hp::QCollection<dim>> quad_hp;
2946  quad_hp.emplace_back(quad);
2947 
2948  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(
2950  dof_handlers,
2951  constraints,
2952  locally_owned_sets,
2953  quad_hp,
2954  additional_data);
2955 }
2956 
2957 
2958 
2959 template <int dim, typename Number, typename VectorizedArrayType>
2960 template <typename QuadratureType, typename number2, typename MappingType>
2961 void
2963  const MappingType & mapping,
2964  const DoFHandler<dim> & dof_handler,
2965  const AffineConstraints<number2> &constraints_in,
2966  const QuadratureType & quad,
2968  &additional_data)
2969 {
2970  std::vector<const DoFHandler<dim, dim> *> dof_handlers;
2971  std::vector<const AffineConstraints<number2> *> constraints;
2972 
2973  dof_handlers.push_back(&dof_handler);
2974  constraints.push_back(&constraints_in);
2975 
2976  std::vector<IndexSet> locally_owned_sets =
2977  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2978  dof_handlers, additional_data.mg_level);
2979 
2980  std::vector<hp::QCollection<dim>> quad_hp;
2981  quad_hp.emplace_back(quad);
2982 
2983  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
2984  dof_handlers,
2985  constraints,
2986  locally_owned_sets,
2987  quad_hp,
2988  additional_data);
2989 }
2990 
2991 
2992 
2993 template <int dim, typename Number, typename VectorizedArrayType>
2994 template <typename QuadratureType, typename number2>
2995 void
2997  const std::vector<const DoFHandler<dim> *> & dof_handler,
2998  const std::vector<const AffineConstraints<number2> *> &constraint,
2999  const std::vector<QuadratureType> & quad,
3001  &additional_data)
3002 {
3003  std::vector<IndexSet> locally_owned_set =
3004  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3005  dof_handler, additional_data.mg_level);
3006  std::vector<hp::QCollection<dim>> quad_hp;
3007  for (unsigned int q = 0; q < quad.size(); ++q)
3008  quad_hp.emplace_back(quad[q]);
3009 
3010  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(
3012  dof_handler,
3013  constraint,
3014  locally_owned_set,
3015  quad_hp,
3016  additional_data);
3017 }
3018 
3019 
3020 
3021 template <int dim, typename Number, typename VectorizedArrayType>
3022 template <typename QuadratureType,
3023  typename number2,
3024  typename DoFHandlerType,
3025  typename MappingType>
3026 void
3028  const MappingType & mapping,
3029  const std::vector<const DoFHandlerType *> & dof_handler,
3030  const std::vector<const AffineConstraints<number2> *> &constraint,
3031  const std::vector<QuadratureType> & quad,
3032  const AdditionalData & additional_data)
3033 {
3034  static_assert(dim == DoFHandlerType::dimension,
3035  "Dimension dim not equal to DoFHandlerType::dimension.");
3036 
3037  std::vector<const DoFHandler<dim> *> dof_handlers;
3038 
3039  for (const auto dh : dof_handler)
3040  dof_handlers.push_back(dh);
3041 
3042  this->reinit(mapping, dof_handlers, constraint, quad, additional_data);
3043 }
3044 
3045 
3046 
3047 template <int dim, typename Number, typename VectorizedArrayType>
3048 template <typename QuadratureType, typename number2>
3049 void
3051  const std::vector<const DoFHandler<dim> *> & dof_handler,
3052  const std::vector<const AffineConstraints<number2> *> &constraint,
3053  const QuadratureType & quad,
3055  &additional_data)
3056 {
3057  std::vector<IndexSet> locally_owned_set =
3058  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3059  dof_handler, additional_data.mg_level);
3060  std::vector<hp::QCollection<dim>> quad_hp;
3061  quad_hp.emplace_back(quad);
3062 
3063  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(
3065  dof_handler,
3066  constraint,
3067  locally_owned_set,
3068  quad_hp,
3069  additional_data);
3070 }
3071 
3072 
3073 
3074 template <int dim, typename Number, typename VectorizedArrayType>
3075 template <typename QuadratureType, typename number2, typename DoFHandlerType>
3076 void
3078  const std::vector<const DoFHandlerType *> & dof_handler,
3079  const std::vector<const AffineConstraints<number2> *> &constraint,
3080  const std::vector<QuadratureType> & quad,
3081  const AdditionalData & additional_data)
3082 {
3083  static_assert(dim == DoFHandlerType::dimension,
3084  "Dimension dim not equal to DoFHandlerType::dimension.");
3085 
3086  std::vector<const DoFHandler<dim> *> dof_handlers;
3087 
3088  for (const auto dh : dof_handler)
3089  dof_handlers.push_back(dof_handler);
3090 
3091  this->reinit(dof_handlers, constraint, quad, additional_data);
3092 }
3093 
3094 
3095 
3096 template <int dim, typename Number, typename VectorizedArrayType>
3097 template <typename QuadratureType, typename number2, typename MappingType>
3098 void
3100  const MappingType & mapping,
3101  const std::vector<const DoFHandler<dim> *> & dof_handler,
3102  const std::vector<const AffineConstraints<number2> *> &constraint,
3103  const QuadratureType & quad,
3105  &additional_data)
3106 {
3107  std::vector<IndexSet> locally_owned_set =
3108  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3109  dof_handler, additional_data.mg_level);
3110  std::vector<hp::QCollection<dim>> quad_hp;
3111  quad_hp.emplace_back(quad);
3112 
3113  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3114  dof_handler,
3115  constraint,
3116  locally_owned_set,
3117  quad_hp,
3118  additional_data);
3119 }
3120 
3121 
3122 
3123 template <int dim, typename Number, typename VectorizedArrayType>
3124 template <typename QuadratureType, typename number2, typename MappingType>
3125 void
3127  const MappingType & mapping,
3128  const std::vector<const DoFHandler<dim> *> & dof_handler,
3129  const std::vector<const AffineConstraints<number2> *> &constraint,
3130  const std::vector<QuadratureType> & quad,
3132  &additional_data)
3133 {
3134  std::vector<IndexSet> locally_owned_set =
3135  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3136  dof_handler, additional_data.mg_level);
3137  std::vector<hp::QCollection<dim>> quad_hp;
3138  for (unsigned int q = 0; q < quad.size(); ++q)
3139  quad_hp.emplace_back(quad[q]);
3140 
3141  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3142  dof_handler,
3143  constraint,
3144  locally_owned_set,
3145  quad_hp,
3146  additional_data);
3147 }
3148 
3149 
3150 
3151 template <int dim, typename Number, typename VectorizedArrayType>
3152 template <typename QuadratureType,
3153  typename number2,
3154  typename DoFHandlerType,
3155  typename MappingType>
3156 void
3158  const MappingType & mapping,
3159  const std::vector<const DoFHandlerType *> & dof_handler,
3160  const std::vector<const AffineConstraints<number2> *> &constraint,
3161  const QuadratureType & quad,
3162  const AdditionalData & additional_data)
3163 {
3164  static_assert(dim == DoFHandlerType::dimension,
3165  "Dimension dim not equal to DoFHandlerType::dimension.");
3166 
3167  std::vector<const DoFHandler<dim> *> dof_handlers;
3168 
3169  for (const auto dh : dof_handler)
3170  dof_handlers.push_back(dof_handler);
3171 
3172  this->reinit(mapping, dof_handlers, constraint, quad, additional_data);
3173 }
3174 
3175 
3176 
3177 template <int dim, typename Number, typename VectorizedArrayType>
3178 template <typename QuadratureType, typename number2, typename DoFHandlerType>
3179 void
3181  const std::vector<const DoFHandlerType *> & dof_handler,
3182  const std::vector<const AffineConstraints<number2> *> &constraint,
3183  const QuadratureType & quad,
3184  const AdditionalData & additional_data)
3185 {
3186  static_assert(dim == DoFHandlerType::dimension,
3187  "Dimension dim not equal to DoFHandlerType::dimension.");
3188 
3189  std::vector<const DoFHandler<dim> *> dof_handlers;
3190 
3191  for (const auto dh : dof_handler)
3192  dof_handlers.push_back(dof_handler);
3193 
3194  this->reinit(dof_handlers, constraint, quad, additional_data);
3195 }
3196 
3197 
3198 
3199 // ------------------------------ implementation of loops --------------------
3200 
3201 // internal helper functions that define how to call MPI data exchange
3202 // functions: for generic vectors, do nothing at all. For distributed vectors,
3203 // call update_ghost_values_start function and so on. If we have collections
3204 // of vectors, just do the individual functions of the components. In order to
3205 // keep ghost values consistent (whether we are in read or write mode), we
3206 // also reset the values at the end. the whole situation is a bit complicated
3207 // by the fact that we need to treat block vectors differently, which use some
3208 // additional helper functions to select the blocks and template magic.
3209 namespace internal
3210 {
3214  template <int dim, typename Number, typename VectorizedArrayType>
3215  struct VectorDataExchange
3216  {
3217  // A shift for the MPI messages to reduce the risk for accidental
3218  // interaction with other open communications that a user program might
3219  // set up (parallel vectors support unfinished communication). We let
3220  // the other vectors use the first 20 assigned numbers and start the
3221  // matrix-free communication.
3222  static constexpr unsigned int channel_shift = 20;
3223 
3224 
3225 
3230  VectorDataExchange(
3231  const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3232  const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3233  DataAccessOnFaces vector_face_access,
3234  const unsigned int n_components)
3235  : matrix_free(matrix_free)
3236  , vector_face_access(
3237  matrix_free.get_task_info().face_partition_data.empty() ?
3239  DataAccessOnFaces::unspecified :
3240  vector_face_access)
3241  , ghosts_were_set(false)
3242 # ifdef DEAL_II_WITH_MPI
3243  , tmp_data(n_components)
3244  , requests(n_components)
3245 # endif
3246  {
3247  (void)n_components;
3248  if (this->vector_face_access !=
3251  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3253  matrix_free.get_dof_info(c).vector_exchanger_face_variants.size(),
3254  5);
3255  }
3256 
3257 
3258 
3262  ~VectorDataExchange() // NOLINT
3263  {
3264 # ifdef DEAL_II_WITH_MPI
3265  for (unsigned int i = 0; i < tmp_data.size(); ++i)
3266  if (tmp_data[i] != nullptr)
3267  matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3268 # endif
3269  }
3270 
3271 
3272 
3277  template <typename VectorType>
3278  unsigned int
3279  find_vector_in_mf(const VectorType &vec,
3280  const bool check_global_compatibility = true) const
3281  {
3282  // case 1: vector was set up with MatrixFree::initialize_dof_vector()
3283  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3284  if (vec.get_partitioner().get() ==
3285  matrix_free.get_dof_info(c).vector_partitioner.get())
3286  return c;
3287 
3288  // case 2: user provided own partitioner (compatibility mode)
3289  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3290  if (check_global_compatibility ?
3291  vec.get_partitioner()->is_globally_compatible(
3292  *matrix_free.get_dof_info(c).vector_partitioner) :
3293  vec.get_partitioner()->is_compatible(
3294  *matrix_free.get_dof_info(c).vector_partitioner))
3295  return c;
3296 
3297  Assert(false,
3298  ExcNotImplemented("Could not find partitioner that fits vector"));
3299 
3301  }
3302 
3303 
3304 
3310  get_partitioner(const unsigned int mf_component) const
3311  {
3312  AssertDimension(matrix_free.get_dof_info(mf_component)
3313  .vector_exchanger_face_variants.size(),
3314  5);
3315  if (vector_face_access ==
3318  return *matrix_free.get_dof_info(mf_component)
3319  .vector_exchanger_face_variants[0];
3320  else if (vector_face_access ==
3323  return *matrix_free.get_dof_info(mf_component)
3324  .vector_exchanger_face_variants[1];
3325  else if (vector_face_access ==
3328  return *matrix_free.get_dof_info(mf_component)
3329  .vector_exchanger_face_variants[2];
3330  else if (vector_face_access ==
3333  return *matrix_free.get_dof_info(mf_component)
3334  .vector_exchanger_face_variants[3];
3335  else if (vector_face_access ==
3338  return *matrix_free.get_dof_info(mf_component)
3339  .vector_exchanger_face_variants[4];
3340  else
3341  return *matrix_free.get_dof_info(mf_component).vector_exchanger.get();
3342  }
3343 
3344 
3345 
3349  template <typename VectorType,
3350  typename std::enable_if<is_serial_or_dummy<VectorType>::value,
3351  VectorType>::type * = nullptr>
3352  void
3353  update_ghost_values_start(const unsigned int /*component_in_block_vector*/,
3354  const VectorType & /*vec*/)
3355  {}
3356 
3357 
3362  template <typename VectorType,
3363  typename std::enable_if<
3364  !has_update_ghost_values_start<VectorType>::value &&
3365  !is_serial_or_dummy<VectorType>::value,
3366  VectorType>::type * = nullptr>
3367  void
3368  update_ghost_values_start(const unsigned int component_in_block_vector,
3369  const VectorType & vec)
3370  {
3371  (void)component_in_block_vector;
3372  bool ghosts_set = vec.has_ghost_elements();
3373  if (ghosts_set)
3374  ghosts_were_set = true;
3375 
3376  vec.update_ghost_values();
3377  }
3378 
3379 
3380 
3386  template <typename VectorType,
3387  typename std::enable_if<
3388  has_update_ghost_values_start<VectorType>::value &&
3389  !has_exchange_on_subset<VectorType>::value,
3390  VectorType>::type * = nullptr>
3391  void
3392  update_ghost_values_start(const unsigned int component_in_block_vector,
3393  const VectorType & vec)
3394  {
3395  (void)component_in_block_vector;
3396  bool ghosts_set = vec.has_ghost_elements();
3397  if (ghosts_set)
3398  ghosts_were_set = true;
3399 
3400  vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3401  }
3402 
3403 
3404 
3411  template <typename VectorType,
3412  typename std::enable_if<
3413  has_update_ghost_values_start<VectorType>::value &&
3414  has_exchange_on_subset<VectorType>::value,
3415  VectorType>::type * = nullptr>
3416  void
3417  update_ghost_values_start(const unsigned int component_in_block_vector,
3418  const VectorType & vec)
3419  {
3420  static_assert(
3421  std::is_same<Number, typename VectorType::value_type>::value,
3422  "Type mismatch between VectorType and VectorDataExchange");
3423  (void)component_in_block_vector;
3424  bool ghosts_set = vec.has_ghost_elements();
3425  if (ghosts_set)
3426  ghosts_were_set = true;
3427 
3428  if (vec.size() != 0)
3429  {
3430 # ifdef DEAL_II_WITH_MPI
3431  const unsigned int mf_component = find_vector_in_mf(vec);
3432 
3433  const auto &part = get_partitioner(mf_component);
3434 
3435  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3436  part.n_import_sm_procs() == 0)
3437  return;
3438 
3439  tmp_data[component_in_block_vector] =
3440  matrix_free.acquire_scratch_data_non_threadsafe();
3441  tmp_data[component_in_block_vector]->resize_fast(
3442  part.n_import_indices());
3443  AssertDimension(requests.size(), tmp_data.size());
3444 
3445  part.export_to_ghosted_array_start(
3446  component_in_block_vector * 2 + channel_shift,
3447  ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3448  vec.shared_vector_data(),
3449  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3450  part.locally_owned_size(),
3451  matrix_free.get_dof_info(mf_component)
3452  .vector_partitioner->n_ghost_indices()),
3453  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3454  part.n_import_indices()),
3455  this->requests[component_in_block_vector]);
3456 # endif
3457  }
3458  }
3459 
3460 
3461 
3466  template <
3467  typename VectorType,
3468  typename std::enable_if<!has_update_ghost_values_start<VectorType>::value,
3469  VectorType>::type * = nullptr>
3470  void
3471  update_ghost_values_finish(const unsigned int /*component_in_block_vector*/,
3472  const VectorType & /*vec*/)
3473  {}
3474 
3475 
3476 
3482  template <typename VectorType,
3483  typename std::enable_if<
3484  has_update_ghost_values_start<VectorType>::value &&
3485  !has_exchange_on_subset<VectorType>::value,
3486  VectorType>::type * = nullptr>
3487  void
3488  update_ghost_values_finish(const unsigned int component_in_block_vector,
3489  const VectorType & vec)
3490  {
3491  (void)component_in_block_vector;
3492  vec.update_ghost_values_finish();
3493  }
3494 
3495 
3496 
3503  template <typename VectorType,
3504  typename std::enable_if<
3505  has_update_ghost_values_start<VectorType>::value &&
3506  has_exchange_on_subset<VectorType>::value,
3507  VectorType>::type * = nullptr>
3508  void
3509  update_ghost_values_finish(const unsigned int component_in_block_vector,
3510  const VectorType & vec)
3511  {
3512  static_assert(
3513  std::is_same<Number, typename VectorType::value_type>::value,
3514  "Type mismatch between VectorType and VectorDataExchange");
3515  (void)component_in_block_vector;
3516 
3517  if (vec.size() != 0)
3518  {
3519 # ifdef DEAL_II_WITH_MPI
3520  AssertIndexRange(component_in_block_vector, tmp_data.size());
3521  AssertDimension(requests.size(), tmp_data.size());
3522 
3523  const unsigned int mf_component = find_vector_in_mf(vec);
3524 
3525  const auto &part = get_partitioner(mf_component);
3526 
3527  if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3528  part.n_import_sm_procs() != 0)
3529  {
3530  part.export_to_ghosted_array_finish(
3531  ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3532  vec.shared_vector_data(),
3533  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3534  part.locally_owned_size(),
3535  matrix_free.get_dof_info(mf_component)
3536  .vector_partitioner->n_ghost_indices()),
3537  this->requests[component_in_block_vector]);
3538 
3539  matrix_free.release_scratch_data_non_threadsafe(
3540  tmp_data[component_in_block_vector]);
3541  tmp_data[component_in_block_vector] = nullptr;
3542  }
3543 # endif
3544  }
3545  // let vector know that ghosts are being updated and we can read from
3546  // them
3547  vec.set_ghost_state(true);
3548  }
3549 
3550 
3551 
3555  template <typename VectorType,
3556  typename std::enable_if<is_serial_or_dummy<VectorType>::value,
3557  VectorType>::type * = nullptr>
3558  void
3559  compress_start(const unsigned int /*component_in_block_vector*/,
3560  VectorType & /*vec*/)
3561  {}
3562 
3563 
3564 
3569  template <typename VectorType,
3570  typename std::enable_if<!has_compress_start<VectorType>::value &&
3571  !is_serial_or_dummy<VectorType>::value,
3572  VectorType>::type * = nullptr>
3573  void
3574  compress_start(const unsigned int component_in_block_vector,
3575  VectorType & vec)
3576  {
3577  (void)component_in_block_vector;
3578  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3579  vec.compress(::VectorOperation::add);
3580  }
3581 
3582 
3583 
3589  template <
3590  typename VectorType,
3591  typename std::enable_if<has_compress_start<VectorType>::value &&
3592  !has_exchange_on_subset<VectorType>::value,
3593  VectorType>::type * = nullptr>
3594  void
3595  compress_start(const unsigned int component_in_block_vector,
3596  VectorType & vec)
3597  {
3598  (void)component_in_block_vector;
3599  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3600  vec.compress_start(component_in_block_vector + channel_shift);
3601  }
3602 
3603 
3604 
3611  template <
3612  typename VectorType,
3613  typename std::enable_if<has_compress_start<VectorType>::value &&
3614  has_exchange_on_subset<VectorType>::value,
3615  VectorType>::type * = nullptr>
3616  void
3617  compress_start(const unsigned int component_in_block_vector,
3618  VectorType & vec)
3619  {
3620  static_assert(
3621  std::is_same<Number, typename VectorType::value_type>::value,
3622  "Type mismatch between VectorType and VectorDataExchange");
3623  (void)component_in_block_vector;
3624  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3625 
3626  if (vec.size() != 0)
3627  {
3628 # ifdef DEAL_II_WITH_MPI
3629  const unsigned int mf_component = find_vector_in_mf(vec);
3630 
3631  const auto &part = get_partitioner(mf_component);
3632 
3633  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3634  part.n_import_sm_procs() == 0)
3635  return;
3636 
3637  tmp_data[component_in_block_vector] =
3638  matrix_free.acquire_scratch_data_non_threadsafe();
3639  tmp_data[component_in_block_vector]->resize_fast(
3640  part.n_import_indices());
3641  AssertDimension(requests.size(), tmp_data.size());
3642 
3643  part.import_from_ghosted_array_start(
3645  component_in_block_vector * 2 + channel_shift,
3646  ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3647  vec.shared_vector_data(),
3648  ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3649  matrix_free.get_dof_info(mf_component)
3650  .vector_partitioner->n_ghost_indices()),
3651  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3652  part.n_import_indices()),
3653  this->requests[component_in_block_vector]);
3654 # endif
3655  }
3656  }
3657 
3658 
3659 
3664  template <typename VectorType,
3665  typename std::enable_if<!has_compress_start<VectorType>::value,
3666  VectorType>::type * = nullptr>
3667  void
3668  compress_finish(const unsigned int /*component_in_block_vector*/,
3669  VectorType & /*vec*/)
3670  {}
3671 
3672 
3673 
3679  template <
3680  typename VectorType,
3681  typename std::enable_if<has_compress_start<VectorType>::value &&
3682  !has_exchange_on_subset<VectorType>::value,
3683  VectorType>::type * = nullptr>
3684  void
3685  compress_finish(const unsigned int component_in_block_vector,
3686  VectorType & vec)
3687  {
3688  (void)component_in_block_vector;
3689  vec.compress_finish(::VectorOperation::add);
3690  }
3691 
3692 
3693 
3700  template <
3701  typename VectorType,
3702  typename std::enable_if<has_compress_start<VectorType>::value &&
3703  has_exchange_on_subset<VectorType>::value,
3704  VectorType>::type * = nullptr>
3705  void
3706  compress_finish(const unsigned int component_in_block_vector,
3707  VectorType & vec)
3708  {
3709  static_assert(
3710  std::is_same<Number, typename VectorType::value_type>::value,
3711  "Type mismatch between VectorType and VectorDataExchange");
3712  (void)component_in_block_vector;
3713  if (vec.size() != 0)
3714  {
3715 # ifdef DEAL_II_WITH_MPI
3716  AssertIndexRange(component_in_block_vector, tmp_data.size());
3717  AssertDimension(requests.size(), tmp_data.size());
3718 
3719  const unsigned int mf_component = find_vector_in_mf(vec);
3720 
3721  const auto &part = get_partitioner(mf_component);
3722 
3723  if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3724  part.n_import_sm_procs() != 0)
3725  {
3726  part.import_from_ghosted_array_finish(
3728  ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3729  vec.shared_vector_data(),
3730  ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3731  matrix_free.get_dof_info(mf_component)
3732  .vector_partitioner->n_ghost_indices()),
3734  tmp_data[component_in_block_vector]->begin(),
3735  part.n_import_indices()),
3736  this->requests[component_in_block_vector]);
3737 
3738  matrix_free.release_scratch_data_non_threadsafe(
3739  tmp_data[component_in_block_vector]);
3740  tmp_data[component_in_block_vector] = nullptr;
3741  }
3742 
3744  MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3745 # endif
3746  }
3747  }
3748 
3749 
3750 
3754  template <typename VectorType,
3755  typename std::enable_if<is_serial_or_dummy<VectorType>::value,
3756  VectorType>::type * = nullptr>
3757  void
3758  reset_ghost_values(const VectorType & /*vec*/) const
3759  {}
3760 
3761 
3762 
3767  template <
3768  typename VectorType,
3769  typename std::enable_if<!has_exchange_on_subset<VectorType>::value &&
3770  !is_serial_or_dummy<VectorType>::value,
3771  VectorType>::type * = nullptr>
3772  void
3773  reset_ghost_values(const VectorType &vec) const
3774  {
3775  if (ghosts_were_set == true)
3776  return;
3777 
3778  vec.zero_out_ghost_values();
3779  }
3780 
3781 
3782 
3788  template <typename VectorType,
3789  typename std::enable_if<has_exchange_on_subset<VectorType>::value,
3790  VectorType>::type * = nullptr>
3791  void
3792  reset_ghost_values(const VectorType &vec) const
3793  {
3794  static_assert(
3795  std::is_same<Number, typename VectorType::value_type>::value,
3796  "Type mismatch between VectorType and VectorDataExchange");
3797  if (ghosts_were_set == true)
3798  return;
3799 
3800  if (vec.size() != 0)
3801  {
3802 # ifdef DEAL_II_WITH_MPI
3803  AssertDimension(requests.size(), tmp_data.size());
3804 
3805  const unsigned int mf_component = find_vector_in_mf(vec);
3806 
3807  const auto &part = get_partitioner(mf_component);
3808 
3809  if (part.n_ghost_indices() > 0)
3810  {
3811  part.reset_ghost_values(ArrayView<Number>(
3813  .begin() +
3814  part.locally_owned_size(),
3815  matrix_free.get_dof_info(mf_component)
3816  .vector_partitioner->n_ghost_indices()));
3817  }
3818 
3819 # endif
3820  }
3821  // let vector know that it's not ghosted anymore
3822  vec.set_ghost_state(false);
3823  }
3824 
3825 
3826 
3832  template <typename VectorType,
3833  typename std::enable_if<has_exchange_on_subset<VectorType>::value,
3834  VectorType>::type * = nullptr>
3835  void
3836  zero_vector_region(const unsigned int range_index, VectorType &vec) const
3837  {
3838  static_assert(
3839  std::is_same<Number, typename VectorType::value_type>::value,
3840  "Type mismatch between VectorType and VectorDataExchange");
3841  if (range_index == numbers::invalid_unsigned_int)
3842  vec = Number();
3843  else
3844  {
3845  const unsigned int mf_component = find_vector_in_mf(vec, false);
3847  matrix_free.get_dof_info(mf_component);
3848  Assert(dof_info.vector_zero_range_list_index.empty() == false,
3849  ExcNotInitialized());
3850 
3851  Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3852  ExcInternalError());
3853  AssertIndexRange(range_index,
3854  dof_info.vector_zero_range_list_index.size() - 1);
3855  for (unsigned int id =
3856  dof_info.vector_zero_range_list_index[range_index];
3857  id != dof_info.vector_zero_range_list_index[range_index + 1];
3858  ++id)
3859  std::memset(vec.begin() + dof_info.vector_zero_range_list[id].first,
3860  0,
3861  (dof_info.vector_zero_range_list[id].second -
3862  dof_info.vector_zero_range_list[id].first) *
3863  sizeof(Number));
3864  }
3865  }
3866 
3867 
3868 
3874  template <
3875  typename VectorType,
3876  typename std::enable_if<!has_exchange_on_subset<VectorType>::value,
3877  VectorType>::type * = nullptr,
3878  typename VectorType::value_type * = nullptr>
3879  void
3880  zero_vector_region(const unsigned int range_index, VectorType &vec) const
3881  {
3882  if (range_index == numbers::invalid_unsigned_int || range_index == 0)
3883  vec = typename VectorType::value_type();
3884  }
3885 
3886 
3887 
3892  void
3893  zero_vector_region(const unsigned int, ...) const
3894  {
3895  Assert(false,
3896  ExcNotImplemented("Zeroing is only implemented for vector types "
3897  "which provide VectorType::value_type"));
3898  }
3899 
3900 
3901 
3902  const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
3903  const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3904  DataAccessOnFaces vector_face_access;
3905  bool ghosts_were_set;
3906 # ifdef DEAL_II_WITH_MPI
3907  std::vector<AlignedVector<Number> *> tmp_data;
3908  std::vector<std::vector<MPI_Request>> requests;
3909 # endif
3910  }; // VectorDataExchange
3911 
3912  template <typename VectorStruct>
3913  unsigned int
3914  n_components(const VectorStruct &vec);
3915 
3916  template <typename VectorStruct>
3917  unsigned int
3918  n_components_block(const VectorStruct &vec,
3919  std::integral_constant<bool, true>)
3920  {
3921  unsigned int components = 0;
3922  for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
3923  components += n_components(vec.block(bl));
3924  return components;
3925  }
3926 
3927  template <typename VectorStruct>
3928  unsigned int
3929  n_components_block(const VectorStruct &, std::integral_constant<bool, false>)
3930  {
3931  return 1;
3932  }
3933 
3934  template <typename VectorStruct>
3935  unsigned int
3936  n_components(const VectorStruct &vec)
3937  {
3938  return n_components_block(
3939  vec, std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3940  }
3941 
3942  template <typename VectorStruct>
3943  inline unsigned int
3944  n_components(const std::vector<VectorStruct> &vec)
3945  {
3946  unsigned int components = 0;
3947  for (unsigned int comp = 0; comp < vec.size(); comp++)
3948  components += n_components_block(
3949  vec[comp],
3950  std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3951  return components;
3952  }
3953 
3954  template <typename VectorStruct>
3955  inline unsigned int
3956  n_components(const std::vector<VectorStruct *> &vec)
3957  {
3958  unsigned int components = 0;
3959  for (unsigned int comp = 0; comp < vec.size(); comp++)
3960  components += n_components_block(
3961  *vec[comp],
3962  std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3963  return components;
3964  }
3965 
3966 
3967 
3968  // A helper function to identify block vectors with many components where we
3969  // should not try to overlap computations and communication because there
3970  // would be too many outstanding communication requests.
3971 
3972  // default value for vectors that do not have communication_block_size
3973  template <
3974  typename VectorStruct,
3975  typename std::enable_if<!has_communication_block_size<VectorStruct>::value,
3976  VectorStruct>::type * = nullptr>
3977  constexpr unsigned int
3978  get_communication_block_size(const VectorStruct &)
3979  {
3981  }
3982 
3983 
3984 
3985  template <
3986  typename VectorStruct,
3987  typename std::enable_if<has_communication_block_size<VectorStruct>::value,
3988  VectorStruct>::type * = nullptr>
3989  constexpr unsigned int
3990  get_communication_block_size(const VectorStruct &)
3991  {
3992  return VectorStruct::communication_block_size;
3993  }
3994 
3995 
3996 
3997  // --------------------------------------------------------------------------
3998  // below we have wrappers to distinguish between block and non-block vectors.
3999  // --------------------------------------------------------------------------
4000 
4001  //
4002  // update_ghost_values_start
4003  //
4004 
4005  // update_ghost_values for block vectors
4006  template <int dim,
4007  typename VectorStruct,
4008  typename Number,
4009  typename VectorizedArrayType,
4010  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4011  VectorStruct>::type * = nullptr>
4012  void
4013  update_ghost_values_start(
4014  const VectorStruct & vec,
4015  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4016  const unsigned int channel = 0)
4017  {
4018  if (get_communication_block_size(vec) < vec.n_blocks())
4019  {
4020  // don't forget to set ghosts_were_set, that otherwise happens
4021  // inside VectorDataExchange::update_ghost_values_start()
4022  exchanger.ghosts_were_set = vec.has_ghost_elements();
4023  vec.update_ghost_values();
4024  }
4025  else
4026  {
4027  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4028  update_ghost_values_start(vec.block(i), exchanger, channel + i);
4029  }
4030  }
4031 
4032 
4033 
4034  // update_ghost_values for non-block vectors
4035  template <int dim,
4036  typename VectorStruct,
4037  typename Number,
4038  typename VectorizedArrayType,
4039  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4040  VectorStruct>::type * = nullptr>
4041  void
4042  update_ghost_values_start(
4043  const VectorStruct & vec,
4044  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4045  const unsigned int channel = 0)
4046  {
4047  exchanger.update_ghost_values_start(channel, vec);
4048  }
4049 
4050 
4051 
4052  // update_ghost_values_start() for vector of vectors
4053  template <int dim,
4054  typename VectorStruct,
4055  typename Number,
4056  typename VectorizedArrayType>
4057  inline void
4058  update_ghost_values_start(
4059  const std::vector<VectorStruct> & vec,
4060  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4061  {
4062  unsigned int component_index = 0;
4063  for (unsigned int comp = 0; comp < vec.size(); comp++)
4064  {
4065  update_ghost_values_start(vec[comp], exchanger, component_index);
4066  component_index += n_components(vec[comp]);
4067  }
4068  }
4069 
4070 
4071 
4072  // update_ghost_values_start() for vector of pointers to vectors
4073  template <int dim,
4074  typename VectorStruct,
4075  typename Number,
4076  typename VectorizedArrayType>
4077  inline void
4078  update_ghost_values_start(
4079  const std::vector<VectorStruct *> & vec,
4080  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4081  {
4082  unsigned int component_index = 0;
4083  for (unsigned int comp = 0; comp < vec.size(); comp++)
4084  {
4085  update_ghost_values_start(*vec[comp], exchanger, component_index);
4086  component_index += n_components(*vec[comp]);
4087  }
4088  }
4089 
4090 
4091 
4092  //
4093  // update_ghost_values_finish
4094  //
4095 
4096  // for block vectors
4097  template <int dim,
4098  typename VectorStruct,
4099  typename Number,
4100  typename VectorizedArrayType,
4101  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4102  VectorStruct>::type * = nullptr>
4103  void
4104  update_ghost_values_finish(
4105  const VectorStruct & vec,
4106  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4107  const unsigned int channel = 0)
4108  {
4109  if (get_communication_block_size(vec) < vec.n_blocks())
4110  {
4111  // do nothing, everything has already been completed in the _start()
4112  // call
4113  }
4114  else
4115  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4116  update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4117  }
4118 
4119 
4120 
4121  // for non-block vectors
4122  template <int dim,
4123  typename VectorStruct,
4124  typename Number,
4125  typename VectorizedArrayType,
4126  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4127  VectorStruct>::type * = nullptr>
4128  void
4129  update_ghost_values_finish(
4130  const VectorStruct & vec,
4131  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4132  const unsigned int channel = 0)
4133  {
4134  exchanger.update_ghost_values_finish(channel, vec);
4135  }
4136 
4137 
4138 
4139  // for vector of vectors
4140  template <int dim,
4141  typename VectorStruct,
4142  typename Number,
4143  typename VectorizedArrayType>
4144  inline void
4145  update_ghost_values_finish(
4146  const std::vector<VectorStruct> & vec,
4147  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4148  {
4149  unsigned int component_index = 0;
4150  for (unsigned int comp = 0; comp < vec.size(); comp++)
4151  {
4152  update_ghost_values_finish(vec[comp], exchanger, component_index);
4153  component_index += n_components(vec[comp]);
4154  }
4155  }
4156 
4157 
4158 
4159  // for vector of pointers to vectors
4160  template <int dim,
4161  typename VectorStruct,
4162  typename Number,
4163  typename VectorizedArrayType>
4164  inline void
4165  update_ghost_values_finish(
4166  const std::vector<VectorStruct *> & vec,
4167  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4168  {
4169  unsigned int component_index = 0;
4170  for (unsigned int comp = 0; comp < vec.size(); comp++)
4171  {
4172  update_ghost_values_finish(*vec[comp], exchanger, component_index);
4173  component_index += n_components(*vec[comp]);
4174  }
4175  }
4176 
4177 
4178 
4179  //
4180  // compress_start
4181  //
4182 
4183  // for block vectors
4184  template <int dim,
4185  typename VectorStruct,
4186  typename Number,
4187  typename VectorizedArrayType,
4188  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4189  VectorStruct>::type * = nullptr>
4190  inline void
4191  compress_start(
4192  VectorStruct & vec,
4193  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4194  const unsigned int channel = 0)
4195  {
4196  if (get_communication_block_size(vec) < vec.n_blocks())
4197  vec.compress(::VectorOperation::add);
4198  else
4199  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4200  compress_start(vec.block(i), exchanger, channel + i);
4201  }
4202 
4203 
4204 
4205  // for non-block vectors
4206  template <int dim,
4207  typename VectorStruct,
4208  typename Number,
4209  typename VectorizedArrayType,
4210  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4211  VectorStruct>::type * = nullptr>
4212  inline void
4213  compress_start(
4214  VectorStruct & vec,
4215  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4216  const unsigned int channel = 0)
4217  {
4218  exchanger.compress_start(channel, vec);
4219  }
4220 
4221 
4222 
4223  // for std::vector of vectors
4224  template <int dim,
4225  typename VectorStruct,
4226  typename Number,
4227  typename VectorizedArrayType>
4228  inline void
4229  compress_start(
4230  std::vector<VectorStruct> & vec,
4231  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4232  {
4233  unsigned int component_index = 0;
4234  for (unsigned int comp = 0; comp < vec.size(); comp++)
4235  {
4236  compress_start(vec[comp], exchanger, component_index);
4237  component_index += n_components(vec[comp]);
4238  }
4239  }
4240 
4241 
4242 
4243  // for std::vector of pointer to vectors
4244  template <int dim,
4245  typename VectorStruct,
4246  typename Number,
4247  typename VectorizedArrayType>
4248  inline void
4249  compress_start(
4250  std::vector<VectorStruct *> & vec,
4251  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4252  {
4253  unsigned int component_index = 0;
4254  for (unsigned int comp = 0; comp < vec.size(); comp++)
4255  {
4256  compress_start(*vec[comp], exchanger, component_index);
4257  component_index += n_components(*vec[comp]);
4258  }
4259  }
4260 
4261 
4262 
4263  //
4264  // compress_finish
4265  //
4266 
4267  // for block vectors
4268  template <int dim,
4269  typename VectorStruct,
4270  typename Number,
4271  typename VectorizedArrayType,
4272  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4273  VectorStruct>::type * = nullptr>
4274  inline void
4275  compress_finish(
4276  VectorStruct & vec,
4277  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4278  const unsigned int channel = 0)
4279  {
4280  if (get_communication_block_size(vec) < vec.n_blocks())
4281  {
4282  // do nothing, everything has already been completed in the _start()
4283  // call
4284  }
4285  else
4286  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4287  compress_finish(vec.block(i), exchanger, channel + i);
4288  }
4289 
4290 
4291 
4292  // for non-block vectors
4293  template <int dim,
4294  typename VectorStruct,
4295  typename Number,
4296  typename VectorizedArrayType,
4297  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4298  VectorStruct>::type * = nullptr>
4299  inline void
4300  compress_finish(
4301  VectorStruct & vec,
4302  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4303  const unsigned int channel = 0)
4304  {
4305  exchanger.compress_finish(channel, vec);
4306  }
4307 
4308 
4309 
4310  // for std::vector of vectors
4311  template <int dim,
4312  typename VectorStruct,
4313  typename Number,
4314  typename VectorizedArrayType>
4315  inline void
4316  compress_finish(
4317  std::vector<VectorStruct> & vec,
4318  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4319  {
4320  unsigned int component_index = 0;
4321  for (unsigned int comp = 0; comp < vec.size(); comp++)
4322  {
4323  compress_finish(vec[comp], exchanger, component_index);
4324  component_index += n_components(vec[comp]);
4325  }
4326  }
4327 
4328 
4329 
4330  // for std::vector of pointer to vectors
4331  template <int dim,
4332  typename VectorStruct,
4333  typename Number,
4334  typename VectorizedArrayType>
4335  inline void
4336  compress_finish(
4337  std::vector<VectorStruct *> & vec,
4338  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4339  {
4340  unsigned int component_index = 0;
4341  for (unsigned int comp = 0; comp < vec.size(); comp++)
4342  {
4343  compress_finish(*vec[comp], exchanger, component_index);
4344  component_index += n_components(*vec[comp]);
4345  }
4346  }
4347 
4348 
4349 
4350  //
4351  // reset_ghost_values:
4352  //
4353  // if the input vector did not have ghosts imported, clear them here again
4354  // in order to avoid subsequent operations e.g. in linear solvers to work
4355  // with ghosts all the time
4356  //
4357 
4358  // for block vectors
4359  template <int dim,
4360  typename VectorStruct,
4361  typename Number,
4362  typename VectorizedArrayType,
4363  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4364  VectorStruct>::type * = nullptr>
4365  inline void
4366  reset_ghost_values(
4367  const VectorStruct & vec,
4368  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4369  {
4370  // return immediately if there is nothing to do.
4371  if (exchanger.ghosts_were_set == true)
4372  return;
4373 
4374  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4375  reset_ghost_values(vec.block(i), exchanger);
4376  }
4377 
4378 
4379 
4380  // for non-block vectors
4381  template <int dim,
4382  typename VectorStruct,
4383  typename Number,
4384  typename VectorizedArrayType,
4385  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4386  VectorStruct>::type * = nullptr>
4387  inline void
4388  reset_ghost_values(
4389  const VectorStruct & vec,
4390  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4391  {
4392  exchanger.reset_ghost_values(vec);
4393  }
4394 
4395 
4396 
4397  // for std::vector of vectors
4398  template <int dim,
4399  typename VectorStruct,
4400  typename Number,
4401  typename VectorizedArrayType>
4402  inline void
4403  reset_ghost_values(
4404  const std::vector<VectorStruct> & vec,
4405  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4406  {
4407  // return immediately if there is nothing to do.
4408  if (exchanger.ghosts_were_set == true)
4409  return;
4410 
4411  for (unsigned int comp = 0; comp < vec.size(); comp++)
4412  reset_ghost_values(vec[comp], exchanger);
4413  }
4414 
4415 
4416 
4417  // for std::vector of pointer to vectors
4418  template <int dim,
4419  typename VectorStruct,
4420  typename Number,
4421  typename VectorizedArrayType>
4422  inline void
4423  reset_ghost_values(
4424  const std::vector<VectorStruct *> & vec,
4425  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4426  {
4427  // return immediately if there is nothing to do.
4428  if (exchanger.ghosts_were_set == true)
4429  return;
4430 
4431  for (unsigned int comp = 0; comp < vec.size(); comp++)
4432  reset_ghost_values(*vec[comp], exchanger);
4433  }
4434 
4435 
4436 
4437  //
4438  // zero_vector_region
4439  //
4440 
4441  // for block vectors
4442  template <int dim,
4443  typename VectorStruct,
4444  typename Number,
4445  typename VectorizedArrayType,
4446  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4447  VectorStruct>::type * = nullptr>
4448  inline void
4449  zero_vector_region(
4450  const unsigned int range_index,
4451  VectorStruct & vec,
4452  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4453  {
4454  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4455  exchanger.zero_vector_region(range_index, vec.block(i));
4456  }
4457 
4458 
4459 
4460  // for non-block vectors
4461  template <int dim,
4462  typename VectorStruct,
4463  typename Number,
4464  typename VectorizedArrayType,
4465  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4466  VectorStruct>::type * = nullptr>
4467  inline void
4468  zero_vector_region(
4469  const unsigned int range_index,
4470  VectorStruct & vec,
4471  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4472  {
4473  exchanger.zero_vector_region(range_index, vec);
4474  }
4475 
4476 
4477 
4478  // for std::vector of vectors
4479  template <int dim,
4480  typename VectorStruct,
4481  typename Number,
4482  typename VectorizedArrayType>
4483  inline void
4484  zero_vector_region(
4485  const unsigned int range_index,
4486  std::vector<VectorStruct> & vec,
4487  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4488  {
4489  for (unsigned int comp = 0; comp < vec.size(); comp++)
4490  zero_vector_region(range_index, vec[comp], exchanger);
4491  }
4492 
4493 
4494 
4495  // for std::vector of pointers to vectors
4496  template <int dim,
4497  typename VectorStruct,
4498  typename Number,
4499  typename VectorizedArrayType>
4500  inline void
4501  zero_vector_region(
4502  const unsigned int range_index,
4503  std::vector<VectorStruct *> & vec,
4504  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4505  {
4506  for (unsigned int comp = 0; comp < vec.size(); comp++)
4507  zero_vector_region(range_index, *vec[comp], exchanger);
4508  }
4509 
4510 
4511 
4512  namespace MatrixFreeFunctions
4513  {
4514  // struct to select between a const interface and a non-const interface
4515  // for MFWorker
4516  template <typename, typename, typename, typename, bool>
4517  struct InterfaceSelector
4518  {};
4519 
4520  // Version for constant functions
4521  template <typename MF,
4522  typename InVector,
4523  typename OutVector,
4524  typename Container>
4525  struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4526  {
4527  using function_type = void (Container::*)(
4528  const MF &,
4529  OutVector &,
4530  const InVector &,
4531  const std::pair<unsigned int, unsigned int> &) const;
4532  };
4533 
4534  // Version for non-constant functions
4535  template <typename MF,
4536  typename InVector,
4537  typename OutVector,
4538  typename Container>
4539  struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4540  {
4541  using function_type =
4542  void (Container::*)(const MF &,
4543  OutVector &,
4544  const InVector &,
4545  const std::pair<unsigned int, unsigned int> &);
4546  };
4547  } // namespace MatrixFreeFunctions
4548 
4549 
4550 
4551  // A implementation class for the worker object that runs the various
4552  // operations we want to perform during the matrix-free loop
4553  template <typename MF,
4554  typename InVector,
4555  typename OutVector,
4556  typename Container,
4557  bool is_constant>
4558  class MFWorker : public MFWorkerInterface
4559  {
4560  public:
4561  // An alias to make the arguments further down more readable
4562  using function_type = typename MatrixFreeFunctions::
4563  InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4564  function_type;
4565 
4566  // constructor, binds all the arguments to this class
4567  MFWorker(const MF & matrix_free,
4568  const InVector & src,
4569  OutVector & dst,
4570  const bool zero_dst_vector_setting,
4571  const Container & container,
4572  function_type cell_function,
4573  function_type face_function,
4574  function_type boundary_function,
4575  const typename MF::DataAccessOnFaces src_vector_face_access =
4577  const typename MF::DataAccessOnFaces dst_vector_face_access =
4579  const std::function<void(const unsigned int, const unsigned int)>
4580  &operation_before_loop = {},
4581  const std::function<void(const unsigned int, const unsigned int)>
4582  & operation_after_loop = {},
4583  const unsigned int dof_handler_index_pre_post = 0)
4584  : matrix_free(matrix_free)
4585  , container(const_cast<Container &>(container))
4586  , cell_function(cell_function)
4587  , face_function(face_function)
4588  , boundary_function(boundary_function)
4589  , src(src)
4590  , dst(dst)
4591  , src_data_exchanger(matrix_free,
4592  src_vector_face_access,
4593  n_components(src))
4594  , dst_data_exchanger(matrix_free,
4595  dst_vector_face_access,
4596  n_components(dst))
4597  , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4598  , zero_dst_vector_setting(zero_dst_vector_setting &&
4599  !src_and_dst_are_same)
4600  , operation_before_loop(operation_before_loop)
4601  , operation_after_loop(operation_after_loop)
4602  , dof_handler_index_pre_post(dof_handler_index_pre_post)
4603  {}
4604 
4605  // Runs the cell work. If no function is given, nothing is done
4606  virtual void
4607  cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4608  {
4609  if (cell_function != nullptr && cell_range.second > cell_range.first)
4610  for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4611  {
4612  const auto cell_subrange =
4613  matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4614 
4615  if (cell_subrange.second <= cell_subrange.first)
4616  continue;
4617 
4618  (container.*
4619  cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4620  }
4621  }
4622 
4623  virtual void
4624  cell(const unsigned int range_index) override
4625  {
4626  process_range(cell_function,
4627  matrix_free.get_task_info().cell_partition_data_hp_ptr,
4628  matrix_free.get_task_info().cell_partition_data_hp,
4629  range_index);
4630  }
4631 
4632  virtual void
4633  face(const unsigned int range_index) override
4634  {
4635  process_range(face_function,
4636  matrix_free.get_task_info().face_partition_data_hp_ptr,
4637  matrix_free.get_task_info().face_partition_data_hp,
4638  range_index);
4639  }
4640 
4641  virtual void
4642  boundary(const unsigned int range_index) override
4643  {
4644  process_range(boundary_function,
4645  matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4646  matrix_free.get_task_info().boundary_partition_data_hp,
4647  range_index);
4648  }
4649 
4650  private:
4651  void
4652  process_range(const function_type & fu,
4653  const std::vector<unsigned int> &ptr,
4654  const std::vector<unsigned int> &data,
4655  const unsigned int range_index)
4656  {
4657  if (fu == nullptr)
4658  return;
4659 
4660  for (unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4661  (container.*fu)(matrix_free,
4662  this->dst,
4663  this->src,
4664  std::make_pair(data[2 * i], data[2 * i + 1]));
4665  }
4666 
4667  public:
4668  // Starts the communication for the update ghost values operation. We
4669  // cannot call this update if ghost and destination are the same because
4670  // that would introduce spurious entries in the destination (there is also
4671  // the problem that reading from a vector that we also write to is usually
4672  // not intended in case there is overlap, but this is up to the
4673  // application code to decide and we cannot catch this case here).
4674  virtual void
4675  vector_update_ghosts_start() override
4676  {
4677  if (!src_and_dst_are_same)
4678  internal::update_ghost_values_start(src, src_data_exchanger);
4679  }
4680 
4681  // Finishes the communication for the update ghost values operation
4682  virtual void
4683  vector_update_ghosts_finish() override
4684  {
4685  if (!src_and_dst_are_same)
4686  internal::update_ghost_values_finish(src, src_data_exchanger);
4687  }
4688 
4689  // Starts the communication for the vector compress operation
4690  virtual void
4691  vector_compress_start() override
4692  {
4693  internal::compress_start(dst, dst_data_exchanger);
4694  }
4695 
4696  // Finishes the communication for the vector compress operation
4697  virtual void
4698  vector_compress_finish() override
4699  {
4700  internal::compress_finish(dst, dst_data_exchanger);
4701  if (!src_and_dst_are_same)
4702  internal::reset_ghost_values(src, src_data_exchanger);
4703  }
4704 
4705  // Zeros the given input vector
4706  virtual void
4707  zero_dst_vector_range(const unsigned int range_index) override
4708  {
4709  if (zero_dst_vector_setting)
4710  internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4711  }
4712 
4713  virtual void
4714  cell_loop_pre_range(const unsigned int range_index) override
4715  {
4716  if (operation_before_loop)
4717  {
4719  matrix_free.get_dof_info(dof_handler_index_pre_post);
4720  if (range_index == numbers::invalid_unsigned_int)
4721  {
4722  // Case with threaded loop -> currently no overlap implemented
4724  0U,
4725  dof_info.vector_partitioner->locally_owned_size(),
4726  operation_before_loop,
4728  }
4729  else
4730  {
4731  AssertIndexRange(range_index,
4732  dof_info.cell_loop_pre_list_index.size() - 1);
4733  for (unsigned int id =
4734  dof_info.cell_loop_pre_list_index[range_index];
4735  id != dof_info.cell_loop_pre_list_index[range_index + 1];
4736  ++id)
4737  operation_before_loop(dof_info.cell_loop_pre_list[id].first,
4738  dof_info.cell_loop_pre_list[id].second);
4739  }
4740  }
4741  }
4742 
4743  virtual void
4744  cell_loop_post_range(const unsigned int range_index) override
4745  {
4746  if (operation_after_loop)
4747  {
4749  matrix_free.get_dof_info(dof_handler_index_pre_post);
4750  if (range_index == numbers::invalid_unsigned_int)
4751  {
4752  // Case with threaded loop -> currently no overlap implemented
4754  0U,
4755  dof_info.vector_partitioner->locally_owned_size(),
4756  operation_after_loop,
4758  }
4759  else
4760  {
4761  AssertIndexRange(range_index,
4762  dof_info.cell_loop_post_list_index.size() - 1);
4763  for (unsigned int id =
4764  dof_info.cell_loop_post_list_index[range_index];
4765  id != dof_info.cell_loop_post_list_index[range_index + 1];
4766  ++id)
4767  operation_after_loop(dof_info.cell_loop_post_list[id].first,
4768  dof_info.cell_loop_post_list[id].second);
4769  }
4770  }
4771  }
4772 
4773  private:
4774  const MF & matrix_free;
4775  Container & container;
4776  function_type cell_function;
4777  function_type face_function;
4778  function_type boundary_function;
4779 
4780  const InVector &src;
4781  OutVector & dst;
4782  VectorDataExchange<MF::dimension,
4783  typename MF::value_type,
4784  typename MF::vectorized_value_type>
4785  src_data_exchanger;
4786  VectorDataExchange<MF::dimension,
4787  typename MF::value_type,
4788  typename MF::vectorized_value_type>
4789  dst_data_exchanger;
4790  const bool src_and_dst_are_same;
4791  const bool zero_dst_vector_setting;
4792  const std::function<void(const unsigned int, const unsigned int)>
4793  operation_before_loop;
4794  const std::function<void(const unsigned int, const unsigned int)>
4795  operation_after_loop;
4796  const unsigned int dof_handler_index_pre_post;
4797  };
4798 
4799 
4800 
4805  template <class MF, typename InVector, typename OutVector>
4806  struct MFClassWrapper
4807  {
4808  using function_type =
4809  std::function<void(const MF &,
4810  OutVector &,
4811  const InVector &,
4812  const std::pair<unsigned int, unsigned int> &)>;
4813 
4814  MFClassWrapper(const function_type cell,
4815  const function_type face,
4816  const function_type boundary)
4817  : cell(cell)
4818  , face(face)
4819  , boundary(boundary)
4820  {}
4821 
4822  void
4823  cell_integrator(const MF & mf,
4824  OutVector & dst,
4825  const InVector & src,
4826  const std::pair<unsigned int, unsigned int> &range) const
4827  {
4828  if (cell)
4829  cell(mf, dst, src, range);
4830  }
4831 
4832  void
4833  face_integrator(const MF & mf,
4834  OutVector & dst,
4835  const InVector & src,
4836  const std::pair<unsigned int, unsigned int> &range) const
4837  {
4838  if (face)
4839  face(mf, dst, src, range);
4840  }
4841 
4842  void
4843  boundary_integrator(
4844  const MF & mf,
4845  OutVector & dst,
4846  const InVector & src,
4847  const std::pair<unsigned int, unsigned int> &range) const
4848  {
4849  if (boundary)
4850  boundary(mf, dst, src, range);
4851  }
4852 
4853  const function_type cell;
4854  const function_type face;
4855  const function_type boundary;
4856  };
4857 
4858 } // end of namespace internal
4859 
4860 
4861 
4862 template <int dim, typename Number, typename VectorizedArrayType>
4863 template <typename OutVector, typename InVector>
4864 inline void
4866  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4867  OutVector &,
4868  const InVector &,
4869  const std::pair<unsigned int, unsigned int> &)>
4870  & cell_operation,
4871  OutVector & dst,
4872  const InVector &src,
4873  const bool zero_dst_vector) const
4874 {
4875  using Wrapper =
4876  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4877  InVector,
4878  OutVector>;
4879  Wrapper wrap(cell_operation, nullptr, nullptr);
4880  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4881  InVector,
4882  OutVector,
4883  Wrapper,
4884  true>
4885  worker(*this,
4886  src,
4887  dst,
4888  zero_dst_vector,
4889  wrap,
4890  &Wrapper::cell_integrator,
4891  &Wrapper::face_integrator,
4892  &Wrapper::boundary_integrator);
4893 
4894  task_info.loop(worker);
4895 }
4896 
4897 
4898 
4899 template <int dim, typename Number, typename VectorizedArrayType>
4900 template <typename OutVector, typename InVector>
4901 inline void
4903  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4904  OutVector &,
4905  const InVector &,
4906  const std::pair<unsigned int, unsigned int> &)>
4907  & cell_operation,
4908  OutVector & dst,
4909  const InVector &src,
4910  const std::function<void(const unsigned int, const unsigned int)>
4911  &operation_before_loop,
4912  const std::function<void(const unsigned int, const unsigned int)>
4913  & operation_after_loop,
4914  const unsigned int dof_handler_index_pre_post) const
4915 {
4916  using Wrapper =
4917  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4918  InVector,
4919  OutVector>;
4920  Wrapper wrap(cell_operation, nullptr, nullptr);
4921  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4922  InVector,
4923  OutVector,
4924  Wrapper,
4925  true>
4926  worker(*this,
4927  src,
4928  dst,
4929  false,
4930  wrap,
4931  &Wrapper::cell_integrator,
4932  &Wrapper::face_integrator,
4933  &Wrapper::boundary_integrator,
4936  operation_before_loop,
4937  operation_after_loop,
4938  dof_handler_index_pre_post);
4939 
4940  task_info.loop(worker);
4941 }
4942 
4943 
4944 
4945 template <int dim, typename Number, typename VectorizedArrayType>
4946 template <typename OutVector, typename InVector>
4947 inline void
4949  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4950  OutVector &,
4951  const InVector &,
4952  const std::pair<unsigned int, unsigned int> &)>
4953  &cell_operation,
4954  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4955  OutVector &,
4956  const InVector &,
4957  const std::pair<unsigned int, unsigned int> &)>
4958  &face_operation,
4959  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4960  OutVector &,
4961  const InVector &,
4962  const std::pair<unsigned int, unsigned int> &)>
4963  & boundary_operation,
4964  OutVector & dst,
4965  const InVector & src,
4966  const bool zero_dst_vector,
4967  const DataAccessOnFaces dst_vector_face_access,
4968  const DataAccessOnFaces src_vector_face_access) const
4969 {
4970  using Wrapper =
4971  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4972  InVector,
4973  OutVector>;
4974  Wrapper wrap(cell_operation, face_operation, boundary_operation);
4975  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4976  InVector,
4977  OutVector,
4978  Wrapper,
4979  true>
4980  worker(*this,
4981  src,
4982  dst,
4983  zero_dst_vector,
4984  wrap,
4985  &Wrapper::cell_integrator,
4986  &Wrapper::face_integrator,
4987  &Wrapper::boundary_integrator,
4988  src_vector_face_access,
4989  dst_vector_face_access);
4990 
4991  task_info.loop(worker);
4992 }
4993 
4994 
4995 
4996 template <int dim, typename Number, typename VectorizedArrayType>
4997 template <typename CLASS, typename OutVector, typename InVector>
4998 inline void
5000  void (CLASS::*function_pointer)(
5002  OutVector &,
5003  const InVector &,
5004  const std::pair<unsigned int, unsigned int> &) const,
5005  const CLASS * owning_class,
5006  OutVector & dst,
5007  const InVector &src,
5008  const bool zero_dst_vector) const
5009 {
5010  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5011  InVector,
5012  OutVector,
5013  CLASS,
5014  true>
5015  worker(*this,
5016  src,
5017  dst,
5018  zero_dst_vector,
5019  *owning_class,
5020  function_pointer,
5021  nullptr,
5022  nullptr);
5023  task_info.loop(worker);
5024 }
5025 
5026 
5027 
5028 template <int dim, typename Number, typename VectorizedArrayType>
5029 template <typename CLASS, typename OutVector, typename InVector>
5030 inline void
5032  void (CLASS::*function_pointer)(
5034  OutVector &,
5035  const InVector &,
5036  const std::pair<unsigned int, unsigned int> &) const,
5037  const CLASS * owning_class,
5038  OutVector & dst,
5039  const InVector &src,
5040  const std::function<void(const unsigned int, const unsigned int)>
5041  &operation_before_loop,
5042  const std::function<void(const unsigned int, const unsigned int)>
5043  & operation_after_loop,
5044  const unsigned int dof_handler_index_pre_post) const
5045 {
5046  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5047  InVector,
5048  OutVector,
5049  CLASS,
5050  true>
5051  worker(*this,
5052  src,
5053  dst,
5054  false,
5055  *owning_class,
5056  function_pointer,
5057  nullptr,
5058  nullptr,
5061  operation_before_loop,
5062  operation_after_loop,
5063  dof_handler_index_pre_post);
5064  task_info.loop(worker);
5065 }
5066 
5067 
5068 
5069 template <int dim, typename Number, typename VectorizedArrayType>
5070 template <typename CLASS, typename OutVector, typename InVector>
5071 inline void
5073  void (CLASS::*cell_operation)(
5075  OutVector &,
5076  const InVector &,
5077  const std::pair<unsigned int, unsigned int> &) const,
5078  void (CLASS::*face_operation)(
5080  OutVector &,
5081  const InVector &,
5082  const std::pair<unsigned int, unsigned int> &) const,
5083  void (CLASS::*boundary_operation)(
5085  OutVector &,
5086  const InVector &,
5087  const std::pair<unsigned int, unsigned int> &) const,
5088  const CLASS * owning_class,
5089  OutVector & dst,
5090  const InVector & src,
5091  const bool zero_dst_vector,
5092  const DataAccessOnFaces dst_vector_face_access,
5093  const DataAccessOnFaces src_vector_face_access) const
5094 {
5095  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5096  InVector,
5097  OutVector,
5098  CLASS,
5099  true>
5100  worker(*this,
5101  src,
5102  dst,
5103  zero_dst_vector,
5104  *owning_class,
5105  cell_operation,
5106  face_operation,
5107  boundary_operation,
5108  src_vector_face_access,
5109  dst_vector_face_access);
5110  task_info.loop(worker);
5111 }
5112 
5113 
5114 
5115 template <int dim, typename Number, typename VectorizedArrayType>
5116 template <typename CLASS, typename OutVector, typename InVector>
5117 inline void
5119  void (CLASS::*function_pointer)(
5121  OutVector &,
5122  const InVector &,
5123  const std::pair<unsigned int, unsigned int> &),
5124  CLASS * owning_class,
5125  OutVector & dst,
5126  const InVector &src,
5127  const bool zero_dst_vector) const
5128 {
5129  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5130  InVector,
5131  OutVector,
5132  CLASS,
5133  false>
5134  worker(*this,
5135  src,
5136  dst,
5137  zero_dst_vector,
5138  *owning_class,
5139  function_pointer,
5140  nullptr,
5141  nullptr);
5142  task_info.loop(worker);
5143 }
5144 
5145 
5146 
5147 template <int dim, typename Number, typename VectorizedArrayType>
5148 template <typename CLASS, typename OutVector, typename InVector>
5149 inline void
5151  void (CLASS::*function_pointer)(
5153  OutVector &,
5154  const InVector &,
5155  const std::pair<unsigned int, unsigned int> &),
5156  CLASS * owning_class,
5157  OutVector & dst,
5158  const InVector &src,
5159  const std::function<void(const unsigned int, const unsigned int)>
5160  &operation_before_loop,
5161  const std::function<void(const unsigned int, const unsigned int)>
5162  & operation_after_loop,
5163  const unsigned int dof_handler_index_pre_post) const
5164 {
5165  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5166  InVector,
5167  OutVector,
5168  CLASS,
5169  false>
5170  worker(*this,
5171  src,
5172  dst,
5173  false,
5174  *owning_class,
5175  function_pointer,
5176  nullptr,
5177  nullptr,
5180  operation_before_loop,
5181  operation_after_loop,
5182  dof_handler_index_pre_post);
5183  task_info.loop(worker);
5184 }
5185 
5186 
5187 
5188 template <int dim, typename Number, typename VectorizedArrayType>
5189 template <typename CLASS, typename OutVector, typename InVector>
5190 inline void
5192  void (CLASS::*cell_operation)(
5194  OutVector &,
5195  const InVector &,
5196  const std::pair<unsigned int, unsigned int> &),
5197  void (CLASS::*face_operation)(
5199  OutVector &,
5200  const InVector &,
5201  const std::pair<unsigned int, unsigned int> &),
5202  void (CLASS::*boundary_operation)(
5204  OutVector &,
5205  const InVector &,
5206  const std::pair<unsigned int, unsigned int> &),
5207  CLASS * owning_class,
5208  OutVector & dst,
5209  const InVector & src,
5210  const bool zero_dst_vector,
5211  const DataAccessOnFaces dst_vector_face_access,
5212  const DataAccessOnFaces src_vector_face_access) const
5213 {
5214  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5215  InVector,
5216  OutVector,
5217  CLASS,
5218  false>
5219  worker(*this,
5220  src,
5221  dst,
5222  zero_dst_vector,
5223  *owning_class,
5224  cell_operation,
5225  face_operation,
5226  boundary_operation,
5227  src_vector_face_access,
5228  dst_vector_face_access);
5229  task_info.loop(worker);
5230 }
5231 
5232 
5233 
5234 template <int dim, typename Number, typename VectorizedArrayType>
5235 template <typename CLASS, typename OutVector, typename InVector>
5236 inline void
5238  void (CLASS::*function_pointer)(
5240  OutVector &,
5241  const InVector &,
5242  const std::pair<unsigned int, unsigned int> &) const,
5243  const CLASS * owning_class,
5244  OutVector & dst,
5245  const InVector & src,
5246  const bool zero_dst_vector,
5247  const DataAccessOnFaces src_vector_face_access) const
5248 {
5249  auto src_vector_face_access_temp = src_vector_face_access;
5250  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5251  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5252  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5253  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5254 
5255  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5256  InVector,
5257  OutVector,
5258  CLASS,
5259  true>
5260  worker(*this,
5261  src,
5262  dst,
5263  zero_dst_vector,
5264  *owning_class,
5265  function_pointer,
5266  nullptr,
5267  nullptr,
5268  src_vector_face_access_temp,
5270  task_info.loop(worker);
5271 }
5272 
5273 
5274 
5275 template <int dim, typename Number, typename VectorizedArrayType>
5276 template <typename CLASS, typename OutVector, typename InVector>
5277 inline void
5279  void (CLASS::*function_pointer)(
5281  OutVector &,
5282  const InVector &,
5283  const std::pair<unsigned int, unsigned int> &),
5284  CLASS * owning_class,
5285  OutVector & dst,
5286  const InVector & src,
5287  const bool zero_dst_vector,
5288  const DataAccessOnFaces src_vector_face_access) const
5289 {
5290  auto src_vector_face_access_temp = src_vector_face_access;
5291  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5292  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5293  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5294  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5295 
5296  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5297  InVector,
5298  OutVector,
5299  CLASS,
5300  false>
5301  worker(*this,
5302  src,
5303  dst,
5304  zero_dst_vector,
5305  *owning_class,
5306  function_pointer,
5307  nullptr,
5308  nullptr,
5309  src_vector_face_access_temp,
5311  task_info.loop(worker);
5312 }
5313 
5314 
5315 
5316 template <int dim, typename Number, typename VectorizedArrayType>
5317 template <typename OutVector, typename InVector>
5318 inline void
5320  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5321  OutVector &,
5322  const InVector &,
5323  const std::pair<unsigned int, unsigned int> &)>
5324  & cell_operation,
5325  OutVector & dst,
5326  const InVector & src,
5327  const bool zero_dst_vector,
5328  const DataAccessOnFaces src_vector_face_access) const
5329 {
5330  auto src_vector_face_access_temp = src_vector_face_access;
5331  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5332  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5333  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5334  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5335 
5336  using Wrapper =
5337  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5338  InVector,
5339  OutVector>;
5340  Wrapper wrap(cell_operation, nullptr, nullptr);
5341 
5342  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5343  InVector,
5344  OutVector,
5345  Wrapper,
5346  true>
5347  worker(*this,
5348  src,
5349  dst,
5350  zero_dst_vector,
5351  wrap,
5352  &Wrapper::cell_integrator,
5353  &Wrapper::face_integrator,
5354  &Wrapper::boundary_integrator,
5355  src_vector_face_access_temp,
5357  task_info.loop(worker);
5358 }
5359 
5360 
5361 #endif // ifndef DOXYGEN
5362 
5363 
5364 
5366 
5367 #endif
Transformed quadrature weights.
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > > shape_info
Definition: matrix_free.h:2120
void print_memory_consumption(StreamType &out) 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
static const unsigned int invalid_unsigned_int
Definition: types.h:196
unsigned int n_components_filled(const unsigned int cell_batch_number) const
std::vector< std::pair< unsigned int, unsigned int > > vector_zero_range_list
Definition: dof_info.h:714
AlignedVector< Number > * acquire_scratch_data_non_threadsafe() const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
~MatrixFree() override=default
unsigned int n_boundary_face_batches() const
A class that provides a separate storage location on each thread that accesses the object...
const Quadrature< dim - 1 > & get_face_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int get_dofs_per_cell(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
std::vector< unsigned int > cell_loop_pre_list_index
Definition: dof_info.h:721
unsigned int get_mg_level() const
unsigned int n_ghost_cell_batches() const
unsigned int n_components() 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
unsigned int get_n_q_points_face(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
void loop(MFWorkerInterface &worker) const
Definition: task_info.cc:330
unsigned int n_inner_face_batches() const
bool mapping_is_initialized
Definition: matrix_free.h:2161
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
AdditionalData(const AdditionalData &other)
Definition: matrix_free.h:253
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition: dof_info.h:565
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
std::vector< unsigned int > vector_zero_range_list_index
Definition: dof_info.h:709
std::list< std::pair< bool, AlignedVector< Number > > > scratch_pad_non_threadsafe
Definition: matrix_free.h:2179
std::pair< unsigned int, unsigned int > get_face_category(const unsigned int macro_face) const
bool mapping_initialized() const
void reinit(const size_type size, const bool omit_zeroing_entries=false)
Number value_type
Definition: matrix_free.h:127
static const char U
UpdateFlags mapping_update_flags_boundary_faces
Definition: matrix_free.h:389
static ::ExceptionBase & ExcNotInitialized()
std::vector< unsigned int > cell_partition_data
Definition: task_info.h:474
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
unsigned int get_n_q_points(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
UpdateFlags mapping_update_flags
Definition: matrix_free.h:368
std::size_t memory_consumption() const
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
const Number * constraint_pool_end(const unsigned int pool_index) const
void clear()
const IndexSet & get_ghost_set(const unsigned int dof_handler_index=0) const
::Table< 3, unsigned int > cell_and_face_to_plain_faces
Definition: face_info.h:172
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:408
static const unsigned int dimension
Definition: matrix_free.h:133
std::vector< unsigned int > cell_loop_post_list_index
Definition: dof_info.h:734
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
void print(std::ostream &out) const
const Number * constraint_pool_begin(const unsigned int pool_index) const
void reinit(const MappingType &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< number2 > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
static ::ExceptionBase & ExcMessage(std::string arg1)
bool indices_initialized() const
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() 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
const IndexSet & get_locally_owned_set(const unsigned int dof_handler_index=0) const
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
No update.
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_base_elements(const unsigned int dof_handler_index) const
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_pre_list
Definition: dof_info.h:727
unsigned int n_ghost_inner_face_batches() const
#define Assert(cond, exc)
Definition: exceptions.h:1465
UpdateFlags
void apply_to_subranges(const tbb::blocked_range< RangeType > &range, const Function &f)
Definition: parallel.h:367
internal::MatrixFreeFunctions::FaceInfo< VectorizedArrayType::size()> face_info
Definition: matrix_free.h:2151
const types::boundary_id invalid_boundary_id
Definition: types.h:239
unsigned int mg_level
Definition: matrix_free.h:2184
std::vector< std::pair< unsigned int, unsigned int > > cell_level_index
Definition: matrix_free.h:2128
unsigned int get_cell_active_fe_index(const std::pair< unsigned int, unsigned int > range) const
UpdateFlags mapping_update_flags_inner_faces
Definition: matrix_free.h:410
std::vector< unsigned int > boundary_partition_data
Definition: task_info.h:518
std::vector< internal::MatrixFreeFunctions::DoFInfo > dof_info
Definition: matrix_free.h:2093
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
unsigned int level
Definition: grid_out.cc:4590
std::pair< int, int > get_cell_level_and_index(const unsigned int cell_batch_index, const unsigned int lane_index) const
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_index) const
unsigned int get_dofs_per_face(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
UpdateFlags mapping_update_flags_faces_by_cells
Definition: matrix_free.h:438
unsigned int & level_mg_handler
Definition: matrix_free.h:455
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
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
unsigned int n_constraint_pool_entries() const
TasksParallelScheme tasks_parallel_scheme
Definition: matrix_free.h:344
std::vector< unsigned int > cell_vectorization_category
Definition: matrix_free.h:524
unsigned int n_cell_batches() const
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner(const unsigned int dof_handler_index=0) const
std::vector< FaceToCellTopology< vectorization_width > > faces
Definition: face_info.h:165
std::vector< Number > constraint_pool_data
Definition: matrix_free.h:2101
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
unsigned int n_active_fe_indices() 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)
types::boundary_id get_boundary_id(const unsigned int macro_face) const
std::vector< unsigned int > face_partition_data
Definition: task_info.h:496
unsigned int n_macro_cells() const
::Table< 3, types::boundary_id > cell_and_face_boundary_id
Definition: face_info.h:178
const std::vector< unsigned int > & get_constrained_dofs(const unsigned int dof_handler_index=0) const
internal::MatrixFreeFunctions::TaskInfo task_info
Definition: matrix_free.h:2144
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)
Definition: matrix_free.h:217
void initialize_dof_handlers(const std::vector< const DoFHandler< dim, dim > *> &dof_handlers, const AdditionalData &additional_data)
unsigned int get_cell_category(const unsigned int cell_batch_index) const
void update_mapping(const Mapping< dim > &mapping)
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
void initialize_dof_vector(VectorType &vec, const unsigned int dof_handler_index=0) const
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:445
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
VectorType::value_type * begin(VectorType &V)
Shape function gradients.
bool indices_are_initialized
Definition: matrix_free.h:2156
void renumber_dofs(std::vector< types::global_dof_index > &renumbering, const unsigned int dof_handler_index=0)
Threads::ThreadLocalStorage< std::list< std::pair< bool, AlignedVector< VectorizedArrayType > > > > scratch_pad
Definition: matrix_free.h:2172
unsigned int minimum_parallel_grain_size
Definition: parallel.cc:34
VectorizedArrayType vectorized_value_type
Definition: matrix_free.h:128
internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > mapping_info
Definition: matrix_free.h:2114
AdditionalData & operator=(const AdditionalData &other)
Definition: matrix_free.h:283
static ::ExceptionBase & ExcNotImplemented()
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 release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
void release_scratch_data_non_threadsafe(const AlignedVector< Number > *memory) const
std::vector< unsigned int > constraint_pool_row_index
Definition: matrix_free.h:2107
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_post_list
Definition: dof_info.h:740
static bool is_supported(const FiniteElement< dim, spacedim > &fe)
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
std::vector< SmartPointer< const DoFHandler< dim > > > dof_handlers
Definition: matrix_free.h:2087
unsigned int tasks_block_size
Definition: matrix_free.h:355
const Quadrature< dim > & get_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int get_face_active_fe_index(const std::pair< unsigned int, unsigned int > range, const bool is_interior_face=true) const
bool job_supports_mpi()
Definition: mpi.cc:1097
#define DEAL_II_DEPRECATED
Definition: config.h:157
unsigned int boundary_id
Definition: types.h:129
unsigned int cell_level_index_end_local
Definition: matrix_free.h:2137
unsigned int n_physical_cells() const
bool at_irregular_cell(const unsigned int cell_batch_index) const
void copy_from(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free_base)
static bool equal(const T *p1, const T *p2)
static ::ExceptionBase & ExcInternalError()