Reference documentation for deal.II version Git 2276b94619 2020-02-28 10:11:55 +0100
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
matrix_free.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2011 - 2019 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 
22 #include <deal.II/base/aligned_vector.h>
23 #include <deal.II/base/exceptions.h>
24 #include <deal.II/base/quadrature.h>
25 #include <deal.II/base/template_constraints.h>
26 #include <deal.II/base/thread_local_storage.h>
27 #include <deal.II/base/vectorization.h>
28 
29 #include <deal.II/dofs/dof_handler.h>
30 
31 #include <deal.II/fe/fe.h>
32 #include <deal.II/fe/mapping.h>
33 #include <deal.II/fe/mapping_q1.h>
34 
35 #include <deal.II/grid/grid_tools.h>
36 
37 #include <deal.II/hp/dof_handler.h>
38 #include <deal.II/hp/q_collection.h>
39 
40 #include <deal.II/lac/affine_constraints.h>
41 #include <deal.II/lac/block_vector_base.h>
42 #include <deal.II/lac/la_parallel_vector.h>
43 #include <deal.II/lac/vector_operation.h>
44 
45 #include <deal.II/matrix_free/dof_info.h>
46 #include <deal.II/matrix_free/mapping_info.h>
47 #include <deal.II/matrix_free/shape_info.h>
48 #include <deal.II/matrix_free/task_info.h>
49 #include <deal.II/matrix_free/type_traits.h>
50 
51 #include <cstdlib>
52 #include <limits>
53 #include <list>
54 #include <memory>
55 
56 
57 DEAL_II_NAMESPACE_OPEN
58 
59 
60 
114 template <int dim,
115  typename Number = double,
116  typename VectorizedArrayType = VectorizedArray<Number>>
117 class MatrixFree : public Subscriptor
118 {
119  static_assert(
120  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
121  "Type of Number and of VectorizedArrayType do not match.");
122 
123 public:
128  using value_type = Number;
129  using vectorized_value_type = VectorizedArrayType;
130 
134  static const unsigned int dimension = dim;
135 
183  {
190  {
194  none = internal::MatrixFreeFunctions::TaskInfo::none,
199  internal::MatrixFreeFunctions::TaskInfo::partition_partition,
204  internal::MatrixFreeFunctions::TaskInfo::partition_color,
209  color = internal::MatrixFreeFunctions::TaskInfo::color
210  };
211 
217  const unsigned int tasks_block_size = 0,
223  const unsigned int mg_level = numbers::invalid_unsigned_int,
224  const bool store_plain_indices = true,
225  const bool initialize_indices = true,
226  const bool initialize_mapping = true,
227  const bool overlap_communication_computation = true,
228  const bool hold_all_faces_to_owned_cells = false,
229  const bool cell_vectorization_categories_strict = false)
236  , mg_level(mg_level)
237  , level_mg_handler(this->mg_level)
245  {}
246 
259  , mg_level(other.mg_level)
260  , level_mg_handler(this->mg_level)
270  {}
271 
276  operator=(const AdditionalData &other)
277  {
286  mg_level = other.mg_level;
296 
297  return *this;
298  }
299 
337 
347  unsigned int tasks_block_size;
348 
361 
382 
403 
431 
440  unsigned int mg_level;
441 
447  DEAL_II_DEPRECATED unsigned int &level_mg_handler;
448 
456 
467 
476 
490 
499 
516  std::vector<unsigned int> cell_vectorization_category;
517 
528  };
529 
537  MatrixFree();
538 
543 
547  ~MatrixFree() override = default;
548 
561  template <typename DoFHandlerType, typename QuadratureType, typename number2>
562  void
563  reinit(const Mapping<dim> & mapping,
564  const DoFHandlerType & dof_handler,
565  const AffineConstraints<number2> &constraint,
566  const QuadratureType & quad,
567  const AdditionalData & additional_data = AdditionalData());
568 
573  template <typename DoFHandlerType, typename QuadratureType, typename number2>
574  void
575  reinit(const DoFHandlerType & dof_handler,
576  const AffineConstraints<number2> &constraint,
577  const QuadratureType & quad,
578  const AdditionalData & additional_data = AdditionalData());
579 
587  template <typename DoFHandlerType, typename QuadratureType, typename number2>
588  DEAL_II_DEPRECATED void
589  reinit(const Mapping<dim> & mapping,
590  const DoFHandlerType & dof_handler,
591  const AffineConstraints<number2> &constraint,
592  const IndexSet & locally_owned_dofs,
593  const QuadratureType & quad,
594  const AdditionalData & additional_data = AdditionalData());
595 
617  template <typename DoFHandlerType, typename QuadratureType, typename number2>
618  void
619  reinit(const Mapping<dim> & mapping,
620  const std::vector<const DoFHandlerType *> & dof_handler,
621  const std::vector<const AffineConstraints<number2> *> &constraint,
622  const std::vector<QuadratureType> & quad,
623  const AdditionalData &additional_data = AdditionalData());
624 
629  template <typename DoFHandlerType, typename QuadratureType, typename number2>
630  void
631  reinit(const std::vector<const DoFHandlerType *> & dof_handler,
632  const std::vector<const AffineConstraints<number2> *> &constraint,
633  const std::vector<QuadratureType> & quad,
634  const AdditionalData &additional_data = AdditionalData());
635 
643  template <typename DoFHandlerType, typename QuadratureType, typename number2>
644  DEAL_II_DEPRECATED void
645  reinit(const Mapping<dim> & mapping,
646  const std::vector<const DoFHandlerType *> & dof_handler,
647  const std::vector<const AffineConstraints<number2> *> &constraint,
648  const std::vector<IndexSet> & locally_owned_set,
649  const std::vector<QuadratureType> &quad,
650  const AdditionalData & additional_data = AdditionalData());
651 
659  template <typename DoFHandlerType, typename QuadratureType, typename number2>
660  void
661  reinit(const Mapping<dim> & mapping,
662  const std::vector<const DoFHandlerType *> & dof_handler,
663  const std::vector<const AffineConstraints<number2> *> &constraint,
664  const QuadratureType & quad,
665  const AdditionalData &additional_data = AdditionalData());
666 
671  template <typename DoFHandlerType, typename QuadratureType, typename number2>
672  void
673  reinit(const std::vector<const DoFHandlerType *> & dof_handler,
674  const std::vector<const AffineConstraints<number2> *> &constraint,
675  const QuadratureType & quad,
676  const AdditionalData &additional_data = AdditionalData());
677 
683  void
684  copy_from(
685  const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free_base);
686 
691  void
692  clear();
693 
695 
707  enum class DataAccessOnFaces
708  {
715  none,
716 
726  values,
727 
736  values_all_faces,
737 
748  gradients,
749 
758  gradients_all_faces,
759 
766  unspecified
767  };
768 
813  template <typename OutVector, typename InVector>
814  void
815  cell_loop(const std::function<void(
817  OutVector &,
818  const InVector &,
819  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
820  OutVector & dst,
821  const InVector & src,
822  const bool zero_dst_vector = false) const;
823 
870  template <typename CLASS, typename OutVector, typename InVector>
871  void
872  cell_loop(void (CLASS::*cell_operation)(
873  const MatrixFree &,
874  OutVector &,
875  const InVector &,
876  const std::pair<unsigned int, unsigned int> &) const,
877  const CLASS * owning_class,
878  OutVector & dst,
879  const InVector &src,
880  const bool zero_dst_vector = false) const;
881 
885  template <typename CLASS, typename OutVector, typename InVector>
886  void
887  cell_loop(void (CLASS::*cell_operation)(
888  const MatrixFree &,
889  OutVector &,
890  const InVector &,
891  const std::pair<unsigned int, unsigned int> &),
892  CLASS * owning_class,
893  OutVector & dst,
894  const InVector &src,
895  const bool zero_dst_vector = false) const;
896 
981  template <typename CLASS, typename OutVector, typename InVector>
982  void
983  cell_loop(void (CLASS::*cell_operation)(
984  const MatrixFree &,
985  OutVector &,
986  const InVector &,
987  const std::pair<unsigned int, unsigned int> &) const,
988  const CLASS * owning_class,
989  OutVector & dst,
990  const InVector &src,
991  const std::function<void(const unsigned int, const unsigned int)>
992  &operation_before_loop,
993  const std::function<void(const unsigned int, const unsigned int)>
994  & operation_after_loop,
995  const unsigned int dof_handler_index_pre_post = 0) const;
996 
1000  template <typename CLASS, typename OutVector, typename InVector>
1001  void
1002  cell_loop(void (CLASS::*cell_operation)(
1003  const MatrixFree &,
1004  OutVector &,
1005  const InVector &,
1006  const std::pair<unsigned int, unsigned int> &),
1007  CLASS * owning_class,
1008  OutVector & dst,
1009  const InVector &src,
1010  const std::function<void(const unsigned int, const unsigned int)>
1011  &operation_before_loop,
1012  const std::function<void(const unsigned int, const unsigned int)>
1013  & operation_after_loop,
1014  const unsigned int dof_handler_index_pre_post = 0) const;
1015 
1020  template <typename OutVector, typename InVector>
1021  void
1022  cell_loop(const std::function<void(
1024  OutVector &,
1025  const InVector &,
1026  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1027  OutVector & dst,
1028  const InVector & src,
1029  const std::function<void(const unsigned int, const unsigned int)>
1030  &operation_before_loop,
1031  const std::function<void(const unsigned int, const unsigned int)>
1032  & operation_after_loop,
1033  const unsigned int dof_handler_index_pre_post = 0) const;
1034 
1110  template <typename OutVector, typename InVector>
1111  void
1112  loop(const std::function<
1114  OutVector &,
1115  const InVector &,
1116  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1117  const std::function<
1119  OutVector &,
1120  const InVector &,
1121  const std::pair<unsigned int, unsigned int> &)> &face_operation,
1122  const std::function<void(
1124  OutVector &,
1125  const InVector &,
1126  const std::pair<unsigned int, unsigned int> &)> &boundary_operation,
1127  OutVector & dst,
1128  const InVector & src,
1129  const bool zero_dst_vector = false,
1130  const DataAccessOnFaces dst_vector_face_access =
1132  const DataAccessOnFaces src_vector_face_access =
1134 
1220  template <typename CLASS, typename OutVector, typename InVector>
1221  void
1222  loop(
1223  void (CLASS::*cell_operation)(const MatrixFree &,
1224  OutVector &,
1225  const InVector &,
1226  const std::pair<unsigned int, unsigned int> &)
1227  const,
1228  void (CLASS::*face_operation)(const MatrixFree &,
1229  OutVector &,
1230  const InVector &,
1231  const std::pair<unsigned int, unsigned int> &)
1232  const,
1233  void (CLASS::*boundary_operation)(
1234  const MatrixFree &,
1235  OutVector &,
1236  const InVector &,
1237  const std::pair<unsigned int, unsigned int> &) const,
1238  const CLASS * owning_class,
1239  OutVector & dst,
1240  const InVector & src,
1241  const bool zero_dst_vector = false,
1242  const DataAccessOnFaces dst_vector_face_access =
1244  const DataAccessOnFaces src_vector_face_access =
1246 
1250  template <typename CLASS, typename OutVector, typename InVector>
1251  void
1252  loop(void (CLASS::*cell_operation)(
1253  const MatrixFree &,
1254  OutVector &,
1255  const InVector &,
1256  const std::pair<unsigned int, unsigned int> &),
1257  void (CLASS::*face_operation)(
1258  const MatrixFree &,
1259  OutVector &,
1260  const InVector &,
1261  const std::pair<unsigned int, unsigned int> &),
1262  void (CLASS::*boundary_operation)(
1263  const MatrixFree &,
1264  OutVector &,
1265  const InVector &,
1266  const std::pair<unsigned int, unsigned int> &),
1267  CLASS * owning_class,
1268  OutVector & dst,
1269  const InVector & src,
1270  const bool zero_dst_vector = false,
1271  const DataAccessOnFaces dst_vector_face_access =
1273  const DataAccessOnFaces src_vector_face_access =
1275 
1340  template <typename CLASS, typename OutVector, typename InVector>
1341  void
1342  loop_cell_centric(void (CLASS::*cell_operation)(
1343  const MatrixFree &,
1344  OutVector &,
1345  const InVector &,
1346  const std::pair<unsigned int, unsigned int> &) const,
1347  const CLASS * owning_class,
1348  OutVector & dst,
1349  const InVector & src,
1350  const bool zero_dst_vector = false,
1351  const DataAccessOnFaces src_vector_face_access =
1353 
1357  template <typename CLASS, typename OutVector, typename InVector>
1358  void
1359  loop_cell_centric(void (CLASS::*cell_operation)(
1360  const MatrixFree &,
1361  OutVector &,
1362  const InVector &,
1363  const std::pair<unsigned int, unsigned int> &),
1364  CLASS * owning_class,
1365  OutVector & dst,
1366  const InVector & src,
1367  const bool zero_dst_vector = false,
1368  const DataAccessOnFaces src_vector_face_access =
1370 
1378  std::pair<unsigned int, unsigned int>
1379  create_cell_subrange_hp(const std::pair<unsigned int, unsigned int> &range,
1380  const unsigned int fe_degree,
1381  const unsigned int dof_handler_index = 0) const;
1382 
1389  std::pair<unsigned int, unsigned int>
1391  const std::pair<unsigned int, unsigned int> &range,
1392  const unsigned int fe_index,
1393  const unsigned int dof_handler_index = 0) const;
1394 
1396 
1421  template <typename VectorType>
1422  void
1423  initialize_dof_vector(VectorType & vec,
1424  const unsigned int dof_handler_index = 0) const;
1425 
1446  template <typename Number2>
1447  void
1449  const unsigned int dof_handler_index = 0) const;
1450 
1461  const std::shared_ptr<const Utilities::MPI::Partitioner> &
1462  get_vector_partitioner(const unsigned int dof_handler_index = 0) const;
1463 
1467  const IndexSet &
1468  get_locally_owned_set(const unsigned int dof_handler_index = 0) const;
1469 
1473  const IndexSet &
1474  get_ghost_set(const unsigned int dof_handler_index = 0) const;
1475 
1485  const std::vector<unsigned int> &
1486  get_constrained_dofs(const unsigned int dof_handler_index = 0) const;
1487 
1498  void
1499  renumber_dofs(std::vector<types::global_dof_index> &renumbering,
1500  const unsigned int dof_handler_index = 0);
1501 
1503 
1511  template <int spacedim>
1512  static bool
1514 
1518  unsigned int
1519  n_components() const;
1520 
1525  unsigned int
1526  n_base_elements(const unsigned int dof_handler_index) const;
1527 
1535  unsigned int
1536  n_physical_cells() const;
1537 
1547  unsigned int
1548  n_macro_cells() const;
1549 
1559  unsigned int
1560  n_cell_batches() const;
1561 
1568  unsigned int
1569  n_ghost_cell_batches() const;
1570 
1578  unsigned int
1579  n_inner_face_batches() const;
1580 
1589  unsigned int
1590  n_boundary_face_batches() const;
1591 
1596  unsigned int
1598 
1606  get_boundary_id(const unsigned int macro_face) const;
1607 
1612  std::array<types::boundary_id, VectorizedArrayType::n_array_elements>
1613  get_faces_by_cells_boundary_id(const unsigned int macro_cell,
1614  const unsigned int face_number) const;
1615 
1624  template <typename DoFHandlerType = DoFHandler<dim>>
1625  const DoFHandlerType &
1626  get_dof_handler(const unsigned int dof_handler_index = 0) const;
1627 
1639  get_cell_iterator(const unsigned int macro_cell_number,
1640  const unsigned int vector_number,
1641  const unsigned int dof_handler_index = 0) const;
1642 
1648  std::pair<int, int>
1649  get_cell_level_and_index(const unsigned int macro_cell_number,
1650  const unsigned int vector_number) const;
1651 
1664  std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>
1665  get_face_iterator(const unsigned int face_batch_number,
1666  const unsigned int vector_number,
1667  const bool interior = true,
1668  const unsigned int fe_component = 0) const;
1669 
1682  get_hp_cell_iterator(const unsigned int macro_cell_number,
1683  const unsigned int vector_number,
1684  const unsigned int dof_handler_index = 0) const;
1685 
1697  bool
1698  at_irregular_cell(const unsigned int macro_cell_number) const;
1699 
1708  unsigned int
1709  n_components_filled(const unsigned int cell_batch_number) const;
1710 
1719  unsigned int
1720  n_active_entries_per_cell_batch(const unsigned int cell_batch_number) const;
1721 
1731  unsigned int
1732  n_active_entries_per_face_batch(const unsigned int face_batch_number) const;
1733 
1737  unsigned int
1738  get_dofs_per_cell(const unsigned int dof_handler_index = 0,
1739  const unsigned int hp_active_fe_index = 0) const;
1740 
1744  unsigned int
1745  get_n_q_points(const unsigned int quad_index = 0,
1746  const unsigned int hp_active_fe_index = 0) const;
1747 
1752  unsigned int
1753  get_dofs_per_face(const unsigned int dof_handler_index = 0,
1754  const unsigned int hp_active_fe_index = 0) const;
1755 
1760  unsigned int
1761  get_n_q_points_face(const unsigned int quad_index = 0,
1762  const unsigned int hp_active_fe_index = 0) const;
1763 
1767  const Quadrature<dim> &
1768  get_quadrature(const unsigned int quad_index = 0,
1769  const unsigned int hp_active_fe_index = 0) const;
1770 
1774  const Quadrature<dim - 1> &
1775  get_face_quadrature(const unsigned int quad_index = 0,
1776  const unsigned int hp_active_fe_index = 0) const;
1777 
1784  unsigned int
1785  get_cell_category(const unsigned int macro_cell) const;
1786 
1791  std::pair<unsigned int, unsigned int>
1792  get_face_category(const unsigned int macro_face) const;
1793 
1797  bool
1798  indices_initialized() const;
1799 
1804  bool
1805  mapping_initialized() const;
1806 
1811  unsigned int
1812  get_mg_level() const;
1813 
1818  std::size_t
1819  memory_consumption() const;
1820 
1825  template <typename StreamType>
1826  void
1827  print_memory_consumption(StreamType &out) const;
1828 
1833  void
1834  print(std::ostream &out) const;
1835 
1837 
1848  get_task_info() const;
1849 
1850  /*
1851  * Return geometry-dependent information on the cells.
1852  */
1853  const internal::MatrixFreeFunctions::
1854  MappingInfo<dim, Number, VectorizedArrayType> &
1855  get_mapping_info() const;
1856 
1861  get_dof_info(const unsigned int dof_handler_index_component = 0) const;
1862 
1866  unsigned int
1867  n_constraint_pool_entries() const;
1868 
1873  const Number *
1874  constraint_pool_begin(const unsigned int pool_index) const;
1875 
1881  const Number *
1882  constraint_pool_end(const unsigned int pool_index) const;
1883 
1888  get_shape_info(const unsigned int dof_handler_index_component = 0,
1889  const unsigned int quad_index = 0,
1890  const unsigned int fe_base_element = 0,
1891  const unsigned int hp_active_fe_index = 0,
1892  const unsigned int hp_active_quad_index = 0) const;
1893 
1898  VectorizedArrayType::n_array_elements> &
1899  get_face_info(const unsigned int face_batch_number) const;
1900 
1901 
1907  const Table<3, unsigned int> &
1909 
1924  acquire_scratch_data() const;
1925 
1929  void
1931 
1943 
1947  void
1949  const AlignedVector<Number> *memory) const;
1950 
1952 
1953 private:
1958  template <typename number2>
1959  void
1961  const Mapping<dim> & mapping,
1962  const std::vector<const DoFHandler<dim> *> & dof_handler,
1963  const std::vector<const AffineConstraints<number2> *> &constraint,
1964  const std::vector<IndexSet> & locally_owned_set,
1965  const std::vector<hp::QCollection<1>> & quad,
1966  const AdditionalData & additional_data);
1967 
1971  template <typename number2>
1972  void
1974  const Mapping<dim> & mapping,
1975  const std::vector<const hp::DoFHandler<dim> *> & dof_handler,
1976  const std::vector<const AffineConstraints<number2> *> &constraint,
1977  const std::vector<IndexSet> & locally_owned_set,
1978  const std::vector<hp::QCollection<1>> & quad,
1979  const AdditionalData & additional_data);
1980 
1987  template <typename number2>
1988  void
1990  const std::vector<const AffineConstraints<number2> *> &constraint,
1991  const std::vector<IndexSet> & locally_owned_set,
1992  const AdditionalData & additional_data);
1993 
1997  void
1999  const std::vector<const DoFHandler<dim> *> &dof_handlers,
2000  const AdditionalData & additional_data);
2001 
2005  void
2007  const std::vector<const hp::DoFHandler<dim> *> &dof_handlers,
2008  const AdditionalData & additional_data);
2009 
2014  void
2016 
2023  {
2024  DoFHandlers()
2025  : active_dof_handler(usual)
2026  , n_dof_handlers(0)
2028  {}
2029 
2030  std::vector<SmartPointer<const DoFHandler<dim>>> dof_handler;
2031  std::vector<SmartPointer<const hp::DoFHandler<dim>>> hp_dof_handler;
2033  {
2042  } active_dof_handler;
2043  unsigned int n_dof_handlers;
2044  unsigned int level;
2045  };
2046 
2051 
2056  std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
2057 
2064  std::vector<Number> constraint_pool_data;
2065 
2070  std::vector<unsigned int> constraint_pool_row_index;
2071 
2078 
2084 
2091  std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
2092 
2093 
2101 
2108 
2115 
2120 
2125 
2134  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2136 
2141  mutable std::list<std::pair<bool, AlignedVector<Number>>>
2143 
2147  unsigned int mg_level;
2148 };
2149 
2150 
2151 
2152 /*----------------------- Inline functions ----------------------------------*/
2153 
2154 #ifndef DOXYGEN
2155 
2156 
2157 
2158 template <int dim, typename Number, typename VectorizedArrayType>
2159 template <typename VectorType>
2160 inline void
2162  VectorType & vec,
2163  const unsigned int comp) const
2164 {
2165  AssertIndexRange(comp, n_components());
2166  vec.reinit(dof_info[comp].vector_partitioner->size());
2167 }
2168 
2169 
2170 
2171 template <int dim, typename Number, typename VectorizedArrayType>
2172 template <typename Number2>
2173 inline void
2176  const unsigned int comp) const
2177 {
2178  AssertIndexRange(comp, n_components());
2179  vec.reinit(dof_info[comp].vector_partitioner);
2180 }
2181 
2182 
2183 
2184 template <int dim, typename Number, typename VectorizedArrayType>
2185 inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2187  const unsigned int comp) const
2188 {
2189  AssertIndexRange(comp, n_components());
2190  return dof_info[comp].vector_partitioner;
2191 }
2192 
2193 
2194 
2195 template <int dim, typename Number, typename VectorizedArrayType>
2196 inline const std::vector<unsigned int> &
2198  const unsigned int comp) const
2199 {
2200  AssertIndexRange(comp, n_components());
2201  return dof_info[comp].constrained_dofs;
2202 }
2203 
2204 
2205 
2206 template <int dim, typename Number, typename VectorizedArrayType>
2207 inline unsigned int
2209 {
2210  AssertDimension(dof_handlers.n_dof_handlers, dof_info.size());
2211  return dof_handlers.n_dof_handlers;
2212 }
2213 
2214 
2215 
2216 template <int dim, typename Number, typename VectorizedArrayType>
2217 inline unsigned int
2219  const unsigned int dof_no) const
2220 {
2221  AssertDimension(dof_handlers.n_dof_handlers, dof_info.size());
2222  AssertIndexRange(dof_no, dof_handlers.n_dof_handlers);
2223  return dof_handlers.dof_handler[dof_no]->get_fe().n_base_elements();
2224 }
2225 
2226 
2227 
2228 template <int dim, typename Number, typename VectorizedArrayType>
2231 {
2232  return task_info;
2233 }
2234 
2235 
2236 
2237 template <int dim, typename Number, typename VectorizedArrayType>
2238 inline unsigned int
2240 {
2241  return *(task_info.cell_partition_data.end() - 2);
2242 }
2243 
2244 
2245 
2246 template <int dim, typename Number, typename VectorizedArrayType>
2247 inline unsigned int
2249 {
2250  return task_info.n_active_cells;
2251 }
2252 
2253 
2254 
2255 template <int dim, typename Number, typename VectorizedArrayType>
2256 inline unsigned int
2258 {
2259  return *(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  return *(task_info.cell_partition_data.end() - 1) -
2269  *(task_info.cell_partition_data.end() - 2);
2270 }
2271 
2272 
2273 
2274 template <int dim, typename Number, typename VectorizedArrayType>
2275 inline unsigned int
2277 {
2278  if (task_info.face_partition_data.size() == 0)
2279  return 0;
2280  return task_info.face_partition_data.back();
2281 }
2282 
2283 
2284 
2285 template <int dim, typename Number, typename VectorizedArrayType>
2286 inline unsigned int
2288 {
2289  if (task_info.face_partition_data.size() == 0)
2290  return 0;
2291  return task_info.boundary_partition_data.back() -
2293 }
2294 
2295 
2296 
2297 template <int dim, typename Number, typename VectorizedArrayType>
2298 inline unsigned int
2300 {
2301  if (task_info.face_partition_data.size() == 0)
2302  return 0;
2303  return face_info.faces.size() - task_info.boundary_partition_data.back();
2304 }
2305 
2306 
2307 
2308 template <int dim, typename Number, typename VectorizedArrayType>
2309 inline types::boundary_id
2311  const unsigned int macro_face) const
2312 {
2313  Assert(macro_face >= task_info.boundary_partition_data[0] &&
2314  macro_face < task_info.boundary_partition_data.back(),
2315  ExcIndexRange(macro_face,
2318  return types::boundary_id(face_info.faces[macro_face].exterior_face_no);
2319 }
2320 
2321 
2322 
2323 template <int dim, typename Number, typename VectorizedArrayType>
2324 inline std::array<types::boundary_id, VectorizedArrayType::n_array_elements>
2326  const unsigned int macro_cell,
2327  const unsigned int face_number) const
2328 {
2329  AssertIndexRange(macro_cell, n_macro_cells());
2332  ExcNotInitialized());
2333  std::array<types::boundary_id, VectorizedArrayType::n_array_elements> result;
2334  result.fill(numbers::invalid_boundary_id);
2335  for (unsigned int v = 0; v < n_active_entries_per_cell_batch(macro_cell); ++v)
2336  result[v] = face_info.cell_and_face_boundary_id(macro_cell, face_number, v);
2337  return result;
2338 }
2339 
2340 
2341 
2342 template <int dim, typename Number, typename VectorizedArrayType>
2343 inline const internal::MatrixFreeFunctions::
2344  MappingInfo<dim, Number, VectorizedArrayType> &
2346 {
2347  return mapping_info;
2348 }
2349 
2350 
2351 
2352 template <int dim, typename Number, typename VectorizedArrayType>
2355  const unsigned int dof_index) const
2356 {
2357  AssertIndexRange(dof_index, n_components());
2358  return dof_info[dof_index];
2359 }
2360 
2361 
2362 
2363 template <int dim, typename Number, typename VectorizedArrayType>
2364 inline unsigned int
2366 {
2367  return constraint_pool_row_index.size() - 1;
2368 }
2369 
2370 
2371 
2372 template <int dim, typename Number, typename VectorizedArrayType>
2373 inline const Number *
2375  const unsigned int row) const
2376 {
2378  return constraint_pool_data.empty() ?
2379  nullptr :
2381 }
2382 
2383 
2384 
2385 template <int dim, typename Number, typename VectorizedArrayType>
2386 inline const Number *
2388  const unsigned int row) const
2389 {
2391  return constraint_pool_data.empty() ?
2392  nullptr :
2394 }
2395 
2396 
2397 
2398 template <int dim, typename Number, typename VectorizedArrayType>
2399 inline std::pair<unsigned int, unsigned int>
2401  const std::pair<unsigned int, unsigned int> &range,
2402  const unsigned int degree,
2403  const unsigned int dof_handler_component) const
2404 {
2405  if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2406  {
2408  dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2410  dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2411  if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2412  return range;
2413  else
2414  return {range.second, range.second};
2415  }
2416 
2417  const unsigned int fe_index =
2418  dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2419  if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2420  return {range.second, range.second};
2421  else
2422  return create_cell_subrange_hp_by_index(range,
2423  fe_index,
2424  dof_handler_component);
2425 }
2426 
2427 
2428 
2429 template <int dim, typename Number, typename VectorizedArrayType>
2430 inline bool
2432  const unsigned int macro_cell) const
2433 {
2434  AssertIndexRange(macro_cell, task_info.cell_partition_data.back());
2435  return VectorizedArrayType::n_array_elements > 1 &&
2436  cell_level_index[(macro_cell + 1) *
2437  VectorizedArrayType::n_array_elements -
2438  1] ==
2439  cell_level_index[(macro_cell + 1) *
2440  VectorizedArrayType::n_array_elements -
2441  2];
2442 }
2443 
2444 
2445 
2446 template <int dim, typename Number, typename VectorizedArrayType>
2447 inline unsigned int
2449  const unsigned int cell_batch_number) const
2450 {
2451  return n_active_entries_per_cell_batch(cell_batch_number);
2452 }
2453 
2454 
2455 
2456 template <int dim, typename Number, typename VectorizedArrayType>
2457 inline unsigned int
2459  const unsigned int cell_batch_number) const
2460 {
2461  AssertIndexRange(cell_batch_number, task_info.cell_partition_data.back());
2462  unsigned int n_components = VectorizedArrayType::n_array_elements;
2463  while (
2464  n_components > 1 &&
2465  cell_level_index[cell_batch_number * VectorizedArrayType::n_array_elements +
2466  n_components - 1] ==
2467  cell_level_index[cell_batch_number *
2468  VectorizedArrayType::n_array_elements +
2469  n_components - 2])
2470  --n_components;
2471  AssertIndexRange(n_components - 1, VectorizedArrayType::n_array_elements);
2472  return n_components;
2473 }
2474 
2475 
2476 
2477 template <int dim, typename Number, typename VectorizedArrayType>
2478 inline unsigned int
2480  const unsigned int face_batch_number) const
2481 {
2482  AssertIndexRange(face_batch_number, face_info.faces.size());
2483  unsigned int n_components = VectorizedArrayType::n_array_elements;
2484  while (n_components > 1 &&
2485  face_info.faces[face_batch_number].cells_interior[n_components - 1] ==
2487  --n_components;
2488  AssertIndexRange(n_components - 1, VectorizedArrayType::n_array_elements);
2489  return n_components;
2490 }
2491 
2492 
2493 
2494 template <int dim, typename Number, typename VectorizedArrayType>
2495 inline unsigned int
2497  const unsigned int dof_handler_index,
2498  const unsigned int active_fe_index) const
2499 {
2500  return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2501 }
2502 
2503 
2504 
2505 template <int dim, typename Number, typename VectorizedArrayType>
2506 inline unsigned int
2508  const unsigned int quad_index,
2509  const unsigned int active_fe_index) const
2510 {
2511  AssertIndexRange(quad_index, mapping_info.cell_data.size());
2512  return mapping_info.cell_data[quad_index]
2513  .descriptor[active_fe_index]
2514  .n_q_points;
2515 }
2516 
2517 
2518 
2519 template <int dim, typename Number, typename VectorizedArrayType>
2520 inline unsigned int
2522  const unsigned int dof_handler_index,
2523  const unsigned int active_fe_index) const
2524 {
2525  return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2526 }
2527 
2528 
2529 
2530 template <int dim, typename Number, typename VectorizedArrayType>
2531 inline unsigned int
2533  const unsigned int quad_index,
2534  const unsigned int active_fe_index) const
2535 {
2536  AssertIndexRange(quad_index, mapping_info.face_data.size());
2537  return mapping_info.face_data[quad_index]
2538  .descriptor[active_fe_index]
2539  .n_q_points;
2540 }
2541 
2542 
2543 
2544 template <int dim, typename Number, typename VectorizedArrayType>
2545 inline const IndexSet &
2547  const unsigned int dof_handler_index) const
2548 {
2549  return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2550 }
2551 
2552 
2553 
2554 template <int dim, typename Number, typename VectorizedArrayType>
2555 inline const IndexSet &
2557  const unsigned int dof_handler_index) const
2558 {
2559  return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2560 }
2561 
2562 
2563 
2564 template <int dim, typename Number, typename VectorizedArrayType>
2567  const unsigned int dof_handler_index,
2568  const unsigned int index_quad,
2569  const unsigned int index_fe,
2570  const unsigned int active_fe_index,
2571  const unsigned int active_quad_index) const
2572 {
2573  AssertIndexRange(dof_handler_index, dof_info.size());
2574  const unsigned int ind =
2575  dof_info[dof_handler_index].global_base_element_offset + index_fe;
2576  AssertIndexRange(ind, shape_info.size(0));
2577  AssertIndexRange(index_quad, shape_info.size(1));
2578  AssertIndexRange(active_fe_index, shape_info.size(2));
2579  AssertIndexRange(active_quad_index, shape_info.size(3));
2580  return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2581 }
2582 
2583 
2584 
2585 template <int dim, typename Number, typename VectorizedArrayType>
2587  VectorizedArrayType::n_array_elements> &
2589  const unsigned int macro_face) const
2590 {
2591  AssertIndexRange(macro_face, face_info.faces.size());
2592  return face_info.faces[macro_face];
2593 }
2594 
2595 
2596 
2597 template <int dim, typename Number, typename VectorizedArrayType>
2598 inline const Table<3, unsigned int> &
2600  const
2601 {
2603 }
2604 
2605 
2606 
2607 template <int dim, typename Number, typename VectorizedArrayType>
2608 inline const Quadrature<dim> &
2610  const unsigned int quad_index,
2611  const unsigned int active_fe_index) const
2612 {
2613  AssertIndexRange(quad_index, mapping_info.cell_data.size());
2614  return mapping_info.cell_data[quad_index]
2615  .descriptor[active_fe_index]
2616  .quadrature;
2617 }
2618 
2619 
2620 
2621 template <int dim, typename Number, typename VectorizedArrayType>
2622 inline const Quadrature<dim - 1> &
2624  const unsigned int quad_index,
2625  const unsigned int active_fe_index) const
2626 {
2627  AssertIndexRange(quad_index, mapping_info.face_data.size());
2628  return mapping_info.face_data[quad_index]
2629  .descriptor[active_fe_index]
2630  .quadrature;
2631 }
2632 
2633 
2634 
2635 template <int dim, typename Number, typename VectorizedArrayType>
2636 inline unsigned int
2638  const unsigned int macro_cell) const
2639 {
2640  AssertIndexRange(0, dof_info.size());
2641  AssertIndexRange(macro_cell, dof_info[0].cell_active_fe_index.size());
2642  if (dof_info[0].cell_active_fe_index.empty())
2643  return 0;
2644  else
2645  return dof_info[0].cell_active_fe_index[macro_cell];
2646 }
2647 
2648 
2649 
2650 template <int dim, typename Number, typename VectorizedArrayType>
2651 inline std::pair<unsigned int, unsigned int>
2653  const unsigned int macro_face) const
2654 {
2655  AssertIndexRange(macro_face, face_info.faces.size());
2656  if (dof_info[0].cell_active_fe_index.empty())
2657  return std::make_pair(0U, 0U);
2658 
2659  std::pair<unsigned int, unsigned int> result;
2660  for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements &&
2661  face_info.faces[macro_face].cells_interior[v] !=
2663  ++v)
2664  result.first = std::max(
2665  result.first,
2666  dof_info[0]
2667  .cell_active_fe_index[face_info.faces[macro_face].cells_interior[v]]);
2668  if (face_info.faces[macro_face].cells_exterior[0] !=
2670  for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements &&
2671  face_info.faces[macro_face].cells_exterior[v] !=
2673  ++v)
2674  result.second = std::max(
2675  result.first,
2676  dof_info[0]
2677  .cell_active_fe_index[face_info.faces[macro_face].cells_exterior[v]]);
2678  else
2679  result.second = numbers::invalid_unsigned_int;
2680  return result;
2681 }
2682 
2683 
2684 
2685 template <int dim, typename Number, typename VectorizedArrayType>
2686 inline bool
2688 {
2689  return indices_are_initialized;
2690 }
2691 
2692 
2693 
2694 template <int dim, typename Number, typename VectorizedArrayType>
2695 inline bool
2697 {
2698  return mapping_is_initialized;
2699 }
2700 
2701 
2702 template <int dim, typename Number, typename VectorizedArrayType>
2703 inline unsigned int
2705 {
2706  return mg_level;
2707 }
2708 
2709 
2710 
2711 template <int dim, typename Number, typename VectorizedArrayType>
2714 {
2715  using list_type =
2716  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
2717  list_type &data = scratch_pad.get();
2718  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
2719  if (it->first == false)
2720  {
2721  it->first = true;
2722  return &it->second;
2723  }
2724  data.emplace_front(true, AlignedVector<VectorizedArrayType>());
2725  return &data.front().second;
2726 }
2727 
2728 
2729 
2730 template <int dim, typename Number, typename VectorizedArrayType>
2731 void
2733  const AlignedVector<VectorizedArrayType> *scratch) const
2734 {
2735  using list_type =
2736  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
2737  list_type &data = scratch_pad.get();
2738  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
2739  if (&it->second == scratch)
2740  {
2741  Assert(it->first == true, ExcInternalError());
2742  it->first = false;
2743  return;
2744  }
2745  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2746 }
2747 
2748 
2749 
2750 template <int dim, typename Number, typename VectorizedArrayType>
2754 {
2755  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
2757  it != scratch_pad_non_threadsafe.end();
2758  ++it)
2759  if (it->first == false)
2760  {
2761  it->first = true;
2762  return &it->second;
2763  }
2764  scratch_pad_non_threadsafe.push_front(
2765  std::make_pair(true, AlignedVector<Number>()));
2766  return &scratch_pad_non_threadsafe.front().second;
2767 }
2768 
2769 
2770 
2771 template <int dim, typename Number, typename VectorizedArrayType>
2772 void
2775  const AlignedVector<Number> *scratch) const
2776 {
2777  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
2779  it != scratch_pad_non_threadsafe.end();
2780  ++it)
2781  if (&it->second == scratch)
2782  {
2783  Assert(it->first == true, ExcInternalError());
2784  it->first = false;
2785  return;
2786  }
2787  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2788 }
2789 
2790 
2791 
2792 // ------------------------------ reinit functions ---------------------------
2793 
2794 namespace internal
2795 {
2796  namespace MatrixFreeImplementation
2797  {
2798  template <typename DoFHandlerType>
2799  inline std::vector<IndexSet>
2800  extract_locally_owned_index_sets(
2801  const std::vector<const DoFHandlerType *> &dofh,
2802  const unsigned int level)
2803  {
2804  std::vector<IndexSet> locally_owned_set;
2805  locally_owned_set.reserve(dofh.size());
2806  for (unsigned int j = 0; j < dofh.size(); j++)
2807  if (level == numbers::invalid_unsigned_int)
2808  locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
2809  else
2810  AssertThrow(false, ExcNotImplemented());
2811  return locally_owned_set;
2812  }
2813 
2814  template <int dim, int spacedim>
2815  inline std::vector<IndexSet>
2816  extract_locally_owned_index_sets(
2817  const std::vector<const ::DoFHandler<dim, spacedim> *> &dofh,
2818  const unsigned int level)
2819  {
2820  std::vector<IndexSet> locally_owned_set;
2821  locally_owned_set.reserve(dofh.size());
2822  for (unsigned int j = 0; j < dofh.size(); j++)
2823  if (level == numbers::invalid_unsigned_int)
2824  locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
2825  else
2826  locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
2827  return locally_owned_set;
2828  }
2829  } // namespace MatrixFreeImplementation
2830 } // namespace internal
2831 
2832 
2833 
2834 template <int dim, typename Number, typename VectorizedArrayType>
2835 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2836 void
2838  const DoFHandlerType & dof_handler,
2839  const AffineConstraints<number2> &constraints_in,
2840  const QuadratureType & quad,
2842  &additional_data)
2843 {
2844  std::vector<const DoFHandlerType *> dof_handlers;
2845  std::vector<const AffineConstraints<number2> *> constraints;
2846  std::vector<QuadratureType> quads;
2847 
2848  dof_handlers.push_back(&dof_handler);
2849  constraints.push_back(&constraints_in);
2850  quads.push_back(quad);
2851 
2852  std::vector<IndexSet> locally_owned_sets =
2853  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2854  dof_handlers, additional_data.mg_level);
2855 
2856  std::vector<hp::QCollection<1>> quad_hp;
2857  quad_hp.emplace_back(quad);
2858 
2860  dof_handlers,
2861  constraints,
2862  locally_owned_sets,
2863  quad_hp,
2864  additional_data);
2865 }
2866 
2867 
2868 
2869 template <int dim, typename Number, typename VectorizedArrayType>
2870 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2871 void
2873  const Mapping<dim> & mapping,
2874  const DoFHandlerType & dof_handler,
2875  const AffineConstraints<number2> &constraints_in,
2876  const QuadratureType & quad,
2878  &additional_data)
2879 {
2880  std::vector<const DoFHandlerType *> dof_handlers;
2881  std::vector<const AffineConstraints<number2> *> constraints;
2882 
2883  dof_handlers.push_back(&dof_handler);
2884  constraints.push_back(&constraints_in);
2885 
2886  std::vector<IndexSet> locally_owned_sets =
2887  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2888  dof_handlers, additional_data.mg_level);
2889 
2890  std::vector<hp::QCollection<1>> quad_hp;
2891  quad_hp.emplace_back(quad);
2892 
2893  internal_reinit(mapping,
2894  dof_handlers,
2895  constraints,
2896  locally_owned_sets,
2897  quad_hp,
2898  additional_data);
2899 }
2900 
2901 
2902 
2903 template <int dim, typename Number, typename VectorizedArrayType>
2904 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2905 void
2907  const std::vector<const DoFHandlerType *> & dof_handler,
2908  const std::vector<const AffineConstraints<number2> *> &constraint,
2909  const std::vector<QuadratureType> & quad,
2911  &additional_data)
2912 {
2913  std::vector<IndexSet> locally_owned_set =
2914  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2915  dof_handler, additional_data.mg_level);
2916  std::vector<hp::QCollection<1>> quad_hp;
2917  for (unsigned int q = 0; q < quad.size(); ++q)
2918  quad_hp.emplace_back(quad[q]);
2920  dof_handler,
2921  constraint,
2922  locally_owned_set,
2923  quad_hp,
2924  additional_data);
2925 }
2926 
2927 
2928 
2929 template <int dim, typename Number, typename VectorizedArrayType>
2930 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2931 void
2933  const std::vector<const DoFHandlerType *> & dof_handler,
2934  const std::vector<const AffineConstraints<number2> *> &constraint,
2935  const QuadratureType & quad,
2937  &additional_data)
2938 {
2939  std::vector<IndexSet> locally_owned_set =
2940  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2941  dof_handler, additional_data.mg_level);
2942  std::vector<hp::QCollection<1>> quad_hp;
2943  quad_hp.emplace_back(quad);
2945  dof_handler,
2946  constraint,
2947  locally_owned_set,
2948  quad_hp,
2949  additional_data);
2950 }
2951 
2952 
2953 
2954 template <int dim, typename Number, typename VectorizedArrayType>
2955 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2956 void
2958  const Mapping<dim> & mapping,
2959  const std::vector<const DoFHandlerType *> & dof_handler,
2960  const std::vector<const AffineConstraints<number2> *> &constraint,
2961  const QuadratureType & quad,
2963  &additional_data)
2964 {
2965  std::vector<IndexSet> locally_owned_set =
2966  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2967  dof_handler, additional_data.mg_level);
2968  std::vector<hp::QCollection<1>> quad_hp;
2969  quad_hp.emplace_back(quad);
2970  internal_reinit(mapping,
2971  dof_handler,
2972  constraint,
2973  locally_owned_set,
2974  quad_hp,
2975  additional_data);
2976 }
2977 
2978 
2979 
2980 template <int dim, typename Number, typename VectorizedArrayType>
2981 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2982 void
2984  const Mapping<dim> & mapping,
2985  const std::vector<const DoFHandlerType *> & dof_handler,
2986  const std::vector<const AffineConstraints<number2> *> &constraint,
2987  const std::vector<QuadratureType> & quad,
2989  &additional_data)
2990 {
2991  std::vector<IndexSet> locally_owned_set =
2992  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2993  dof_handler, additional_data.mg_level);
2994  std::vector<hp::QCollection<1>> quad_hp;
2995  for (unsigned int q = 0; q < quad.size(); ++q)
2996  quad_hp.emplace_back(quad[q]);
2997  internal_reinit(mapping,
2998  dof_handler,
2999  constraint,
3000  locally_owned_set,
3001  quad_hp,
3002  additional_data);
3003 }
3004 
3005 
3006 
3007 template <int dim, typename Number, typename VectorizedArrayType>
3008 template <typename DoFHandlerType, typename QuadratureType, typename number2>
3009 void
3011  const Mapping<dim> & mapping,
3012  const std::vector<const DoFHandlerType *> & dof_handler,
3013  const std::vector<const AffineConstraints<number2> *> &constraint,
3014  const std::vector<IndexSet> & locally_owned_set,
3015  const std::vector<QuadratureType> & quad,
3017  &additional_data)
3018 {
3019  // find out whether we use a hp Quadrature or a standard quadrature
3020  std::vector<hp::QCollection<1>> quad_hp;
3021  for (unsigned int q = 0; q < quad.size(); ++q)
3022  quad_hp.emplace_back(quad[q]);
3023  internal_reinit(mapping,
3024  dof_handler,
3025  constraint,
3026  locally_owned_set,
3027  quad_hp,
3028  additional_data);
3029 }
3030 
3031 
3032 
3033 // ------------------------------ implementation of loops --------------------
3034 
3035 // internal helper functions that define how to call MPI data exchange
3036 // functions: for generic vectors, do nothing at all. For distributed vectors,
3037 // call update_ghost_values_start function and so on. If we have collections
3038 // of vectors, just do the individual functions of the components. In order to
3039 // keep ghost values consistent (whether we are in read or write mode), we
3040 // also reset the values at the end. the whole situation is a bit complicated
3041 // by the fact that we need to treat block vectors differently, which use some
3042 // additional helper functions to select the blocks and template magic.
3043 namespace internal
3044 {
3048  template <int dim, typename Number, typename VectorizedArrayType>
3049  struct VectorDataExchange
3050  {
3051  // A shift for the MPI messages to reduce the risk for accidental
3052  // interaction with other open communications that a user program might
3053  // set up (parallel vectors support unfinished communication). We let
3054  // the other vectors use the first 20 assigned numbers and start the
3055  // matrix-free communication.
3056  static constexpr unsigned int channel_shift = 20;
3057 
3058 
3059 
3064  VectorDataExchange(
3065  const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3066  const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3067  DataAccessOnFaces vector_face_access,
3068  const unsigned int n_components)
3069  : matrix_free(matrix_free)
3070  , vector_face_access(
3071  matrix_free.get_task_info().face_partition_data.empty() ?
3073  DataAccessOnFaces::unspecified :
3074  vector_face_access)
3075  , ghosts_were_set(false)
3076 # ifdef DEAL_II_WITH_MPI
3077  , tmp_data(n_components)
3078  , requests(n_components)
3079 # endif
3080  {
3081  (void)n_components;
3082  if (this->vector_face_access !=
3085  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3087  matrix_free.get_dof_info(c).vector_partitioner_face_variants.size(),
3088  5);
3089  }
3090 
3091 
3092 
3096  ~VectorDataExchange() // NOLINT
3097  {
3098 # ifdef DEAL_II_WITH_MPI
3099  for (unsigned int i = 0; i < tmp_data.size(); ++i)
3100  if (tmp_data[i] != nullptr)
3101  matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3102 # endif
3103  }
3104 
3105 
3106 
3111  template <typename VectorType>
3112  unsigned int
3113  find_vector_in_mf(const VectorType &vec,
3114  const bool check_global_compatibility = true) const
3115  {
3116  unsigned int mf_component = numbers::invalid_unsigned_int;
3117  (void)check_global_compatibility;
3118  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3119  if (
3120 # ifdef DEBUG
3121  check_global_compatibility ?
3122  vec.get_partitioner()->is_globally_compatible(
3123  *matrix_free.get_dof_info(c).vector_partitioner) :
3124 # endif
3125  vec.get_partitioner()->is_compatible(
3126  *matrix_free.get_dof_info(c).vector_partitioner))
3127  {
3128  mf_component = c;
3129  break;
3130  }
3131  return mf_component;
3132  }
3133 
3134 
3135 
3141  get_partitioner(const unsigned int mf_component) const
3142  {
3143  AssertDimension(matrix_free.get_dof_info(mf_component)
3144  .vector_partitioner_face_variants.size(),
3145  5);
3146  if (vector_face_access ==
3149  return *matrix_free.get_dof_info(mf_component)
3150  .vector_partitioner_face_variants[0];
3151  else if (vector_face_access ==
3154  return *matrix_free.get_dof_info(mf_component)
3155  .vector_partitioner_face_variants[1];
3156  else if (vector_face_access ==
3159  return *matrix_free.get_dof_info(mf_component)
3160  .vector_partitioner_face_variants[2];
3161  else if (vector_face_access ==
3164  return *matrix_free.get_dof_info(mf_component)
3165  .vector_partitioner_face_variants[3];
3166  else /*if (vector_face_access ==
3167  ::MatrixFree<dim,
3168  Number>::DataAccessOnFaces::gradients_all_faces)*/
3169  return *matrix_free.get_dof_info(mf_component)
3170  .vector_partitioner_face_variants[4];
3171  }
3172 
3173 
3174 
3178  template <typename VectorType,
3179  typename std::enable_if<is_serial_or_dummy<VectorType>::value,
3180  VectorType>::type * = nullptr>
3181  void
3182  update_ghost_values_start(const unsigned int /*component_in_block_vector*/,
3183  const VectorType & /*vec*/)
3184  {}
3185 
3186 
3191  template <typename VectorType,
3192  typename std::enable_if<
3193  !has_update_ghost_values_start<VectorType>::value &&
3194  !is_serial_or_dummy<VectorType>::value,
3195  VectorType>::type * = nullptr>
3196  void
3197  update_ghost_values_start(const unsigned int component_in_block_vector,
3198  const VectorType & vec)
3199  {
3200  (void)component_in_block_vector;
3201  bool ghosts_set = vec.has_ghost_elements();
3202  if (ghosts_set)
3203  ghosts_were_set = true;
3204 
3205  vec.update_ghost_values();
3206  }
3207 
3208 
3209 
3215  template <typename VectorType,
3216  typename std::enable_if<
3217  has_update_ghost_values_start<VectorType>::value &&
3218  !has_exchange_on_subset<VectorType>::value,
3219  VectorType>::type * = nullptr>
3220  void
3221  update_ghost_values_start(const unsigned int component_in_block_vector,
3222  const VectorType & vec)
3223  {
3224  (void)component_in_block_vector;
3225  bool ghosts_set = vec.has_ghost_elements();
3226  if (ghosts_set)
3227  ghosts_were_set = true;
3228 
3229  vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3230  }
3231 
3232 
3233 
3240  template <typename VectorType,
3241  typename std::enable_if<
3242  has_update_ghost_values_start<VectorType>::value &&
3243  has_exchange_on_subset<VectorType>::value,
3244  VectorType>::type * = nullptr>
3245  void
3246  update_ghost_values_start(const unsigned int component_in_block_vector,
3247  const VectorType & vec)
3248  {
3249  static_assert(
3250  std::is_same<Number, typename VectorType::value_type>::value,
3251  "Type mismatch between VectorType and VectorDataExchange");
3252  (void)component_in_block_vector;
3253  bool ghosts_set = vec.has_ghost_elements();
3254  if (ghosts_set)
3255  ghosts_were_set = true;
3256  if (vector_face_access ==
3259  vec.size() == 0)
3260  vec.update_ghost_values_start(component_in_block_vector +
3261  channel_shift);
3262  else
3263  {
3264 # ifdef DEAL_II_WITH_MPI
3265  const unsigned int mf_component = find_vector_in_mf(vec);
3266  if (&get_partitioner(mf_component) ==
3267  matrix_free.get_dof_info(mf_component).vector_partitioner.get())
3268  {
3269  vec.update_ghost_values_start(component_in_block_vector +
3270  channel_shift);
3271  return;
3272  }
3273 
3274  const Utilities::MPI::Partitioner &part =
3275  get_partitioner(mf_component);
3276  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0)
3277  return;
3278 
3279  tmp_data[component_in_block_vector] =
3280  matrix_free.acquire_scratch_data_non_threadsafe();
3281  tmp_data[component_in_block_vector]->resize_fast(
3282  part.n_import_indices());
3283  AssertDimension(requests.size(), tmp_data.size());
3284 
3286  component_in_block_vector + channel_shift,
3287  ArrayView<const Number>(vec.begin(), part.local_size()),
3288  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3289  part.n_import_indices()),
3290  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3291  vec.get_partitioner()->local_size(),
3292  vec.get_partitioner()->n_ghost_indices()),
3293  this->requests[component_in_block_vector]);
3294 # endif
3295  }
3296  }
3297 
3298 
3299 
3304  template <
3305  typename VectorType,
3306  typename std::enable_if<!has_update_ghost_values_start<VectorType>::value,
3307  VectorType>::type * = nullptr>
3308  void
3309  update_ghost_values_finish(const unsigned int /*component_in_block_vector*/,
3310  const VectorType & /*vec*/)
3311  {}
3312 
3313 
3314 
3320  template <typename VectorType,
3321  typename std::enable_if<
3322  has_update_ghost_values_start<VectorType>::value &&
3323  !has_exchange_on_subset<VectorType>::value,
3324  VectorType>::type * = nullptr>
3325  void
3326  update_ghost_values_finish(const unsigned int component_in_block_vector,
3327  const VectorType & vec)
3328  {
3329  (void)component_in_block_vector;
3330  vec.update_ghost_values_finish();
3331  }
3332 
3333 
3334 
3341  template <typename VectorType,
3342  typename std::enable_if<
3343  has_update_ghost_values_start<VectorType>::value &&
3344  has_exchange_on_subset<VectorType>::value,
3345  VectorType>::type * = nullptr>
3346  void
3347  update_ghost_values_finish(const unsigned int component_in_block_vector,
3348  const VectorType & vec)
3349  {
3350  static_assert(
3351  std::is_same<Number, typename VectorType::value_type>::value,
3352  "Type mismatch between VectorType and VectorDataExchange");
3353  (void)component_in_block_vector;
3354  if (vector_face_access ==
3357  vec.size() == 0)
3358  vec.update_ghost_values_finish();
3359  else
3360  {
3361 # ifdef DEAL_II_WITH_MPI
3362 
3363  AssertIndexRange(component_in_block_vector, tmp_data.size());
3364  AssertDimension(requests.size(), tmp_data.size());
3365 
3366  const unsigned int mf_component = find_vector_in_mf(vec);
3367  const Utilities::MPI::Partitioner &part =
3368  get_partitioner(mf_component);
3369  if (&part ==
3370  matrix_free.get_dof_info(mf_component).vector_partitioner.get())
3371  {
3372  vec.update_ghost_values_finish();
3373  return;
3374  }
3375 
3376  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0)
3377  return;
3378 
3380  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3381  vec.get_partitioner()->local_size(),
3382  vec.get_partitioner()->n_ghost_indices()),
3383  this->requests[component_in_block_vector]);
3384 
3385  matrix_free.release_scratch_data_non_threadsafe(
3386  tmp_data[component_in_block_vector]);
3387  tmp_data[component_in_block_vector] = nullptr;
3388 
3389  // let vector know that ghosts are being updated and we can read from
3390  // them
3391  vec.set_ghost_state(true);
3392 # endif
3393  }
3394  }
3395 
3396 
3397 
3401  template <typename VectorType,
3402  typename std::enable_if<is_serial_or_dummy<VectorType>::value,
3403  VectorType>::type * = nullptr>
3404  void
3405  compress_start(const unsigned int /*component_in_block_vector*/,
3406  VectorType & /*vec*/)
3407  {}
3408 
3409 
3410 
3415  template <typename VectorType,
3416  typename std::enable_if<!has_compress_start<VectorType>::value &&
3417  !is_serial_or_dummy<VectorType>::value,
3418  VectorType>::type * = nullptr>
3419  void
3420  compress_start(const unsigned int component_in_block_vector,
3421  VectorType & vec)
3422  {
3423  (void)component_in_block_vector;
3424  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3425  vec.compress(::VectorOperation::add);
3426  }
3427 
3428 
3429 
3435  template <
3436  typename VectorType,
3437  typename std::enable_if<has_compress_start<VectorType>::value &&
3438  !has_exchange_on_subset<VectorType>::value,
3439  VectorType>::type * = nullptr>
3440  void
3441  compress_start(const unsigned int component_in_block_vector,
3442  VectorType & vec)
3443  {
3444  (void)component_in_block_vector;
3445  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3446  vec.compress_start(component_in_block_vector + channel_shift);
3447  }
3448 
3449 
3450 
3457  template <
3458  typename VectorType,
3459  typename std::enable_if<has_compress_start<VectorType>::value &&
3460  has_exchange_on_subset<VectorType>::value,
3461  VectorType>::type * = nullptr>
3462  void
3463  compress_start(const unsigned int component_in_block_vector,
3464  VectorType & vec)
3465  {
3466  static_assert(
3467  std::is_same<Number, typename VectorType::value_type>::value,
3468  "Type mismatch between VectorType and VectorDataExchange");
3469  (void)component_in_block_vector;
3470  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3471  if (vector_face_access ==
3474  vec.size() == 0)
3475  vec.compress_start(component_in_block_vector + channel_shift);
3476  else
3477  {
3478 # ifdef DEAL_II_WITH_MPI
3479 
3480  const unsigned int mf_component = find_vector_in_mf(vec);
3481  const Utilities::MPI::Partitioner &part =
3482  get_partitioner(mf_component);
3483  if (&part ==
3484  matrix_free.get_dof_info(mf_component).vector_partitioner.get())
3485  {
3486  vec.compress_start(component_in_block_vector + channel_shift);
3487  return;
3488  }
3489 
3490  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0)
3491  return;
3492 
3493  tmp_data[component_in_block_vector] =
3494  matrix_free.acquire_scratch_data_non_threadsafe();
3495  tmp_data[component_in_block_vector]->resize_fast(
3496  part.n_import_indices());
3497  AssertDimension(requests.size(), tmp_data.size());
3498 
3501  component_in_block_vector + channel_shift,
3502  ArrayView<Number>(vec.begin() + vec.get_partitioner()->local_size(),
3503  vec.get_partitioner()->n_ghost_indices()),
3504  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3505  part.n_import_indices()),
3506  this->requests[component_in_block_vector]);
3507 # endif
3508  }
3509  }
3510 
3511 
3512 
3517  template <typename VectorType,
3518  typename std::enable_if<!has_compress_start<VectorType>::value,
3519  VectorType>::type * = nullptr>
3520  void
3521  compress_finish(const unsigned int /*component_in_block_vector*/,
3522  VectorType & /*vec*/)
3523  {}
3524 
3525 
3526 
3532  template <
3533  typename VectorType,
3534  typename std::enable_if<has_compress_start<VectorType>::value &&
3535  !has_exchange_on_subset<VectorType>::value,
3536  VectorType>::type * = nullptr>
3537  void
3538  compress_finish(const unsigned int component_in_block_vector,
3539  VectorType & vec)
3540  {
3541  (void)component_in_block_vector;
3542  vec.compress_finish(::VectorOperation::add);
3543  }
3544 
3545 
3546 
3553  template <
3554  typename VectorType,
3555  typename std::enable_if<has_compress_start<VectorType>::value &&
3556  has_exchange_on_subset<VectorType>::value,
3557  VectorType>::type * = nullptr>
3558  void
3559  compress_finish(const unsigned int component_in_block_vector,
3560  VectorType & vec)
3561  {
3562  static_assert(
3563  std::is_same<Number, typename VectorType::value_type>::value,
3564  "Type mismatch between VectorType and VectorDataExchange");
3565  (void)component_in_block_vector;
3566  if (vector_face_access ==
3569  vec.size() == 0)
3570  vec.compress_finish(::VectorOperation::add);
3571  else
3572  {
3573 # ifdef DEAL_II_WITH_MPI
3574  AssertIndexRange(component_in_block_vector, tmp_data.size());
3575  AssertDimension(requests.size(), tmp_data.size());
3576 
3577  const unsigned int mf_component = find_vector_in_mf(vec);
3578 
3579  const Utilities::MPI::Partitioner &part =
3580  get_partitioner(mf_component);
3581  if (&part ==
3582  matrix_free.get_dof_info(mf_component).vector_partitioner.get())
3583  {
3584  vec.compress_finish(::VectorOperation::add);
3585  return;
3586  }
3587 
3588  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0)
3589  return;
3590 
3594  tmp_data[component_in_block_vector]->begin(),
3595  part.n_import_indices()),
3596  ArrayView<Number>(vec.begin(), part.local_size()),
3597  ArrayView<Number>(vec.begin() + vec.get_partitioner()->local_size(),
3598  vec.get_partitioner()->n_ghost_indices()),
3599  this->requests[component_in_block_vector]);
3600 
3601  matrix_free.release_scratch_data_non_threadsafe(
3602  tmp_data[component_in_block_vector]);
3603  tmp_data[component_in_block_vector] = nullptr;
3604 # endif
3605  }
3606  }
3607 
3608 
3609 
3613  template <typename VectorType,
3614  typename std::enable_if<is_serial_or_dummy<VectorType>::value,
3615  VectorType>::type * = nullptr>
3616  void
3617  reset_ghost_values(const VectorType & /*vec*/) const
3618  {}
3619 
3620 
3621 
3626  template <
3627  typename VectorType,
3628  typename std::enable_if<!has_exchange_on_subset<VectorType>::value &&
3629  !is_serial_or_dummy<VectorType>::value,
3630  VectorType>::type * = nullptr>
3631  void
3632  reset_ghost_values(const VectorType &vec) const
3633  {
3634  if (ghosts_were_set == true)
3635  return;
3636 
3637  vec.zero_out_ghosts();
3638  }
3639 
3640 
3641 
3647  template <typename VectorType,
3648  typename std::enable_if<has_exchange_on_subset<VectorType>::value,
3649  VectorType>::type * = nullptr>
3650  void
3651  reset_ghost_values(const VectorType &vec) const
3652  {
3653  static_assert(
3654  std::is_same<Number, typename VectorType::value_type>::value,
3655  "Type mismatch between VectorType and VectorDataExchange");
3656  if (ghosts_were_set == true)
3657  return;
3658 
3659  if (vector_face_access ==
3662  vec.size() == 0)
3663  vec.zero_out_ghosts();
3664  else
3665  {
3666 # ifdef DEAL_II_WITH_MPI
3667  AssertDimension(requests.size(), tmp_data.size());
3668 
3669  const unsigned int mf_component = find_vector_in_mf(vec);
3670  const Utilities::MPI::Partitioner &part =
3671  get_partitioner(mf_component);
3672  if (&part ==
3673  matrix_free.get_dof_info(mf_component).vector_partitioner.get())
3674  vec.zero_out_ghosts();
3675  else if (part.n_ghost_indices() > 0)
3676  {
3677  for (std::vector<std::pair<unsigned int, unsigned int>>::
3678  const_iterator my_ghosts =
3680  my_ghosts !=
3682  ++my_ghosts)
3683  for (unsigned int j = my_ghosts->first; j < my_ghosts->second;
3684  j++)
3685  {
3687  vec)
3688  .local_element(j + part.local_size()) = 0.;
3689  }
3690  }
3691 
3692  // let vector know that it's not ghosted anymore
3693  vec.set_ghost_state(false);
3694 # endif
3695  }
3696  }
3697 
3698 
3699 
3705  template <typename VectorType,
3706  typename std::enable_if<has_exchange_on_subset<VectorType>::value,
3707  VectorType>::type * = nullptr>
3708  void
3709  zero_vector_region(const unsigned int range_index, VectorType &vec) const
3710  {
3711  static_assert(
3712  std::is_same<Number, typename VectorType::value_type>::value,
3713  "Type mismatch between VectorType and VectorDataExchange");
3714  if (range_index == numbers::invalid_unsigned_int)
3715  vec = Number();
3716  else
3717  {
3718  const unsigned int mf_component = find_vector_in_mf(vec, false);
3720  matrix_free.get_dof_info(mf_component);
3721  Assert(dof_info.vector_zero_range_list_index.empty() == false,
3722  ExcNotInitialized());
3723 
3724  Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3725  ExcInternalError());
3726  AssertIndexRange(range_index,
3727  dof_info.vector_zero_range_list_index.size() - 1);
3728  for (unsigned int id =
3729  dof_info.vector_zero_range_list_index[range_index];
3730  id != dof_info.vector_zero_range_list_index[range_index + 1];
3731  ++id)
3732  std::memset(vec.begin() + dof_info.vector_zero_range_list[id].first,
3733  0,
3734  (dof_info.vector_zero_range_list[id].second -
3735  dof_info.vector_zero_range_list[id].first) *
3736  sizeof(Number));
3737  }
3738  }
3739 
3740 
3741 
3747  template <
3748  typename VectorType,
3749  typename std::enable_if<!has_exchange_on_subset<VectorType>::value,
3750  VectorType>::type * = nullptr,
3751  typename VectorType::value_type * = nullptr>
3752  void
3753  zero_vector_region(const unsigned int range_index, VectorType &vec) const
3754  {
3755  if (range_index == numbers::invalid_unsigned_int)
3756  vec = typename VectorType::value_type();
3757  else
3758  {
3759  Assert(false, ExcNotImplemented());
3760  }
3761  }
3762 
3763 
3764 
3769  void
3770  zero_vector_region(const unsigned int, ...) const
3771  {
3772  Assert(false,
3773  ExcNotImplemented("Zeroing is only implemented for vector types "
3774  "which provide VectorType::value_type"));
3775  }
3776 
3777 
3778 
3779  const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
3780  const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3781  DataAccessOnFaces vector_face_access;
3782  bool ghosts_were_set;
3783 # ifdef DEAL_II_WITH_MPI
3784  std::vector<AlignedVector<Number> *> tmp_data;
3785  std::vector<std::vector<MPI_Request>> requests;
3786 # endif
3787  }; // VectorDataExchange
3788 
3789  template <typename VectorStruct>
3790  unsigned int
3791  n_components(const VectorStruct &vec);
3792 
3793  template <typename VectorStruct>
3794  unsigned int
3795  n_components_block(const VectorStruct &vec,
3796  std::integral_constant<bool, true>)
3797  {
3798  unsigned int components = 0;
3799  for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
3800  components += n_components(vec.block(bl));
3801  return components;
3802  }
3803 
3804  template <typename VectorStruct>
3805  unsigned int
3806  n_components_block(const VectorStruct &, std::integral_constant<bool, false>)
3807  {
3808  return 1;
3809  }
3810 
3811  template <typename VectorStruct>
3812  unsigned int
3813  n_components(const VectorStruct &vec)
3814  {
3815  return n_components_block(
3816  vec, std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3817  }
3818 
3819  template <typename VectorStruct>
3820  inline unsigned int
3821  n_components(const std::vector<VectorStruct> &vec)
3822  {
3823  unsigned int components = 0;
3824  for (unsigned int comp = 0; comp < vec.size(); comp++)
3825  components += n_components_block(
3826  vec[comp],
3827  std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3828  return components;
3829  }
3830 
3831  template <typename VectorStruct>
3832  inline unsigned int
3833  n_components(const std::vector<VectorStruct *> &vec)
3834  {
3835  unsigned int components = 0;
3836  for (unsigned int comp = 0; comp < vec.size(); comp++)
3837  components += n_components_block(
3838  *vec[comp],
3839  std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3840  return components;
3841  }
3842 
3843 
3844 
3845  // A helper function to identify block vectors with many components where we
3846  // should not try to overlap computations and communication because there
3847  // would be too many outstanding communication requests.
3848 
3849  // default value for vectors that do not have communication_block_size
3850  template <
3851  typename VectorStruct,
3852  typename std::enable_if<!has_communication_block_size<VectorStruct>::value,
3853  VectorStruct>::type * = nullptr>
3854  constexpr unsigned int
3855  get_communication_block_size(const VectorStruct &)
3856  {
3858  }
3859 
3860 
3861 
3862  template <
3863  typename VectorStruct,
3864  typename std::enable_if<has_communication_block_size<VectorStruct>::value,
3865  VectorStruct>::type * = nullptr>
3866  constexpr unsigned int
3867  get_communication_block_size(const VectorStruct &)
3868  {
3869  return VectorStruct::communication_block_size;
3870  }
3871 
3872 
3873 
3874  // --------------------------------------------------------------------------
3875  // below we have wrappers to distinguish between block and non-block vectors.
3876  // --------------------------------------------------------------------------
3877 
3878  //
3879  // update_ghost_values_start
3880  //
3881 
3882  // update_ghost_values for block vectors
3883  template <int dim,
3884  typename VectorStruct,
3885  typename Number,
3886  typename VectorizedArrayType,
3887  typename std::enable_if<IsBlockVector<VectorStruct>::value,
3888  VectorStruct>::type * = nullptr>
3889  void
3890  update_ghost_values_start(
3891  const VectorStruct & vec,
3892  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3893  const unsigned int channel = 0)
3894  {
3895  if (get_communication_block_size(vec) < vec.n_blocks())
3896  {
3897  // don't forget to set ghosts_were_set, that otherwise happens
3898  // inside VectorDataExchange::update_ghost_values_start()
3899  exchanger.ghosts_were_set = vec.has_ghost_elements();
3900  vec.update_ghost_values();
3901  }
3902  else
3903  {
3904  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
3905  update_ghost_values_start(vec.block(i), exchanger, channel + i);
3906  }
3907  }
3908 
3909 
3910 
3911  // update_ghost_values for non-block vectors
3912  template <int dim,
3913  typename VectorStruct,
3914  typename Number,
3915  typename VectorizedArrayType,
3916  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
3917  VectorStruct>::type * = nullptr>
3918  void
3919  update_ghost_values_start(
3920  const VectorStruct & vec,
3921  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3922  const unsigned int channel = 0)
3923  {
3924  exchanger.update_ghost_values_start(channel, vec);
3925  }
3926 
3927 
3928 
3929  // update_ghost_values_start() for vector of vectors
3930  template <int dim,
3931  typename VectorStruct,
3932  typename Number,
3933  typename VectorizedArrayType>
3934  inline void
3935  update_ghost_values_start(
3936  const std::vector<VectorStruct> & vec,
3937  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3938  {
3939  unsigned int component_index = 0;
3940  for (unsigned int comp = 0; comp < vec.size(); comp++)
3941  {
3942  update_ghost_values_start(vec[comp], exchanger, component_index);
3943  component_index += n_components(vec[comp]);
3944  }
3945  }
3946 
3947 
3948 
3949  // update_ghost_values_start() for vector of pointers to vectors
3950  template <int dim,
3951  typename VectorStruct,
3952  typename Number,
3953  typename VectorizedArrayType>
3954  inline void
3955  update_ghost_values_start(
3956  const std::vector<VectorStruct *> & vec,
3957  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3958  {
3959  unsigned int component_index = 0;
3960  for (unsigned int comp = 0; comp < vec.size(); comp++)
3961  {
3962  update_ghost_values_start(*vec[comp], exchanger, component_index);
3963  component_index += n_components(*vec[comp]);
3964  }
3965  }
3966 
3967 
3968 
3969  //
3970  // update_ghost_values_finish
3971  //
3972 
3973  // for block vectors
3974  template <int dim,
3975  typename VectorStruct,
3976  typename Number,
3977  typename VectorizedArrayType,
3978  typename std::enable_if<IsBlockVector<VectorStruct>::value,
3979  VectorStruct>::type * = nullptr>
3980  void
3981  update_ghost_values_finish(
3982  const VectorStruct & vec,
3983  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3984  const unsigned int channel = 0)
3985  {
3986  if (get_communication_block_size(vec) < vec.n_blocks())
3987  {
3988  // do nothing, everything has already been completed in the _start()
3989  // call
3990  }
3991  else
3992  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
3993  update_ghost_values_finish(vec.block(i), exchanger, channel + i);
3994  }
3995 
3996 
3997 
3998  // for non-block vectors
3999  template <int dim,
4000  typename VectorStruct,
4001  typename Number,
4002  typename VectorizedArrayType,
4003  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4004  VectorStruct>::type * = nullptr>
4005  void
4006  update_ghost_values_finish(
4007  const VectorStruct & vec,
4008  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4009  const unsigned int channel = 0)
4010  {
4011  exchanger.update_ghost_values_finish(channel, vec);
4012  }
4013 
4014 
4015 
4016  // for vector of vectors
4017  template <int dim,
4018  typename VectorStruct,
4019  typename Number,
4020  typename VectorizedArrayType>
4021  inline void
4022  update_ghost_values_finish(
4023  const std::vector<VectorStruct> & vec,
4024  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4025  {
4026  unsigned int component_index = 0;
4027  for (unsigned int comp = 0; comp < vec.size(); comp++)
4028  {
4029  update_ghost_values_finish(vec[comp], exchanger, component_index);
4030  component_index += n_components(vec[comp]);
4031  }
4032  }
4033 
4034 
4035 
4036  // for vector of pointers to vectors
4037  template <int dim,
4038  typename VectorStruct,
4039  typename Number,
4040  typename VectorizedArrayType>
4041  inline void
4042  update_ghost_values_finish(
4043  const std::vector<VectorStruct *> & vec,
4044  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4045  {
4046  unsigned int component_index = 0;
4047  for (unsigned int comp = 0; comp < vec.size(); comp++)
4048  {
4049  update_ghost_values_finish(*vec[comp], exchanger, component_index);
4050  component_index += n_components(*vec[comp]);
4051  }
4052  }
4053 
4054 
4055 
4056  //
4057  // compress_start
4058  //
4059 
4060  // for block vectors
4061  template <int dim,
4062  typename VectorStruct,
4063  typename Number,
4064  typename VectorizedArrayType,
4065  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4066  VectorStruct>::type * = nullptr>
4067  inline void
4068  compress_start(
4069  VectorStruct & vec,
4070  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4071  const unsigned int channel = 0)
4072  {
4073  if (get_communication_block_size(vec) < vec.n_blocks())
4074  vec.compress(::VectorOperation::add);
4075  else
4076  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4077  compress_start(vec.block(i), exchanger, channel + i);
4078  }
4079 
4080 
4081 
4082  // for non-block vectors
4083  template <int dim,
4084  typename VectorStruct,
4085  typename Number,
4086  typename VectorizedArrayType,
4087  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4088  VectorStruct>::type * = nullptr>
4089  inline void
4090  compress_start(
4091  VectorStruct & vec,
4092  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4093  const unsigned int channel = 0)
4094  {
4095  exchanger.compress_start(channel, vec);
4096  }
4097 
4098 
4099 
4100  // for std::vector of vectors
4101  template <int dim,
4102  typename VectorStruct,
4103  typename Number,
4104  typename VectorizedArrayType>
4105  inline void
4106  compress_start(
4107  std::vector<VectorStruct> & vec,
4108  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4109  {
4110  unsigned int component_index = 0;
4111  for (unsigned int comp = 0; comp < vec.size(); comp++)
4112  {
4113  compress_start(vec[comp], exchanger, component_index);
4114  component_index += n_components(vec[comp]);
4115  }
4116  }
4117 
4118 
4119 
4120  // for std::vector of pointer to vectors
4121  template <int dim,
4122  typename VectorStruct,
4123  typename Number,
4124  typename VectorizedArrayType>
4125  inline void
4126  compress_start(
4127  std::vector<VectorStruct *> & vec,
4128  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4129  {
4130  unsigned int component_index = 0;
4131  for (unsigned int comp = 0; comp < vec.size(); comp++)
4132  {
4133  compress_start(*vec[comp], exchanger, component_index);
4134  component_index += n_components(*vec[comp]);
4135  }
4136  }
4137 
4138 
4139 
4140  //
4141  // compress_finish
4142  //
4143 
4144  // for block vectors
4145  template <int dim,
4146  typename VectorStruct,
4147  typename Number,
4148  typename VectorizedArrayType,
4149  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4150  VectorStruct>::type * = nullptr>
4151  inline void
4152  compress_finish(
4153  VectorStruct & vec,
4154  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4155  const unsigned int channel = 0)
4156  {
4157  if (get_communication_block_size(vec) < vec.n_blocks())
4158  {
4159  // do nothing, everything has already been completed in the _start()
4160  // call
4161  }
4162  else
4163  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4164  compress_finish(vec.block(i), exchanger, channel + i);
4165  }
4166 
4167 
4168 
4169  // for non-block vectors
4170  template <int dim,
4171  typename VectorStruct,
4172  typename Number,
4173  typename VectorizedArrayType,
4174  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4175  VectorStruct>::type * = nullptr>
4176  inline void
4177  compress_finish(
4178  VectorStruct & vec,
4179  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4180  const unsigned int channel = 0)
4181  {
4182  exchanger.compress_finish(channel, vec);
4183  }
4184 
4185 
4186 
4187  // for std::vector of vectors
4188  template <int dim,
4189  typename VectorStruct,
4190  typename Number,
4191  typename VectorizedArrayType>
4192  inline void
4193  compress_finish(
4194  std::vector<VectorStruct> & vec,
4195  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4196  {
4197  unsigned int component_index = 0;
4198  for (unsigned int comp = 0; comp < vec.size(); comp++)
4199  {
4200  compress_finish(vec[comp], exchanger, component_index);
4201  component_index += n_components(vec[comp]);
4202  }
4203  }
4204 
4205 
4206 
4207  // for std::vector of pointer to vectors
4208  template <int dim,
4209  typename VectorStruct,
4210  typename Number,
4211  typename VectorizedArrayType>
4212  inline void
4213  compress_finish(
4214  std::vector<VectorStruct *> & vec,
4215  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4216  {
4217  unsigned int component_index = 0;
4218  for (unsigned int comp = 0; comp < vec.size(); comp++)
4219  {
4220  compress_finish(*vec[comp], exchanger, component_index);
4221  component_index += n_components(*vec[comp]);
4222  }
4223  }
4224 
4225 
4226 
4227  //
4228  // reset_ghost_values:
4229  //
4230  // if the input vector did not have ghosts imported, clear them here again
4231  // in order to avoid subsequent operations e.g. in linear solvers to work
4232  // with ghosts all the time
4233  //
4234 
4235  // for block vectors
4236  template <int dim,
4237  typename VectorStruct,
4238  typename Number,
4239  typename VectorizedArrayType,
4240  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4241  VectorStruct>::type * = nullptr>
4242  inline void
4243  reset_ghost_values(
4244  const VectorStruct & vec,
4245  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4246  {
4247  // return immediately if there is nothing to do.
4248  if (exchanger.ghosts_were_set == true)
4249  return;
4250 
4251  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4252  reset_ghost_values(vec.block(i), exchanger);
4253  }
4254 
4255 
4256 
4257  // for non-block vectors
4258  template <int dim,
4259  typename VectorStruct,
4260  typename Number,
4261  typename VectorizedArrayType,
4262  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4263  VectorStruct>::type * = nullptr>
4264  inline void
4265  reset_ghost_values(
4266  const VectorStruct & vec,
4267  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4268  {
4269  exchanger.reset_ghost_values(vec);
4270  }
4271 
4272 
4273 
4274  // for std::vector of vectors
4275  template <int dim,
4276  typename VectorStruct,
4277  typename Number,
4278  typename VectorizedArrayType>
4279  inline void
4280  reset_ghost_values(
4281  const std::vector<VectorStruct> & vec,
4282  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4283  {
4284  // return immediately if there is nothing to do.
4285  if (exchanger.ghosts_were_set == true)
4286  return;
4287 
4288  for (unsigned int comp = 0; comp < vec.size(); comp++)
4289  reset_ghost_values(vec[comp], exchanger);
4290  }
4291 
4292 
4293 
4294  // for std::vector of pointer to vectors
4295  template <int dim,
4296  typename VectorStruct,
4297  typename Number,
4298  typename VectorizedArrayType>
4299  inline void
4300  reset_ghost_values(
4301  const std::vector<VectorStruct *> & vec,
4302  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4303  {
4304  // return immediately if there is nothing to do.
4305  if (exchanger.ghosts_were_set == true)
4306  return;
4307 
4308  for (unsigned int comp = 0; comp < vec.size(); comp++)
4309  reset_ghost_values(*vec[comp], exchanger);
4310  }
4311 
4312 
4313 
4314  //
4315  // zero_vector_region
4316  //
4317 
4318  // for block vectors
4319  template <int dim,
4320  typename VectorStruct,
4321  typename Number,
4322  typename VectorizedArrayType,
4323  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4324  VectorStruct>::type * = nullptr>
4325  inline void
4326  zero_vector_region(
4327  const unsigned int range_index,
4328  VectorStruct & vec,
4329  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4330  {
4331  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4332  exchanger.zero_vector_region(range_index, vec.block(i));
4333  }
4334 
4335 
4336 
4337  // for non-block vectors
4338  template <int dim,
4339  typename VectorStruct,
4340  typename Number,
4341  typename VectorizedArrayType,
4342  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4343  VectorStruct>::type * = nullptr>
4344  inline void
4345  zero_vector_region(
4346  const unsigned int range_index,
4347  VectorStruct & vec,
4348  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4349  {
4350  exchanger.zero_vector_region(range_index, vec);
4351  }
4352 
4353 
4354 
4355  // for std::vector of vectors
4356  template <int dim,
4357  typename VectorStruct,
4358  typename Number,
4359  typename VectorizedArrayType>
4360  inline void
4361  zero_vector_region(
4362  const unsigned int range_index,
4363  std::vector<VectorStruct> & vec,
4364  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4365  {
4366  for (unsigned int comp = 0; comp < vec.size(); comp++)
4367  zero_vector_region(range_index, vec[comp], exchanger);
4368  }
4369 
4370 
4371 
4372  // for std::vector of pointers to vectors
4373  template <int dim,
4374  typename VectorStruct,
4375  typename Number,
4376  typename VectorizedArrayType>
4377  inline void
4378  zero_vector_region(
4379  const unsigned int range_index,
4380  std::vector<VectorStruct *> & vec,
4381  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4382  {
4383  for (unsigned int comp = 0; comp < vec.size(); comp++)
4384  zero_vector_region(range_index, *vec[comp], exchanger);
4385  }
4386 
4387 
4388 
4389  namespace MatrixFreeFunctions
4390  {
4391  // struct to select between a const interface and a non-const interface
4392  // for MFWorker
4393  template <typename, typename, typename, typename, bool>
4394  struct InterfaceSelector
4395  {};
4396 
4397  // Version for constant functions
4398  template <typename MF,
4399  typename InVector,
4400  typename OutVector,
4401  typename Container>
4402  struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4403  {
4404  using function_type = void (Container::*)(
4405  const MF &,
4406  OutVector &,
4407  const InVector &,
4408  const std::pair<unsigned int, unsigned int> &) const;
4409  };
4410 
4411  // Version for non-constant functions
4412  template <typename MF,
4413  typename InVector,
4414  typename OutVector,
4415  typename Container>
4416  struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4417  {
4418  using function_type =
4419  void (Container::*)(const MF &,
4420  OutVector &,
4421  const InVector &,
4422  const std::pair<unsigned int, unsigned int> &);
4423  };
4424  } // namespace MatrixFreeFunctions
4425 
4426 
4427 
4428  // A implementation class for the worker object that runs the various
4429  // operations we want to perform during the matrix-free loop
4430  template <typename MF,
4431  typename InVector,
4432  typename OutVector,
4433  typename Container,
4434  bool is_constant>
4435  class MFWorker : public MFWorkerInterface
4436  {
4437  public:
4438  // An alias to make the arguments further down more readable
4439  using function_type = typename MatrixFreeFunctions::
4440  InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4441  function_type;
4442 
4443  // constructor, binds all the arguments to this class
4444  MFWorker(const MF & matrix_free,
4445  const InVector & src,
4446  OutVector & dst,
4447  const bool zero_dst_vector_setting,
4448  const Container & container,
4449  function_type cell_function,
4450  function_type face_function,
4451  function_type boundary_function,
4452  const typename MF::DataAccessOnFaces src_vector_face_access =
4453  MF::DataAccessOnFaces::none,
4454  const typename MF::DataAccessOnFaces dst_vector_face_access =
4455  MF::DataAccessOnFaces::none,
4456  const std::function<void(const unsigned int, const unsigned int)>
4457  &operation_before_loop = {},
4458  const std::function<void(const unsigned int, const unsigned int)>
4459  & operation_after_loop = {},
4460  const unsigned int dof_handler_index_pre_post = 0)
4461  : matrix_free(matrix_free)
4462  , container(const_cast<Container &>(container))
4463  , cell_function(cell_function)
4464  , face_function(face_function)
4465  , boundary_function(boundary_function)
4466  , src(src)
4467  , dst(dst)
4468  , src_data_exchanger(matrix_free,
4469  src_vector_face_access,
4470  n_components(src))
4471  , dst_data_exchanger(matrix_free,
4472  dst_vector_face_access,
4473  n_components(dst))
4474  , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4475  , zero_dst_vector_setting(zero_dst_vector_setting &&
4476  !src_and_dst_are_same)
4477  , operation_before_loop(operation_before_loop)
4478  , operation_after_loop(operation_after_loop)
4479  , dof_handler_index_pre_post(dof_handler_index_pre_post)
4480  {}
4481 
4482  // Runs the cell work. If no function is given, nothing is done
4483  virtual void
4484  cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4485  {
4486  if (cell_function != nullptr && cell_range.second > cell_range.first)
4487  (container.*
4488  cell_function)(matrix_free, this->dst, this->src, cell_range);
4489  }
4490 
4491  // Runs the assembler on interior faces. If no function is given, nothing
4492  // is done
4493  virtual void
4494  face(const std::pair<unsigned int, unsigned int> &face_range) override
4495  {
4496  if (face_function != nullptr && face_range.second > face_range.first)
4497  (container.*
4498  face_function)(matrix_free, this->dst, this->src, face_range);
4499  }
4500 
4501  // Runs the assembler on boundary faces. If no function is given, nothing
4502  // is done
4503  virtual void
4504  boundary(const std::pair<unsigned int, unsigned int> &face_range) override
4505  {
4506  if (boundary_function != nullptr && face_range.second > face_range.first)
4507  (container.*
4508  boundary_function)(matrix_free, this->dst, this->src, face_range);
4509  }
4510 
4511  // Starts the communication for the update ghost values operation. We
4512  // cannot call this update if ghost and destination are the same because
4513  // that would introduce spurious entries in the destination (there is also
4514  // the problem that reading from a vector that we also write to is usually
4515  // not intended in case there is overlap, but this is up to the
4516  // application code to decide and we cannot catch this case here).
4517  virtual void
4518  vector_update_ghosts_start() override
4519  {
4520  if (!src_and_dst_are_same)
4521  internal::update_ghost_values_start(src, src_data_exchanger);
4522  }
4523 
4524  // Finishes the communication for the update ghost values operation
4525  virtual void
4526  vector_update_ghosts_finish() override
4527  {
4528  if (!src_and_dst_are_same)
4529  internal::update_ghost_values_finish(src, src_data_exchanger);
4530  }
4531 
4532  // Starts the communication for the vector compress operation
4533  virtual void
4534  vector_compress_start() override
4535  {
4536  internal::compress_start(dst, dst_data_exchanger);
4537  }
4538 
4539  // Finishes the communication for the vector compress operation
4540  virtual void
4541  vector_compress_finish() override
4542  {
4543  internal::compress_finish(dst, dst_data_exchanger);
4544  if (!src_and_dst_are_same)
4545  internal::reset_ghost_values(src, src_data_exchanger);
4546  }
4547 
4548  // Zeros the given input vector
4549  virtual void
4550  zero_dst_vector_range(const unsigned int range_index) override
4551  {
4552  if (zero_dst_vector_setting)
4553  internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4554  }
4555 
4556  virtual void
4557  cell_loop_pre_range(const unsigned int range_index) override
4558  {
4559  if (operation_before_loop)
4560  {
4562  matrix_free.get_dof_info(dof_handler_index_pre_post);
4563  AssertIndexRange(range_index,
4564  dof_info.cell_loop_pre_list_index.size() - 1);
4565  for (unsigned int id = dof_info.cell_loop_pre_list_index[range_index];
4566  id != dof_info.cell_loop_pre_list_index[range_index + 1];
4567  ++id)
4568  operation_before_loop(dof_info.cell_loop_pre_list[id].first,
4569  dof_info.cell_loop_pre_list[id].second);
4570  }
4571  }
4572 
4573  virtual void
4574  cell_loop_post_range(const unsigned int range_index) override
4575  {
4576  if (operation_after_loop)
4577  {
4579  matrix_free.get_dof_info(dof_handler_index_pre_post);
4580  AssertIndexRange(range_index,
4581  dof_info.cell_loop_post_list_index.size() - 1);
4582  for (unsigned int id =
4583  dof_info.cell_loop_post_list_index[range_index];
4584  id != dof_info.cell_loop_post_list_index[range_index + 1];
4585  ++id)
4586  operation_after_loop(dof_info.cell_loop_post_list[id].first,
4587  dof_info.cell_loop_post_list[id].second);
4588  }
4589  }
4590 
4591  private:
4592  const MF & matrix_free;
4593  Container & container;
4594  function_type cell_function;
4595  function_type face_function;
4596  function_type boundary_function;
4597 
4598  const InVector &src;
4599  OutVector & dst;
4600  VectorDataExchange<MF::dimension,
4601  typename MF::value_type,
4602  typename MF::vectorized_value_type>
4603  src_data_exchanger;
4604  VectorDataExchange<MF::dimension,
4605  typename MF::value_type,
4606  typename MF::vectorized_value_type>
4607  dst_data_exchanger;
4608  const bool src_and_dst_are_same;
4609  const bool zero_dst_vector_setting;
4610  const std::function<void(const unsigned int, const unsigned int)>
4611  operation_before_loop;
4612  const std::function<void(const unsigned int, const unsigned int)>
4613  operation_after_loop;
4614  const unsigned int dof_handler_index_pre_post;
4615  };
4616 
4617 
4618 
4623  template <class MF, typename InVector, typename OutVector>
4624  struct MFClassWrapper
4625  {
4626  using function_type =
4627  std::function<void(const MF &,
4628  OutVector &,
4629  const InVector &,
4630  const std::pair<unsigned int, unsigned int> &)>;
4631 
4632  MFClassWrapper(const function_type cell,
4633  const function_type face,
4634  const function_type boundary)
4635  : cell(cell)
4636  , face(face)
4637  , boundary(boundary)
4638  {}
4639 
4640  void
4641  cell_integrator(const MF & mf,
4642  OutVector & dst,
4643  const InVector & src,
4644  const std::pair<unsigned int, unsigned int> &range) const
4645  {
4646  if (cell)
4647  cell(mf, dst, src, range);
4648  }
4649 
4650  void
4651  face_integrator(const MF & mf,
4652  OutVector & dst,
4653  const InVector & src,
4654  const std::pair<unsigned int, unsigned int> &range) const
4655  {
4656  if (face)
4657  face(mf, dst, src, range);
4658  }
4659 
4660  void
4661  boundary_integrator(
4662  const MF & mf,
4663  OutVector & dst,
4664  const InVector & src,
4665  const std::pair<unsigned int, unsigned int> &range) const
4666  {
4667  if (boundary)
4668  boundary(mf, dst, src, range);
4669  }
4670 
4671  const function_type cell;
4672  const function_type face;
4673  const function_type boundary;
4674  };
4675 
4676 } // end of namespace internal
4677 
4678 
4679 
4680 template <int dim, typename Number, typename VectorizedArrayType>
4681 template <typename OutVector, typename InVector>
4682 inline void
4684  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4685  OutVector &,
4686  const InVector &,
4687  const std::pair<unsigned int, unsigned int> &)>
4688  & cell_operation,
4689  OutVector & dst,
4690  const InVector &src,
4691  const bool zero_dst_vector) const
4692 {
4693  using Wrapper =
4694  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4695  InVector,
4696  OutVector>;
4697  Wrapper wrap(cell_operation, nullptr, nullptr);
4698  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4699  InVector,
4700  OutVector,
4701  Wrapper,
4702  true>
4703  worker(*this,
4704  src,
4705  dst,
4706  zero_dst_vector,
4707  wrap,
4708  &Wrapper::cell_integrator,
4709  &Wrapper::face_integrator,
4710  &Wrapper::boundary_integrator);
4711 
4712  task_info.loop(worker);
4713 }
4714 
4715 
4716 
4717 template <int dim, typename Number, typename VectorizedArrayType>
4718 template <typename OutVector, typename InVector>
4719 inline void
4721  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4722  OutVector &,
4723  const InVector &,
4724  const std::pair<unsigned int, unsigned int> &)>
4725  & cell_operation,
4726  OutVector & dst,
4727  const InVector &src,
4728  const std::function<void(const unsigned int, const unsigned int)>
4729  &operation_before_loop,
4730  const std::function<void(const unsigned int, const unsigned int)>
4731  & operation_after_loop,
4732  const unsigned int dof_handler_index_pre_post) const
4733 {
4734  using Wrapper =
4735  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4736  InVector,
4737  OutVector>;
4738  Wrapper wrap(cell_operation, nullptr, nullptr);
4739  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4740  InVector,
4741  OutVector,
4742  Wrapper,
4743  true>
4744  worker(*this,
4745  src,
4746  dst,
4747  false,
4748  wrap,
4749  &Wrapper::cell_integrator,
4750  &Wrapper::face_integrator,
4751  &Wrapper::boundary_integrator,
4754  operation_before_loop,
4755  operation_after_loop,
4756  dof_handler_index_pre_post);
4757 
4758  task_info.loop(worker);
4759 }
4760 
4761 
4762 
4763 template <int dim, typename Number, typename VectorizedArrayType>
4764 template <typename OutVector, typename InVector>
4765 inline void
4767  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4768  OutVector &,
4769  const InVector &,
4770  const std::pair<unsigned int, unsigned int> &)>
4771  &cell_operation,
4772  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4773  OutVector &,
4774  const InVector &,
4775  const std::pair<unsigned int, unsigned int> &)>
4776  &face_operation,
4777  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4778  OutVector &,
4779  const InVector &,
4780  const std::pair<unsigned int, unsigned int> &)>
4781  & boundary_operation,
4782  OutVector & dst,
4783  const InVector & src,
4784  const bool zero_dst_vector,
4785  const DataAccessOnFaces dst_vector_face_access,
4786  const DataAccessOnFaces src_vector_face_access) const
4787 {
4788  using Wrapper =
4789  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4790  InVector,
4791  OutVector>;
4792  Wrapper wrap(cell_operation, face_operation, boundary_operation);
4793  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4794  InVector,
4795  OutVector,
4796  Wrapper,
4797  true>
4798  worker(*this,
4799  src,
4800  dst,
4801  zero_dst_vector,
4802  wrap,
4803  &Wrapper::cell_integrator,
4804  &Wrapper::face_integrator,
4805  &Wrapper::boundary_integrator,
4806  src_vector_face_access,
4807  dst_vector_face_access);
4808 
4809  task_info.loop(worker);
4810 }
4811 
4812 
4813 
4814 template <int dim, typename Number, typename VectorizedArrayType>
4815 template <typename CLASS, typename OutVector, typename InVector>
4816 inline void
4818  void (CLASS::*function_pointer)(
4820  OutVector &,
4821  const InVector &,
4822  const std::pair<unsigned int, unsigned int> &) const,
4823  const CLASS * owning_class,
4824  OutVector & dst,
4825  const InVector &src,
4826  const bool zero_dst_vector) const
4827 {
4828  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4829  InVector,
4830  OutVector,
4831  CLASS,
4832  true>
4833  worker(*this,
4834  src,
4835  dst,
4836  zero_dst_vector,
4837  *owning_class,
4838  function_pointer,
4839  nullptr,
4840  nullptr);
4841  task_info.loop(worker);
4842 }
4843 
4844 
4845 
4846 template <int dim, typename Number, typename VectorizedArrayType>
4847 template <typename CLASS, typename OutVector, typename InVector>
4848 inline void
4850  void (CLASS::*function_pointer)(
4852  OutVector &,
4853  const InVector &,
4854  const std::pair<unsigned int, unsigned int> &) const,
4855  const CLASS * owning_class,
4856  OutVector & dst,
4857  const InVector &src,
4858  const std::function<void(const unsigned int, const unsigned int)>
4859  &operation_before_loop,
4860  const std::function<void(const unsigned int, const unsigned int)>
4861  & operation_after_loop,
4862  const unsigned int dof_handler_index_pre_post) const
4863 {
4864  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4865  InVector,
4866  OutVector,
4867  CLASS,
4868  true>
4869  worker(*this,
4870  src,
4871  dst,
4872  false,
4873  *owning_class,
4874  function_pointer,
4875  nullptr,
4876  nullptr,
4879  operation_before_loop,
4880  operation_after_loop,
4881  dof_handler_index_pre_post);
4882  task_info.loop(worker);
4883 }
4884 
4885 
4886 
4887 template <int dim, typename Number, typename VectorizedArrayType>
4888 template <typename CLASS, typename OutVector, typename InVector>
4889 inline void
4891  void (CLASS::*cell_operation)(
4893  OutVector &,
4894  const InVector &,
4895  const std::pair<unsigned int, unsigned int> &) const,
4896  void (CLASS::*face_operation)(
4898  OutVector &,
4899  const InVector &,
4900  const std::pair<unsigned int, unsigned int> &) const,
4901  void (CLASS::*boundary_operation)(
4903  OutVector &,
4904  const InVector &,
4905  const std::pair<unsigned int, unsigned int> &) const,
4906  const CLASS * owning_class,
4907  OutVector & dst,
4908  const InVector & src,
4909  const bool zero_dst_vector,
4910  const DataAccessOnFaces dst_vector_face_access,
4911  const DataAccessOnFaces src_vector_face_access) const
4912 {
4913  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4914  InVector,
4915  OutVector,
4916  CLASS,
4917  true>
4918  worker(*this,
4919  src,
4920  dst,
4921  zero_dst_vector,
4922  *owning_class,
4923  cell_operation,
4924  face_operation,
4925  boundary_operation,
4926  src_vector_face_access,
4927  dst_vector_face_access);
4928  task_info.loop(worker);
4929 }
4930 
4931 
4932 
4933 template <int dim, typename Number, typename VectorizedArrayType>
4934 template <typename CLASS, typename OutVector, typename InVector>
4935 inline void
4937  void (CLASS::*function_pointer)(
4939  OutVector &,
4940  const InVector &,
4941  const std::pair<unsigned int, unsigned int> &),
4942  CLASS * owning_class,
4943  OutVector & dst,
4944  const InVector &src,
4945  const bool zero_dst_vector) const
4946 {
4947  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4948  InVector,
4949  OutVector,
4950  CLASS,
4951  false>
4952  worker(*this,
4953  src,
4954  dst,
4955  zero_dst_vector,
4956  *owning_class,
4957  function_pointer,
4958  nullptr,
4959  nullptr);
4960  task_info.loop(worker);
4961 }
4962 
4963 
4964 
4965 template <int dim, typename Number, typename VectorizedArrayType>
4966 template <typename CLASS, typename OutVector, typename InVector>
4967 inline void
4969  void (CLASS::*function_pointer)(
4971  OutVector &,
4972  const InVector &,
4973  const std::pair<unsigned int, unsigned int> &),
4974  CLASS * owning_class,
4975  OutVector & dst,
4976  const InVector &src,
4977  const std::function<void(const unsigned int, const unsigned int)>
4978  &operation_before_loop,
4979  const std::function<void(const unsigned int, const unsigned int)>
4980  & operation_after_loop,
4981  const unsigned int dof_handler_index_pre_post) const
4982 {
4983  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4984  InVector,
4985  OutVector,
4986  CLASS,
4987  false>
4988  worker(*this,
4989  src,
4990  dst,
4991  false,
4992  *owning_class,
4993  function_pointer,
4994  nullptr,
4995  nullptr,
4998  operation_before_loop,
4999  operation_after_loop,
5000  dof_handler_index_pre_post);
5001  task_info.loop(worker);
5002 }
5003 
5004 
5005 
5006 template <int dim, typename Number, typename VectorizedArrayType>
5007 template <typename CLASS, typename OutVector, typename InVector>
5008 inline void
5010  void (CLASS::*cell_operation)(
5012  OutVector &,
5013  const InVector &,
5014  const std::pair<unsigned int, unsigned int> &),
5015  void (CLASS::*face_operation)(
5017  OutVector &,
5018  const InVector &,
5019  const std::pair<unsigned int, unsigned int> &),
5020  void (CLASS::*boundary_operation)(
5022  OutVector &,
5023  const InVector &,
5024  const std::pair<unsigned int, unsigned int> &),
5025  CLASS * owning_class,
5026  OutVector & dst,
5027  const InVector & src,
5028  const bool zero_dst_vector,
5029  const DataAccessOnFaces dst_vector_face_access,
5030  const DataAccessOnFaces src_vector_face_access) const
5031 {
5032  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5033  InVector,
5034  OutVector,
5035  CLASS,
5036  false>
5037  worker(*this,
5038  src,
5039  dst,
5040  zero_dst_vector,
5041  *owning_class,
5042  cell_operation,
5043  face_operation,
5044  boundary_operation,
5045  src_vector_face_access,
5046  dst_vector_face_access);
5047  task_info.loop(worker);
5048 }
5049 
5050 
5051 
5052 template <int dim, typename Number, typename VectorizedArrayType>
5053 template <typename CLASS, typename OutVector, typename InVector>
5054 inline void
5056  void (CLASS::*function_pointer)(
5058  OutVector &,
5059  const InVector &,
5060  const std::pair<unsigned int, unsigned int> &) const,
5061  const CLASS * owning_class,
5062  OutVector & dst,
5063  const InVector & src,
5064  const bool zero_dst_vector,
5065  const DataAccessOnFaces src_vector_face_access) const
5066 {
5067  auto src_vector_face_access_temp = src_vector_face_access;
5068  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5069  src_vector_face_access_temp = DataAccessOnFaces::gradients_cell;
5070  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5071  src_vector_face_access_temp = DataAccessOnFaces::values_cell;
5072 
5073  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5074  InVector,
5075  OutVector,
5076  CLASS,
5077  true>
5078  worker(*this,
5079  src,
5080  dst,
5081  zero_dst_vector,
5082  *owning_class,
5083  function_pointer,
5084  nullptr,
5085  nullptr,
5086  src_vector_face_access_temp,
5088  task_info.loop(worker);
5089 }
5090 
5091 
5092 
5093 template <int dim, typename Number, typename VectorizedArrayType>
5094 template <typename CLASS, typename OutVector, typename InVector>
5095 inline void
5097  void (CLASS::*function_pointer)(
5099  OutVector &,
5100  const InVector &,
5101  const std::pair<unsigned int, unsigned int> &),
5102  CLASS * owning_class,
5103  OutVector & dst,
5104  const InVector & src,
5105  const bool zero_dst_vector,
5106  const DataAccessOnFaces src_vector_face_access) const
5107 {
5108  auto src_vector_face_access_temp = src_vector_face_access;
5109  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5110  src_vector_face_access_temp = DataAccessOnFaces::gradients_cell;
5111  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5112  src_vector_face_access_temp = DataAccessOnFaces::values_cell;
5113 
5114  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5115  InVector,
5116  OutVector,
5117  CLASS,
5118  false>
5119  worker(*this,
5120  src,
5121  dst,
5122  zero_dst_vector,
5123  *owning_class,
5124  function_pointer,
5125  nullptr,
5126  nullptr,
5127  src_vector_face_access_temp,
5129  task_info.loop(worker);
5130 }
5131 
5132 
5133 #endif // ifndef DOXYGEN
5134 
5135 
5136 
5137 DEAL_II_NAMESPACE_CLOSE
5138 
5139 #endif
Transformed quadrature weights.
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > > shape_info
Definition: matrix_free.h:2083
void print_memory_consumption(StreamType &out) const
std::pair< unsigned int, unsigned int > create_cell_subrange_hp(const std::pair< unsigned int, unsigned int > &range, const unsigned int fe_degree, const unsigned int dof_handler_index=0) const
static const unsigned int invalid_unsigned_int
Definition: types.h:187
unsigned int n_components_filled(const unsigned int cell_batch_number) const
std::vector< std::pair< unsigned int, unsigned int > > vector_zero_range_list
Definition: dof_info.h:614
AlignedVector< Number > * acquire_scratch_data_non_threadsafe() const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1571
~MatrixFree() override=default
unsigned int n_boundary_face_batches() const
A class that provides a separate storage location on each thread that accesses the object...
unsigned int local_size() const
void internal_reinit(const Mapping< dim > &mapping, const std::vector< const DoFHandler< dim > *> &dof_handler, const std::vector< const AffineConstraints< number2 > *> &constraint, const std::vector< IndexSet > &locally_owned_set, const std::vector< hp::QCollection< 1 >> &quad, const AdditionalData &additional_data)
void import_from_ghosted_array_finish(const VectorOperation::values vector_operation, const ArrayView< const Number, MemorySpaceType > &temporary_storage, const ArrayView< Number, MemorySpaceType > &locally_owned_storage, const ArrayView< Number, MemorySpaceType > &ghost_array, std::vector< MPI_Request > &requests) 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 get_dofs_per_cell(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
std::vector< unsigned int > cell_loop_pre_list_index
Definition: dof_info.h:621
unsigned int get_cell_category(const unsigned int macro_cell) const
unsigned int get_mg_level() const
void import_from_ghosted_array_start(const VectorOperation::values vector_operation, const unsigned int communication_channel, const ArrayView< Number, MemorySpaceType > &ghost_array, const ArrayView< Number, MemorySpaceType > &temporary_storage, std::vector< MPI_Request > &requests) const
unsigned int n_ghost_cell_batches() const
unsigned int n_components() const
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
unsigned int get_n_q_points_face(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
void loop(MFWorkerInterface &worker) const
Definition: task_info.cc:337
unsigned int n_inner_face_batches() const
DoFHandlers dof_handlers
Definition: matrix_free.h:2050
bool mapping_is_initialized
Definition: matrix_free.h:2124
#define AssertIndexRange(index, range)
Definition: exceptions.h:1641
AdditionalData(const AdditionalData &other)
Definition: matrix_free.h:250
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition: dof_info.h:476
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
std::vector< unsigned int > vector_zero_range_list_index
Definition: dof_info.h:609
std::list< std::pair< bool, AlignedVector< Number > > > scratch_pad_non_threadsafe
Definition: matrix_free.h:2142
void make_connectivity_graph_faces(DynamicSparsityPattern &connectivity)
std::pair< unsigned int, unsigned int > get_face_category(const unsigned int macro_face) const
const DoFHandlerType & get_dof_handler(const unsigned int dof_handler_index=0) const
const std::vector< std::pair< unsigned int, unsigned int > > & ghost_indices_within_larger_ghost_set() const
hp::DoFHandler< dim >::active_cell_iterator get_hp_cell_iterator(const unsigned int macro_cell_number, const unsigned int vector_number, const unsigned int dof_handler_index=0) const
bool mapping_initialized() const
void reinit(const size_type size, const bool omit_zeroing_entries=false)
Number value_type
Definition: matrix_free.h:128
UpdateFlags mapping_update_flags_boundary_faces
Definition: matrix_free.h:381
static ::ExceptionBase & ExcNotInitialized()
std::vector< unsigned int > cell_partition_data
Definition: task_info.h:473
#define AssertThrow(cond, exc)
Definition: exceptions.h:1523
unsigned int get_n_q_points(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
UpdateFlags mapping_update_flags
Definition: matrix_free.h:360
std::size_t memory_consumption() const
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
const Number * constraint_pool_end(const unsigned int pool_index) const
void clear()
const IndexSet & get_ghost_set(const unsigned int dof_handler_index=0) const
::Table< 3, unsigned int > cell_and_face_to_plain_faces
Definition: face_info.h:163
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_number) const
static const unsigned int dimension
Definition: matrix_free.h:134
std::vector< unsigned int > cell_loop_post_list_index
Definition: dof_info.h:634
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:312
std::pair< unsigned int, unsigned int > create_cell_subrange_hp_by_index(const std::pair< unsigned int, unsigned int > &range, const unsigned int fe_index, const unsigned int dof_handler_index=0) const
void print(std::ostream &out) const
const Number * constraint_pool_begin(const unsigned int pool_index) const
DoFHandler< dim >::cell_iterator get_cell_iterator(const unsigned int macro_cell_number, const unsigned int vector_number, const unsigned int dof_handler_index=0) const
void set_ghost_state(const bool ghosted) const
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
static ::ExceptionBase & ExcMessage(std::string arg1)
bool indices_initialized() const
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
const IndexSet & get_locally_owned_set(const unsigned int dof_handler_index=0) const
void initialize_dof_handlers(const std::vector< const DoFHandler< dim > *> &dof_handlers, const AdditionalData &additional_data)
No update.
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
unsigned int n_base_elements(const unsigned int dof_handler_index) const
unsigned int n_ghost_indices() const
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_pre_list
Definition: dof_info.h:627
unsigned int n_ghost_inner_face_batches() const
#define Assert(cond, exc)
Definition: exceptions.h:1411
UpdateFlags
const types::boundary_id invalid_boundary_id
Definition: types.h:228
unsigned int mg_level
Definition: matrix_free.h:2147
std::vector< std::pair< unsigned int, unsigned int > > cell_level_index
Definition: matrix_free.h:2091
UpdateFlags mapping_update_flags_inner_faces
Definition: matrix_free.h:402
std::vector< unsigned int > boundary_partition_data
Definition: task_info.h:491
std::vector< internal::MatrixFreeFunctions::DoFInfo > dof_info
Definition: matrix_free.h:2056
unsigned int get_dofs_per_face(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::n_array_elements > & get_face_info(const unsigned int face_batch_number) const
UpdateFlags mapping_update_flags_faces_by_cells
Definition: matrix_free.h:430
unsigned int & level_mg_handler
Definition: matrix_free.h:447
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_number) const
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
unsigned int n_constraint_pool_entries() const
TasksParallelScheme tasks_parallel_scheme
Definition: matrix_free.h:336
std::vector< unsigned int > cell_vectorization_category
Definition: matrix_free.h:516
unsigned int n_cell_batches() const
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner(const unsigned int dof_handler_index=0) const
Definition: hp.h:117
std::vector< FaceToCellTopology< vectorization_width > > faces
Definition: face_info.h:156
std::vector< Number > constraint_pool_data
Definition: matrix_free.h:2064
types::boundary_id get_boundary_id(const unsigned int macro_face) const
std::vector< unsigned int > face_partition_data
Definition: task_info.h:482
unsigned int n_macro_cells() const
::Table< 3, types::boundary_id > cell_and_face_boundary_id
Definition: face_info.h:169
const std::vector< unsigned int > & get_constrained_dofs(const unsigned int dof_handler_index=0) const
internal::MatrixFreeFunctions::TaskInfo task_info
Definition: matrix_free.h:2107
AdditionalData(const TasksParallelScheme tasks_parallel_scheme=partition_partition, const unsigned int tasks_block_size=0, const UpdateFlags mapping_update_flags=update_gradients|update_JxW_values, const UpdateFlags mapping_update_flags_boundary_faces=update_default, const UpdateFlags mapping_update_flags_inner_faces=update_default, const UpdateFlags mapping_update_flags_faces_by_cells=update_default, const unsigned int mg_level=numbers::invalid_unsigned_int, const bool store_plain_indices=true, const bool initialize_indices=true, const bool initialize_mapping=true, const bool overlap_communication_computation=true, const bool hold_all_faces_to_owned_cells=false, const bool cell_vectorization_categories_strict=false)
Definition: matrix_free.h:215
size_type size(const unsigned int i) const
std::array< types::boundary_id, VectorizedArrayType::n_array_elements > get_faces_by_cells_boundary_id(const unsigned int macro_cell, const unsigned int face_number) const
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
void initialize_dof_vector(VectorType &vec, const unsigned int dof_handler_index=0) const
Shape function gradients.
bool indices_are_initialized
Definition: matrix_free.h:2119
std::pair< int, int > get_cell_level_and_index(const unsigned int macro_cell_number, const unsigned int vector_number) const
void renumber_dofs(std::vector< types::global_dof_index > &renumbering, const unsigned int dof_handler_index=0)
std::pair< typename DoFHandler< dim >::cell_iterator, unsigned int > get_face_iterator(const unsigned int face_batch_number, const unsigned int vector_number, const bool interior=true, const unsigned int fe_component=0) const
Threads::ThreadLocalStorage< std::list< std::pair< bool, AlignedVector< VectorizedArrayType > > > > scratch_pad
Definition: matrix_free.h:2135
void reinit(const Mapping< dim > &mapping, const DoFHandlerType &dof_handler, const AffineConstraints< number2 > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
void export_to_ghosted_array_start(const unsigned int communication_channel, const ArrayView< const Number, MemorySpaceType > &locally_owned_array, const ArrayView< Number, MemorySpaceType > &temporary_storage, const ArrayView< Number, MemorySpaceType > &ghost_array, std::vector< MPI_Request > &requests) const
void export_to_ghosted_array_finish(const ArrayView< Number, MemorySpaceType > &ghost_array, std::vector< MPI_Request > &requests) const
internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > mapping_info
Definition: matrix_free.h:2077
AdditionalData & operator=(const AdditionalData &other)
Definition: matrix_free.h:276
static ::ExceptionBase & ExcNotImplemented()
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
void release_scratch_data_non_threadsafe(const AlignedVector< Number > *memory) const
bool at_irregular_cell(const unsigned int macro_cell_number) const
std::vector< unsigned int > constraint_pool_row_index
Definition: matrix_free.h:2070
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_post_list
Definition: dof_info.h:640
static bool is_supported(const FiniteElement< dim, spacedim > &fe)
void loop_cell_centric(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
internal::MatrixFreeFunctions::FaceInfo< VectorizedArrayType::n_array_elements > face_info
Definition: matrix_free.h:2114
unsigned int tasks_block_size
Definition: matrix_free.h:347
const Quadrature< dim > & get_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int n_import_indices() const
unsigned int boundary_id
Definition: types.h:125
unsigned int cell_level_index_end_local
Definition: matrix_free.h:2100
unsigned int n_physical_cells() const
void copy_from(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free_base)
static bool equal(const T *p1, const T *p2)
static ::ExceptionBase & ExcInternalError()