Reference documentation for deal.II version Git cb0bd54b52 2019-09-21 16:31:22 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
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/aligned_vector.h>
21 #include <deal.II/base/exceptions.h>
22 #include <deal.II/base/quadrature.h>
23 #include <deal.II/base/template_constraints.h>
24 #include <deal.II/base/thread_local_storage.h>
25 #include <deal.II/base/vectorization.h>
26 
27 #include <deal.II/dofs/dof_handler.h>
28 
29 #include <deal.II/fe/fe.h>
30 #include <deal.II/fe/mapping.h>
31 #include <deal.II/fe/mapping_q1.h>
32 
33 #include <deal.II/grid/grid_tools.h>
34 
35 #include <deal.II/hp/dof_handler.h>
36 #include <deal.II/hp/q_collection.h>
37 
38 #include <deal.II/lac/affine_constraints.h>
39 #include <deal.II/lac/block_vector_base.h>
40 #include <deal.II/lac/la_parallel_vector.h>
41 #include <deal.II/lac/vector_operation.h>
42 
43 #include <deal.II/matrix_free/dof_info.h>
44 #include <deal.II/matrix_free/mapping_info.h>
45 #include <deal.II/matrix_free/shape_info.h>
46 #include <deal.II/matrix_free/task_info.h>
47 #include <deal.II/matrix_free/type_traits.h>
48 
49 #include <cstdlib>
50 #include <limits>
51 #include <list>
52 #include <memory>
53 
54 
55 DEAL_II_NAMESPACE_OPEN
56 
57 
58 
112 template <int dim,
113  typename Number = double,
114  typename VectorizedArrayType = VectorizedArray<Number>>
115 class MatrixFree : public Subscriptor
116 {
117  static_assert(
118  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
119  "Type of Number and of VectorizedArrayType do not match.");
120 
121 public:
126  using value_type = Number;
127  using vectorized_value_type = VectorizedArrayType;
128 
132  static const unsigned int dimension = dim;
133 
181  {
188  {
192  none = internal::MatrixFreeFunctions::TaskInfo::none,
197  internal::MatrixFreeFunctions::TaskInfo::partition_partition,
202  internal::MatrixFreeFunctions::TaskInfo::partition_color,
207  color = internal::MatrixFreeFunctions::TaskInfo::color
208  };
209 
215  const unsigned int tasks_block_size = 0,
221  const unsigned int level_mg_handler = numbers::invalid_unsigned_int,
222  const bool store_plain_indices = true,
223  const bool initialize_indices = true,
224  const bool initialize_mapping = true,
225  const bool overlap_communication_computation = true,
226  const bool hold_all_faces_to_owned_cells = false,
227  const bool cell_vectorization_categories_strict = false)
242  {}
243 
281 
291  unsigned int tasks_block_size;
292 
305 
326 
347 
375 
384  unsigned int level_mg_handler;
385 
393 
404 
413 
427 
436 
453  std::vector<unsigned int> cell_vectorization_category;
454 
465  };
466 
474  MatrixFree();
475 
480 
484  ~MatrixFree() override = default;
485 
498  template <typename DoFHandlerType, typename QuadratureType, typename number2>
499  void
500  reinit(const Mapping<dim> & mapping,
501  const DoFHandlerType & dof_handler,
502  const AffineConstraints<number2> &constraint,
503  const QuadratureType & quad,
504  const AdditionalData & additional_data = AdditionalData());
505 
510  template <typename DoFHandlerType, typename QuadratureType, typename number2>
511  void
512  reinit(const DoFHandlerType & dof_handler,
513  const AffineConstraints<number2> &constraint,
514  const QuadratureType & quad,
515  const AdditionalData & additional_data = AdditionalData());
516 
524  template <typename DoFHandlerType, typename QuadratureType, typename number2>
525  DEAL_II_DEPRECATED void
526  reinit(const Mapping<dim> & mapping,
527  const DoFHandlerType & dof_handler,
528  const AffineConstraints<number2> &constraint,
529  const IndexSet & locally_owned_dofs,
530  const QuadratureType & quad,
531  const AdditionalData & additional_data = AdditionalData());
532 
554  template <typename DoFHandlerType, typename QuadratureType, typename number2>
555  void
556  reinit(const Mapping<dim> & mapping,
557  const std::vector<const DoFHandlerType *> & dof_handler,
558  const std::vector<const AffineConstraints<number2> *> &constraint,
559  const std::vector<QuadratureType> & quad,
560  const AdditionalData &additional_data = AdditionalData());
561 
566  template <typename DoFHandlerType, typename QuadratureType, typename number2>
567  void
568  reinit(const std::vector<const DoFHandlerType *> & dof_handler,
569  const std::vector<const AffineConstraints<number2> *> &constraint,
570  const std::vector<QuadratureType> & quad,
571  const AdditionalData &additional_data = AdditionalData());
572 
580  template <typename DoFHandlerType, typename QuadratureType, typename number2>
581  DEAL_II_DEPRECATED void
582  reinit(const Mapping<dim> & mapping,
583  const std::vector<const DoFHandlerType *> & dof_handler,
584  const std::vector<const AffineConstraints<number2> *> &constraint,
585  const std::vector<IndexSet> & locally_owned_set,
586  const std::vector<QuadratureType> &quad,
587  const AdditionalData & additional_data = AdditionalData());
588 
596  template <typename DoFHandlerType, typename QuadratureType, typename number2>
597  void
598  reinit(const Mapping<dim> & mapping,
599  const std::vector<const DoFHandlerType *> & dof_handler,
600  const std::vector<const AffineConstraints<number2> *> &constraint,
601  const QuadratureType & quad,
602  const AdditionalData &additional_data = AdditionalData());
603 
608  template <typename DoFHandlerType, typename QuadratureType, typename number2>
609  void
610  reinit(const std::vector<const DoFHandlerType *> & dof_handler,
611  const std::vector<const AffineConstraints<number2> *> &constraint,
612  const QuadratureType & quad,
613  const AdditionalData &additional_data = AdditionalData());
614 
620  void
621  copy_from(
622  const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free_base);
623 
628  void
629  clear();
630 
632 
644  enum class DataAccessOnFaces
645  {
652  none,
653 
663  values,
664 
675  gradients,
676 
683  unspecified
684  };
685 
726  template <typename OutVector, typename InVector>
727  void
728  cell_loop(const std::function<void(
730  OutVector &,
731  const InVector &,
732  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
733  OutVector & dst,
734  const InVector & src,
735  const bool zero_dst_vector = false) const;
736 
778  template <typename CLASS, typename OutVector, typename InVector>
779  void
780  cell_loop(void (CLASS::*cell_operation)(
781  const MatrixFree &,
782  OutVector &,
783  const InVector &,
784  const std::pair<unsigned int, unsigned int> &) const,
785  const CLASS * owning_class,
786  OutVector & dst,
787  const InVector &src,
788  const bool zero_dst_vector = false) const;
789 
793  template <typename CLASS, typename OutVector, typename InVector>
794  void
795  cell_loop(void (CLASS::*cell_operation)(
796  const MatrixFree &,
797  OutVector &,
798  const InVector &,
799  const std::pair<unsigned int, unsigned int> &),
800  CLASS * owning_class,
801  OutVector & dst,
802  const InVector &src,
803  const bool zero_dst_vector = false) const;
804 
880  template <typename OutVector, typename InVector>
881  void
882  loop(const std::function<
884  OutVector &,
885  const InVector &,
886  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
887  const std::function<
889  OutVector &,
890  const InVector &,
891  const std::pair<unsigned int, unsigned int> &)> &face_operation,
892  const std::function<void(
894  OutVector &,
895  const InVector &,
896  const std::pair<unsigned int, unsigned int> &)> &boundary_operation,
897  OutVector & dst,
898  const InVector & src,
899  const bool zero_dst_vector = false,
900  const DataAccessOnFaces dst_vector_face_access =
902  const DataAccessOnFaces src_vector_face_access =
904 
989  template <typename CLASS, typename OutVector, typename InVector>
990  void
991  loop(
992  void (CLASS::*cell_operation)(const MatrixFree &,
993  OutVector &,
994  const InVector &,
995  const std::pair<unsigned int, unsigned int> &)
996  const,
997  void (CLASS::*face_operation)(const MatrixFree &,
998  OutVector &,
999  const InVector &,
1000  const std::pair<unsigned int, unsigned int> &)
1001  const,
1002  void (CLASS::*boundary_operation)(
1003  const MatrixFree &,
1004  OutVector &,
1005  const InVector &,
1006  const std::pair<unsigned int, unsigned int> &) const,
1007  const CLASS * owning_class,
1008  OutVector & dst,
1009  const InVector & src,
1010  const bool zero_dst_vector = false,
1011  const DataAccessOnFaces dst_vector_face_access =
1013  const DataAccessOnFaces src_vector_face_access =
1015 
1019  template <typename CLASS, typename OutVector, typename InVector>
1020  void
1021  loop(void (CLASS::*cell_operation)(
1022  const MatrixFree &,
1023  OutVector &,
1024  const InVector &,
1025  const std::pair<unsigned int, unsigned int> &),
1026  void (CLASS::*face_operation)(
1027  const MatrixFree &,
1028  OutVector &,
1029  const InVector &,
1030  const std::pair<unsigned int, unsigned int> &),
1031  void (CLASS::*boundary_operation)(
1032  const MatrixFree &,
1033  OutVector &,
1034  const InVector &,
1035  const std::pair<unsigned int, unsigned int> &),
1036  CLASS * owning_class,
1037  OutVector & dst,
1038  const InVector & src,
1039  const bool zero_dst_vector = false,
1040  const DataAccessOnFaces dst_vector_face_access =
1042  const DataAccessOnFaces src_vector_face_access =
1044 
1052  std::pair<unsigned int, unsigned int>
1053  create_cell_subrange_hp(const std::pair<unsigned int, unsigned int> &range,
1054  const unsigned int fe_degree,
1055  const unsigned int dof_handler_index = 0) const;
1056 
1063  std::pair<unsigned int, unsigned int>
1065  const std::pair<unsigned int, unsigned int> &range,
1066  const unsigned int fe_index,
1067  const unsigned int dof_handler_index = 0) const;
1068 
1070 
1095  template <typename VectorType>
1096  void
1097  initialize_dof_vector(VectorType & vec,
1098  const unsigned int dof_handler_index = 0) const;
1099 
1120  template <typename Number2>
1121  void
1123  const unsigned int dof_handler_index = 0) const;
1124 
1135  const std::shared_ptr<const Utilities::MPI::Partitioner> &
1136  get_vector_partitioner(const unsigned int dof_handler_index = 0) const;
1137 
1141  const IndexSet &
1142  get_locally_owned_set(const unsigned int dof_handler_index = 0) const;
1143 
1147  const IndexSet &
1148  get_ghost_set(const unsigned int dof_handler_index = 0) const;
1149 
1159  const std::vector<unsigned int> &
1160  get_constrained_dofs(const unsigned int dof_handler_index = 0) const;
1161 
1172  void
1173  renumber_dofs(std::vector<types::global_dof_index> &renumbering,
1174  const unsigned int dof_handler_index = 0);
1175 
1177 
1185  template <int spacedim>
1186  static bool
1188 
1192  unsigned int
1193  n_components() const;
1194 
1199  unsigned int
1200  n_base_elements(const unsigned int dof_handler_index) const;
1201 
1209  unsigned int
1210  n_physical_cells() const;
1211 
1221  unsigned int
1222  n_macro_cells() const;
1223 
1233  unsigned int
1234  n_cell_batches() const;
1235 
1242  unsigned int
1243  n_ghost_cell_batches() const;
1244 
1252  unsigned int
1253  n_inner_face_batches() const;
1254 
1263  unsigned int
1264  n_boundary_face_batches() const;
1265 
1270  unsigned int
1272 
1280  get_boundary_id(const unsigned int macro_face) const;
1281 
1286  std::array<types::boundary_id, VectorizedArrayType::n_array_elements>
1287  get_faces_by_cells_boundary_id(const unsigned int macro_cell,
1288  const unsigned int face_number) const;
1289 
1294  const DoFHandler<dim> &
1295  get_dof_handler(const unsigned int dof_handler_index = 0) const;
1296 
1308  get_cell_iterator(const unsigned int macro_cell_number,
1309  const unsigned int vector_number,
1310  const unsigned int dof_handler_index = 0) const;
1311 
1324  get_hp_cell_iterator(const unsigned int macro_cell_number,
1325  const unsigned int vector_number,
1326  const unsigned int dof_handler_index = 0) const;
1327 
1339  bool
1340  at_irregular_cell(const unsigned int macro_cell_number) const;
1341 
1350  unsigned int
1351  n_components_filled(const unsigned int cell_batch_number) const;
1352 
1361  unsigned int
1362  n_active_entries_per_cell_batch(const unsigned int cell_batch_number) const;
1363 
1373  unsigned int
1374  n_active_entries_per_face_batch(const unsigned int face_batch_number) const;
1375 
1379  unsigned int
1380  get_dofs_per_cell(const unsigned int dof_handler_index = 0,
1381  const unsigned int hp_active_fe_index = 0) const;
1382 
1386  unsigned int
1387  get_n_q_points(const unsigned int quad_index = 0,
1388  const unsigned int hp_active_fe_index = 0) const;
1389 
1394  unsigned int
1395  get_dofs_per_face(const unsigned int dof_handler_index = 0,
1396  const unsigned int hp_active_fe_index = 0) const;
1397 
1402  unsigned int
1403  get_n_q_points_face(const unsigned int quad_index = 0,
1404  const unsigned int hp_active_fe_index = 0) const;
1405 
1409  const Quadrature<dim> &
1410  get_quadrature(const unsigned int quad_index = 0,
1411  const unsigned int hp_active_fe_index = 0) const;
1412 
1416  const Quadrature<dim - 1> &
1417  get_face_quadrature(const unsigned int quad_index = 0,
1418  const unsigned int hp_active_fe_index = 0) const;
1419 
1426  unsigned int
1427  get_cell_category(const unsigned int macro_cell) const;
1428 
1433  std::pair<unsigned int, unsigned int>
1434  get_face_category(const unsigned int macro_face) const;
1435 
1439  bool
1440  indices_initialized() const;
1441 
1446  bool
1447  mapping_initialized() const;
1448 
1453  std::size_t
1454  memory_consumption() const;
1455 
1460  template <typename StreamType>
1461  void
1462  print_memory_consumption(StreamType &out) const;
1463 
1468  void
1469  print(std::ostream &out) const;
1470 
1472 
1483  get_task_info() const;
1484 
1488  DEAL_II_DEPRECATED
1490  get_size_info() const;
1491 
1492  /*
1493  * Return geometry-dependent information on the cells.
1494  */
1495  const internal::MatrixFreeFunctions::
1496  MappingInfo<dim, Number, VectorizedArrayType> &
1497  get_mapping_info() const;
1498 
1503  get_dof_info(const unsigned int dof_handler_index_component = 0) const;
1504 
1508  unsigned int
1509  n_constraint_pool_entries() const;
1510 
1515  const Number *
1516  constraint_pool_begin(const unsigned int pool_index) const;
1517 
1523  const Number *
1524  constraint_pool_end(const unsigned int pool_index) const;
1525 
1530  get_shape_info(const unsigned int dof_handler_index_component = 0,
1531  const unsigned int quad_index = 0,
1532  const unsigned int fe_base_element = 0,
1533  const unsigned int hp_active_fe_index = 0,
1534  const unsigned int hp_active_quad_index = 0) const;
1535 
1540  VectorizedArrayType::n_array_elements> &
1541  get_face_info(const unsigned int face_batch_number) const;
1542 
1557  acquire_scratch_data() const;
1558 
1562  void
1564 
1576 
1580  void
1582  const AlignedVector<Number> *memory) const;
1583 
1585 
1586 private:
1591  template <typename number2>
1592  void
1594  const Mapping<dim> & mapping,
1595  const std::vector<const DoFHandler<dim> *> & dof_handler,
1596  const std::vector<const AffineConstraints<number2> *> &constraint,
1597  const std::vector<IndexSet> & locally_owned_set,
1598  const std::vector<hp::QCollection<1>> & quad,
1599  const AdditionalData & additional_data);
1600 
1604  template <typename number2>
1605  void
1607  const Mapping<dim> & mapping,
1608  const std::vector<const hp::DoFHandler<dim> *> & dof_handler,
1609  const std::vector<const AffineConstraints<number2> *> &constraint,
1610  const std::vector<IndexSet> & locally_owned_set,
1611  const std::vector<hp::QCollection<1>> & quad,
1612  const AdditionalData & additional_data);
1613 
1620  template <typename number2>
1621  void
1623  const std::vector<const AffineConstraints<number2> *> &constraint,
1624  const std::vector<IndexSet> & locally_owned_set,
1625  const AdditionalData & additional_data);
1626 
1630  void
1632  const std::vector<const DoFHandler<dim> *> &dof_handlers,
1633  const AdditionalData & additional_data);
1634 
1638  void
1640  const std::vector<const hp::DoFHandler<dim> *> &dof_handlers,
1641  const AdditionalData & additional_data);
1642 
1647  void
1649 
1656  {
1657  DoFHandlers()
1658  : active_dof_handler(usual)
1659  , n_dof_handlers(0)
1661  {}
1662 
1663  std::vector<SmartPointer<const DoFHandler<dim>>> dof_handler;
1664  std::vector<SmartPointer<const hp::DoFHandler<dim>>> hp_dof_handler;
1666  {
1675  } active_dof_handler;
1676  unsigned int n_dof_handlers;
1677  unsigned int level;
1678  };
1679 
1684 
1689  std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
1690 
1697  std::vector<Number> constraint_pool_data;
1698 
1703  std::vector<unsigned int> constraint_pool_row_index;
1704 
1711 
1717 
1724  std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
1725 
1726 
1734 
1741 
1748 
1753 
1758 
1767  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
1769 
1774  mutable std::list<std::pair<bool, AlignedVector<Number>>>
1776 };
1777 
1778 
1779 
1780 /*----------------------- Inline functions ----------------------------------*/
1781 
1782 #ifndef DOXYGEN
1783 
1784 
1785 
1786 template <int dim, typename Number, typename VectorizedArrayType>
1787 template <typename VectorType>
1788 inline void
1790  VectorType & vec,
1791  const unsigned int comp) const
1792 {
1793  AssertIndexRange(comp, n_components());
1794  vec.reinit(dof_info[comp].vector_partitioner->size());
1795 }
1796 
1797 
1798 
1799 template <int dim, typename Number, typename VectorizedArrayType>
1800 template <typename Number2>
1801 inline void
1804  const unsigned int comp) const
1805 {
1806  AssertIndexRange(comp, n_components());
1807  vec.reinit(dof_info[comp].vector_partitioner);
1808 }
1809 
1810 
1811 
1812 template <int dim, typename Number, typename VectorizedArrayType>
1813 inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
1815  const unsigned int comp) const
1816 {
1817  AssertIndexRange(comp, n_components());
1818  return dof_info[comp].vector_partitioner;
1819 }
1820 
1821 
1822 
1823 template <int dim, typename Number, typename VectorizedArrayType>
1824 inline const std::vector<unsigned int> &
1826  const unsigned int comp) const
1827 {
1828  AssertIndexRange(comp, n_components());
1829  return dof_info[comp].constrained_dofs;
1830 }
1831 
1832 
1833 
1834 template <int dim, typename Number, typename VectorizedArrayType>
1835 inline unsigned int
1837 {
1838  AssertDimension(dof_handlers.n_dof_handlers, dof_info.size());
1839  return dof_handlers.n_dof_handlers;
1840 }
1841 
1842 
1843 
1844 template <int dim, typename Number, typename VectorizedArrayType>
1845 inline unsigned int
1847  const unsigned int dof_no) const
1848 {
1849  AssertDimension(dof_handlers.n_dof_handlers, dof_info.size());
1850  AssertIndexRange(dof_no, dof_handlers.n_dof_handlers);
1851  return dof_handlers.dof_handler[dof_no]->get_fe().n_base_elements();
1852 }
1853 
1854 
1855 
1856 template <int dim, typename Number, typename VectorizedArrayType>
1859 {
1860  return task_info;
1861 }
1862 
1863 
1864 
1865 template <int dim, typename Number, typename VectorizedArrayType>
1868 {
1869  return task_info;
1870 }
1871 
1872 
1873 
1874 template <int dim, typename Number, typename VectorizedArrayType>
1875 inline unsigned int
1877 {
1878  return *(task_info.cell_partition_data.end() - 2);
1879 }
1880 
1881 
1882 
1883 template <int dim, typename Number, typename VectorizedArrayType>
1884 inline unsigned int
1886 {
1887  return task_info.n_active_cells;
1888 }
1889 
1890 
1891 
1892 template <int dim, typename Number, typename VectorizedArrayType>
1893 inline unsigned int
1895 {
1896  return *(task_info.cell_partition_data.end() - 2);
1897 }
1898 
1899 
1900 
1901 template <int dim, typename Number, typename VectorizedArrayType>
1902 inline unsigned int
1904 {
1905  return *(task_info.cell_partition_data.end() - 1) -
1906  *(task_info.cell_partition_data.end() - 2);
1907 }
1908 
1909 
1910 
1911 template <int dim, typename Number, typename VectorizedArrayType>
1912 inline unsigned int
1914 {
1915  if (task_info.face_partition_data.size() == 0)
1916  return 0;
1917  return task_info.face_partition_data.back();
1918 }
1919 
1920 
1921 
1922 template <int dim, typename Number, typename VectorizedArrayType>
1923 inline unsigned int
1925 {
1926  if (task_info.face_partition_data.size() == 0)
1927  return 0;
1928  return task_info.boundary_partition_data.back() -
1930 }
1931 
1932 
1933 
1934 template <int dim, typename Number, typename VectorizedArrayType>
1935 inline unsigned int
1937 {
1938  if (task_info.face_partition_data.size() == 0)
1939  return 0;
1940  return face_info.faces.size() - task_info.boundary_partition_data.back();
1941 }
1942 
1943 
1944 
1945 template <int dim, typename Number, typename VectorizedArrayType>
1946 inline types::boundary_id
1948  const unsigned int macro_face) const
1949 {
1950  Assert(macro_face >= task_info.boundary_partition_data[0] &&
1951  macro_face < task_info.boundary_partition_data.back(),
1952  ExcIndexRange(macro_face,
1955  return types::boundary_id(face_info.faces[macro_face].exterior_face_no);
1956 }
1957 
1958 
1959 
1960 template <int dim, typename Number, typename VectorizedArrayType>
1961 inline std::array<types::boundary_id, VectorizedArrayType::n_array_elements>
1963  const unsigned int macro_cell,
1964  const unsigned int face_number) const
1965 {
1966  AssertIndexRange(macro_cell, n_macro_cells());
1969  ExcNotInitialized());
1970  std::array<types::boundary_id, VectorizedArrayType::n_array_elements> result;
1971  result.fill(numbers::invalid_boundary_id);
1972  for (unsigned int v = 0; v < n_active_entries_per_cell_batch(macro_cell); ++v)
1973  result[v] = face_info.cell_and_face_boundary_id(macro_cell, face_number, v);
1974  return result;
1975 }
1976 
1977 
1978 
1979 template <int dim, typename Number, typename VectorizedArrayType>
1980 inline const internal::MatrixFreeFunctions::
1981  MappingInfo<dim, Number, VectorizedArrayType> &
1983 {
1984  return mapping_info;
1985 }
1986 
1987 
1988 
1989 template <int dim, typename Number, typename VectorizedArrayType>
1992  const unsigned int dof_index) const
1993 {
1994  AssertIndexRange(dof_index, n_components());
1995  return dof_info[dof_index];
1996 }
1997 
1998 
1999 
2000 template <int dim, typename Number, typename VectorizedArrayType>
2001 inline unsigned int
2003 {
2004  return constraint_pool_row_index.size() - 1;
2005 }
2006 
2007 
2008 
2009 template <int dim, typename Number, typename VectorizedArrayType>
2010 inline const Number *
2012  const unsigned int row) const
2013 {
2015  return constraint_pool_data.empty() ?
2016  nullptr :
2018 }
2019 
2020 
2021 
2022 template <int dim, typename Number, typename VectorizedArrayType>
2023 inline const Number *
2025  const unsigned int row) const
2026 {
2028  return constraint_pool_data.empty() ?
2029  nullptr :
2031 }
2032 
2033 
2034 
2035 template <int dim, typename Number, typename VectorizedArrayType>
2036 inline std::pair<unsigned int, unsigned int>
2038  const std::pair<unsigned int, unsigned int> &range,
2039  const unsigned int degree,
2040  const unsigned int dof_handler_component) const
2041 {
2042  if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2043  {
2045  dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2047  dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2048  if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2049  return range;
2050  else
2051  return {range.second, range.second};
2052  }
2053 
2054  const unsigned int fe_index =
2055  dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2056  if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2057  return {range.second, range.second};
2058  else
2059  return create_cell_subrange_hp_by_index(range,
2060  fe_index,
2061  dof_handler_component);
2062 }
2063 
2064 
2065 
2066 template <int dim, typename Number, typename VectorizedArrayType>
2067 inline bool
2069  const unsigned int macro_cell) const
2070 {
2071  AssertIndexRange(macro_cell, task_info.cell_partition_data.back());
2072  return VectorizedArrayType::n_array_elements > 1 &&
2073  cell_level_index[(macro_cell + 1) *
2074  VectorizedArrayType::n_array_elements -
2075  1] ==
2076  cell_level_index[(macro_cell + 1) *
2077  VectorizedArrayType::n_array_elements -
2078  2];
2079 }
2080 
2081 
2082 
2083 template <int dim, typename Number, typename VectorizedArrayType>
2084 inline unsigned int
2086  const unsigned int cell_batch_number) const
2087 {
2088  return n_active_entries_per_cell_batch(cell_batch_number);
2089 }
2090 
2091 
2092 
2093 template <int dim, typename Number, typename VectorizedArrayType>
2094 inline unsigned int
2096  const unsigned int cell_batch_number) const
2097 {
2098  AssertIndexRange(cell_batch_number, task_info.cell_partition_data.back());
2099  unsigned int n_components = VectorizedArrayType::n_array_elements;
2100  while (
2101  n_components > 1 &&
2102  cell_level_index[cell_batch_number * VectorizedArrayType::n_array_elements +
2103  n_components - 1] ==
2104  cell_level_index[cell_batch_number *
2105  VectorizedArrayType::n_array_elements +
2106  n_components - 2])
2107  --n_components;
2108  AssertIndexRange(n_components - 1, VectorizedArrayType::n_array_elements);
2109  return n_components;
2110 }
2111 
2112 
2113 
2114 template <int dim, typename Number, typename VectorizedArrayType>
2115 inline unsigned int
2117  const unsigned int face_batch_number) const
2118 {
2119  AssertIndexRange(face_batch_number, face_info.faces.size());
2120  unsigned int n_components = VectorizedArrayType::n_array_elements;
2121  while (n_components > 1 &&
2122  face_info.faces[face_batch_number].cells_interior[n_components - 1] ==
2124  --n_components;
2125  AssertIndexRange(n_components - 1, VectorizedArrayType::n_array_elements);
2126  return n_components;
2127 }
2128 
2129 
2130 
2131 template <int dim, typename Number, typename VectorizedArrayType>
2132 inline unsigned int
2134  const unsigned int dof_handler_index,
2135  const unsigned int active_fe_index) const
2136 {
2137  return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2138 }
2139 
2140 
2141 
2142 template <int dim, typename Number, typename VectorizedArrayType>
2143 inline unsigned int
2145  const unsigned int quad_index,
2146  const unsigned int active_fe_index) const
2147 {
2148  AssertIndexRange(quad_index, mapping_info.cell_data.size());
2149  return mapping_info.cell_data[quad_index]
2150  .descriptor[active_fe_index]
2151  .n_q_points;
2152 }
2153 
2154 
2155 
2156 template <int dim, typename Number, typename VectorizedArrayType>
2157 inline unsigned int
2159  const unsigned int dof_handler_index,
2160  const unsigned int active_fe_index) const
2161 {
2162  return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2163 }
2164 
2165 
2166 
2167 template <int dim, typename Number, typename VectorizedArrayType>
2168 inline unsigned int
2170  const unsigned int quad_index,
2171  const unsigned int active_fe_index) const
2172 {
2173  AssertIndexRange(quad_index, mapping_info.face_data.size());
2174  return mapping_info.face_data[quad_index]
2175  .descriptor[active_fe_index]
2176  .n_q_points;
2177 }
2178 
2179 
2180 
2181 template <int dim, typename Number, typename VectorizedArrayType>
2182 inline const IndexSet &
2184  const unsigned int dof_handler_index) const
2185 {
2186  return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2187 }
2188 
2189 
2190 
2191 template <int dim, typename Number, typename VectorizedArrayType>
2192 inline const IndexSet &
2194  const unsigned int dof_handler_index) const
2195 {
2196  return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2197 }
2198 
2199 
2200 
2201 template <int dim, typename Number, typename VectorizedArrayType>
2204  const unsigned int dof_handler_index,
2205  const unsigned int index_quad,
2206  const unsigned int index_fe,
2207  const unsigned int active_fe_index,
2208  const unsigned int active_quad_index) const
2209 {
2210  AssertIndexRange(dof_handler_index, dof_info.size());
2211  const unsigned int ind =
2212  dof_info[dof_handler_index].global_base_element_offset + index_fe;
2213  AssertIndexRange(ind, shape_info.size(0));
2214  AssertIndexRange(index_quad, shape_info.size(1));
2215  AssertIndexRange(active_fe_index, shape_info.size(2));
2216  AssertIndexRange(active_quad_index, shape_info.size(3));
2217  return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2218 }
2219 
2220 
2221 
2222 template <int dim, typename Number, typename VectorizedArrayType>
2224  VectorizedArrayType::n_array_elements> &
2226  const unsigned int macro_face) const
2227 {
2228  AssertIndexRange(macro_face, face_info.faces.size());
2229  return face_info.faces[macro_face];
2230 }
2231 
2232 
2233 
2234 template <int dim, typename Number, typename VectorizedArrayType>
2235 inline const Quadrature<dim> &
2237  const unsigned int quad_index,
2238  const unsigned int active_fe_index) const
2239 {
2240  AssertIndexRange(quad_index, mapping_info.cell_data.size());
2241  return mapping_info.cell_data[quad_index]
2242  .descriptor[active_fe_index]
2243  .quadrature;
2244 }
2245 
2246 
2247 
2248 template <int dim, typename Number, typename VectorizedArrayType>
2249 inline const Quadrature<dim - 1> &
2251  const unsigned int quad_index,
2252  const unsigned int active_fe_index) const
2253 {
2254  AssertIndexRange(quad_index, mapping_info.face_data.size());
2255  return mapping_info.face_data[quad_index]
2256  .descriptor[active_fe_index]
2257  .quadrature;
2258 }
2259 
2260 
2261 
2262 template <int dim, typename Number, typename VectorizedArrayType>
2263 inline unsigned int
2265  const unsigned int macro_cell) const
2266 {
2267  AssertIndexRange(0, dof_info.size());
2268  AssertIndexRange(macro_cell, dof_info[0].cell_active_fe_index.size());
2269  if (dof_info[0].cell_active_fe_index.empty())
2270  return 0;
2271  else
2272  return dof_info[0].cell_active_fe_index[macro_cell];
2273 }
2274 
2275 
2276 
2277 template <int dim, typename Number, typename VectorizedArrayType>
2278 inline std::pair<unsigned int, unsigned int>
2280  const unsigned int macro_face) const
2281 {
2282  AssertIndexRange(macro_face, face_info.faces.size());
2283  if (dof_info[0].cell_active_fe_index.empty())
2284  return std::make_pair(0U, 0U);
2285 
2286  std::pair<unsigned int, unsigned int> result;
2287  for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements &&
2288  face_info.faces[macro_face].cells_interior[v] !=
2290  ++v)
2291  result.first = std::max(
2292  result.first,
2293  dof_info[0]
2294  .cell_active_fe_index[face_info.faces[macro_face].cells_interior[v]]);
2295  if (face_info.faces[macro_face].cells_exterior[0] !=
2297  for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements &&
2298  face_info.faces[macro_face].cells_exterior[v] !=
2300  ++v)
2301  result.second = std::max(
2302  result.first,
2303  dof_info[0]
2304  .cell_active_fe_index[face_info.faces[macro_face].cells_exterior[v]]);
2305  else
2306  result.second = numbers::invalid_unsigned_int;
2307  return result;
2308 }
2309 
2310 
2311 
2312 template <int dim, typename Number, typename VectorizedArrayType>
2313 inline bool
2315 {
2316  return indices_are_initialized;
2317 }
2318 
2319 
2320 
2321 template <int dim, typename Number, typename VectorizedArrayType>
2322 inline bool
2324 {
2325  return mapping_is_initialized;
2326 }
2327 
2328 
2329 
2330 template <int dim, typename Number, typename VectorizedArrayType>
2333 {
2334  using list_type =
2335  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
2336  list_type &data = scratch_pad.get();
2337  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
2338  if (it->first == false)
2339  {
2340  it->first = true;
2341  return &it->second;
2342  }
2343  data.emplace_front(true, AlignedVector<VectorizedArrayType>());
2344  return &data.front().second;
2345 }
2346 
2347 
2348 
2349 template <int dim, typename Number, typename VectorizedArrayType>
2350 void
2352  const AlignedVector<VectorizedArrayType> *scratch) const
2353 {
2354  using list_type =
2355  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
2356  list_type &data = scratch_pad.get();
2357  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
2358  if (&it->second == scratch)
2359  {
2360  Assert(it->first == true, ExcInternalError());
2361  it->first = false;
2362  return;
2363  }
2364  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2365 }
2366 
2367 
2368 
2369 template <int dim, typename Number, typename VectorizedArrayType>
2373 {
2374  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
2376  it != scratch_pad_non_threadsafe.end();
2377  ++it)
2378  if (it->first == false)
2379  {
2380  it->first = true;
2381  return &it->second;
2382  }
2383  scratch_pad_non_threadsafe.push_front(
2384  std::make_pair(true, AlignedVector<Number>()));
2385  return &scratch_pad_non_threadsafe.front().second;
2386 }
2387 
2388 
2389 
2390 template <int dim, typename Number, typename VectorizedArrayType>
2391 void
2394  const AlignedVector<Number> *scratch) const
2395 {
2396  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
2398  it != scratch_pad_non_threadsafe.end();
2399  ++it)
2400  if (&it->second == scratch)
2401  {
2402  Assert(it->first == true, ExcInternalError());
2403  it->first = false;
2404  return;
2405  }
2406  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2407 }
2408 
2409 
2410 
2411 // ------------------------------ reinit functions ---------------------------
2412 
2413 namespace internal
2414 {
2415  namespace MatrixFreeImplementation
2416  {
2417  template <typename DoFHandlerType>
2418  inline std::vector<IndexSet>
2419  extract_locally_owned_index_sets(
2420  const std::vector<const DoFHandlerType *> &dofh,
2421  const unsigned int level)
2422  {
2423  std::vector<IndexSet> locally_owned_set;
2424  locally_owned_set.reserve(dofh.size());
2425  for (unsigned int j = 0; j < dofh.size(); j++)
2426  if (level == numbers::invalid_unsigned_int)
2427  locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
2428  else
2429  AssertThrow(false, ExcNotImplemented());
2430  return locally_owned_set;
2431  }
2432 
2433  template <int dim, int spacedim>
2434  inline std::vector<IndexSet>
2435  extract_locally_owned_index_sets(
2436  const std::vector<const ::DoFHandler<dim, spacedim> *> &dofh,
2437  const unsigned int level)
2438  {
2439  std::vector<IndexSet> locally_owned_set;
2440  locally_owned_set.reserve(dofh.size());
2441  for (unsigned int j = 0; j < dofh.size(); j++)
2442  if (level == numbers::invalid_unsigned_int)
2443  locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
2444  else
2445  locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
2446  return locally_owned_set;
2447  }
2448  } // namespace MatrixFreeImplementation
2449 } // namespace internal
2450 
2451 
2452 
2453 template <int dim, typename Number, typename VectorizedArrayType>
2454 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2455 void
2457  const DoFHandlerType & dof_handler,
2458  const AffineConstraints<number2> &constraints_in,
2459  const QuadratureType & quad,
2461  &additional_data)
2462 {
2463  std::vector<const DoFHandlerType *> dof_handlers;
2464  std::vector<const AffineConstraints<number2> *> constraints;
2465  std::vector<QuadratureType> quads;
2466 
2467  dof_handlers.push_back(&dof_handler);
2468  constraints.push_back(&constraints_in);
2469  quads.push_back(quad);
2470 
2471  std::vector<IndexSet> locally_owned_sets =
2472  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2473  dof_handlers, additional_data.level_mg_handler);
2474 
2475  std::vector<hp::QCollection<1>> quad_hp;
2476  quad_hp.emplace_back(quad);
2477 
2479  dof_handlers,
2480  constraints,
2481  locally_owned_sets,
2482  quad_hp,
2483  additional_data);
2484 }
2485 
2486 
2487 
2488 template <int dim, typename Number, typename VectorizedArrayType>
2489 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2490 void
2492  const Mapping<dim> & mapping,
2493  const DoFHandlerType & dof_handler,
2494  const AffineConstraints<number2> &constraints_in,
2495  const QuadratureType & quad,
2497  &additional_data)
2498 {
2499  std::vector<const DoFHandlerType *> dof_handlers;
2500  std::vector<const AffineConstraints<number2> *> constraints;
2501 
2502  dof_handlers.push_back(&dof_handler);
2503  constraints.push_back(&constraints_in);
2504 
2505  std::vector<IndexSet> locally_owned_sets =
2506  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2507  dof_handlers, additional_data.level_mg_handler);
2508 
2509  std::vector<hp::QCollection<1>> quad_hp;
2510  quad_hp.emplace_back(quad);
2511 
2512  internal_reinit(mapping,
2513  dof_handlers,
2514  constraints,
2515  locally_owned_sets,
2516  quad_hp,
2517  additional_data);
2518 }
2519 
2520 
2521 
2522 template <int dim, typename Number, typename VectorizedArrayType>
2523 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2524 void
2526  const std::vector<const DoFHandlerType *> & dof_handler,
2527  const std::vector<const AffineConstraints<number2> *> &constraint,
2528  const std::vector<QuadratureType> & quad,
2530  &additional_data)
2531 {
2532  std::vector<IndexSet> locally_owned_set =
2533  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2534  dof_handler, additional_data.level_mg_handler);
2535  std::vector<hp::QCollection<1>> quad_hp;
2536  for (unsigned int q = 0; q < quad.size(); ++q)
2537  quad_hp.emplace_back(quad[q]);
2539  dof_handler,
2540  constraint,
2541  locally_owned_set,
2542  quad_hp,
2543  additional_data);
2544 }
2545 
2546 
2547 
2548 template <int dim, typename Number, typename VectorizedArrayType>
2549 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2550 void
2552  const std::vector<const DoFHandlerType *> & dof_handler,
2553  const std::vector<const AffineConstraints<number2> *> &constraint,
2554  const QuadratureType & quad,
2556  &additional_data)
2557 {
2558  std::vector<IndexSet> locally_owned_set =
2559  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2560  dof_handler, additional_data.level_mg_handler);
2561  std::vector<hp::QCollection<1>> quad_hp;
2562  quad_hp.emplace_back(quad);
2564  dof_handler,
2565  constraint,
2566  locally_owned_set,
2567  quad_hp,
2568  additional_data);
2569 }
2570 
2571 
2572 
2573 template <int dim, typename Number, typename VectorizedArrayType>
2574 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2575 void
2577  const Mapping<dim> & mapping,
2578  const std::vector<const DoFHandlerType *> & dof_handler,
2579  const std::vector<const AffineConstraints<number2> *> &constraint,
2580  const QuadratureType & quad,
2582  &additional_data)
2583 {
2584  std::vector<IndexSet> locally_owned_set =
2585  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2586  dof_handler, additional_data.level_mg_handler);
2587  std::vector<hp::QCollection<1>> quad_hp;
2588  quad_hp.emplace_back(quad);
2589  internal_reinit(mapping,
2590  dof_handler,
2591  constraint,
2592  locally_owned_set,
2593  quad_hp,
2594  additional_data);
2595 }
2596 
2597 
2598 
2599 template <int dim, typename Number, typename VectorizedArrayType>
2600 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2601 void
2603  const Mapping<dim> & mapping,
2604  const std::vector<const DoFHandlerType *> & dof_handler,
2605  const std::vector<const AffineConstraints<number2> *> &constraint,
2606  const std::vector<QuadratureType> & quad,
2608  &additional_data)
2609 {
2610  std::vector<IndexSet> locally_owned_set =
2611  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2612  dof_handler, additional_data.level_mg_handler);
2613  std::vector<hp::QCollection<1>> quad_hp;
2614  for (unsigned int q = 0; q < quad.size(); ++q)
2615  quad_hp.emplace_back(quad[q]);
2616  internal_reinit(mapping,
2617  dof_handler,
2618  constraint,
2619  locally_owned_set,
2620  quad_hp,
2621  additional_data);
2622 }
2623 
2624 
2625 
2626 template <int dim, typename Number, typename VectorizedArrayType>
2627 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2628 void
2630  const Mapping<dim> & mapping,
2631  const std::vector<const DoFHandlerType *> & dof_handler,
2632  const std::vector<const AffineConstraints<number2> *> &constraint,
2633  const std::vector<IndexSet> & locally_owned_set,
2634  const std::vector<QuadratureType> & quad,
2636  &additional_data)
2637 {
2638  // find out whether we use a hp Quadrature or a standard quadrature
2639  std::vector<hp::QCollection<1>> quad_hp;
2640  for (unsigned int q = 0; q < quad.size(); ++q)
2641  quad_hp.emplace_back(quad[q]);
2642  internal_reinit(mapping,
2643  dof_handler,
2644  constraint,
2645  locally_owned_set,
2646  quad_hp,
2647  additional_data);
2648 }
2649 
2650 
2651 
2652 // ------------------------------ implementation of loops --------------------
2653 
2654 // internal helper functions that define how to call MPI data exchange
2655 // functions: for generic vectors, do nothing at all. For distributed vectors,
2656 // call update_ghost_values_start function and so on. If we have collections
2657 // of vectors, just do the individual functions of the components. In order to
2658 // keep ghost values consistent (whether we are in read or write mode), we
2659 // also reset the values at the end. the whole situation is a bit complicated
2660 // by the fact that we need to treat block vectors differently, which use some
2661 // additional helper functions to select the blocks and template magic.
2662 namespace internal
2663 {
2667  template <int dim, typename Number, typename VectorizedArrayType>
2668  struct VectorDataExchange
2669  {
2670  // An arbitrary shift for communication to reduce the risk for accidental
2671  // interaction with other open communications that a user program might
2672  // set up
2673  static constexpr unsigned int channel_shift = 103;
2674 
2675 
2676 
2681  VectorDataExchange(
2682  const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
2683  const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
2684  DataAccessOnFaces vector_face_access,
2685  const unsigned int n_components)
2686  : matrix_free(matrix_free)
2687  , vector_face_access(
2688  matrix_free.get_task_info().face_partition_data.empty() ?
2690  DataAccessOnFaces::unspecified :
2691  vector_face_access)
2692  , ghosts_were_set(false)
2693 # ifdef DEAL_II_WITH_MPI
2694  , tmp_data(n_components)
2695  , requests(n_components)
2696 # endif
2697  {
2698  (void)n_components;
2699  if (this->vector_face_access !=
2702  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
2704  matrix_free.get_dof_info(c).vector_partitioner_face_variants.size(),
2705  3);
2706  }
2707 
2708 
2709 
2713  ~VectorDataExchange() // NOLINT
2714  {
2715 # ifdef DEAL_II_WITH_MPI
2716  for (unsigned int i = 0; i < tmp_data.size(); ++i)
2717  if (tmp_data[i] != nullptr)
2718  matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
2719 # endif
2720  }
2721 
2722 
2723 
2728  template <typename VectorType>
2729  unsigned int
2730  find_vector_in_mf(const VectorType &vec,
2731  const bool check_global_compatibility = true) const
2732  {
2733  unsigned int mf_component = numbers::invalid_unsigned_int;
2734  (void)check_global_compatibility;
2735  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
2736  if (
2737 # ifdef DEBUG
2738  check_global_compatibility ?
2739  vec.get_partitioner()->is_globally_compatible(
2740  *matrix_free.get_dof_info(c).vector_partitioner) :
2741 # endif
2742  vec.get_partitioner()->is_compatible(
2743  *matrix_free.get_dof_info(c).vector_partitioner))
2744  {
2745  mf_component = c;
2746  break;
2747  }
2748  return mf_component;
2749  }
2750 
2751 
2752 
2758  get_partitioner(const unsigned int mf_component) const
2759  {
2760  AssertDimension(matrix_free.get_dof_info(mf_component)
2761  .vector_partitioner_face_variants.size(),
2762  3);
2763  if (vector_face_access ==
2766  return *matrix_free.get_dof_info(mf_component)
2767  .vector_partitioner_face_variants[0];
2768  else if (vector_face_access ==
2771  return *matrix_free.get_dof_info(mf_component)
2772  .vector_partitioner_face_variants[1];
2773  else
2774  return *matrix_free.get_dof_info(mf_component)
2775  .vector_partitioner_face_variants[2];
2776  }
2777 
2778 
2779 
2783  template <typename VectorType,
2784  typename std::enable_if<is_serial_or_dummy<VectorType>::value,
2785  VectorType>::type * = nullptr>
2786  void
2787  update_ghost_values_start(const unsigned int /*component_in_block_vector*/,
2788  const VectorType & /*vec*/)
2789  {}
2790 
2791 
2796  template <typename VectorType,
2797  typename std::enable_if<
2798  !has_update_ghost_values_start<VectorType>::value &&
2799  !is_serial_or_dummy<VectorType>::value,
2800  VectorType>::type * = nullptr>
2801  void
2802  update_ghost_values_start(const unsigned int component_in_block_vector,
2803  const VectorType & vec)
2804  {
2805  (void)component_in_block_vector;
2806  bool ghosts_set = vec.has_ghost_elements();
2807  if (ghosts_set)
2808  ghosts_were_set = true;
2809 
2810  vec.update_ghost_values();
2811  }
2812 
2813 
2814 
2820  template <typename VectorType,
2821  typename std::enable_if<
2822  has_update_ghost_values_start<VectorType>::value &&
2823  !has_exchange_on_subset<VectorType>::value,
2824  VectorType>::type * = nullptr>
2825  void
2826  update_ghost_values_start(const unsigned int component_in_block_vector,
2827  const VectorType & vec)
2828  {
2829  (void)component_in_block_vector;
2830  bool ghosts_set = vec.has_ghost_elements();
2831  if (ghosts_set)
2832  ghosts_were_set = true;
2833 
2834  vec.update_ghost_values_start(component_in_block_vector + channel_shift);
2835  }
2836 
2837 
2838 
2845  template <typename VectorType,
2846  typename std::enable_if<
2847  has_update_ghost_values_start<VectorType>::value &&
2848  has_exchange_on_subset<VectorType>::value,
2849  VectorType>::type * = nullptr>
2850  void
2851  update_ghost_values_start(const unsigned int component_in_block_vector,
2852  const VectorType & vec)
2853  {
2854  static_assert(
2855  std::is_same<Number, typename VectorType::value_type>::value,
2856  "Type mismatch between VectorType and VectorDataExchange");
2857  (void)component_in_block_vector;
2858  bool ghosts_set = vec.has_ghost_elements();
2859  if (ghosts_set)
2860  ghosts_were_set = true;
2861  if (vector_face_access ==
2864  vec.size() == 0)
2865  vec.update_ghost_values_start(component_in_block_vector +
2866  channel_shift);
2867  else
2868  {
2869 # ifdef DEAL_II_WITH_MPI
2870  const unsigned int mf_component = find_vector_in_mf(vec);
2871  if (&get_partitioner(mf_component) ==
2872  matrix_free.get_dof_info(mf_component).vector_partitioner.get())
2873  {
2874  vec.update_ghost_values_start(component_in_block_vector +
2875  channel_shift);
2876  return;
2877  }
2878 
2879  const Utilities::MPI::Partitioner &part =
2880  get_partitioner(mf_component);
2881  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0)
2882  return;
2883 
2884  tmp_data[component_in_block_vector] =
2885  matrix_free.acquire_scratch_data_non_threadsafe();
2886  tmp_data[component_in_block_vector]->resize_fast(
2887  part.n_import_indices());
2888  AssertDimension(requests.size(), tmp_data.size());
2889 
2891  component_in_block_vector + channel_shift,
2892  ArrayView<const Number>(vec.begin(), part.local_size()),
2893  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
2894  part.n_import_indices()),
2895  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
2896  vec.get_partitioner()->local_size(),
2897  vec.get_partitioner()->n_ghost_indices()),
2898  this->requests[component_in_block_vector]);
2899 # endif
2900  }
2901  }
2902 
2903 
2904 
2909  template <
2910  typename VectorType,
2911  typename std::enable_if<!has_update_ghost_values_start<VectorType>::value,
2912  VectorType>::type * = nullptr>
2913  void
2914  update_ghost_values_finish(const unsigned int /*component_in_block_vector*/,
2915  const VectorType & /*vec*/)
2916  {}
2917 
2918 
2919 
2925  template <typename VectorType,
2926  typename std::enable_if<
2927  has_update_ghost_values_start<VectorType>::value &&
2928  !has_exchange_on_subset<VectorType>::value,
2929  VectorType>::type * = nullptr>
2930  void
2931  update_ghost_values_finish(const unsigned int component_in_block_vector,
2932  const VectorType & vec)
2933  {
2934  (void)component_in_block_vector;
2935  vec.update_ghost_values_finish();
2936  }
2937 
2938 
2939 
2946  template <typename VectorType,
2947  typename std::enable_if<
2948  has_update_ghost_values_start<VectorType>::value &&
2949  has_exchange_on_subset<VectorType>::value,
2950  VectorType>::type * = nullptr>
2951  void
2952  update_ghost_values_finish(const unsigned int component_in_block_vector,
2953  const VectorType & vec)
2954  {
2955  static_assert(
2956  std::is_same<Number, typename VectorType::value_type>::value,
2957  "Type mismatch between VectorType and VectorDataExchange");
2958  (void)component_in_block_vector;
2959  if (vector_face_access ==
2962  vec.size() == 0)
2963  vec.update_ghost_values_finish();
2964  else
2965  {
2966 # ifdef DEAL_II_WITH_MPI
2967 
2968  AssertIndexRange(component_in_block_vector, tmp_data.size());
2969  AssertDimension(requests.size(), tmp_data.size());
2970 
2971  const unsigned int mf_component = find_vector_in_mf(vec);
2972  const Utilities::MPI::Partitioner &part =
2973  get_partitioner(mf_component);
2974  if (&part ==
2975  matrix_free.get_dof_info(mf_component).vector_partitioner.get())
2976  {
2977  vec.update_ghost_values_finish();
2978  return;
2979  }
2980 
2981  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0)
2982  return;
2983 
2985  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
2986  vec.get_partitioner()->local_size(),
2987  vec.get_partitioner()->n_ghost_indices()),
2988  this->requests[component_in_block_vector]);
2989 
2990  matrix_free.release_scratch_data_non_threadsafe(
2991  tmp_data[component_in_block_vector]);
2992  tmp_data[component_in_block_vector] = nullptr;
2993 
2994  // let vector know that ghosts are being updated and we can read from
2995  // them
2996  vec.set_ghost_state(true);
2997 # endif
2998  }
2999  }
3000 
3001 
3002 
3006  template <typename VectorType,
3007  typename std::enable_if<is_serial_or_dummy<VectorType>::value,
3008  VectorType>::type * = nullptr>
3009  void
3010  compress_start(const unsigned int /*component_in_block_vector*/,
3011  VectorType & /*vec*/)
3012  {}
3013 
3014 
3015 
3020  template <typename VectorType,
3021  typename std::enable_if<!has_compress_start<VectorType>::value &&
3022  !is_serial_or_dummy<VectorType>::value,
3023  VectorType>::type * = nullptr>
3024  void
3025  compress_start(const unsigned int component_in_block_vector,
3026  VectorType & vec)
3027  {
3028  (void)component_in_block_vector;
3029  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3030  vec.compress(::VectorOperation::add);
3031  }
3032 
3033 
3034 
3040  template <
3041  typename VectorType,
3042  typename std::enable_if<has_compress_start<VectorType>::value &&
3043  !has_exchange_on_subset<VectorType>::value,
3044  VectorType>::type * = nullptr>
3045  void
3046  compress_start(const unsigned int component_in_block_vector,
3047  VectorType & vec)
3048  {
3049  (void)component_in_block_vector;
3050  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3051  vec.compress_start(component_in_block_vector + channel_shift);
3052  }
3053 
3054 
3055 
3062  template <
3063  typename VectorType,
3064  typename std::enable_if<has_compress_start<VectorType>::value &&
3065  has_exchange_on_subset<VectorType>::value,
3066  VectorType>::type * = nullptr>
3067  void
3068  compress_start(const unsigned int component_in_block_vector,
3069  VectorType & vec)
3070  {
3071  static_assert(
3072  std::is_same<Number, typename VectorType::value_type>::value,
3073  "Type mismatch between VectorType and VectorDataExchange");
3074  (void)component_in_block_vector;
3075  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3076  if (vector_face_access ==
3079  vec.size() == 0)
3080  vec.compress_start(component_in_block_vector + channel_shift);
3081  else
3082  {
3083 # ifdef DEAL_II_WITH_MPI
3084 
3085  const unsigned int mf_component = find_vector_in_mf(vec);
3086  const Utilities::MPI::Partitioner &part =
3087  get_partitioner(mf_component);
3088  if (&part ==
3089  matrix_free.get_dof_info(mf_component).vector_partitioner.get())
3090  {
3091  vec.compress_start(component_in_block_vector + channel_shift);
3092  return;
3093  }
3094 
3095  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0)
3096  return;
3097 
3098  tmp_data[component_in_block_vector] =
3099  matrix_free.acquire_scratch_data_non_threadsafe();
3100  tmp_data[component_in_block_vector]->resize_fast(
3101  part.n_import_indices());
3102  AssertDimension(requests.size(), tmp_data.size());
3103 
3106  component_in_block_vector + channel_shift,
3107  ArrayView<Number>(vec.begin() + vec.get_partitioner()->local_size(),
3108  vec.get_partitioner()->n_ghost_indices()),
3109  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3110  part.n_import_indices()),
3111  this->requests[component_in_block_vector]);
3112 # endif
3113  }
3114  }
3115 
3116 
3117 
3122  template <typename VectorType,
3123  typename std::enable_if<!has_compress_start<VectorType>::value,
3124  VectorType>::type * = nullptr>
3125  void
3126  compress_finish(const unsigned int /*component_in_block_vector*/,
3127  VectorType & /*vec*/)
3128  {}
3129 
3130 
3131 
3137  template <
3138  typename VectorType,
3139  typename std::enable_if<has_compress_start<VectorType>::value &&
3140  !has_exchange_on_subset<VectorType>::value,
3141  VectorType>::type * = nullptr>
3142  void
3143  compress_finish(const unsigned int component_in_block_vector,
3144  VectorType & vec)
3145  {
3146  (void)component_in_block_vector;
3147  vec.compress_finish(::VectorOperation::add);
3148  }
3149 
3150 
3151 
3158  template <
3159  typename VectorType,
3160  typename std::enable_if<has_compress_start<VectorType>::value &&
3161  has_exchange_on_subset<VectorType>::value,
3162  VectorType>::type * = nullptr>
3163  void
3164  compress_finish(const unsigned int component_in_block_vector,
3165  VectorType & vec)
3166  {
3167  static_assert(
3168  std::is_same<Number, typename VectorType::value_type>::value,
3169  "Type mismatch between VectorType and VectorDataExchange");
3170  (void)component_in_block_vector;
3171  if (vector_face_access ==
3174  vec.size() == 0)
3175  vec.compress_finish(::VectorOperation::add);
3176  else
3177  {
3178 # ifdef DEAL_II_WITH_MPI
3179  AssertIndexRange(component_in_block_vector, tmp_data.size());
3180  AssertDimension(requests.size(), tmp_data.size());
3181 
3182  const unsigned int mf_component = find_vector_in_mf(vec);
3183 
3184  const Utilities::MPI::Partitioner &part =
3185  get_partitioner(mf_component);
3186  if (&part ==
3187  matrix_free.get_dof_info(mf_component).vector_partitioner.get())
3188  {
3189  vec.compress_finish(::VectorOperation::add);
3190  return;
3191  }
3192 
3193  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0)
3194  return;
3195 
3199  tmp_data[component_in_block_vector]->begin(),
3200  part.n_import_indices()),
3201  ArrayView<Number>(vec.begin(), part.local_size()),
3202  ArrayView<Number>(vec.begin() + vec.get_partitioner()->local_size(),
3203  vec.get_partitioner()->n_ghost_indices()),
3204  this->requests[component_in_block_vector]);
3205 
3206  matrix_free.release_scratch_data_non_threadsafe(
3207  tmp_data[component_in_block_vector]);
3208  tmp_data[component_in_block_vector] = nullptr;
3209 # endif
3210  }
3211  }
3212 
3213 
3214 
3218  template <typename VectorType,
3219  typename std::enable_if<is_serial_or_dummy<VectorType>::value,
3220  VectorType>::type * = nullptr>
3221  void
3222  reset_ghost_values(const VectorType & /*vec*/) const
3223  {}
3224 
3225 
3226 
3231  template <
3232  typename VectorType,
3233  typename std::enable_if<!has_exchange_on_subset<VectorType>::value &&
3234  !is_serial_or_dummy<VectorType>::value,
3235  VectorType>::type * = nullptr>
3236  void
3237  reset_ghost_values(const VectorType &vec) const
3238  {
3239  if (ghosts_were_set == true)
3240  return;
3241 
3242  vec.zero_out_ghosts();
3243  }
3244 
3245 
3246 
3252  template <typename VectorType,
3253  typename std::enable_if<has_exchange_on_subset<VectorType>::value,
3254  VectorType>::type * = nullptr>
3255  void
3256  reset_ghost_values(const VectorType &vec) const
3257  {
3258  static_assert(
3259  std::is_same<Number, typename VectorType::value_type>::value,
3260  "Type mismatch between VectorType and VectorDataExchange");
3261  if (ghosts_were_set == true)
3262  return;
3263 
3264  if (vector_face_access ==
3267  vec.size() == 0)
3268  vec.zero_out_ghosts();
3269  else
3270  {
3271 # ifdef DEAL_II_WITH_MPI
3272  AssertDimension(requests.size(), tmp_data.size());
3273 
3274  const unsigned int mf_component = find_vector_in_mf(vec);
3275  const Utilities::MPI::Partitioner &part =
3276  get_partitioner(mf_component);
3277  if (&part ==
3278  matrix_free.get_dof_info(mf_component).vector_partitioner.get())
3279  vec.zero_out_ghosts();
3280  else if (part.n_ghost_indices() > 0)
3281  {
3282  for (std::vector<std::pair<unsigned int, unsigned int>>::
3283  const_iterator my_ghosts =
3285  my_ghosts !=
3287  ++my_ghosts)
3288  for (unsigned int j = my_ghosts->first; j < my_ghosts->second;
3289  j++)
3290  {
3292  vec)
3293  .local_element(j + part.local_size()) = 0.;
3294  }
3295  }
3296 
3297  // let vector know that it's not ghosted anymore
3298  vec.set_ghost_state(false);
3299 # endif
3300  }
3301  }
3302 
3303 
3304 
3310  template <typename VectorType,
3311  typename std::enable_if<has_exchange_on_subset<VectorType>::value,
3312  VectorType>::type * = nullptr>
3313  void
3314  zero_vector_region(const unsigned int range_index, VectorType &vec) const
3315  {
3316  static_assert(
3317  std::is_same<Number, typename VectorType::value_type>::value,
3318  "Type mismatch between VectorType and VectorDataExchange");
3319  if (range_index == numbers::invalid_unsigned_int)
3320  vec = Number();
3321  else
3322  {
3323  const unsigned int mf_component = find_vector_in_mf(vec, false);
3325  matrix_free.get_dof_info(mf_component);
3326  Assert(dof_info.vector_zero_range_list_index.empty() == false,
3327  ExcNotInitialized());
3328 
3329  Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3330  ExcInternalError());
3331  AssertIndexRange(range_index,
3332  dof_info.vector_zero_range_list_index.size() - 1);
3333  for (unsigned int id =
3334  dof_info.vector_zero_range_list_index[range_index];
3335  id != dof_info.vector_zero_range_list_index[range_index + 1];
3336  ++id)
3337  {
3338  const unsigned int start_pos =
3339  dof_info.vector_zero_range_list[id] *
3341  const unsigned int end_pos =
3342  std::min((dof_info.vector_zero_range_list[id] + 1) *
3344  chunk_size_zero_vector,
3345  dof_info.vector_partitioner->local_size() +
3346  dof_info.vector_partitioner->n_ghost_indices());
3347  std::memset(vec.begin() + start_pos,
3348  0,
3349  (end_pos - start_pos) * sizeof(Number));
3350  }
3351  }
3352  }
3353 
3354 
3355 
3361  template <
3362  typename VectorType,
3363  typename std::enable_if<!has_exchange_on_subset<VectorType>::value,
3364  VectorType>::type * = nullptr,
3365  typename VectorType::value_type * = nullptr>
3366  void
3367  zero_vector_region(const unsigned int range_index, VectorType &vec) const
3368  {
3369  if (range_index == numbers::invalid_unsigned_int)
3370  vec = typename VectorType::value_type();
3371  else
3372  {
3373  Assert(false, ExcNotImplemented());
3374  }
3375  }
3376 
3377 
3378 
3383  void
3384  zero_vector_region(const unsigned int, ...) const
3385  {
3386  Assert(false,
3387  ExcNotImplemented("Zeroing is only implemented for vector types "
3388  "which provide VectorType::value_type"));
3389  }
3390 
3391 
3392 
3393  const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
3394  const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3395  DataAccessOnFaces vector_face_access;
3396  bool ghosts_were_set;
3397 # ifdef DEAL_II_WITH_MPI
3398  std::vector<AlignedVector<Number> *> tmp_data;
3399  std::vector<std::vector<MPI_Request>> requests;
3400 # endif
3401  }; // VectorDataExchange
3402 
3403  template <typename VectorStruct>
3404  unsigned int
3405  n_components(const VectorStruct &vec);
3406 
3407  template <typename VectorStruct>
3408  unsigned int
3409  n_components_block(const VectorStruct &vec,
3410  std::integral_constant<bool, true>)
3411  {
3412  unsigned int components = 0;
3413  for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
3414  components += n_components(vec.block(bl));
3415  return components;
3416  }
3417 
3418  template <typename VectorStruct>
3419  unsigned int
3420  n_components_block(const VectorStruct &, std::integral_constant<bool, false>)
3421  {
3422  return 1;
3423  }
3424 
3425  template <typename VectorStruct>
3426  unsigned int
3427  n_components(const VectorStruct &vec)
3428  {
3429  return n_components_block(
3430  vec, std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3431  }
3432 
3433  template <typename VectorStruct>
3434  inline unsigned int
3435  n_components(const std::vector<VectorStruct> &vec)
3436  {
3437  unsigned int components = 0;
3438  for (unsigned int comp = 0; comp < vec.size(); comp++)
3439  components += n_components_block(
3440  vec[comp],
3441  std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3442  return components;
3443  }
3444 
3445  template <typename VectorStruct>
3446  inline unsigned int
3447  n_components(const std::vector<VectorStruct *> &vec)
3448  {
3449  unsigned int components = 0;
3450  for (unsigned int comp = 0; comp < vec.size(); comp++)
3451  components += n_components_block(
3452  *vec[comp],
3453  std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3454  return components;
3455  }
3456 
3457 
3458 
3459  // A helper function to identify block vectors with many components where we
3460  // should not try to overlap computations and communication because there
3461  // would be too many outstanding communication requests.
3462 
3463  // default value for vectors that do not have communication_block_size
3464  template <
3465  typename VectorStruct,
3466  typename std::enable_if<!has_communication_block_size<VectorStruct>::value,
3467  VectorStruct>::type * = nullptr>
3468  constexpr unsigned int
3469  get_communication_block_size(const VectorStruct &)
3470  {
3472  }
3473 
3474 
3475 
3476  template <
3477  typename VectorStruct,
3478  typename std::enable_if<has_communication_block_size<VectorStruct>::value,
3479  VectorStruct>::type * = nullptr>
3480  constexpr unsigned int
3481  get_communication_block_size(const VectorStruct &)
3482  {
3483  return VectorStruct::communication_block_size;
3484  }
3485 
3486 
3487 
3488  // --------------------------------------------------------------------------
3489  // below we have wrappers to distinguish between block and non-block vectors.
3490  // --------------------------------------------------------------------------
3491 
3492  //
3493  // update_ghost_values_start
3494  //
3495 
3496  // update_ghost_values for block vectors
3497  template <int dim,
3498  typename VectorStruct,
3499  typename Number,
3500  typename VectorizedArrayType,
3501  typename std::enable_if<IsBlockVector<VectorStruct>::value,
3502  VectorStruct>::type * = nullptr>
3503  void
3504  update_ghost_values_start(
3505  const VectorStruct & vec,
3506  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3507  const unsigned int channel = 0)
3508  {
3509  if (get_communication_block_size(vec) < vec.n_blocks())
3510  {
3511  // don't forget to set ghosts_were_set, that otherwise happens
3512  // inside VectorDataExchange::update_ghost_values_start()
3513  exchanger.ghosts_were_set = vec.has_ghost_elements();
3514  vec.update_ghost_values();
3515  }
3516  else
3517  {
3518  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
3519  update_ghost_values_start(vec.block(i), exchanger, channel + i);
3520  }
3521  }
3522 
3523 
3524 
3525  // update_ghost_values for non-block vectors
3526  template <int dim,
3527  typename VectorStruct,
3528  typename Number,
3529  typename VectorizedArrayType,
3530  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
3531  VectorStruct>::type * = nullptr>
3532  void
3533  update_ghost_values_start(
3534  const VectorStruct & vec,
3535  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3536  const unsigned int channel = 0)
3537  {
3538  exchanger.update_ghost_values_start(channel, vec);
3539  }
3540 
3541 
3542 
3543  // update_ghost_values_start() for vector of vectors
3544  template <int dim,
3545  typename VectorStruct,
3546  typename Number,
3547  typename VectorizedArrayType>
3548  inline void
3549  update_ghost_values_start(
3550  const std::vector<VectorStruct> & vec,
3551  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3552  {
3553  unsigned int component_index = 0;
3554  for (unsigned int comp = 0; comp < vec.size(); comp++)
3555  {
3556  update_ghost_values_start(vec[comp], exchanger, component_index);
3557  component_index += n_components(vec[comp]);
3558  }
3559  }
3560 
3561 
3562 
3563  // update_ghost_values_start() for vector of pointers to vectors
3564  template <int dim,
3565  typename VectorStruct,
3566  typename Number,
3567  typename VectorizedArrayType>
3568  inline void
3569  update_ghost_values_start(
3570  const std::vector<VectorStruct *> & vec,
3571  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3572  {
3573  unsigned int component_index = 0;
3574  for (unsigned int comp = 0; comp < vec.size(); comp++)
3575  {
3576  update_ghost_values_start(*vec[comp], exchanger, component_index);
3577  component_index += n_components(*vec[comp]);
3578  }
3579  }
3580 
3581 
3582 
3583  //
3584  // update_ghost_values_finish
3585  //
3586 
3587  // for block vectors
3588  template <int dim,
3589  typename VectorStruct,
3590  typename Number,
3591  typename VectorizedArrayType,
3592  typename std::enable_if<IsBlockVector<VectorStruct>::value,
3593  VectorStruct>::type * = nullptr>
3594  void
3595  update_ghost_values_finish(
3596  const VectorStruct & vec,
3597  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3598  const unsigned int channel = 0)
3599  {
3600  if (get_communication_block_size(vec) < vec.n_blocks())
3601  {
3602  // do nothing, everything has already been completed in the _start()
3603  // call
3604  }
3605  else
3606  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
3607  update_ghost_values_finish(vec.block(i), exchanger, channel + i);
3608  }
3609 
3610 
3611 
3612  // for non-block vectors
3613  template <int dim,
3614  typename VectorStruct,
3615  typename Number,
3616  typename VectorizedArrayType,
3617  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
3618  VectorStruct>::type * = nullptr>
3619  void
3620  update_ghost_values_finish(
3621  const VectorStruct & vec,
3622  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3623  const unsigned int channel = 0)
3624  {
3625  exchanger.update_ghost_values_finish(channel, vec);
3626  }
3627 
3628 
3629 
3630  // for vector of vectors
3631  template <int dim,
3632  typename VectorStruct,
3633  typename Number,
3634  typename VectorizedArrayType>
3635  inline void
3636  update_ghost_values_finish(
3637  const std::vector<VectorStruct> & vec,
3638  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3639  {
3640  unsigned int component_index = 0;
3641  for (unsigned int comp = 0; comp < vec.size(); comp++)
3642  {
3643  update_ghost_values_finish(vec[comp], exchanger, component_index);
3644  component_index += n_components(vec[comp]);
3645  }
3646  }
3647 
3648 
3649 
3650  // for vector of pointers to vectors
3651  template <int dim,
3652  typename VectorStruct,
3653  typename Number,
3654  typename VectorizedArrayType>
3655  inline void
3656  update_ghost_values_finish(
3657  const std::vector<VectorStruct *> & vec,
3658  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3659  {
3660  unsigned int component_index = 0;
3661  for (unsigned int comp = 0; comp < vec.size(); comp++)
3662  {
3663  update_ghost_values_finish(*vec[comp], exchanger, component_index);
3664  component_index += n_components(*vec[comp]);
3665  }
3666  }
3667 
3668 
3669 
3670  //
3671  // compress_start
3672  //
3673 
3674  // for block vectors
3675  template <int dim,
3676  typename VectorStruct,
3677  typename Number,
3678  typename VectorizedArrayType,
3679  typename std::enable_if<IsBlockVector<VectorStruct>::value,
3680  VectorStruct>::type * = nullptr>
3681  inline void
3682  compress_start(
3683  VectorStruct & vec,
3684  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3685  const unsigned int channel = 0)
3686  {
3687  if (get_communication_block_size(vec) < vec.n_blocks())
3688  vec.compress(::VectorOperation::add);
3689  else
3690  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
3691  compress_start(vec.block(i), exchanger, channel + i);
3692  }
3693 
3694 
3695 
3696  // for non-block vectors
3697  template <int dim,
3698  typename VectorStruct,
3699  typename Number,
3700  typename VectorizedArrayType,
3701  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
3702  VectorStruct>::type * = nullptr>
3703  inline void
3704  compress_start(
3705  VectorStruct & vec,
3706  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3707  const unsigned int channel = 0)
3708  {
3709  exchanger.compress_start(channel, vec);
3710  }
3711 
3712 
3713 
3714  // for std::vector of vectors
3715  template <int dim,
3716  typename VectorStruct,
3717  typename Number,
3718  typename VectorizedArrayType>
3719  inline void
3720  compress_start(
3721  std::vector<VectorStruct> & vec,
3722  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3723  {
3724  unsigned int component_index = 0;
3725  for (unsigned int comp = 0; comp < vec.size(); comp++)
3726  {
3727  compress_start(vec[comp], exchanger, component_index);
3728  component_index += n_components(vec[comp]);
3729  }
3730  }
3731 
3732 
3733 
3734  // for std::vector of pointer to vectors
3735  template <int dim,
3736  typename VectorStruct,
3737  typename Number,
3738  typename VectorizedArrayType>
3739  inline void
3740  compress_start(
3741  std::vector<VectorStruct *> & vec,
3742  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3743  {
3744  unsigned int component_index = 0;
3745  for (unsigned int comp = 0; comp < vec.size(); comp++)
3746  {
3747  compress_start(*vec[comp], exchanger, component_index);
3748  component_index += n_components(*vec[comp]);
3749  }
3750  }
3751 
3752 
3753 
3754  //
3755  // compress_finish
3756  //
3757 
3758  // for block vectors
3759  template <int dim,
3760  typename VectorStruct,
3761  typename Number,
3762  typename VectorizedArrayType,
3763  typename std::enable_if<IsBlockVector<VectorStruct>::value,
3764  VectorStruct>::type * = nullptr>
3765  inline void
3766  compress_finish(
3767  VectorStruct & vec,
3768  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3769  const unsigned int channel = 0)
3770  {
3771  if (get_communication_block_size(vec) < vec.n_blocks())
3772  {
3773  // do nothing, everything has already been completed in the _start()
3774  // call
3775  }
3776  else
3777  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
3778  compress_finish(vec.block(i), exchanger, channel + i);
3779  }
3780 
3781 
3782 
3783  // for non-block vectors
3784  template <int dim,
3785  typename VectorStruct,
3786  typename Number,
3787  typename VectorizedArrayType,
3788  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
3789  VectorStruct>::type * = nullptr>
3790  inline void
3791  compress_finish(
3792  VectorStruct & vec,
3793  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3794  const unsigned int channel = 0)
3795  {
3796  exchanger.compress_finish(channel, vec);
3797  }
3798 
3799 
3800 
3801  // for std::vector of vectors
3802  template <int dim,
3803  typename VectorStruct,
3804  typename Number,
3805  typename VectorizedArrayType>
3806  inline void
3807  compress_finish(
3808  std::vector<VectorStruct> & vec,
3809  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3810  {
3811  unsigned int component_index = 0;
3812  for (unsigned int comp = 0; comp < vec.size(); comp++)
3813  {
3814  compress_finish(vec[comp], exchanger, component_index);
3815  component_index += n_components(vec[comp]);
3816  }
3817  }
3818 
3819 
3820 
3821  // for std::vector of pointer to vectors
3822  template <int dim,
3823  typename VectorStruct,
3824  typename Number,
3825  typename VectorizedArrayType>
3826  inline void
3827  compress_finish(
3828  std::vector<VectorStruct *> & vec,
3829  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3830  {
3831  unsigned int component_index = 0;
3832  for (unsigned int comp = 0; comp < vec.size(); comp++)
3833  {
3834  compress_finish(*vec[comp], exchanger, component_index);
3835  component_index += n_components(*vec[comp]);
3836  }
3837  }
3838 
3839 
3840 
3841  //
3842  // reset_ghost_values:
3843  //
3844  // if the input vector did not have ghosts imported, clear them here again
3845  // in order to avoid subsequent operations e.g. in linear solvers to work
3846  // with ghosts all the time
3847  //
3848 
3849  // for block vectors
3850  template <int dim,
3851  typename VectorStruct,
3852  typename Number,
3853  typename VectorizedArrayType,
3854  typename std::enable_if<IsBlockVector<VectorStruct>::value,
3855  VectorStruct>::type * = nullptr>
3856  inline void
3857  reset_ghost_values(
3858  const VectorStruct & vec,
3859  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3860  {
3861  // return immediately if there is nothing to do.
3862  if (exchanger.ghosts_were_set == true)
3863  return;
3864 
3865  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
3866  reset_ghost_values(vec.block(i), exchanger);
3867  }
3868 
3869 
3870 
3871  // for non-block vectors
3872  template <int dim,
3873  typename VectorStruct,
3874  typename Number,
3875  typename VectorizedArrayType,
3876  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
3877  VectorStruct>::type * = nullptr>
3878  inline void
3879  reset_ghost_values(
3880  const VectorStruct & vec,
3881  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3882  {
3883  exchanger.reset_ghost_values(vec);
3884  }
3885 
3886 
3887 
3888  // for std::vector of vectors
3889  template <int dim,
3890  typename VectorStruct,
3891  typename Number,
3892  typename VectorizedArrayType>
3893  inline void
3894  reset_ghost_values(
3895  const std::vector<VectorStruct> & vec,
3896  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3897  {
3898  // return immediately if there is nothing to do.
3899  if (exchanger.ghosts_were_set == true)
3900  return;
3901 
3902  for (unsigned int comp = 0; comp < vec.size(); comp++)
3903  reset_ghost_values(vec[comp], exchanger);
3904  }
3905 
3906 
3907 
3908  // for std::vector of pointer to vectors
3909  template <int dim,
3910  typename VectorStruct,
3911  typename Number,
3912  typename VectorizedArrayType>
3913  inline void
3914  reset_ghost_values(
3915  const std::vector<VectorStruct *> & vec,
3916  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3917  {
3918  // return immediately if there is nothing to do.
3919  if (exchanger.ghosts_were_set == true)
3920  return;
3921 
3922  for (unsigned int comp = 0; comp < vec.size(); comp++)
3923  reset_ghost_values(*vec[comp], exchanger);
3924  }
3925 
3926 
3927 
3928  //
3929  // zero_vector_region
3930  //
3931 
3932  // for block vectors
3933  template <int dim,
3934  typename VectorStruct,
3935  typename Number,
3936  typename VectorizedArrayType,
3937  typename std::enable_if<IsBlockVector<VectorStruct>::value,
3938  VectorStruct>::type * = nullptr>
3939  inline void
3940  zero_vector_region(
3941  const unsigned int range_index,
3942  VectorStruct & vec,
3943  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3944  {
3945  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
3946  exchanger.zero_vector_region(range_index, vec.block(i));
3947  }
3948 
3949 
3950 
3951  // for non-block vectors
3952  template <int dim,
3953  typename VectorStruct,
3954  typename Number,
3955  typename VectorizedArrayType,
3956  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
3957  VectorStruct>::type * = nullptr>
3958  inline void
3959  zero_vector_region(
3960  const unsigned int range_index,
3961  VectorStruct & vec,
3962  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3963  {
3964  exchanger.zero_vector_region(range_index, vec);
3965  }
3966 
3967 
3968 
3969  // for std::vector of vectors
3970  template <int dim,
3971  typename VectorStruct,
3972  typename Number,
3973  typename VectorizedArrayType>
3974  inline void
3975  zero_vector_region(
3976  const unsigned int range_index,
3977  std::vector<VectorStruct> & vec,
3978  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3979  {
3980  for (unsigned int comp = 0; comp < vec.size(); comp++)
3981  zero_vector_region(range_index, vec[comp], exchanger);
3982  }
3983 
3984 
3985 
3986  // for std::vector of pointers to vectors
3987  template <int dim,
3988  typename VectorStruct,
3989  typename Number,
3990  typename VectorizedArrayType>
3991  inline void
3992  zero_vector_region(
3993  const unsigned int range_index,
3994  std::vector<VectorStruct *> & vec,
3995  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3996  {
3997  for (unsigned int comp = 0; comp < vec.size(); comp++)
3998  zero_vector_region(range_index, *vec[comp], exchanger);
3999  }
4000 
4001 
4002 
4003  namespace MatrixFreeFunctions
4004  {
4005  // struct to select between a const interface and a non-const interface
4006  // for MFWorker
4007  template <typename, typename, typename, typename, bool>
4008  struct InterfaceSelector
4009  {};
4010 
4011  // Version for constant functions
4012  template <typename MF,
4013  typename InVector,
4014  typename OutVector,
4015  typename Container>
4016  struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4017  {
4018  using function_type = void (Container::*)(
4019  const MF &,
4020  OutVector &,
4021  const InVector &,
4022  const std::pair<unsigned int, unsigned int> &) const;
4023  };
4024 
4025  // Version for non-constant functions
4026  template <typename MF,
4027  typename InVector,
4028  typename OutVector,
4029  typename Container>
4030  struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4031  {
4032  using function_type =
4033  void (Container::*)(const MF &,
4034  OutVector &,
4035  const InVector &,
4036  const std::pair<unsigned int, unsigned int> &);
4037  };
4038  } // namespace MatrixFreeFunctions
4039 
4040 
4041 
4042  // A implementation class for the worker object that runs the various
4043  // operations we want to perform during the matrix-free loop
4044  template <typename MF,
4045  typename InVector,
4046  typename OutVector,
4047  typename Container,
4048  bool is_constant>
4049  class MFWorker : public MFWorkerInterface
4050  {
4051  public:
4052  // An alias to make the arguments further down more readable
4053  using function_type = typename MatrixFreeFunctions::
4054  InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4055  function_type;
4056 
4057  // constructor, binds all the arguments to this class
4058  MFWorker(const MF & matrix_free,
4059  const InVector & src,
4060  OutVector & dst,
4061  const bool zero_dst_vector_setting,
4062  const Container & container,
4063  function_type cell_function,
4064  function_type face_function,
4065  function_type boundary_function,
4066  const typename MF::DataAccessOnFaces src_vector_face_access =
4067  MF::DataAccessOnFaces::none,
4068  const typename MF::DataAccessOnFaces dst_vector_face_access =
4069  MF::DataAccessOnFaces::none)
4070  : matrix_free(matrix_free)
4071  , container(const_cast<Container &>(container))
4072  , cell_function(cell_function)
4073  , face_function(face_function)
4074  , boundary_function(boundary_function)
4075  , src(src)
4076  , dst(dst)
4077  , src_data_exchanger(matrix_free,
4078  src_vector_face_access,
4079  n_components(src))
4080  , dst_data_exchanger(matrix_free,
4081  dst_vector_face_access,
4082  n_components(dst))
4083  , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4084  , zero_dst_vector_setting(zero_dst_vector_setting &&
4085  !src_and_dst_are_same)
4086  {}
4087 
4088  // Runs the cell work. If no function is given, nothing is done
4089  virtual void
4090  cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4091  {
4092  if (cell_function != nullptr && cell_range.second > cell_range.first)
4093  (container.*
4094  cell_function)(matrix_free, this->dst, this->src, cell_range);
4095  }
4096 
4097  // Runs the assembler on interior faces. If no function is given, nothing
4098  // is done
4099  virtual void
4100  face(const std::pair<unsigned int, unsigned int> &face_range) override
4101  {
4102  if (face_function != nullptr && face_range.second > face_range.first)
4103  (container.*
4104  face_function)(matrix_free, this->dst, this->src, face_range);
4105  }
4106 
4107  // Runs the assembler on boundary faces. If no function is given, nothing
4108  // is done
4109  virtual void
4110  boundary(const std::pair<unsigned int, unsigned int> &face_range) override
4111  {
4112  if (boundary_function != nullptr && face_range.second > face_range.first)
4113  (container.*
4114  boundary_function)(matrix_free, this->dst, this->src, face_range);
4115  }
4116 
4117  // Starts the communication for the update ghost values operation. We
4118  // cannot call this update if ghost and destination are the same because
4119  // that would introduce spurious entries in the destination (there is also
4120  // the problem that reading from a vector that we also write to is usually
4121  // not intended in case there is overlap, but this is up to the
4122  // application code to decide and we cannot catch this case here).
4123  virtual void
4124  vector_update_ghosts_start() override
4125  {
4126  if (!src_and_dst_are_same)
4127  internal::update_ghost_values_start(src, src_data_exchanger);
4128  }
4129 
4130  // Finishes the communication for the update ghost values operation
4131  virtual void
4132  vector_update_ghosts_finish() override
4133  {
4134  if (!src_and_dst_are_same)
4135  internal::update_ghost_values_finish(src, src_data_exchanger);
4136  }
4137 
4138  // Starts the communication for the vector compress operation
4139  virtual void
4140  vector_compress_start() override
4141  {
4142  internal::compress_start(dst, dst_data_exchanger);
4143  }
4144 
4145  // Finishes the communication for the vector compress operation
4146  virtual void
4147  vector_compress_finish() override
4148  {
4149  internal::compress_finish(dst, dst_data_exchanger);
4150  if (!src_and_dst_are_same)
4151  internal::reset_ghost_values(src, src_data_exchanger);
4152  }
4153 
4154  // Zeros the given input vector
4155  virtual void
4156  zero_dst_vector_range(const unsigned int range_index) override
4157  {
4158  if (zero_dst_vector_setting)
4159  internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4160  }
4161 
4162  private:
4163  const MF & matrix_free;
4164  Container & container;
4165  function_type cell_function;
4166  function_type face_function;
4167  function_type boundary_function;
4168 
4169  const InVector &src;
4170  OutVector & dst;
4171  VectorDataExchange<MF::dimension,
4172  typename MF::value_type,
4173  typename MF::vectorized_value_type>
4174  src_data_exchanger;
4175  VectorDataExchange<MF::dimension,
4176  typename MF::value_type,
4177  typename MF::vectorized_value_type>
4178  dst_data_exchanger;
4179  const bool src_and_dst_are_same;
4180  const bool zero_dst_vector_setting;
4181  };
4182 
4183 
4184 
4189  template <class MF, typename InVector, typename OutVector>
4190  struct MFClassWrapper
4191  {
4192  using function_type =
4193  std::function<void(const MF &,
4194  OutVector &,
4195  const InVector &,
4196  const std::pair<unsigned int, unsigned int> &)>;
4197 
4198  MFClassWrapper(const function_type cell,
4199  const function_type face,
4200  const function_type boundary)
4201  : cell(cell)
4202  , face(face)
4203  , boundary(boundary)
4204  {}
4205 
4206  void
4207  cell_integrator(const MF & mf,
4208  OutVector & dst,
4209  const InVector & src,
4210  const std::pair<unsigned int, unsigned int> &range) const
4211  {
4212  if (cell)
4213  cell(mf, dst, src, range);
4214  }
4215 
4216  void
4217  face_integrator(const MF & mf,
4218  OutVector & dst,
4219  const InVector & src,
4220  const std::pair<unsigned int, unsigned int> &range) const
4221  {
4222  if (face)
4223  face(mf, dst, src, range);
4224  }
4225 
4226  void
4227  boundary_integrator(
4228  const MF & mf,
4229  OutVector & dst,
4230  const InVector & src,
4231  const std::pair<unsigned int, unsigned int> &range) const
4232  {
4233  if (boundary)
4234  boundary(mf, dst, src, range);
4235  }
4236 
4237  const function_type cell;
4238  const function_type face;
4239  const function_type boundary;
4240  };
4241 
4242 } // end of namespace internal
4243 
4244 
4245 
4246 template <int dim, typename Number, typename VectorizedArrayType>
4247 template <typename OutVector, typename InVector>
4248 inline void
4250  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4251  OutVector &,
4252  const InVector &,
4253  const std::pair<unsigned int, unsigned int> &)>
4254  & cell_operation,
4255  OutVector & dst,
4256  const InVector &src,
4257  const bool zero_dst_vector) const
4258 {
4259  using Wrapper =
4260  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4261  InVector,
4262  OutVector>;
4263  Wrapper wrap(cell_operation, nullptr, nullptr);
4264  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4265  InVector,
4266  OutVector,
4267  Wrapper,
4268  true>
4269  worker(*this,
4270  src,
4271  dst,
4272  zero_dst_vector,
4273  wrap,
4274  &Wrapper::cell_integrator,
4275  &Wrapper::face_integrator,
4276  &Wrapper::boundary_integrator);
4277 
4278  task_info.loop(worker);
4279 }
4280 
4281 
4282 
4283 template <int dim, typename Number, typename VectorizedArrayType>
4284 template <typename OutVector, typename InVector>
4285 inline void
4287  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4288  OutVector &,
4289  const InVector &,
4290  const std::pair<unsigned int, unsigned int> &)>
4291  &cell_operation,
4292  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4293  OutVector &,
4294  const InVector &,
4295  const std::pair<unsigned int, unsigned int> &)>
4296  &face_operation,
4297  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4298  OutVector &,
4299  const InVector &,
4300  const std::pair<unsigned int, unsigned int> &)>
4301  & boundary_operation,
4302  OutVector & dst,
4303  const InVector & src,
4304  const bool zero_dst_vector,
4305  const DataAccessOnFaces dst_vector_face_access,
4306  const DataAccessOnFaces src_vector_face_access) const
4307 {
4308  using Wrapper =
4309  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4310  InVector,
4311  OutVector>;
4312  Wrapper wrap(cell_operation, face_operation, boundary_operation);
4313  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4314  InVector,
4315  OutVector,
4316  Wrapper,
4317  true>
4318  worker(*this,
4319  src,
4320  dst,
4321  zero_dst_vector,
4322  wrap,
4323  &Wrapper::cell_integrator,
4324  &Wrapper::face_integrator,
4325  &Wrapper::boundary_integrator,
4326  src_vector_face_access,
4327  dst_vector_face_access);
4328 
4329  task_info.loop(worker);
4330 }
4331 
4332 
4333 
4334 template <int dim, typename Number, typename VectorizedArrayType>
4335 template <typename CLASS, typename OutVector, typename InVector>
4336 inline void
4338  void (CLASS::*function_pointer)(
4340  OutVector &,
4341  const InVector &,
4342  const std::pair<unsigned int, unsigned int> &) const,
4343  const CLASS * owning_class,
4344  OutVector & dst,
4345  const InVector &src,
4346  const bool zero_dst_vector) const
4347 {
4348  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4349  InVector,
4350  OutVector,
4351  CLASS,
4352  true>
4353  worker(*this,
4354  src,
4355  dst,
4356  zero_dst_vector,
4357  *owning_class,
4358  function_pointer,
4359  nullptr,
4360  nullptr);
4361  task_info.loop(worker);
4362 }
4363 
4364 
4365 
4366 template <int dim, typename Number, typename VectorizedArrayType>
4367 template <typename CLASS, typename OutVector, typename InVector>
4368 inline void
4370  void (CLASS::*cell_operation)(
4372  OutVector &,
4373  const InVector &,
4374  const std::pair<unsigned int, unsigned int> &) const,
4375  void (CLASS::*face_operation)(
4377  OutVector &,
4378  const InVector &,
4379  const std::pair<unsigned int, unsigned int> &) const,
4380  void (CLASS::*boundary_operation)(
4382  OutVector &,
4383  const InVector &,
4384  const std::pair<unsigned int, unsigned int> &) const,
4385  const CLASS * owning_class,
4386  OutVector & dst,
4387  const InVector & src,
4388  const bool zero_dst_vector,
4389  const DataAccessOnFaces dst_vector_face_access,
4390  const DataAccessOnFaces src_vector_face_access) const
4391 {
4392  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4393  InVector,
4394  OutVector,
4395  CLASS,
4396  true>
4397  worker(*this,
4398  src,
4399  dst,
4400  zero_dst_vector,
4401  *owning_class,
4402  cell_operation,
4403  face_operation,
4404  boundary_operation,
4405  src_vector_face_access,
4406  dst_vector_face_access);
4407  task_info.loop(worker);
4408 }
4409 
4410 
4411 
4412 template <int dim, typename Number, typename VectorizedArrayType>
4413 template <typename CLASS, typename OutVector, typename InVector>
4414 inline void
4416  void (CLASS::*function_pointer)(
4418  OutVector &,
4419  const InVector &,
4420  const std::pair<unsigned int, unsigned int> &),
4421  CLASS * owning_class,
4422  OutVector & dst,
4423  const InVector &src,
4424  const bool zero_dst_vector) const
4425 {
4426  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4427  InVector,
4428  OutVector,
4429  CLASS,
4430  false>
4431  worker(*this,
4432  src,
4433  dst,
4434  zero_dst_vector,
4435  *owning_class,
4436  function_pointer,
4437  nullptr,
4438  nullptr);
4439  task_info.loop(worker);
4440 }
4441 
4442 
4443 
4444 template <int dim, typename Number, typename VectorizedArrayType>
4445 template <typename CLASS, typename OutVector, typename InVector>
4446 inline void
4448  void (CLASS::*cell_operation)(
4450  OutVector &,
4451  const InVector &,
4452  const std::pair<unsigned int, unsigned int> &),
4453  void (CLASS::*face_operation)(
4455  OutVector &,
4456  const InVector &,
4457  const std::pair<unsigned int, unsigned int> &),
4458  void (CLASS::*boundary_operation)(
4460  OutVector &,
4461  const InVector &,
4462  const std::pair<unsigned int, unsigned int> &),
4463  CLASS * owning_class,
4464  OutVector & dst,
4465  const InVector & src,
4466  const bool zero_dst_vector,
4467  const DataAccessOnFaces dst_vector_face_access,
4468  const DataAccessOnFaces src_vector_face_access) const
4469 {
4470  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4471  InVector,
4472  OutVector,
4473  CLASS,
4474  false>
4475  worker(*this,
4476  src,
4477  dst,
4478  zero_dst_vector,
4479  *owning_class,
4480  cell_operation,
4481  face_operation,
4482  boundary_operation,
4483  src_vector_face_access,
4484  dst_vector_face_access);
4485  task_info.loop(worker);
4486 }
4487 
4488 
4489 #endif // ifndef DOXYGEN
4490 
4491 
4492 
4493 DEAL_II_NAMESPACE_CLOSE
4494 
4495 #endif
Transformed quadrature weights.
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > > shape_info
Definition: matrix_free.h:1716
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
AlignedVector< Number > * acquire_scratch_data_non_threadsafe() const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
~MatrixFree() override=default
unsigned int n_boundary_face_batches() const
A class that provides a separate storage location on each thread that accesses the object...
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
unsigned int get_cell_category(const unsigned int macro_cell) 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:1683
bool mapping_is_initialized
Definition: matrix_free.h:1757
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition: dof_info.h:472
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:598
std::list< std::pair< bool, AlignedVector< Number > > > scratch_pad_non_threadsafe
Definition: matrix_free.h:1775
void make_connectivity_graph_faces(DynamicSparsityPattern &connectivity)
std::pair< unsigned int, unsigned int > get_face_category(const unsigned int macro_face) 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:126
UpdateFlags mapping_update_flags_boundary_faces
Definition: matrix_free.h:325
static ::ExceptionBase & ExcNotInitialized()
std::vector< unsigned int > cell_partition_data
Definition: task_info.h:465
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
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:304
std::size_t memory_consumption() const
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
unsigned int level_mg_handler
Definition: matrix_free.h:384
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
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_number) const
static const unsigned int dimension
Definition: matrix_free.h:132
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
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
unsigned int n_ghost_inner_face_batches() const
#define Assert(cond, exc)
Definition: exceptions.h:1407
UpdateFlags
const types::boundary_id invalid_boundary_id
Definition: types.h:228
std::vector< std::pair< unsigned int, unsigned int > > cell_level_index
Definition: matrix_free.h:1724
UpdateFlags mapping_update_flags_inner_faces
Definition: matrix_free.h:346
std::vector< unsigned int > boundary_partition_data
Definition: task_info.h:483
const internal::MatrixFreeFunctions::TaskInfo & get_size_info() const
std::vector< internal::MatrixFreeFunctions::DoFInfo > dof_info
Definition: matrix_free.h:1689
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:374
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 level_mg_handler=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:213
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:280
std::vector< unsigned int > cell_vectorization_category
Definition: matrix_free.h:453
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:154
std::vector< Number > constraint_pool_data
Definition: matrix_free.h:1697
types::boundary_id get_boundary_id(const unsigned int macro_face) const
std::vector< unsigned int > face_partition_data
Definition: task_info.h:474
unsigned int n_macro_cells() const
::Table< 3, types::boundary_id > cell_and_face_boundary_id
Definition: face_info.h:167
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:1740
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:1752
void renumber_dofs(std::vector< types::global_dof_index > &renumbering, const unsigned int dof_handler_index=0)
Threads::ThreadLocalStorage< std::list< std::pair< bool, AlignedVector< VectorizedArrayType > > > > scratch_pad
Definition: matrix_free.h:1768
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:1710
static ::ExceptionBase & ExcNotImplemented()
std::vector< unsigned int > vector_zero_range_list
Definition: dof_info.h:603
Definition: table.h:699
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:1703
static bool is_supported(const FiniteElement< dim, spacedim > &fe)
internal::MatrixFreeFunctions::FaceInfo< VectorizedArrayType::n_array_elements > face_info
Definition: matrix_free.h:1747
unsigned int tasks_block_size
Definition: matrix_free.h:291
const Quadrature< dim > & get_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
static const unsigned int chunk_size_zero_vector
Definition: dof_info.h:81
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:1733
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()