Reference documentation for deal.II version 9.1.1
\(\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, typename Number = double>
113 class MatrixFree : public Subscriptor
114 {
115 public:
120  using value_type = Number;
121 
125  static const unsigned int dimension = dim;
126 
173  struct AdditionalData
174  {
180  enum TasksParallelScheme
181  {
185  none = internal::MatrixFreeFunctions::TaskInfo::none,
189  partition_partition =
190  internal::MatrixFreeFunctions::TaskInfo::partition_partition,
194  partition_color =
195  internal::MatrixFreeFunctions::TaskInfo::partition_color,
200  color = internal::MatrixFreeFunctions::TaskInfo::color
201  };
202 
206  AdditionalData(
207  const TasksParallelScheme tasks_parallel_scheme = partition_partition,
208  const unsigned int tasks_block_size = 0,
209  const UpdateFlags mapping_update_flags = update_gradients |
210  update_JxW_values,
211  const UpdateFlags mapping_update_flags_boundary_faces = update_default,
212  const UpdateFlags mapping_update_flags_inner_faces = update_default,
213  const UpdateFlags mapping_update_flags_faces_by_cells = update_default,
214  const unsigned int level_mg_handler = numbers::invalid_unsigned_int,
215  const bool store_plain_indices = true,
216  const bool initialize_indices = true,
217  const bool initialize_mapping = true,
218  const bool overlap_communication_computation = true,
219  const bool hold_all_faces_to_owned_cells = false,
220  const bool cell_vectorization_categories_strict = false)
221  : tasks_parallel_scheme(tasks_parallel_scheme)
222  , tasks_block_size(tasks_block_size)
223  , mapping_update_flags(mapping_update_flags)
224  , mapping_update_flags_boundary_faces(mapping_update_flags_boundary_faces)
225  , mapping_update_flags_inner_faces(mapping_update_flags_inner_faces)
226  , mapping_update_flags_faces_by_cells(mapping_update_flags_faces_by_cells)
227  , level_mg_handler(level_mg_handler)
228  , store_plain_indices(store_plain_indices)
229  , initialize_indices(initialize_indices)
230  , initialize_mapping(initialize_mapping)
231  , overlap_communication_computation(overlap_communication_computation)
232  , hold_all_faces_to_owned_cells(hold_all_faces_to_owned_cells)
233  , cell_vectorization_categories_strict(
234  cell_vectorization_categories_strict)
235  {}
236 
273  TasksParallelScheme tasks_parallel_scheme;
274 
284  unsigned int tasks_block_size;
285 
297  UpdateFlags mapping_update_flags;
298 
318  UpdateFlags mapping_update_flags_boundary_faces;
319 
339  UpdateFlags mapping_update_flags_inner_faces;
340 
367  UpdateFlags mapping_update_flags_faces_by_cells;
368 
377  unsigned int level_mg_handler;
378 
385  bool store_plain_indices;
386 
396  bool initialize_indices;
397 
405  bool initialize_mapping;
406 
419  bool overlap_communication_computation;
420 
428  bool hold_all_faces_to_owned_cells;
429 
446  std::vector<unsigned int> cell_vectorization_category;
447 
457  bool cell_vectorization_categories_strict;
458  };
459 
467  MatrixFree();
468 
472  MatrixFree(const MatrixFree<dim, Number> &other);
473 
477  ~MatrixFree() override = default;
478 
491  template <typename DoFHandlerType, typename QuadratureType, typename number2>
492  void
493  reinit(const Mapping<dim> & mapping,
494  const DoFHandlerType & dof_handler,
495  const AffineConstraints<number2> &constraint,
496  const QuadratureType & quad,
497  const AdditionalData additional_data = AdditionalData());
498 
503  template <typename DoFHandlerType, typename QuadratureType, typename number2>
504  void
505  reinit(const DoFHandlerType & dof_handler,
506  const AffineConstraints<number2> &constraint,
507  const QuadratureType & quad,
508  const AdditionalData additional_data = AdditionalData());
509 
517  template <typename DoFHandlerType, typename QuadratureType, typename number2>
518  DEAL_II_DEPRECATED void
519  reinit(const Mapping<dim> & mapping,
520  const DoFHandlerType & dof_handler,
521  const AffineConstraints<number2> &constraint,
522  const IndexSet & locally_owned_dofs,
523  const QuadratureType & quad,
524  const AdditionalData additional_data = AdditionalData());
525 
547  template <typename DoFHandlerType, typename QuadratureType, typename number2>
548  void
549  reinit(const Mapping<dim> & mapping,
550  const std::vector<const DoFHandlerType *> & dof_handler,
551  const std::vector<const AffineConstraints<number2> *> &constraint,
552  const std::vector<QuadratureType> & quad,
553  const AdditionalData additional_data = AdditionalData());
554 
559  template <typename DoFHandlerType, typename QuadratureType, typename number2>
560  void
561  reinit(const std::vector<const DoFHandlerType *> & dof_handler,
562  const std::vector<const AffineConstraints<number2> *> &constraint,
563  const std::vector<QuadratureType> & quad,
564  const AdditionalData additional_data = AdditionalData());
565 
573  template <typename DoFHandlerType, typename QuadratureType, typename number2>
574  DEAL_II_DEPRECATED void
575  reinit(const Mapping<dim> & mapping,
576  const std::vector<const DoFHandlerType *> & dof_handler,
577  const std::vector<const AffineConstraints<number2> *> &constraint,
578  const std::vector<IndexSet> & locally_owned_set,
579  const std::vector<QuadratureType> &quad,
580  const AdditionalData additional_data = AdditionalData());
581 
589  template <typename DoFHandlerType, typename QuadratureType, typename number2>
590  void
591  reinit(const Mapping<dim> & mapping,
592  const std::vector<const DoFHandlerType *> & dof_handler,
593  const std::vector<const AffineConstraints<number2> *> &constraint,
594  const QuadratureType & quad,
595  const AdditionalData additional_data = AdditionalData());
596 
601  template <typename DoFHandlerType, typename QuadratureType, typename number2>
602  void
603  reinit(const std::vector<const DoFHandlerType *> & dof_handler,
604  const std::vector<const AffineConstraints<number2> *> &constraint,
605  const QuadratureType & quad,
606  const AdditionalData additional_data = AdditionalData());
607 
613  void
614  copy_from(const MatrixFree<dim, Number> &matrix_free_base);
615 
620  void
621  clear();
622 
624 
636  enum class DataAccessOnFaces
637  {
644  none,
645 
655  values,
656 
667  gradients,
668 
675  unspecified
676  };
677 
718  template <typename OutVector, typename InVector>
719  void
720  cell_loop(
721  const std::function<void(const MatrixFree<dim, Number> &,
722  OutVector &,
723  const InVector &,
724  const std::pair<unsigned int, unsigned int> &)>
725  & cell_operation,
726  OutVector & dst,
727  const InVector &src,
728  const bool zero_dst_vector = false) const;
729 
771  template <typename CLASS, typename OutVector, typename InVector>
772  void
773  cell_loop(void (CLASS::*cell_operation)(
774  const MatrixFree &,
775  OutVector &,
776  const InVector &,
777  const std::pair<unsigned int, unsigned int> &) const,
778  const CLASS * owning_class,
779  OutVector & dst,
780  const InVector &src,
781  const bool zero_dst_vector = false) const;
782 
786  template <typename CLASS, typename OutVector, typename InVector>
787  void
788  cell_loop(void (CLASS::*cell_operation)(
789  const MatrixFree &,
790  OutVector &,
791  const InVector &,
792  const std::pair<unsigned int, unsigned int> &),
793  CLASS * owning_class,
794  OutVector & dst,
795  const InVector &src,
796  const bool zero_dst_vector = false) const;
797 
873  template <typename OutVector, typename InVector>
874  void
875  loop(const std::function<void(const MatrixFree<dim, Number> &,
876  OutVector &,
877  const InVector &,
878  const std::pair<unsigned int, unsigned int> &)>
879  &cell_operation,
880  const std::function<void(const MatrixFree<dim, Number> &,
881  OutVector &,
882  const InVector &,
883  const std::pair<unsigned int, unsigned int> &)>
884  &face_operation,
885  const std::function<void(const MatrixFree<dim, Number> &,
886  OutVector &,
887  const InVector &,
888  const std::pair<unsigned int, unsigned int> &)>
889  & boundary_operation,
890  OutVector & dst,
891  const InVector & src,
892  const bool zero_dst_vector = false,
893  const DataAccessOnFaces dst_vector_face_access =
894  DataAccessOnFaces::unspecified,
895  const DataAccessOnFaces src_vector_face_access =
896  DataAccessOnFaces::unspecified) const;
897 
982  template <typename CLASS, typename OutVector, typename InVector>
983  void
984  loop(
985  void (CLASS::*cell_operation)(const MatrixFree &,
986  OutVector &,
987  const InVector &,
988  const std::pair<unsigned int, unsigned int> &)
989  const,
990  void (CLASS::*face_operation)(const MatrixFree &,
991  OutVector &,
992  const InVector &,
993  const std::pair<unsigned int, unsigned int> &)
994  const,
995  void (CLASS::*boundary_operation)(
996  const MatrixFree &,
997  OutVector &,
998  const InVector &,
999  const std::pair<unsigned int, unsigned int> &) const,
1000  const CLASS * owning_class,
1001  OutVector & dst,
1002  const InVector & src,
1003  const bool zero_dst_vector = false,
1004  const DataAccessOnFaces dst_vector_face_access =
1005  DataAccessOnFaces::unspecified,
1006  const DataAccessOnFaces src_vector_face_access =
1007  DataAccessOnFaces::unspecified) const;
1008 
1012  template <typename CLASS, typename OutVector, typename InVector>
1013  void
1014  loop(void (CLASS::*cell_operation)(
1015  const MatrixFree &,
1016  OutVector &,
1017  const InVector &,
1018  const std::pair<unsigned int, unsigned int> &),
1019  void (CLASS::*face_operation)(
1020  const MatrixFree &,
1021  OutVector &,
1022  const InVector &,
1023  const std::pair<unsigned int, unsigned int> &),
1024  void (CLASS::*boundary_operation)(
1025  const MatrixFree &,
1026  OutVector &,
1027  const InVector &,
1028  const std::pair<unsigned int, unsigned int> &),
1029  CLASS * owning_class,
1030  OutVector & dst,
1031  const InVector & src,
1032  const bool zero_dst_vector = false,
1033  const DataAccessOnFaces dst_vector_face_access =
1034  DataAccessOnFaces::unspecified,
1035  const DataAccessOnFaces src_vector_face_access =
1036  DataAccessOnFaces::unspecified) const;
1037 
1045  std::pair<unsigned int, unsigned int>
1046  create_cell_subrange_hp(const std::pair<unsigned int, unsigned int> &range,
1047  const unsigned int fe_degree,
1048  const unsigned int dof_handler_index = 0) const;
1049 
1056  std::pair<unsigned int, unsigned int>
1057  create_cell_subrange_hp_by_index(
1058  const std::pair<unsigned int, unsigned int> &range,
1059  const unsigned int fe_index,
1060  const unsigned int dof_handler_index = 0) const;
1061 
1063 
1088  template <typename VectorType>
1089  void
1090  initialize_dof_vector(VectorType & vec,
1091  const unsigned int dof_handler_index = 0) const;
1092 
1113  template <typename Number2>
1114  void
1115  initialize_dof_vector(LinearAlgebra::distributed::Vector<Number2> &vec,
1116  const unsigned int dof_handler_index = 0) const;
1117 
1128  const std::shared_ptr<const Utilities::MPI::Partitioner> &
1129  get_vector_partitioner(const unsigned int dof_handler_index = 0) const;
1130 
1134  const IndexSet &
1135  get_locally_owned_set(const unsigned int dof_handler_index = 0) const;
1136 
1140  const IndexSet &
1141  get_ghost_set(const unsigned int dof_handler_index = 0) const;
1142 
1152  const std::vector<unsigned int> &
1153  get_constrained_dofs(const unsigned int dof_handler_index = 0) const;
1154 
1165  void
1166  renumber_dofs(std::vector<types::global_dof_index> &renumbering,
1167  const unsigned int dof_handler_index = 0);
1168 
1170 
1178  template <int spacedim>
1179  static bool
1180  is_supported(const FiniteElement<dim, spacedim> &fe);
1181 
1185  unsigned int
1186  n_components() const;
1187 
1192  unsigned int
1193  n_base_elements(const unsigned int dof_handler_index) const;
1194 
1202  unsigned int
1203  n_physical_cells() const;
1204 
1214  unsigned int
1215  n_macro_cells() const;
1216 
1226  unsigned int
1227  n_cell_batches() const;
1228 
1235  unsigned int
1236  n_ghost_cell_batches() const;
1237 
1245  unsigned int
1246  n_inner_face_batches() const;
1247 
1256  unsigned int
1257  n_boundary_face_batches() const;
1258 
1263  unsigned int
1264  n_ghost_inner_face_batches() const;
1265 
1272  types::boundary_id
1273  get_boundary_id(const unsigned int macro_face) const;
1274 
1279  std::array<types::boundary_id, VectorizedArray<Number>::n_array_elements>
1280  get_faces_by_cells_boundary_id(const unsigned int macro_cell,
1281  const unsigned int face_number) const;
1282 
1287  const DoFHandler<dim> &
1288  get_dof_handler(const unsigned int dof_handler_index = 0) const;
1289 
1300  typename DoFHandler<dim>::cell_iterator
1301  get_cell_iterator(const unsigned int macro_cell_number,
1302  const unsigned int vector_number,
1303  const unsigned int fe_component = 0) const;
1304 
1316  typename hp::DoFHandler<dim>::active_cell_iterator
1317  get_hp_cell_iterator(const unsigned int macro_cell_number,
1318  const unsigned int vector_number,
1319  const unsigned int dof_handler_index = 0) const;
1320 
1332  bool
1333  at_irregular_cell(const unsigned int macro_cell_number) const;
1334 
1343  unsigned int
1344  n_components_filled(const unsigned int cell_batch_number) const;
1345 
1354  unsigned int
1355  n_active_entries_per_cell_batch(const unsigned int cell_batch_number) const;
1356 
1366  unsigned int
1367  n_active_entries_per_face_batch(const unsigned int face_batch_number) const;
1368 
1372  unsigned int
1373  get_dofs_per_cell(const unsigned int dof_handler_index = 0,
1374  const unsigned int hp_active_fe_index = 0) const;
1375 
1379  unsigned int
1380  get_n_q_points(const unsigned int quad_index = 0,
1381  const unsigned int hp_active_fe_index = 0) const;
1382 
1387  unsigned int
1388  get_dofs_per_face(const unsigned int fe_component = 0,
1389  const unsigned int hp_active_fe_index = 0) const;
1390 
1395  unsigned int
1396  get_n_q_points_face(const unsigned int quad_index = 0,
1397  const unsigned int hp_active_fe_index = 0) const;
1398 
1402  const Quadrature<dim> &
1403  get_quadrature(const unsigned int quad_index = 0,
1404  const unsigned int hp_active_fe_index = 0) const;
1405 
1409  const Quadrature<dim - 1> &
1410  get_face_quadrature(const unsigned int quad_index = 0,
1411  const unsigned int hp_active_fe_index = 0) const;
1412 
1419  unsigned int
1420  get_cell_category(const unsigned int macro_cell) const;
1421 
1426  std::pair<unsigned int, unsigned int>
1427  get_face_category(const unsigned int macro_face) const;
1428 
1432  bool
1433  indices_initialized() const;
1434 
1439  bool
1440  mapping_initialized() const;
1441 
1446  std::size_t
1447  memory_consumption() const;
1448 
1453  template <typename StreamType>
1454  void
1455  print_memory_consumption(StreamType &out) const;
1456 
1461  void
1462  print(std::ostream &out) const;
1463 
1465 
1475  const internal::MatrixFreeFunctions::TaskInfo &
1476  get_task_info() const;
1477 
1481  DEAL_II_DEPRECATED
1482  const internal::MatrixFreeFunctions::TaskInfo &
1483  get_size_info() const;
1484 
1485  /*
1486  * Return geometry-dependent information on the cells.
1487  */
1488  const internal::MatrixFreeFunctions::MappingInfo<dim, Number> &
1489  get_mapping_info() const;
1490 
1494  const internal::MatrixFreeFunctions::DoFInfo &
1495  get_dof_info(const unsigned int dof_handler_index_component = 0) const;
1496 
1500  unsigned int
1501  n_constraint_pool_entries() const;
1502 
1507  const Number *
1508  constraint_pool_begin(const unsigned int pool_index) const;
1509 
1515  const Number *
1516  constraint_pool_end(const unsigned int pool_index) const;
1517 
1521  const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &
1522  get_shape_info(const unsigned int dof_handler_index_component = 0,
1523  const unsigned int quad_index = 0,
1524  const unsigned int fe_base_element = 0,
1525  const unsigned int hp_active_fe_index = 0,
1526  const unsigned int hp_active_quad_index = 0) const;
1527 
1531  const internal::MatrixFreeFunctions::FaceToCellTopology<
1532  VectorizedArray<Number>::n_array_elements> &
1533  get_face_info(const unsigned int face_batch_number) const;
1534 
1548  AlignedVector<VectorizedArray<Number>> *
1549  acquire_scratch_data() const;
1550 
1554  void
1555  release_scratch_data(
1556  const AlignedVector<VectorizedArray<Number>> *memory) const;
1557 
1567  AlignedVector<Number> *
1568  acquire_scratch_data_non_threadsafe() const;
1569 
1573  void
1574  release_scratch_data_non_threadsafe(
1575  const AlignedVector<Number> *memory) const;
1576 
1578 
1579 private:
1584  template <typename number2>
1585  void
1586  internal_reinit(
1587  const Mapping<dim> & mapping,
1588  const std::vector<const DoFHandler<dim> *> & dof_handler,
1589  const std::vector<const AffineConstraints<number2> *> &constraint,
1590  const std::vector<IndexSet> & locally_owned_set,
1591  const std::vector<hp::QCollection<1>> & quad,
1592  const AdditionalData & additional_data);
1593 
1597  template <typename number2>
1598  void
1599  internal_reinit(
1600  const Mapping<dim> & mapping,
1601  const std::vector<const hp::DoFHandler<dim> *> & dof_handler,
1602  const std::vector<const AffineConstraints<number2> *> &constraint,
1603  const std::vector<IndexSet> & locally_owned_set,
1604  const std::vector<hp::QCollection<1>> & quad,
1605  const AdditionalData & additional_data);
1606 
1613  template <typename number2>
1614  void
1615  initialize_indices(
1616  const std::vector<const AffineConstraints<number2> *> &constraint,
1617  const std::vector<IndexSet> & locally_owned_set,
1618  const AdditionalData & additional_data);
1619 
1623  void
1624  initialize_dof_handlers(
1625  const std::vector<const DoFHandler<dim> *> &dof_handlers,
1626  const AdditionalData & additional_data);
1627 
1631  void
1632  initialize_dof_handlers(
1633  const std::vector<const hp::DoFHandler<dim> *> &dof_handlers,
1634  const AdditionalData & additional_data);
1635 
1640  void
1641  make_connectivity_graph_faces(DynamicSparsityPattern &connectivity);
1642 
1648  struct DoFHandlers
1649  {
1650  DoFHandlers()
1651  : active_dof_handler(usual)
1652  , n_dof_handlers(0)
1653  , level(numbers::invalid_unsigned_int)
1654  {}
1655 
1656  std::vector<SmartPointer<const DoFHandler<dim>>> dof_handler;
1657  std::vector<SmartPointer<const hp::DoFHandler<dim>>> hp_dof_handler;
1658  enum ActiveDoFHandler
1659  {
1663  usual,
1667  hp
1668  } active_dof_handler;
1669  unsigned int n_dof_handlers;
1670  unsigned int level;
1671  };
1672 
1676  DoFHandlers dof_handlers;
1677 
1682  std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
1683 
1690  std::vector<Number> constraint_pool_data;
1691 
1696  std::vector<unsigned int> constraint_pool_row_index;
1697 
1702  internal::MatrixFreeFunctions::MappingInfo<dim, Number> mapping_info;
1703 
1707  Table<4, internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>>>
1708  shape_info;
1709 
1716  std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
1717 
1718 
1725  unsigned int cell_level_index_end_local;
1726 
1732  internal::MatrixFreeFunctions::TaskInfo task_info;
1733 
1738  internal::MatrixFreeFunctions::FaceInfo<
1739  VectorizedArray<Number>::n_array_elements>
1740  face_info;
1741 
1745  bool indices_are_initialized;
1746 
1750  bool mapping_is_initialized;
1751 
1759  mutable Threads::ThreadLocalStorage<
1760  std::list<std::pair<bool, AlignedVector<VectorizedArray<Number>>>>>
1761  scratch_pad;
1762 
1767  mutable std::list<std::pair<bool, AlignedVector<Number>>>
1768  scratch_pad_non_threadsafe;
1769 };
1770 
1771 
1772 
1773 /*----------------------- Inline functions ----------------------------------*/
1774 
1775 #ifndef DOXYGEN
1776 
1777 
1778 
1779 template <int dim, typename Number>
1780 template <typename VectorType>
1781 inline void
1782 MatrixFree<dim, Number>::initialize_dof_vector(VectorType & vec,
1783  const unsigned int comp) const
1784 {
1785  AssertIndexRange(comp, n_components());
1786  vec.reinit(dof_info[comp].vector_partitioner->size());
1787 }
1788 
1789 
1790 
1791 template <int dim, typename Number>
1792 template <typename Number2>
1793 inline void
1794 MatrixFree<dim, Number>::initialize_dof_vector(
1795  LinearAlgebra::distributed::Vector<Number2> &vec,
1796  const unsigned int comp) const
1797 {
1798  AssertIndexRange(comp, n_components());
1799  vec.reinit(dof_info[comp].vector_partitioner);
1800 }
1801 
1802 
1803 
1804 template <int dim, typename Number>
1805 inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
1806 MatrixFree<dim, Number>::get_vector_partitioner(const unsigned int comp) const
1807 {
1808  AssertIndexRange(comp, n_components());
1809  return dof_info[comp].vector_partitioner;
1810 }
1811 
1812 
1813 
1814 template <int dim, typename Number>
1815 inline const std::vector<unsigned int> &
1816 MatrixFree<dim, Number>::get_constrained_dofs(const unsigned int comp) const
1817 {
1818  AssertIndexRange(comp, n_components());
1819  return dof_info[comp].constrained_dofs;
1820 }
1821 
1822 
1823 
1824 template <int dim, typename Number>
1825 inline unsigned int
1826 MatrixFree<dim, Number>::n_components() const
1827 {
1828  AssertDimension(dof_handlers.n_dof_handlers, dof_info.size());
1829  return dof_handlers.n_dof_handlers;
1830 }
1831 
1832 
1833 
1834 template <int dim, typename Number>
1835 inline unsigned int
1836 MatrixFree<dim, Number>::n_base_elements(const unsigned int dof_no) const
1837 {
1838  AssertDimension(dof_handlers.n_dof_handlers, dof_info.size());
1839  AssertIndexRange(dof_no, dof_handlers.n_dof_handlers);
1840  return dof_handlers.dof_handler[dof_no]->get_fe().n_base_elements();
1841 }
1842 
1843 
1844 
1845 template <int dim, typename Number>
1846 inline const internal::MatrixFreeFunctions::TaskInfo &
1847 MatrixFree<dim, Number>::get_task_info() const
1848 {
1849  return task_info;
1850 }
1851 
1852 
1853 
1854 template <int dim, typename Number>
1855 inline const internal::MatrixFreeFunctions::TaskInfo &
1856 MatrixFree<dim, Number>::get_size_info() const
1857 {
1858  return task_info;
1859 }
1860 
1861 
1862 
1863 template <int dim, typename Number>
1864 inline unsigned int
1865 MatrixFree<dim, Number>::n_macro_cells() const
1866 {
1867  return *(task_info.cell_partition_data.end() - 2);
1868 }
1869 
1870 
1871 
1872 template <int dim, typename Number>
1873 inline unsigned int
1874 MatrixFree<dim, Number>::n_physical_cells() const
1875 {
1876  return task_info.n_active_cells;
1877 }
1878 
1879 
1880 
1881 template <int dim, typename Number>
1882 inline unsigned int
1883 MatrixFree<dim, Number>::n_cell_batches() const
1884 {
1885  return *(task_info.cell_partition_data.end() - 2);
1886 }
1887 
1888 
1889 
1890 template <int dim, typename Number>
1891 inline unsigned int
1892 MatrixFree<dim, Number>::n_ghost_cell_batches() const
1893 {
1894  return *(task_info.cell_partition_data.end() - 1) -
1895  *(task_info.cell_partition_data.end() - 2);
1896 }
1897 
1898 
1899 
1900 template <int dim, typename Number>
1901 inline unsigned int
1902 MatrixFree<dim, Number>::n_inner_face_batches() const
1903 {
1904  if (task_info.face_partition_data.size() == 0)
1905  return 0;
1906  return task_info.face_partition_data.back();
1907 }
1908 
1909 
1910 
1911 template <int dim, typename Number>
1912 inline unsigned int
1913 MatrixFree<dim, Number>::n_boundary_face_batches() const
1914 {
1915  if (task_info.face_partition_data.size() == 0)
1916  return 0;
1917  return task_info.boundary_partition_data.back() -
1918  task_info.face_partition_data.back();
1919 }
1920 
1921 
1922 
1923 template <int dim, typename Number>
1924 inline unsigned int
1925 MatrixFree<dim, Number>::n_ghost_inner_face_batches() const
1926 {
1927  if (task_info.face_partition_data.size() == 0)
1928  return 0;
1929  return face_info.faces.size() - task_info.boundary_partition_data.back();
1930 }
1931 
1932 
1933 
1934 template <int dim, typename Number>
1935 inline types::boundary_id
1936 MatrixFree<dim, Number>::get_boundary_id(const unsigned int macro_face) const
1937 {
1938  Assert(macro_face >= task_info.boundary_partition_data[0] &&
1939  macro_face < task_info.boundary_partition_data.back(),
1940  ExcIndexRange(macro_face,
1941  task_info.boundary_partition_data[0],
1942  task_info.boundary_partition_data.back()));
1943  return types::boundary_id(face_info.faces[macro_face].exterior_face_no);
1944 }
1945 
1946 
1947 
1948 template <int dim, typename Number>
1949 inline std::array<types::boundary_id, VectorizedArray<Number>::n_array_elements>
1950 MatrixFree<dim, Number>::get_faces_by_cells_boundary_id(
1951  const unsigned int macro_cell,
1952  const unsigned int face_number) const
1953 {
1954  AssertIndexRange(macro_cell, n_macro_cells());
1955  AssertIndexRange(face_number, GeometryInfo<dim>::faces_per_cell);
1956  Assert(face_info.cell_and_face_boundary_id.size(0) >= n_macro_cells(),
1957  ExcNotInitialized());
1958  std::array<types::boundary_id, VectorizedArray<Number>::n_array_elements>
1959  result;
1960  result.fill(numbers::invalid_boundary_id);
1961  for (unsigned int v = 0; v < n_active_entries_per_cell_batch(macro_cell); ++v)
1962  result[v] = face_info.cell_and_face_boundary_id(macro_cell, face_number, v);
1963  return result;
1964 }
1965 
1966 
1967 
1968 template <int dim, typename Number>
1969 inline const internal::MatrixFreeFunctions::MappingInfo<dim, Number> &
1970 MatrixFree<dim, Number>::get_mapping_info() const
1971 {
1972  return mapping_info;
1973 }
1974 
1975 
1976 
1977 template <int dim, typename Number>
1978 inline const internal::MatrixFreeFunctions::DoFInfo &
1979 MatrixFree<dim, Number>::get_dof_info(const unsigned int dof_index) const
1980 {
1981  AssertIndexRange(dof_index, n_components());
1982  return dof_info[dof_index];
1983 }
1984 
1985 
1986 
1987 template <int dim, typename Number>
1988 inline unsigned int
1989 MatrixFree<dim, Number>::n_constraint_pool_entries() const
1990 {
1991  return constraint_pool_row_index.size() - 1;
1992 }
1993 
1994 
1995 
1996 template <int dim, typename Number>
1997 inline const Number *
1998 MatrixFree<dim, Number>::constraint_pool_begin(const unsigned int row) const
1999 {
2000  AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2001  return constraint_pool_data.empty() ?
2002  nullptr :
2003  constraint_pool_data.data() + constraint_pool_row_index[row];
2004 }
2005 
2006 
2007 
2008 template <int dim, typename Number>
2009 inline const Number *
2010 MatrixFree<dim, Number>::constraint_pool_end(const unsigned int row) const
2011 {
2012  AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2013  return constraint_pool_data.empty() ?
2014  nullptr :
2015  constraint_pool_data.data() + constraint_pool_row_index[row + 1];
2016 }
2017 
2018 
2019 
2020 template <int dim, typename Number>
2021 inline std::pair<unsigned int, unsigned int>
2022 MatrixFree<dim, Number>::create_cell_subrange_hp(
2023  const std::pair<unsigned int, unsigned int> &range,
2024  const unsigned int degree,
2025  const unsigned int dof_handler_component) const
2026 {
2027  if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2028  {
2029  AssertDimension(
2030  dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2031  AssertDimension(
2032  dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2033  if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2034  return range;
2035  else
2036  return {range.second, range.second};
2037  }
2038 
2039  const unsigned int fe_index =
2040  dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2041  if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2042  return {range.second, range.second};
2043  else
2044  return create_cell_subrange_hp_by_index(range,
2045  fe_index,
2046  dof_handler_component);
2047 }
2048 
2049 
2050 
2051 template <int dim, typename Number>
2052 inline bool
2053 MatrixFree<dim, Number>::at_irregular_cell(const unsigned int macro_cell) const
2054 {
2055  AssertIndexRange(macro_cell, task_info.cell_partition_data.back());
2056  return VectorizedArray<Number>::n_array_elements > 1 &&
2057  cell_level_index[(macro_cell + 1) *
2058  VectorizedArray<Number>::n_array_elements -
2059  1] ==
2060  cell_level_index[(macro_cell + 1) *
2061  VectorizedArray<Number>::n_array_elements -
2062  2];
2063 }
2064 
2065 
2066 
2067 template <int dim, typename Number>
2068 inline unsigned int
2069 MatrixFree<dim, Number>::n_components_filled(
2070  const unsigned int cell_batch_number) const
2071 {
2072  return n_active_entries_per_cell_batch(cell_batch_number);
2073 }
2074 
2075 
2076 
2077 template <int dim, typename Number>
2078 inline unsigned int
2079 MatrixFree<dim, Number>::n_active_entries_per_cell_batch(
2080  const unsigned int cell_batch_number) const
2081 {
2082  AssertIndexRange(cell_batch_number, task_info.cell_partition_data.back());
2083  unsigned int n_components = VectorizedArray<Number>::n_array_elements;
2084  while (n_components > 1 &&
2085  cell_level_index[cell_batch_number *
2086  VectorizedArray<Number>::n_array_elements +
2087  n_components - 1] ==
2088  cell_level_index[cell_batch_number *
2089  VectorizedArray<Number>::n_array_elements +
2090  n_components - 2])
2091  --n_components;
2092  AssertIndexRange(n_components - 1, VectorizedArray<Number>::n_array_elements);
2093  return n_components;
2094 }
2095 
2096 
2097 
2098 template <int dim, typename Number>
2099 inline unsigned int
2100 MatrixFree<dim, Number>::n_active_entries_per_face_batch(
2101  const unsigned int face_batch_number) const
2102 {
2103  AssertIndexRange(face_batch_number, face_info.faces.size());
2104  unsigned int n_components = VectorizedArray<Number>::n_array_elements;
2105  while (n_components > 1 &&
2106  face_info.faces[face_batch_number].cells_interior[n_components - 1] ==
2107  numbers::invalid_unsigned_int)
2108  --n_components;
2109  AssertIndexRange(n_components - 1, VectorizedArray<Number>::n_array_elements);
2110  return n_components;
2111 }
2112 
2113 
2114 
2115 template <int dim, typename Number>
2116 inline unsigned int
2117 MatrixFree<dim, Number>::get_dofs_per_cell(
2118  const unsigned int dof_handler_index,
2119  const unsigned int active_fe_index) const
2120 {
2121  return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2122 }
2123 
2124 
2125 
2126 template <int dim, typename Number>
2127 inline unsigned int
2128 MatrixFree<dim, Number>::get_n_q_points(
2129  const unsigned int quad_index,
2130  const unsigned int active_fe_index) const
2131 {
2132  AssertIndexRange(quad_index, mapping_info.cell_data.size());
2133  return mapping_info.cell_data[quad_index]
2134  .descriptor[active_fe_index]
2135  .n_q_points;
2136 }
2137 
2138 
2139 
2140 template <int dim, typename Number>
2141 inline unsigned int
2142 MatrixFree<dim, Number>::get_dofs_per_face(
2143  const unsigned int dof_handler_index,
2144  const unsigned int active_fe_index) const
2145 {
2146  return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2147 }
2148 
2149 
2150 
2151 template <int dim, typename Number>
2152 inline unsigned int
2153 MatrixFree<dim, Number>::get_n_q_points_face(
2154  const unsigned int quad_index,
2155  const unsigned int active_fe_index) const
2156 {
2157  AssertIndexRange(quad_index, mapping_info.face_data.size());
2158  return mapping_info.face_data[quad_index]
2159  .descriptor[active_fe_index]
2160  .n_q_points;
2161 }
2162 
2163 
2164 
2165 template <int dim, typename Number>
2166 inline const IndexSet &
2167 MatrixFree<dim, Number>::get_locally_owned_set(
2168  const unsigned int dof_handler_index) const
2169 {
2170  return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2171 }
2172 
2173 
2174 
2175 template <int dim, typename Number>
2176 inline const IndexSet &
2177 MatrixFree<dim, Number>::get_ghost_set(
2178  const unsigned int dof_handler_index) const
2179 {
2180  return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2181 }
2182 
2183 
2184 
2185 template <int dim, typename Number>
2186 inline const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &
2187 MatrixFree<dim, Number>::get_shape_info(
2188  const unsigned int dof_handler_index,
2189  const unsigned int index_quad,
2190  const unsigned int index_fe,
2191  const unsigned int active_fe_index,
2192  const unsigned int active_quad_index) const
2193 {
2194  AssertIndexRange(dof_handler_index, dof_info.size());
2195  const unsigned int ind =
2196  dof_info[dof_handler_index].global_base_element_offset + index_fe;
2197  AssertIndexRange(ind, shape_info.size(0));
2198  AssertIndexRange(index_quad, shape_info.size(1));
2199  AssertIndexRange(active_fe_index, shape_info.size(2));
2200  AssertIndexRange(active_quad_index, shape_info.size(3));
2201  return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2202 }
2203 
2204 
2205 
2206 template <int dim, typename Number>
2207 inline const internal::MatrixFreeFunctions::FaceToCellTopology<
2208  VectorizedArray<Number>::n_array_elements> &
2209 MatrixFree<dim, Number>::get_face_info(const unsigned int macro_face) const
2210 {
2211  AssertIndexRange(macro_face, face_info.faces.size());
2212  return face_info.faces[macro_face];
2213 }
2214 
2215 
2216 
2217 template <int dim, typename Number>
2218 inline const Quadrature<dim> &
2219 MatrixFree<dim, Number>::get_quadrature(
2220  const unsigned int quad_index,
2221  const unsigned int active_fe_index) const
2222 {
2223  AssertIndexRange(quad_index, mapping_info.cell_data.size());
2224  return mapping_info.cell_data[quad_index]
2225  .descriptor[active_fe_index]
2226  .quadrature;
2227 }
2228 
2229 
2230 
2231 template <int dim, typename Number>
2232 inline const Quadrature<dim - 1> &
2233 MatrixFree<dim, Number>::get_face_quadrature(
2234  const unsigned int quad_index,
2235  const unsigned int active_fe_index) const
2236 {
2237  AssertIndexRange(quad_index, mapping_info.face_data.size());
2238  return mapping_info.face_data[quad_index]
2239  .descriptor[active_fe_index]
2240  .quadrature;
2241 }
2242 
2243 
2244 
2245 template <int dim, typename Number>
2246 inline unsigned int
2247 MatrixFree<dim, Number>::get_cell_category(const unsigned int macro_cell) const
2248 {
2249  AssertIndexRange(0, dof_info.size());
2250  AssertIndexRange(macro_cell, dof_info[0].cell_active_fe_index.size());
2251  if (dof_info[0].cell_active_fe_index.empty())
2252  return 0;
2253  else
2254  return dof_info[0].cell_active_fe_index[macro_cell];
2255 }
2256 
2257 
2258 
2259 template <int dim, typename Number>
2260 inline std::pair<unsigned int, unsigned int>
2261 MatrixFree<dim, Number>::get_face_category(const unsigned int macro_face) const
2262 {
2263  AssertIndexRange(macro_face, face_info.faces.size());
2264  if (dof_info[0].cell_active_fe_index.empty())
2265  return std::make_pair(0U, 0U);
2266 
2267  std::pair<unsigned int, unsigned int> result;
2268  for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements &&
2269  face_info.faces[macro_face].cells_interior[v] !=
2270  numbers::invalid_unsigned_int;
2271  ++v)
2272  result.first = std::max(
2273  result.first,
2274  dof_info[0]
2275  .cell_active_fe_index[face_info.faces[macro_face].cells_interior[v]]);
2276  if (face_info.faces[macro_face].cells_exterior[0] !=
2277  numbers::invalid_unsigned_int)
2278  for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements &&
2279  face_info.faces[macro_face].cells_exterior[v] !=
2280  numbers::invalid_unsigned_int;
2281  ++v)
2282  result.second = std::max(
2283  result.first,
2284  dof_info[0]
2285  .cell_active_fe_index[face_info.faces[macro_face].cells_exterior[v]]);
2286  else
2287  result.second = numbers::invalid_unsigned_int;
2288  return result;
2289 }
2290 
2291 
2292 
2293 template <int dim, typename Number>
2294 inline bool
2295 MatrixFree<dim, Number>::indices_initialized() const
2296 {
2297  return indices_are_initialized;
2298 }
2299 
2300 
2301 
2302 template <int dim, typename Number>
2303 inline bool
2304 MatrixFree<dim, Number>::mapping_initialized() const
2305 {
2306  return mapping_is_initialized;
2307 }
2308 
2309 
2310 
2311 template <int dim, typename Number>
2312 AlignedVector<VectorizedArray<Number>> *
2313 MatrixFree<dim, Number>::acquire_scratch_data() const
2314 {
2315  using list_type =
2316  std::list<std::pair<bool, AlignedVector<VectorizedArray<Number>>>>;
2317  list_type &data = scratch_pad.get();
2318  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
2319  if (it->first == false)
2320  {
2321  it->first = true;
2322  return &it->second;
2323  }
2324  data.push_front(
2325  std::make_pair(true, AlignedVector<VectorizedArray<Number>>()));
2326  return &data.front().second;
2327 }
2328 
2329 
2330 
2331 template <int dim, typename Number>
2332 void
2333 MatrixFree<dim, Number>::release_scratch_data(
2334  const AlignedVector<VectorizedArray<Number>> *scratch) const
2335 {
2336  using list_type =
2337  std::list<std::pair<bool, AlignedVector<VectorizedArray<Number>>>>;
2338  list_type &data = scratch_pad.get();
2339  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
2340  if (&it->second == scratch)
2341  {
2342  Assert(it->first == true, ExcInternalError());
2343  it->first = false;
2344  return;
2345  }
2346  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2347 }
2348 
2349 
2350 
2351 template <int dim, typename Number>
2352 AlignedVector<Number> *
2353 MatrixFree<dim, Number>::acquire_scratch_data_non_threadsafe() const
2354 {
2355  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
2356  scratch_pad_non_threadsafe.begin();
2357  it != scratch_pad_non_threadsafe.end();
2358  ++it)
2359  if (it->first == false)
2360  {
2361  it->first = true;
2362  return &it->second;
2363  }
2364  scratch_pad_non_threadsafe.push_front(
2365  std::make_pair(true, AlignedVector<Number>()));
2366  return &scratch_pad_non_threadsafe.front().second;
2367 }
2368 
2369 
2370 
2371 template <int dim, typename Number>
2372 void
2373 MatrixFree<dim, Number>::release_scratch_data_non_threadsafe(
2374  const AlignedVector<Number> *scratch) const
2375 {
2376  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
2377  scratch_pad_non_threadsafe.begin();
2378  it != scratch_pad_non_threadsafe.end();
2379  ++it)
2380  if (&it->second == scratch)
2381  {
2382  Assert(it->first == true, ExcInternalError());
2383  it->first = false;
2384  return;
2385  }
2386  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2387 }
2388 
2389 
2390 
2391 // ------------------------------ reinit functions ---------------------------
2392 
2393 namespace internal
2394 {
2395  namespace MatrixFreeImplementation
2396  {
2397  template <typename DoFHandlerType>
2398  inline std::vector<IndexSet>
2399  extract_locally_owned_index_sets(
2400  const std::vector<const DoFHandlerType *> &dofh,
2401  const unsigned int level)
2402  {
2403  std::vector<IndexSet> locally_owned_set;
2404  locally_owned_set.reserve(dofh.size());
2405  for (unsigned int j = 0; j < dofh.size(); j++)
2406  if (level == numbers::invalid_unsigned_int)
2407  locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
2408  else
2409  AssertThrow(false, ExcNotImplemented());
2410  return locally_owned_set;
2411  }
2412 
2413  template <int dim, int spacedim>
2414  inline std::vector<IndexSet>
2415  extract_locally_owned_index_sets(
2416  const std::vector<const ::DoFHandler<dim, spacedim> *> &dofh,
2417  const unsigned int level)
2418  {
2419  std::vector<IndexSet> locally_owned_set;
2420  locally_owned_set.reserve(dofh.size());
2421  for (unsigned int j = 0; j < dofh.size(); j++)
2422  if (level == numbers::invalid_unsigned_int)
2423  locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
2424  else
2425  locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
2426  return locally_owned_set;
2427  }
2428  } // namespace MatrixFreeImplementation
2429 } // namespace internal
2430 
2431 
2432 
2433 template <int dim, typename Number>
2434 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2435 void
2436 MatrixFree<dim, Number>::reinit(
2437  const DoFHandlerType & dof_handler,
2438  const AffineConstraints<number2> & constraints_in,
2439  const QuadratureType & quad,
2440  const typename MatrixFree<dim, Number>::AdditionalData additional_data)
2441 {
2442  std::vector<const DoFHandlerType *> dof_handlers;
2443  std::vector<const AffineConstraints<number2> *> constraints;
2444  std::vector<QuadratureType> quads;
2445 
2446  dof_handlers.push_back(&dof_handler);
2447  constraints.push_back(&constraints_in);
2448  quads.push_back(quad);
2449 
2450  std::vector<IndexSet> locally_owned_sets =
2451  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2452  dof_handlers, additional_data.level_mg_handler);
2453 
2454  std::vector<hp::QCollection<1>> quad_hp;
2455  quad_hp.emplace_back(quad);
2456 
2457  internal_reinit(StaticMappingQ1<dim>::mapping,
2458  dof_handlers,
2459  constraints,
2460  locally_owned_sets,
2461  quad_hp,
2462  additional_data);
2463 }
2464 
2465 
2466 
2467 template <int dim, typename Number>
2468 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2469 void
2470 MatrixFree<dim, Number>::reinit(
2471  const Mapping<dim> & mapping,
2472  const DoFHandlerType & dof_handler,
2473  const AffineConstraints<number2> & constraints_in,
2474  const QuadratureType & quad,
2475  const typename MatrixFree<dim, Number>::AdditionalData additional_data)
2476 {
2477  std::vector<const DoFHandlerType *> dof_handlers;
2478  std::vector<const AffineConstraints<number2> *> constraints;
2479 
2480  dof_handlers.push_back(&dof_handler);
2481  constraints.push_back(&constraints_in);
2482 
2483  std::vector<IndexSet> locally_owned_sets =
2484  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2485  dof_handlers, additional_data.level_mg_handler);
2486 
2487  std::vector<hp::QCollection<1>> quad_hp;
2488  quad_hp.emplace_back(quad);
2489 
2490  internal_reinit(mapping,
2491  dof_handlers,
2492  constraints,
2493  locally_owned_sets,
2494  quad_hp,
2495  additional_data);
2496 }
2497 
2498 
2499 
2500 template <int dim, typename Number>
2501 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2502 void
2503 MatrixFree<dim, Number>::reinit(
2504  const std::vector<const DoFHandlerType *> & dof_handler,
2505  const std::vector<const AffineConstraints<number2> *> &constraint,
2506  const std::vector<QuadratureType> & quad,
2507  const typename MatrixFree<dim, Number>::AdditionalData additional_data)
2508 {
2509  std::vector<IndexSet> locally_owned_set =
2510  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2511  dof_handler, additional_data.level_mg_handler);
2512  std::vector<hp::QCollection<1>> quad_hp;
2513  for (unsigned int q = 0; q < quad.size(); ++q)
2514  quad_hp.emplace_back(quad[q]);
2515  internal_reinit(StaticMappingQ1<dim>::mapping,
2516  dof_handler,
2517  constraint,
2518  locally_owned_set,
2519  quad_hp,
2520  additional_data);
2521 }
2522 
2523 
2524 
2525 template <int dim, typename Number>
2526 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2527 void
2528 MatrixFree<dim, Number>::reinit(
2529  const std::vector<const DoFHandlerType *> & dof_handler,
2530  const std::vector<const AffineConstraints<number2> *> &constraint,
2531  const QuadratureType & quad,
2532  const typename MatrixFree<dim, Number>::AdditionalData additional_data)
2533 {
2534  std::vector<IndexSet> locally_owned_set =
2535  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2536  dof_handler, additional_data.level_mg_handler);
2537  std::vector<hp::QCollection<1>> quad_hp;
2538  quad_hp.emplace_back(quad);
2539  internal_reinit(StaticMappingQ1<dim>::mapping,
2540  dof_handler,
2541  constraint,
2542  locally_owned_set,
2543  quad_hp,
2544  additional_data);
2545 }
2546 
2547 
2548 
2549 template <int dim, typename Number>
2550 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2551 void
2552 MatrixFree<dim, Number>::reinit(
2553  const Mapping<dim> & mapping,
2554  const std::vector<const DoFHandlerType *> & dof_handler,
2555  const std::vector<const AffineConstraints<number2> *> &constraint,
2556  const QuadratureType & quad,
2557  const typename MatrixFree<dim, Number>::AdditionalData additional_data)
2558 {
2559  std::vector<IndexSet> locally_owned_set =
2560  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2561  dof_handler, additional_data.level_mg_handler);
2562  std::vector<hp::QCollection<1>> quad_hp;
2563  quad_hp.emplace_back(quad);
2564  internal_reinit(mapping,
2565  dof_handler,
2566  constraint,
2567  locally_owned_set,
2568  quad_hp,
2569  additional_data);
2570 }
2571 
2572 
2573 
2574 template <int dim, typename Number>
2575 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2576 void
2577 MatrixFree<dim, Number>::reinit(
2578  const Mapping<dim> & mapping,
2579  const std::vector<const DoFHandlerType *> & dof_handler,
2580  const std::vector<const AffineConstraints<number2> *> &constraint,
2581  const std::vector<QuadratureType> & quad,
2582  const typename MatrixFree<dim, Number>::AdditionalData 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  for (unsigned int q = 0; q < quad.size(); ++q)
2589  quad_hp.emplace_back(quad[q]);
2590  internal_reinit(mapping,
2591  dof_handler,
2592  constraint,
2593  locally_owned_set,
2594  quad_hp,
2595  additional_data);
2596 }
2597 
2598 
2599 
2600 template <int dim, typename Number>
2601 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2602 void
2603 MatrixFree<dim, Number>::reinit(
2604  const Mapping<dim> & mapping,
2605  const std::vector<const DoFHandlerType *> & dof_handler,
2606  const std::vector<const AffineConstraints<number2> *> &constraint,
2607  const std::vector<IndexSet> & locally_owned_set,
2608  const std::vector<QuadratureType> & quad,
2609  const typename MatrixFree<dim, Number>::AdditionalData additional_data)
2610 {
2611  // find out whether we use a hp Quadrature or a standard quadrature
2612  std::vector<hp::QCollection<1>> quad_hp;
2613  for (unsigned int q = 0; q < quad.size(); ++q)
2614  quad_hp.emplace_back(quad[q]);
2615  internal_reinit(mapping,
2616  dof_handler,
2617  constraint,
2618  locally_owned_set,
2619  quad_hp,
2620  additional_data);
2621 }
2622 
2623 
2624 
2625 // ------------------------------ implementation of loops --------------------
2626 
2627 // internal helper functions that define how to call MPI data exchange
2628 // functions: for generic vectors, do nothing at all. For distributed vectors,
2629 // call update_ghost_values_start function and so on. If we have collections
2630 // of vectors, just do the individual functions of the components. In order to
2631 // keep ghost values consistent (whether we are in read or write mode), we
2632 // also reset the values at the end. the whole situation is a bit complicated
2633 // by the fact that we need to treat block vectors differently, which use some
2634 // additional helper functions to select the blocks and template magic.
2635 namespace internal
2636 {
2640  template <int dim, typename Number>
2641  struct VectorDataExchange
2642  {
2643  // An arbitrary shift for communication to reduce the risk for accidental
2644  // interaction with other open communications that a user program might
2645  // set up
2646  static constexpr unsigned int channel_shift = 103;
2647 
2648 
2649 
2654  VectorDataExchange(
2655  const ::MatrixFree<dim, Number> &matrix_free,
2656  const typename ::MatrixFree<dim, Number>::DataAccessOnFaces
2657  vector_face_access,
2658  const unsigned int n_components)
2659  : matrix_free(matrix_free)
2660  , vector_face_access(
2661  matrix_free.get_task_info().face_partition_data.empty() ?
2662  ::MatrixFree<dim, Number>::DataAccessOnFaces::unspecified :
2663  vector_face_access)
2664  , ghosts_were_set(false)
2665 # ifdef DEAL_II_WITH_MPI
2666  , tmp_data(n_components)
2667  , requests(n_components)
2668 # endif
2669  {
2670  (void)n_components;
2671  if (this->vector_face_access !=
2672  ::MatrixFree<dim, Number>::DataAccessOnFaces::unspecified)
2673  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
2674  AssertDimension(
2675  matrix_free.get_dof_info(c).vector_partitioner_face_variants.size(),
2676  3);
2677  }
2678 
2679 
2680 
2684  ~VectorDataExchange() // NOLINT
2685  {
2686 # ifdef DEAL_II_WITH_MPI
2687  for (unsigned int i = 0; i < tmp_data.size(); ++i)
2688  if (tmp_data[i] != nullptr)
2689  matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
2690 # endif
2691  }
2692 
2693 
2694 
2699  template <typename VectorType>
2700  unsigned int
2701  find_vector_in_mf(const VectorType &vec,
2702  const bool check_global_compatibility = true) const
2703  {
2704  unsigned int mf_component = numbers::invalid_unsigned_int;
2705  (void)check_global_compatibility;
2706  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
2707  if (
2708 # ifdef DEBUG
2709  check_global_compatibility ?
2710  vec.get_partitioner()->is_globally_compatible(
2711  *matrix_free.get_dof_info(c).vector_partitioner) :
2712 # endif
2713  vec.get_partitioner()->is_compatible(
2714  *matrix_free.get_dof_info(c).vector_partitioner))
2715  {
2716  mf_component = c;
2717  break;
2718  }
2719  return mf_component;
2720  }
2721 
2722 
2723 
2728  const Utilities::MPI::Partitioner &
2729  get_partitioner(const unsigned int mf_component) const
2730  {
2731  AssertDimension(matrix_free.get_dof_info(mf_component)
2732  .vector_partitioner_face_variants.size(),
2733  3);
2734  if (vector_face_access ==
2735  ::MatrixFree<dim, Number>::DataAccessOnFaces::none)
2736  return *matrix_free.get_dof_info(mf_component)
2737  .vector_partitioner_face_variants[0];
2738  else if (vector_face_access ==
2739  ::MatrixFree<dim, Number>::DataAccessOnFaces::values)
2740  return *matrix_free.get_dof_info(mf_component)
2741  .vector_partitioner_face_variants[1];
2742  else
2743  return *matrix_free.get_dof_info(mf_component)
2744  .vector_partitioner_face_variants[2];
2745  }
2746 
2747 
2748 
2752  template <typename VectorType,
2753  typename std::enable_if<is_serial_or_dummy<VectorType>::value,
2754  VectorType>::type * = nullptr>
2755  void
2756  update_ghost_values_start(const unsigned int /*component_in_block_vector*/,
2757  const VectorType & /*vec*/)
2758  {}
2759 
2760 
2765  template <typename VectorType,
2766  typename std::enable_if<
2767  !has_update_ghost_values_start<VectorType>::value &&
2768  !is_serial_or_dummy<VectorType>::value,
2769  VectorType>::type * = nullptr>
2770  void
2771  update_ghost_values_start(const unsigned int component_in_block_vector,
2772  const VectorType & vec)
2773  {
2774  (void)component_in_block_vector;
2775  bool ghosts_set = vec.has_ghost_elements();
2776  if (ghosts_set)
2777  ghosts_were_set = true;
2778 
2779  vec.update_ghost_values();
2780  }
2781 
2782 
2783 
2789  template <typename VectorType,
2790  typename std::enable_if<
2791  has_update_ghost_values_start<VectorType>::value &&
2792  !has_exchange_on_subset<VectorType>::value,
2793  VectorType>::type * = nullptr>
2794  void
2795  update_ghost_values_start(const unsigned int component_in_block_vector,
2796  const VectorType & vec)
2797  {
2798  (void)component_in_block_vector;
2799  bool ghosts_set = vec.has_ghost_elements();
2800  if (ghosts_set)
2801  ghosts_were_set = true;
2802 
2803  vec.update_ghost_values_start(component_in_block_vector + channel_shift);
2804  }
2805 
2806 
2807 
2814  template <typename VectorType,
2815  typename std::enable_if<
2816  has_update_ghost_values_start<VectorType>::value &&
2817  has_exchange_on_subset<VectorType>::value,
2818  VectorType>::type * = nullptr>
2819  void
2820  update_ghost_values_start(const unsigned int component_in_block_vector,
2821  const VectorType & vec)
2822  {
2823  static_assert(
2824  std::is_same<Number, typename VectorType::value_type>::value,
2825  "Type mismatch between VectorType and VectorDataExchange");
2826  (void)component_in_block_vector;
2827  bool ghosts_set = vec.has_ghost_elements();
2828  if (ghosts_set)
2829  ghosts_were_set = true;
2830  if (vector_face_access ==
2831  ::MatrixFree<dim, Number>::DataAccessOnFaces::unspecified ||
2832  vec.size() == 0)
2833  vec.update_ghost_values_start(component_in_block_vector +
2834  channel_shift);
2835  else
2836  {
2837 # ifdef DEAL_II_WITH_MPI
2838  const unsigned int mf_component = find_vector_in_mf(vec);
2839  if (&get_partitioner(mf_component) ==
2840  matrix_free.get_dof_info(mf_component).vector_partitioner.get())
2841  {
2842  vec.update_ghost_values_start(component_in_block_vector +
2843  channel_shift);
2844  return;
2845  }
2846 
2847  const Utilities::MPI::Partitioner &part =
2848  get_partitioner(mf_component);
2849  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0)
2850  return;
2851 
2852  tmp_data[component_in_block_vector] =
2853  matrix_free.acquire_scratch_data_non_threadsafe();
2854  tmp_data[component_in_block_vector]->resize_fast(
2855  part.n_import_indices());
2856  AssertDimension(requests.size(), tmp_data.size());
2857 
2858  part.export_to_ghosted_array_start(
2859  component_in_block_vector + channel_shift,
2860  ArrayView<const Number>(vec.begin(), part.local_size()),
2861  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
2862  part.n_import_indices()),
2863  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
2864  vec.get_partitioner()->local_size(),
2865  vec.get_partitioner()->n_ghost_indices()),
2866  this->requests[component_in_block_vector]);
2867 # endif
2868  }
2869  }
2870 
2871 
2872 
2877  template <
2878  typename VectorType,
2879  typename std::enable_if<!has_update_ghost_values_start<VectorType>::value,
2880  VectorType>::type * = nullptr>
2881  void
2882  update_ghost_values_finish(const unsigned int /*component_in_block_vector*/,
2883  const VectorType & /*vec*/)
2884  {}
2885 
2886 
2887 
2893  template <typename VectorType,
2894  typename std::enable_if<
2895  has_update_ghost_values_start<VectorType>::value &&
2896  !has_exchange_on_subset<VectorType>::value,
2897  VectorType>::type * = nullptr>
2898  void
2899  update_ghost_values_finish(const unsigned int component_in_block_vector,
2900  const VectorType & vec)
2901  {
2902  (void)component_in_block_vector;
2903  vec.update_ghost_values_finish();
2904  }
2905 
2906 
2907 
2914  template <typename VectorType,
2915  typename std::enable_if<
2916  has_update_ghost_values_start<VectorType>::value &&
2917  has_exchange_on_subset<VectorType>::value,
2918  VectorType>::type * = nullptr>
2919  void
2920  update_ghost_values_finish(const unsigned int component_in_block_vector,
2921  const VectorType & vec)
2922  {
2923  static_assert(
2924  std::is_same<Number, typename VectorType::value_type>::value,
2925  "Type mismatch between VectorType and VectorDataExchange");
2926  (void)component_in_block_vector;
2927  if (vector_face_access ==
2928  ::MatrixFree<dim, Number>::DataAccessOnFaces::unspecified ||
2929  vec.size() == 0)
2930  vec.update_ghost_values_finish();
2931  else
2932  {
2933 # ifdef DEAL_II_WITH_MPI
2934 
2935  AssertIndexRange(component_in_block_vector, tmp_data.size());
2936  AssertDimension(requests.size(), tmp_data.size());
2937 
2938  const unsigned int mf_component = find_vector_in_mf(vec);
2939  const Utilities::MPI::Partitioner &part =
2940  get_partitioner(mf_component);
2941  if (&part ==
2942  matrix_free.get_dof_info(mf_component).vector_partitioner.get())
2943  {
2944  vec.update_ghost_values_finish();
2945  return;
2946  }
2947 
2948  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0)
2949  return;
2950 
2951  part.export_to_ghosted_array_finish(
2952  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
2953  vec.get_partitioner()->local_size(),
2954  vec.get_partitioner()->n_ghost_indices()),
2955  this->requests[component_in_block_vector]);
2956 
2957  matrix_free.release_scratch_data_non_threadsafe(
2958  tmp_data[component_in_block_vector]);
2959  tmp_data[component_in_block_vector] = nullptr;
2960 
2961  // let vector know that ghosts are being updated and we can read from
2962  // them
2963  vec.set_ghost_state(true);
2964 # endif
2965  }
2966  }
2967 
2968 
2969 
2973  template <typename VectorType,
2974  typename std::enable_if<is_serial_or_dummy<VectorType>::value,
2975  VectorType>::type * = nullptr>
2976  void
2977  compress_start(const unsigned int /*component_in_block_vector*/,
2978  VectorType & /*vec*/)
2979  {}
2980 
2981 
2982 
2987  template <typename VectorType,
2988  typename std::enable_if<!has_compress_start<VectorType>::value &&
2989  !is_serial_or_dummy<VectorType>::value,
2990  VectorType>::type * = nullptr>
2991  void
2992  compress_start(const unsigned int component_in_block_vector,
2993  VectorType & vec)
2994  {
2995  (void)component_in_block_vector;
2996  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
2997  vec.compress(::VectorOperation::add);
2998  }
2999 
3000 
3001 
3007  template <
3008  typename VectorType,
3009  typename std::enable_if<has_compress_start<VectorType>::value &&
3010  !has_exchange_on_subset<VectorType>::value,
3011  VectorType>::type * = nullptr>
3012  void
3013  compress_start(const unsigned int component_in_block_vector,
3014  VectorType & vec)
3015  {
3016  (void)component_in_block_vector;
3017  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3018  vec.compress_start(component_in_block_vector + channel_shift);
3019  }
3020 
3021 
3022 
3029  template <
3030  typename VectorType,
3031  typename std::enable_if<has_compress_start<VectorType>::value &&
3032  has_exchange_on_subset<VectorType>::value,
3033  VectorType>::type * = nullptr>
3034  void
3035  compress_start(const unsigned int component_in_block_vector,
3036  VectorType & vec)
3037  {
3038  static_assert(
3039  std::is_same<Number, typename VectorType::value_type>::value,
3040  "Type mismatch between VectorType and VectorDataExchange");
3041  (void)component_in_block_vector;
3042  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3043  if (vector_face_access ==
3044  ::MatrixFree<dim, Number>::DataAccessOnFaces::unspecified ||
3045  vec.size() == 0)
3046  vec.compress_start(component_in_block_vector + channel_shift);
3047  else
3048  {
3049 # ifdef DEAL_II_WITH_MPI
3050 
3051  const unsigned int mf_component = find_vector_in_mf(vec);
3052  const Utilities::MPI::Partitioner &part =
3053  get_partitioner(mf_component);
3054  if (&part ==
3055  matrix_free.get_dof_info(mf_component).vector_partitioner.get())
3056  {
3057  vec.compress_start(component_in_block_vector + channel_shift);
3058  return;
3059  }
3060 
3061  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0)
3062  return;
3063 
3064  tmp_data[component_in_block_vector] =
3065  matrix_free.acquire_scratch_data_non_threadsafe();
3066  tmp_data[component_in_block_vector]->resize_fast(
3067  part.n_import_indices());
3068  AssertDimension(requests.size(), tmp_data.size());
3069 
3070  part.import_from_ghosted_array_start(
3071  ::VectorOperation::add,
3072  component_in_block_vector + channel_shift,
3073  ArrayView<Number>(vec.begin() + vec.get_partitioner()->local_size(),
3074  vec.get_partitioner()->n_ghost_indices()),
3075  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3076  part.n_import_indices()),
3077  this->requests[component_in_block_vector]);
3078 # endif
3079  }
3080  }
3081 
3082 
3083 
3088  template <typename VectorType,
3089  typename std::enable_if<!has_compress_start<VectorType>::value,
3090  VectorType>::type * = nullptr>
3091  void
3092  compress_finish(const unsigned int /*component_in_block_vector*/,
3093  VectorType & /*vec*/)
3094  {}
3095 
3096 
3097 
3103  template <
3104  typename VectorType,
3105  typename std::enable_if<has_compress_start<VectorType>::value &&
3106  !has_exchange_on_subset<VectorType>::value,
3107  VectorType>::type * = nullptr>
3108  void
3109  compress_finish(const unsigned int component_in_block_vector,
3110  VectorType & vec)
3111  {
3112  (void)component_in_block_vector;
3113  vec.compress_finish(::VectorOperation::add);
3114  }
3115 
3116 
3117 
3124  template <
3125  typename VectorType,
3126  typename std::enable_if<has_compress_start<VectorType>::value &&
3127  has_exchange_on_subset<VectorType>::value,
3128  VectorType>::type * = nullptr>
3129  void
3130  compress_finish(const unsigned int component_in_block_vector,
3131  VectorType & vec)
3132  {
3133  static_assert(
3134  std::is_same<Number, typename VectorType::value_type>::value,
3135  "Type mismatch between VectorType and VectorDataExchange");
3136  (void)component_in_block_vector;
3137  if (vector_face_access ==
3138  ::MatrixFree<dim, Number>::DataAccessOnFaces::unspecified ||
3139  vec.size() == 0)
3140  vec.compress_finish(::VectorOperation::add);
3141  else
3142  {
3143 # ifdef DEAL_II_WITH_MPI
3144  AssertIndexRange(component_in_block_vector, tmp_data.size());
3145  AssertDimension(requests.size(), tmp_data.size());
3146 
3147  const unsigned int mf_component = find_vector_in_mf(vec);
3148 
3149  const Utilities::MPI::Partitioner &part =
3150  get_partitioner(mf_component);
3151  if (&part ==
3152  matrix_free.get_dof_info(mf_component).vector_partitioner.get())
3153  {
3154  vec.compress_finish(::VectorOperation::add);
3155  return;
3156  }
3157 
3158  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0)
3159  return;
3160 
3161  part.import_from_ghosted_array_finish(
3162  VectorOperation::add,
3163  ArrayView<const Number>(
3164  tmp_data[component_in_block_vector]->begin(),
3165  part.n_import_indices()),
3166  ArrayView<Number>(vec.begin(), part.local_size()),
3167  ArrayView<Number>(vec.begin() + vec.get_partitioner()->local_size(),
3168  vec.get_partitioner()->n_ghost_indices()),
3169  this->requests[component_in_block_vector]);
3170 
3171  matrix_free.release_scratch_data_non_threadsafe(
3172  tmp_data[component_in_block_vector]);
3173  tmp_data[component_in_block_vector] = nullptr;
3174 # endif
3175  }
3176  }
3177 
3178 
3179 
3183  template <typename VectorType,
3184  typename std::enable_if<is_serial_or_dummy<VectorType>::value,
3185  VectorType>::type * = nullptr>
3186  void
3187  reset_ghost_values(const VectorType & /*vec*/) const
3188  {}
3189 
3190 
3191 
3196  template <
3197  typename VectorType,
3198  typename std::enable_if<!has_exchange_on_subset<VectorType>::value &&
3199  !is_serial_or_dummy<VectorType>::value,
3200  VectorType>::type * = nullptr>
3201  void
3202  reset_ghost_values(const VectorType &vec) const
3203  {
3204  if (ghosts_were_set == true)
3205  return;
3206 
3207  vec.zero_out_ghosts();
3208  }
3209 
3210 
3211 
3217  template <typename VectorType,
3218  typename std::enable_if<has_exchange_on_subset<VectorType>::value,
3219  VectorType>::type * = nullptr>
3220  void
3221  reset_ghost_values(const VectorType &vec) const
3222  {
3223  static_assert(
3224  std::is_same<Number, typename VectorType::value_type>::value,
3225  "Type mismatch between VectorType and VectorDataExchange");
3226  if (ghosts_were_set == true)
3227  return;
3228 
3229  if (vector_face_access ==
3230  ::MatrixFree<dim, Number>::DataAccessOnFaces::unspecified ||
3231  vec.size() == 0)
3232  vec.zero_out_ghosts();
3233  else
3234  {
3235 # ifdef DEAL_II_WITH_MPI
3236  AssertDimension(requests.size(), tmp_data.size());
3237 
3238  const unsigned int mf_component = find_vector_in_mf(vec);
3239  const Utilities::MPI::Partitioner &part =
3240  get_partitioner(mf_component);
3241  if (&part ==
3242  matrix_free.get_dof_info(mf_component).vector_partitioner.get())
3243  vec.zero_out_ghosts();
3244  else if (part.n_ghost_indices() > 0)
3245  {
3246  for (std::vector<std::pair<unsigned int, unsigned int>>::
3247  const_iterator my_ghosts =
3248  part.ghost_indices_within_larger_ghost_set().begin();
3249  my_ghosts !=
3250  part.ghost_indices_within_larger_ghost_set().end();
3251  ++my_ghosts)
3252  for (unsigned int j = my_ghosts->first; j < my_ghosts->second;
3253  j++)
3254  {
3255  const_cast<LinearAlgebra::distributed::Vector<Number> &>(
3256  vec)
3257  .local_element(j + part.local_size()) = 0.;
3258  }
3259  }
3260 
3261  // let vector know that it's not ghosted anymore
3262  vec.set_ghost_state(false);
3263 # endif
3264  }
3265  }
3266 
3267 
3268 
3274  template <typename VectorType,
3275  typename std::enable_if<has_exchange_on_subset<VectorType>::value,
3276  VectorType>::type * = nullptr>
3277  void
3278  zero_vector_region(const unsigned int range_index, VectorType &vec) const
3279  {
3280  static_assert(
3281  std::is_same<Number, typename VectorType::value_type>::value,
3282  "Type mismatch between VectorType and VectorDataExchange");
3283  if (range_index == numbers::invalid_unsigned_int)
3284  vec = Number();
3285  else
3286  {
3287  const unsigned int mf_component = find_vector_in_mf(vec, false);
3288  const internal::MatrixFreeFunctions::DoFInfo &dof_info =
3289  matrix_free.get_dof_info(mf_component);
3290  Assert(dof_info.vector_zero_range_list_index.empty() == false,
3291  ExcNotInitialized());
3292 
3293  Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3294  ExcInternalError());
3295  AssertIndexRange(range_index,
3296  dof_info.vector_zero_range_list_index.size() - 1);
3297  for (unsigned int id =
3298  dof_info.vector_zero_range_list_index[range_index];
3299  id != dof_info.vector_zero_range_list_index[range_index + 1];
3300  ++id)
3301  {
3302  const unsigned int start_pos =
3303  dof_info.vector_zero_range_list[id] *
3304  internal::MatrixFreeFunctions::DoFInfo::chunk_size_zero_vector;
3305  const unsigned int end_pos =
3306  std::min((dof_info.vector_zero_range_list[id] + 1) *
3307  internal::MatrixFreeFunctions::DoFInfo::
3308  chunk_size_zero_vector,
3309  dof_info.vector_partitioner->local_size() +
3310  dof_info.vector_partitioner->n_ghost_indices());
3311  std::memset(vec.begin() + start_pos,
3312  0,
3313  (end_pos - start_pos) * sizeof(Number));
3314  }
3315  }
3316  }
3317 
3318 
3319 
3325  template <
3326  typename VectorType,
3327  typename std::enable_if<!has_exchange_on_subset<VectorType>::value,
3328  VectorType>::type * = nullptr,
3329  typename VectorType::value_type * = nullptr>
3330  void
3331  zero_vector_region(const unsigned int range_index, VectorType &vec) const
3332  {
3333  if (range_index == numbers::invalid_unsigned_int)
3334  vec = typename VectorType::value_type();
3335  else
3336  {
3337  Assert(false, ExcNotImplemented());
3338  }
3339  }
3340 
3341 
3342 
3347  void
3348  zero_vector_region(const unsigned int, ...) const
3349  {
3350  Assert(false,
3351  ExcNotImplemented("Zeroing is only implemented for vector types "
3352  "which provide VectorType::value_type"));
3353  }
3354 
3355 
3356 
3357  const ::MatrixFree<dim, Number> &matrix_free;
3358  const typename ::MatrixFree<dim, Number>::DataAccessOnFaces
3359  vector_face_access;
3360  bool ghosts_were_set;
3361 # ifdef DEAL_II_WITH_MPI
3362  std::vector<AlignedVector<Number> *> tmp_data;
3363  std::vector<std::vector<MPI_Request>> requests;
3364 # endif
3365  }; // VectorDataExchange
3366 
3367  template <typename VectorStruct>
3368  unsigned int
3369  n_components(const VectorStruct &vec);
3370 
3371  template <typename VectorStruct>
3372  unsigned int
3373  n_components_block(const VectorStruct &vec,
3374  std::integral_constant<bool, true>)
3375  {
3376  unsigned int components = 0;
3377  for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
3378  components += n_components(vec.block(bl));
3379  return components;
3380  }
3381 
3382  template <typename VectorStruct>
3383  unsigned int
3384  n_components_block(const VectorStruct &, std::integral_constant<bool, false>)
3385  {
3386  return 1;
3387  }
3388 
3389  template <typename VectorStruct>
3390  unsigned int
3391  n_components(const VectorStruct &vec)
3392  {
3393  return n_components_block(
3394  vec, std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3395  }
3396 
3397  template <typename VectorStruct>
3398  inline unsigned int
3399  n_components(const std::vector<VectorStruct> &vec)
3400  {
3401  unsigned int components = 0;
3402  for (unsigned int comp = 0; comp < vec.size(); comp++)
3403  components += n_components_block(
3404  vec[comp],
3405  std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3406  return components;
3407  }
3408 
3409  template <typename VectorStruct>
3410  inline unsigned int
3411  n_components(const std::vector<VectorStruct *> &vec)
3412  {
3413  unsigned int components = 0;
3414  for (unsigned int comp = 0; comp < vec.size(); comp++)
3415  components += n_components_block(
3416  *vec[comp],
3417  std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3418  return components;
3419  }
3420 
3421 
3422 
3423  // A helper function to identify block vectors with many components where we
3424  // should not try to overlap computations and communication because there
3425  // would be too many outstanding communication requests.
3426 
3427  // default value for vectors that do not have communication_block_size
3428  template <
3429  typename VectorStruct,
3430  typename std::enable_if<!has_communication_block_size<VectorStruct>::value,
3431  VectorStruct>::type * = nullptr>
3432  constexpr unsigned int
3433  get_communication_block_size(const VectorStruct &)
3434  {
3435  return numbers::invalid_unsigned_int;
3436  }
3437 
3438 
3439 
3440  template <
3441  typename VectorStruct,
3442  typename std::enable_if<has_communication_block_size<VectorStruct>::value,
3443  VectorStruct>::type * = nullptr>
3444  constexpr unsigned int
3445  get_communication_block_size(const VectorStruct &)
3446  {
3447  return VectorStruct::communication_block_size;
3448  }
3449 
3450 
3451 
3452  // --------------------------------------------------------------------------
3453  // below we have wrappers to distinguish between block and non-block vectors.
3454  // --------------------------------------------------------------------------
3455 
3456  //
3457  // update_ghost_values_start
3458  //
3459 
3460  // update_ghost_values for block vectors
3461  template <int dim,
3462  typename VectorStruct,
3463  typename Number,
3464  typename std::enable_if<IsBlockVector<VectorStruct>::value,
3465  VectorStruct>::type * = nullptr>
3466  void
3467  update_ghost_values_start(const VectorStruct & vec,
3468  VectorDataExchange<dim, Number> &exchanger,
3469  const unsigned int channel = 0)
3470  {
3471  if (get_communication_block_size(vec) < vec.n_blocks())
3472  {
3473  // don't forget to set ghosts_were_set, that otherwise happens
3474  // inside VectorDataExchange::update_ghost_values_start()
3475  exchanger.ghosts_were_set = vec.has_ghost_elements();
3476  vec.update_ghost_values();
3477  }
3478  else
3479  {
3480  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
3481  update_ghost_values_start(vec.block(i), exchanger, channel + i);
3482  }
3483  }
3484 
3485 
3486 
3487  // update_ghost_values for non-block vectors
3488  template <int dim,
3489  typename VectorStruct,
3490  typename Number,
3491  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
3492  VectorStruct>::type * = nullptr>
3493  void
3494  update_ghost_values_start(const VectorStruct & vec,
3495  VectorDataExchange<dim, Number> &exchanger,
3496  const unsigned int channel = 0)
3497  {
3498  exchanger.update_ghost_values_start(channel, vec);
3499  }
3500 
3501 
3502 
3503  // update_ghost_values_start() for vector of vectors
3504  template <int dim, typename VectorStruct, typename Number>
3505  inline void
3506  update_ghost_values_start(const std::vector<VectorStruct> &vec,
3507  VectorDataExchange<dim, Number> &exchanger)
3508  {
3509  unsigned int component_index = 0;
3510  for (unsigned int comp = 0; comp < vec.size(); comp++)
3511  {
3512  update_ghost_values_start(vec[comp], exchanger, component_index);
3513  component_index += n_components(vec[comp]);
3514  }
3515  }
3516 
3517 
3518 
3519  // update_ghost_values_start() for vector of pointers to vectors
3520  template <int dim, typename VectorStruct, typename Number>
3521  inline void
3522  update_ghost_values_start(const std::vector<VectorStruct *> &vec,
3523  VectorDataExchange<dim, Number> & exchanger)
3524  {
3525  unsigned int component_index = 0;
3526  for (unsigned int comp = 0; comp < vec.size(); comp++)
3527  {
3528  update_ghost_values_start(*vec[comp], exchanger, component_index);
3529  component_index += n_components(*vec[comp]);
3530  }
3531  }
3532 
3533 
3534 
3535  //
3536  // update_ghost_values_finish
3537  //
3538 
3539  // for block vectors
3540  template <int dim,
3541  typename VectorStruct,
3542  typename Number,
3543  typename std::enable_if<IsBlockVector<VectorStruct>::value,
3544  VectorStruct>::type * = nullptr>
3545  void
3546  update_ghost_values_finish(const VectorStruct & vec,
3547  VectorDataExchange<dim, Number> &exchanger,
3548  const unsigned int channel = 0)
3549  {
3550  if (get_communication_block_size(vec) < vec.n_blocks())
3551  {
3552  // do nothing, everything has already been completed in the _start()
3553  // call
3554  }
3555  else
3556  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
3557  update_ghost_values_finish(vec.block(i), exchanger, channel + i);
3558  }
3559 
3560 
3561 
3562  // for non-block vectors
3563  template <int dim,
3564  typename VectorStruct,
3565  typename Number,
3566  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
3567  VectorStruct>::type * = nullptr>
3568  void
3569  update_ghost_values_finish(const VectorStruct & vec,
3570  VectorDataExchange<dim, Number> &exchanger,
3571  const unsigned int channel = 0)
3572  {
3573  exchanger.update_ghost_values_finish(channel, vec);
3574  }
3575 
3576 
3577 
3578  // for vector of vectors
3579  template <int dim, typename VectorStruct, typename Number>
3580  inline void
3581  update_ghost_values_finish(const std::vector<VectorStruct> &vec,
3582  VectorDataExchange<dim, Number> &exchanger)
3583  {
3584  unsigned int component_index = 0;
3585  for (unsigned int comp = 0; comp < vec.size(); comp++)
3586  {
3587  update_ghost_values_finish(vec[comp], exchanger, component_index);
3588  component_index += n_components(vec[comp]);
3589  }
3590  }
3591 
3592 
3593 
3594  // for vector of pointers to vectors
3595  template <int dim, typename VectorStruct, typename Number>
3596  inline void
3597  update_ghost_values_finish(const std::vector<VectorStruct *> &vec,
3598  VectorDataExchange<dim, Number> & exchanger)
3599  {
3600  unsigned int component_index = 0;
3601  for (unsigned int comp = 0; comp < vec.size(); comp++)
3602  {
3603  update_ghost_values_finish(*vec[comp], exchanger, component_index);
3604  component_index += n_components(*vec[comp]);
3605  }
3606  }
3607 
3608 
3609 
3610  //
3611  // compress_start
3612  //
3613 
3614  // for block vectors
3615  template <int dim,
3616  typename VectorStruct,
3617  typename Number,
3618  typename std::enable_if<IsBlockVector<VectorStruct>::value,
3619  VectorStruct>::type * = nullptr>
3620  inline void
3621  compress_start(VectorStruct & vec,
3622  VectorDataExchange<dim, Number> &exchanger,
3623  const unsigned int channel = 0)
3624  {
3625  if (get_communication_block_size(vec) < vec.n_blocks())
3626  vec.compress(::VectorOperation::add);
3627  else
3628  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
3629  compress_start(vec.block(i), exchanger, channel + i);
3630  }
3631 
3632 
3633 
3634  // for non-block vectors
3635  template <int dim,
3636  typename VectorStruct,
3637  typename Number,
3638  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
3639  VectorStruct>::type * = nullptr>
3640  inline void
3641  compress_start(VectorStruct & vec,
3642  VectorDataExchange<dim, Number> &exchanger,
3643  const unsigned int channel = 0)
3644  {
3645  exchanger.compress_start(channel, vec);
3646  }
3647 
3648 
3649 
3650  // for std::vector of vectors
3651  template <int dim, typename VectorStruct, typename Number>
3652  inline void
3653  compress_start(std::vector<VectorStruct> & vec,
3654  VectorDataExchange<dim, Number> &exchanger)
3655  {
3656  unsigned int component_index = 0;
3657  for (unsigned int comp = 0; comp < vec.size(); comp++)
3658  {
3659  compress_start(vec[comp], exchanger, component_index);
3660  component_index += n_components(vec[comp]);
3661  }
3662  }
3663 
3664 
3665 
3666  // for std::vector of pointer to vectors
3667  template <int dim, typename VectorStruct, typename Number>
3668  inline void
3669  compress_start(std::vector<VectorStruct *> & vec,
3670  VectorDataExchange<dim, Number> &exchanger)
3671  {
3672  unsigned int component_index = 0;
3673  for (unsigned int comp = 0; comp < vec.size(); comp++)
3674  {
3675  compress_start(*vec[comp], exchanger, component_index);
3676  component_index += n_components(*vec[comp]);
3677  }
3678  }
3679 
3680 
3681 
3682  //
3683  // compress_finish
3684  //
3685 
3686  // for block vectors
3687  template <int dim,
3688  typename VectorStruct,
3689  typename Number,
3690  typename std::enable_if<IsBlockVector<VectorStruct>::value,
3691  VectorStruct>::type * = nullptr>
3692  inline void
3693  compress_finish(VectorStruct & vec,
3694  VectorDataExchange<dim, Number> &exchanger,
3695  const unsigned int channel = 0)
3696  {
3697  if (get_communication_block_size(vec) < vec.n_blocks())
3698  {
3699  // do nothing, everything has already been completed in the _start()
3700  // call
3701  }
3702  else
3703  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
3704  compress_finish(vec.block(i), exchanger, channel + i);
3705  }
3706 
3707 
3708 
3709  // for non-block vectors
3710  template <int dim,
3711  typename VectorStruct,
3712  typename Number,
3713  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
3714  VectorStruct>::type * = nullptr>
3715  inline void
3716  compress_finish(VectorStruct & vec,
3717  VectorDataExchange<dim, Number> &exchanger,
3718  const unsigned int channel = 0)
3719  {
3720  exchanger.compress_finish(channel, vec);
3721  }
3722 
3723 
3724 
3725  // for std::vector of vectors
3726  template <int dim, typename VectorStruct, typename Number>
3727  inline void
3728  compress_finish(std::vector<VectorStruct> & vec,
3729  VectorDataExchange<dim, Number> &exchanger)
3730  {
3731  unsigned int component_index = 0;
3732  for (unsigned int comp = 0; comp < vec.size(); comp++)
3733  {
3734  compress_finish(vec[comp], exchanger, component_index);
3735  component_index += n_components(vec[comp]);
3736  }
3737  }
3738 
3739 
3740 
3741  // for std::vector of pointer to vectors
3742  template <int dim, typename VectorStruct, typename Number>
3743  inline void
3744  compress_finish(std::vector<VectorStruct *> & vec,
3745  VectorDataExchange<dim, Number> &exchanger)
3746  {
3747  unsigned int component_index = 0;
3748  for (unsigned int comp = 0; comp < vec.size(); comp++)
3749  {
3750  compress_finish(*vec[comp], exchanger, component_index);
3751  component_index += n_components(*vec[comp]);
3752  }
3753  }
3754 
3755 
3756 
3757  //
3758  // reset_ghost_values:
3759  //
3760  // if the input vector did not have ghosts imported, clear them here again
3761  // in order to avoid subsequent operations e.g. in linear solvers to work
3762  // with ghosts all the time
3763  //
3764 
3765  // for block vectors
3766  template <int dim,
3767  typename VectorStruct,
3768  typename Number,
3769  typename std::enable_if<IsBlockVector<VectorStruct>::value,
3770  VectorStruct>::type * = nullptr>
3771  inline void
3772  reset_ghost_values(const VectorStruct & vec,
3773  VectorDataExchange<dim, Number> &exchanger)
3774  {
3775  // return immediately if there is nothing to do.
3776  if (exchanger.ghosts_were_set == true)
3777  return;
3778 
3779  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
3780  reset_ghost_values(vec.block(i), exchanger);
3781  }
3782 
3783 
3784 
3785  // for non-block vectors
3786  template <int dim,
3787  typename VectorStruct,
3788  typename Number,
3789  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
3790  VectorStruct>::type * = nullptr>
3791  inline void
3792  reset_ghost_values(const VectorStruct & vec,
3793  VectorDataExchange<dim, Number> &exchanger)
3794  {
3795  exchanger.reset_ghost_values(vec);
3796  }
3797 
3798 
3799 
3800  // for std::vector of vectors
3801  template <int dim, typename VectorStruct, typename Number>
3802  inline void
3803  reset_ghost_values(const std::vector<VectorStruct> &vec,
3804  VectorDataExchange<dim, Number> &exchanger)
3805  {
3806  // return immediately if there is nothing to do.
3807  if (exchanger.ghosts_were_set == true)
3808  return;
3809 
3810  for (unsigned int comp = 0; comp < vec.size(); comp++)
3811  reset_ghost_values(vec[comp], exchanger);
3812  }
3813 
3814 
3815 
3816  // for std::vector of pointer to vectors
3817  template <int dim, typename VectorStruct, typename Number>
3818  inline void
3819  reset_ghost_values(const std::vector<VectorStruct *> &vec,
3820  VectorDataExchange<dim, Number> & exchanger)
3821  {
3822  // return immediately if there is nothing to do.
3823  if (exchanger.ghosts_were_set == true)
3824  return;
3825 
3826  for (unsigned int comp = 0; comp < vec.size(); comp++)
3827  reset_ghost_values(*vec[comp], exchanger);
3828  }
3829 
3830 
3831 
3832  //
3833  // zero_vector_region
3834  //
3835 
3836  // for block vectors
3837  template <int dim,
3838  typename VectorStruct,
3839  typename Number,
3840  typename std::enable_if<IsBlockVector<VectorStruct>::value,
3841  VectorStruct>::type * = nullptr>
3842  inline void
3843  zero_vector_region(const unsigned int range_index,
3844  VectorStruct & vec,
3845  VectorDataExchange<dim, Number> &exchanger)
3846  {
3847  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
3848  exchanger.zero_vector_region(range_index, vec.block(i));
3849  }
3850 
3851 
3852 
3853  // for non-block vectors
3854  template <int dim,
3855  typename VectorStruct,
3856  typename Number,
3857  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
3858  VectorStruct>::type * = nullptr>
3859  inline void
3860  zero_vector_region(const unsigned int range_index,
3861  VectorStruct & vec,
3862  VectorDataExchange<dim, Number> &exchanger)
3863  {
3864  exchanger.zero_vector_region(range_index, vec);
3865  }
3866 
3867 
3868 
3869  // for std::vector of vectors
3870  template <int dim, typename VectorStruct, typename Number>
3871  inline void
3872  zero_vector_region(const unsigned int range_index,
3873  std::vector<VectorStruct> & vec,
3874  VectorDataExchange<dim, Number> &exchanger)
3875  {
3876  for (unsigned int comp = 0; comp < vec.size(); comp++)
3877  zero_vector_region(range_index, vec[comp], exchanger);
3878  }
3879 
3880 
3881 
3882  // for std::vector of pointers to vectors
3883  template <int dim, typename VectorStruct, typename Number>
3884  inline void
3885  zero_vector_region(const unsigned int range_index,
3886  std::vector<VectorStruct *> & vec,
3887  VectorDataExchange<dim, Number> &exchanger)
3888  {
3889  for (unsigned int comp = 0; comp < vec.size(); comp++)
3890  zero_vector_region(range_index, *vec[comp], exchanger);
3891  }
3892 
3893 
3894 
3895  namespace MatrixFreeFunctions
3896  {
3897  // struct to select between a const interface and a non-const interface
3898  // for MFWorker
3899  template <typename, typename, typename, typename, bool>
3900  struct InterfaceSelector
3901  {};
3902 
3903  // Version for constant functions
3904  template <typename MF,
3905  typename InVector,
3906  typename OutVector,
3907  typename Container>
3908  struct InterfaceSelector<MF, InVector, OutVector, Container, true>
3909  {
3910  using function_type = void (Container::*)(
3911  const MF &,
3912  OutVector &,
3913  const InVector &,
3914  const std::pair<unsigned int, unsigned int> &) const;
3915  };
3916 
3917  // Version for non-constant functions
3918  template <typename MF,
3919  typename InVector,
3920  typename OutVector,
3921  typename Container>
3922  struct InterfaceSelector<MF, InVector, OutVector, Container, false>
3923  {
3924  using function_type =
3925  void (Container::*)(const MF &,
3926  OutVector &,
3927  const InVector &,
3928  const std::pair<unsigned int, unsigned int> &);
3929  };
3930  } // namespace MatrixFreeFunctions
3931 
3932 
3933 
3934  // A implementation class for the worker object that runs the various
3935  // operations we want to perform during the matrix-free loop
3936  template <typename MF,
3937  typename InVector,
3938  typename OutVector,
3939  typename Container,
3940  bool is_constant>
3941  class MFWorker : public MFWorkerInterface
3942  {
3943  public:
3944  // An alias to make the arguments further down more readable
3945  using function_type = typename MatrixFreeFunctions::
3946  InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
3947  function_type;
3948 
3949  // constructor, binds all the arguments to this class
3950  MFWorker(const MF & matrix_free,
3951  const InVector & src,
3952  OutVector & dst,
3953  const bool zero_dst_vector_setting,
3954  const Container & container,
3955  function_type cell_function,
3956  function_type face_function,
3957  function_type boundary_function,
3958  const typename MF::DataAccessOnFaces src_vector_face_access =
3959  MF::DataAccessOnFaces::none,
3960  const typename MF::DataAccessOnFaces dst_vector_face_access =
3961  MF::DataAccessOnFaces::none)
3962  : matrix_free(matrix_free)
3963  , container(const_cast<Container &>(container))
3964  , cell_function(cell_function)
3965  , face_function(face_function)
3966  , boundary_function(boundary_function)
3967  , src(src)
3968  , dst(dst)
3969  , src_data_exchanger(matrix_free,
3970  src_vector_face_access,
3971  n_components(src))
3972  , dst_data_exchanger(matrix_free,
3973  dst_vector_face_access,
3974  n_components(dst))
3975  , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
3976  , zero_dst_vector_setting(zero_dst_vector_setting &&
3977  !src_and_dst_are_same)
3978  {}
3979 
3980  // Runs the cell work. If no function is given, nothing is done
3981  virtual void
3982  cell(const std::pair<unsigned int, unsigned int> &cell_range) override
3983  {
3984  if (cell_function != nullptr && cell_range.second > cell_range.first)
3985  (container.*
3986  cell_function)(matrix_free, this->dst, this->src, cell_range);
3987  }
3988 
3989  // Runs the assembler on interior faces. If no function is given, nothing
3990  // is done
3991  virtual void
3992  face(const std::pair<unsigned int, unsigned int> &face_range) override
3993  {
3994  if (face_function != nullptr && face_range.second > face_range.first)
3995  (container.*
3996  face_function)(matrix_free, this->dst, this->src, face_range);
3997  }
3998 
3999  // Runs the assembler on boundary faces. If no function is given, nothing
4000  // is done
4001  virtual void
4002  boundary(const std::pair<unsigned int, unsigned int> &face_range) override
4003  {
4004  if (boundary_function != nullptr && face_range.second > face_range.first)
4005  (container.*
4006  boundary_function)(matrix_free, this->dst, this->src, face_range);
4007  }
4008 
4009  // Starts the communication for the update ghost values operation. We
4010  // cannot call this update if ghost and destination are the same because
4011  // that would introduce spurious entries in the destination (there is also
4012  // the problem that reading from a vector that we also write to is usually
4013  // not intended in case there is overlap, but this is up to the
4014  // application code to decide and we cannot catch this case here).
4015  virtual void
4016  vector_update_ghosts_start() override
4017  {
4018  if (!src_and_dst_are_same)
4019  internal::update_ghost_values_start(src, src_data_exchanger);
4020  }
4021 
4022  // Finishes the communication for the update ghost values operation
4023  virtual void
4024  vector_update_ghosts_finish() override
4025  {
4026  if (!src_and_dst_are_same)
4027  internal::update_ghost_values_finish(src, src_data_exchanger);
4028  }
4029 
4030  // Starts the communication for the vector compress operation
4031  virtual void
4032  vector_compress_start() override
4033  {
4034  internal::compress_start(dst, dst_data_exchanger);
4035  }
4036 
4037  // Finishes the communication for the vector compress operation
4038  virtual void
4039  vector_compress_finish() override
4040  {
4041  internal::compress_finish(dst, dst_data_exchanger);
4042  if (!src_and_dst_are_same)
4043  internal::reset_ghost_values(src, src_data_exchanger);
4044  }
4045 
4046  // Zeros the given input vector
4047  virtual void
4048  zero_dst_vector_range(const unsigned int range_index) override
4049  {
4050  if (zero_dst_vector_setting)
4051  internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4052  }
4053 
4054  private:
4055  const MF & matrix_free;
4056  Container & container;
4057  function_type cell_function;
4058  function_type face_function;
4059  function_type boundary_function;
4060 
4061  const InVector &src;
4062  OutVector & dst;
4063  VectorDataExchange<MF::dimension, typename MF::value_type>
4064  src_data_exchanger;
4065  VectorDataExchange<MF::dimension, typename MF::value_type>
4066  dst_data_exchanger;
4067  const bool src_and_dst_are_same;
4068  const bool zero_dst_vector_setting;
4069  };
4070 
4071 
4072 
4077  template <class MF, typename InVector, typename OutVector>
4078  struct MFClassWrapper
4079  {
4080  using function_type =
4081  std::function<void(const MF &,
4082  OutVector &,
4083  const InVector &,
4084  const std::pair<unsigned int, unsigned int> &)>;
4085 
4086  MFClassWrapper(const function_type cell,
4087  const function_type face,
4088  const function_type boundary)
4089  : cell(cell)
4090  , face(face)
4091  , boundary(boundary)
4092  {}
4093 
4094  void
4095  cell_integrator(const MF & mf,
4096  OutVector & dst,
4097  const InVector & src,
4098  const std::pair<unsigned int, unsigned int> &range) const
4099  {
4100  if (cell)
4101  cell(mf, dst, src, range);
4102  }
4103 
4104  void
4105  face_integrator(const MF & mf,
4106  OutVector & dst,
4107  const InVector & src,
4108  const std::pair<unsigned int, unsigned int> &range) const
4109  {
4110  if (face)
4111  face(mf, dst, src, range);
4112  }
4113 
4114  void
4115  boundary_integrator(
4116  const MF & mf,
4117  OutVector & dst,
4118  const InVector & src,
4119  const std::pair<unsigned int, unsigned int> &range) const
4120  {
4121  if (boundary)
4122  boundary(mf, dst, src, range);
4123  }
4124 
4125  const function_type cell;
4126  const function_type face;
4127  const function_type boundary;
4128  };
4129 
4130 } // end of namespace internal
4131 
4132 
4133 
4134 template <int dim, typename Number>
4135 template <typename OutVector, typename InVector>
4136 inline void
4137 MatrixFree<dim, Number>::cell_loop(
4138  const std::function<void(const MatrixFree<dim, Number> &,
4139  OutVector &,
4140  const InVector &,
4141  const std::pair<unsigned int, unsigned int> &)>
4142  & cell_operation,
4143  OutVector & dst,
4144  const InVector &src,
4145  const bool zero_dst_vector) const
4146 {
4147  using Wrapper =
4148  internal::MFClassWrapper<MatrixFree<dim, Number>, InVector, OutVector>;
4149  Wrapper wrap(cell_operation, nullptr, nullptr);
4150  internal::
4151  MFWorker<MatrixFree<dim, Number>, InVector, OutVector, Wrapper, true>
4152  worker(*this,
4153  src,
4154  dst,
4155  zero_dst_vector,
4156  wrap,
4157  &Wrapper::cell_integrator,
4158  &Wrapper::face_integrator,
4159  &Wrapper::boundary_integrator);
4160 
4161  task_info.loop(worker);
4162 }
4163 
4164 
4165 
4166 template <int dim, typename Number>
4167 template <typename OutVector, typename InVector>
4168 inline void
4169 MatrixFree<dim, Number>::loop(
4170  const std::function<void(const MatrixFree<dim, Number> &,
4171  OutVector &,
4172  const InVector &,
4173  const std::pair<unsigned int, unsigned int> &)>
4174  &cell_operation,
4175  const std::function<void(const MatrixFree<dim, Number> &,
4176  OutVector &,
4177  const InVector &,
4178  const std::pair<unsigned int, unsigned int> &)>
4179  &face_operation,
4180  const std::function<void(const MatrixFree<dim, Number> &,
4181  OutVector &,
4182  const InVector &,
4183  const std::pair<unsigned int, unsigned int> &)>
4184  & boundary_operation,
4185  OutVector & dst,
4186  const InVector & src,
4187  const bool zero_dst_vector,
4188  const DataAccessOnFaces dst_vector_face_access,
4189  const DataAccessOnFaces src_vector_face_access) const
4190 {
4191  using Wrapper =
4192  internal::MFClassWrapper<MatrixFree<dim, Number>, InVector, OutVector>;
4193  Wrapper wrap(cell_operation, face_operation, boundary_operation);
4194  internal::
4195  MFWorker<MatrixFree<dim, Number>, InVector, OutVector, Wrapper, true>
4196  worker(*this,
4197  src,
4198  dst,
4199  zero_dst_vector,
4200  wrap,
4201  &Wrapper::cell_integrator,
4202  &Wrapper::face_integrator,
4203  &Wrapper::boundary_integrator,
4204  src_vector_face_access,
4205  dst_vector_face_access);
4206 
4207  task_info.loop(worker);
4208 }
4209 
4210 
4211 
4212 template <int dim, typename Number>
4213 template <typename CLASS, typename OutVector, typename InVector>
4214 inline void
4215 MatrixFree<dim, Number>::cell_loop(
4216  void (CLASS::*function_pointer)(const MatrixFree<dim, Number> &,
4217  OutVector &,
4218  const InVector &,
4219  const std::pair<unsigned int, unsigned int> &)
4220  const,
4221  const CLASS * owning_class,
4222  OutVector & dst,
4223  const InVector &src,
4224  const bool zero_dst_vector) const
4225 {
4226  internal::MFWorker<MatrixFree<dim, Number>, InVector, OutVector, CLASS, true>
4227  worker(*this,
4228  src,
4229  dst,
4230  zero_dst_vector,
4231  *owning_class,
4232  function_pointer,
4233  nullptr,
4234  nullptr);
4235  task_info.loop(worker);
4236 }
4237 
4238 
4239 
4240 template <int dim, typename Number>
4241 template <typename CLASS, typename OutVector, typename InVector>
4242 inline void
4243 MatrixFree<dim, Number>::loop(
4244  void (CLASS::*cell_operation)(const MatrixFree<dim, Number> &,
4245  OutVector &,
4246  const InVector &,
4247  const std::pair<unsigned int, unsigned int> &)
4248  const,
4249  void (CLASS::*face_operation)(const MatrixFree<dim, Number> &,
4250  OutVector &,
4251  const InVector &,
4252  const std::pair<unsigned int, unsigned int> &)
4253  const,
4254  void (CLASS::*boundary_operation)(
4255  const MatrixFree<dim, Number> &,
4256  OutVector &,
4257  const InVector &,
4258  const std::pair<unsigned int, unsigned int> &) const,
4259  const CLASS * owning_class,
4260  OutVector & dst,
4261  const InVector & src,
4262  const bool zero_dst_vector,
4263  const DataAccessOnFaces dst_vector_face_access,
4264  const DataAccessOnFaces src_vector_face_access) const
4265 {
4266  internal::MFWorker<MatrixFree<dim, Number>, InVector, OutVector, CLASS, true>
4267  worker(*this,
4268  src,
4269  dst,
4270  zero_dst_vector,
4271  *owning_class,
4272  cell_operation,
4273  face_operation,
4274  boundary_operation,
4275  src_vector_face_access,
4276  dst_vector_face_access);
4277  task_info.loop(worker);
4278 }
4279 
4280 
4281 
4282 template <int dim, typename Number>
4283 template <typename CLASS, typename OutVector, typename InVector>
4284 inline void
4285 MatrixFree<dim, Number>::cell_loop(
4286  void (CLASS::*function_pointer)(
4287  const MatrixFree<dim, Number> &,
4288  OutVector &,
4289  const InVector &,
4290  const std::pair<unsigned int, unsigned int> &),
4291  CLASS * owning_class,
4292  OutVector & dst,
4293  const InVector &src,
4294  const bool zero_dst_vector) const
4295 {
4296  internal::MFWorker<MatrixFree<dim, Number>, InVector, OutVector, CLASS, false>
4297  worker(*this,
4298  src,
4299  dst,
4300  zero_dst_vector,
4301  *owning_class,
4302  function_pointer,
4303  nullptr,
4304  nullptr);
4305  task_info.loop(worker);
4306 }
4307 
4308 
4309 
4310 template <int dim, typename Number>
4311 template <typename CLASS, typename OutVector, typename InVector>
4312 inline void
4313 MatrixFree<dim, Number>::loop(
4314  void (CLASS::*cell_operation)(const MatrixFree<dim, Number> &,
4315  OutVector &,
4316  const InVector &,
4317  const std::pair<unsigned int, unsigned int> &),
4318  void (CLASS::*face_operation)(const MatrixFree<dim, Number> &,
4319  OutVector &,
4320  const InVector &,
4321  const std::pair<unsigned int, unsigned int> &),
4322  void (CLASS::*boundary_operation)(
4323  const MatrixFree<dim, Number> &,
4324  OutVector &,
4325  const InVector &,
4326  const std::pair<unsigned int, unsigned int> &),
4327  CLASS * owning_class,
4328  OutVector & dst,
4329  const InVector & src,
4330  const bool zero_dst_vector,
4331  const DataAccessOnFaces dst_vector_face_access,
4332  const DataAccessOnFaces src_vector_face_access) const
4333 {
4334  internal::MFWorker<MatrixFree<dim, Number>, InVector, OutVector, CLASS, false>
4335  worker(*this,
4336  src,
4337  dst,
4338  zero_dst_vector,
4339  *owning_class,
4340  cell_operation,
4341  face_operation,
4342  boundary_operation,
4343  src_vector_face_access,
4344  dst_vector_face_access);
4345  task_info.loop(worker);
4346 }
4347 
4348 
4349 #endif // ifndef DOXYGEN
4350 
4351 
4352 
4353 DEAL_II_NAMESPACE_CLOSE
4354 
4355 #endif
MatrixFree::value_type
Number value_type
Definition: matrix_free.h:120
MatrixFree::get_ghost_set
const IndexSet & get_ghost_set(const unsigned int dof_handler_index=0) const
update_JxW_values
Transformed quadrature weights.
Definition: fe_update_flags.h:124
MatrixFree::AdditionalData::mapping_update_flags_inner_faces
UpdateFlags mapping_update_flags_inner_faces
Definition: matrix_free.h:339
MatrixFree::n_physical_cells
unsigned int n_physical_cells() const
MatrixFree::get_cell_iterator
DoFHandler< dim >::cell_iterator get_cell_iterator(const unsigned int macro_cell_number, const unsigned int vector_number, const unsigned int fe_component=0) const
MatrixFree::loop
void loop(const std::function< void(const MatrixFree< dim, Number > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &face_operation, const std::function< void(const MatrixFree< dim, Number > &, 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
MatrixFree::task_info
internal::MatrixFreeFunctions::TaskInfo task_info
Definition: matrix_free.h:1732
MatrixFree::indices_are_initialized
bool indices_are_initialized
Definition: matrix_free.h:1745
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:173
MatrixFree::get_dofs_per_face
unsigned int get_dofs_per_face(const unsigned int fe_component=0, const unsigned int hp_active_fe_index=0) const
MatrixFree::AdditionalData::mapping_update_flags_boundary_faces
UpdateFlags mapping_update_flags_boundary_faces
Definition: matrix_free.h:318
MatrixFree::constraint_pool_begin
const Number * constraint_pool_begin(const unsigned int pool_index) const
MatrixFree::constraint_pool_row_index
std::vector< unsigned int > constraint_pool_row_index
Definition: matrix_free.h:1696
MatrixFree::get_n_q_points_face
unsigned int get_n_q_points_face(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
Threads::ThreadLocalStorage
A class that provides a separate storage location on each thread that accesses the object...
Definition: thread_local_storage.h:72
MatrixFree::n_ghost_cell_batches
unsigned int n_ghost_cell_batches() const
MatrixFree::dimension
static const unsigned int dimension
Definition: matrix_free.h:125
MatrixFree::scratch_pad_non_threadsafe
std::list< std::pair< bool, AlignedVector< Number > > > scratch_pad_non_threadsafe
Definition: matrix_free.h:1768
Utilities::MPI::Partitioner::local_size
unsigned int local_size() const
MatrixFree::dof_handlers
DoFHandlers dof_handlers
Definition: matrix_free.h:1676
MatrixFree::initialize_dof_handlers
void initialize_dof_handlers(const std::vector< const DoFHandler< dim > *> &dof_handlers, const AdditionalData &additional_data)
MatrixFree::get_size_info
const internal::MatrixFreeFunctions::TaskInfo & get_size_info() const
Utilities::MPI::Partitioner::import_from_ghosted_array_finish
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
MatrixFree::n_components
unsigned int n_components() const
Utilities::MPI::Partitioner::import_from_ghosted_array_start
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
MatrixFree::cell_level_index
std::vector< std::pair< unsigned int, unsigned int > > cell_level_index
Definition: matrix_free.h:1716
MatrixFree::get_locally_owned_set
const IndexSet & get_locally_owned_set(const unsigned int dof_handler_index=0) const
MatrixFree::n_active_entries_per_face_batch
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_number) const
MatrixFree::constraint_pool_data
std::vector< Number > constraint_pool_data
Definition: matrix_free.h:1690
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
internal::MatrixFreeFunctions::DoFInfo::vector_partitioner
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition: dof_info.h:469
internal::MatrixFreeFunctions::DoFInfo::vector_zero_range_list_index
std::vector< unsigned int > vector_zero_range_list_index
Definition: dof_info.h:595
MatrixFree::copy_from
void copy_from(const MatrixFree< dim, Number > &matrix_free_base)
MatrixFree::get_faces_by_cells_boundary_id
std::array< types::boundary_id, VectorizedArray< Number >::n_array_elements > get_faces_by_cells_boundary_id(const unsigned int macro_cell, const unsigned int face_number) const
Utilities::MPI::Partitioner::ghost_indices_within_larger_ghost_set
const std::vector< std::pair< unsigned int, unsigned int > > & ghost_indices_within_larger_ghost_set() const
MatrixFree::make_connectivity_graph_faces
void make_connectivity_graph_faces(DynamicSparsityPattern &connectivity)
LinearAlgebra::distributed::Vector::reinit
void reinit(const size_type size, const bool omit_zeroing_entries=false)
StandardExceptions::ExcNotInitialized
static ::ExceptionBase & ExcNotInitialized()
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
MatrixFree::initialize_indices
void initialize_indices(const std::vector< const AffineConstraints< number2 > *> &constraint, const std::vector< IndexSet > &locally_owned_set, const AdditionalData &additional_data)
MatrixFree::get_hp_cell_iterator
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
StandardExceptions::ExcIndexRange
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
MatrixFree::get_dofs_per_cell
unsigned int get_dofs_per_cell(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
MatrixFree::n_boundary_face_batches
unsigned int n_boundary_face_batches() const
MatrixFree::n_inner_face_batches
unsigned int n_inner_face_batches() const
MatrixFree::mapping_is_initialized
bool mapping_is_initialized
Definition: matrix_free.h:1750
MatrixFree::~MatrixFree
~MatrixFree() override=default
MatrixFree::acquire_scratch_data
AlignedVector< VectorizedArray< Number > > * acquire_scratch_data() const
hp::DoFHandler::active_cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:303
MatrixFree::initialize_dof_vector
void initialize_dof_vector(VectorType &vec, const unsigned int dof_handler_index=0) const
MatrixFree::create_cell_subrange_hp
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
MatrixFree::n_components_filled
unsigned int n_components_filled(const unsigned int cell_batch_number) const
LinearAlgebra::distributed::Vector::set_ghost_state
void set_ghost_state(const bool ghosted) const
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
MatrixFree::clear
void clear()
MatrixFree::renumber_dofs
void renumber_dofs(std::vector< types::global_dof_index > &renumbering, const unsigned int dof_handler_index=0)
update_default
No update.
Definition: fe_update_flags.h:67
MatrixFree::get_vector_partitioner
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner(const unsigned int dof_handler_index=0) const
Utilities::MPI::Partitioner::n_ghost_indices
unsigned int n_ghost_indices() const
MatrixFree::at_irregular_cell
bool at_irregular_cell(const unsigned int macro_cell_number) const
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1407
MatrixFree::is_supported
static bool is_supported(const FiniteElement< dim, spacedim > &fe)
UpdateFlags
UpdateFlags
Definition: fe_update_flags.h:64
MatrixFree::print_memory_consumption
void print_memory_consumption(StreamType &out) const
numbers::invalid_boundary_id
const types::boundary_id invalid_boundary_id
Definition: types.h:207
MatrixFree::constraint_pool_end
const Number * constraint_pool_end(const unsigned int pool_index) const
MatrixFree::get_shape_info
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArray< Number > > & 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
MatrixFree::n_constraint_pool_entries
unsigned int n_constraint_pool_entries() const
MatrixFree::AdditionalData::AdditionalData
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:206
MatrixFree::release_scratch_data_non_threadsafe
void release_scratch_data_non_threadsafe(const AlignedVector< Number > *memory) const
MatrixFree::n_macro_cells
unsigned int n_macro_cells() const
MatrixFree::AdditionalData::tasks_parallel_scheme
TasksParallelScheme tasks_parallel_scheme
Definition: matrix_free.h:273
MatrixFree::acquire_scratch_data_non_threadsafe
AlignedVector< Number > * acquire_scratch_data_non_threadsafe() const
MatrixFree::get_boundary_id
types::boundary_id get_boundary_id(const unsigned int macro_face) const
MatrixFree::mapping_initialized
bool mapping_initialized() const
hp
Definition: hp.h:117
MatrixFree::n_ghost_inner_face_batches
unsigned int n_ghost_inner_face_batches() const
MatrixFree::face_info
internal::MatrixFreeFunctions::FaceInfo< VectorizedArray< Number >::n_array_elements > face_info
Definition: matrix_free.h:1740
MatrixFree::shape_info
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< VectorizedArray< Number > > > shape_info
Definition: matrix_free.h:1708
MatrixFree::internal_reinit
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)
MatrixFree::reinit
void reinit(const Mapping< dim > &mapping, const DoFHandlerType &dof_handler, const AffineConstraints< number2 > &constraint, const QuadratureType &quad, const AdditionalData additional_data=AdditionalData())
MatrixFree::AdditionalData::tasks_block_size
unsigned int tasks_block_size
Definition: matrix_free.h:284
MatrixFree::mapping_info
internal::MatrixFreeFunctions::MappingInfo< dim, Number > mapping_info
Definition: matrix_free.h:1702
MatrixFree::cell_level_index_end_local
unsigned int cell_level_index_end_local
Definition: matrix_free.h:1725
MatrixFree::AdditionalData::level_mg_handler
unsigned int level_mg_handler
Definition: matrix_free.h:377
MatrixFree::get_dof_handler
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
update_gradients
Shape function gradients.
Definition: fe_update_flags.h:82
MatrixFree::get_n_q_points
unsigned int get_n_q_points(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
AlignedVector::begin
iterator begin()
MatrixFree::get_face_info
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArray< Number >::n_array_elements > & get_face_info(const unsigned int face_batch_number) const
MatrixFree::n_active_entries_per_cell_batch
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_number) const
MatrixFree::create_cell_subrange_hp_by_index
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
Utilities::MPI::Partitioner::export_to_ghosted_array_start
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
MatrixFree::n_base_elements
unsigned int n_base_elements(const unsigned int dof_handler_index) const
MatrixFree::AdditionalData::mapping_update_flags
UpdateFlags mapping_update_flags
Definition: matrix_free.h:297
Utilities::MPI::Partitioner::export_to_ghosted_array_finish
void export_to_ghosted_array_finish(const ArrayView< Number, MemorySpaceType > &ghost_array, std::vector< MPI_Request > &requests) const
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
MatrixFree::memory_consumption
std::size_t memory_consumption() const
MatrixFree::get_face_quadrature
const Quadrature< dim - 1 > & get_face_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
MatrixFree::scratch_pad
Threads::ThreadLocalStorage< std::list< std::pair< bool, AlignedVector< VectorizedArray< Number > > > > > scratch_pad
Definition: matrix_free.h:1761
internal::MatrixFreeFunctions::DoFInfo::vector_zero_range_list
std::vector< unsigned int > vector_zero_range_list
Definition: dof_info.h:600
Table
Definition: table.h:37
MatrixFree::indices_initialized
bool indices_initialized() const
MatrixFree::get_cell_category
unsigned int get_cell_category(const unsigned int macro_cell) const
MatrixFree::get_face_category
std::pair< unsigned int, unsigned int > get_face_category(const unsigned int macro_face) const
MatrixFree::get_dof_info
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
MatrixFree::get_quadrature
const Quadrature< dim > & get_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
DoFHandler::n_components
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
MatrixFree::dof_info
std::vector< internal::MatrixFreeFunctions::DoFInfo > dof_info
Definition: matrix_free.h:1682
internal::MatrixFreeFunctions::DoFInfo::chunk_size_zero_vector
static const unsigned int chunk_size_zero_vector
Definition: dof_info.h:78
DoFTools::n_components
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
Utilities::MPI::Partitioner::n_import_indices
unsigned int n_import_indices() const
types::boundary_id
unsigned int boundary_id
Definition: types.h:111
MatrixFree::n_cell_batches
unsigned int n_cell_batches() const
MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells
UpdateFlags mapping_update_flags_faces_by_cells
Definition: matrix_free.h:367
MatrixFree::get_task_info
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
MatrixFree::AdditionalData::cell_vectorization_category
std::vector< unsigned int > cell_vectorization_category
Definition: matrix_free.h:446
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
MatrixFree::cell_loop
void cell_loop(const std::function< void(const MatrixFree< dim, Number > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
MatrixFree::print
void print(std::ostream &out) const
MatrixFree::release_scratch_data
void release_scratch_data(const AlignedVector< VectorizedArray< Number >> *memory) const
MatrixFree::get_constrained_dofs
const std::vector< unsigned int > & get_constrained_dofs(const unsigned int dof_handler_index=0) const
Generated by   doxygen 1.8.14