Reference documentation for deal.II version Git 0b4b57890a 2021-10-16 12:51:50 +0200
\(\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 
216  const unsigned int tasks_block_size = 0,
222  const unsigned int mg_level = numbers::invalid_unsigned_int,
223  const bool store_plain_indices = true,
224  const bool initialize_indices = true,
225  const bool initialize_mapping = true,
226  const bool overlap_communication_computation = true,
227  const bool hold_all_faces_to_owned_cells = false,
228  const bool cell_vectorization_categories_strict = false,
229  const bool allow_ghosted_vectors_in_loops = true,
230  const bool use_fast_hanging_node_algorithm = true)
237  , mg_level(mg_level)
247  , communicator_sm(MPI_COMM_SELF)
248  {}
249 
262  , mg_level(other.mg_level)
275  {}
276 
281  operator=(const AdditionalData &other)
282  {
291  mg_level = other.mg_level;
304 
305  return *this;
306  }
307 
345 
355  unsigned int tasks_block_size;
356 
369 
389 
409 
437 
446  unsigned int mg_level;
447 
455 
466 
475 
489 
498 
515  std::vector<unsigned int> cell_vectorization_category;
516 
527 
539 
544 
549  };
550 
558  MatrixFree();
559 
564 
568  ~MatrixFree() override = default;
569 
582  template <typename QuadratureType, typename number2, typename MappingType>
583  void
584  reinit(const MappingType & mapping,
585  const DoFHandler<dim> & dof_handler,
586  const AffineConstraints<number2> &constraint,
587  const QuadratureType & quad,
588  const AdditionalData & additional_data = AdditionalData());
589 
596  template <typename QuadratureType, typename number2>
597  DEAL_II_DEPRECATED void
598  reinit(const DoFHandler<dim> & dof_handler,
599  const AffineConstraints<number2> &constraint,
600  const QuadratureType & quad,
601  const AdditionalData & additional_data = AdditionalData());
602 
624  template <typename QuadratureType, typename number2, typename MappingType>
625  void
626  reinit(const MappingType & mapping,
627  const std::vector<const DoFHandler<dim> *> & dof_handler,
628  const std::vector<const AffineConstraints<number2> *> &constraint,
629  const std::vector<QuadratureType> & quad,
630  const AdditionalData &additional_data = AdditionalData());
631 
637  template <typename QuadratureType,
638  typename number2,
639  typename DoFHandlerType,
640  typename MappingType>
641  DEAL_II_DEPRECATED void
642  reinit(const MappingType & mapping,
643  const std::vector<const DoFHandlerType *> & dof_handler,
644  const std::vector<const AffineConstraints<number2> *> &constraint,
645  const std::vector<QuadratureType> & quad,
646  const AdditionalData &additional_data = AdditionalData());
647 
654  template <typename QuadratureType, typename number2>
655  DEAL_II_DEPRECATED void
656  reinit(const std::vector<const DoFHandler<dim> *> & dof_handler,
657  const std::vector<const AffineConstraints<number2> *> &constraint,
658  const std::vector<QuadratureType> & quad,
659  const AdditionalData &additional_data = AdditionalData());
660 
666  template <typename QuadratureType, typename number2, typename DoFHandlerType>
667  DEAL_II_DEPRECATED void
668  reinit(const std::vector<const DoFHandlerType *> & dof_handler,
669  const std::vector<const AffineConstraints<number2> *> &constraint,
670  const std::vector<QuadratureType> & quad,
671  const AdditionalData &additional_data = AdditionalData());
672 
680  template <typename QuadratureType, typename number2, typename MappingType>
681  void
682  reinit(const MappingType & mapping,
683  const std::vector<const DoFHandler<dim> *> & dof_handler,
684  const std::vector<const AffineConstraints<number2> *> &constraint,
685  const QuadratureType & quad,
686  const AdditionalData &additional_data = AdditionalData());
687 
693  template <typename QuadratureType,
694  typename number2,
695  typename DoFHandlerType,
696  typename MappingType>
697  DEAL_II_DEPRECATED void
698  reinit(const MappingType & mapping,
699  const std::vector<const DoFHandlerType *> & dof_handler,
700  const std::vector<const AffineConstraints<number2> *> &constraint,
701  const QuadratureType & quad,
702  const AdditionalData &additional_data = AdditionalData());
703 
710  template <typename QuadratureType, typename number2>
711  DEAL_II_DEPRECATED void
712  reinit(const std::vector<const DoFHandler<dim> *> & dof_handler,
713  const std::vector<const AffineConstraints<number2> *> &constraint,
714  const QuadratureType & quad,
715  const AdditionalData &additional_data = AdditionalData());
716 
722  template <typename QuadratureType, typename number2, typename DoFHandlerType>
723  DEAL_II_DEPRECATED void
724  reinit(const std::vector<const DoFHandlerType *> & dof_handler,
725  const std::vector<const AffineConstraints<number2> *> &constraint,
726  const QuadratureType & quad,
727  const AdditionalData &additional_data = AdditionalData());
728 
734  void
735  copy_from(
736  const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free_base);
737 
747  void
748  update_mapping(const Mapping<dim> &mapping);
749 
753  void
754  update_mapping(const std::shared_ptr<hp::MappingCollection<dim>> &mapping);
755 
760  void
761  clear();
762 
764 
776  enum class DataAccessOnFaces
777  {
784  none,
785 
796  values,
797 
806  values_all_faces,
807 
818  gradients,
819 
828  gradients_all_faces,
829 
836  unspecified
837  };
838 
883  template <typename OutVector, typename InVector>
884  void
885  cell_loop(const std::function<void(
887  OutVector &,
888  const InVector &,
889  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
890  OutVector & dst,
891  const InVector & src,
892  const bool zero_dst_vector = false) const;
893 
940  template <typename CLASS, typename OutVector, typename InVector>
941  void
942  cell_loop(void (CLASS::*cell_operation)(
943  const MatrixFree &,
944  OutVector &,
945  const InVector &,
946  const std::pair<unsigned int, unsigned int> &) const,
947  const CLASS * owning_class,
948  OutVector & dst,
949  const InVector &src,
950  const bool zero_dst_vector = false) const;
951 
955  template <typename CLASS, typename OutVector, typename InVector>
956  void
957  cell_loop(void (CLASS::*cell_operation)(
958  const MatrixFree &,
959  OutVector &,
960  const InVector &,
961  const std::pair<unsigned int, unsigned int> &),
962  CLASS * owning_class,
963  OutVector & dst,
964  const InVector &src,
965  const bool zero_dst_vector = false) const;
966 
1051  template <typename CLASS, typename OutVector, typename InVector>
1052  void
1053  cell_loop(void (CLASS::*cell_operation)(
1054  const MatrixFree &,
1055  OutVector &,
1056  const InVector &,
1057  const std::pair<unsigned int, unsigned int> &) const,
1058  const CLASS * owning_class,
1059  OutVector & dst,
1060  const InVector &src,
1061  const std::function<void(const unsigned int, const unsigned int)>
1062  &operation_before_loop,
1063  const std::function<void(const unsigned int, const unsigned int)>
1064  & operation_after_loop,
1065  const unsigned int dof_handler_index_pre_post = 0) const;
1066 
1070  template <typename CLASS, typename OutVector, typename InVector>
1071  void
1072  cell_loop(void (CLASS::*cell_operation)(
1073  const MatrixFree &,
1074  OutVector &,
1075  const InVector &,
1076  const std::pair<unsigned int, unsigned int> &),
1077  CLASS * owning_class,
1078  OutVector & dst,
1079  const InVector &src,
1080  const std::function<void(const unsigned int, const unsigned int)>
1081  &operation_before_loop,
1082  const std::function<void(const unsigned int, const unsigned int)>
1083  & operation_after_loop,
1084  const unsigned int dof_handler_index_pre_post = 0) const;
1085 
1090  template <typename OutVector, typename InVector>
1091  void
1092  cell_loop(const std::function<void(
1094  OutVector &,
1095  const InVector &,
1096  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1097  OutVector & dst,
1098  const InVector & src,
1099  const std::function<void(const unsigned int, const unsigned int)>
1100  &operation_before_loop,
1101  const std::function<void(const unsigned int, const unsigned int)>
1102  & operation_after_loop,
1103  const unsigned int dof_handler_index_pre_post = 0) const;
1104 
1180  template <typename OutVector, typename InVector>
1181  void
1182  loop(const std::function<
1184  OutVector &,
1185  const InVector &,
1186  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1187  const std::function<
1189  OutVector &,
1190  const InVector &,
1191  const std::pair<unsigned int, unsigned int> &)> &face_operation,
1192  const std::function<void(
1194  OutVector &,
1195  const InVector &,
1196  const std::pair<unsigned int, unsigned int> &)> &boundary_operation,
1197  OutVector & dst,
1198  const InVector & src,
1199  const bool zero_dst_vector = false,
1200  const DataAccessOnFaces dst_vector_face_access =
1202  const DataAccessOnFaces src_vector_face_access =
1204 
1290  template <typename CLASS, typename OutVector, typename InVector>
1291  void
1292  loop(
1293  void (CLASS::*cell_operation)(const MatrixFree &,
1294  OutVector &,
1295  const InVector &,
1296  const std::pair<unsigned int, unsigned int> &)
1297  const,
1298  void (CLASS::*face_operation)(const MatrixFree &,
1299  OutVector &,
1300  const InVector &,
1301  const std::pair<unsigned int, unsigned int> &)
1302  const,
1303  void (CLASS::*boundary_operation)(
1304  const MatrixFree &,
1305  OutVector &,
1306  const InVector &,
1307  const std::pair<unsigned int, unsigned int> &) const,
1308  const CLASS * owning_class,
1309  OutVector & dst,
1310  const InVector & src,
1311  const bool zero_dst_vector = false,
1312  const DataAccessOnFaces dst_vector_face_access =
1314  const DataAccessOnFaces src_vector_face_access =
1316 
1320  template <typename CLASS, typename OutVector, typename InVector>
1321  void
1322  loop(void (CLASS::*cell_operation)(
1323  const MatrixFree &,
1324  OutVector &,
1325  const InVector &,
1326  const std::pair<unsigned int, unsigned int> &),
1327  void (CLASS::*face_operation)(
1328  const MatrixFree &,
1329  OutVector &,
1330  const InVector &,
1331  const std::pair<unsigned int, unsigned int> &),
1332  void (CLASS::*boundary_operation)(
1333  const MatrixFree &,
1334  OutVector &,
1335  const InVector &,
1336  const std::pair<unsigned int, unsigned int> &),
1337  CLASS * owning_class,
1338  OutVector & dst,
1339  const InVector & src,
1340  const bool zero_dst_vector = false,
1341  const DataAccessOnFaces dst_vector_face_access =
1343  const DataAccessOnFaces src_vector_face_access =
1345 
1410  template <typename CLASS, typename OutVector, typename InVector>
1411  void
1412  loop_cell_centric(void (CLASS::*cell_operation)(
1413  const MatrixFree &,
1414  OutVector &,
1415  const InVector &,
1416  const std::pair<unsigned int, unsigned int> &) const,
1417  const CLASS * owning_class,
1418  OutVector & dst,
1419  const InVector & src,
1420  const bool zero_dst_vector = false,
1421  const DataAccessOnFaces src_vector_face_access =
1423 
1427  template <typename CLASS, typename OutVector, typename InVector>
1428  void
1429  loop_cell_centric(void (CLASS::*cell_operation)(
1430  const MatrixFree &,
1431  OutVector &,
1432  const InVector &,
1433  const std::pair<unsigned int, unsigned int> &),
1434  CLASS * owning_class,
1435  OutVector & dst,
1436  const InVector & src,
1437  const bool zero_dst_vector = false,
1438  const DataAccessOnFaces src_vector_face_access =
1440 
1444  template <typename OutVector, typename InVector>
1445  void
1447  const std::function<void(const MatrixFree &,
1448  OutVector &,
1449  const InVector &,
1450  const std::pair<unsigned int, unsigned int> &)>
1451  & cell_operation,
1452  OutVector & dst,
1453  const InVector & src,
1454  const bool zero_dst_vector = false,
1455  const DataAccessOnFaces src_vector_face_access =
1457 
1465  std::pair<unsigned int, unsigned int>
1466  create_cell_subrange_hp(const std::pair<unsigned int, unsigned int> &range,
1467  const unsigned int fe_degree,
1468  const unsigned int dof_handler_index = 0) const;
1469 
1476  std::pair<unsigned int, unsigned int>
1478  const std::pair<unsigned int, unsigned int> &range,
1479  const unsigned int fe_index,
1480  const unsigned int dof_handler_index = 0) const;
1481 
1485  unsigned int
1486  n_active_fe_indices() const;
1487 
1491  unsigned int
1493  const std::pair<unsigned int, unsigned int> range) const;
1494 
1498  unsigned int
1499  get_face_active_fe_index(const std::pair<unsigned int, unsigned int> range,
1500  const bool is_interior_face = true) const;
1501 
1503 
1528  template <typename VectorType>
1529  void
1531  const unsigned int dof_handler_index = 0) const;
1532 
1553  template <typename Number2>
1554  void
1556  const unsigned int dof_handler_index = 0) const;
1557 
1568  const std::shared_ptr<const Utilities::MPI::Partitioner> &
1569  get_vector_partitioner(const unsigned int dof_handler_index = 0) const;
1570 
1574  const IndexSet &
1575  get_locally_owned_set(const unsigned int dof_handler_index = 0) const;
1576 
1580  const IndexSet &
1581  get_ghost_set(const unsigned int dof_handler_index = 0) const;
1582 
1592  const std::vector<unsigned int> &
1593  get_constrained_dofs(const unsigned int dof_handler_index = 0) const;
1594 
1605  void
1606  renumber_dofs(std::vector<types::global_dof_index> &renumbering,
1607  const unsigned int dof_handler_index = 0);
1608 
1610 
1618  template <int spacedim>
1619  static bool
1621 
1625  unsigned int
1626  n_components() const;
1627 
1632  unsigned int
1633  n_base_elements(const unsigned int dof_handler_index) const;
1634 
1642  unsigned int
1643  n_physical_cells() const;
1644 
1648  DEAL_II_DEPRECATED unsigned int
1649  n_macro_cells() const;
1650 
1660  unsigned int
1661  n_cell_batches() const;
1662 
1669  unsigned int
1670  n_ghost_cell_batches() const;
1671 
1681  unsigned int
1682  n_inner_face_batches() const;
1683 
1694  unsigned int
1695  n_boundary_face_batches() const;
1696 
1701  unsigned int
1703 
1711  get_boundary_id(const unsigned int macro_face) const;
1712 
1717  std::array<types::boundary_id, VectorizedArrayType::size()>
1718  get_faces_by_cells_boundary_id(const unsigned int cell_batch_index,
1719  const unsigned int face_number) const;
1720 
1725  const DoFHandler<dim> &
1726  get_dof_handler(const unsigned int dof_handler_index = 0) const;
1727 
1738  template <typename DoFHandlerType>
1739  DEAL_II_DEPRECATED const DoFHandlerType &
1740  get_dof_handler(const unsigned int dof_handler_index = 0) const;
1741 
1755  get_cell_iterator(const unsigned int cell_batch_index,
1756  const unsigned int lane_index,
1757  const unsigned int dof_handler_index = 0) const;
1758 
1764  std::pair<int, int>
1765  get_cell_level_and_index(const unsigned int cell_batch_index,
1766  const unsigned int lane_index) const;
1767 
1780  std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>
1781  get_face_iterator(const unsigned int face_batch_index,
1782  const unsigned int lane_index,
1783  const bool interior = true,
1784  const unsigned int fe_component = 0) const;
1785 
1792  get_hp_cell_iterator(const unsigned int cell_batch_index,
1793  const unsigned int lane_index,
1794  const unsigned int dof_handler_index = 0) const;
1795 
1808  bool
1809  at_irregular_cell(const unsigned int cell_batch_index) const;
1810 
1814  DEAL_II_DEPRECATED unsigned int
1815  n_components_filled(const unsigned int cell_batch_number) const;
1816 
1826  unsigned int
1827  n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const;
1828 
1838  unsigned int
1839  n_active_entries_per_face_batch(const unsigned int face_batch_index) const;
1840 
1844  unsigned int
1845  get_dofs_per_cell(const unsigned int dof_handler_index = 0,
1846  const unsigned int hp_active_fe_index = 0) const;
1847 
1851  unsigned int
1852  get_n_q_points(const unsigned int quad_index = 0,
1853  const unsigned int hp_active_fe_index = 0) const;
1854 
1859  unsigned int
1860  get_dofs_per_face(const unsigned int dof_handler_index = 0,
1861  const unsigned int hp_active_fe_index = 0) const;
1862 
1867  unsigned int
1868  get_n_q_points_face(const unsigned int quad_index = 0,
1869  const unsigned int hp_active_fe_index = 0) const;
1870 
1874  const Quadrature<dim> &
1875  get_quadrature(const unsigned int quad_index = 0,
1876  const unsigned int hp_active_fe_index = 0) const;
1877 
1881  const Quadrature<dim - 1> &
1882  get_face_quadrature(const unsigned int quad_index = 0,
1883  const unsigned int hp_active_fe_index = 0) const;
1884 
1891  unsigned int
1892  get_cell_category(const unsigned int cell_batch_index) const;
1893 
1898  std::pair<unsigned int, unsigned int>
1899  get_face_category(const unsigned int macro_face) const;
1900 
1904  bool
1905  indices_initialized() const;
1906 
1911  bool
1912  mapping_initialized() const;
1913 
1918  unsigned int
1919  get_mg_level() const;
1920 
1925  std::size_t
1926  memory_consumption() const;
1927 
1932  template <typename StreamType>
1933  void
1934  print_memory_consumption(StreamType &out) const;
1935 
1940  void
1941  print(std::ostream &out) const;
1942 
1944 
1955  get_task_info() const;
1956 
1957  /*
1958  * Return geometry-dependent information on the cells.
1959  */
1960  const internal::MatrixFreeFunctions::
1961  MappingInfo<dim, Number, VectorizedArrayType> &
1962  get_mapping_info() const;
1963 
1968  get_dof_info(const unsigned int dof_handler_index_component = 0) const;
1969 
1973  unsigned int
1974  n_constraint_pool_entries() const;
1975 
1980  const Number *
1981  constraint_pool_begin(const unsigned int pool_index) const;
1982 
1988  const Number *
1989  constraint_pool_end(const unsigned int pool_index) const;
1990 
1995  get_shape_info(const unsigned int dof_handler_index_component = 0,
1996  const unsigned int quad_index = 0,
1997  const unsigned int fe_base_element = 0,
1998  const unsigned int hp_active_fe_index = 0,
1999  const unsigned int hp_active_quad_index = 0) const;
2000 
2005  VectorizedArrayType::size()> &
2006  get_face_info(const unsigned int face_batch_index) const;
2007 
2008 
2014  const Table<3, unsigned int> &
2016 
2031  acquire_scratch_data() const;
2032 
2036  void
2038 
2050 
2054  void
2056  const AlignedVector<Number> *memory) const;
2057 
2059 
2060 private:
2065  template <typename number2, int q_dim>
2066  void
2068  const std::shared_ptr<hp::MappingCollection<dim>> & mapping,
2069  const std::vector<const DoFHandler<dim, dim> *> & dof_handlers,
2070  const std::vector<const AffineConstraints<number2> *> &constraint,
2071  const std::vector<IndexSet> & locally_owned_set,
2072  const std::vector<hp::QCollection<q_dim>> & quad,
2073  const AdditionalData & additional_data);
2074 
2081  template <typename number2>
2082  void
2084  const std::vector<const AffineConstraints<number2> *> &constraint,
2085  const std::vector<IndexSet> & locally_owned_set,
2086  const AdditionalData & additional_data);
2087 
2091  void
2093  const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
2094  const AdditionalData & additional_data);
2095 
2099  std::vector<SmartPointer<const DoFHandler<dim>>> dof_handlers;
2100 
2105  std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
2106 
2113  std::vector<Number> constraint_pool_data;
2114 
2119  std::vector<unsigned int> constraint_pool_row_index;
2120 
2127 
2133 
2140  std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
2141 
2142 
2150 
2157 
2164 
2169 
2174 
2183  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2185 
2190  mutable std::list<std::pair<bool, AlignedVector<Number>>>
2192 
2196  unsigned int mg_level;
2197 };
2198 
2199 
2200 
2201 /*----------------------- Inline functions ----------------------------------*/
2202 
2203 #ifndef DOXYGEN
2204 
2205 
2206 
2207 template <int dim, typename Number, typename VectorizedArrayType>
2208 template <typename VectorType>
2209 inline void
2211  VectorType & vec,
2212  const unsigned int comp) const
2213 {
2214  AssertIndexRange(comp, n_components());
2215  vec.reinit(dof_info[comp].vector_partitioner->size());
2216 }
2217 
2218 
2219 
2220 template <int dim, typename Number, typename VectorizedArrayType>
2221 template <typename Number2>
2222 inline void
2225  const unsigned int comp) const
2226 {
2227  AssertIndexRange(comp, n_components());
2228  vec.reinit(dof_info[comp].vector_partitioner, task_info.communicator_sm);
2229 }
2230 
2231 
2232 
2233 template <int dim, typename Number, typename VectorizedArrayType>
2234 inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2236  const unsigned int comp) const
2237 {
2238  AssertIndexRange(comp, n_components());
2239  return dof_info[comp].vector_partitioner;
2240 }
2241 
2242 
2243 
2244 template <int dim, typename Number, typename VectorizedArrayType>
2245 inline const std::vector<unsigned int> &
2247  const unsigned int comp) const
2248 {
2249  AssertIndexRange(comp, n_components());
2250  return dof_info[comp].constrained_dofs;
2251 }
2252 
2253 
2254 
2255 template <int dim, typename Number, typename VectorizedArrayType>
2256 inline unsigned int
2258 {
2259  AssertDimension(dof_handlers.size(), dof_info.size());
2260  return dof_handlers.size();
2261 }
2262 
2263 
2264 
2265 template <int dim, typename Number, typename VectorizedArrayType>
2266 inline unsigned int
2268  const unsigned int dof_no) const
2269 {
2270  AssertDimension(dof_handlers.size(), dof_info.size());
2271  AssertIndexRange(dof_no, dof_handlers.size());
2272  return dof_handlers[dof_no]->get_fe().n_base_elements();
2273 }
2274 
2275 
2276 
2277 template <int dim, typename Number, typename VectorizedArrayType>
2280 {
2281  return task_info;
2282 }
2283 
2284 
2285 
2286 template <int dim, typename Number, typename VectorizedArrayType>
2287 inline unsigned int
2289 {
2290  return *(task_info.cell_partition_data.end() - 2);
2291 }
2292 
2293 
2294 
2295 template <int dim, typename Number, typename VectorizedArrayType>
2296 inline unsigned int
2298 {
2299  return task_info.n_active_cells;
2300 }
2301 
2302 
2303 
2304 template <int dim, typename Number, typename VectorizedArrayType>
2305 inline unsigned int
2307 {
2308  return *(task_info.cell_partition_data.end() - 2);
2309 }
2310 
2311 
2312 
2313 template <int dim, typename Number, typename VectorizedArrayType>
2314 inline unsigned int
2316 {
2317  return *(task_info.cell_partition_data.end() - 1) -
2318  *(task_info.cell_partition_data.end() - 2);
2319 }
2320 
2321 
2322 
2323 template <int dim, typename Number, typename VectorizedArrayType>
2324 inline unsigned int
2326 {
2327  if (task_info.face_partition_data.size() == 0)
2328  return 0;
2329  return task_info.face_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 task_info.boundary_partition_data.back() -
2342 }
2343 
2344 
2345 
2346 template <int dim, typename Number, typename VectorizedArrayType>
2347 inline unsigned int
2349 {
2350  if (task_info.face_partition_data.size() == 0)
2351  return 0;
2352  return face_info.faces.size() - task_info.boundary_partition_data.back();
2353 }
2354 
2355 
2356 
2357 template <int dim, typename Number, typename VectorizedArrayType>
2358 inline types::boundary_id
2360  const unsigned int macro_face) const
2361 {
2362  Assert(macro_face >= task_info.boundary_partition_data[0] &&
2363  macro_face < task_info.boundary_partition_data.back(),
2364  ExcIndexRange(macro_face,
2367  return types::boundary_id(face_info.faces[macro_face].exterior_face_no);
2368 }
2369 
2370 
2371 
2372 template <int dim, typename Number, typename VectorizedArrayType>
2373 inline std::array<types::boundary_id, VectorizedArrayType::size()>
2375  const unsigned int cell_batch_index,
2376  const unsigned int face_number) const
2377 {
2378  AssertIndexRange(cell_batch_index, n_cell_batches());
2381  ExcNotInitialized());
2382  std::array<types::boundary_id, VectorizedArrayType::size()> result;
2383  result.fill(numbers::invalid_boundary_id);
2384  for (unsigned int v = 0;
2385  v < n_active_entries_per_cell_batch(cell_batch_index);
2386  ++v)
2387  result[v] =
2388  face_info.cell_and_face_boundary_id(cell_batch_index, face_number, v);
2389  return result;
2390 }
2391 
2392 
2393 
2394 template <int dim, typename Number, typename VectorizedArrayType>
2395 inline const internal::MatrixFreeFunctions::
2396  MappingInfo<dim, Number, VectorizedArrayType> &
2398 {
2399  return mapping_info;
2400 }
2401 
2402 
2403 
2404 template <int dim, typename Number, typename VectorizedArrayType>
2407  const unsigned int dof_index) const
2408 {
2409  AssertIndexRange(dof_index, n_components());
2410  return dof_info[dof_index];
2411 }
2412 
2413 
2414 
2415 template <int dim, typename Number, typename VectorizedArrayType>
2416 inline unsigned int
2418 {
2419  return constraint_pool_row_index.size() - 1;
2420 }
2421 
2422 
2423 
2424 template <int dim, typename Number, typename VectorizedArrayType>
2425 inline const Number *
2427  const unsigned int row) const
2428 {
2430  return constraint_pool_data.empty() ?
2431  nullptr :
2433 }
2434 
2435 
2436 
2437 template <int dim, typename Number, typename VectorizedArrayType>
2438 inline const Number *
2440  const unsigned int row) const
2441 {
2443  return constraint_pool_data.empty() ?
2444  nullptr :
2446 }
2447 
2448 
2449 
2450 template <int dim, typename Number, typename VectorizedArrayType>
2451 inline std::pair<unsigned int, unsigned int>
2453  const std::pair<unsigned int, unsigned int> &range,
2454  const unsigned int degree,
2455  const unsigned int dof_handler_component) const
2456 {
2457  if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2458  {
2460  dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2462  dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2463  if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2464  return range;
2465  else
2466  return {range.second, range.second};
2467  }
2468 
2469  const unsigned int fe_index =
2470  dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2471  if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2472  return {range.second, range.second};
2473  else
2474  return create_cell_subrange_hp_by_index(range,
2475  fe_index,
2476  dof_handler_component);
2477 }
2478 
2479 
2480 
2481 template <int dim, typename Number, typename VectorizedArrayType>
2482 inline bool
2484  const unsigned int cell_batch_index) const
2485 {
2486  AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2487  return VectorizedArrayType::size() > 1 &&
2488  cell_level_index[(cell_batch_index + 1) * VectorizedArrayType::size() -
2489  1] == cell_level_index[(cell_batch_index + 1) *
2490  VectorizedArrayType::size() -
2491  2];
2492 }
2493 
2494 
2495 
2496 template <int dim, typename Number, typename VectorizedArrayType>
2497 unsigned int
2499 {
2500  return shape_info.size(2);
2501 }
2502 
2503 
2504 template <int dim, typename Number, typename VectorizedArrayType>
2505 unsigned int
2507  const std::pair<unsigned int, unsigned int> range) const
2508 {
2509  const auto &fe_indices = dof_info[0].cell_active_fe_index;
2510 
2511  if (fe_indices.empty() == true)
2512  return 0;
2513 
2514  const auto index = fe_indices[range.first];
2515 
2516  for (unsigned int i = range.first; i < range.second; ++i)
2517  AssertDimension(index, fe_indices[i]);
2518 
2519  return index;
2520 }
2521 
2522 
2523 
2524 template <int dim, typename Number, typename VectorizedArrayType>
2525 unsigned int
2527  const std::pair<unsigned int, unsigned int> range,
2528  const bool is_interior_face) const
2529 {
2530  const auto &fe_indices = dof_info[0].cell_active_fe_index;
2531 
2532  if (fe_indices.empty() == true)
2533  return 0;
2534 
2535  if (is_interior_face)
2536  {
2537  const unsigned int index =
2538  fe_indices[face_info.faces[range.first].cells_interior[0] /
2539  VectorizedArrayType::size()];
2540 
2541  for (unsigned int i = range.first; i < range.second; ++i)
2542  AssertDimension(index,
2543  fe_indices[face_info.faces[i].cells_interior[0] /
2544  VectorizedArrayType::size()]);
2545 
2546  return index;
2547  }
2548  else
2549  {
2550  const unsigned int index =
2551  fe_indices[face_info.faces[range.first].cells_exterior[0] /
2552  VectorizedArrayType::size()];
2553 
2554  for (unsigned int i = range.first; i < range.second; ++i)
2555  AssertDimension(index,
2556  fe_indices[face_info.faces[i].cells_exterior[0] /
2557  VectorizedArrayType::size()]);
2558 
2559  return index;
2560  }
2561 }
2562 
2563 
2564 
2565 template <int dim, typename Number, typename VectorizedArrayType>
2566 inline unsigned int
2568  const unsigned int cell_batch_index) const
2569 {
2570  return n_active_entries_per_cell_batch(cell_batch_index);
2571 }
2572 
2573 
2574 
2575 template <int dim, typename Number, typename VectorizedArrayType>
2576 inline unsigned int
2578  const unsigned int cell_batch_index) const
2579 {
2580  AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2581  unsigned int n_lanes = VectorizedArrayType::size();
2582  while (n_lanes > 1 &&
2583  cell_level_index[cell_batch_index * VectorizedArrayType::size() +
2584  n_lanes - 1] ==
2585  cell_level_index[cell_batch_index * VectorizedArrayType::size() +
2586  n_lanes - 2])
2587  --n_lanes;
2588  AssertIndexRange(n_lanes - 1, VectorizedArrayType::size());
2589  return n_lanes;
2590 }
2591 
2592 
2593 
2594 template <int dim, typename Number, typename VectorizedArrayType>
2595 inline unsigned int
2597  const unsigned int face_batch_index) const
2598 {
2599  AssertIndexRange(face_batch_index, face_info.faces.size());
2600  unsigned int n_lanes = VectorizedArrayType::size();
2601  while (n_lanes > 1 &&
2602  face_info.faces[face_batch_index].cells_interior[n_lanes - 1] ==
2604  --n_lanes;
2605  AssertIndexRange(n_lanes - 1, VectorizedArrayType::size());
2606  return n_lanes;
2607 }
2608 
2609 
2610 
2611 template <int dim, typename Number, typename VectorizedArrayType>
2612 inline unsigned int
2614  const unsigned int dof_handler_index,
2615  const unsigned int active_fe_index) const
2616 {
2617  return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2618 }
2619 
2620 
2621 
2622 template <int dim, typename Number, typename VectorizedArrayType>
2623 inline unsigned int
2625  const unsigned int quad_index,
2626  const unsigned int active_fe_index) const
2627 {
2628  AssertIndexRange(quad_index, mapping_info.cell_data.size());
2629  return mapping_info.cell_data[quad_index]
2630  .descriptor[active_fe_index]
2631  .n_q_points;
2632 }
2633 
2634 
2635 
2636 template <int dim, typename Number, typename VectorizedArrayType>
2637 inline unsigned int
2639  const unsigned int dof_handler_index,
2640  const unsigned int active_fe_index) const
2641 {
2642  return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2643 }
2644 
2645 
2646 
2647 template <int dim, typename Number, typename VectorizedArrayType>
2648 inline unsigned int
2650  const unsigned int quad_index,
2651  const unsigned int active_fe_index) const
2652 {
2653  AssertIndexRange(quad_index, mapping_info.face_data.size());
2654  return mapping_info.face_data[quad_index]
2655  .descriptor[active_fe_index]
2656  .n_q_points;
2657 }
2658 
2659 
2660 
2661 template <int dim, typename Number, typename VectorizedArrayType>
2662 inline const IndexSet &
2664  const unsigned int dof_handler_index) const
2665 {
2666  return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2667 }
2668 
2669 
2670 
2671 template <int dim, typename Number, typename VectorizedArrayType>
2672 inline const IndexSet &
2674  const unsigned int dof_handler_index) const
2675 {
2676  return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2677 }
2678 
2679 
2680 
2681 template <int dim, typename Number, typename VectorizedArrayType>
2684  const unsigned int dof_handler_index,
2685  const unsigned int index_quad,
2686  const unsigned int index_fe,
2687  const unsigned int active_fe_index,
2688  const unsigned int active_quad_index) const
2689 {
2690  AssertIndexRange(dof_handler_index, dof_info.size());
2691  const unsigned int ind =
2692  dof_info[dof_handler_index].global_base_element_offset + index_fe;
2693  AssertIndexRange(ind, shape_info.size(0));
2694  AssertIndexRange(index_quad, shape_info.size(1));
2695  AssertIndexRange(active_fe_index, shape_info.size(2));
2696  AssertIndexRange(active_quad_index, shape_info.size(3));
2697  return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2698 }
2699 
2700 
2701 
2702 template <int dim, typename Number, typename VectorizedArrayType>
2704  VectorizedArrayType::size()> &
2706  const unsigned int macro_face) const
2707 {
2708  AssertIndexRange(macro_face, face_info.faces.size());
2709  return face_info.faces[macro_face];
2710 }
2711 
2712 
2713 
2714 template <int dim, typename Number, typename VectorizedArrayType>
2715 inline const Table<3, unsigned int> &
2717  const
2718 {
2720 }
2721 
2722 
2723 
2724 template <int dim, typename Number, typename VectorizedArrayType>
2725 inline const Quadrature<dim> &
2727  const unsigned int quad_index,
2728  const unsigned int active_fe_index) const
2729 {
2730  AssertIndexRange(quad_index, mapping_info.cell_data.size());
2731  return mapping_info.cell_data[quad_index]
2732  .descriptor[active_fe_index]
2733  .quadrature;
2734 }
2735 
2736 
2737 
2738 template <int dim, typename Number, typename VectorizedArrayType>
2739 inline const Quadrature<dim - 1> &
2741  const unsigned int quad_index,
2742  const unsigned int active_fe_index) const
2743 {
2744  AssertIndexRange(quad_index, mapping_info.face_data.size());
2745  return mapping_info.face_data[quad_index]
2746  .descriptor[active_fe_index]
2747  .quadrature;
2748 }
2749 
2750 
2751 
2752 template <int dim, typename Number, typename VectorizedArrayType>
2753 inline unsigned int
2755  const unsigned int cell_batch_index) const
2756 {
2757  AssertIndexRange(0, dof_info.size());
2758  AssertIndexRange(cell_batch_index, dof_info[0].cell_active_fe_index.size());
2759  if (dof_info[0].cell_active_fe_index.empty())
2760  return 0;
2761  else
2762  return dof_info[0].cell_active_fe_index[cell_batch_index];
2763 }
2764 
2765 
2766 
2767 template <int dim, typename Number, typename VectorizedArrayType>
2768 inline std::pair<unsigned int, unsigned int>
2770  const unsigned int macro_face) const
2771 {
2772  AssertIndexRange(macro_face, face_info.faces.size());
2773  if (dof_info[0].cell_active_fe_index.empty())
2774  return std::make_pair(0U, 0U);
2775 
2776  std::pair<unsigned int, unsigned int> result;
2777  for (unsigned int v = 0; v < VectorizedArrayType::size() &&
2778  face_info.faces[macro_face].cells_interior[v] !=
2780  ++v)
2781  result.first = std::max(
2782  result.first,
2783  dof_info[0]
2784  .cell_active_fe_index[face_info.faces[macro_face].cells_interior[v]]);
2785  if (face_info.faces[macro_face].cells_exterior[0] !=
2787  for (unsigned int v = 0; v < VectorizedArrayType::size() &&
2788  face_info.faces[macro_face].cells_exterior[v] !=
2790  ++v)
2791  result.second = std::max(
2792  result.first,
2793  dof_info[0]
2794  .cell_active_fe_index[face_info.faces[macro_face].cells_exterior[v]]);
2795  else
2796  result.second = numbers::invalid_unsigned_int;
2797  return result;
2798 }
2799 
2800 
2801 
2802 template <int dim, typename Number, typename VectorizedArrayType>
2803 inline bool
2805 {
2806  return indices_are_initialized;
2807 }
2808 
2809 
2810 
2811 template <int dim, typename Number, typename VectorizedArrayType>
2812 inline bool
2814 {
2815  return mapping_is_initialized;
2816 }
2817 
2818 
2819 template <int dim, typename Number, typename VectorizedArrayType>
2820 inline unsigned int
2822 {
2823  return mg_level;
2824 }
2825 
2826 
2827 
2828 template <int dim, typename Number, typename VectorizedArrayType>
2831 {
2832  using list_type =
2833  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
2834  list_type &data = scratch_pad.get();
2835  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
2836  if (it->first == false)
2837  {
2838  it->first = true;
2839  return &it->second;
2840  }
2841  data.emplace_front(true, AlignedVector<VectorizedArrayType>());
2842  return &data.front().second;
2843 }
2844 
2845 
2846 
2847 template <int dim, typename Number, typename VectorizedArrayType>
2848 void
2850  const AlignedVector<VectorizedArrayType> *scratch) const
2851 {
2852  using list_type =
2853  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
2854  list_type &data = scratch_pad.get();
2855  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
2856  if (&it->second == scratch)
2857  {
2858  Assert(it->first == true, ExcInternalError());
2859  it->first = false;
2860  return;
2861  }
2862  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2863 }
2864 
2865 
2866 
2867 template <int dim, typename Number, typename VectorizedArrayType>
2871 {
2872  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
2874  it != scratch_pad_non_threadsafe.end();
2875  ++it)
2876  if (it->first == false)
2877  {
2878  it->first = true;
2879  return &it->second;
2880  }
2881  scratch_pad_non_threadsafe.push_front(
2882  std::make_pair(true, AlignedVector<Number>()));
2883  return &scratch_pad_non_threadsafe.front().second;
2884 }
2885 
2886 
2887 
2888 template <int dim, typename Number, typename VectorizedArrayType>
2889 void
2892  const AlignedVector<Number> *scratch) const
2893 {
2894  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
2896  it != scratch_pad_non_threadsafe.end();
2897  ++it)
2898  if (&it->second == scratch)
2899  {
2900  Assert(it->first == true, ExcInternalError());
2901  it->first = false;
2902  return;
2903  }
2904  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2905 }
2906 
2907 
2908 
2909 // ------------------------------ reinit functions ---------------------------
2910 
2911 namespace internal
2912 {
2913  namespace MatrixFreeImplementation
2914  {
2915  template <int dim, int spacedim>
2916  inline std::vector<IndexSet>
2917  extract_locally_owned_index_sets(
2918  const std::vector<const ::DoFHandler<dim, spacedim> *> &dofh,
2919  const unsigned int level)
2920  {
2921  std::vector<IndexSet> locally_owned_set;
2922  locally_owned_set.reserve(dofh.size());
2923  for (unsigned int j = 0; j < dofh.size(); ++j)
2924  if (level == numbers::invalid_unsigned_int)
2925  locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
2926  else
2927  locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
2928  return locally_owned_set;
2929  }
2930  } // namespace MatrixFreeImplementation
2931 } // namespace internal
2932 
2933 
2934 
2935 template <int dim, typename Number, typename VectorizedArrayType>
2936 template <typename QuadratureType, typename number2>
2937 void
2939  const DoFHandler<dim> & dof_handler,
2940  const AffineConstraints<number2> &constraints_in,
2941  const QuadratureType & quad,
2943  &additional_data)
2944 {
2945  std::vector<const DoFHandler<dim, dim> *> dof_handlers;
2946  std::vector<const AffineConstraints<number2> *> constraints;
2947  std::vector<QuadratureType> quads;
2948 
2949  dof_handlers.push_back(&dof_handler);
2950  constraints.push_back(&constraints_in);
2951  quads.push_back(quad);
2952 
2953  std::vector<IndexSet> locally_owned_sets =
2954  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2955  dof_handlers, additional_data.mg_level);
2956 
2957  std::vector<hp::QCollection<dim>> quad_hp;
2958  quad_hp.emplace_back(quad);
2959 
2960  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(
2962  dof_handlers,
2963  constraints,
2964  locally_owned_sets,
2965  quad_hp,
2966  additional_data);
2967 }
2968 
2969 
2970 
2971 template <int dim, typename Number, typename VectorizedArrayType>
2972 template <typename QuadratureType, typename number2, typename MappingType>
2973 void
2975  const MappingType & mapping,
2976  const DoFHandler<dim> & dof_handler,
2977  const AffineConstraints<number2> &constraints_in,
2978  const QuadratureType & quad,
2980  &additional_data)
2981 {
2982  std::vector<const DoFHandler<dim, dim> *> dof_handlers;
2983  std::vector<const AffineConstraints<number2> *> constraints;
2984 
2985  dof_handlers.push_back(&dof_handler);
2986  constraints.push_back(&constraints_in);
2987 
2988  std::vector<IndexSet> locally_owned_sets =
2989  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2990  dof_handlers, additional_data.mg_level);
2991 
2992  std::vector<hp::QCollection<dim>> quad_hp;
2993  quad_hp.emplace_back(quad);
2994 
2995  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
2996  dof_handlers,
2997  constraints,
2998  locally_owned_sets,
2999  quad_hp,
3000  additional_data);
3001 }
3002 
3003 
3004 
3005 template <int dim, typename Number, typename VectorizedArrayType>
3006 template <typename QuadratureType, typename number2>
3007 void
3009  const std::vector<const DoFHandler<dim> *> & dof_handler,
3010  const std::vector<const AffineConstraints<number2> *> &constraint,
3011  const std::vector<QuadratureType> & quad,
3013  &additional_data)
3014 {
3015  std::vector<IndexSet> locally_owned_set =
3016  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3017  dof_handler, additional_data.mg_level);
3018  std::vector<hp::QCollection<dim>> quad_hp;
3019  for (unsigned int q = 0; q < quad.size(); ++q)
3020  quad_hp.emplace_back(quad[q]);
3021 
3022  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(
3024  dof_handler,
3025  constraint,
3026  locally_owned_set,
3027  quad_hp,
3028  additional_data);
3029 }
3030 
3031 
3032 
3033 template <int dim, typename Number, typename VectorizedArrayType>
3034 template <typename QuadratureType,
3035  typename number2,
3036  typename DoFHandlerType,
3037  typename MappingType>
3038 void
3040  const MappingType & mapping,
3041  const std::vector<const DoFHandlerType *> & dof_handler,
3042  const std::vector<const AffineConstraints<number2> *> &constraint,
3043  const std::vector<QuadratureType> & quad,
3044  const AdditionalData & additional_data)
3045 {
3046  static_assert(dim == DoFHandlerType::dimension,
3047  "Dimension dim not equal to DoFHandlerType::dimension.");
3048 
3049  std::vector<const DoFHandler<dim> *> dof_handlers;
3050 
3051  for (const auto dh : dof_handler)
3052  dof_handlers.push_back(dh);
3053 
3054  this->reinit(mapping, dof_handlers, constraint, quad, additional_data);
3055 }
3056 
3057 
3058 
3059 template <int dim, typename Number, typename VectorizedArrayType>
3060 template <typename QuadratureType, typename number2>
3061 void
3063  const std::vector<const DoFHandler<dim> *> & dof_handler,
3064  const std::vector<const AffineConstraints<number2> *> &constraint,
3065  const QuadratureType & quad,
3067  &additional_data)
3068 {
3069  std::vector<IndexSet> locally_owned_set =
3070  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3071  dof_handler, additional_data.mg_level);
3072  std::vector<hp::QCollection<dim>> quad_hp;
3073  quad_hp.emplace_back(quad);
3074 
3075  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(
3077  dof_handler,
3078  constraint,
3079  locally_owned_set,
3080  quad_hp,
3081  additional_data);
3082 }
3083 
3084 
3085 
3086 template <int dim, typename Number, typename VectorizedArrayType>
3087 template <typename QuadratureType, typename number2, typename DoFHandlerType>
3088 void
3090  const std::vector<const DoFHandlerType *> & dof_handler,
3091  const std::vector<const AffineConstraints<number2> *> &constraint,
3092  const std::vector<QuadratureType> & quad,
3093  const AdditionalData & additional_data)
3094 {
3095  static_assert(dim == DoFHandlerType::dimension,
3096  "Dimension dim not equal to DoFHandlerType::dimension.");
3097 
3098  std::vector<const DoFHandler<dim> *> dof_handlers;
3099 
3100  for (const auto dh : dof_handler)
3101  dof_handlers.push_back(dof_handler);
3102 
3103  this->reinit(dof_handlers, constraint, quad, additional_data);
3104 }
3105 
3106 
3107 
3108 template <int dim, typename Number, typename VectorizedArrayType>
3109 template <typename QuadratureType, typename number2, typename MappingType>
3110 void
3112  const MappingType & mapping,
3113  const std::vector<const DoFHandler<dim> *> & dof_handler,
3114  const std::vector<const AffineConstraints<number2> *> &constraint,
3115  const QuadratureType & quad,
3117  &additional_data)
3118 {
3119  std::vector<IndexSet> locally_owned_set =
3120  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3121  dof_handler, additional_data.mg_level);
3122  std::vector<hp::QCollection<dim>> quad_hp;
3123  quad_hp.emplace_back(quad);
3124 
3125  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3126  dof_handler,
3127  constraint,
3128  locally_owned_set,
3129  quad_hp,
3130  additional_data);
3131 }
3132 
3133 
3134 
3135 template <int dim, typename Number, typename VectorizedArrayType>
3136 template <typename QuadratureType, typename number2, typename MappingType>
3137 void
3139  const MappingType & mapping,
3140  const std::vector<const DoFHandler<dim> *> & dof_handler,
3141  const std::vector<const AffineConstraints<number2> *> &constraint,
3142  const std::vector<QuadratureType> & quad,
3144  &additional_data)
3145 {
3146  std::vector<IndexSet> locally_owned_set =
3147  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3148  dof_handler, additional_data.mg_level);
3149  std::vector<hp::QCollection<dim>> quad_hp;
3150  for (unsigned int q = 0; q < quad.size(); ++q)
3151  quad_hp.emplace_back(quad[q]);
3152 
3153  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3154  dof_handler,
3155  constraint,
3156  locally_owned_set,
3157  quad_hp,
3158  additional_data);
3159 }
3160 
3161 
3162 
3163 template <int dim, typename Number, typename VectorizedArrayType>
3164 template <typename QuadratureType,
3165  typename number2,
3166  typename DoFHandlerType,
3167  typename MappingType>
3168 void
3170  const MappingType & mapping,
3171  const std::vector<const DoFHandlerType *> & dof_handler,
3172  const std::vector<const AffineConstraints<number2> *> &constraint,
3173  const QuadratureType & quad,
3174  const AdditionalData & additional_data)
3175 {
3176  static_assert(dim == DoFHandlerType::dimension,
3177  "Dimension dim not equal to DoFHandlerType::dimension.");
3178 
3179  std::vector<const DoFHandler<dim> *> dof_handlers;
3180 
3181  for (const auto dh : dof_handler)
3182  dof_handlers.push_back(dof_handler);
3183 
3184  this->reinit(mapping, dof_handlers, constraint, quad, additional_data);
3185 }
3186 
3187 
3188 
3189 template <int dim, typename Number, typename VectorizedArrayType>
3190 template <typename QuadratureType, typename number2, typename DoFHandlerType>
3191 void
3193  const std::vector<const DoFHandlerType *> & dof_handler,
3194  const std::vector<const AffineConstraints<number2> *> &constraint,
3195  const QuadratureType & quad,
3196  const AdditionalData & additional_data)
3197 {
3198  static_assert(dim == DoFHandlerType::dimension,
3199  "Dimension dim not equal to DoFHandlerType::dimension.");
3200 
3201  std::vector<const DoFHandler<dim> *> dof_handlers;
3202 
3203  for (const auto dh : dof_handler)
3204  dof_handlers.push_back(dof_handler);
3205 
3206  this->reinit(dof_handlers, constraint, quad, additional_data);
3207 }
3208 
3209 
3210 
3211 // ------------------------------ implementation of loops --------------------
3212 
3213 // internal helper functions that define how to call MPI data exchange
3214 // functions: for generic vectors, do nothing at all. For distributed vectors,
3215 // call update_ghost_values_start function and so on. If we have collections
3216 // of vectors, just do the individual functions of the components. In order to
3217 // keep ghost values consistent (whether we are in read or write mode), we
3218 // also reset the values at the end. the whole situation is a bit complicated
3219 // by the fact that we need to treat block vectors differently, which use some
3220 // additional helper functions to select the blocks and template magic.
3221 namespace internal
3222 {
3226  template <int dim, typename Number, typename VectorizedArrayType>
3227  struct VectorDataExchange
3228  {
3229  // A shift for the MPI messages to reduce the risk for accidental
3230  // interaction with other open communications that a user program might
3231  // set up (parallel vectors support unfinished communication). We let
3232  // the other vectors use the first 20 assigned numbers and start the
3233  // matrix-free communication.
3234  static constexpr unsigned int channel_shift = 20;
3235 
3236 
3237 
3242  VectorDataExchange(
3243  const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3244  const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3245  DataAccessOnFaces vector_face_access,
3246  const unsigned int n_components)
3247  : matrix_free(matrix_free)
3248  , vector_face_access(
3249  matrix_free.get_task_info().face_partition_data.empty() ?
3251  DataAccessOnFaces::unspecified :
3252  vector_face_access)
3253  , ghosts_were_set(false)
3254 # ifdef DEAL_II_WITH_MPI
3255  , tmp_data(n_components)
3256  , requests(n_components)
3257 # endif
3258  {
3259  (void)n_components;
3260  if (this->vector_face_access !=
3263  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3265  matrix_free.get_dof_info(c).vector_exchanger_face_variants.size(),
3266  5);
3267  }
3268 
3269 
3270 
3274  ~VectorDataExchange() // NOLINT
3275  {
3276 # ifdef DEAL_II_WITH_MPI
3277  for (unsigned int i = 0; i < tmp_data.size(); ++i)
3278  if (tmp_data[i] != nullptr)
3279  matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3280 # endif
3281  }
3282 
3283 
3284 
3289  template <typename VectorType>
3290  unsigned int
3291  find_vector_in_mf(const VectorType &vec,
3292  const bool check_global_compatibility = true) const
3293  {
3294  // case 1: vector was set up with MatrixFree::initialize_dof_vector()
3295  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3296  if (vec.get_partitioner().get() ==
3297  matrix_free.get_dof_info(c).vector_partitioner.get())
3298  return c;
3299 
3300  // case 2: user provided own partitioner (compatibility mode)
3301  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3302  if (check_global_compatibility ?
3303  vec.get_partitioner()->is_globally_compatible(
3304  *matrix_free.get_dof_info(c).vector_partitioner) :
3305  vec.get_partitioner()->is_compatible(
3306  *matrix_free.get_dof_info(c).vector_partitioner))
3307  return c;
3308 
3309  Assert(false,
3310  ExcNotImplemented("Could not find partitioner that fits vector"));
3311 
3313  }
3314 
3315 
3316 
3322  get_partitioner(const unsigned int mf_component) const
3323  {
3324  AssertDimension(matrix_free.get_dof_info(mf_component)
3325  .vector_exchanger_face_variants.size(),
3326  5);
3327  if (vector_face_access ==
3330  return *matrix_free.get_dof_info(mf_component)
3331  .vector_exchanger_face_variants[0];
3332  else if (vector_face_access ==
3335  return *matrix_free.get_dof_info(mf_component)
3336  .vector_exchanger_face_variants[1];
3337  else if (vector_face_access ==
3340  return *matrix_free.get_dof_info(mf_component)
3341  .vector_exchanger_face_variants[2];
3342  else if (vector_face_access ==
3345  return *matrix_free.get_dof_info(mf_component)
3346  .vector_exchanger_face_variants[3];
3347  else if (vector_face_access ==
3350  return *matrix_free.get_dof_info(mf_component)
3351  .vector_exchanger_face_variants[4];
3352  else
3353  return *matrix_free.get_dof_info(mf_component).vector_exchanger.get();
3354  }
3355 
3356 
3357 
3361  template <typename VectorType,
3362  typename std::enable_if<is_serial_or_dummy<VectorType>::value,
3363  VectorType>::type * = nullptr>
3364  void
3365  update_ghost_values_start(const unsigned int /*component_in_block_vector*/,
3366  const VectorType & /*vec*/)
3367  {}
3368 
3369 
3374  template <typename VectorType,
3375  typename std::enable_if<
3376  !has_update_ghost_values_start<VectorType>::value &&
3377  !is_serial_or_dummy<VectorType>::value,
3378  VectorType>::type * = nullptr>
3379  void
3380  update_ghost_values_start(const unsigned int component_in_block_vector,
3381  const VectorType & vec)
3382  {
3383  (void)component_in_block_vector;
3384  bool ghosts_set = vec.has_ghost_elements();
3385 
3386  Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3387  ghosts_set == false,
3388  ExcNotImplemented());
3389 
3390  if (ghosts_set)
3391  ghosts_were_set = true;
3392 
3393  vec.update_ghost_values();
3394  }
3395 
3396 
3397 
3403  template <typename VectorType,
3404  typename std::enable_if<
3405  has_update_ghost_values_start<VectorType>::value &&
3406  !has_exchange_on_subset<VectorType>::value,
3407  VectorType>::type * = nullptr>
3408  void
3409  update_ghost_values_start(const unsigned int component_in_block_vector,
3410  const VectorType & vec)
3411  {
3412  (void)component_in_block_vector;
3413  bool ghosts_set = vec.has_ghost_elements();
3414 
3415  Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3416  ghosts_set == false,
3417  ExcNotImplemented());
3418 
3419  if (ghosts_set)
3420  ghosts_were_set = true;
3421 
3422  vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3423  }
3424 
3425 
3426 
3433  template <typename VectorType,
3434  typename std::enable_if<
3435  has_update_ghost_values_start<VectorType>::value &&
3436  has_exchange_on_subset<VectorType>::value,
3437  VectorType>::type * = nullptr>
3438  void
3439  update_ghost_values_start(const unsigned int component_in_block_vector,
3440  const VectorType & vec)
3441  {
3442  static_assert(
3443  std::is_same<Number, typename VectorType::value_type>::value,
3444  "Type mismatch between VectorType and VectorDataExchange");
3445  (void)component_in_block_vector;
3446  bool ghosts_set = vec.has_ghost_elements();
3447 
3448  Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3449  ghosts_set == false,
3450  ExcNotImplemented());
3451 
3452  if (ghosts_set)
3453  ghosts_were_set = true;
3454 
3455  if (vec.size() != 0)
3456  {
3457 # ifdef DEAL_II_WITH_MPI
3458  const unsigned int mf_component = find_vector_in_mf(vec);
3459 
3460  const auto &part = get_partitioner(mf_component);
3461 
3462  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3463  part.n_import_sm_procs() == 0)
3464  return;
3465 
3466  tmp_data[component_in_block_vector] =
3467  matrix_free.acquire_scratch_data_non_threadsafe();
3468  tmp_data[component_in_block_vector]->resize_fast(
3469  part.n_import_indices());
3470  AssertDimension(requests.size(), tmp_data.size());
3471 
3472  part.export_to_ghosted_array_start(
3473  component_in_block_vector * 2 + channel_shift,
3474  ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3475  vec.shared_vector_data(),
3476  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3477  part.locally_owned_size(),
3478  matrix_free.get_dof_info(mf_component)
3479  .vector_partitioner->n_ghost_indices()),
3480  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3481  part.n_import_indices()),
3482  this->requests[component_in_block_vector]);
3483 # endif
3484  }
3485  }
3486 
3487 
3488 
3493  template <
3494  typename VectorType,
3495  typename std::enable_if<!has_update_ghost_values_start<VectorType>::value,
3496  VectorType>::type * = nullptr>
3497  void
3498  update_ghost_values_finish(const unsigned int /*component_in_block_vector*/,
3499  const VectorType & /*vec*/)
3500  {}
3501 
3502 
3503 
3509  template <typename VectorType,
3510  typename std::enable_if<
3511  has_update_ghost_values_start<VectorType>::value &&
3512  !has_exchange_on_subset<VectorType>::value,
3513  VectorType>::type * = nullptr>
3514  void
3515  update_ghost_values_finish(const unsigned int component_in_block_vector,
3516  const VectorType & vec)
3517  {
3518  (void)component_in_block_vector;
3519  vec.update_ghost_values_finish();
3520  }
3521 
3522 
3523 
3530  template <typename VectorType,
3531  typename std::enable_if<
3532  has_update_ghost_values_start<VectorType>::value &&
3533  has_exchange_on_subset<VectorType>::value,
3534  VectorType>::type * = nullptr>
3535  void
3536  update_ghost_values_finish(const unsigned int component_in_block_vector,
3537  const VectorType & vec)
3538  {
3539  static_assert(
3540  std::is_same<Number, typename VectorType::value_type>::value,
3541  "Type mismatch between VectorType and VectorDataExchange");
3542  (void)component_in_block_vector;
3543 
3544  if (vec.size() != 0)
3545  {
3546 # ifdef DEAL_II_WITH_MPI
3547  AssertIndexRange(component_in_block_vector, tmp_data.size());
3548  AssertDimension(requests.size(), tmp_data.size());
3549 
3550  const unsigned int mf_component = find_vector_in_mf(vec);
3551 
3552  const auto &part = get_partitioner(mf_component);
3553 
3554  if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3555  part.n_import_sm_procs() != 0)
3556  {
3557  part.export_to_ghosted_array_finish(
3558  ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3559  vec.shared_vector_data(),
3560  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3561  part.locally_owned_size(),
3562  matrix_free.get_dof_info(mf_component)
3563  .vector_partitioner->n_ghost_indices()),
3564  this->requests[component_in_block_vector]);
3565 
3566  matrix_free.release_scratch_data_non_threadsafe(
3567  tmp_data[component_in_block_vector]);
3568  tmp_data[component_in_block_vector] = nullptr;
3569  }
3570 # endif
3571  }
3572  // let vector know that ghosts are being updated and we can read from
3573  // them
3574  vec.set_ghost_state(true);
3575  }
3576 
3577 
3578 
3582  template <typename VectorType,
3583  typename std::enable_if<is_serial_or_dummy<VectorType>::value,
3584  VectorType>::type * = nullptr>
3585  void
3586  compress_start(const unsigned int /*component_in_block_vector*/,
3587  VectorType & /*vec*/)
3588  {}
3589 
3590 
3591 
3596  template <typename VectorType,
3597  typename std::enable_if<!has_compress_start<VectorType>::value &&
3598  !is_serial_or_dummy<VectorType>::value,
3599  VectorType>::type * = nullptr>
3600  void
3601  compress_start(const unsigned int component_in_block_vector,
3602  VectorType & vec)
3603  {
3604  (void)component_in_block_vector;
3605  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3606  vec.compress(::VectorOperation::add);
3607  }
3608 
3609 
3610 
3616  template <
3617  typename VectorType,
3618  typename std::enable_if<has_compress_start<VectorType>::value &&
3619  !has_exchange_on_subset<VectorType>::value,
3620  VectorType>::type * = nullptr>
3621  void
3622  compress_start(const unsigned int component_in_block_vector,
3623  VectorType & vec)
3624  {
3625  (void)component_in_block_vector;
3626  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3627  vec.compress_start(component_in_block_vector + channel_shift);
3628  }
3629 
3630 
3631 
3638  template <
3639  typename VectorType,
3640  typename std::enable_if<has_compress_start<VectorType>::value &&
3641  has_exchange_on_subset<VectorType>::value,
3642  VectorType>::type * = nullptr>
3643  void
3644  compress_start(const unsigned int component_in_block_vector,
3645  VectorType & vec)
3646  {
3647  static_assert(
3648  std::is_same<Number, typename VectorType::value_type>::value,
3649  "Type mismatch between VectorType and VectorDataExchange");
3650  (void)component_in_block_vector;
3651  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3652 
3653  if (vec.size() != 0)
3654  {
3655 # ifdef DEAL_II_WITH_MPI
3656  const unsigned int mf_component = find_vector_in_mf(vec);
3657 
3658  const auto &part = get_partitioner(mf_component);
3659 
3660  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3661  part.n_import_sm_procs() == 0)
3662  return;
3663 
3664  tmp_data[component_in_block_vector] =
3665  matrix_free.acquire_scratch_data_non_threadsafe();
3666  tmp_data[component_in_block_vector]->resize_fast(
3667  part.n_import_indices());
3668  AssertDimension(requests.size(), tmp_data.size());
3669 
3670  part.import_from_ghosted_array_start(
3672  component_in_block_vector * 2 + channel_shift,
3673  ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3674  vec.shared_vector_data(),
3675  ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3676  matrix_free.get_dof_info(mf_component)
3677  .vector_partitioner->n_ghost_indices()),
3678  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3679  part.n_import_indices()),
3680  this->requests[component_in_block_vector]);
3681 # endif
3682  }
3683  }
3684 
3685 
3686 
3691  template <typename VectorType,
3692  typename std::enable_if<!has_compress_start<VectorType>::value,
3693  VectorType>::type * = nullptr>
3694  void
3695  compress_finish(const unsigned int /*component_in_block_vector*/,
3696  VectorType & /*vec*/)
3697  {}
3698 
3699 
3700 
3706  template <
3707  typename VectorType,
3708  typename std::enable_if<has_compress_start<VectorType>::value &&
3709  !has_exchange_on_subset<VectorType>::value,
3710  VectorType>::type * = nullptr>
3711  void
3712  compress_finish(const unsigned int component_in_block_vector,
3713  VectorType & vec)
3714  {
3715  (void)component_in_block_vector;
3716  vec.compress_finish(::VectorOperation::add);
3717  }
3718 
3719 
3720 
3727  template <
3728  typename VectorType,
3729  typename std::enable_if<has_compress_start<VectorType>::value &&
3730  has_exchange_on_subset<VectorType>::value,
3731  VectorType>::type * = nullptr>
3732  void
3733  compress_finish(const unsigned int component_in_block_vector,
3734  VectorType & vec)
3735  {
3736  static_assert(
3737  std::is_same<Number, typename VectorType::value_type>::value,
3738  "Type mismatch between VectorType and VectorDataExchange");
3739  (void)component_in_block_vector;
3740  if (vec.size() != 0)
3741  {
3742 # ifdef DEAL_II_WITH_MPI
3743  AssertIndexRange(component_in_block_vector, tmp_data.size());
3744  AssertDimension(requests.size(), tmp_data.size());
3745 
3746  const unsigned int mf_component = find_vector_in_mf(vec);
3747 
3748  const auto &part = get_partitioner(mf_component);
3749 
3750  if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3751  part.n_import_sm_procs() != 0)
3752  {
3753  part.import_from_ghosted_array_finish(
3755  ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3756  vec.shared_vector_data(),
3757  ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3758  matrix_free.get_dof_info(mf_component)
3759  .vector_partitioner->n_ghost_indices()),
3761  tmp_data[component_in_block_vector]->begin(),
3762  part.n_import_indices()),
3763  this->requests[component_in_block_vector]);
3764 
3765  matrix_free.release_scratch_data_non_threadsafe(
3766  tmp_data[component_in_block_vector]);
3767  tmp_data[component_in_block_vector] = nullptr;
3768  }
3769 
3771  {
3772  const int ierr =
3773  MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3774  AssertThrowMPI(ierr);
3775  }
3776 # endif
3777  }
3778  }
3779 
3780 
3781 
3785  template <typename VectorType,
3786  typename std::enable_if<is_serial_or_dummy<VectorType>::value,
3787  VectorType>::type * = nullptr>
3788  void
3789  reset_ghost_values(const VectorType & /*vec*/) const
3790  {}
3791 
3792 
3793 
3798  template <
3799  typename VectorType,
3800  typename std::enable_if<!has_exchange_on_subset<VectorType>::value &&
3801  !is_serial_or_dummy<VectorType>::value,
3802  VectorType>::type * = nullptr>
3803  void
3804  reset_ghost_values(const VectorType &vec) const
3805  {
3806  if (ghosts_were_set == true)
3807  return;
3808 
3809  vec.zero_out_ghost_values();
3810  }
3811 
3812 
3813 
3819  template <typename VectorType,
3820  typename std::enable_if<has_exchange_on_subset<VectorType>::value,
3821  VectorType>::type * = nullptr>
3822  void
3823  reset_ghost_values(const VectorType &vec) const
3824  {
3825  static_assert(
3826  std::is_same<Number, typename VectorType::value_type>::value,
3827  "Type mismatch between VectorType and VectorDataExchange");
3828  if (ghosts_were_set == true)
3829  return;
3830 
3831  if (vec.size() != 0)
3832  {
3833 # ifdef DEAL_II_WITH_MPI
3834  AssertDimension(requests.size(), tmp_data.size());
3835 
3836  const unsigned int mf_component = find_vector_in_mf(vec);
3837 
3838  const auto &part = get_partitioner(mf_component);
3839 
3840  if (part.n_ghost_indices() > 0)
3841  {
3842  part.reset_ghost_values(ArrayView<Number>(
3844  .begin() +
3845  part.locally_owned_size(),
3846  matrix_free.get_dof_info(mf_component)
3847  .vector_partitioner->n_ghost_indices()));
3848  }
3849 
3850 # endif
3851  }
3852  // let vector know that it's not ghosted anymore
3853  vec.set_ghost_state(false);
3854  }
3855 
3856 
3857 
3863  template <typename VectorType,
3864  typename std::enable_if<has_exchange_on_subset<VectorType>::value,
3865  VectorType>::type * = nullptr>
3866  void
3867  zero_vector_region(const unsigned int range_index, VectorType &vec) const
3868  {
3869  static_assert(
3870  std::is_same<Number, typename VectorType::value_type>::value,
3871  "Type mismatch between VectorType and VectorDataExchange");
3872  if (range_index == numbers::invalid_unsigned_int)
3873  vec = Number();
3874  else
3875  {
3876  const unsigned int mf_component = find_vector_in_mf(vec, false);
3878  matrix_free.get_dof_info(mf_component);
3879  Assert(dof_info.vector_zero_range_list_index.empty() == false,
3880  ExcNotInitialized());
3881 
3882  Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3883  ExcInternalError());
3884  AssertIndexRange(range_index,
3885  dof_info.vector_zero_range_list_index.size() - 1);
3886  for (unsigned int id =
3887  dof_info.vector_zero_range_list_index[range_index];
3888  id != dof_info.vector_zero_range_list_index[range_index + 1];
3889  ++id)
3890  std::memset(vec.begin() + dof_info.vector_zero_range_list[id].first,
3891  0,
3892  (dof_info.vector_zero_range_list[id].second -
3893  dof_info.vector_zero_range_list[id].first) *
3894  sizeof(Number));
3895  }
3896  }
3897 
3898 
3899 
3905  template <
3906  typename VectorType,
3907  typename std::enable_if<!has_exchange_on_subset<VectorType>::value,
3908  VectorType>::type * = nullptr,
3909  typename VectorType::value_type * = nullptr>
3910  void
3911  zero_vector_region(const unsigned int range_index, VectorType &vec) const
3912  {
3913  if (range_index == numbers::invalid_unsigned_int || range_index == 0)
3914  vec = typename VectorType::value_type();
3915  }
3916 
3917 
3918 
3923  void
3924  zero_vector_region(const unsigned int, ...) const
3925  {
3926  Assert(false,
3927  ExcNotImplemented("Zeroing is only implemented for vector types "
3928  "which provide VectorType::value_type"));
3929  }
3930 
3931 
3932 
3933  const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
3934  const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3935  DataAccessOnFaces vector_face_access;
3936  bool ghosts_were_set;
3937 # ifdef DEAL_II_WITH_MPI
3938  std::vector<AlignedVector<Number> *> tmp_data;
3939  std::vector<std::vector<MPI_Request>> requests;
3940 # endif
3941  }; // VectorDataExchange
3942 
3943  template <typename VectorStruct>
3944  unsigned int
3945  n_components(const VectorStruct &vec);
3946 
3947  template <typename VectorStruct>
3948  unsigned int
3949  n_components_block(const VectorStruct &vec,
3950  std::integral_constant<bool, true>)
3951  {
3952  unsigned int components = 0;
3953  for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
3954  components += n_components(vec.block(bl));
3955  return components;
3956  }
3957 
3958  template <typename VectorStruct>
3959  unsigned int
3960  n_components_block(const VectorStruct &, std::integral_constant<bool, false>)
3961  {
3962  return 1;
3963  }
3964 
3965  template <typename VectorStruct>
3966  unsigned int
3967  n_components(const VectorStruct &vec)
3968  {
3969  return n_components_block(
3970  vec, std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3971  }
3972 
3973  template <typename VectorStruct>
3974  inline unsigned int
3975  n_components(const std::vector<VectorStruct> &vec)
3976  {
3977  unsigned int components = 0;
3978  for (unsigned int comp = 0; comp < vec.size(); ++comp)
3979  components += n_components_block(
3980  vec[comp],
3981  std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3982  return components;
3983  }
3984 
3985  template <typename VectorStruct>
3986  inline unsigned int
3987  n_components(const std::vector<VectorStruct *> &vec)
3988  {
3989  unsigned int components = 0;
3990  for (unsigned int comp = 0; comp < vec.size(); ++comp)
3991  components += n_components_block(
3992  *vec[comp],
3993  std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3994  return components;
3995  }
3996 
3997 
3998 
3999  // A helper function to identify block vectors with many components where we
4000  // should not try to overlap computations and communication because there
4001  // would be too many outstanding communication requests.
4002 
4003  // default value for vectors that do not have communication_block_size
4004  template <
4005  typename VectorStruct,
4006  typename std::enable_if<!has_communication_block_size<VectorStruct>::value,
4007  VectorStruct>::type * = nullptr>
4008  constexpr unsigned int
4009  get_communication_block_size(const VectorStruct &)
4010  {
4012  }
4013 
4014 
4015 
4016  template <
4017  typename VectorStruct,
4018  typename std::enable_if<has_communication_block_size<VectorStruct>::value,
4019  VectorStruct>::type * = nullptr>
4020  constexpr unsigned int
4021  get_communication_block_size(const VectorStruct &)
4022  {
4023  return VectorStruct::communication_block_size;
4024  }
4025 
4026 
4027 
4028  // --------------------------------------------------------------------------
4029  // below we have wrappers to distinguish between block and non-block vectors.
4030  // --------------------------------------------------------------------------
4031 
4032  //
4033  // update_ghost_values_start
4034  //
4035 
4036  // update_ghost_values for block vectors
4037  template <int dim,
4038  typename VectorStruct,
4039  typename Number,
4040  typename VectorizedArrayType,
4041  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4042  VectorStruct>::type * = nullptr>
4043  void
4044  update_ghost_values_start(
4045  const VectorStruct & vec,
4046  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4047  const unsigned int channel = 0)
4048  {
4049  if (get_communication_block_size(vec) < vec.n_blocks())
4050  {
4051  // don't forget to set ghosts_were_set, that otherwise happens
4052  // inside VectorDataExchange::update_ghost_values_start()
4053  exchanger.ghosts_were_set = vec.has_ghost_elements();
4054  vec.update_ghost_values();
4055  }
4056  else
4057  {
4058  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4059  update_ghost_values_start(vec.block(i), exchanger, channel + i);
4060  }
4061  }
4062 
4063 
4064 
4065  // update_ghost_values for non-block vectors
4066  template <int dim,
4067  typename VectorStruct,
4068  typename Number,
4069  typename VectorizedArrayType,
4070  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4071  VectorStruct>::type * = nullptr>
4072  void
4073  update_ghost_values_start(
4074  const VectorStruct & vec,
4075  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4076  const unsigned int channel = 0)
4077  {
4078  exchanger.update_ghost_values_start(channel, vec);
4079  }
4080 
4081 
4082 
4083  // update_ghost_values_start() for vector of vectors
4084  template <int dim,
4085  typename VectorStruct,
4086  typename Number,
4087  typename VectorizedArrayType>
4088  inline void
4089  update_ghost_values_start(
4090  const std::vector<VectorStruct> & vec,
4091  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4092  {
4093  unsigned int component_index = 0;
4094  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4095  {
4096  update_ghost_values_start(vec[comp], exchanger, component_index);
4097  component_index += n_components(vec[comp]);
4098  }
4099  }
4100 
4101 
4102 
4103  // update_ghost_values_start() for vector of pointers to vectors
4104  template <int dim,
4105  typename VectorStruct,
4106  typename Number,
4107  typename VectorizedArrayType>
4108  inline void
4109  update_ghost_values_start(
4110  const std::vector<VectorStruct *> & vec,
4111  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4112  {
4113  unsigned int component_index = 0;
4114  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4115  {
4116  update_ghost_values_start(*vec[comp], exchanger, component_index);
4117  component_index += n_components(*vec[comp]);
4118  }
4119  }
4120 
4121 
4122 
4123  //
4124  // update_ghost_values_finish
4125  //
4126 
4127  // for block vectors
4128  template <int dim,
4129  typename VectorStruct,
4130  typename Number,
4131  typename VectorizedArrayType,
4132  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4133  VectorStruct>::type * = nullptr>
4134  void
4135  update_ghost_values_finish(
4136  const VectorStruct & vec,
4137  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4138  const unsigned int channel = 0)
4139  {
4140  if (get_communication_block_size(vec) < vec.n_blocks())
4141  {
4142  // do nothing, everything has already been completed in the _start()
4143  // call
4144  }
4145  else
4146  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4147  update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4148  }
4149 
4150 
4151 
4152  // for non-block vectors
4153  template <int dim,
4154  typename VectorStruct,
4155  typename Number,
4156  typename VectorizedArrayType,
4157  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4158  VectorStruct>::type * = nullptr>
4159  void
4160  update_ghost_values_finish(
4161  const VectorStruct & vec,
4162  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4163  const unsigned int channel = 0)
4164  {
4165  exchanger.update_ghost_values_finish(channel, vec);
4166  }
4167 
4168 
4169 
4170  // for vector of vectors
4171  template <int dim,
4172  typename VectorStruct,
4173  typename Number,
4174  typename VectorizedArrayType>
4175  inline void
4176  update_ghost_values_finish(
4177  const std::vector<VectorStruct> & vec,
4178  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4179  {
4180  unsigned int component_index = 0;
4181  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4182  {
4183  update_ghost_values_finish(vec[comp], exchanger, component_index);
4184  component_index += n_components(vec[comp]);
4185  }
4186  }
4187 
4188 
4189 
4190  // for vector of pointers to vectors
4191  template <int dim,
4192  typename VectorStruct,
4193  typename Number,
4194  typename VectorizedArrayType>
4195  inline void
4196  update_ghost_values_finish(
4197  const std::vector<VectorStruct *> & vec,
4198  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4199  {
4200  unsigned int component_index = 0;
4201  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4202  {
4203  update_ghost_values_finish(*vec[comp], exchanger, component_index);
4204  component_index += n_components(*vec[comp]);
4205  }
4206  }
4207 
4208 
4209 
4210  //
4211  // compress_start
4212  //
4213 
4214  // for block vectors
4215  template <int dim,
4216  typename VectorStruct,
4217  typename Number,
4218  typename VectorizedArrayType,
4219  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4220  VectorStruct>::type * = nullptr>
4221  inline void
4222  compress_start(
4223  VectorStruct & vec,
4224  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4225  const unsigned int channel = 0)
4226  {
4227  if (get_communication_block_size(vec) < vec.n_blocks())
4228  vec.compress(::VectorOperation::add);
4229  else
4230  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4231  compress_start(vec.block(i), exchanger, channel + i);
4232  }
4233 
4234 
4235 
4236  // for non-block vectors
4237  template <int dim,
4238  typename VectorStruct,
4239  typename Number,
4240  typename VectorizedArrayType,
4241  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4242  VectorStruct>::type * = nullptr>
4243  inline void
4244  compress_start(
4245  VectorStruct & vec,
4246  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4247  const unsigned int channel = 0)
4248  {
4249  exchanger.compress_start(channel, vec);
4250  }
4251 
4252 
4253 
4254  // for std::vector of vectors
4255  template <int dim,
4256  typename VectorStruct,
4257  typename Number,
4258  typename VectorizedArrayType>
4259  inline void
4260  compress_start(
4261  std::vector<VectorStruct> & vec,
4262  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4263  {
4264  unsigned int component_index = 0;
4265  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4266  {
4267  compress_start(vec[comp], exchanger, component_index);
4268  component_index += n_components(vec[comp]);
4269  }
4270  }
4271 
4272 
4273 
4274  // for std::vector of pointer to vectors
4275  template <int dim,
4276  typename VectorStruct,
4277  typename Number,
4278  typename VectorizedArrayType>
4279  inline void
4280  compress_start(
4281  std::vector<VectorStruct *> & vec,
4282  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4283  {
4284  unsigned int component_index = 0;
4285  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4286  {
4287  compress_start(*vec[comp], exchanger, component_index);
4288  component_index += n_components(*vec[comp]);
4289  }
4290  }
4291 
4292 
4293 
4294  //
4295  // compress_finish
4296  //
4297 
4298  // for block vectors
4299  template <int dim,
4300  typename VectorStruct,
4301  typename Number,
4302  typename VectorizedArrayType,
4303  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4304  VectorStruct>::type * = nullptr>
4305  inline void
4306  compress_finish(
4307  VectorStruct & vec,
4308  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4309  const unsigned int channel = 0)
4310  {
4311  if (get_communication_block_size(vec) < vec.n_blocks())
4312  {
4313  // do nothing, everything has already been completed in the _start()
4314  // call
4315  }
4316  else
4317  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4318  compress_finish(vec.block(i), exchanger, channel + i);
4319  }
4320 
4321 
4322 
4323  // for non-block vectors
4324  template <int dim,
4325  typename VectorStruct,
4326  typename Number,
4327  typename VectorizedArrayType,
4328  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4329  VectorStruct>::type * = nullptr>
4330  inline void
4331  compress_finish(
4332  VectorStruct & vec,
4333  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4334  const unsigned int channel = 0)
4335  {
4336  exchanger.compress_finish(channel, vec);
4337  }
4338 
4339 
4340 
4341  // for std::vector of vectors
4342  template <int dim,
4343  typename VectorStruct,
4344  typename Number,
4345  typename VectorizedArrayType>
4346  inline void
4347  compress_finish(
4348  std::vector<VectorStruct> & vec,
4349  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4350  {
4351  unsigned int component_index = 0;
4352  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4353  {
4354  compress_finish(vec[comp], exchanger, component_index);
4355  component_index += n_components(vec[comp]);
4356  }
4357  }
4358 
4359 
4360 
4361  // for std::vector of pointer to vectors
4362  template <int dim,
4363  typename VectorStruct,
4364  typename Number,
4365  typename VectorizedArrayType>
4366  inline void
4367  compress_finish(
4368  std::vector<VectorStruct *> & vec,
4369  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4370  {
4371  unsigned int component_index = 0;
4372  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4373  {
4374  compress_finish(*vec[comp], exchanger, component_index);
4375  component_index += n_components(*vec[comp]);
4376  }
4377  }
4378 
4379 
4380 
4381  //
4382  // reset_ghost_values:
4383  //
4384  // if the input vector did not have ghosts imported, clear them here again
4385  // in order to avoid subsequent operations e.g. in linear solvers to work
4386  // with ghosts all the time
4387  //
4388 
4389  // for block vectors
4390  template <int dim,
4391  typename VectorStruct,
4392  typename Number,
4393  typename VectorizedArrayType,
4394  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4395  VectorStruct>::type * = nullptr>
4396  inline void
4397  reset_ghost_values(
4398  const VectorStruct & vec,
4399  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4400  {
4401  // return immediately if there is nothing to do.
4402  if (exchanger.ghosts_were_set == true)
4403  return;
4404 
4405  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4406  reset_ghost_values(vec.block(i), exchanger);
4407  }
4408 
4409 
4410 
4411  // for non-block vectors
4412  template <int dim,
4413  typename VectorStruct,
4414  typename Number,
4415  typename VectorizedArrayType,
4416  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4417  VectorStruct>::type * = nullptr>
4418  inline void
4419  reset_ghost_values(
4420  const VectorStruct & vec,
4421  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4422  {
4423  exchanger.reset_ghost_values(vec);
4424  }
4425 
4426 
4427 
4428  // for std::vector of vectors
4429  template <int dim,
4430  typename VectorStruct,
4431  typename Number,
4432  typename VectorizedArrayType>
4433  inline void
4434  reset_ghost_values(
4435  const std::vector<VectorStruct> & vec,
4436  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4437  {
4438  // return immediately if there is nothing to do.
4439  if (exchanger.ghosts_were_set == true)
4440  return;
4441 
4442  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4443  reset_ghost_values(vec[comp], exchanger);
4444  }
4445 
4446 
4447 
4448  // for std::vector of pointer to vectors
4449  template <int dim,
4450  typename VectorStruct,
4451  typename Number,
4452  typename VectorizedArrayType>
4453  inline void
4454  reset_ghost_values(
4455  const std::vector<VectorStruct *> & vec,
4456  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4457  {
4458  // return immediately if there is nothing to do.
4459  if (exchanger.ghosts_were_set == true)
4460  return;
4461 
4462  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4463  reset_ghost_values(*vec[comp], exchanger);
4464  }
4465 
4466 
4467 
4468  //
4469  // zero_vector_region
4470  //
4471 
4472  // for block vectors
4473  template <int dim,
4474  typename VectorStruct,
4475  typename Number,
4476  typename VectorizedArrayType,
4477  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4478  VectorStruct>::type * = nullptr>
4479  inline void
4480  zero_vector_region(
4481  const unsigned int range_index,
4482  VectorStruct & vec,
4483  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4484  {
4485  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4486  exchanger.zero_vector_region(range_index, vec.block(i));
4487  }
4488 
4489 
4490 
4491  // for non-block vectors
4492  template <int dim,
4493  typename VectorStruct,
4494  typename Number,
4495  typename VectorizedArrayType,
4496  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4497  VectorStruct>::type * = nullptr>
4498  inline void
4499  zero_vector_region(
4500  const unsigned int range_index,
4501  VectorStruct & vec,
4502  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4503  {
4504  exchanger.zero_vector_region(range_index, vec);
4505  }
4506 
4507 
4508 
4509  // for std::vector of vectors
4510  template <int dim,
4511  typename VectorStruct,
4512  typename Number,
4513  typename VectorizedArrayType>
4514  inline void
4515  zero_vector_region(
4516  const unsigned int range_index,
4517  std::vector<VectorStruct> & vec,
4518  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4519  {
4520  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4521  zero_vector_region(range_index, vec[comp], exchanger);
4522  }
4523 
4524 
4525 
4526  // for std::vector of pointers to vectors
4527  template <int dim,
4528  typename VectorStruct,
4529  typename Number,
4530  typename VectorizedArrayType>
4531  inline void
4532  zero_vector_region(
4533  const unsigned int range_index,
4534  std::vector<VectorStruct *> & vec,
4535  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4536  {
4537  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4538  zero_vector_region(range_index, *vec[comp], exchanger);
4539  }
4540 
4541 
4542 
4543  namespace MatrixFreeFunctions
4544  {
4545  // struct to select between a const interface and a non-const interface
4546  // for MFWorker
4547  template <typename, typename, typename, typename, bool>
4548  struct InterfaceSelector
4549  {};
4550 
4551  // Version for constant functions
4552  template <typename MF,
4553  typename InVector,
4554  typename OutVector,
4555  typename Container>
4556  struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4557  {
4558  using function_type = void (Container::*)(
4559  const MF &,
4560  OutVector &,
4561  const InVector &,
4562  const std::pair<unsigned int, unsigned int> &) const;
4563  };
4564 
4565  // Version for non-constant functions
4566  template <typename MF,
4567  typename InVector,
4568  typename OutVector,
4569  typename Container>
4570  struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4571  {
4572  using function_type =
4573  void (Container::*)(const MF &,
4574  OutVector &,
4575  const InVector &,
4576  const std::pair<unsigned int, unsigned int> &);
4577  };
4578  } // namespace MatrixFreeFunctions
4579 
4580 
4581 
4582  // A implementation class for the worker object that runs the various
4583  // operations we want to perform during the matrix-free loop
4584  template <typename MF,
4585  typename InVector,
4586  typename OutVector,
4587  typename Container,
4588  bool is_constant>
4589  class MFWorker : public MFWorkerInterface
4590  {
4591  public:
4592  // An alias to make the arguments further down more readable
4593  using function_type = typename MatrixFreeFunctions::
4594  InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4595  function_type;
4596 
4597  // constructor, binds all the arguments to this class
4598  MFWorker(const MF & matrix_free,
4599  const InVector & src,
4600  OutVector & dst,
4601  const bool zero_dst_vector_setting,
4602  const Container & container,
4603  function_type cell_function,
4604  function_type face_function,
4605  function_type boundary_function,
4606  const typename MF::DataAccessOnFaces src_vector_face_access =
4608  const typename MF::DataAccessOnFaces dst_vector_face_access =
4610  const std::function<void(const unsigned int, const unsigned int)>
4611  &operation_before_loop = {},
4612  const std::function<void(const unsigned int, const unsigned int)>
4613  & operation_after_loop = {},
4614  const unsigned int dof_handler_index_pre_post = 0)
4615  : matrix_free(matrix_free)
4616  , container(const_cast<Container &>(container))
4617  , cell_function(cell_function)
4618  , face_function(face_function)
4619  , boundary_function(boundary_function)
4620  , src(src)
4621  , dst(dst)
4622  , src_data_exchanger(matrix_free,
4623  src_vector_face_access,
4624  n_components(src))
4625  , dst_data_exchanger(matrix_free,
4626  dst_vector_face_access,
4627  n_components(dst))
4628  , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4629  , zero_dst_vector_setting(zero_dst_vector_setting &&
4630  !src_and_dst_are_same)
4631  , operation_before_loop(operation_before_loop)
4632  , operation_after_loop(operation_after_loop)
4633  , dof_handler_index_pre_post(dof_handler_index_pre_post)
4634  {}
4635 
4636  // Runs the cell work. If no function is given, nothing is done
4637  virtual void
4638  cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4639  {
4640  if (cell_function != nullptr && cell_range.second > cell_range.first)
4641  for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4642  {
4643  const auto cell_subrange =
4644  matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4645 
4646  if (cell_subrange.second <= cell_subrange.first)
4647  continue;
4648 
4649  (container.*
4650  cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4651  }
4652  }
4653 
4654  virtual void
4655  cell(const unsigned int range_index) override
4656  {
4657  process_range(cell_function,
4658  matrix_free.get_task_info().cell_partition_data_hp_ptr,
4659  matrix_free.get_task_info().cell_partition_data_hp,
4660  range_index);
4661  }
4662 
4663  virtual void
4664  face(const unsigned int range_index) override
4665  {
4666  process_range(face_function,
4667  matrix_free.get_task_info().face_partition_data_hp_ptr,
4668  matrix_free.get_task_info().face_partition_data_hp,
4669  range_index);
4670  }
4671 
4672  virtual void
4673  boundary(const unsigned int range_index) override
4674  {
4675  process_range(boundary_function,
4676  matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4677  matrix_free.get_task_info().boundary_partition_data_hp,
4678  range_index);
4679  }
4680 
4681  private:
4682  void
4683  process_range(const function_type & fu,
4684  const std::vector<unsigned int> &ptr,
4685  const std::vector<unsigned int> &data,
4686  const unsigned int range_index)
4687  {
4688  if (fu == nullptr)
4689  return;
4690 
4691  for (unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4692  (container.*fu)(matrix_free,
4693  this->dst,
4694  this->src,
4695  std::make_pair(data[2 * i], data[2 * i + 1]));
4696  }
4697 
4698  public:
4699  // Starts the communication for the update ghost values operation. We
4700  // cannot call this update if ghost and destination are the same because
4701  // that would introduce spurious entries in the destination (there is also
4702  // the problem that reading from a vector that we also write to is usually
4703  // not intended in case there is overlap, but this is up to the
4704  // application code to decide and we cannot catch this case here).
4705  virtual void
4706  vector_update_ghosts_start() override
4707  {
4708  if (!src_and_dst_are_same)
4709  internal::update_ghost_values_start(src, src_data_exchanger);
4710  }
4711 
4712  // Finishes the communication for the update ghost values operation
4713  virtual void
4714  vector_update_ghosts_finish() override
4715  {
4716  if (!src_and_dst_are_same)
4717  internal::update_ghost_values_finish(src, src_data_exchanger);
4718  }
4719 
4720  // Starts the communication for the vector compress operation
4721  virtual void
4722  vector_compress_start() override
4723  {
4724  internal::compress_start(dst, dst_data_exchanger);
4725  }
4726 
4727  // Finishes the communication for the vector compress operation
4728  virtual void
4729  vector_compress_finish() override
4730  {
4731  internal::compress_finish(dst, dst_data_exchanger);
4732  if (!src_and_dst_are_same)
4733  internal::reset_ghost_values(src, src_data_exchanger);
4734  }
4735 
4736  // Zeros the given input vector
4737  virtual void
4738  zero_dst_vector_range(const unsigned int range_index) override
4739  {
4740  if (zero_dst_vector_setting)
4741  internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4742  }
4743 
4744  virtual void
4745  cell_loop_pre_range(const unsigned int range_index) override
4746  {
4747  if (operation_before_loop)
4748  {
4750  matrix_free.get_dof_info(dof_handler_index_pre_post);
4751  if (range_index == numbers::invalid_unsigned_int)
4752  {
4753  // Case with threaded loop -> currently no overlap implemented
4755  0U,
4756  dof_info.vector_partitioner->locally_owned_size(),
4757  operation_before_loop,
4759  }
4760  else
4761  {
4762  AssertIndexRange(range_index,
4763  dof_info.cell_loop_pre_list_index.size() - 1);
4764  for (unsigned int id =
4765  dof_info.cell_loop_pre_list_index[range_index];
4766  id != dof_info.cell_loop_pre_list_index[range_index + 1];
4767  ++id)
4768  operation_before_loop(dof_info.cell_loop_pre_list[id].first,
4769  dof_info.cell_loop_pre_list[id].second);
4770  }
4771  }
4772  }
4773 
4774  virtual void
4775  cell_loop_post_range(const unsigned int range_index) override
4776  {
4777  if (operation_after_loop)
4778  {
4780  matrix_free.get_dof_info(dof_handler_index_pre_post);
4781  if (range_index == numbers::invalid_unsigned_int)
4782  {
4783  // Case with threaded loop -> currently no overlap implemented
4785  0U,
4786  dof_info.vector_partitioner->locally_owned_size(),
4787  operation_after_loop,
4789  }
4790  else
4791  {
4792  AssertIndexRange(range_index,
4793  dof_info.cell_loop_post_list_index.size() - 1);
4794  for (unsigned int id =
4795  dof_info.cell_loop_post_list_index[range_index];
4796  id != dof_info.cell_loop_post_list_index[range_index + 1];
4797  ++id)
4798  operation_after_loop(dof_info.cell_loop_post_list[id].first,
4799  dof_info.cell_loop_post_list[id].second);
4800  }
4801  }
4802  }
4803 
4804  private:
4805  const MF & matrix_free;
4806  Container & container;
4807  function_type cell_function;
4808  function_type face_function;
4809  function_type boundary_function;
4810 
4811  const InVector &src;
4812  OutVector & dst;
4813  VectorDataExchange<MF::dimension,
4814  typename MF::value_type,
4815  typename MF::vectorized_value_type>
4816  src_data_exchanger;
4817  VectorDataExchange<MF::dimension,
4818  typename MF::value_type,
4819  typename MF::vectorized_value_type>
4820  dst_data_exchanger;
4821  const bool src_and_dst_are_same;
4822  const bool zero_dst_vector_setting;
4823  const std::function<void(const unsigned int, const unsigned int)>
4824  operation_before_loop;
4825  const std::function<void(const unsigned int, const unsigned int)>
4826  operation_after_loop;
4827  const unsigned int dof_handler_index_pre_post;
4828  };
4829 
4830 
4831 
4836  template <class MF, typename InVector, typename OutVector>
4837  struct MFClassWrapper
4838  {
4839  using function_type =
4840  std::function<void(const MF &,
4841  OutVector &,
4842  const InVector &,
4843  const std::pair<unsigned int, unsigned int> &)>;
4844 
4845  MFClassWrapper(const function_type cell,
4846  const function_type face,
4847  const function_type boundary)
4848  : cell(cell)
4849  , face(face)
4850  , boundary(boundary)
4851  {}
4852 
4853  void
4854  cell_integrator(const MF & mf,
4855  OutVector & dst,
4856  const InVector & src,
4857  const std::pair<unsigned int, unsigned int> &range) const
4858  {
4859  if (cell)
4860  cell(mf, dst, src, range);
4861  }
4862 
4863  void
4864  face_integrator(const MF & mf,
4865  OutVector & dst,
4866  const InVector & src,
4867  const std::pair<unsigned int, unsigned int> &range) const
4868  {
4869  if (face)
4870  face(mf, dst, src, range);
4871  }
4872 
4873  void
4874  boundary_integrator(
4875  const MF & mf,
4876  OutVector & dst,
4877  const InVector & src,
4878  const std::pair<unsigned int, unsigned int> &range) const
4879  {
4880  if (boundary)
4881  boundary(mf, dst, src, range);
4882  }
4883 
4884  const function_type cell;
4885  const function_type face;
4886  const function_type boundary;
4887  };
4888 
4889 } // end of namespace internal
4890 
4891 
4892 
4893 template <int dim, typename Number, typename VectorizedArrayType>
4894 template <typename OutVector, typename InVector>
4895 inline void
4897  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4898  OutVector &,
4899  const InVector &,
4900  const std::pair<unsigned int, unsigned int> &)>
4901  & cell_operation,
4902  OutVector & dst,
4903  const InVector &src,
4904  const bool zero_dst_vector) const
4905 {
4906  using Wrapper =
4907  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4908  InVector,
4909  OutVector>;
4910  Wrapper wrap(cell_operation, nullptr, nullptr);
4911  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4912  InVector,
4913  OutVector,
4914  Wrapper,
4915  true>
4916  worker(*this,
4917  src,
4918  dst,
4919  zero_dst_vector,
4920  wrap,
4921  &Wrapper::cell_integrator,
4922  &Wrapper::face_integrator,
4923  &Wrapper::boundary_integrator);
4924 
4925  task_info.loop(worker);
4926 }
4927 
4928 
4929 
4930 template <int dim, typename Number, typename VectorizedArrayType>
4931 template <typename OutVector, typename InVector>
4932 inline void
4934  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4935  OutVector &,
4936  const InVector &,
4937  const std::pair<unsigned int, unsigned int> &)>
4938  & cell_operation,
4939  OutVector & dst,
4940  const InVector &src,
4941  const std::function<void(const unsigned int, const unsigned int)>
4942  &operation_before_loop,
4943  const std::function<void(const unsigned int, const unsigned int)>
4944  & operation_after_loop,
4945  const unsigned int dof_handler_index_pre_post) const
4946 {
4947  using Wrapper =
4948  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4949  InVector,
4950  OutVector>;
4951  Wrapper wrap(cell_operation, nullptr, nullptr);
4952  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4953  InVector,
4954  OutVector,
4955  Wrapper,
4956  true>
4957  worker(*this,
4958  src,
4959  dst,
4960  false,
4961  wrap,
4962  &Wrapper::cell_integrator,
4963  &Wrapper::face_integrator,
4964  &Wrapper::boundary_integrator,
4967  operation_before_loop,
4968  operation_after_loop,
4969  dof_handler_index_pre_post);
4970 
4971  task_info.loop(worker);
4972 }
4973 
4974 
4975 
4976 template <int dim, typename Number, typename VectorizedArrayType>
4977 template <typename OutVector, typename InVector>
4978 inline void
4980  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4981  OutVector &,
4982  const InVector &,
4983  const std::pair<unsigned int, unsigned int> &)>
4984  &cell_operation,
4985  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4986  OutVector &,
4987  const InVector &,
4988  const std::pair<unsigned int, unsigned int> &)>
4989  &face_operation,
4990  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4991  OutVector &,
4992  const InVector &,
4993  const std::pair<unsigned int, unsigned int> &)>
4994  & boundary_operation,
4995  OutVector & dst,
4996  const InVector & src,
4997  const bool zero_dst_vector,
4998  const DataAccessOnFaces dst_vector_face_access,
4999  const DataAccessOnFaces src_vector_face_access) const
5000 {
5001  using Wrapper =
5002  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5003  InVector,
5004  OutVector>;
5005  Wrapper wrap(cell_operation, face_operation, boundary_operation);
5006  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5007  InVector,
5008  OutVector,
5009  Wrapper,
5010  true>
5011  worker(*this,
5012  src,
5013  dst,
5014  zero_dst_vector,
5015  wrap,
5016  &Wrapper::cell_integrator,
5017  &Wrapper::face_integrator,
5018  &Wrapper::boundary_integrator,
5019  src_vector_face_access,
5020  dst_vector_face_access);
5021 
5022  task_info.loop(worker);
5023 }
5024 
5025 
5026 
5027 template <int dim, typename Number, typename VectorizedArrayType>
5028 template <typename CLASS, typename OutVector, typename InVector>
5029 inline void
5031  void (CLASS::*function_pointer)(
5033  OutVector &,
5034  const InVector &,
5035  const std::pair<unsigned int, unsigned int> &) const,
5036  const CLASS * owning_class,
5037  OutVector & dst,
5038  const InVector &src,
5039  const bool zero_dst_vector) const
5040 {
5041  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5042  InVector,
5043  OutVector,
5044  CLASS,
5045  true>
5046  worker(*this,
5047  src,
5048  dst,
5049  zero_dst_vector,
5050  *owning_class,
5051  function_pointer,
5052  nullptr,
5053  nullptr);
5054  task_info.loop(worker);
5055 }
5056 
5057 
5058 
5059 template <int dim, typename Number, typename VectorizedArrayType>
5060 template <typename CLASS, typename OutVector, typename InVector>
5061 inline void
5063  void (CLASS::*function_pointer)(
5065  OutVector &,
5066  const InVector &,
5067  const std::pair<unsigned int, unsigned int> &) const,
5068  const CLASS * owning_class,
5069  OutVector & dst,
5070  const InVector &src,
5071  const std::function<void(const unsigned int, const unsigned int)>
5072  &operation_before_loop,
5073  const std::function<void(const unsigned int, const unsigned int)>
5074  & operation_after_loop,
5075  const unsigned int dof_handler_index_pre_post) const
5076 {
5077  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5078  InVector,
5079  OutVector,
5080  CLASS,
5081  true>
5082  worker(*this,
5083  src,
5084  dst,
5085  false,
5086  *owning_class,
5087  function_pointer,
5088  nullptr,
5089  nullptr,
5092  operation_before_loop,
5093  operation_after_loop,
5094  dof_handler_index_pre_post);
5095  task_info.loop(worker);
5096 }
5097 
5098 
5099 
5100 template <int dim, typename Number, typename VectorizedArrayType>
5101 template <typename CLASS, typename OutVector, typename InVector>
5102 inline void
5104  void (CLASS::*cell_operation)(
5106  OutVector &,
5107  const InVector &,
5108  const std::pair<unsigned int, unsigned int> &) const,
5109  void (CLASS::*face_operation)(
5111  OutVector &,
5112  const InVector &,
5113  const std::pair<unsigned int, unsigned int> &) const,
5114  void (CLASS::*boundary_operation)(
5116  OutVector &,
5117  const InVector &,
5118  const std::pair<unsigned int, unsigned int> &) const,
5119  const CLASS * owning_class,
5120  OutVector & dst,
5121  const InVector & src,
5122  const bool zero_dst_vector,
5123  const DataAccessOnFaces dst_vector_face_access,
5124  const DataAccessOnFaces src_vector_face_access) const
5125 {
5126  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5127  InVector,
5128  OutVector,
5129  CLASS,
5130  true>
5131  worker(*this,
5132  src,
5133  dst,
5134  zero_dst_vector,
5135  *owning_class,
5136  cell_operation,
5137  face_operation,
5138  boundary_operation,
5139  src_vector_face_access,
5140  dst_vector_face_access);
5141  task_info.loop(worker);
5142 }
5143 
5144 
5145 
5146 template <int dim, typename Number, typename VectorizedArrayType>
5147 template <typename CLASS, typename OutVector, typename InVector>
5148 inline void
5150  void (CLASS::*function_pointer)(
5152  OutVector &,
5153  const InVector &,
5154  const std::pair<unsigned int, unsigned int> &),
5155  CLASS * owning_class,
5156  OutVector & dst,
5157  const InVector &src,
5158  const bool zero_dst_vector) const
5159 {
5160  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5161  InVector,
5162  OutVector,
5163  CLASS,
5164  false>
5165  worker(*this,
5166  src,
5167  dst,
5168  zero_dst_vector,
5169  *owning_class,
5170  function_pointer,
5171  nullptr,
5172  nullptr);
5173  task_info.loop(worker);
5174 }
5175 
5176 
5177 
5178 template <int dim, typename Number, typename VectorizedArrayType>
5179 template <typename CLASS, typename OutVector, typename InVector>
5180 inline void
5182  void (CLASS::*function_pointer)(
5184  OutVector &,
5185  const InVector &,
5186  const std::pair<unsigned int, unsigned int> &),
5187  CLASS * owning_class,
5188  OutVector & dst,
5189  const InVector &src,
5190  const std::function<void(const unsigned int, const unsigned int)>
5191  &operation_before_loop,
5192  const std::function<void(const unsigned int, const unsigned int)>
5193  & operation_after_loop,
5194  const unsigned int dof_handler_index_pre_post) const
5195 {
5196  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5197  InVector,
5198  OutVector,
5199  CLASS,
5200  false>
5201  worker(*this,
5202  src,
5203  dst,
5204  false,
5205  *owning_class,
5206  function_pointer,
5207  nullptr,
5208  nullptr,
5211  operation_before_loop,
5212  operation_after_loop,
5213  dof_handler_index_pre_post);
5214  task_info.loop(worker);
5215 }
5216 
5217 
5218 
5219 template <int dim, typename Number, typename VectorizedArrayType>
5220 template <typename CLASS, typename OutVector, typename InVector>
5221 inline void
5223  void (CLASS::*cell_operation)(
5225  OutVector &,
5226  const InVector &,
5227  const std::pair<unsigned int, unsigned int> &),
5228  void (CLASS::*face_operation)(
5230  OutVector &,
5231  const InVector &,
5232  const std::pair<unsigned int, unsigned int> &),
5233  void (CLASS::*boundary_operation)(
5235  OutVector &,
5236  const InVector &,
5237  const std::pair<unsigned int, unsigned int> &),
5238  CLASS * owning_class,
5239  OutVector & dst,
5240  const InVector & src,
5241  const bool zero_dst_vector,
5242  const DataAccessOnFaces dst_vector_face_access,
5243  const DataAccessOnFaces src_vector_face_access) const
5244 {
5245  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5246  InVector,
5247  OutVector,
5248  CLASS,
5249  false>
5250  worker(*this,
5251  src,
5252  dst,
5253  zero_dst_vector,
5254  *owning_class,
5255  cell_operation,
5256  face_operation,
5257  boundary_operation,
5258  src_vector_face_access,
5259  dst_vector_face_access);
5260  task_info.loop(worker);
5261 }
5262 
5263 
5264 
5265 template <int dim, typename Number, typename VectorizedArrayType>
5266 template <typename CLASS, typename OutVector, typename InVector>
5267 inline void
5269  void (CLASS::*function_pointer)(
5271  OutVector &,
5272  const InVector &,
5273  const std::pair<unsigned int, unsigned int> &) const,
5274  const CLASS * owning_class,
5275  OutVector & dst,
5276  const InVector & src,
5277  const bool zero_dst_vector,
5278  const DataAccessOnFaces src_vector_face_access) const
5279 {
5280  auto src_vector_face_access_temp = src_vector_face_access;
5281  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5282  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5283  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5284  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5285 
5286  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5287  InVector,
5288  OutVector,
5289  CLASS,
5290  true>
5291  worker(*this,
5292  src,
5293  dst,
5294  zero_dst_vector,
5295  *owning_class,
5296  function_pointer,
5297  nullptr,
5298  nullptr,
5299  src_vector_face_access_temp,
5301  task_info.loop(worker);
5302 }
5303 
5304 
5305 
5306 template <int dim, typename Number, typename VectorizedArrayType>
5307 template <typename CLASS, typename OutVector, typename InVector>
5308 inline void
5310  void (CLASS::*function_pointer)(
5312  OutVector &,
5313  const InVector &,
5314  const std::pair<unsigned int, unsigned int> &),
5315  CLASS * owning_class,
5316  OutVector & dst,
5317  const InVector & src,
5318  const bool zero_dst_vector,
5319  const DataAccessOnFaces src_vector_face_access) const
5320 {
5321  auto src_vector_face_access_temp = src_vector_face_access;
5322  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5323  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5324  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5325  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5326 
5327  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5328  InVector,
5329  OutVector,
5330  CLASS,
5331  false>
5332  worker(*this,
5333  src,
5334  dst,
5335  zero_dst_vector,
5336  *owning_class,
5337  function_pointer,
5338  nullptr,
5339  nullptr,
5340  src_vector_face_access_temp,
5342  task_info.loop(worker);
5343 }
5344 
5345 
5346 
5347 template <int dim, typename Number, typename VectorizedArrayType>
5348 template <typename OutVector, typename InVector>
5349 inline void
5351  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5352  OutVector &,
5353  const InVector &,
5354  const std::pair<unsigned int, unsigned int> &)>
5355  & cell_operation,
5356  OutVector & dst,
5357  const InVector & src,
5358  const bool zero_dst_vector,
5359  const DataAccessOnFaces src_vector_face_access) const
5360 {
5361  auto src_vector_face_access_temp = src_vector_face_access;
5362  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5363  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5364  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5365  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5366 
5367  using Wrapper =
5368  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5369  InVector,
5370  OutVector>;
5371  Wrapper wrap(cell_operation, nullptr, nullptr);
5372 
5373  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5374  InVector,
5375  OutVector,
5376  Wrapper,
5377  true>
5378  worker(*this,
5379  src,
5380  dst,
5381  zero_dst_vector,
5382  wrap,
5383  &Wrapper::cell_integrator,
5384  &Wrapper::face_integrator,
5385  &Wrapper::boundary_integrator,
5386  src_vector_face_access_temp,
5388  task_info.loop(worker);
5389 }
5390 
5391 
5392 #endif // ifndef DOXYGEN
5393 
5394 
5395 
5397 
5398 #endif
Transformed quadrature weights.
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > > shape_info
Definition: matrix_free.h:2132
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:747
AlignedVector< Number > * acquire_scratch_data_non_threadsafe() const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1655
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:754
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:2173
#define AssertIndexRange(index, range)
Definition: exceptions.h:1720
AdditionalData(const AdditionalData &other)
Definition: matrix_free.h:253
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition: dof_info.h:598
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:742
std::list< std::pair< bool, AlignedVector< Number > > > scratch_pad_non_threadsafe
Definition: matrix_free.h:2191
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:388
static ::ExceptionBase & ExcNotInitialized()
std::vector< unsigned int > cell_partition_data
Definition: task_info.h:474
#define AssertThrow(cond, exc)
Definition: exceptions.h:1571
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
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
static const unsigned int dimension
Definition: matrix_free.h:133
std::vector< unsigned int > cell_loop_post_list_index
Definition: dof_info.h:767
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
AdditionalData(const TasksParallelScheme tasks_parallel_scheme=partition_partition, const unsigned int tasks_block_size=0, const UpdateFlags mapping_update_flags=update_gradients|update_JxW_values, const UpdateFlags mapping_update_flags_boundary_faces=update_default, const UpdateFlags mapping_update_flags_inner_faces=update_default, const UpdateFlags mapping_update_flags_faces_by_cells=update_default, const unsigned int mg_level=numbers::invalid_unsigned_int, const bool store_plain_indices=true, const bool initialize_indices=true, const bool initialize_mapping=true, const bool overlap_communication_computation=true, const bool hold_all_faces_to_owned_cells=false, const bool cell_vectorization_categories_strict=false, const bool allow_ghosted_vectors_in_loops=true, const bool use_fast_hanging_node_algorithm=true)
Definition: matrix_free.h:214
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_pre_list
Definition: dof_info.h:760
unsigned int n_ghost_inner_face_batches() const
#define Assert(cond, exc)
Definition: exceptions.h:1461
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:2163
const types::boundary_id invalid_boundary_id
Definition: types.h:239
unsigned int mg_level
Definition: matrix_free.h:2196
std::vector< std::pair< unsigned int, unsigned int > > cell_level_index
Definition: matrix_free.h:2140
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:408
std::vector< unsigned int > boundary_partition_data
Definition: task_info.h:518
std::vector< internal::MatrixFreeFunctions::DoFInfo > dof_info
Definition: matrix_free.h:2105
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
unsigned int level
Definition: grid_out.cc:4589
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:436
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:515
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:2113
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:2156
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
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1778
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_NAMESPACE_OPEN
Definition: config.h:400
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
VectorType::value_type * begin(VectorType &V)
Shape function gradients.
bool indices_are_initialized
Definition: matrix_free.h:2168
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:2184
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:2126
AdditionalData & operator=(const AdditionalData &other)
Definition: matrix_free.h:281
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:2119
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_post_list
Definition: dof_info.h:773
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:2099
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:985
#define DEAL_II_DEPRECATED
Definition: config.h:160
unsigned int boundary_id
Definition: types.h:129
unsigned int cell_level_index_end_local
Definition: matrix_free.h:2149
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()