Reference documentation for deal.II version GIT c415589cf0 2022-08-14 18:50:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
matrix_free.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2011 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii_matrix_free_h
18 #define dealii_matrix_free_h
19 
20 #include <deal.II/base/config.h>
21 
28 
30 
31 #include <deal.II/fe/fe.h>
32 #include <deal.II/fe/mapping.h>
33 #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 
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 
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 
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 
1640  get_affine_constraints(const unsigned int dof_handler_index = 0) const;
1641 
1655  get_cell_iterator(const unsigned int cell_batch_index,
1656  const unsigned int lane_index,
1657  const unsigned int dof_handler_index = 0) const;
1658 
1664  std::pair<int, int>
1665  get_cell_level_and_index(const unsigned int cell_batch_index,
1666  const unsigned int lane_index) const;
1667 
1680  std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>
1681  get_face_iterator(const unsigned int face_batch_index,
1682  const unsigned int lane_index,
1683  const bool interior = true,
1684  const unsigned int fe_component = 0) const;
1685 
1698  bool
1699  at_irregular_cell(const unsigned int cell_batch_index) const;
1700 
1710  unsigned int
1711  n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const;
1712 
1722  unsigned int
1723  n_active_entries_per_face_batch(const unsigned int face_batch_index) const;
1724 
1728  unsigned int
1729  get_dofs_per_cell(const unsigned int dof_handler_index = 0,
1730  const unsigned int hp_active_fe_index = 0) const;
1731 
1735  unsigned int
1736  get_n_q_points(const unsigned int quad_index = 0,
1737  const unsigned int hp_active_fe_index = 0) const;
1738 
1743  unsigned int
1744  get_dofs_per_face(const unsigned int dof_handler_index = 0,
1745  const unsigned int hp_active_fe_index = 0) const;
1746 
1751  unsigned int
1752  get_n_q_points_face(const unsigned int quad_index = 0,
1753  const unsigned int hp_active_fe_index = 0) const;
1754 
1758  const Quadrature<dim> &
1759  get_quadrature(const unsigned int quad_index = 0,
1760  const unsigned int hp_active_fe_index = 0) const;
1761 
1765  const Quadrature<dim - 1> &
1766  get_face_quadrature(const unsigned int quad_index = 0,
1767  const unsigned int hp_active_fe_index = 0) const;
1768 
1781  unsigned int
1783  const std::pair<unsigned int, unsigned int> cell_batch_range) const;
1784 
1789  std::pair<unsigned int, unsigned int>
1791  const std::pair<unsigned int, unsigned int> face_batch_range) const;
1792 
1804  unsigned int
1805  get_cell_category(const unsigned int cell_batch_index) const;
1806 
1811  std::pair<unsigned int, unsigned int>
1812  get_face_category(const unsigned int face_batch_index) const;
1813 
1817  bool
1819 
1824  bool
1826 
1831  unsigned int
1832  get_mg_level() const;
1833 
1838  std::size_t
1840 
1845  template <typename StreamType>
1846  void
1847  print_memory_consumption(StreamType &out) const;
1848 
1853  void
1854  print(std::ostream &out) const;
1855 
1868  get_task_info() const;
1869 
1870  /*
1871  * Return geometry-dependent information on the cells.
1872  */
1873  const internal::MatrixFreeFunctions::
1874  MappingInfo<dim, Number, VectorizedArrayType> &
1876 
1881  get_dof_info(const unsigned int dof_handler_index_component = 0) const;
1882 
1886  unsigned int
1888 
1893  const Number *
1894  constraint_pool_begin(const unsigned int pool_index) const;
1895 
1901  const Number *
1902  constraint_pool_end(const unsigned int pool_index) const;
1903 
1908  get_shape_info(const unsigned int dof_handler_index_component = 0,
1909  const unsigned int quad_index = 0,
1910  const unsigned int fe_base_element = 0,
1911  const unsigned int hp_active_fe_index = 0,
1912  const unsigned int hp_active_quad_index = 0) const;
1913 
1918  VectorizedArrayType::size()> &
1919  get_face_info(const unsigned int face_batch_index) const;
1920 
1921 
1927  const Table<3, unsigned int> &
1929 
1945 
1949  void
1951 
1963 
1967  void
1969  const AlignedVector<Number> *memory) const;
1970 
1973 private:
1978  template <typename number2, int q_dim>
1979  void
1981  const std::shared_ptr<hp::MappingCollection<dim>> & mapping,
1982  const std::vector<const DoFHandler<dim, dim> *> & dof_handlers,
1983  const std::vector<const AffineConstraints<number2> *> &constraint,
1984  const std::vector<IndexSet> & locally_owned_set,
1985  const std::vector<hp::QCollection<q_dim>> & quad,
1986  const AdditionalData & additional_data);
1987 
1994  template <typename number2>
1995  void
1997  const std::vector<const AffineConstraints<number2> *> &constraint,
1998  const std::vector<IndexSet> & locally_owned_set,
1999  const AdditionalData & additional_data);
2000 
2004  void
2006  const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
2007  const AdditionalData & additional_data);
2008 
2012  std::vector<SmartPointer<const DoFHandler<dim>>> dof_handlers;
2013 
2020  std::vector<SmartPointer<const AffineConstraints<Number>>> affine_constraints;
2021 
2026  std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
2027 
2034  std::vector<Number> constraint_pool_data;
2035 
2040  std::vector<unsigned int> constraint_pool_row_index;
2041 
2048 
2054 
2061  std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
2062 
2063 
2071 
2078 
2083  internal::MatrixFreeFunctions::FaceInfo<VectorizedArrayType::size()>
2085 
2090 
2095 
2104  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2106 
2111  mutable std::list<std::pair<bool, AlignedVector<Number>>>
2113 
2117  unsigned int mg_level;
2118 };
2119 
2120 
2121 
2122 /*----------------------- Inline functions ----------------------------------*/
2123 
2124 #ifndef DOXYGEN
2125 
2126 
2127 
2128 template <int dim, typename Number, typename VectorizedArrayType>
2129 template <typename T>
2130 inline void
2132  AlignedVector<T> &vec) const
2133 {
2134  vec.resize(this->n_cell_batches() + this->n_ghost_cell_batches());
2135 }
2136 
2137 
2138 
2139 template <int dim, typename Number, typename VectorizedArrayType>
2140 template <typename T>
2141 inline void
2143  AlignedVector<T> &vec) const
2144 {
2145  vec.resize(this->n_inner_face_batches() + this->n_boundary_face_batches() +
2146  this->n_ghost_inner_face_batches());
2147 }
2148 
2149 
2150 
2151 template <int dim, typename Number, typename VectorizedArrayType>
2152 template <typename VectorType>
2153 inline void
2155  VectorType & vec,
2156  const unsigned int comp) const
2157 {
2158  static_assert(IsBlockVector<VectorType>::value == false,
2159  "This function is not supported for block vectors.");
2160 
2161  Assert(task_info.n_procs == 1,
2162  ExcMessage("This function can only be used in serial."));
2163 
2164  AssertIndexRange(comp, n_components());
2165  vec.reinit(dof_info[comp].vector_partitioner->size());
2166 }
2167 
2168 
2169 
2170 template <int dim, typename Number, typename VectorizedArrayType>
2171 template <typename Number2>
2172 inline void
2175  const unsigned int comp) const
2176 {
2177  AssertIndexRange(comp, n_components());
2178  vec.reinit(dof_info[comp].vector_partitioner, task_info.communicator_sm);
2179 }
2180 
2181 
2182 
2183 template <int dim, typename Number, typename VectorizedArrayType>
2184 inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2186  const unsigned int comp) const
2187 {
2188  AssertIndexRange(comp, n_components());
2189  return dof_info[comp].vector_partitioner;
2190 }
2191 
2192 
2193 
2194 template <int dim, typename Number, typename VectorizedArrayType>
2195 inline const std::vector<unsigned int> &
2197  const unsigned int comp) const
2198 {
2199  AssertIndexRange(comp, n_components());
2200  return dof_info[comp].constrained_dofs;
2201 }
2202 
2203 
2204 
2205 template <int dim, typename Number, typename VectorizedArrayType>
2206 inline unsigned int
2208 {
2209  AssertDimension(dof_handlers.size(), dof_info.size());
2210  return dof_handlers.size();
2211 }
2212 
2213 
2214 
2215 template <int dim, typename Number, typename VectorizedArrayType>
2216 inline unsigned int
2218  const unsigned int dof_no) const
2219 {
2220  AssertDimension(dof_handlers.size(), dof_info.size());
2221  AssertIndexRange(dof_no, dof_handlers.size());
2222  return dof_handlers[dof_no]->get_fe().n_base_elements();
2223 }
2224 
2225 
2226 
2227 template <int dim, typename Number, typename VectorizedArrayType>
2230 {
2231  return task_info;
2232 }
2233 
2234 
2235 
2236 template <int dim, typename Number, typename VectorizedArrayType>
2237 inline unsigned int
2239 {
2240  return task_info.n_active_cells;
2241 }
2242 
2243 
2244 
2245 template <int dim, typename Number, typename VectorizedArrayType>
2246 inline unsigned int
2248 {
2249  return *(task_info.cell_partition_data.end() - 2);
2250 }
2251 
2252 
2253 
2254 template <int dim, typename Number, typename VectorizedArrayType>
2255 inline unsigned int
2257 {
2258  return *(task_info.cell_partition_data.end() - 1) -
2259  *(task_info.cell_partition_data.end() - 2);
2260 }
2261 
2262 
2263 
2264 template <int dim, typename Number, typename VectorizedArrayType>
2265 inline unsigned int
2267 {
2268  if (task_info.face_partition_data.size() == 0)
2269  return 0;
2270  return task_info.face_partition_data.back();
2271 }
2272 
2273 
2274 
2275 template <int dim, typename Number, typename VectorizedArrayType>
2276 inline unsigned int
2278 {
2279  if (task_info.face_partition_data.size() == 0)
2280  return 0;
2281  return task_info.boundary_partition_data.back() -
2283 }
2284 
2285 
2286 
2287 template <int dim, typename Number, typename VectorizedArrayType>
2288 inline unsigned int
2290 {
2291  if (task_info.face_partition_data.size() == 0)
2292  return 0;
2293  return face_info.faces.size() - task_info.boundary_partition_data.back();
2294 }
2295 
2296 
2297 
2298 template <int dim, typename Number, typename VectorizedArrayType>
2299 inline types::boundary_id
2301  const unsigned int face_batch_index) const
2302 {
2303  Assert(face_batch_index >= task_info.boundary_partition_data[0] &&
2304  face_batch_index < task_info.boundary_partition_data.back(),
2305  ExcIndexRange(face_batch_index,
2308  return types::boundary_id(face_info.faces[face_batch_index].exterior_face_no);
2309 }
2310 
2311 
2312 
2313 template <int dim, typename Number, typename VectorizedArrayType>
2314 inline std::array<types::boundary_id, VectorizedArrayType::size()>
2316  const unsigned int cell_batch_index,
2317  const unsigned int face_number) const
2318 {
2319  AssertIndexRange(cell_batch_index, n_cell_batches());
2322  ExcNotInitialized());
2323  std::array<types::boundary_id, VectorizedArrayType::size()> result;
2324  result.fill(numbers::invalid_boundary_id);
2325  for (unsigned int v = 0;
2326  v < n_active_entries_per_cell_batch(cell_batch_index);
2327  ++v)
2328  result[v] =
2329  face_info.cell_and_face_boundary_id(cell_batch_index, face_number, v);
2330  return result;
2331 }
2332 
2333 
2334 
2335 template <int dim, typename Number, typename VectorizedArrayType>
2336 inline const internal::MatrixFreeFunctions::
2337  MappingInfo<dim, Number, VectorizedArrayType> &
2339 {
2340  return mapping_info;
2341 }
2342 
2343 
2344 
2345 template <int dim, typename Number, typename VectorizedArrayType>
2348  const unsigned int dof_index) const
2349 {
2350  AssertIndexRange(dof_index, n_components());
2351  return dof_info[dof_index];
2352 }
2353 
2354 
2355 
2356 template <int dim, typename Number, typename VectorizedArrayType>
2357 inline unsigned int
2359 {
2360  return constraint_pool_row_index.size() - 1;
2361 }
2362 
2363 
2364 
2365 template <int dim, typename Number, typename VectorizedArrayType>
2366 inline const Number *
2368  const unsigned int row) const
2369 {
2371  return constraint_pool_data.empty() ?
2372  nullptr :
2374 }
2375 
2376 
2377 
2378 template <int dim, typename Number, typename VectorizedArrayType>
2379 inline const Number *
2381  const unsigned int row) const
2382 {
2384  return constraint_pool_data.empty() ?
2385  nullptr :
2387 }
2388 
2389 
2390 
2391 template <int dim, typename Number, typename VectorizedArrayType>
2392 inline std::pair<unsigned int, unsigned int>
2394  const std::pair<unsigned int, unsigned int> &range,
2395  const unsigned int degree,
2396  const unsigned int dof_handler_component) const
2397 {
2398  if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2399  {
2401  dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2403  dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2404  if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2405  return range;
2406  else
2407  return {range.second, range.second};
2408  }
2409 
2410  const unsigned int fe_index =
2411  dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2412  if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2413  return {range.second, range.second};
2414  else
2415  return create_cell_subrange_hp_by_index(range,
2416  fe_index,
2417  dof_handler_component);
2418 }
2419 
2420 
2421 
2422 template <int dim, typename Number, typename VectorizedArrayType>
2423 inline bool
2425  const unsigned int cell_batch_index) const
2426 {
2427  AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2428  return VectorizedArrayType::size() > 1 &&
2429  cell_level_index[(cell_batch_index + 1) * VectorizedArrayType::size() -
2430  1] == cell_level_index[(cell_batch_index + 1) *
2431  VectorizedArrayType::size() -
2432  2];
2433 }
2434 
2435 
2436 
2437 template <int dim, typename Number, typename VectorizedArrayType>
2438 unsigned int
2440 {
2441  return shape_info.size(2);
2442 }
2443 
2444 
2445 template <int dim, typename Number, typename VectorizedArrayType>
2446 unsigned int
2448  const std::pair<unsigned int, unsigned int> range) const
2449 {
2450  const auto &fe_indices = dof_info[0].cell_active_fe_index;
2451 
2452  if (fe_indices.empty() == true)
2453  return 0;
2454 
2455  const auto index = fe_indices[range.first];
2456 
2457  for (unsigned int i = range.first; i < range.second; ++i)
2458  AssertDimension(index, fe_indices[i]);
2459 
2460  return index;
2461 }
2462 
2463 
2464 
2465 template <int dim, typename Number, typename VectorizedArrayType>
2466 unsigned int
2468  const std::pair<unsigned int, unsigned int> range,
2469  const bool is_interior_face) const
2470 {
2471  const auto &fe_indices = dof_info[0].cell_active_fe_index;
2472 
2473  if (fe_indices.empty() == true)
2474  return 0;
2475 
2476  if (is_interior_face)
2477  {
2478  const unsigned int index =
2479  fe_indices[face_info.faces[range.first].cells_interior[0] /
2480  VectorizedArrayType::size()];
2481 
2482  for (unsigned int i = range.first; i < range.second; ++i)
2483  AssertDimension(index,
2484  fe_indices[face_info.faces[i].cells_interior[0] /
2485  VectorizedArrayType::size()]);
2486 
2487  return index;
2488  }
2489  else
2490  {
2491  const unsigned int index =
2492  fe_indices[face_info.faces[range.first].cells_exterior[0] /
2493  VectorizedArrayType::size()];
2494 
2495  for (unsigned int i = range.first; i < range.second; ++i)
2496  AssertDimension(index,
2497  fe_indices[face_info.faces[i].cells_exterior[0] /
2498  VectorizedArrayType::size()]);
2499 
2500  return index;
2501  }
2502 }
2503 
2504 
2505 
2506 template <int dim, typename Number, typename VectorizedArrayType>
2507 inline unsigned int
2509  const unsigned int cell_batch_index) const
2510 {
2511  Assert(!dof_info.empty(), ExcNotInitialized());
2512  AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2513  const std::vector<unsigned char> &n_lanes_filled =
2514  dof_info[0].n_vectorization_lanes_filled
2516  AssertIndexRange(cell_batch_index, n_lanes_filled.size());
2517 
2518  return n_lanes_filled[cell_batch_index];
2519 }
2520 
2521 
2522 
2523 template <int dim, typename Number, typename VectorizedArrayType>
2524 inline unsigned int
2526  const unsigned int face_batch_index) const
2527 {
2528  AssertIndexRange(face_batch_index, face_info.faces.size());
2529  Assert(!dof_info.empty(), ExcNotInitialized());
2530  const std::vector<unsigned char> &n_lanes_filled =
2531  dof_info[0].n_vectorization_lanes_filled
2533  AssertIndexRange(face_batch_index, n_lanes_filled.size());
2534  return n_lanes_filled[face_batch_index];
2535 }
2536 
2537 
2538 
2539 template <int dim, typename Number, typename VectorizedArrayType>
2540 inline unsigned int
2542  const unsigned int dof_handler_index,
2543  const unsigned int active_fe_index) const
2544 {
2545  return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2546 }
2547 
2548 
2549 
2550 template <int dim, typename Number, typename VectorizedArrayType>
2551 inline unsigned int
2553  const unsigned int quad_index,
2554  const unsigned int active_fe_index) const
2555 {
2556  AssertIndexRange(quad_index, mapping_info.cell_data.size());
2557  return mapping_info.cell_data[quad_index]
2558  .descriptor[active_fe_index]
2559  .n_q_points;
2560 }
2561 
2562 
2563 
2564 template <int dim, typename Number, typename VectorizedArrayType>
2565 inline unsigned int
2567  const unsigned int dof_handler_index,
2568  const unsigned int active_fe_index) const
2569 {
2570  return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2571 }
2572 
2573 
2574 
2575 template <int dim, typename Number, typename VectorizedArrayType>
2576 inline unsigned int
2578  const unsigned int quad_index,
2579  const unsigned int active_fe_index) const
2580 {
2581  AssertIndexRange(quad_index, mapping_info.face_data.size());
2582  return mapping_info.face_data[quad_index]
2583  .descriptor[active_fe_index]
2584  .n_q_points;
2585 }
2586 
2587 
2588 
2589 template <int dim, typename Number, typename VectorizedArrayType>
2590 inline const IndexSet &
2592  const unsigned int dof_handler_index) const
2593 {
2594  return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2595 }
2596 
2597 
2598 
2599 template <int dim, typename Number, typename VectorizedArrayType>
2600 inline const IndexSet &
2602  const unsigned int dof_handler_index) const
2603 {
2604  return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2605 }
2606 
2607 
2608 
2609 template <int dim, typename Number, typename VectorizedArrayType>
2612  const unsigned int dof_handler_index,
2613  const unsigned int index_quad,
2614  const unsigned int index_fe,
2615  const unsigned int active_fe_index,
2616  const unsigned int active_quad_index) const
2617 {
2618  AssertIndexRange(dof_handler_index, dof_info.size());
2619  const unsigned int ind =
2620  dof_info[dof_handler_index].global_base_element_offset + index_fe;
2621  AssertIndexRange(ind, shape_info.size(0));
2622  AssertIndexRange(index_quad, shape_info.size(1));
2623  AssertIndexRange(active_fe_index, shape_info.size(2));
2624  AssertIndexRange(active_quad_index, shape_info.size(3));
2625  return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2626 }
2627 
2628 
2629 
2630 template <int dim, typename Number, typename VectorizedArrayType>
2632  VectorizedArrayType::size()> &
2634  const unsigned int face_batch_index) const
2635 {
2636  AssertIndexRange(face_batch_index, face_info.faces.size());
2637  return face_info.faces[face_batch_index];
2638 }
2639 
2640 
2641 
2642 template <int dim, typename Number, typename VectorizedArrayType>
2643 inline const Table<3, unsigned int> &
2645  const
2646 {
2648 }
2649 
2650 
2651 
2652 template <int dim, typename Number, typename VectorizedArrayType>
2653 inline const Quadrature<dim> &
2655  const unsigned int quad_index,
2656  const unsigned int active_fe_index) const
2657 {
2658  AssertIndexRange(quad_index, mapping_info.cell_data.size());
2659  return mapping_info.cell_data[quad_index]
2660  .descriptor[active_fe_index]
2661  .quadrature;
2662 }
2663 
2664 
2665 
2666 template <int dim, typename Number, typename VectorizedArrayType>
2667 inline const Quadrature<dim - 1> &
2669  const unsigned int quad_index,
2670  const unsigned int active_fe_index) const
2671 {
2672  AssertIndexRange(quad_index, mapping_info.face_data.size());
2673  return mapping_info.face_data[quad_index]
2674  .descriptor[active_fe_index]
2675  .quadrature;
2676 }
2677 
2678 
2679 
2680 template <int dim, typename Number, typename VectorizedArrayType>
2681 inline unsigned int
2683  const std::pair<unsigned int, unsigned int> range) const
2684 {
2685  auto result = get_cell_category(range.first);
2686 
2687  for (unsigned int i = range.first; i < range.second; ++i)
2688  result = std::max(result, get_cell_category(i));
2689 
2690  return result;
2691 }
2692 
2693 
2694 
2695 template <int dim, typename Number, typename VectorizedArrayType>
2696 inline std::pair<unsigned int, unsigned int>
2698  const std::pair<unsigned int, unsigned int> range) const
2699 {
2700  auto result = get_face_category(range.first);
2701 
2702  for (unsigned int i = range.first; i < range.second; ++i)
2703  {
2704  result.first = std::max(result.first, get_face_category(i).first);
2705  result.second = std::max(result.second, get_face_category(i).second);
2706  }
2707 
2708  return result;
2709 }
2710 
2711 
2712 
2713 template <int dim, typename Number, typename VectorizedArrayType>
2714 inline unsigned int
2716  const unsigned int cell_batch_index) const
2717 {
2718  AssertIndexRange(0, dof_info.size());
2719  AssertIndexRange(cell_batch_index, dof_info[0].cell_active_fe_index.size());
2720  if (dof_info[0].cell_active_fe_index.empty())
2721  return 0;
2722  else
2723  return dof_info[0].cell_active_fe_index[cell_batch_index];
2724 }
2725 
2726 
2727 
2728 template <int dim, typename Number, typename VectorizedArrayType>
2729 inline std::pair<unsigned int, unsigned int>
2731  const unsigned int face_batch_index) const
2732 {
2733  AssertIndexRange(face_batch_index, face_info.faces.size());
2734  if (dof_info[0].cell_active_fe_index.empty())
2735  return std::make_pair(0U, 0U);
2736 
2737  std::pair<unsigned int, unsigned int> result = std::make_pair(0U, 0U);
2738  for (unsigned int v = 0;
2739  v < VectorizedArrayType::size() &&
2740  face_info.faces[face_batch_index].cells_interior[v] !=
2742  ++v)
2743  result.first = std::max(
2744  result.first,
2745  dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
2746  .cells_interior[v] /
2747  VectorizedArrayType::size()]);
2748  if (face_info.faces[face_batch_index].cells_exterior[0] !=
2750  for (unsigned int v = 0;
2751  v < VectorizedArrayType::size() &&
2752  face_info.faces[face_batch_index].cells_exterior[v] !=
2754  ++v)
2755  result.second = std::max(
2756  result.second,
2757  dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
2758  .cells_exterior[v] /
2759  VectorizedArrayType::size()]);
2760  else
2761  result.second = numbers::invalid_unsigned_int;
2762  return result;
2763 }
2764 
2765 
2766 
2767 template <int dim, typename Number, typename VectorizedArrayType>
2768 inline bool
2770 {
2771  return indices_are_initialized;
2772 }
2773 
2774 
2775 
2776 template <int dim, typename Number, typename VectorizedArrayType>
2777 inline bool
2779 {
2780  return mapping_is_initialized;
2781 }
2782 
2783 
2784 template <int dim, typename Number, typename VectorizedArrayType>
2785 inline unsigned int
2787 {
2788  return mg_level;
2789 }
2790 
2791 
2792 
2793 template <int dim, typename Number, typename VectorizedArrayType>
2796 {
2797  using list_type =
2798  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
2799  list_type &data = scratch_pad.get();
2800  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
2801  if (it->first == false)
2802  {
2803  it->first = true;
2804  return &it->second;
2805  }
2806  data.emplace_front(true, AlignedVector<VectorizedArrayType>());
2807  return &data.front().second;
2808 }
2809 
2810 
2811 
2812 template <int dim, typename Number, typename VectorizedArrayType>
2813 void
2815  const AlignedVector<VectorizedArrayType> *scratch) const
2816 {
2817  using list_type =
2818  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
2819  list_type &data = scratch_pad.get();
2820  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
2821  if (&it->second == scratch)
2822  {
2823  Assert(it->first == true, ExcInternalError());
2824  it->first = false;
2825  return;
2826  }
2827  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2828 }
2829 
2830 
2831 
2832 template <int dim, typename Number, typename VectorizedArrayType>
2836 {
2837  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
2839  it != scratch_pad_non_threadsafe.end();
2840  ++it)
2841  if (it->first == false)
2842  {
2843  it->first = true;
2844  return &it->second;
2845  }
2846  scratch_pad_non_threadsafe.push_front(
2847  std::make_pair(true, AlignedVector<Number>()));
2848  return &scratch_pad_non_threadsafe.front().second;
2849 }
2850 
2851 
2852 
2853 template <int dim, typename Number, typename VectorizedArrayType>
2854 void
2857  const AlignedVector<Number> *scratch) const
2858 {
2859  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
2861  it != scratch_pad_non_threadsafe.end();
2862  ++it)
2863  if (&it->second == scratch)
2864  {
2865  Assert(it->first == true, ExcInternalError());
2866  it->first = false;
2867  return;
2868  }
2869  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2870 }
2871 
2872 
2873 
2874 // ------------------------------ reinit functions ---------------------------
2875 
2876 namespace internal
2877 {
2878  namespace MatrixFreeImplementation
2879  {
2880  template <int dim, int spacedim>
2881  inline std::vector<IndexSet>
2882  extract_locally_owned_index_sets(
2883  const std::vector<const ::DoFHandler<dim, spacedim> *> &dofh,
2884  const unsigned int level)
2885  {
2886  std::vector<IndexSet> locally_owned_set;
2887  locally_owned_set.reserve(dofh.size());
2888  for (unsigned int j = 0; j < dofh.size(); ++j)
2890  locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
2891  else
2892  locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
2893  return locally_owned_set;
2894  }
2895  } // namespace MatrixFreeImplementation
2896 } // namespace internal
2897 
2898 
2899 
2900 template <int dim, typename Number, typename VectorizedArrayType>
2901 template <typename QuadratureType, typename number2, typename MappingType>
2902 void
2904  const MappingType & mapping,
2905  const DoFHandler<dim> & dof_handler,
2906  const AffineConstraints<number2> &constraints_in,
2907  const QuadratureType & quad,
2909  &additional_data)
2910 {
2911  std::vector<const DoFHandler<dim, dim> *> dof_handlers;
2912  std::vector<const AffineConstraints<number2> *> constraints;
2913 
2914  dof_handlers.push_back(&dof_handler);
2915  constraints.push_back(&constraints_in);
2916 
2917  std::vector<IndexSet> locally_owned_sets =
2918  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2919  dof_handlers, additional_data.mg_level);
2920 
2921  std::vector<hp::QCollection<dim>> quad_hp;
2922  quad_hp.emplace_back(quad);
2923 
2924  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
2925  dof_handlers,
2926  constraints,
2927  locally_owned_sets,
2928  quad_hp,
2929  additional_data);
2930 }
2931 
2932 
2933 
2934 template <int dim, typename Number, typename VectorizedArrayType>
2935 template <typename QuadratureType, typename number2, typename MappingType>
2936 void
2938  const MappingType & mapping,
2939  const std::vector<const DoFHandler<dim> *> & dof_handler,
2940  const std::vector<const AffineConstraints<number2> *> &constraint,
2941  const QuadratureType & quad,
2943  &additional_data)
2944 {
2945  std::vector<IndexSet> locally_owned_set =
2946  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2947  dof_handler, additional_data.mg_level);
2948  std::vector<hp::QCollection<dim>> quad_hp;
2949  quad_hp.emplace_back(quad);
2950 
2951  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
2952  dof_handler,
2953  constraint,
2954  locally_owned_set,
2955  quad_hp,
2956  additional_data);
2957 }
2958 
2959 
2960 
2961 template <int dim, typename Number, typename VectorizedArrayType>
2962 template <typename QuadratureType, typename number2, typename MappingType>
2963 void
2965  const MappingType & mapping,
2966  const std::vector<const DoFHandler<dim> *> & dof_handler,
2967  const std::vector<const AffineConstraints<number2> *> &constraint,
2968  const std::vector<QuadratureType> & quad,
2970  &additional_data)
2971 {
2972  std::vector<IndexSet> locally_owned_set =
2973  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2974  dof_handler, additional_data.mg_level);
2975  std::vector<hp::QCollection<dim>> quad_hp;
2976  for (unsigned int q = 0; q < quad.size(); ++q)
2977  quad_hp.emplace_back(quad[q]);
2978 
2979  internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
2980  dof_handler,
2981  constraint,
2982  locally_owned_set,
2983  quad_hp,
2984  additional_data);
2985 }
2986 
2987 
2988 
2989 // ------------------------------ implementation of loops --------------------
2990 
2991 // internal helper functions that define how to call MPI data exchange
2992 // functions: for generic vectors, do nothing at all. For distributed vectors,
2993 // call update_ghost_values_start function and so on. If we have collections
2994 // of vectors, just do the individual functions of the components. In order to
2995 // keep ghost values consistent (whether we are in read or write mode), we
2996 // also reset the values at the end. the whole situation is a bit complicated
2997 // by the fact that we need to treat block vectors differently, which use some
2998 // additional helper functions to select the blocks and template magic.
2999 namespace internal
3000 {
3004  template <int dim, typename Number, typename VectorizedArrayType>
3005  struct VectorDataExchange
3006  {
3007  // A shift for the MPI messages to reduce the risk for accidental
3008  // interaction with other open communications that a user program might
3009  // set up (parallel vectors support unfinished communication). We let
3010  // the other vectors use the first 20 assigned numbers and start the
3011  // matrix-free communication.
3012  static constexpr unsigned int channel_shift = 20;
3013 
3014 
3015 
3020  VectorDataExchange(
3021  const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3022  const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3023  DataAccessOnFaces vector_face_access,
3024  const unsigned int n_components)
3025  : matrix_free(matrix_free)
3026  , vector_face_access(
3027  matrix_free.get_task_info().face_partition_data.empty() ?
3028  ::MatrixFree<dim, Number, VectorizedArrayType>::
3029  DataAccessOnFaces::unspecified :
3030  vector_face_access)
3031  , ghosts_were_set(false)
3032 # ifdef DEAL_II_WITH_MPI
3033  , tmp_data(n_components)
3034  , requests(n_components)
3035 # endif
3036  {
3037  (void)n_components;
3038  if (this->vector_face_access !=
3040  DataAccessOnFaces::unspecified)
3041  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3043  matrix_free.get_dof_info(c).vector_exchanger_face_variants.size(),
3044  5);
3045  }
3046 
3047 
3048 
3052  ~VectorDataExchange() // NOLINT
3053  {
3054 # ifdef DEAL_II_WITH_MPI
3055  for (unsigned int i = 0; i < tmp_data.size(); ++i)
3056  if (tmp_data[i] != nullptr)
3057  matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3058 # endif
3059  }
3060 
3061 
3062 
3067  template <typename VectorType>
3068  unsigned int
3069  find_vector_in_mf(const VectorType &vec,
3070  const bool check_global_compatibility = true) const
3071  {
3072  // case 1: vector was set up with MatrixFree::initialize_dof_vector()
3073  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3074  if (vec.get_partitioner().get() ==
3075  matrix_free.get_dof_info(c).vector_partitioner.get())
3076  return c;
3077 
3078  // case 2: user provided own partitioner (compatibility mode)
3079  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3080  if (check_global_compatibility ?
3081  vec.get_partitioner()->is_globally_compatible(
3082  *matrix_free.get_dof_info(c).vector_partitioner) :
3083  vec.get_partitioner()->is_compatible(
3084  *matrix_free.get_dof_info(c).vector_partitioner))
3085  return c;
3086 
3087  Assert(false,
3088  ExcNotImplemented("Could not find partitioner that fits vector"));
3089 
3091  }
3092 
3093 
3094 
3100  get_partitioner(const unsigned int mf_component) const
3101  {
3102  AssertDimension(matrix_free.get_dof_info(mf_component)
3103  .vector_exchanger_face_variants.size(),
3104  5);
3105  if (vector_face_access ==
3108  return *matrix_free.get_dof_info(mf_component)
3109  .vector_exchanger_face_variants[0];
3110  else if (vector_face_access ==
3113  return *matrix_free.get_dof_info(mf_component)
3114  .vector_exchanger_face_variants[1];
3115  else if (vector_face_access ==
3118  return *matrix_free.get_dof_info(mf_component)
3119  .vector_exchanger_face_variants[2];
3120  else if (vector_face_access ==
3122  DataAccessOnFaces::values_all_faces)
3123  return *matrix_free.get_dof_info(mf_component)
3124  .vector_exchanger_face_variants[3];
3125  else if (vector_face_access ==
3127  DataAccessOnFaces::gradients_all_faces)
3128  return *matrix_free.get_dof_info(mf_component)
3129  .vector_exchanger_face_variants[4];
3130  else
3131  return *matrix_free.get_dof_info(mf_component).vector_exchanger.get();
3132  }
3133 
3134 
3135 
3139  template <typename VectorType,
3140  std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3141  * = nullptr>
3142  void
3143  update_ghost_values_start(const unsigned int /*component_in_block_vector*/,
3144  const VectorType & /*vec*/)
3145  {}
3146 
3147 
3152  template <typename VectorType,
3153  std::enable_if_t<!has_update_ghost_values_start<VectorType> &&
3154  !is_not_parallel_vector<VectorType>,
3155  VectorType> * = nullptr>
3156  void
3157  update_ghost_values_start(const unsigned int component_in_block_vector,
3158  const VectorType & vec)
3159  {
3160  (void)component_in_block_vector;
3161  bool ghosts_set = vec.has_ghost_elements();
3162 
3163  Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3164  ghosts_set == false,
3165  ExcNotImplemented());
3166 
3167  if (ghosts_set)
3168  ghosts_were_set = true;
3169 
3170  vec.update_ghost_values();
3171  }
3172 
3173 
3174 
3180  template <typename VectorType,
3181  std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3182  !has_exchange_on_subset<VectorType>,
3183  VectorType> * = nullptr>
3184  void
3185  update_ghost_values_start(const unsigned int component_in_block_vector,
3186  const VectorType & vec)
3187  {
3188  (void)component_in_block_vector;
3189  bool ghosts_set = vec.has_ghost_elements();
3190 
3191  Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3192  ghosts_set == false,
3193  ExcNotImplemented());
3194 
3195  if (ghosts_set)
3196  ghosts_were_set = true;
3197 
3198  vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3199  }
3200 
3201 
3202 
3209  template <typename VectorType,
3210  std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3211  has_exchange_on_subset<VectorType>,
3212  VectorType> * = nullptr>
3213  void
3214  update_ghost_values_start(const unsigned int component_in_block_vector,
3215  const VectorType & vec)
3216  {
3217  static_assert(
3218  std::is_same<Number, typename VectorType::value_type>::value,
3219  "Type mismatch between VectorType and VectorDataExchange");
3220  (void)component_in_block_vector;
3221  bool ghosts_set = vec.has_ghost_elements();
3222 
3223  Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3224  ghosts_set == false,
3225  ExcNotImplemented());
3226 
3227  if (ghosts_set)
3228  ghosts_were_set = true;
3229 
3230  if (vec.size() != 0)
3231  {
3232 # ifdef DEAL_II_WITH_MPI
3233  const unsigned int mf_component = find_vector_in_mf(vec);
3234 
3235  const auto &part = get_partitioner(mf_component);
3236 
3237  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3238  part.n_import_sm_procs() == 0)
3239  return;
3240 
3241  tmp_data[component_in_block_vector] =
3242  matrix_free.acquire_scratch_data_non_threadsafe();
3243  tmp_data[component_in_block_vector]->resize_fast(
3244  part.n_import_indices());
3245  AssertDimension(requests.size(), tmp_data.size());
3246 
3247  part.export_to_ghosted_array_start(
3248  component_in_block_vector * 2 + channel_shift,
3249  ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3250  vec.shared_vector_data(),
3251  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3252  part.locally_owned_size(),
3253  matrix_free.get_dof_info(mf_component)
3254  .vector_partitioner->n_ghost_indices()),
3255  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3256  part.n_import_indices()),
3257  this->requests[component_in_block_vector]);
3258 # endif
3259  }
3260  }
3261 
3262 
3263 
3268  template <typename VectorType,
3269  std::enable_if_t<!has_update_ghost_values_start<VectorType>,
3270  VectorType> * = nullptr>
3271  void
3272  update_ghost_values_finish(const unsigned int /*component_in_block_vector*/,
3273  const VectorType & /*vec*/)
3274  {}
3275 
3276 
3277 
3283  template <typename VectorType,
3284  std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3285  !has_exchange_on_subset<VectorType>,
3286  VectorType> * = nullptr>
3287  void
3288  update_ghost_values_finish(const unsigned int component_in_block_vector,
3289  const VectorType & vec)
3290  {
3291  (void)component_in_block_vector;
3292  vec.update_ghost_values_finish();
3293  }
3294 
3295 
3296 
3303  template <typename VectorType,
3304  std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3305  has_exchange_on_subset<VectorType>,
3306  VectorType> * = nullptr>
3307  void
3308  update_ghost_values_finish(const unsigned int component_in_block_vector,
3309  const VectorType & vec)
3310  {
3311  static_assert(
3312  std::is_same<Number, typename VectorType::value_type>::value,
3313  "Type mismatch between VectorType and VectorDataExchange");
3314  (void)component_in_block_vector;
3315 
3316  if (vec.size() != 0)
3317  {
3318 # ifdef DEAL_II_WITH_MPI
3319  AssertIndexRange(component_in_block_vector, tmp_data.size());
3320  AssertDimension(requests.size(), tmp_data.size());
3321 
3322  const unsigned int mf_component = find_vector_in_mf(vec);
3323 
3324  const auto &part = get_partitioner(mf_component);
3325 
3326  if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3327  part.n_import_sm_procs() != 0)
3328  {
3329  part.export_to_ghosted_array_finish(
3330  ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3331  vec.shared_vector_data(),
3332  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3333  part.locally_owned_size(),
3334  matrix_free.get_dof_info(mf_component)
3335  .vector_partitioner->n_ghost_indices()),
3336  this->requests[component_in_block_vector]);
3337 
3338  matrix_free.release_scratch_data_non_threadsafe(
3339  tmp_data[component_in_block_vector]);
3340  tmp_data[component_in_block_vector] = nullptr;
3341  }
3342 # endif
3343  }
3344  // let vector know that ghosts are being updated and we can read from
3345  // them
3346  vec.set_ghost_state(true);
3347  }
3348 
3349 
3350 
3354  template <typename VectorType,
3355  std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3356  * = nullptr>
3357  void
3358  compress_start(const unsigned int /*component_in_block_vector*/,
3359  VectorType & /*vec*/)
3360  {}
3361 
3362 
3363 
3368  template <typename VectorType,
3369  std::enable_if_t<!has_compress_start<VectorType> &&
3370  !is_not_parallel_vector<VectorType>,
3371  VectorType> * = nullptr>
3372  void
3373  compress_start(const unsigned int component_in_block_vector,
3374  VectorType & vec)
3375  {
3376  (void)component_in_block_vector;
3377  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3378  vec.compress(::VectorOperation::add);
3379  }
3380 
3381 
3382 
3388  template <typename VectorType,
3389  std::enable_if_t<has_compress_start<VectorType> &&
3390  !has_exchange_on_subset<VectorType>,
3391  VectorType> * = nullptr>
3392  void
3393  compress_start(const unsigned int component_in_block_vector,
3394  VectorType & vec)
3395  {
3396  (void)component_in_block_vector;
3397  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3398  vec.compress_start(component_in_block_vector + channel_shift);
3399  }
3400 
3401 
3402 
3409  template <typename VectorType,
3410  std::enable_if_t<has_compress_start<VectorType> &&
3411  has_exchange_on_subset<VectorType>,
3412  VectorType> * = nullptr>
3413  void
3414  compress_start(const unsigned int component_in_block_vector,
3415  VectorType & vec)
3416  {
3417  static_assert(
3418  std::is_same<Number, typename VectorType::value_type>::value,
3419  "Type mismatch between VectorType and VectorDataExchange");
3420  (void)component_in_block_vector;
3421  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3422 
3423  if (vec.size() != 0)
3424  {
3425 # ifdef DEAL_II_WITH_MPI
3426  const unsigned int mf_component = find_vector_in_mf(vec);
3427 
3428  const auto &part = get_partitioner(mf_component);
3429 
3430  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3431  part.n_import_sm_procs() == 0)
3432  return;
3433 
3434  tmp_data[component_in_block_vector] =
3435  matrix_free.acquire_scratch_data_non_threadsafe();
3436  tmp_data[component_in_block_vector]->resize_fast(
3437  part.n_import_indices());
3438  AssertDimension(requests.size(), tmp_data.size());
3439 
3440  part.import_from_ghosted_array_start(
3442  component_in_block_vector * 2 + channel_shift,
3443  ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3444  vec.shared_vector_data(),
3445  ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3446  matrix_free.get_dof_info(mf_component)
3447  .vector_partitioner->n_ghost_indices()),
3448  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3449  part.n_import_indices()),
3450  this->requests[component_in_block_vector]);
3451 # endif
3452  }
3453  }
3454 
3455 
3456 
3461  template <
3462  typename VectorType,
3463  std::enable_if_t<!has_compress_start<VectorType>, VectorType> * = nullptr>
3464  void
3465  compress_finish(const unsigned int /*component_in_block_vector*/,
3466  VectorType & /*vec*/)
3467  {}
3468 
3469 
3470 
3476  template <typename VectorType,
3477  std::enable_if_t<has_compress_start<VectorType> &&
3478  !has_exchange_on_subset<VectorType>,
3479  VectorType> * = nullptr>
3480  void
3481  compress_finish(const unsigned int component_in_block_vector,
3482  VectorType & vec)
3483  {
3484  (void)component_in_block_vector;
3485  vec.compress_finish(::VectorOperation::add);
3486  }
3487 
3488 
3489 
3496  template <typename VectorType,
3497  std::enable_if_t<has_compress_start<VectorType> &&
3498  has_exchange_on_subset<VectorType>,
3499  VectorType> * = nullptr>
3500  void
3501  compress_finish(const unsigned int component_in_block_vector,
3502  VectorType & vec)
3503  {
3504  static_assert(
3505  std::is_same<Number, typename VectorType::value_type>::value,
3506  "Type mismatch between VectorType and VectorDataExchange");
3507  (void)component_in_block_vector;
3508  if (vec.size() != 0)
3509  {
3510 # ifdef DEAL_II_WITH_MPI
3511  AssertIndexRange(component_in_block_vector, tmp_data.size());
3512  AssertDimension(requests.size(), tmp_data.size());
3513 
3514  const unsigned int mf_component = find_vector_in_mf(vec);
3515 
3516  const auto &part = get_partitioner(mf_component);
3517 
3518  if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3519  part.n_import_sm_procs() != 0)
3520  {
3521  part.import_from_ghosted_array_finish(
3523  ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3524  vec.shared_vector_data(),
3525  ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3526  matrix_free.get_dof_info(mf_component)
3527  .vector_partitioner->n_ghost_indices()),
3529  tmp_data[component_in_block_vector]->begin(),
3530  part.n_import_indices()),
3531  this->requests[component_in_block_vector]);
3532 
3533  matrix_free.release_scratch_data_non_threadsafe(
3534  tmp_data[component_in_block_vector]);
3535  tmp_data[component_in_block_vector] = nullptr;
3536  }
3537 
3539  {
3540  const int ierr =
3541  MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3542  AssertThrowMPI(ierr);
3543  }
3544 # endif
3545  }
3546  }
3547 
3548 
3549 
3553  template <typename VectorType,
3554  std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3555  * = nullptr>
3556  void
3557  reset_ghost_values(const VectorType & /*vec*/) const
3558  {}
3559 
3560 
3561 
3566  template <typename VectorType,
3567  std::enable_if_t<!has_exchange_on_subset<VectorType> &&
3568  !is_not_parallel_vector<VectorType>,
3569  VectorType> * = nullptr>
3570  void
3571  reset_ghost_values(const VectorType &vec) const
3572  {
3573  if (ghosts_were_set == true)
3574  return;
3575 
3576  vec.zero_out_ghost_values();
3577  }
3578 
3579 
3580 
3586  template <typename VectorType,
3587  std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3588  * = nullptr>
3589  void
3590  reset_ghost_values(const VectorType &vec) const
3591  {
3592  static_assert(
3593  std::is_same<Number, typename VectorType::value_type>::value,
3594  "Type mismatch between VectorType and VectorDataExchange");
3595  if (ghosts_were_set == true)
3596  return;
3597 
3598  if (vec.size() != 0)
3599  {
3600 # ifdef DEAL_II_WITH_MPI
3601  AssertDimension(requests.size(), tmp_data.size());
3602 
3603  const unsigned int mf_component = find_vector_in_mf(vec);
3604 
3605  const auto &part = get_partitioner(mf_component);
3606 
3607  if (part.n_ghost_indices() > 0)
3608  {
3609  part.reset_ghost_values(ArrayView<Number>(
3611  .begin() +
3612  part.locally_owned_size(),
3613  matrix_free.get_dof_info(mf_component)
3614  .vector_partitioner->n_ghost_indices()));
3615  }
3616 
3617 # endif
3618  }
3619  // let vector know that it's not ghosted anymore
3620  vec.set_ghost_state(false);
3621  }
3622 
3623 
3624 
3630  template <typename VectorType,
3631  std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3632  * = nullptr>
3633  void
3634  zero_vector_region(const unsigned int range_index, VectorType &vec) const
3635  {
3636  static_assert(
3637  std::is_same<Number, typename VectorType::value_type>::value,
3638  "Type mismatch between VectorType and VectorDataExchange");
3639  if (range_index == numbers::invalid_unsigned_int)
3640  vec = Number();
3641  else
3642  {
3643  const unsigned int mf_component = find_vector_in_mf(vec, false);
3644  const internal::MatrixFreeFunctions::DoFInfo &dof_info =
3645  matrix_free.get_dof_info(mf_component);
3646  Assert(dof_info.vector_zero_range_list_index.empty() == false,
3647  ExcNotInitialized());
3648 
3649  Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3650  ExcInternalError());
3651  AssertIndexRange(range_index,
3652  dof_info.vector_zero_range_list_index.size() - 1);
3653  for (unsigned int id =
3654  dof_info.vector_zero_range_list_index[range_index];
3655  id != dof_info.vector_zero_range_list_index[range_index + 1];
3656  ++id)
3657  std::memset(vec.begin() + dof_info.vector_zero_range_list[id].first,
3658  0,
3659  (dof_info.vector_zero_range_list[id].second -
3660  dof_info.vector_zero_range_list[id].first) *
3661  sizeof(Number));
3662  }
3663  }
3664 
3665 
3666 
3672  template <typename VectorType,
3673  std::enable_if_t<!has_exchange_on_subset<VectorType>, VectorType>
3674  * = nullptr,
3675  typename VectorType::value_type * = nullptr>
3676  void
3677  zero_vector_region(const unsigned int range_index, VectorType &vec) const
3678  {
3679  if (range_index == numbers::invalid_unsigned_int || range_index == 0)
3680  vec = typename VectorType::value_type();
3681  }
3682 
3683 
3684 
3689  void
3690  zero_vector_region(const unsigned int, ...) const
3691  {
3692  Assert(false,
3693  ExcNotImplemented("Zeroing is only implemented for vector types "
3694  "which provide VectorType::value_type"));
3695  }
3696 
3697 
3698 
3699  const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
3700  const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3701  DataAccessOnFaces vector_face_access;
3702  bool ghosts_were_set;
3703 # ifdef DEAL_II_WITH_MPI
3704  std::vector<AlignedVector<Number> *> tmp_data;
3705  std::vector<std::vector<MPI_Request>> requests;
3706 # endif
3707  }; // VectorDataExchange
3708 
3709  template <typename VectorStruct>
3710  unsigned int
3711  n_components(const VectorStruct &vec);
3712 
3713  template <typename VectorStruct>
3714  unsigned int
3715  n_components_block(const VectorStruct &vec,
3716  std::integral_constant<bool, true>)
3717  {
3718  unsigned int components = 0;
3719  for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
3720  components += n_components(vec.block(bl));
3721  return components;
3722  }
3723 
3724  template <typename VectorStruct>
3725  unsigned int
3726  n_components_block(const VectorStruct &, std::integral_constant<bool, false>)
3727  {
3728  return 1;
3729  }
3730 
3731  template <typename VectorStruct>
3732  unsigned int
3733  n_components(const VectorStruct &vec)
3734  {
3735  return n_components_block(
3736  vec, std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3737  }
3738 
3739  template <typename VectorStruct>
3740  inline unsigned int
3741  n_components(const std::vector<VectorStruct> &vec)
3742  {
3743  unsigned int components = 0;
3744  for (unsigned int comp = 0; comp < vec.size(); ++comp)
3745  components += n_components_block(
3746  vec[comp],
3747  std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3748  return components;
3749  }
3750 
3751  template <typename VectorStruct>
3752  inline unsigned int
3753  n_components(const std::vector<VectorStruct *> &vec)
3754  {
3755  unsigned int components = 0;
3756  for (unsigned int comp = 0; comp < vec.size(); ++comp)
3757  components += n_components_block(
3758  *vec[comp],
3759  std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3760  return components;
3761  }
3762 
3763 
3764 
3765  // A helper function to identify block vectors with many components where we
3766  // should not try to overlap computations and communication because there
3767  // would be too many outstanding communication requests.
3768 
3769  // default value for vectors that do not have communication_block_size
3770  template <typename VectorStruct,
3771  std::enable_if_t<!has_communication_block_size<VectorStruct>,
3772  VectorStruct> * = nullptr>
3773  constexpr unsigned int
3774  get_communication_block_size(const VectorStruct &)
3775  {
3777  }
3778 
3779 
3780 
3781  template <typename VectorStruct,
3782  std::enable_if_t<has_communication_block_size<VectorStruct>,
3783  VectorStruct> * = nullptr>
3784  constexpr unsigned int
3785  get_communication_block_size(const VectorStruct &)
3786  {
3787  return VectorStruct::communication_block_size;
3788  }
3789 
3790 
3791 
3792  // --------------------------------------------------------------------------
3793  // below we have wrappers to distinguish between block and non-block vectors.
3794  // --------------------------------------------------------------------------
3795 
3796  //
3797  // update_ghost_values_start
3798  //
3799 
3800  // update_ghost_values for block vectors
3801  template <int dim,
3802  typename VectorStruct,
3803  typename Number,
3804  typename VectorizedArrayType,
3805  std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
3806  * = nullptr>
3807  void
3808  update_ghost_values_start(
3809  const VectorStruct & vec,
3810  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3811  const unsigned int channel = 0)
3812  {
3813  if (get_communication_block_size(vec) < vec.n_blocks())
3814  {
3815  // don't forget to set ghosts_were_set, that otherwise happens
3816  // inside VectorDataExchange::update_ghost_values_start()
3817  exchanger.ghosts_were_set = vec.has_ghost_elements();
3818  vec.update_ghost_values();
3819  }
3820  else
3821  {
3822  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
3823  update_ghost_values_start(vec.block(i), exchanger, channel + i);
3824  }
3825  }
3826 
3827 
3828 
3829  // update_ghost_values for non-block vectors
3830  template <int dim,
3831  typename VectorStruct,
3832  typename Number,
3833  typename VectorizedArrayType,
3834  std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
3835  * = nullptr>
3836  void
3837  update_ghost_values_start(
3838  const VectorStruct & vec,
3839  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3840  const unsigned int channel = 0)
3841  {
3842  exchanger.update_ghost_values_start(channel, vec);
3843  }
3844 
3845 
3846 
3847  // update_ghost_values_start() for vector of vectors
3848  template <int dim,
3849  typename VectorStruct,
3850  typename Number,
3851  typename VectorizedArrayType>
3852  inline void
3853  update_ghost_values_start(
3854  const std::vector<VectorStruct> & vec,
3855  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3856  {
3857  unsigned int component_index = 0;
3858  for (unsigned int comp = 0; comp < vec.size(); ++comp)
3859  {
3860  update_ghost_values_start(vec[comp], exchanger, component_index);
3861  component_index += n_components(vec[comp]);
3862  }
3863  }
3864 
3865 
3866 
3867  // update_ghost_values_start() for vector of pointers to vectors
3868  template <int dim,
3869  typename VectorStruct,
3870  typename Number,
3871  typename VectorizedArrayType>
3872  inline void
3873  update_ghost_values_start(
3874  const std::vector<VectorStruct *> & vec,
3875  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3876  {
3877  unsigned int component_index = 0;
3878  for (unsigned int comp = 0; comp < vec.size(); ++comp)
3879  {
3880  update_ghost_values_start(*vec[comp], exchanger, component_index);
3881  component_index += n_components(*vec[comp]);
3882  }
3883  }
3884 
3885 
3886 
3887  //
3888  // update_ghost_values_finish
3889  //
3890 
3891  // for block vectors
3892  template <int dim,
3893  typename VectorStruct,
3894  typename Number,
3895  typename VectorizedArrayType,
3896  std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
3897  * = nullptr>
3898  void
3899  update_ghost_values_finish(
3900  const VectorStruct & vec,
3901  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3902  const unsigned int channel = 0)
3903  {
3904  if (get_communication_block_size(vec) < vec.n_blocks())
3905  {
3906  // do nothing, everything has already been completed in the _start()
3907  // call
3908  }
3909  else
3910  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
3911  update_ghost_values_finish(vec.block(i), exchanger, channel + i);
3912  }
3913 
3914 
3915 
3916  // for non-block vectors
3917  template <int dim,
3918  typename VectorStruct,
3919  typename Number,
3920  typename VectorizedArrayType,
3921  std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
3922  * = nullptr>
3923  void
3924  update_ghost_values_finish(
3925  const VectorStruct & vec,
3926  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3927  const unsigned int channel = 0)
3928  {
3929  exchanger.update_ghost_values_finish(channel, vec);
3930  }
3931 
3932 
3933 
3934  // for vector of vectors
3935  template <int dim,
3936  typename VectorStruct,
3937  typename Number,
3938  typename VectorizedArrayType>
3939  inline void
3940  update_ghost_values_finish(
3941  const std::vector<VectorStruct> & vec,
3942  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3943  {
3944  unsigned int component_index = 0;
3945  for (unsigned int comp = 0; comp < vec.size(); ++comp)
3946  {
3947  update_ghost_values_finish(vec[comp], exchanger, component_index);
3948  component_index += n_components(vec[comp]);
3949  }
3950  }
3951 
3952 
3953 
3954  // for vector of pointers to vectors
3955  template <int dim,
3956  typename VectorStruct,
3957  typename Number,
3958  typename VectorizedArrayType>
3959  inline void
3960  update_ghost_values_finish(
3961  const std::vector<VectorStruct *> & vec,
3962  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3963  {
3964  unsigned int component_index = 0;
3965  for (unsigned int comp = 0; comp < vec.size(); ++comp)
3966  {
3967  update_ghost_values_finish(*vec[comp], exchanger, component_index);
3968  component_index += n_components(*vec[comp]);
3969  }
3970  }
3971 
3972 
3973 
3974  //
3975  // compress_start
3976  //
3977 
3978  // for block vectors
3979  template <int dim,
3980  typename VectorStruct,
3981  typename Number,
3982  typename VectorizedArrayType,
3983  std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
3984  * = nullptr>
3985  inline void
3986  compress_start(
3987  VectorStruct & vec,
3988  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3989  const unsigned int channel = 0)
3990  {
3991  if (get_communication_block_size(vec) < vec.n_blocks())
3992  vec.compress(::VectorOperation::add);
3993  else
3994  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
3995  compress_start(vec.block(i), exchanger, channel + i);
3996  }
3997 
3998 
3999 
4000  // for non-block vectors
4001  template <int dim,
4002  typename VectorStruct,
4003  typename Number,
4004  typename VectorizedArrayType,
4005  std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4006  * = nullptr>
4007  inline void
4008  compress_start(
4009  VectorStruct & vec,
4010  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4011  const unsigned int channel = 0)
4012  {
4013  exchanger.compress_start(channel, vec);
4014  }
4015 
4016 
4017 
4018  // for std::vector of vectors
4019  template <int dim,
4020  typename VectorStruct,
4021  typename Number,
4022  typename VectorizedArrayType>
4023  inline void
4024  compress_start(
4025  std::vector<VectorStruct> & vec,
4026  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4027  {
4028  unsigned int component_index = 0;
4029  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4030  {
4031  compress_start(vec[comp], exchanger, component_index);
4032  component_index += n_components(vec[comp]);
4033  }
4034  }
4035 
4036 
4037 
4038  // for std::vector of pointer to vectors
4039  template <int dim,
4040  typename VectorStruct,
4041  typename Number,
4042  typename VectorizedArrayType>
4043  inline void
4044  compress_start(
4045  std::vector<VectorStruct *> & vec,
4046  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4047  {
4048  unsigned int component_index = 0;
4049  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4050  {
4051  compress_start(*vec[comp], exchanger, component_index);
4052  component_index += n_components(*vec[comp]);
4053  }
4054  }
4055 
4056 
4057 
4058  //
4059  // compress_finish
4060  //
4061 
4062  // for block vectors
4063  template <int dim,
4064  typename VectorStruct,
4065  typename Number,
4066  typename VectorizedArrayType,
4067  std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4068  * = nullptr>
4069  inline void
4070  compress_finish(
4071  VectorStruct & vec,
4072  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4073  const unsigned int channel = 0)
4074  {
4075  if (get_communication_block_size(vec) < vec.n_blocks())
4076  {
4077  // do nothing, everything has already been completed in the _start()
4078  // call
4079  }
4080  else
4081  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4082  compress_finish(vec.block(i), exchanger, channel + i);
4083  }
4084 
4085 
4086 
4087  // for non-block vectors
4088  template <int dim,
4089  typename VectorStruct,
4090  typename Number,
4091  typename VectorizedArrayType,
4092  std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4093  * = nullptr>
4094  inline void
4095  compress_finish(
4096  VectorStruct & vec,
4097  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4098  const unsigned int channel = 0)
4099  {
4100  exchanger.compress_finish(channel, vec);
4101  }
4102 
4103 
4104 
4105  // for std::vector of vectors
4106  template <int dim,
4107  typename VectorStruct,
4108  typename Number,
4109  typename VectorizedArrayType>
4110  inline void
4111  compress_finish(
4112  std::vector<VectorStruct> & vec,
4113  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4114  {
4115  unsigned int component_index = 0;
4116  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4117  {
4118  compress_finish(vec[comp], exchanger, component_index);
4119  component_index += n_components(vec[comp]);
4120  }
4121  }
4122 
4123 
4124 
4125  // for std::vector of pointer to vectors
4126  template <int dim,
4127  typename VectorStruct,
4128  typename Number,
4129  typename VectorizedArrayType>
4130  inline void
4131  compress_finish(
4132  std::vector<VectorStruct *> & vec,
4133  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4134  {
4135  unsigned int component_index = 0;
4136  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4137  {
4138  compress_finish(*vec[comp], exchanger, component_index);
4139  component_index += n_components(*vec[comp]);
4140  }
4141  }
4142 
4143 
4144 
4145  //
4146  // reset_ghost_values:
4147  //
4148  // if the input vector did not have ghosts imported, clear them here again
4149  // in order to avoid subsequent operations e.g. in linear solvers to work
4150  // with ghosts all the time
4151  //
4152 
4153  // for block vectors
4154  template <int dim,
4155  typename VectorStruct,
4156  typename Number,
4157  typename VectorizedArrayType,
4158  std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4159  * = nullptr>
4160  inline void
4161  reset_ghost_values(
4162  const VectorStruct & vec,
4163  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4164  {
4165  // return immediately if there is nothing to do.
4166  if (exchanger.ghosts_were_set == true)
4167  return;
4168 
4169  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4170  reset_ghost_values(vec.block(i), exchanger);
4171  }
4172 
4173 
4174 
4175  // for non-block vectors
4176  template <int dim,
4177  typename VectorStruct,
4178  typename Number,
4179  typename VectorizedArrayType,
4180  std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4181  * = nullptr>
4182  inline void
4183  reset_ghost_values(
4184  const VectorStruct & vec,
4185  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4186  {
4187  exchanger.reset_ghost_values(vec);
4188  }
4189 
4190 
4191 
4192  // for std::vector of vectors
4193  template <int dim,
4194  typename VectorStruct,
4195  typename Number,
4196  typename VectorizedArrayType>
4197  inline void
4198  reset_ghost_values(
4199  const std::vector<VectorStruct> & vec,
4200  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4201  {
4202  // return immediately if there is nothing to do.
4203  if (exchanger.ghosts_were_set == true)
4204  return;
4205 
4206  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4207  reset_ghost_values(vec[comp], exchanger);
4208  }
4209 
4210 
4211 
4212  // for std::vector of pointer to vectors
4213  template <int dim,
4214  typename VectorStruct,
4215  typename Number,
4216  typename VectorizedArrayType>
4217  inline void
4218  reset_ghost_values(
4219  const std::vector<VectorStruct *> & vec,
4220  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4221  {
4222  // return immediately if there is nothing to do.
4223  if (exchanger.ghosts_were_set == true)
4224  return;
4225 
4226  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4227  reset_ghost_values(*vec[comp], exchanger);
4228  }
4229 
4230 
4231 
4232  //
4233  // zero_vector_region
4234  //
4235 
4236  // for block vectors
4237  template <int dim,
4238  typename VectorStruct,
4239  typename Number,
4240  typename VectorizedArrayType,
4241  std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4242  * = nullptr>
4243  inline void
4244  zero_vector_region(
4245  const unsigned int range_index,
4246  VectorStruct & vec,
4247  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4248  {
4249  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4250  exchanger.zero_vector_region(range_index, vec.block(i));
4251  }
4252 
4253 
4254 
4255  // for non-block vectors
4256  template <int dim,
4257  typename VectorStruct,
4258  typename Number,
4259  typename VectorizedArrayType,
4260  std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4261  * = nullptr>
4262  inline void
4263  zero_vector_region(
4264  const unsigned int range_index,
4265  VectorStruct & vec,
4266  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4267  {
4268  exchanger.zero_vector_region(range_index, vec);
4269  }
4270 
4271 
4272 
4273  // for std::vector of vectors
4274  template <int dim,
4275  typename VectorStruct,
4276  typename Number,
4277  typename VectorizedArrayType>
4278  inline void
4279  zero_vector_region(
4280  const unsigned int range_index,
4281  std::vector<VectorStruct> & vec,
4282  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4283  {
4284  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4285  zero_vector_region(range_index, vec[comp], exchanger);
4286  }
4287 
4288 
4289 
4290  // for std::vector of pointers to vectors
4291  template <int dim,
4292  typename VectorStruct,
4293  typename Number,
4294  typename VectorizedArrayType>
4295  inline void
4296  zero_vector_region(
4297  const unsigned int range_index,
4298  std::vector<VectorStruct *> & vec,
4299  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4300  {
4301  for (unsigned int comp = 0; comp < vec.size(); ++comp)
4302  zero_vector_region(range_index, *vec[comp], exchanger);
4303  }
4304 
4305 
4306 
4307  // Apply a unit matrix operation to constrained DoFs: Default cases where we
4308  // cannot detect a LinearAlgebra::distributed::Vector, we do not do
4309  // anything, else we apply the constraints as a unit operation
4310  template <typename VectorStruct1, typename VectorStruct2>
4311  inline void
4312  apply_operation_to_constrained_dofs(const std::vector<unsigned int> &,
4313  const VectorStruct1 &,
4314  VectorStruct2 &)
4315  {}
4316 
4317  template <typename Number>
4318  inline void
4319  apply_operation_to_constrained_dofs(
4320  const std::vector<unsigned int> & constrained_dofs,
4323  {
4324  for (const unsigned int i : constrained_dofs)
4325  dst.local_element(i) = src.local_element(i);
4326  }
4327 
4328 
4329  namespace MatrixFreeFunctions
4330  {
4331  // struct to select between a const interface and a non-const interface
4332  // for MFWorker
4333  template <typename, typename, typename, typename, bool>
4334  struct InterfaceSelector
4335  {};
4336 
4337  // Version for constant functions
4338  template <typename MF,
4339  typename InVector,
4340  typename OutVector,
4341  typename Container>
4342  struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4343  {
4344  using function_type = void (Container::*)(
4345  const MF &,
4346  OutVector &,
4347  const InVector &,
4348  const std::pair<unsigned int, unsigned int> &) const;
4349  };
4350 
4351  // Version for non-constant functions
4352  template <typename MF,
4353  typename InVector,
4354  typename OutVector,
4355  typename Container>
4356  struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4357  {
4358  using function_type =
4359  void (Container::*)(const MF &,
4360  OutVector &,
4361  const InVector &,
4362  const std::pair<unsigned int, unsigned int> &);
4363  };
4364  } // namespace MatrixFreeFunctions
4365 
4366 
4367 
4368  // A implementation class for the worker object that runs the various
4369  // operations we want to perform during the matrix-free loop
4370  template <typename MF,
4371  typename InVector,
4372  typename OutVector,
4373  typename Container,
4374  bool is_constant>
4375  class MFWorker : public MFWorkerInterface
4376  {
4377  public:
4378  // An alias to make the arguments further down more readable
4379  using function_type = typename MatrixFreeFunctions::
4380  InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4381  function_type;
4382 
4383  // constructor, binds all the arguments to this class
4384  MFWorker(const MF & matrix_free,
4385  const InVector & src,
4386  OutVector & dst,
4387  const bool zero_dst_vector_setting,
4388  const Container & container,
4389  function_type cell_function,
4390  function_type face_function,
4391  function_type boundary_function,
4392  const typename MF::DataAccessOnFaces src_vector_face_access =
4394  const typename MF::DataAccessOnFaces dst_vector_face_access =
4396  const std::function<void(const unsigned int, const unsigned int)>
4397  &operation_before_loop = {},
4398  const std::function<void(const unsigned int, const unsigned int)>
4399  & operation_after_loop = {},
4400  const unsigned int dof_handler_index_pre_post = 0)
4401  : matrix_free(matrix_free)
4402  , container(const_cast<Container &>(container))
4403  , cell_function(cell_function)
4404  , face_function(face_function)
4405  , boundary_function(boundary_function)
4406  , src(src)
4407  , dst(dst)
4408  , src_data_exchanger(matrix_free,
4409  src_vector_face_access,
4410  n_components(src))
4411  , dst_data_exchanger(matrix_free,
4412  dst_vector_face_access,
4413  n_components(dst))
4414  , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4415  , zero_dst_vector_setting(zero_dst_vector_setting &&
4416  !src_and_dst_are_same)
4417  , operation_before_loop(operation_before_loop)
4418  , operation_after_loop(operation_after_loop)
4419  , dof_handler_index_pre_post(dof_handler_index_pre_post)
4420  {}
4421 
4422  // Runs the cell work. If no function is given, nothing is done
4423  virtual void
4424  cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4425  {
4426  if (cell_function != nullptr && cell_range.second > cell_range.first)
4427  for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4428  {
4429  const auto cell_subrange =
4430  matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4431 
4432  if (cell_subrange.second <= cell_subrange.first)
4433  continue;
4434 
4435  (container.*
4436  cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4437  }
4438  }
4439 
4440  virtual void
4441  cell(const unsigned int range_index) override
4442  {
4443  process_range(cell_function,
4444  matrix_free.get_task_info().cell_partition_data_hp_ptr,
4445  matrix_free.get_task_info().cell_partition_data_hp,
4446  range_index);
4447  }
4448 
4449  virtual void
4450  face(const unsigned int range_index) override
4451  {
4452  process_range(face_function,
4453  matrix_free.get_task_info().face_partition_data_hp_ptr,
4454  matrix_free.get_task_info().face_partition_data_hp,
4455  range_index);
4456  }
4457 
4458  virtual void
4459  boundary(const unsigned int range_index) override
4460  {
4461  process_range(boundary_function,
4462  matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4463  matrix_free.get_task_info().boundary_partition_data_hp,
4464  range_index);
4465  }
4466 
4467  private:
4468  void
4469  process_range(const function_type & fu,
4470  const std::vector<unsigned int> &ptr,
4471  const std::vector<unsigned int> &data,
4472  const unsigned int range_index)
4473  {
4474  if (fu == nullptr)
4475  return;
4476 
4477  AssertIndexRange(range_index + 1, ptr.size());
4478  for (unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4479  {
4480  AssertIndexRange(2 * i + 1, data.size());
4481  (container.*fu)(matrix_free,
4482  this->dst,
4483  this->src,
4484  std::make_pair(data[2 * i], data[2 * i + 1]));
4485  }
4486  }
4487 
4488  public:
4489  // Starts the communication for the update ghost values operation. We
4490  // cannot call this update if ghost and destination are the same because
4491  // that would introduce spurious entries in the destination (there is also
4492  // the problem that reading from a vector that we also write to is usually
4493  // not intended in case there is overlap, but this is up to the
4494  // application code to decide and we cannot catch this case here).
4495  virtual void
4496  vector_update_ghosts_start() override
4497  {
4498  if (!src_and_dst_are_same)
4499  internal::update_ghost_values_start(src, src_data_exchanger);
4500  }
4501 
4502  // Finishes the communication for the update ghost values operation
4503  virtual void
4504  vector_update_ghosts_finish() override
4505  {
4506  if (!src_and_dst_are_same)
4507  internal::update_ghost_values_finish(src, src_data_exchanger);
4508  }
4509 
4510  // Starts the communication for the vector compress operation
4511  virtual void
4512  vector_compress_start() override
4513  {
4514  internal::compress_start(dst, dst_data_exchanger);
4515  }
4516 
4517  // Finishes the communication for the vector compress operation
4518  virtual void
4519  vector_compress_finish() override
4520  {
4521  internal::compress_finish(dst, dst_data_exchanger);
4522  if (!src_and_dst_are_same)
4523  internal::reset_ghost_values(src, src_data_exchanger);
4524  }
4525 
4526  // Zeros the given input vector
4527  virtual void
4528  zero_dst_vector_range(const unsigned int range_index) override
4529  {
4530  if (zero_dst_vector_setting)
4531  internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4532  }
4533 
4534  virtual void
4535  cell_loop_pre_range(const unsigned int range_index) override
4536  {
4537  if (operation_before_loop)
4538  {
4539  const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4540  matrix_free.get_dof_info(dof_handler_index_pre_post);
4541  if (range_index == numbers::invalid_unsigned_int)
4542  {
4543  // Case with threaded loop -> currently no overlap implemented
4545  0U,
4546  dof_info.vector_partitioner->locally_owned_size(),
4547  operation_before_loop,
4549  }
4550  else
4551  {
4552  AssertIndexRange(range_index,
4553  dof_info.cell_loop_pre_list_index.size() - 1);
4554  for (unsigned int id =
4555  dof_info.cell_loop_pre_list_index[range_index];
4556  id != dof_info.cell_loop_pre_list_index[range_index + 1];
4557  ++id)
4558  operation_before_loop(dof_info.cell_loop_pre_list[id].first,
4559  dof_info.cell_loop_pre_list[id].second);
4560  }
4561  }
4562  }
4563 
4564  virtual void
4565  cell_loop_post_range(const unsigned int range_index) override
4566  {
4567  if (operation_after_loop)
4568  {
4569  // Run unit matrix operation on constrained dofs if we are at the
4570  // last range
4571  const std::vector<unsigned int> &partition_row_index =
4572  matrix_free.get_task_info().partition_row_index;
4573  if (range_index ==
4574  partition_row_index[partition_row_index.size() - 2] - 1)
4575  apply_operation_to_constrained_dofs(
4576  matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
4577  src,
4578  dst);
4579 
4580  const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4581  matrix_free.get_dof_info(dof_handler_index_pre_post);
4582  if (range_index == numbers::invalid_unsigned_int)
4583  {
4584  // Case with threaded loop -> currently no overlap implemented
4586  0U,
4587  dof_info.vector_partitioner->locally_owned_size(),
4588  operation_after_loop,
4590  }
4591  else
4592  {
4593  AssertIndexRange(range_index,
4594  dof_info.cell_loop_post_list_index.size() - 1);
4595  for (unsigned int id =
4596  dof_info.cell_loop_post_list_index[range_index];
4597  id != dof_info.cell_loop_post_list_index[range_index + 1];
4598  ++id)
4599  operation_after_loop(dof_info.cell_loop_post_list[id].first,
4600  dof_info.cell_loop_post_list[id].second);
4601  }
4602  }
4603  }
4604 
4605  private:
4606  const MF & matrix_free;
4607  Container & container;
4608  function_type cell_function;
4609  function_type face_function;
4610  function_type boundary_function;
4611 
4612  const InVector &src;
4613  OutVector & dst;
4614  VectorDataExchange<MF::dimension,
4615  typename MF::value_type,
4616  typename MF::vectorized_value_type>
4617  src_data_exchanger;
4618  VectorDataExchange<MF::dimension,
4619  typename MF::value_type,
4620  typename MF::vectorized_value_type>
4621  dst_data_exchanger;
4622  const bool src_and_dst_are_same;
4623  const bool zero_dst_vector_setting;
4624  const std::function<void(const unsigned int, const unsigned int)>
4625  operation_before_loop;
4626  const std::function<void(const unsigned int, const unsigned int)>
4627  operation_after_loop;
4628  const unsigned int dof_handler_index_pre_post;
4629  };
4630 
4631 
4632 
4637  template <class MF, typename InVector, typename OutVector>
4638  struct MFClassWrapper
4639  {
4640  using function_type =
4641  std::function<void(const MF &,
4642  OutVector &,
4643  const InVector &,
4644  const std::pair<unsigned int, unsigned int> &)>;
4645 
4646  MFClassWrapper(const function_type cell,
4647  const function_type face,
4648  const function_type boundary)
4649  : cell(cell)
4650  , face(face)
4651  , boundary(boundary)
4652  {}
4653 
4654  void
4655  cell_integrator(const MF & mf,
4656  OutVector & dst,
4657  const InVector & src,
4658  const std::pair<unsigned int, unsigned int> &range) const
4659  {
4660  if (cell)
4661  cell(mf, dst, src, range);
4662  }
4663 
4664  void
4665  face_integrator(const MF & mf,
4666  OutVector & dst,
4667  const InVector & src,
4668  const std::pair<unsigned int, unsigned int> &range) const
4669  {
4670  if (face)
4671  face(mf, dst, src, range);
4672  }
4673 
4674  void
4675  boundary_integrator(
4676  const MF & mf,
4677  OutVector & dst,
4678  const InVector & src,
4679  const std::pair<unsigned int, unsigned int> &range) const
4680  {
4681  if (boundary)
4682  boundary(mf, dst, src, range);
4683  }
4684 
4685  const function_type cell;
4686  const function_type face;
4687  const function_type boundary;
4688  };
4689 
4690 } // end of namespace internal
4691 
4692 
4693 
4694 template <int dim, typename Number, typename VectorizedArrayType>
4695 template <typename OutVector, typename InVector>
4696 inline void
4698  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4699  OutVector &,
4700  const InVector &,
4701  const std::pair<unsigned int, unsigned int> &)>
4702  & cell_operation,
4703  OutVector & dst,
4704  const InVector &src,
4705  const bool zero_dst_vector) const
4706 {
4707  using Wrapper =
4708  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4709  InVector,
4710  OutVector>;
4711  Wrapper wrap(cell_operation, nullptr, nullptr);
4712  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4713  InVector,
4714  OutVector,
4715  Wrapper,
4716  true>
4717  worker(*this,
4718  src,
4719  dst,
4720  zero_dst_vector,
4721  wrap,
4722  &Wrapper::cell_integrator,
4723  &Wrapper::face_integrator,
4724  &Wrapper::boundary_integrator);
4725 
4726  task_info.loop(worker);
4727 }
4728 
4729 
4730 
4731 template <int dim, typename Number, typename VectorizedArrayType>
4732 template <typename OutVector, typename InVector>
4733 inline void
4735  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4736  OutVector &,
4737  const InVector &,
4738  const std::pair<unsigned int, unsigned int> &)>
4739  & cell_operation,
4740  OutVector & dst,
4741  const InVector &src,
4742  const std::function<void(const unsigned int, const unsigned int)>
4743  &operation_before_loop,
4744  const std::function<void(const unsigned int, const unsigned int)>
4745  & operation_after_loop,
4746  const unsigned int dof_handler_index_pre_post) const
4747 {
4748  using Wrapper =
4749  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4750  InVector,
4751  OutVector>;
4752  Wrapper wrap(cell_operation, nullptr, nullptr);
4753  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4754  InVector,
4755  OutVector,
4756  Wrapper,
4757  true>
4758  worker(*this,
4759  src,
4760  dst,
4761  false,
4762  wrap,
4763  &Wrapper::cell_integrator,
4764  &Wrapper::face_integrator,
4765  &Wrapper::boundary_integrator,
4768  operation_before_loop,
4769  operation_after_loop,
4770  dof_handler_index_pre_post);
4771 
4772  task_info.loop(worker);
4773 }
4774 
4775 
4776 
4777 template <int dim, typename Number, typename VectorizedArrayType>
4778 template <typename OutVector, typename InVector>
4779 inline void
4781  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4782  OutVector &,
4783  const InVector &,
4784  const std::pair<unsigned int, unsigned int> &)>
4785  &cell_operation,
4786  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4787  OutVector &,
4788  const InVector &,
4789  const std::pair<unsigned int, unsigned int> &)>
4790  &face_operation,
4791  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4792  OutVector &,
4793  const InVector &,
4794  const std::pair<unsigned int, unsigned int> &)>
4795  & boundary_operation,
4796  OutVector & dst,
4797  const InVector & src,
4798  const bool zero_dst_vector,
4799  const DataAccessOnFaces dst_vector_face_access,
4800  const DataAccessOnFaces src_vector_face_access) const
4801 {
4802  using Wrapper =
4803  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4804  InVector,
4805  OutVector>;
4806  Wrapper wrap(cell_operation, face_operation, boundary_operation);
4807  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4808  InVector,
4809  OutVector,
4810  Wrapper,
4811  true>
4812  worker(*this,
4813  src,
4814  dst,
4815  zero_dst_vector,
4816  wrap,
4817  &Wrapper::cell_integrator,
4818  &Wrapper::face_integrator,
4819  &Wrapper::boundary_integrator,
4820  src_vector_face_access,
4821  dst_vector_face_access);
4822 
4823  task_info.loop(worker);
4824 }
4825 
4826 
4827 
4828 template <int dim, typename Number, typename VectorizedArrayType>
4829 template <typename CLASS, typename OutVector, typename InVector>
4830 inline void
4832  void (CLASS::*function_pointer)(
4834  OutVector &,
4835  const InVector &,
4836  const std::pair<unsigned int, unsigned int> &) const,
4837  const CLASS * owning_class,
4838  OutVector & dst,
4839  const InVector &src,
4840  const bool zero_dst_vector) const
4841 {
4842  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4843  InVector,
4844  OutVector,
4845  CLASS,
4846  true>
4847  worker(*this,
4848  src,
4849  dst,
4850  zero_dst_vector,
4851  *owning_class,
4852  function_pointer,
4853  nullptr,
4854  nullptr);
4855  task_info.loop(worker);
4856 }
4857 
4858 
4859 
4860 template <int dim, typename Number, typename VectorizedArrayType>
4861 template <typename CLASS, typename OutVector, typename InVector>
4862 inline void
4864  void (CLASS::*function_pointer)(
4866  OutVector &,
4867  const InVector &,
4868  const std::pair<unsigned int, unsigned int> &) const,
4869  const CLASS * owning_class,
4870  OutVector & dst,
4871  const InVector &src,
4872  const std::function<void(const unsigned int, const unsigned int)>
4873  &operation_before_loop,
4874  const std::function<void(const unsigned int, const unsigned int)>
4875  & operation_after_loop,
4876  const unsigned int dof_handler_index_pre_post) const
4877 {
4878  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4879  InVector,
4880  OutVector,
4881  CLASS,
4882  true>
4883  worker(*this,
4884  src,
4885  dst,
4886  false,
4887  *owning_class,
4888  function_pointer,
4889  nullptr,
4890  nullptr,
4893  operation_before_loop,
4894  operation_after_loop,
4895  dof_handler_index_pre_post);
4896  task_info.loop(worker);
4897 }
4898 
4899 
4900 
4901 template <int dim, typename Number, typename VectorizedArrayType>
4902 template <typename CLASS, typename OutVector, typename InVector>
4903 inline void
4905  void (CLASS::*cell_operation)(
4907  OutVector &,
4908  const InVector &,
4909  const std::pair<unsigned int, unsigned int> &) const,
4910  void (CLASS::*face_operation)(
4912  OutVector &,
4913  const InVector &,
4914  const std::pair<unsigned int, unsigned int> &) const,
4915  void (CLASS::*boundary_operation)(
4917  OutVector &,
4918  const InVector &,
4919  const std::pair<unsigned int, unsigned int> &) const,
4920  const CLASS * owning_class,
4921  OutVector & dst,
4922  const InVector & src,
4923  const bool zero_dst_vector,
4924  const DataAccessOnFaces dst_vector_face_access,
4925  const DataAccessOnFaces src_vector_face_access) const
4926 {
4927  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4928  InVector,
4929  OutVector,
4930  CLASS,
4931  true>
4932  worker(*this,
4933  src,
4934  dst,
4935  zero_dst_vector,
4936  *owning_class,
4937  cell_operation,
4938  face_operation,
4939  boundary_operation,
4940  src_vector_face_access,
4941  dst_vector_face_access);
4942  task_info.loop(worker);
4943 }
4944 
4945 
4946 
4947 template <int dim, typename Number, typename VectorizedArrayType>
4948 template <typename CLASS, typename OutVector, typename InVector>
4949 inline void
4951  void (CLASS::*function_pointer)(
4953  OutVector &,
4954  const InVector &,
4955  const std::pair<unsigned int, unsigned int> &),
4956  CLASS * owning_class,
4957  OutVector & dst,
4958  const InVector &src,
4959  const bool zero_dst_vector) const
4960 {
4961  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4962  InVector,
4963  OutVector,
4964  CLASS,
4965  false>
4966  worker(*this,
4967  src,
4968  dst,
4969  zero_dst_vector,
4970  *owning_class,
4971  function_pointer,
4972  nullptr,
4973  nullptr);
4974  task_info.loop(worker);
4975 }
4976 
4977 
4978 
4979 template <int dim, typename Number, typename VectorizedArrayType>
4980 template <typename CLASS, typename OutVector, typename InVector>
4981 inline void
4983  void (CLASS::*function_pointer)(
4985  OutVector &,
4986  const InVector &,
4987  const std::pair<unsigned int, unsigned int> &),
4988  CLASS * owning_class,
4989  OutVector & dst,
4990  const InVector &src,
4991  const std::function<void(const unsigned int, const unsigned int)>
4992  &operation_before_loop,
4993  const std::function<void(const unsigned int, const unsigned int)>
4994  & operation_after_loop,
4995  const unsigned int dof_handler_index_pre_post) const
4996 {
4997  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4998  InVector,
4999  OutVector,
5000  CLASS,
5001  false>
5002  worker(*this,
5003  src,
5004  dst,
5005  false,
5006  *owning_class,
5007  function_pointer,
5008  nullptr,
5009  nullptr,
5012  operation_before_loop,
5013  operation_after_loop,
5014  dof_handler_index_pre_post);
5015  task_info.loop(worker);
5016 }
5017 
5018 
5019 
5020 template <int dim, typename Number, typename VectorizedArrayType>
5021 template <typename CLASS, typename OutVector, typename InVector>
5022 inline void
5024  void (CLASS::*cell_operation)(
5026  OutVector &,
5027  const InVector &,
5028  const std::pair<unsigned int, unsigned int> &),
5029  void (CLASS::*face_operation)(
5031  OutVector &,
5032  const InVector &,
5033  const std::pair<unsigned int, unsigned int> &),
5034  void (CLASS::*boundary_operation)(
5036  OutVector &,
5037  const InVector &,
5038  const std::pair<unsigned int, unsigned int> &),
5039  CLASS * owning_class,
5040  OutVector & dst,
5041  const InVector & src,
5042  const bool zero_dst_vector,
5043  const DataAccessOnFaces dst_vector_face_access,
5044  const DataAccessOnFaces src_vector_face_access) const
5045 {
5046  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5047  InVector,
5048  OutVector,
5049  CLASS,
5050  false>
5051  worker(*this,
5052  src,
5053  dst,
5054  zero_dst_vector,
5055  *owning_class,
5056  cell_operation,
5057  face_operation,
5058  boundary_operation,
5059  src_vector_face_access,
5060  dst_vector_face_access);
5061  task_info.loop(worker);
5062 }
5063 
5064 
5065 
5066 template <int dim, typename Number, typename VectorizedArrayType>
5067 template <typename CLASS, typename OutVector, typename InVector>
5068 inline void
5070  void (CLASS::*function_pointer)(
5072  OutVector &,
5073  const InVector &,
5074  const std::pair<unsigned int, unsigned int> &) const,
5075  const CLASS * owning_class,
5076  OutVector & dst,
5077  const InVector & src,
5078  const bool zero_dst_vector,
5079  const DataAccessOnFaces src_vector_face_access) const
5080 {
5081  auto src_vector_face_access_temp = src_vector_face_access;
5082  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5083  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5084  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5085  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5086 
5087  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5088  InVector,
5089  OutVector,
5090  CLASS,
5091  true>
5092  worker(*this,
5093  src,
5094  dst,
5095  zero_dst_vector,
5096  *owning_class,
5097  function_pointer,
5098  nullptr,
5099  nullptr,
5100  src_vector_face_access_temp,
5102  task_info.loop(worker);
5103 }
5104 
5105 
5106 
5107 template <int dim, typename Number, typename VectorizedArrayType>
5108 template <typename CLASS, typename OutVector, typename InVector>
5109 inline void
5111  void (CLASS::*function_pointer)(
5113  OutVector &,
5114  const InVector &,
5115  const std::pair<unsigned int, unsigned int> &),
5116  CLASS * owning_class,
5117  OutVector & dst,
5118  const InVector & src,
5119  const bool zero_dst_vector,
5120  const DataAccessOnFaces src_vector_face_access) const
5121 {
5122  auto src_vector_face_access_temp = src_vector_face_access;
5123  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5124  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5125  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5126  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5127 
5128  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5129  InVector,
5130  OutVector,
5131  CLASS,
5132  false>
5133  worker(*this,
5134  src,
5135  dst,
5136  zero_dst_vector,
5137  *owning_class,
5138  function_pointer,
5139  nullptr,
5140  nullptr,
5141  src_vector_face_access_temp,
5143  task_info.loop(worker);
5144 }
5145 
5146 
5147 
5148 template <int dim, typename Number, typename VectorizedArrayType>
5149 template <typename OutVector, typename InVector>
5150 inline void
5152  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5153  OutVector &,
5154  const InVector &,
5155  const std::pair<unsigned int, unsigned int> &)>
5156  & cell_operation,
5157  OutVector & dst,
5158  const InVector & src,
5159  const bool zero_dst_vector,
5160  const DataAccessOnFaces src_vector_face_access) const
5161 {
5162  auto src_vector_face_access_temp = src_vector_face_access;
5163  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5164  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5165  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5166  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5167 
5168  using Wrapper =
5169  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5170  InVector,
5171  OutVector>;
5172  Wrapper wrap(cell_operation, nullptr, nullptr);
5173 
5174  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5175  InVector,
5176  OutVector,
5177  Wrapper,
5178  true>
5179  worker(*this,
5180  src,
5181  dst,
5182  zero_dst_vector,
5183  wrap,
5184  &Wrapper::cell_integrator,
5185  &Wrapper::face_integrator,
5186  &Wrapper::boundary_integrator,
5187  src_vector_face_access_temp,
5189  task_info.loop(worker);
5190 }
5191 
5192 
5193 #endif // ifndef DOXYGEN
5194 
5195 
5196 
5198 
5199 #endif
void resize(const size_type new_size)
Number local_element(const size_type local_index) const
void reinit(const size_type size, const bool omit_zeroing_entries=false)
Abstract base class for mapping classes.
Definition: mapping.h:311
unsigned int n_active_fe_indices() const
void loop_cell_centric(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > mapping_info
Definition: matrix_free.h:2047
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:2053
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:2077
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:2094
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:2117
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:2061
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:2034
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:2105
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
const AffineConstraints< Number > & get_affine_constraints(const unsigned int dof_handler_index=0) const
unsigned int n_boundary_face_batches() const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
std::list< std::pair< bool, AlignedVector< Number > > > scratch_pad_non_threadsafe
Definition: matrix_free.h:2112
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:2012
void initialize_indices(const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const AdditionalData &additional_data)
std::vector< SmartPointer< const AffineConstraints< Number > > > affine_constraints
Definition: matrix_free.h:2020
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:2026
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:2089
internal::MatrixFreeFunctions::FaceInfo< VectorizedArrayType::size()> face_info
Definition: matrix_free.h:2084
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:2040
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:2070
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
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
UpdateFlags
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_default
No update.
static const char U
VectorType::value_type * begin(VectorType &V)
bool job_supports_mpi()
Definition: mpi.cc: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:520
void loop(MFWorkerInterface &worker) const
Definition: task_info.cc:350
std::vector< unsigned int > cell_partition_data
Definition: task_info.h:476
std::vector< unsigned int > face_partition_data
Definition: task_info.h:498