Reference documentation for deal.II version GIT 7bd1c95dbb 2022-06-29 20:50:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
matrix_free.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2011 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii_matrix_free_h
18 #define dealii_matrix_free_h
19 
20 #include <deal.II/base/config.h>
21 
28 
30 
31 #include <deal.II/fe/fe.h>
32 #include <deal.II/fe/mapping.h>
33 #include <deal.II/fe/mapping_q1.h>
34 
37 
42 
48 
49 #include <cstdlib>
50 #include <limits>
51 #include <list>
52 #include <memory>
53 
54 
56 
57 
58 
110 template <int dim,
111  typename Number = double,
112  typename VectorizedArrayType = VectorizedArray<Number>>
113 class MatrixFree : public Subscriptor
114 {
115  static_assert(
116  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
117  "Type of Number and of VectorizedArrayType do not match.");
118 
119 public:
124  using value_type = Number;
125  using vectorized_value_type = VectorizedArrayType;
126 
130  static constexpr unsigned int dimension = dim;
131 
179  {
184 
191  {
211  };
212 
218  const unsigned int tasks_block_size = 0,
224  const unsigned int mg_level = numbers::invalid_unsigned_int,
225  const bool store_plain_indices = true,
226  const bool initialize_indices = true,
227  const bool initialize_mapping = true,
228  const bool overlap_communication_computation = true,
229  const bool hold_all_faces_to_owned_cells = false,
230  const bool cell_vectorization_categories_strict = false,
231  const bool allow_ghosted_vectors_in_loops = true)
238  , mg_level(mg_level)
247  , communicator_sm(MPI_COMM_SELF)
248  {}
249 
262  , mg_level(other.mg_level)
274  {}
275 
280  operator=(const AdditionalData &other)
281  {
290  mg_level = other.mg_level;
302 
303  return *this;
304  }
305 
343 
353  unsigned int tasks_block_size;
354 
367 
387 
407 
435 
444  unsigned int mg_level;
445 
453 
464 
473 
487 
496 
518  std::vector<unsigned int> cell_vectorization_category;
519 
530 
542 
546  MPI_Comm communicator_sm;
547  };
548 
557 
562 
566  ~MatrixFree() override = default;
567 
580  template <typename QuadratureType, typename number2, typename MappingType>
581  void
582  reinit(const MappingType & mapping,
583  const DoFHandler<dim> & dof_handler,
584  const AffineConstraints<number2> &constraint,
585  const QuadratureType & quad,
586  const AdditionalData & additional_data = AdditionalData());
587 
609  template <typename QuadratureType, typename number2, typename MappingType>
610  void
611  reinit(const MappingType & mapping,
612  const std::vector<const DoFHandler<dim> *> & dof_handler,
613  const std::vector<const AffineConstraints<number2> *> &constraint,
614  const std::vector<QuadratureType> & quad,
615  const AdditionalData &additional_data = AdditionalData());
616 
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 QuadratureType & quad,
630  const AdditionalData &additional_data = AdditionalData());
631 
637  void
639  const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free_base);
640 
650  void
651  update_mapping(const Mapping<dim> &mapping);
652 
656  void
657  update_mapping(const std::shared_ptr<hp::MappingCollection<dim>> &mapping);
658 
663  void
664  clear();
665 
667 
679  enum class DataAccessOnFaces
680  {
687  none,
688 
699  values,
700 
710 
721  gradients,
722 
732 
740  };
741 
786  template <typename OutVector, typename InVector>
787  void
788  cell_loop(const std::function<void(
790  OutVector &,
791  const InVector &,
792  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
793  OutVector & dst,
794  const InVector & src,
795  const bool zero_dst_vector = false) const;
796 
843  template <typename CLASS, typename OutVector, typename InVector>
844  void
845  cell_loop(void (CLASS::*cell_operation)(
846  const MatrixFree &,
847  OutVector &,
848  const InVector &,
849  const std::pair<unsigned int, unsigned int> &) const,
850  const CLASS * owning_class,
851  OutVector & dst,
852  const InVector &src,
853  const bool zero_dst_vector = false) const;
854 
858  template <typename CLASS, typename OutVector, typename InVector>
859  void
860  cell_loop(void (CLASS::*cell_operation)(
861  const MatrixFree &,
862  OutVector &,
863  const InVector &,
864  const std::pair<unsigned int, unsigned int> &),
865  CLASS * owning_class,
866  OutVector & dst,
867  const InVector &src,
868  const bool zero_dst_vector = false) const;
869 
954  template <typename CLASS, typename OutVector, typename InVector>
955  void
956  cell_loop(void (CLASS::*cell_operation)(
957  const MatrixFree &,
958  OutVector &,
959  const InVector &,
960  const std::pair<unsigned int, unsigned int> &) const,
961  const CLASS * owning_class,
962  OutVector & dst,
963  const InVector &src,
964  const std::function<void(const unsigned int, const unsigned int)>
965  &operation_before_loop,
966  const std::function<void(const unsigned int, const unsigned int)>
967  & operation_after_loop,
968  const unsigned int dof_handler_index_pre_post = 0) const;
969 
973  template <typename CLASS, typename OutVector, typename InVector>
974  void
975  cell_loop(void (CLASS::*cell_operation)(
976  const MatrixFree &,
977  OutVector &,
978  const InVector &,
979  const std::pair<unsigned int, unsigned int> &),
980  CLASS * owning_class,
981  OutVector & dst,
982  const InVector &src,
983  const std::function<void(const unsigned int, const unsigned int)>
984  &operation_before_loop,
985  const std::function<void(const unsigned int, const unsigned int)>
986  & operation_after_loop,
987  const unsigned int dof_handler_index_pre_post = 0) const;
988 
993  template <typename OutVector, typename InVector>
994  void
995  cell_loop(const std::function<void(
997  OutVector &,
998  const InVector &,
999  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1000  OutVector & dst,
1001  const InVector & src,
1002  const std::function<void(const unsigned int, const unsigned int)>
1003  &operation_before_loop,
1004  const std::function<void(const unsigned int, const unsigned int)>
1005  & operation_after_loop,
1006  const unsigned int dof_handler_index_pre_post = 0) const;
1007 
1083  template <typename OutVector, typename InVector>
1084  void
1085  loop(const std::function<
1087  OutVector &,
1088  const InVector &,
1089  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1090  const std::function<
1092  OutVector &,
1093  const InVector &,
1094  const std::pair<unsigned int, unsigned int> &)> &face_operation,
1095  const std::function<void(
1097  OutVector &,
1098  const InVector &,
1099  const std::pair<unsigned int, unsigned int> &)> &boundary_operation,
1100  OutVector & dst,
1101  const InVector & src,
1102  const bool zero_dst_vector = false,
1103  const DataAccessOnFaces dst_vector_face_access =
1105  const DataAccessOnFaces src_vector_face_access =
1107 
1193  template <typename CLASS, typename OutVector, typename InVector>
1194  void
1196  void (CLASS::*cell_operation)(const MatrixFree &,
1197  OutVector &,
1198  const InVector &,
1199  const std::pair<unsigned int, unsigned int> &)
1200  const,
1201  void (CLASS::*face_operation)(const MatrixFree &,
1202  OutVector &,
1203  const InVector &,
1204  const std::pair<unsigned int, unsigned int> &)
1205  const,
1206  void (CLASS::*boundary_operation)(
1207  const MatrixFree &,
1208  OutVector &,
1209  const InVector &,
1210  const std::pair<unsigned int, unsigned int> &) const,
1211  const CLASS * owning_class,
1212  OutVector & dst,
1213  const InVector & src,
1214  const bool zero_dst_vector = false,
1215  const DataAccessOnFaces dst_vector_face_access =
1217  const DataAccessOnFaces src_vector_face_access =
1219 
1223  template <typename CLASS, typename OutVector, typename InVector>
1224  void
1225  loop(void (CLASS::*cell_operation)(
1226  const MatrixFree &,
1227  OutVector &,
1228  const InVector &,
1229  const std::pair<unsigned int, unsigned int> &),
1230  void (CLASS::*face_operation)(
1231  const MatrixFree &,
1232  OutVector &,
1233  const InVector &,
1234  const std::pair<unsigned int, unsigned int> &),
1235  void (CLASS::*boundary_operation)(
1236  const MatrixFree &,
1237  OutVector &,
1238  const InVector &,
1239  const std::pair<unsigned int, unsigned int> &),
1240  CLASS * owning_class,
1241  OutVector & dst,
1242  const InVector & src,
1243  const bool zero_dst_vector = false,
1244  const DataAccessOnFaces dst_vector_face_access =
1246  const DataAccessOnFaces src_vector_face_access =
1248 
1313  template <typename CLASS, typename OutVector, typename InVector>
1314  void
1315  loop_cell_centric(void (CLASS::*cell_operation)(
1316  const MatrixFree &,
1317  OutVector &,
1318  const InVector &,
1319  const std::pair<unsigned int, unsigned int> &) const,
1320  const CLASS * owning_class,
1321  OutVector & dst,
1322  const InVector & src,
1323  const bool zero_dst_vector = false,
1324  const DataAccessOnFaces src_vector_face_access =
1326 
1330  template <typename CLASS, typename OutVector, typename InVector>
1331  void
1332  loop_cell_centric(void (CLASS::*cell_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 src_vector_face_access =
1343 
1347  template <typename OutVector, typename InVector>
1348  void
1350  const std::function<void(const MatrixFree &,
1351  OutVector &,
1352  const InVector &,
1353  const std::pair<unsigned int, unsigned int> &)>
1354  & cell_operation,
1355  OutVector & dst,
1356  const InVector & src,
1357  const bool zero_dst_vector = false,
1358  const DataAccessOnFaces src_vector_face_access =
1360 
1368  std::pair<unsigned int, unsigned int>
1369  create_cell_subrange_hp(const std::pair<unsigned int, unsigned int> &range,
1370  const unsigned int fe_degree,
1371  const unsigned int dof_handler_index = 0) const;
1372 
1379  std::pair<unsigned int, unsigned int>
1381  const std::pair<unsigned int, unsigned int> &range,
1382  const unsigned int fe_index,
1383  const unsigned int dof_handler_index = 0) const;
1384 
1388  unsigned int
1390 
1394  unsigned int
1396  const std::pair<unsigned int, unsigned int> range) const;
1397 
1401  unsigned int
1402  get_face_active_fe_index(const std::pair<unsigned int, unsigned int> range,
1403  const bool is_interior_face = true) const;
1404 
1406 
1416  template <typename T>
1417  void
1419 
1425  template <typename T>
1426  void
1428 
1442  template <typename VectorType>
1443  void
1444  initialize_dof_vector(VectorType & vec,
1445  const unsigned int dof_handler_index = 0) const;
1446 
1464  template <typename Number2>
1465  void
1467  const unsigned int dof_handler_index = 0) const;
1468 
1479  const std::shared_ptr<const Utilities::MPI::Partitioner> &
1480  get_vector_partitioner(const unsigned int dof_handler_index = 0) const;
1481 
1485  const IndexSet &
1486  get_locally_owned_set(const unsigned int dof_handler_index = 0) const;
1487 
1491  const IndexSet &
1492  get_ghost_set(const unsigned int dof_handler_index = 0) const;
1493 
1503  const std::vector<unsigned int> &
1504  get_constrained_dofs(const unsigned int dof_handler_index = 0) const;
1505 
1516  void
1517  renumber_dofs(std::vector<types::global_dof_index> &renumbering,
1518  const unsigned int dof_handler_index = 0);
1519 
1521 
1529  template <int spacedim>
1530  static bool
1532 
1536  unsigned int
1537  n_components() const;
1538 
1543  unsigned int
1544  n_base_elements(const unsigned int dof_handler_index) const;
1545 
1553  unsigned int
1555 
1565  unsigned int
1567 
1574  unsigned int
1576 
1586  unsigned int
1588 
1599  unsigned int
1601 
1606  unsigned int
1608 
1616  get_boundary_id(const unsigned int face_batch_index) const;
1617 
1622  std::array<types::boundary_id, VectorizedArrayType::size()>
1623  get_faces_by_cells_boundary_id(const unsigned int cell_batch_index,
1624  const unsigned int face_number) const;
1625 
1630  const DoFHandler<dim> &
1631  get_dof_handler(const unsigned int dof_handler_index = 0) const;
1632 
1646  get_cell_iterator(const unsigned int cell_batch_index,
1647  const unsigned int lane_index,
1648  const unsigned int dof_handler_index = 0) const;
1649 
1655  std::pair<int, int>
1656  get_cell_level_and_index(const unsigned int cell_batch_index,
1657  const unsigned int lane_index) const;
1658 
1671  std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>
1672  get_face_iterator(const unsigned int face_batch_index,
1673  const unsigned int lane_index,
1674  const bool interior = true,
1675  const unsigned int fe_component = 0) const;
1676 
1689  bool
1690  at_irregular_cell(const unsigned int cell_batch_index) const;
1691 
1701  unsigned int
1702  n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const;
1703 
1713  unsigned int
1714  n_active_entries_per_face_batch(const unsigned int face_batch_index) const;
1715 
1719  unsigned int
1720  get_dofs_per_cell(const unsigned int dof_handler_index = 0,
1721  const unsigned int hp_active_fe_index = 0) const;
1722 
1726  unsigned int
1727  get_n_q_points(const unsigned int quad_index = 0,
1728  const unsigned int hp_active_fe_index = 0) const;
1729 
1734  unsigned int
1735  get_dofs_per_face(const unsigned int dof_handler_index = 0,
1736  const unsigned int hp_active_fe_index = 0) const;
1737 
1742  unsigned int
1743  get_n_q_points_face(const unsigned int quad_index = 0,
1744  const unsigned int hp_active_fe_index = 0) const;
1745 
1749  const Quadrature<dim> &
1750  get_quadrature(const unsigned int quad_index = 0,
1751  const unsigned int hp_active_fe_index = 0) const;
1752 
1756  const Quadrature<dim - 1> &
1757  get_face_quadrature(const unsigned int quad_index = 0,
1758  const unsigned int hp_active_fe_index = 0) const;
1759 
1772  unsigned int
1774  const std::pair<unsigned int, unsigned int> cell_batch_range) const;
1775 
1780  std::pair<unsigned int, unsigned int>
1782  const std::pair<unsigned int, unsigned int> face_batch_range) const;
1783 
1795  unsigned int
1796  get_cell_category(const unsigned int cell_batch_index) const;
1797 
1802  std::pair<unsigned int, unsigned int>
1803  get_face_category(const unsigned int face_batch_index) const;
1804 
1808  bool
1810 
1815  bool
1817 
1822  unsigned int
1823  get_mg_level() const;
1824 
1829  std::size_t
1831 
1836  template <typename StreamType>
1837  void
1838  print_memory_consumption(StreamType &out) const;
1839 
1844  void
1845  print(std::ostream &out) const;
1846 
1848 
1859  get_task_info() const;
1860 
1861  /*
1862  * Return geometry-dependent information on the cells.
1863  */
1864  const internal::MatrixFreeFunctions::
1865  MappingInfo<dim, Number, VectorizedArrayType> &
1867 
1872  get_dof_info(const unsigned int dof_handler_index_component = 0) const;
1873 
1877  unsigned int
1879 
1884  const Number *
1885  constraint_pool_begin(const unsigned int pool_index) const;
1886 
1892  const Number *
1893  constraint_pool_end(const unsigned int pool_index) const;
1894 
1899  get_shape_info(const unsigned int dof_handler_index_component = 0,
1900  const unsigned int quad_index = 0,
1901  const unsigned int fe_base_element = 0,
1902  const unsigned int hp_active_fe_index = 0,
1903  const unsigned int hp_active_quad_index = 0) const;
1904 
1909  VectorizedArrayType::size()> &
1910  get_face_info(const unsigned int face_batch_index) const;
1911 
1912 
1918  const Table<3, unsigned int> &
1920 
1936 
1940  void
1942 
1954 
1958  void
1960  const AlignedVector<Number> *memory) const;
1961 
1963 
1964 private:
1969  template <typename number2, int q_dim>
1970  void
1972  const std::shared_ptr<hp::MappingCollection<dim>> & mapping,
1973  const std::vector<const DoFHandler<dim, dim> *> & dof_handlers,
1974  const std::vector<const AffineConstraints<number2> *> &constraint,
1975  const std::vector<IndexSet> & locally_owned_set,
1976  const std::vector<hp::QCollection<q_dim>> & quad,
1977  const AdditionalData & additional_data);
1978 
1985  template <typename number2>
1986  void
1988  const std::vector<const AffineConstraints<number2> *> &constraint,
1989  const std::vector<IndexSet> & locally_owned_set,
1990  const AdditionalData & additional_data);
1991 
1995  void
1997  const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
1998  const AdditionalData & additional_data);
1999 
2003  std::vector<SmartPointer<const DoFHandler<dim>>> dof_handlers;
2004 
2009  std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
2010 
2017  std::vector<Number> constraint_pool_data;
2018 
2023  std::vector<unsigned int> constraint_pool_row_index;
2024 
2031 
2037 
2044  std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
2045 
2046 
2054 
2061 
2066  internal::MatrixFreeFunctions::FaceInfo<VectorizedArrayType::size()>
2068 
2073 
2078 
2087  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2089 
2094  mutable std::list<std::pair<bool, AlignedVector<Number>>>
2096 
2100  unsigned int mg_level;
2101 };
2102 
2103 
2104 
2105 /*----------------------- Inline functions ----------------------------------*/
2106 
2107 #ifndef DOXYGEN
2108 
2109 
2110 
2111 template <int dim, typename Number, typename VectorizedArrayType>
2112 template <typename T>
2113 inline void
2115  AlignedVector<T> &vec) const
2116 {
2117  vec.resize(this->n_cell_batches() + this->n_ghost_cell_batches());
2118 }
2119 
2120 
2121 
2122 template <int dim, typename Number, typename VectorizedArrayType>
2123 template <typename T>
2124 inline void
2126  AlignedVector<T> &vec) const
2127 {
2128  vec.resize(this->n_inner_face_batches() + this->n_boundary_face_batches() +
2129  this->n_ghost_inner_face_batches());
2130 }
2131 
2132 
2133 
2134 template <int dim, typename Number, typename VectorizedArrayType>
2135 template <typename VectorType>
2136 inline void
2138  VectorType & vec,
2139  const unsigned int comp) const
2140 {
2141  static_assert(IsBlockVector<VectorType>::value == false,
2142  "This function is not supported for block vectors.");
2143 
2144  Assert(task_info.n_procs == 1,
2145  ExcMessage("This function can only be used in serial."));
2146 
2147  AssertIndexRange(comp, n_components());
2148  vec.reinit(dof_info[comp].vector_partitioner->size());
2149 }
2150 
2151 
2152 
2153 template <int dim, typename Number, typename VectorizedArrayType>
2154 template <typename Number2>
2155 inline void
2158  const unsigned int comp) const
2159 {
2160  AssertIndexRange(comp, n_components());
2161  vec.reinit(dof_info[comp].vector_partitioner, task_info.communicator_sm);
2162 }
2163 
2164 
2165 
2166 template <int dim, typename Number, typename VectorizedArrayType>
2167 inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2169  const unsigned int comp) const
2170 {
2171  AssertIndexRange(comp, n_components());
2172  return dof_info[comp].vector_partitioner;
2173 }
2174 
2175 
2176 
2177 template <int dim, typename Number, typename VectorizedArrayType>
2178 inline const std::vector<unsigned int> &
2180  const unsigned int comp) const
2181 {
2182  AssertIndexRange(comp, n_components());
2183  return dof_info[comp].constrained_dofs;
2184 }
2185 
2186 
2187 
2188 template <int dim, typename Number, typename VectorizedArrayType>
2189 inline unsigned int
2191 {
2192  AssertDimension(dof_handlers.size(), dof_info.size());
2193  return dof_handlers.size();
2194 }
2195 
2196 
2197 
2198 template <int dim, typename Number, typename VectorizedArrayType>
2199 inline unsigned int
2201  const unsigned int dof_no) const
2202 {
2203  AssertDimension(dof_handlers.size(), dof_info.size());
2204  AssertIndexRange(dof_no, dof_handlers.size());
2205  return dof_handlers[dof_no]->get_fe().n_base_elements();
2206 }
2207 
2208 
2209 
2210 template <int dim, typename Number, typename VectorizedArrayType>
2213 {
2214  return task_info;
2215 }
2216 
2217 
2218 
2219 template <int dim, typename Number, typename VectorizedArrayType>
2220 inline unsigned int
2222 {
2223  return task_info.n_active_cells;
2224 }
2225 
2226 
2227 
2228 template <int dim, typename Number, typename VectorizedArrayType>
2229 inline unsigned int
2231 {
2232  return *(task_info.cell_partition_data.end() - 2);
2233 }
2234 
2235 
2236 
2237 template <int dim, typename Number, typename VectorizedArrayType>
2238 inline unsigned int
2240 {
2241  return *(task_info.cell_partition_data.end() - 1) -
2242  *(task_info.cell_partition_data.end() - 2);
2243 }
2244 
2245 
2246 
2247 template <int dim, typename Number, typename VectorizedArrayType>
2248 inline unsigned int
2250 {
2251  if (task_info.face_partition_data.size() == 0)
2252  return 0;
2253  return task_info.face_partition_data.back();
2254 }
2255 
2256 
2257 
2258 template <int dim, typename Number, typename VectorizedArrayType>
2259 inline unsigned int
2261 {
2262  if (task_info.face_partition_data.size() == 0)
2263  return 0;
2264  return task_info.boundary_partition_data.back() -
2266 }
2267 
2268 
2269 
2270 template <int dim, typename Number, typename VectorizedArrayType>
2271 inline unsigned int
2273 {
2274  if (task_info.face_partition_data.size() == 0)
2275  return 0;
2276  return face_info.faces.size() - task_info.boundary_partition_data.back();
2277 }
2278 
2279 
2280 
2281 template <int dim, typename Number, typename VectorizedArrayType>
2282 inline types::boundary_id
2284  const unsigned int face_batch_index) const
2285 {
2286  Assert(face_batch_index >= task_info.boundary_partition_data[0] &&
2287  face_batch_index < task_info.boundary_partition_data.back(),
2288  ExcIndexRange(face_batch_index,
2291  return types::boundary_id(face_info.faces[face_batch_index].exterior_face_no);
2292 }
2293 
2294 
2295 
2296 template <int dim, typename Number, typename VectorizedArrayType>
2297 inline std::array<types::boundary_id, VectorizedArrayType::size()>
2299  const unsigned int cell_batch_index,
2300  const unsigned int face_number) const
2301 {
2302  AssertIndexRange(cell_batch_index, n_cell_batches());
2305  ExcNotInitialized());
2306  std::array<types::boundary_id, VectorizedArrayType::size()> result;
2307  result.fill(numbers::invalid_boundary_id);
2308  for (unsigned int v = 0;
2309  v < n_active_entries_per_cell_batch(cell_batch_index);
2310  ++v)
2311  result[v] =
2312  face_info.cell_and_face_boundary_id(cell_batch_index, face_number, v);
2313  return result;
2314 }
2315 
2316 
2317 
2318 template <int dim, typename Number, typename VectorizedArrayType>
2319 inline const internal::MatrixFreeFunctions::
2320  MappingInfo<dim, Number, VectorizedArrayType> &
2322 {
2323  return mapping_info;
2324 }
2325 
2326 
2327 
2328 template <int dim, typename Number, typename VectorizedArrayType>
2331  const unsigned int dof_index) const
2332 {
2333  AssertIndexRange(dof_index, n_components());
2334  return dof_info[dof_index];
2335 }
2336 
2337 
2338 
2339 template <int dim, typename Number, typename VectorizedArrayType>
2340 inline unsigned int
2342 {
2343  return constraint_pool_row_index.size() - 1;
2344 }
2345 
2346 
2347 
2348 template <int dim, typename Number, typename VectorizedArrayType>
2349 inline const Number *
2351  const unsigned int row) const
2352 {
2354  return constraint_pool_data.empty() ?
2355  nullptr :
2357 }
2358 
2359 
2360 
2361 template <int dim, typename Number, typename VectorizedArrayType>
2362 inline const Number *
2364  const unsigned int row) const
2365 {
2367  return constraint_pool_data.empty() ?
2368  nullptr :
2370 }
2371 
2372 
2373 
2374 template <int dim, typename Number, typename VectorizedArrayType>
2375 inline std::pair<unsigned int, unsigned int>
2377  const std::pair<unsigned int, unsigned int> &range,
2378  const unsigned int degree,
2379  const unsigned int dof_handler_component) const
2380 {
2381  if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2382  {
2384  dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2386  dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2387  if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2388  return range;
2389  else
2390  return {range.second, range.second};
2391  }
2392 
2393  const unsigned int fe_index =
2394  dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2395  if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2396  return {range.second, range.second};
2397  else
2398  return create_cell_subrange_hp_by_index(range,
2399  fe_index,
2400  dof_handler_component);
2401 }
2402 
2403 
2404 
2405 template <int dim, typename Number, typename VectorizedArrayType>
2406 inline bool
2408  const unsigned int cell_batch_index) const
2409 {
2410  AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2411  return VectorizedArrayType::size() > 1 &&
2412  cell_level_index[(cell_batch_index + 1) * VectorizedArrayType::size() -
2413  1] == cell_level_index[(cell_batch_index + 1) *
2414  VectorizedArrayType::size() -
2415  2];
2416 }
2417 
2418 
2419 
2420 template <int dim, typename Number, typename VectorizedArrayType>
2421 unsigned int
2423 {
2424  return shape_info.size(2);
2425 }
2426 
2427 
2428 template <int dim, typename Number, typename VectorizedArrayType>
2429 unsigned int
2431  const std::pair<unsigned int, unsigned int> range) const
2432 {
2433  const auto &fe_indices = dof_info[0].cell_active_fe_index;
2434 
2435  if (fe_indices.empty() == true)
2436  return 0;
2437 
2438  const auto index = fe_indices[range.first];
2439 
2440  for (unsigned int i = range.first; i < range.second; ++i)
2441  AssertDimension(index, fe_indices[i]);
2442 
2443  return index;
2444 }
2445 
2446 
2447 
2448 template <int dim, typename Number, typename VectorizedArrayType>
2449 unsigned int
2451  const std::pair<unsigned int, unsigned int> range,
2452  const bool is_interior_face) const
2453 {
2454  const auto &fe_indices = dof_info[0].cell_active_fe_index;
2455 
2456  if (fe_indices.empty() == true)
2457  return 0;
2458 
2459  if (is_interior_face)
2460  {
2461  const unsigned int index =
2462  fe_indices[face_info.faces[range.first].cells_interior[0] /
2463  VectorizedArrayType::size()];
2464 
2465  for (unsigned int i = range.first; i < range.second; ++i)
2466  AssertDimension(index,
2467  fe_indices[face_info.faces[i].cells_interior[0] /
2468  VectorizedArrayType::size()]);
2469 
2470  return index;
2471  }
2472  else
2473  {
2474  const unsigned int index =
2475  fe_indices[face_info.faces[range.first].cells_exterior[0] /
2476  VectorizedArrayType::size()];
2477 
2478  for (unsigned int i = range.first; i < range.second; ++i)
2479  AssertDimension(index,
2480  fe_indices[face_info.faces[i].cells_exterior[0] /
2481  VectorizedArrayType::size()]);
2482 
2483  return index;
2484  }
2485 }
2486 
2487 
2488 
2489 template <int dim, typename Number, typename VectorizedArrayType>
2490 inline unsigned int
2492  const unsigned int cell_batch_index) const
2493 {
2494  Assert(!dof_info.empty(), ExcNotInitialized());
2495  AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2496  const std::vector<unsigned char> &n_lanes_filled =
2497  dof_info[0].n_vectorization_lanes_filled
2499  AssertIndexRange(cell_batch_index, n_lanes_filled.size());
2500 
2501  return n_lanes_filled[cell_batch_index];
2502 }
2503 
2504 
2505 
2506 template <int dim, typename Number, typename VectorizedArrayType>
2507 inline unsigned int
2509  const unsigned int face_batch_index) const
2510 {
2511  AssertIndexRange(face_batch_index, face_info.faces.size());
2512  Assert(!dof_info.empty(), ExcNotInitialized());
2513  const std::vector<unsigned char> &n_lanes_filled =
2514  dof_info[0].n_vectorization_lanes_filled
2516  AssertIndexRange(face_batch_index, n_lanes_filled.size());
2517  return n_lanes_filled[face_batch_index];
2518 }
2519 
2520 
2521 
2522 template <int dim, typename Number, typename VectorizedArrayType>
2523 inline unsigned int
2525  const unsigned int dof_handler_index,
2526  const unsigned int active_fe_index) const
2527 {
2528  return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2529 }
2530 
2531 
2532 
2533 template <int dim, typename Number, typename VectorizedArrayType>
2534 inline unsigned int
2536  const unsigned int quad_index,
2537  const unsigned int active_fe_index) const
2538 {
2539  AssertIndexRange(quad_index, mapping_info.cell_data.size());
2540  return mapping_info.cell_data[quad_index]
2541  .descriptor[active_fe_index]
2542  .n_q_points;
2543 }
2544 
2545 
2546 
2547 template <int dim, typename Number, typename VectorizedArrayType>
2548 inline unsigned int
2550  const unsigned int dof_handler_index,
2551  const unsigned int active_fe_index) const
2552 {
2553  return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2554 }
2555 
2556 
2557 
2558 template <int dim, typename Number, typename VectorizedArrayType>
2559 inline unsigned int
2561  const unsigned int quad_index,
2562  const unsigned int active_fe_index) const
2563 {
2564  AssertIndexRange(quad_index, mapping_info.face_data.size());
2565  return mapping_info.face_data[quad_index]
2566  .descriptor[active_fe_index]
2567  .n_q_points;
2568 }
2569 
2570 
2571 
2572 template <int dim, typename Number, typename VectorizedArrayType>
2573 inline const IndexSet &
2575  const unsigned int dof_handler_index) const
2576 {
2577  return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2578 }
2579 
2580 
2581 
2582 template <int dim, typename Number, typename VectorizedArrayType>
2583 inline const IndexSet &
2585  const unsigned int dof_handler_index) const
2586 {
2587  return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2588 }
2589 
2590 
2591 
2592 template <int dim, typename Number, typename VectorizedArrayType>
2595  const unsigned int dof_handler_index,
2596  const unsigned int index_quad,
2597  const unsigned int index_fe,
2598  const unsigned int active_fe_index,
2599  const unsigned int active_quad_index) const
2600 {
2601  AssertIndexRange(dof_handler_index, dof_info.size());
2602  const unsigned int ind =
2603  dof_info[dof_handler_index].global_base_element_offset + index_fe;
2604  AssertIndexRange(ind, shape_info.size(0));
2605  AssertIndexRange(index_quad, shape_info.size(1));
2606  AssertIndexRange(active_fe_index, shape_info.size(2));
2607  AssertIndexRange(active_quad_index, shape_info.size(3));
2608  return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2609 }
2610 
2611 
2612 
2613 template <int dim, typename Number, typename VectorizedArrayType>
2615  VectorizedArrayType::size()> &
2617  const unsigned int face_batch_index) const
2618 {
2619  AssertIndexRange(face_batch_index, face_info.faces.size());
2620  return face_info.faces[face_batch_index];
2621 }
2622 
2623 
2624 
2625 template <int dim, typename Number, typename VectorizedArrayType>
2626 inline const Table<3, unsigned int> &
2628  const
2629 {
2631 }
2632 
2633 
2634 
2635 template <int dim, typename Number, typename VectorizedArrayType>
2636 inline const Quadrature<dim> &
2638  const unsigned int quad_index,
2639  const unsigned int active_fe_index) const
2640 {
2641  AssertIndexRange(quad_index, mapping_info.cell_data.size());
2642  return mapping_info.cell_data[quad_index]
2643  .descriptor[active_fe_index]
2644  .quadrature;
2645 }
2646 
2647 
2648 
2649 template <int dim, typename Number, typename VectorizedArrayType>
2650 inline const Quadrature<dim - 1> &
2652  const unsigned int quad_index,
2653  const unsigned int active_fe_index) const
2654 {
2655  AssertIndexRange(quad_index, mapping_info.face_data.size());
2656  return mapping_info.face_data[quad_index]
2657  .descriptor[active_fe_index]
2658  .quadrature;
2659 }
2660 
2661 
2662 
2663 template <int dim, typename Number, typename VectorizedArrayType>
2664 inline unsigned int
2666  const std::pair<unsigned int, unsigned int> range) const
2667 {
2668  auto result = get_cell_category(range.first);
2669 
2670  for (unsigned int i = range.first; i < range.second; ++i)
2671  result = std::max(result, get_cell_category(i));
2672 
2673  return result;
2674 }
2675 
2676 
2677 
2678 template <int dim, typename Number, typename VectorizedArrayType>
2679 inline std::pair<unsigned int, unsigned int>
2681  const std::pair<unsigned int, unsigned int> range) const
2682 {
2683  auto result = get_face_category(range.first);
2684 
2685  for (unsigned int i = range.first; i < range.second; ++i)
2686  {
2687  result.first = std::max(result.first, get_face_category(i).first);
2688  result.second = std::max(result.second, get_face_category(i).second);
2689  }
2690 
2691  return result;
2692 }
2693 
2694 
2695 
2696 template <int dim, typename Number, typename VectorizedArrayType>
2697 inline unsigned int
2699  const unsigned int cell_batch_index) const
2700 {
2701  AssertIndexRange(0, dof_info.size());
2702  AssertIndexRange(cell_batch_index, dof_info[0].cell_active_fe_index.size());
2703  if (dof_info[0].cell_active_fe_index.empty())
2704  return 0;
2705  else
2706  return dof_info[0].cell_active_fe_index[cell_batch_index];
2707 }
2708 
2709 
2710 
2711 template <int dim, typename Number, typename VectorizedArrayType>
2712 inline std::pair<unsigned int, unsigned int>
2714  const unsigned int face_batch_index) const
2715 {
2716  AssertIndexRange(face_batch_index, face_info.faces.size());
2717  if (dof_info[0].cell_active_fe_index.empty())
2718  return std::make_pair(0U, 0U);
2719 
2720  std::pair<unsigned int, unsigned int> result = std::make_pair(0U, 0U);
2721  for (unsigned int v = 0;
2722  v < VectorizedArrayType::size() &&
2723  face_info.faces[face_batch_index].cells_interior[v] !=
2725  ++v)
2726  result.first = std::max(
2727  result.first,
2728  dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
2729  .cells_interior[v] /
2730  VectorizedArrayType::size()]);
2731  if (face_info.faces[face_batch_index].cells_exterior[0] !=
2733  for (unsigned int v = 0;
2734  v < VectorizedArrayType::size() &&
2735  face_info.faces[face_batch_index].cells_exterior[v] !=
2737  ++v)
2738  result.second = std::max(
2739  result.second,
2740  dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
2741  .cells_exterior[v] /
2742  VectorizedArrayType::size()]);
2743  else
2744  result.second = numbers::invalid_unsigned_int;
2745  return result;
2746 }
2747 
2748 
2749 
2750 template <int dim, typename Number, typename VectorizedArrayType>
2751 inline bool
2753 {
2754  return indices_are_initialized;
2755 }
2756 
2757 
2758 
2759 template <int dim, typename Number, typename VectorizedArrayType>
2760 inline bool
2762 {
2763  return mapping_is_initialized;
2764 }
2765 
2766 
2767 template <int dim, typename Number, typename VectorizedArrayType>
2768 inline unsigned int
2770 {
2771  return mg_level;
2772 }
2773 
2774 
2775 
2776 template <int dim, typename Number, typename VectorizedArrayType>
2779 {
2780  using list_type =
2781  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
2782  list_type &data = scratch_pad.get();
2783  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
2784  if (it->first == false)
2785  {
2786  it->first = true;
2787  return &it->second;
2788  }
2789  data.emplace_front(true, AlignedVector<VectorizedArrayType>());
2790  return &data.front().second;
2791 }
2792 
2793 
2794 
2795 template <int dim, typename Number, typename VectorizedArrayType>
2796 void
2798  const AlignedVector<VectorizedArrayType> *scratch) const
2799 {
2800  using list_type =
2801  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
2802  list_type &data = scratch_pad.get();
2803  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
2804  if (&it->second == scratch)
2805  {
2806  Assert(it->first == true, ExcInternalError());
2807  it->first = false;
2808  return;
2809  }
2810  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2811 }
2812 
2813 
2814 
2815 template <int dim, typename Number, typename VectorizedArrayType>
2819 {
2820  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
2822  it != scratch_pad_non_threadsafe.end();
2823  ++it)
2824  if (it->first == false)
2825  {
2826  it->first = true;
2827  return &it->second;
2828  }
2829  scratch_pad_non_threadsafe.push_front(
2830  std::make_pair(true, AlignedVector<Number>()));
2831  return &scratch_pad_non_threadsafe.front().second;
2832 }
2833 
2834 
2835 
2836 template <int dim, typename Number, typename VectorizedArrayType>
2837 void
2840  const AlignedVector<Number> *scratch) const
2841 {
2842  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
2844  it != scratch_pad_non_threadsafe.end();
2845  ++it)
2846  if (&it->second == scratch)
2847  {
2848  Assert(it->first == true, ExcInternalError());
2849  it->first = false;
2850  return;
2851  }
2852  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2853 }
2854 
2855 
2856 
2857 // ------------------------------ reinit functions ---------------------------
2858 
2859 namespace internal
2860 {
2861  namespace MatrixFreeImplementation
2862  {
2863  template <int dim, int spacedim>
2864  inline std::vector<IndexSet>
2865  extract_locally_owned_index_sets(
2866  const std::vector<const ::DoFHandler<dim, spacedim> *> &dofh,
2867  const unsigned int level)
2868  {
2869  std::vector<IndexSet> locally_owned_set;
2870  locally_owned_set.reserve(dofh.size());
2871  for (unsigned int j = 0; j < dofh.size(); ++j)
2873  locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
2874  else
2875  locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
2876  return locally_owned_set;
2877  }
2878  } // namespace MatrixFreeImplementation
2879 } // namespace internal
2880 
2881 
2882 
2883 template <int dim, typename Number, typename VectorizedArrayType>
2884 template <typename QuadratureType, typename number2, typename MappingType>
2885 void
2887  const MappingType & mapping,
2888  const DoFHandler<dim> & dof_handler,
2889  const AffineConstraints<number2> &constraints_in,
2890  const QuadratureType & quad,
2892  &additional_data)
2893 {
2894  std::vector<const DoFHandler<dim, dim> *> dof_handlers;
2895  std::vector<const AffineConstraints<number2> *> constraints;
2896 
2897  dof_handlers.push_back(&dof_handler);
2898  constraints.push_back(&constraints_in);
2899 
2900  std::vector<IndexSet> locally_owned_sets =
2901  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2902  dof_handlers, additional_data.mg_level);
2903 
2904  std::vector<hp::QCollection<dim>> quad_hp;
2905  quad_hp.emplace_back(quad);
2906 
2907  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
2908  dof_handlers,
2909  constraints,
2910  locally_owned_sets,
2911  quad_hp,
2912  additional_data);
2913 }
2914 
2915 
2916 
2917 template <int dim, typename Number, typename VectorizedArrayType>
2918 template <typename QuadratureType, typename number2, typename MappingType>
2919 void
2921  const MappingType & mapping,
2922  const std::vector<const DoFHandler<dim> *> & dof_handler,
2923  const std::vector<const AffineConstraints<number2> *> &constraint,
2924  const QuadratureType & quad,
2926  &additional_data)
2927 {
2928  std::vector<IndexSet> locally_owned_set =
2929  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2930  dof_handler, additional_data.mg_level);
2931  std::vector<hp::QCollection<dim>> quad_hp;
2932  quad_hp.emplace_back(quad);
2933 
2934  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
2935  dof_handler,
2936  constraint,
2937  locally_owned_set,
2938  quad_hp,
2939  additional_data);
2940 }
2941 
2942 
2943 
2944 template <int dim, typename Number, typename VectorizedArrayType>
2945 template <typename QuadratureType, typename number2, typename MappingType>
2946 void
2948  const MappingType & mapping,
2949  const std::vector<const DoFHandler<dim> *> & dof_handler,
2950  const std::vector<const AffineConstraints<number2> *> &constraint,
2951  const std::vector<QuadratureType> & quad,
2953  &additional_data)
2954 {
2955  std::vector<IndexSet> locally_owned_set =
2956  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2957  dof_handler, additional_data.mg_level);
2958  std::vector<hp::QCollection<dim>> quad_hp;
2959  for (unsigned int q = 0; q < quad.size(); ++q)
2960  quad_hp.emplace_back(quad[q]);
2961 
2962  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
2963  dof_handler,
2964  constraint,
2965  locally_owned_set,
2966  quad_hp,
2967  additional_data);
2968 }
2969 
2970 
2971 
2972 // ------------------------------ implementation of loops --------------------
2973 
2974 // internal helper functions that define how to call MPI data exchange
2975 // functions: for generic vectors, do nothing at all. For distributed vectors,
2976 // call update_ghost_values_start function and so on. If we have collections
2977 // of vectors, just do the individual functions of the components. In order to
2978 // keep ghost values consistent (whether we are in read or write mode), we
2979 // also reset the values at the end. the whole situation is a bit complicated
2980 // by the fact that we need to treat block vectors differently, which use some
2981 // additional helper functions to select the blocks and template magic.
2982 namespace internal
2983 {
2987  template <int dim, typename Number, typename VectorizedArrayType>
2988  struct VectorDataExchange
2989  {
2990  // A shift for the MPI messages to reduce the risk for accidental
2991  // interaction with other open communications that a user program might
2992  // set up (parallel vectors support unfinished communication). We let
2993  // the other vectors use the first 20 assigned numbers and start the
2994  // matrix-free communication.
2995  static constexpr unsigned int channel_shift = 20;
2996 
2997 
2998 
3003  VectorDataExchange(
3004  const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3005  const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3006  DataAccessOnFaces vector_face_access,
3007  const unsigned int n_components)
3008  : matrix_free(matrix_free)
3009  , vector_face_access(
3010  matrix_free.get_task_info().face_partition_data.empty() ?
3011  ::MatrixFree<dim, Number, VectorizedArrayType>::
3012  DataAccessOnFaces::unspecified :
3013  vector_face_access)
3014  , ghosts_were_set(false)
3015 # ifdef DEAL_II_WITH_MPI
3016  , tmp_data(n_components)
3017  , requests(n_components)
3018 # endif
3019  {
3020  (void)n_components;
3021  if (this->vector_face_access !=
3023  DataAccessOnFaces::unspecified)
3024  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3026  matrix_free.get_dof_info(c).vector_exchanger_face_variants.size(),
3027  5);
3028  }
3029 
3030 
3031 
3035  ~VectorDataExchange() // NOLINT
3036  {
3037 # ifdef DEAL_II_WITH_MPI
3038  for (unsigned int i = 0; i < tmp_data.size(); ++i)
3039  if (tmp_data[i] != nullptr)
3040  matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3041 # endif
3042  }
3043 
3044 
3045 
3050  template <typename VectorType>
3051  unsigned int
3052  find_vector_in_mf(const VectorType &vec,
3053  const bool check_global_compatibility = true) const
3054  {
3055  // case 1: vector was set up with MatrixFree::initialize_dof_vector()
3056  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3057  if (vec.get_partitioner().get() ==
3058  matrix_free.get_dof_info(c).vector_partitioner.get())
3059  return c;
3060 
3061  // case 2: user provided own partitioner (compatibility mode)
3062  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3063  if (check_global_compatibility ?
3064  vec.get_partitioner()->is_globally_compatible(
3065  *matrix_free.get_dof_info(c).vector_partitioner) :
3066  vec.get_partitioner()->is_compatible(
3067  *matrix_free.get_dof_info(c).vector_partitioner))
3068  return c;
3069 
3070  Assert(false,
3071  ExcNotImplemented("Could not find partitioner that fits vector"));
3072 
3074  }
3075 
3076 
3077 
3083  get_partitioner(const unsigned int mf_component) const
3084  {
3085  AssertDimension(matrix_free.get_dof_info(mf_component)
3086  .vector_exchanger_face_variants.size(),
3087  5);
3088  if (vector_face_access ==
3091  return *matrix_free.get_dof_info(mf_component)
3092  .vector_exchanger_face_variants[0];
3093  else if (vector_face_access ==
3096  return *matrix_free.get_dof_info(mf_component)
3097  .vector_exchanger_face_variants[1];
3098  else if (vector_face_access ==
3101  return *matrix_free.get_dof_info(mf_component)
3102  .vector_exchanger_face_variants[2];
3103  else if (vector_face_access ==
3105  DataAccessOnFaces::values_all_faces)
3106  return *matrix_free.get_dof_info(mf_component)
3107  .vector_exchanger_face_variants[3];
3108  else if (vector_face_access ==
3110  DataAccessOnFaces::gradients_all_faces)
3111  return *matrix_free.get_dof_info(mf_component)
3112  .vector_exchanger_face_variants[4];
3113  else
3114  return *matrix_free.get_dof_info(mf_component).vector_exchanger.get();
3115  }
3116 
3117 
3118 
3122  template <typename VectorType,
3123  typename std::enable_if<is_not_parallel_vector<VectorType>,
3124  VectorType>::type * = nullptr>
3125  void
3126  update_ghost_values_start(const unsigned int /*component_in_block_vector*/,
3127  const VectorType & /*vec*/)
3128  {}
3129 
3130 
3135  template <
3136  typename VectorType,
3137  typename std::enable_if<!has_update_ghost_values_start<VectorType> &&
3138  !is_not_parallel_vector<VectorType>,
3139  VectorType>::type * = nullptr>
3140  void
3141  update_ghost_values_start(const unsigned int component_in_block_vector,
3142  const VectorType & vec)
3143  {
3144  (void)component_in_block_vector;
3145  bool ghosts_set = vec.has_ghost_elements();
3146 
3147  Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3148  ghosts_set == false,
3149  ExcNotImplemented());
3150 
3151  if (ghosts_set)
3152  ghosts_were_set = true;
3153 
3154  vec.update_ghost_values();
3155  }
3156 
3157 
3158 
3164  template <
3165  typename VectorType,
3166  typename std::enable_if<has_update_ghost_values_start<VectorType> &&
3167  !has_exchange_on_subset<VectorType>,
3168  VectorType>::type * = nullptr>
3169  void
3170  update_ghost_values_start(const unsigned int component_in_block_vector,
3171  const VectorType & vec)
3172  {
3173  (void)component_in_block_vector;
3174  bool ghosts_set = vec.has_ghost_elements();
3175 
3176  Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3177  ghosts_set == false,
3178  ExcNotImplemented());
3179 
3180  if (ghosts_set)
3181  ghosts_were_set = true;
3182 
3183  vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3184  }
3185 
3186 
3187 
3194  template <
3195  typename VectorType,
3196  typename std::enable_if<has_update_ghost_values_start<VectorType> &&
3197  has_exchange_on_subset<VectorType>,
3198  VectorType>::type * = nullptr>
3199  void
3200  update_ghost_values_start(const unsigned int component_in_block_vector,
3201  const VectorType & vec)
3202  {
3203  static_assert(
3204  std::is_same<Number, typename VectorType::value_type>::value,
3205  "Type mismatch between VectorType and VectorDataExchange");
3206  (void)component_in_block_vector;
3207  bool ghosts_set = vec.has_ghost_elements();
3208 
3209  Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3210  ghosts_set == false,
3211  ExcNotImplemented());
3212 
3213  if (ghosts_set)
3214  ghosts_were_set = true;
3215 
3216  if (vec.size() != 0)
3217  {
3218 # ifdef DEAL_II_WITH_MPI
3219  const unsigned int mf_component = find_vector_in_mf(vec);
3220 
3221  const auto &part = get_partitioner(mf_component);
3222 
3223  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3224  part.n_import_sm_procs() == 0)
3225  return;
3226 
3227  tmp_data[component_in_block_vector] =
3228  matrix_free.acquire_scratch_data_non_threadsafe();
3229  tmp_data[component_in_block_vector]->resize_fast(
3230  part.n_import_indices());
3231  AssertDimension(requests.size(), tmp_data.size());
3232 
3233  part.export_to_ghosted_array_start(
3234  component_in_block_vector * 2 + channel_shift,
3235  ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3236  vec.shared_vector_data(),
3237  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3238  part.locally_owned_size(),
3239  matrix_free.get_dof_info(mf_component)
3240  .vector_partitioner->n_ghost_indices()),
3241  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3242  part.n_import_indices()),
3243  this->requests[component_in_block_vector]);
3244 # endif
3245  }
3246  }
3247 
3248 
3249 
3254  template <
3255  typename VectorType,
3256  typename std::enable_if<!has_update_ghost_values_start<VectorType>,
3257  VectorType>::type * = nullptr>
3258  void
3259  update_ghost_values_finish(const unsigned int /*component_in_block_vector*/,
3260  const VectorType & /*vec*/)
3261  {}
3262 
3263 
3264 
3270  template <
3271  typename VectorType,
3272  typename std::enable_if<has_update_ghost_values_start<VectorType> &&
3273  !has_exchange_on_subset<VectorType>,
3274  VectorType>::type * = nullptr>
3275  void
3276  update_ghost_values_finish(const unsigned int component_in_block_vector,
3277  const VectorType & vec)
3278  {
3279  (void)component_in_block_vector;
3280  vec.update_ghost_values_finish();
3281  }
3282 
3283 
3284 
3291  template <
3292  typename VectorType,
3293  typename std::enable_if<has_update_ghost_values_start<VectorType> &&
3294  has_exchange_on_subset<VectorType>,
3295  VectorType>::type * = nullptr>
3296  void
3297  update_ghost_values_finish(const unsigned int component_in_block_vector,
3298  const VectorType & vec)
3299  {
3300  static_assert(
3301  std::is_same<Number, typename VectorType::value_type>::value,
3302  "Type mismatch between VectorType and VectorDataExchange");
3303  (void)component_in_block_vector;
3304 
3305  if (vec.size() != 0)
3306  {
3307 # ifdef DEAL_II_WITH_MPI
3308  AssertIndexRange(component_in_block_vector, tmp_data.size());
3309  AssertDimension(requests.size(), tmp_data.size());
3310 
3311  const unsigned int mf_component = find_vector_in_mf(vec);
3312 
3313  const auto &part = get_partitioner(mf_component);
3314 
3315  if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3316  part.n_import_sm_procs() != 0)
3317  {
3318  part.export_to_ghosted_array_finish(
3319  ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3320  vec.shared_vector_data(),
3321  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3322  part.locally_owned_size(),
3323  matrix_free.get_dof_info(mf_component)
3324  .vector_partitioner->n_ghost_indices()),
3325  this->requests[component_in_block_vector]);
3326 
3327  matrix_free.release_scratch_data_non_threadsafe(
3328  tmp_data[component_in_block_vector]);
3329  tmp_data[component_in_block_vector] = nullptr;
3330  }
3331 # endif
3332  }
3333  // let vector know that ghosts are being updated and we can read from
3334  // them
3335  vec.set_ghost_state(true);
3336  }
3337 
3338 
3339 
3343  template <typename VectorType,
3344  typename std::enable_if<is_not_parallel_vector<VectorType>,
3345  VectorType>::type * = nullptr>
3346  void
3347  compress_start(const unsigned int /*component_in_block_vector*/,
3348  VectorType & /*vec*/)
3349  {}
3350 
3351 
3352 
3357  template <typename VectorType,
3358  typename std::enable_if<!has_compress_start<VectorType> &&
3359  !is_not_parallel_vector<VectorType>,
3360  VectorType>::type * = nullptr>
3361  void
3362  compress_start(const unsigned int component_in_block_vector,
3363  VectorType & vec)
3364  {
3365  (void)component_in_block_vector;
3366  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3367  vec.compress(::VectorOperation::add);
3368  }
3369 
3370 
3371 
3377  template <typename VectorType,
3378  typename std::enable_if<has_compress_start<VectorType> &&
3379  !has_exchange_on_subset<VectorType>,
3380  VectorType>::type * = nullptr>
3381  void
3382  compress_start(const unsigned int component_in_block_vector,
3383  VectorType & vec)
3384  {
3385  (void)component_in_block_vector;
3386  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3387  vec.compress_start(component_in_block_vector + channel_shift);
3388  }
3389 
3390 
3391 
3398  template <typename VectorType,
3399  typename std::enable_if<has_compress_start<VectorType> &&
3400  has_exchange_on_subset<VectorType>,
3401  VectorType>::type * = nullptr>
3402  void
3403  compress_start(const unsigned int component_in_block_vector,
3404  VectorType & vec)
3405  {
3406  static_assert(
3407  std::is_same<Number, typename VectorType::value_type>::value,
3408  "Type mismatch between VectorType and VectorDataExchange");
3409  (void)component_in_block_vector;
3410  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3411 
3412  if (vec.size() != 0)
3413  {
3414 # ifdef DEAL_II_WITH_MPI
3415  const unsigned int mf_component = find_vector_in_mf(vec);
3416 
3417  const auto &part = get_partitioner(mf_component);
3418 
3419  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3420  part.n_import_sm_procs() == 0)
3421  return;
3422 
3423  tmp_data[component_in_block_vector] =
3424  matrix_free.acquire_scratch_data_non_threadsafe();
3425  tmp_data[component_in_block_vector]->resize_fast(
3426  part.n_import_indices());
3427  AssertDimension(requests.size(), tmp_data.size());
3428 
3429  part.import_from_ghosted_array_start(
3431  component_in_block_vector * 2 + channel_shift,
3432  ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3433  vec.shared_vector_data(),
3434  ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3435  matrix_free.get_dof_info(mf_component)
3436  .vector_partitioner->n_ghost_indices()),
3437  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3438  part.n_import_indices()),
3439  this->requests[component_in_block_vector]);
3440 # endif
3441  }
3442  }
3443 
3444 
3445 
3450  template <typename VectorType,
3451  typename std::enable_if<!has_compress_start<VectorType>,
3452  VectorType>::type * = nullptr>
3453  void
3454  compress_finish(const unsigned int /*component_in_block_vector*/,
3455  VectorType & /*vec*/)
3456  {}
3457 
3458 
3459 
3465  template <typename VectorType,
3466  typename std::enable_if<has_compress_start<VectorType> &&
3467  !has_exchange_on_subset<VectorType>,
3468  VectorType>::type * = nullptr>
3469  void
3470  compress_finish(const unsigned int component_in_block_vector,
3471  VectorType & vec)
3472  {
3473  (void)component_in_block_vector;
3474  vec.compress_finish(::VectorOperation::add);
3475  }
3476 
3477 
3478 
3485  template <typename VectorType,
3486  typename std::enable_if<has_compress_start<VectorType> &&
3487  has_exchange_on_subset<VectorType>,
3488  VectorType>::type * = nullptr>
3489  void
3490  compress_finish(const unsigned int component_in_block_vector,
3491  VectorType & vec)
3492  {
3493  static_assert(
3494  std::is_same<Number, typename VectorType::value_type>::value,
3495  "Type mismatch between VectorType and VectorDataExchange");
3496  (void)component_in_block_vector;
3497  if (vec.size() != 0)
3498  {
3499 # ifdef DEAL_II_WITH_MPI
3500  AssertIndexRange(component_in_block_vector, tmp_data.size());
3501  AssertDimension(requests.size(), tmp_data.size());
3502 
3503  const unsigned int mf_component = find_vector_in_mf(vec);
3504 
3505  const auto &part = get_partitioner(mf_component);
3506 
3507  if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3508  part.n_import_sm_procs() != 0)
3509  {
3510  part.import_from_ghosted_array_finish(
3512  ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3513  vec.shared_vector_data(),
3514  ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3515  matrix_free.get_dof_info(mf_component)
3516  .vector_partitioner->n_ghost_indices()),
3518  tmp_data[component_in_block_vector]->begin(),
3519  part.n_import_indices()),
3520  this->requests[component_in_block_vector]);
3521 
3522  matrix_free.release_scratch_data_non_threadsafe(
3523  tmp_data[component_in_block_vector]);
3524  tmp_data[component_in_block_vector] = nullptr;
3525  }
3526 
3528  {
3529  const int ierr =
3530  MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3531  AssertThrowMPI(ierr);
3532  }
3533 # endif
3534  }
3535  }
3536 
3537 
3538 
3542  template <typename VectorType,
3543  typename std::enable_if<is_not_parallel_vector<VectorType>,
3544  VectorType>::type * = nullptr>
3545  void
3546  reset_ghost_values(const VectorType & /*vec*/) const
3547  {}
3548 
3549 
3550 
3555  template <typename VectorType,
3556  typename std::enable_if<!has_exchange_on_subset<VectorType> &&
3557  !is_not_parallel_vector<VectorType>,
3558  VectorType>::type * = nullptr>
3559  void
3560  reset_ghost_values(const VectorType &vec) const
3561  {
3562  if (ghosts_were_set == true)
3563  return;
3564 
3565  vec.zero_out_ghost_values();
3566  }
3567 
3568 
3569 
3575  template <typename VectorType,
3576  typename std::enable_if<has_exchange_on_subset<VectorType>,
3577  VectorType>::type * = nullptr>
3578  void
3579  reset_ghost_values(const VectorType &vec) const
3580  {
3581  static_assert(
3582  std::is_same<Number, typename VectorType::value_type>::value,
3583  "Type mismatch between VectorType and VectorDataExchange");
3584  if (ghosts_were_set == true)
3585  return;
3586 
3587  if (vec.size() != 0)
3588  {
3589 # ifdef DEAL_II_WITH_MPI
3590  AssertDimension(requests.size(), tmp_data.size());
3591 
3592  const unsigned int mf_component = find_vector_in_mf(vec);
3593 
3594  const auto &part = get_partitioner(mf_component);
3595 
3596  if (part.n_ghost_indices() > 0)
3597  {
3598  part.reset_ghost_values(ArrayView<Number>(
3600  .begin() +
3601  part.locally_owned_size(),
3602  matrix_free.get_dof_info(mf_component)
3603  .vector_partitioner->n_ghost_indices()));
3604  }
3605 
3606 # endif
3607  }
3608  // let vector know that it's not ghosted anymore
3609  vec.set_ghost_state(false);
3610  }
3611 
3612 
3613 
3619  template <typename VectorType,
3620  typename std::enable_if<has_exchange_on_subset<VectorType>,
3621  VectorType>::type * = nullptr>
3622  void
3623  zero_vector_region(const unsigned int range_index, VectorType &vec) const
3624  {
3625  static_assert(
3626  std::is_same<Number, typename VectorType::value_type>::value,
3627  "Type mismatch between VectorType and VectorDataExchange");
3628  if (range_index == numbers::invalid_unsigned_int)
3629  vec = Number();
3630  else
3631  {
3632  const unsigned int mf_component = find_vector_in_mf(vec, false);
3633  const internal::MatrixFreeFunctions::DoFInfo &dof_info =
3634  matrix_free.get_dof_info(mf_component);
3635  Assert(dof_info.vector_zero_range_list_index.empty() == false,
3636  ExcNotInitialized());
3637 
3638  Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3639  ExcInternalError());
3640  AssertIndexRange(range_index,
3641  dof_info.vector_zero_range_list_index.size() - 1);
3642  for (unsigned int id =
3643  dof_info.vector_zero_range_list_index[range_index];
3644  id != dof_info.vector_zero_range_list_index[range_index + 1];
3645  ++id)
3646  std::memset(vec.begin() + dof_info.vector_zero_range_list[id].first,
3647  0,
3648  (dof_info.vector_zero_range_list[id].second -
3649  dof_info.vector_zero_range_list[id].first) *
3650  sizeof(Number));
3651  }
3652  }
3653 
3654 
3655 
3661  template <typename VectorType,
3662  typename std::enable_if<!has_exchange_on_subset<VectorType>,
3663  VectorType>::type * = nullptr,
3664  typename VectorType::value_type * = nullptr>
3665  void
3666  zero_vector_region(const unsigned int range_index, VectorType &vec) const
3667  {
3668  if (range_index == numbers::invalid_unsigned_int || range_index == 0)
3669  vec = typename VectorType::value_type();
3670  }
3671 
3672 
3673 
3678  void
3679  zero_vector_region(const unsigned int, ...) const
3680  {
3681  Assert(false,
3682  ExcNotImplemented("Zeroing is only implemented for vector types "
3683  "which provide VectorType::value_type"));
3684  }
3685 
3686 
3687 
3688  const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
3689  const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3690  DataAccessOnFaces vector_face_access;
3691  bool ghosts_were_set;
3692 # ifdef DEAL_II_WITH_MPI
3693  std::vector<AlignedVector<Number> *> tmp_data;
3694  std::vector<std::vector<MPI_Request>> requests;
3695 # endif
3696  }; // VectorDataExchange
3697 
3698  template <typename VectorStruct>
3699  unsigned int
3700  n_components(const VectorStruct &vec);
3701 
3702  template <typename VectorStruct>
3703  unsigned int
3704  n_components_block(const VectorStruct &vec,
3705  std::integral_constant<bool, true>)
3706  {
3707  unsigned int components = 0;
3708  for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
3709  components += n_components(vec.block(bl));
3710  return components;
3711  }
3712 
3713  template <typename VectorStruct>
3714  unsigned int
3715  n_components_block(const VectorStruct &, std::integral_constant<bool, false>)
3716  {
3717  return 1;
3718  }
3719 
3720  template <typename VectorStruct>
3721  unsigned int
3722  n_components(const VectorStruct &vec)
3723  {
3724  return n_components_block(
3725  vec, std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3726  }
3727 
3728  template <typename VectorStruct>
3729  inline unsigned int
3730  n_components(const std::vector<VectorStruct> &vec)
3731  {
3732  unsigned int components = 0;
3733  for (unsigned int comp = 0; comp < vec.size(); ++comp)
3734  components += n_components_block(
3735  vec[comp],
3736  std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3737  return components;
3738  }
3739 
3740  template <typename VectorStruct>
3741  inline unsigned int
3742  n_components(const std::vector<VectorStruct *> &vec)
3743  {
3744  unsigned int components = 0;
3745  for (unsigned int comp = 0; comp < vec.size(); ++comp)
3746  components += n_components_block(
3747  *vec[comp],
3748  std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3749  return components;
3750  }
3751 
3752 
3753 
3754  // A helper function to identify block vectors with many components where we
3755  // should not try to overlap computations and communication because there
3756  // would be too many outstanding communication requests.
3757 
3758  // default value for vectors that do not have communication_block_size
3759  template <typename VectorStruct,
3760  typename std::enable_if<!has_communication_block_size<VectorStruct>,
3761  VectorStruct>::type * = nullptr>
3762  constexpr unsigned int
3763  get_communication_block_size(const VectorStruct &)
3764  {
3766  }
3767 
3768 
3769 
3770  template <typename VectorStruct,
3771  typename std::enable_if<has_communication_block_size<VectorStruct>,
3772  VectorStruct>::type * = nullptr>
3773  constexpr unsigned int
3774  get_communication_block_size(const VectorStruct &)
3775  {
3776  return VectorStruct::communication_block_size;
3777  }
3778 
3779 
3780 
3781  // --------------------------------------------------------------------------
3782  // below we have wrappers to distinguish between block and non-block vectors.
3783  // --------------------------------------------------------------------------
3784 
3785  //
3786  // update_ghost_values_start
3787  //
3788 
3789  // update_ghost_values for block vectors
3790  template <int dim,
3791  typename VectorStruct,
3792  typename Number,
3793  typename VectorizedArrayType,
3794  typename std::enable_if<IsBlockVector<VectorStruct>::value,
3795  VectorStruct>::type * = nullptr>
3796  void
3797  update_ghost_values_start(
3798  const VectorStruct & vec,
3799  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3800  const unsigned int channel = 0)
3801  {
3802  if (get_communication_block_size(vec) < vec.n_blocks())
3803  {
3804  // don't forget to set ghosts_were_set, that otherwise happens
3805  // inside VectorDataExchange::update_ghost_values_start()
3806  exchanger.ghosts_were_set = vec.has_ghost_elements();
3807  vec.update_ghost_values();
3808  }
3809  else
3810  {
3811  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
3812  update_ghost_values_start(vec.block(i), exchanger, channel + i);
3813  }
3814  }
3815 
3816 
3817 
3818  // update_ghost_values for non-block vectors
3819  template <int dim,
3820  typename VectorStruct,
3821  typename Number,
3822  typename VectorizedArrayType,
3823  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
3824  VectorStruct>::type * = nullptr>
3825  void
3826  update_ghost_values_start(
3827  const VectorStruct & vec,
3828  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3829  const unsigned int channel = 0)
3830  {
3831  exchanger.update_ghost_values_start(channel, vec);
3832  }
3833 
3834 
3835 
3836  // update_ghost_values_start() for vector of vectors
3837  template <int dim,
3838  typename VectorStruct,
3839  typename Number,
3840  typename VectorizedArrayType>
3841  inline void
3842  update_ghost_values_start(
3843  const std::vector<VectorStruct> & vec,
3844  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3845  {
3846  unsigned int component_index = 0;
3847  for (unsigned int comp = 0; comp < vec.size(); ++comp)
3848  {
3849  update_ghost_values_start(vec[comp], exchanger, component_index);
3850  component_index += n_components(vec[comp]);
3851  }
3852  }
3853 
3854 
3855 
3856  // update_ghost_values_start() for vector of pointers to vectors
3857  template <int dim,
3858  typename VectorStruct,
3859  typename Number,
3860  typename VectorizedArrayType>
3861  inline void
3862  update_ghost_values_start(
3863  const std::vector<VectorStruct *> & vec,
3864  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3865  {
3866  unsigned int component_index = 0;
3867  for (unsigned int comp = 0; comp < vec.size(); ++comp)
3868  {
3869  update_ghost_values_start(*vec[comp], exchanger, component_index);
3870  component_index += n_components(*vec[comp]);
3871  }
3872  }
3873 
3874 
3875 
3876  //
3877  // update_ghost_values_finish
3878  //
3879 
3880  // for block vectors
3881  template <int dim,
3882  typename VectorStruct,
3883  typename Number,
3884  typename VectorizedArrayType,
3885  typename std::enable_if<IsBlockVector<VectorStruct>::value,
3886  VectorStruct>::type * = nullptr>
3887  void
3888  update_ghost_values_finish(
3889  const VectorStruct & vec,
3890  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3891  const unsigned int channel = 0)
3892  {
3893  if (get_communication_block_size(vec) < vec.n_blocks())
3894  {
3895  // do nothing, everything has already been completed in the _start()
3896  // call
3897  }
3898  else
3899  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
3900  update_ghost_values_finish(vec.block(i), exchanger, channel + i);
3901  }
3902 
3903 
3904 
3905  // for non-block vectors
3906  template <int dim,
3907  typename VectorStruct,
3908  typename Number,
3909  typename VectorizedArrayType,
3910  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
3911  VectorStruct>::type * = nullptr>
3912  void
3913  update_ghost_values_finish(
3914  const VectorStruct & vec,
3915  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3916  const unsigned int channel = 0)
3917  {
3918  exchanger.update_ghost_values_finish(channel, vec);
3919  }
3920 
3921 
3922 
3923  // for vector of vectors
3924  template <int dim,
3925  typename VectorStruct,
3926  typename Number,
3927  typename VectorizedArrayType>
3928  inline void
3929  update_ghost_values_finish(
3930  const std::vector<VectorStruct> & vec,
3931  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3932  {
3933  unsigned int component_index = 0;
3934  for (unsigned int comp = 0; comp < vec.size(); ++comp)
3935  {
3936  update_ghost_values_finish(vec[comp], exchanger, component_index);
3937  component_index += n_components(vec[comp]);
3938  }
3939  }
3940 
3941 
3942 
3943  // for vector of pointers to vectors
3944  template <int dim,
3945  typename VectorStruct,
3946  typename Number,
3947  typename VectorizedArrayType>
3948  inline void
3949  update_ghost_values_finish(
3950  const std::vector<VectorStruct *> & vec,
3951  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3952  {
3953  unsigned int component_index = 0;
3954  for (unsigned int comp = 0; comp < vec.size(); ++comp)
3955  {
3956  update_ghost_values_finish(*vec[comp], exchanger, component_index);
3957  component_index += n_components(*vec[comp]);
3958  }
3959  }
3960 
3961 
3962 
3963  //
3964  // compress_start
3965  //
3966 
3967  // for block vectors
3968  template <int dim,
3969  typename VectorStruct,
3970  typename Number,
3971  typename VectorizedArrayType,
3972  typename std::enable_if<IsBlockVector<VectorStruct>::value,
3973  VectorStruct>::type * = nullptr>
3974  inline void
3975  compress_start(
3976  VectorStruct & vec,
3977  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3978  const unsigned int channel = 0)
3979  {
3980  if (get_communication_block_size(vec) < vec.n_blocks())
3981  vec.compress(::VectorOperation::add);
3982  else
3983  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
3984  compress_start(vec.block(i), exchanger, channel + i);
3985  }
3986 
3987 
3988 
3989  // for non-block vectors
3990  template <int dim,
3991  typename VectorStruct,
3992  typename Number,
3993  typename VectorizedArrayType,
3994  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
3995  VectorStruct>::type * = nullptr>
3996  inline void
3997  compress_start(
3998  VectorStruct & vec,
3999  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4000  const unsigned int channel = 0)
4001  {
4002  exchanger.compress_start(channel, vec);
4003  }
4004 
4005 
4006 
4007  // for std::vector of vectors
4008  template <int dim,
4009  typename VectorStruct,
4010  typename Number,
4011  typename VectorizedArrayType>
4012  inline void
4013  compress_start(
4014  std::vector<VectorStruct> & vec,
4015  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4016  {
4017  unsigned int component_index = 0;
4018  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4019  {
4020  compress_start(vec[comp], exchanger, component_index);
4021  component_index += n_components(vec[comp]);
4022  }
4023  }
4024 
4025 
4026 
4027  // for std::vector of pointer to vectors
4028  template <int dim,
4029  typename VectorStruct,
4030  typename Number,
4031  typename VectorizedArrayType>
4032  inline void
4033  compress_start(
4034  std::vector<VectorStruct *> & vec,
4035  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4036  {
4037  unsigned int component_index = 0;
4038  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4039  {
4040  compress_start(*vec[comp], exchanger, component_index);
4041  component_index += n_components(*vec[comp]);
4042  }
4043  }
4044 
4045 
4046 
4047  //
4048  // compress_finish
4049  //
4050 
4051  // for block vectors
4052  template <int dim,
4053  typename VectorStruct,
4054  typename Number,
4055  typename VectorizedArrayType,
4056  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4057  VectorStruct>::type * = nullptr>
4058  inline void
4059  compress_finish(
4060  VectorStruct & vec,
4061  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4062  const unsigned int channel = 0)
4063  {
4064  if (get_communication_block_size(vec) < vec.n_blocks())
4065  {
4066  // do nothing, everything has already been completed in the _start()
4067  // call
4068  }
4069  else
4070  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4071  compress_finish(vec.block(i), exchanger, channel + i);
4072  }
4073 
4074 
4075 
4076  // for non-block vectors
4077  template <int dim,
4078  typename VectorStruct,
4079  typename Number,
4080  typename VectorizedArrayType,
4081  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4082  VectorStruct>::type * = nullptr>
4083  inline void
4084  compress_finish(
4085  VectorStruct & vec,
4086  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4087  const unsigned int channel = 0)
4088  {
4089  exchanger.compress_finish(channel, vec);
4090  }
4091 
4092 
4093 
4094  // for std::vector of vectors
4095  template <int dim,
4096  typename VectorStruct,
4097  typename Number,
4098  typename VectorizedArrayType>
4099  inline void
4100  compress_finish(
4101  std::vector<VectorStruct> & vec,
4102  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4103  {
4104  unsigned int component_index = 0;
4105  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4106  {
4107  compress_finish(vec[comp], exchanger, component_index);
4108  component_index += n_components(vec[comp]);
4109  }
4110  }
4111 
4112 
4113 
4114  // for std::vector of pointer to vectors
4115  template <int dim,
4116  typename VectorStruct,
4117  typename Number,
4118  typename VectorizedArrayType>
4119  inline void
4120  compress_finish(
4121  std::vector<VectorStruct *> & vec,
4122  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4123  {
4124  unsigned int component_index = 0;
4125  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4126  {
4127  compress_finish(*vec[comp], exchanger, component_index);
4128  component_index += n_components(*vec[comp]);
4129  }
4130  }
4131 
4132 
4133 
4134  //
4135  // reset_ghost_values:
4136  //
4137  // if the input vector did not have ghosts imported, clear them here again
4138  // in order to avoid subsequent operations e.g. in linear solvers to work
4139  // with ghosts all the time
4140  //
4141 
4142  // for block vectors
4143  template <int dim,
4144  typename VectorStruct,
4145  typename Number,
4146  typename VectorizedArrayType,
4147  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4148  VectorStruct>::type * = nullptr>
4149  inline void
4150  reset_ghost_values(
4151  const VectorStruct & vec,
4152  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4153  {
4154  // return immediately if there is nothing to do.
4155  if (exchanger.ghosts_were_set == true)
4156  return;
4157 
4158  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4159  reset_ghost_values(vec.block(i), exchanger);
4160  }
4161 
4162 
4163 
4164  // for non-block vectors
4165  template <int dim,
4166  typename VectorStruct,
4167  typename Number,
4168  typename VectorizedArrayType,
4169  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4170  VectorStruct>::type * = nullptr>
4171  inline void
4172  reset_ghost_values(
4173  const VectorStruct & vec,
4174  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4175  {
4176  exchanger.reset_ghost_values(vec);
4177  }
4178 
4179 
4180 
4181  // for std::vector of vectors
4182  template <int dim,
4183  typename VectorStruct,
4184  typename Number,
4185  typename VectorizedArrayType>
4186  inline void
4187  reset_ghost_values(
4188  const std::vector<VectorStruct> & vec,
4189  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4190  {
4191  // return immediately if there is nothing to do.
4192  if (exchanger.ghosts_were_set == true)
4193  return;
4194 
4195  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4196  reset_ghost_values(vec[comp], exchanger);
4197  }
4198 
4199 
4200 
4201  // for std::vector of pointer to vectors
4202  template <int dim,
4203  typename VectorStruct,
4204  typename Number,
4205  typename VectorizedArrayType>
4206  inline void
4207  reset_ghost_values(
4208  const std::vector<VectorStruct *> & vec,
4209  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4210  {
4211  // return immediately if there is nothing to do.
4212  if (exchanger.ghosts_were_set == true)
4213  return;
4214 
4215  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4216  reset_ghost_values(*vec[comp], exchanger);
4217  }
4218 
4219 
4220 
4221  //
4222  // zero_vector_region
4223  //
4224 
4225  // for block vectors
4226  template <int dim,
4227  typename VectorStruct,
4228  typename Number,
4229  typename VectorizedArrayType,
4230  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4231  VectorStruct>::type * = nullptr>
4232  inline void
4233  zero_vector_region(
4234  const unsigned int range_index,
4235  VectorStruct & vec,
4236  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4237  {
4238  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4239  exchanger.zero_vector_region(range_index, vec.block(i));
4240  }
4241 
4242 
4243 
4244  // for non-block vectors
4245  template <int dim,
4246  typename VectorStruct,
4247  typename Number,
4248  typename VectorizedArrayType,
4249  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4250  VectorStruct>::type * = nullptr>
4251  inline void
4252  zero_vector_region(
4253  const unsigned int range_index,
4254  VectorStruct & vec,
4255  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4256  {
4257  exchanger.zero_vector_region(range_index, vec);
4258  }
4259 
4260 
4261 
4262  // for std::vector of vectors
4263  template <int dim,
4264  typename VectorStruct,
4265  typename Number,
4266  typename VectorizedArrayType>
4267  inline void
4268  zero_vector_region(
4269  const unsigned int range_index,
4270  std::vector<VectorStruct> & vec,
4271  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4272  {
4273  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4274  zero_vector_region(range_index, vec[comp], exchanger);
4275  }
4276 
4277 
4278 
4279  // for std::vector of pointers to vectors
4280  template <int dim,
4281  typename VectorStruct,
4282  typename Number,
4283  typename VectorizedArrayType>
4284  inline void
4285  zero_vector_region(
4286  const unsigned int range_index,
4287  std::vector<VectorStruct *> & vec,
4288  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4289  {
4290  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4291  zero_vector_region(range_index, *vec[comp], exchanger);
4292  }
4293 
4294 
4295 
4296  // Apply a unit matrix operation to constrained DoFs: Default cases where we
4297  // cannot detect a LinearAlgebra::distributed::Vector, we do not do
4298  // anything, else we apply the constraints as a unit operation
4299  template <typename VectorStruct1, typename VectorStruct2>
4300  inline void
4301  apply_operation_to_constrained_dofs(const std::vector<unsigned int> &,
4302  const VectorStruct1 &,
4303  VectorStruct2 &)
4304  {}
4305 
4306  template <typename Number>
4307  inline void
4308  apply_operation_to_constrained_dofs(
4309  const std::vector<unsigned int> & constrained_dofs,
4312  {
4313  for (const unsigned int i : constrained_dofs)
4314  dst.local_element(i) = src.local_element(i);
4315  }
4316 
4317 
4318  namespace MatrixFreeFunctions
4319  {
4320  // struct to select between a const interface and a non-const interface
4321  // for MFWorker
4322  template <typename, typename, typename, typename, bool>
4323  struct InterfaceSelector
4324  {};
4325 
4326  // Version for constant functions
4327  template <typename MF,
4328  typename InVector,
4329  typename OutVector,
4330  typename Container>
4331  struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4332  {
4333  using function_type = void (Container::*)(
4334  const MF &,
4335  OutVector &,
4336  const InVector &,
4337  const std::pair<unsigned int, unsigned int> &) const;
4338  };
4339 
4340  // Version for non-constant functions
4341  template <typename MF,
4342  typename InVector,
4343  typename OutVector,
4344  typename Container>
4345  struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4346  {
4347  using function_type =
4348  void (Container::*)(const MF &,
4349  OutVector &,
4350  const InVector &,
4351  const std::pair<unsigned int, unsigned int> &);
4352  };
4353  } // namespace MatrixFreeFunctions
4354 
4355 
4356 
4357  // A implementation class for the worker object that runs the various
4358  // operations we want to perform during the matrix-free loop
4359  template <typename MF,
4360  typename InVector,
4361  typename OutVector,
4362  typename Container,
4363  bool is_constant>
4364  class MFWorker : public MFWorkerInterface
4365  {
4366  public:
4367  // An alias to make the arguments further down more readable
4368  using function_type = typename MatrixFreeFunctions::
4369  InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4370  function_type;
4371 
4372  // constructor, binds all the arguments to this class
4373  MFWorker(const MF & matrix_free,
4374  const InVector & src,
4375  OutVector & dst,
4376  const bool zero_dst_vector_setting,
4377  const Container & container,
4378  function_type cell_function,
4379  function_type face_function,
4380  function_type boundary_function,
4381  const typename MF::DataAccessOnFaces src_vector_face_access =
4383  const typename MF::DataAccessOnFaces dst_vector_face_access =
4385  const std::function<void(const unsigned int, const unsigned int)>
4386  &operation_before_loop = {},
4387  const std::function<void(const unsigned int, const unsigned int)>
4388  & operation_after_loop = {},
4389  const unsigned int dof_handler_index_pre_post = 0)
4390  : matrix_free(matrix_free)
4391  , container(const_cast<Container &>(container))
4392  , cell_function(cell_function)
4393  , face_function(face_function)
4394  , boundary_function(boundary_function)
4395  , src(src)
4396  , dst(dst)
4397  , src_data_exchanger(matrix_free,
4398  src_vector_face_access,
4399  n_components(src))
4400  , dst_data_exchanger(matrix_free,
4401  dst_vector_face_access,
4402  n_components(dst))
4403  , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4404  , zero_dst_vector_setting(zero_dst_vector_setting &&
4405  !src_and_dst_are_same)
4406  , operation_before_loop(operation_before_loop)
4407  , operation_after_loop(operation_after_loop)
4408  , dof_handler_index_pre_post(dof_handler_index_pre_post)
4409  {}
4410 
4411  // Runs the cell work. If no function is given, nothing is done
4412  virtual void
4413  cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4414  {
4415  if (cell_function != nullptr && cell_range.second > cell_range.first)
4416  for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4417  {
4418  const auto cell_subrange =
4419  matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4420 
4421  if (cell_subrange.second <= cell_subrange.first)
4422  continue;
4423 
4424  (container.*
4425  cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4426  }
4427  }
4428 
4429  virtual void
4430  cell(const unsigned int range_index) override
4431  {
4432  process_range(cell_function,
4433  matrix_free.get_task_info().cell_partition_data_hp_ptr,
4434  matrix_free.get_task_info().cell_partition_data_hp,
4435  range_index);
4436  }
4437 
4438  virtual void
4439  face(const unsigned int range_index) override
4440  {
4441  process_range(face_function,
4442  matrix_free.get_task_info().face_partition_data_hp_ptr,
4443  matrix_free.get_task_info().face_partition_data_hp,
4444  range_index);
4445  }
4446 
4447  virtual void
4448  boundary(const unsigned int range_index) override
4449  {
4450  process_range(boundary_function,
4451  matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4452  matrix_free.get_task_info().boundary_partition_data_hp,
4453  range_index);
4454  }
4455 
4456  private:
4457  void
4458  process_range(const function_type & fu,
4459  const std::vector<unsigned int> &ptr,
4460  const std::vector<unsigned int> &data,
4461  const unsigned int range_index)
4462  {
4463  if (fu == nullptr)
4464  return;
4465 
4466  AssertIndexRange(range_index + 1, ptr.size());
4467  for (unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4468  {
4469  AssertIndexRange(2 * i + 1, data.size());
4470  (container.*fu)(matrix_free,
4471  this->dst,
4472  this->src,
4473  std::make_pair(data[2 * i], data[2 * i + 1]));
4474  }
4475  }
4476 
4477  public:
4478  // Starts the communication for the update ghost values operation. We
4479  // cannot call this update if ghost and destination are the same because
4480  // that would introduce spurious entries in the destination (there is also
4481  // the problem that reading from a vector that we also write to is usually
4482  // not intended in case there is overlap, but this is up to the
4483  // application code to decide and we cannot catch this case here).
4484  virtual void
4485  vector_update_ghosts_start() override
4486  {
4487  if (!src_and_dst_are_same)
4488  internal::update_ghost_values_start(src, src_data_exchanger);
4489  }
4490 
4491  // Finishes the communication for the update ghost values operation
4492  virtual void
4493  vector_update_ghosts_finish() override
4494  {
4495  if (!src_and_dst_are_same)
4496  internal::update_ghost_values_finish(src, src_data_exchanger);
4497  }
4498 
4499  // Starts the communication for the vector compress operation
4500  virtual void
4501  vector_compress_start() override
4502  {
4503  internal::compress_start(dst, dst_data_exchanger);
4504  }
4505 
4506  // Finishes the communication for the vector compress operation
4507  virtual void
4508  vector_compress_finish() override
4509  {
4510  internal::compress_finish(dst, dst_data_exchanger);
4511  if (!src_and_dst_are_same)
4512  internal::reset_ghost_values(src, src_data_exchanger);
4513  }
4514 
4515  // Zeros the given input vector
4516  virtual void
4517  zero_dst_vector_range(const unsigned int range_index) override
4518  {
4519  if (zero_dst_vector_setting)
4520  internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4521  }
4522 
4523  virtual void
4524  cell_loop_pre_range(const unsigned int range_index) override
4525  {
4526  if (operation_before_loop)
4527  {
4528  const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4529  matrix_free.get_dof_info(dof_handler_index_pre_post);
4530  if (range_index == numbers::invalid_unsigned_int)
4531  {
4532  // Case with threaded loop -> currently no overlap implemented
4534  0U,
4535  dof_info.vector_partitioner->locally_owned_size(),
4536  operation_before_loop,
4538  }
4539  else
4540  {
4541  AssertIndexRange(range_index,
4542  dof_info.cell_loop_pre_list_index.size() - 1);
4543  for (unsigned int id =
4544  dof_info.cell_loop_pre_list_index[range_index];
4545  id != dof_info.cell_loop_pre_list_index[range_index + 1];
4546  ++id)
4547  operation_before_loop(dof_info.cell_loop_pre_list[id].first,
4548  dof_info.cell_loop_pre_list[id].second);
4549  }
4550  }
4551  }
4552 
4553  virtual void
4554  cell_loop_post_range(const unsigned int range_index) override
4555  {
4556  if (operation_after_loop)
4557  {
4558  // Run unit matrix operation on constrained dofs if we are at the
4559  // last range
4560  const std::vector<unsigned int> &partition_row_index =
4561  matrix_free.get_task_info().partition_row_index;
4562  if (range_index ==
4563  partition_row_index[partition_row_index.size() - 2] - 1)
4564  apply_operation_to_constrained_dofs(
4565  matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
4566  src,
4567  dst);
4568 
4569  const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4570  matrix_free.get_dof_info(dof_handler_index_pre_post);
4571  if (range_index == numbers::invalid_unsigned_int)
4572  {
4573  // Case with threaded loop -> currently no overlap implemented
4575  0U,
4576  dof_info.vector_partitioner->locally_owned_size(),
4577  operation_after_loop,
4579  }
4580  else
4581  {
4582  AssertIndexRange(range_index,
4583  dof_info.cell_loop_post_list_index.size() - 1);
4584  for (unsigned int id =
4585  dof_info.cell_loop_post_list_index[range_index];
4586  id != dof_info.cell_loop_post_list_index[range_index + 1];
4587  ++id)
4588  operation_after_loop(dof_info.cell_loop_post_list[id].first,
4589  dof_info.cell_loop_post_list[id].second);
4590  }
4591  }
4592  }
4593 
4594  private:
4595  const MF & matrix_free;
4596  Container & container;
4597  function_type cell_function;
4598  function_type face_function;
4599  function_type boundary_function;
4600 
4601  const InVector &src;
4602  OutVector & dst;
4603  VectorDataExchange<MF::dimension,
4604  typename MF::value_type,
4605  typename MF::vectorized_value_type>
4606  src_data_exchanger;
4607  VectorDataExchange<MF::dimension,
4608  typename MF::value_type,
4609  typename MF::vectorized_value_type>
4610  dst_data_exchanger;
4611  const bool src_and_dst_are_same;
4612  const bool zero_dst_vector_setting;
4613  const std::function<void(const unsigned int, const unsigned int)>
4614  operation_before_loop;
4615  const std::function<void(const unsigned int, const unsigned int)>
4616  operation_after_loop;
4617  const unsigned int dof_handler_index_pre_post;
4618  };
4619 
4620 
4621 
4626  template <class MF, typename InVector, typename OutVector>
4627  struct MFClassWrapper
4628  {
4629  using function_type =
4630  std::function<void(const MF &,
4631  OutVector &,
4632  const InVector &,
4633  const std::pair<unsigned int, unsigned int> &)>;
4634 
4635  MFClassWrapper(const function_type cell,
4636  const function_type face,
4637  const function_type boundary)
4638  : cell(cell)
4639  , face(face)
4640  , boundary(boundary)
4641  {}
4642 
4643  void
4644  cell_integrator(const MF & mf,
4645  OutVector & dst,
4646  const InVector & src,
4647  const std::pair<unsigned int, unsigned int> &range) const
4648  {
4649  if (cell)
4650  cell(mf, dst, src, range);
4651  }
4652 
4653  void
4654  face_integrator(const MF & mf,
4655  OutVector & dst,
4656  const InVector & src,
4657  const std::pair<unsigned int, unsigned int> &range) const
4658  {
4659  if (face)
4660  face(mf, dst, src, range);
4661  }
4662 
4663  void
4664  boundary_integrator(
4665  const MF & mf,
4666  OutVector & dst,
4667  const InVector & src,
4668  const std::pair<unsigned int, unsigned int> &range) const
4669  {
4670  if (boundary)
4671  boundary(mf, dst, src, range);
4672  }
4673 
4674  const function_type cell;
4675  const function_type face;
4676  const function_type boundary;
4677  };
4678 
4679 } // end of namespace internal
4680 
4681 
4682 
4683 template <int dim, typename Number, typename VectorizedArrayType>
4684 template <typename OutVector, typename InVector>
4685 inline void
4687  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4688  OutVector &,
4689  const InVector &,
4690  const std::pair<unsigned int, unsigned int> &)>
4691  & cell_operation,
4692  OutVector & dst,
4693  const InVector &src,
4694  const bool zero_dst_vector) const
4695 {
4696  using Wrapper =
4697  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4698  InVector,
4699  OutVector>;
4700  Wrapper wrap(cell_operation, nullptr, nullptr);
4701  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4702  InVector,
4703  OutVector,
4704  Wrapper,
4705  true>
4706  worker(*this,
4707  src,
4708  dst,
4709  zero_dst_vector,
4710  wrap,
4711  &Wrapper::cell_integrator,
4712  &Wrapper::face_integrator,
4713  &Wrapper::boundary_integrator);
4714 
4715  task_info.loop(worker);
4716 }
4717 
4718 
4719 
4720 template <int dim, typename Number, typename VectorizedArrayType>
4721 template <typename OutVector, typename InVector>
4722 inline void
4724  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4725  OutVector &,
4726  const InVector &,
4727  const std::pair<unsigned int, unsigned int> &)>
4728  & cell_operation,
4729  OutVector & dst,
4730  const InVector &src,
4731  const std::function<void(const unsigned int, const unsigned int)>
4732  &operation_before_loop,
4733  const std::function<void(const unsigned int, const unsigned int)>
4734  & operation_after_loop,
4735  const unsigned int dof_handler_index_pre_post) const
4736 {
4737  using Wrapper =
4738  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4739  InVector,
4740  OutVector>;
4741  Wrapper wrap(cell_operation, nullptr, nullptr);
4742  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4743  InVector,
4744  OutVector,
4745  Wrapper,
4746  true>
4747  worker(*this,
4748  src,
4749  dst,
4750  false,
4751  wrap,
4752  &Wrapper::cell_integrator,
4753  &Wrapper::face_integrator,
4754  &Wrapper::boundary_integrator,
4757  operation_before_loop,
4758  operation_after_loop,
4759  dof_handler_index_pre_post);
4760 
4761  task_info.loop(worker);
4762 }
4763 
4764 
4765 
4766 template <int dim, typename Number, typename VectorizedArrayType>
4767 template <typename OutVector, typename InVector>
4768 inline void
4770  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4771  OutVector &,
4772  const InVector &,
4773  const std::pair<unsigned int, unsigned int> &)>
4774  &cell_operation,
4775  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4776  OutVector &,
4777  const InVector &,
4778  const std::pair<unsigned int, unsigned int> &)>
4779  &face_operation,
4780  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4781  OutVector &,
4782  const InVector &,
4783  const std::pair<unsigned int, unsigned int> &)>
4784  & boundary_operation,
4785  OutVector & dst,
4786  const InVector & src,
4787  const bool zero_dst_vector,
4788  const DataAccessOnFaces dst_vector_face_access,
4789  const DataAccessOnFaces src_vector_face_access) const
4790 {
4791  using Wrapper =
4792  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4793  InVector,
4794  OutVector>;
4795  Wrapper wrap(cell_operation, face_operation, boundary_operation);
4796  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4797  InVector,
4798  OutVector,
4799  Wrapper,
4800  true>
4801  worker(*this,
4802  src,
4803  dst,
4804  zero_dst_vector,
4805  wrap,
4806  &Wrapper::cell_integrator,
4807  &Wrapper::face_integrator,
4808  &Wrapper::boundary_integrator,
4809  src_vector_face_access,
4810  dst_vector_face_access);
4811 
4812  task_info.loop(worker);
4813 }
4814 
4815 
4816 
4817 template <int dim, typename Number, typename VectorizedArrayType>
4818 template <typename CLASS, typename OutVector, typename InVector>
4819 inline void
4821  void (CLASS::*function_pointer)(
4823  OutVector &,
4824  const InVector &,
4825  const std::pair<unsigned int, unsigned int> &) const,
4826  const CLASS * owning_class,
4827  OutVector & dst,
4828  const InVector &src,
4829  const bool zero_dst_vector) const
4830 {
4831  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4832  InVector,
4833  OutVector,
4834  CLASS,
4835  true>
4836  worker(*this,
4837  src,
4838  dst,
4839  zero_dst_vector,
4840  *owning_class,
4841  function_pointer,
4842  nullptr,
4843  nullptr);
4844  task_info.loop(worker);
4845 }
4846 
4847 
4848 
4849 template <int dim, typename Number, typename VectorizedArrayType>
4850 template <typename CLASS, typename OutVector, typename InVector>
4851 inline void
4853  void (CLASS::*function_pointer)(
4855  OutVector &,
4856  const InVector &,
4857  const std::pair<unsigned int, unsigned int> &) const,
4858  const CLASS * owning_class,
4859  OutVector & dst,
4860  const InVector &src,
4861  const std::function<void(const unsigned int, const unsigned int)>
4862  &operation_before_loop,
4863  const std::function<void(const unsigned int, const unsigned int)>
4864  & operation_after_loop,
4865  const unsigned int dof_handler_index_pre_post) const
4866 {
4867  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4868  InVector,
4869  OutVector,
4870  CLASS,
4871  true>
4872  worker(*this,
4873  src,
4874  dst,
4875  false,
4876  *owning_class,
4877  function_pointer,
4878  nullptr,
4879  nullptr,
4882  operation_before_loop,
4883  operation_after_loop,
4884  dof_handler_index_pre_post);
4885  task_info.loop(worker);
4886 }
4887 
4888 
4889 
4890 template <int dim, typename Number, typename VectorizedArrayType>
4891 template <typename CLASS, typename OutVector, typename InVector>
4892 inline void
4894  void (CLASS::*cell_operation)(
4896  OutVector &,
4897  const InVector &,
4898  const std::pair<unsigned int, unsigned int> &) const,
4899  void (CLASS::*face_operation)(
4901  OutVector &,
4902  const InVector &,
4903  const std::pair<unsigned int, unsigned int> &) const,
4904  void (CLASS::*boundary_operation)(
4906  OutVector &,
4907  const InVector &,
4908  const std::pair<unsigned int, unsigned int> &) const,
4909  const CLASS * owning_class,
4910  OutVector & dst,
4911  const InVector & src,
4912  const bool zero_dst_vector,
4913  const DataAccessOnFaces dst_vector_face_access,
4914  const DataAccessOnFaces src_vector_face_access) const
4915 {
4916  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4917  InVector,
4918  OutVector,
4919  CLASS,
4920  true>
4921  worker(*this,
4922  src,
4923  dst,
4924  zero_dst_vector,
4925  *owning_class,
4926  cell_operation,
4927  face_operation,
4928  boundary_operation,
4929  src_vector_face_access,
4930  dst_vector_face_access);
4931  task_info.loop(worker);
4932 }
4933 
4934 
4935 
4936 template <int dim, typename Number, typename VectorizedArrayType>
4937 template <typename CLASS, typename OutVector, typename InVector>
4938 inline void
4940  void (CLASS::*function_pointer)(
4942  OutVector &,
4943  const InVector &,
4944  const std::pair<unsigned int, unsigned int> &),
4945  CLASS * owning_class,
4946  OutVector & dst,
4947  const InVector &src,
4948  const bool zero_dst_vector) const
4949 {
4950  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4951  InVector,
4952  OutVector,
4953  CLASS,
4954  false>
4955  worker(*this,
4956  src,
4957  dst,
4958  zero_dst_vector,
4959  *owning_class,
4960  function_pointer,
4961  nullptr,
4962  nullptr);
4963  task_info.loop(worker);
4964 }
4965 
4966 
4967 
4968 template <int dim, typename Number, typename VectorizedArrayType>
4969 template <typename CLASS, typename OutVector, typename InVector>
4970 inline void
4972  void (CLASS::*function_pointer)(
4974  OutVector &,
4975  const InVector &,
4976  const std::pair<unsigned int, unsigned int> &),
4977  CLASS * owning_class,
4978  OutVector & dst,
4979  const InVector &src,
4980  const std::function<void(const unsigned int, const unsigned int)>
4981  &operation_before_loop,
4982  const std::function<void(const unsigned int, const unsigned int)>
4983  & operation_after_loop,
4984  const unsigned int dof_handler_index_pre_post) const
4985 {
4986  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4987  InVector,
4988  OutVector,
4989  CLASS,
4990  false>
4991  worker(*this,
4992  src,
4993  dst,
4994  false,
4995  *owning_class,
4996  function_pointer,
4997  nullptr,
4998  nullptr,
5001  operation_before_loop,
5002  operation_after_loop,
5003  dof_handler_index_pre_post);
5004  task_info.loop(worker);
5005 }
5006 
5007 
5008 
5009 template <int dim, typename Number, typename VectorizedArrayType>
5010 template <typename CLASS, typename OutVector, typename InVector>
5011 inline void
5013  void (CLASS::*cell_operation)(
5015  OutVector &,
5016  const InVector &,
5017  const std::pair<unsigned int, unsigned int> &),
5018  void (CLASS::*face_operation)(
5020  OutVector &,
5021  const InVector &,
5022  const std::pair<unsigned int, unsigned int> &),
5023  void (CLASS::*boundary_operation)(
5025  OutVector &,
5026  const InVector &,
5027  const std::pair<unsigned int, unsigned int> &),
5028  CLASS * owning_class,
5029  OutVector & dst,
5030  const InVector & src,
5031  const bool zero_dst_vector,
5032  const DataAccessOnFaces dst_vector_face_access,
5033  const DataAccessOnFaces src_vector_face_access) const
5034 {
5035  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5036  InVector,
5037  OutVector,
5038  CLASS,
5039  false>
5040  worker(*this,
5041  src,
5042  dst,
5043  zero_dst_vector,
5044  *owning_class,
5045  cell_operation,
5046  face_operation,
5047  boundary_operation,
5048  src_vector_face_access,
5049  dst_vector_face_access);
5050  task_info.loop(worker);
5051 }
5052 
5053 
5054 
5055 template <int dim, typename Number, typename VectorizedArrayType>
5056 template <typename CLASS, typename OutVector, typename InVector>
5057 inline void
5059  void (CLASS::*function_pointer)(
5061  OutVector &,
5062  const InVector &,
5063  const std::pair<unsigned int, unsigned int> &) const,
5064  const CLASS * owning_class,
5065  OutVector & dst,
5066  const InVector & src,
5067  const bool zero_dst_vector,
5068  const DataAccessOnFaces src_vector_face_access) const
5069 {
5070  auto src_vector_face_access_temp = src_vector_face_access;
5071  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5072  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5073  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5074  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5075 
5076  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5077  InVector,
5078  OutVector,
5079  CLASS,
5080  true>
5081  worker(*this,
5082  src,
5083  dst,
5084  zero_dst_vector,
5085  *owning_class,
5086  function_pointer,
5087  nullptr,
5088  nullptr,
5089  src_vector_face_access_temp,
5091  task_info.loop(worker);
5092 }
5093 
5094 
5095 
5096 template <int dim, typename Number, typename VectorizedArrayType>
5097 template <typename CLASS, typename OutVector, typename InVector>
5098 inline void
5100  void (CLASS::*function_pointer)(
5102  OutVector &,
5103  const InVector &,
5104  const std::pair<unsigned int, unsigned int> &),
5105  CLASS * owning_class,
5106  OutVector & dst,
5107  const InVector & src,
5108  const bool zero_dst_vector,
5109  const DataAccessOnFaces src_vector_face_access) const
5110 {
5111  auto src_vector_face_access_temp = src_vector_face_access;
5112  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5113  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5114  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5115  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5116 
5117  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5118  InVector,
5119  OutVector,
5120  CLASS,
5121  false>
5122  worker(*this,
5123  src,
5124  dst,
5125  zero_dst_vector,
5126  *owning_class,
5127  function_pointer,
5128  nullptr,
5129  nullptr,
5130  src_vector_face_access_temp,
5132  task_info.loop(worker);
5133 }
5134 
5135 
5136 
5137 template <int dim, typename Number, typename VectorizedArrayType>
5138 template <typename OutVector, typename InVector>
5139 inline void
5141  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5142  OutVector &,
5143  const InVector &,
5144  const std::pair<unsigned int, unsigned int> &)>
5145  & cell_operation,
5146  OutVector & dst,
5147  const InVector & src,
5148  const bool zero_dst_vector,
5149  const DataAccessOnFaces src_vector_face_access) const
5150 {
5151  auto src_vector_face_access_temp = src_vector_face_access;
5152  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5153  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5154  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5155  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5156 
5157  using Wrapper =
5158  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5159  InVector,
5160  OutVector>;
5161  Wrapper wrap(cell_operation, nullptr, nullptr);
5162 
5163  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5164  InVector,
5165  OutVector,
5166  Wrapper,
5167  true>
5168  worker(*this,
5169  src,
5170  dst,
5171  zero_dst_vector,
5172  wrap,
5173  &Wrapper::cell_integrator,
5174  &Wrapper::face_integrator,
5175  &Wrapper::boundary_integrator,
5176  src_vector_face_access_temp,
5178  task_info.loop(worker);
5179 }
5180 
5181 
5182 #endif // ifndef DOXYGEN
5183 
5184 
5185 
5187 
5188 #endif
void resize(const size_type new_size)
Abstract base class for mapping classes.
Definition: mapping.h:311
unsigned int n_active_fe_indices() const
void loop_cell_centric(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > mapping_info
Definition: matrix_free.h:2030
unsigned int get_cell_active_fe_index(const std::pair< unsigned int, unsigned int > range) const
AlignedVector< Number > * acquire_scratch_data_non_threadsafe() const
unsigned int n_ghost_cell_batches() const
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void initialize_dof_handlers(const std::vector< const DoFHandler< dim, dim > * > &dof_handlers, const AdditionalData &additional_data)
types::boundary_id get_boundary_id(const unsigned int face_batch_index) const
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > > shape_info
Definition: matrix_free.h:2036
const Quadrature< dim > & get_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
std::pair< unsigned int, unsigned int > create_cell_subrange_hp(const std::pair< unsigned int, unsigned int > &range, const unsigned int fe_degree, const unsigned int dof_handler_index=0) const
const IndexSet & get_ghost_set(const unsigned int dof_handler_index=0) const
void print(std::ostream &out) const
const Quadrature< dim - 1 > & get_face_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int n_inner_face_batches() const
void internal_reinit(const std::shared_ptr< hp::MappingCollection< dim >> &mapping, const std::vector< const DoFHandler< dim, dim > * > &dof_handlers, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const std::vector< hp::QCollection< q_dim >> &quad, const AdditionalData &additional_data)
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
void update_mapping(const Mapping< dim > &mapping)
unsigned int get_mg_level() const
void print_memory_consumption(StreamType &out) const
void update_mapping(const std::shared_ptr< hp::MappingCollection< dim >> &mapping)
bool mapping_initialized() const
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
const Number * constraint_pool_begin(const unsigned int pool_index) const
void clear()
unsigned int get_n_q_points_face(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
std::pair< unsigned int, unsigned int > create_cell_subrange_hp_by_index(const std::pair< unsigned int, unsigned int > &range, const unsigned int fe_index, const unsigned int dof_handler_index=0) const
bool at_irregular_cell(const unsigned int cell_batch_index) const
~MatrixFree() override=default
void copy_from(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free_base)
void initialize_face_data_vector(AlignedVector< T > &vec) const
unsigned int get_dofs_per_cell(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int n_constraint_pool_entries() const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
std::pair< int, int > get_cell_level_and_index(const unsigned int cell_batch_index, const unsigned int lane_index) const
std::pair< unsigned int, unsigned int > get_face_range_category(const std::pair< unsigned int, unsigned int > face_batch_range) const
const std::vector< unsigned int > & get_constrained_dofs(const unsigned int dof_handler_index=0) const
std::pair< unsigned int, unsigned int > get_face_category(const unsigned int face_batch_index) const
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
unsigned int get_dofs_per_face(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
const Number * constraint_pool_end(const unsigned int pool_index) const
internal::MatrixFreeFunctions::TaskInfo task_info
Definition: matrix_free.h:2060
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*boundary_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
std::pair< typename DoFHandler< dim >::cell_iterator, unsigned int > get_face_iterator(const unsigned int face_batch_index, const unsigned int lane_index, const bool interior=true, const unsigned int fe_component=0) const
bool mapping_is_initialized
Definition: matrix_free.h:2077
void loop_cell_centric(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void initialize_dof_vector(LinearAlgebra::distributed::Vector< Number2 > &vec, const unsigned int dof_handler_index=0) const
MatrixFree(const MatrixFree< dim, Number, VectorizedArrayType > &other)
unsigned int mg_level
Definition: matrix_free.h:2100
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
std::vector< std::pair< unsigned int, unsigned int > > cell_level_index
Definition: matrix_free.h:2044
std::array< types::boundary_id, VectorizedArrayType::size()> get_faces_by_cells_boundary_id(const unsigned int cell_batch_index, const unsigned int face_number) const
std::vector< Number > constraint_pool_data
Definition: matrix_free.h:2017
VectorizedArrayType vectorized_value_type
Definition: matrix_free.h:125
unsigned int n_cell_batches() const
bool indices_initialized() const
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
unsigned int get_cell_category(const unsigned int cell_batch_index) const
Threads::ThreadLocalStorage< std::list< std::pair< bool, AlignedVector< VectorizedArrayType > > > > scratch_pad
Definition: matrix_free.h:2088
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_index) const
Number value_type
Definition: matrix_free.h:124
std::size_t memory_consumption() const
unsigned int n_boundary_face_batches() const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
std::list< std::pair< bool, AlignedVector< Number > > > scratch_pad_non_threadsafe
Definition: matrix_free.h:2095
DoFHandler< dim >::cell_iterator get_cell_iterator(const unsigned int cell_batch_index, const unsigned int lane_index, const unsigned int dof_handler_index=0) const
void loop_cell_centric(const std::function< void(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
std::vector< SmartPointer< const DoFHandler< dim > > > dof_handlers
Definition: matrix_free.h:2003
void initialize_indices(const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const AdditionalData &additional_data)
unsigned int n_components() const
unsigned int n_physical_cells() const
void initialize_dof_vector(VectorType &vec, const unsigned int dof_handler_index=0) const
unsigned int get_n_q_points(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
void reinit(const MappingType &mapping, const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
std::vector< internal::MatrixFreeFunctions::DoFInfo > dof_info
Definition: matrix_free.h:2009
void initialize_cell_data_vector(AlignedVector< T > &vec) const
unsigned int n_ghost_inner_face_batches() const
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
void release_scratch_data_non_threadsafe(const AlignedVector< Number > *memory) const
void renumber_dofs(std::vector< types::global_dof_index > &renumbering, const unsigned int dof_handler_index=0)
unsigned int get_face_active_fe_index(const std::pair< unsigned int, unsigned int > range, const bool is_interior_face=true) const
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
bool indices_are_initialized
Definition: matrix_free.h:2072
internal::MatrixFreeFunctions::FaceInfo< VectorizedArrayType::size()> face_info
Definition: matrix_free.h:2067
void reinit(const MappingType &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< number2 > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
unsigned int get_cell_range_category(const std::pair< unsigned int, unsigned int > cell_batch_range) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*boundary_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void reinit(const MappingType &mapping, const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< QuadratureType > &quad, const AdditionalData &additional_data=AdditionalData())
std::vector< unsigned int > constraint_pool_row_index
Definition: matrix_free.h:2023
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner(const unsigned int dof_handler_index=0) const
unsigned int cell_level_index_end_local
Definition: matrix_free.h:2053
unsigned int n_base_elements(const unsigned int dof_handler_index) const
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
const IndexSet & get_locally_owned_set(const unsigned int dof_handler_index=0) const
static bool is_supported(const FiniteElement< dim, spacedim > &fe)
static constexpr unsigned int dimension
Definition: matrix_free.h:130
A class that provides a separate storage location on each thread that accesses the object.
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
UpdateFlags
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_default
No update.
Point< 2 > second
Definition: grid_out.cc:4605
Point< 2 > first
Definition: grid_out.cc:4604
unsigned int level
Definition: grid_out.cc:4607
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1790
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:466
Number local_element(const size_type local_index) const
void reinit(const size_type size, const bool omit_zeroing_entries=false)
static const char U
VectorType::value_type * begin(VectorType &V)
bool job_supports_mpi()
Definition: mpi.cc:1027
unsigned int minimum_parallel_grain_size
Definition: parallel.cc:34
const types::boundary_id invalid_boundary_id
Definition: types.h:244
static const unsigned int invalid_unsigned_int
Definition: types.h:201
void apply_to_subranges(const tbb::blocked_range< RangeType > &range, const Function &f)
Definition: parallel.h:302
unsigned int boundary_id
Definition: types.h:129
TasksParallelScheme tasks_parallel_scheme
Definition: matrix_free.h:342
UpdateFlags mapping_update_flags_inner_faces
Definition: matrix_free.h:406
AdditionalData(const TasksParallelScheme tasks_parallel_scheme=partition_partition, const unsigned int tasks_block_size=0, const UpdateFlags mapping_update_flags=update_gradients|update_JxW_values, const UpdateFlags mapping_update_flags_boundary_faces=update_default, const UpdateFlags mapping_update_flags_inner_faces=update_default, const UpdateFlags mapping_update_flags_faces_by_cells=update_default, const unsigned int mg_level=numbers::invalid_unsigned_int, const bool store_plain_indices=true, const bool initialize_indices=true, const bool initialize_mapping=true, const bool overlap_communication_computation=true, const bool hold_all_faces_to_owned_cells=false, const bool cell_vectorization_categories_strict=false, const bool allow_ghosted_vectors_in_loops=true)
Definition: matrix_free.h:216
AdditionalData & operator=(const AdditionalData &other)
Definition: matrix_free.h:280
std::vector< unsigned int > cell_vectorization_category
Definition: matrix_free.h:518
UpdateFlags mapping_update_flags_boundary_faces
Definition: matrix_free.h:386
UpdateFlags mapping_update_flags
Definition: matrix_free.h:366
UpdateFlags mapping_update_flags_faces_by_cells
Definition: matrix_free.h:434
AdditionalData(const AdditionalData &other)
Definition: matrix_free.h:253
unsigned int tasks_block_size
Definition: matrix_free.h:353
static bool equal(const T *p1, const T *p2)
std::vector< unsigned int > cell_loop_pre_list_index
Definition: dof_info.h:785
std::vector< unsigned int > cell_loop_post_list_index
Definition: dof_info.h:798
std::vector< std::pair< unsigned int, unsigned int > > vector_zero_range_list
Definition: dof_info.h:778
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition: dof_info.h:629
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_pre_list
Definition: dof_info.h:791
std::vector< unsigned int > vector_zero_range_list_index
Definition: dof_info.h:773
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_post_list
Definition: dof_info.h:804
::Table< 3, unsigned int > cell_and_face_to_plain_faces
Definition: face_info.h:172
std::vector< FaceToCellTopology< vectorization_width > > faces
Definition: face_info.h:165
::Table< 3, types::boundary_id > cell_and_face_boundary_id
Definition: face_info.h:178
std::vector< unsigned int > boundary_partition_data
Definition: task_info.h:521
void loop(MFWorkerInterface &worker) const
Definition: task_info.cc:350
std::vector< unsigned int > cell_partition_data
Definition: task_info.h:477
std::vector< unsigned int > face_partition_data
Definition: task_info.h:499