Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
matrix_free.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2011 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii_matrix_free_h
18 #define dealii_matrix_free_h
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/aligned_vector.h>
23 #include <deal.II/base/exceptions.h>
24 #include <deal.II/base/quadrature.h>
25 #include <deal.II/base/template_constraints.h>
26 #include <deal.II/base/thread_local_storage.h>
27 #include <deal.II/base/vectorization.h>
28 
29 #include <deal.II/dofs/dof_handler.h>
30 
31 #include <deal.II/fe/fe.h>
32 #include <deal.II/fe/mapping.h>
33 #include <deal.II/fe/mapping_q1.h>
34 
35 #include <deal.II/grid/grid_tools.h>
36 
37 #include <deal.II/hp/dof_handler.h>
38 #include <deal.II/hp/q_collection.h>
39 
40 #include <deal.II/lac/affine_constraints.h>
41 #include <deal.II/lac/block_vector_base.h>
42 #include <deal.II/lac/la_parallel_vector.h>
43 #include <deal.II/lac/vector_operation.h>
44 
45 #include <deal.II/matrix_free/dof_info.h>
46 #include <deal.II/matrix_free/mapping_info.h>
47 #include <deal.II/matrix_free/shape_info.h>
48 #include <deal.II/matrix_free/task_info.h>
49 #include <deal.II/matrix_free/type_traits.h>
50 
51 #include <cstdlib>
52 #include <limits>
53 #include <list>
54 #include <memory>
55 
56 
57 DEAL_II_NAMESPACE_OPEN
58 
59 
60 
114 template <int dim,
115  typename Number = double,
116  typename VectorizedArrayType = VectorizedArray<Number>>
117 class MatrixFree : public Subscriptor
118 {
119  static_assert(
120  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
121  "Type of Number and of VectorizedArrayType do not match.");
122 
123 public:
128  using value_type = Number;
129  using vectorized_value_type = VectorizedArrayType;
130 
134  static const unsigned int dimension = dim;
135 
182  struct AdditionalData
183  {
189  enum TasksParallelScheme
190  {
194  none = internal::MatrixFreeFunctions::TaskInfo::none,
198  partition_partition =
199  internal::MatrixFreeFunctions::TaskInfo::partition_partition,
203  partition_color =
204  internal::MatrixFreeFunctions::TaskInfo::partition_color,
209  color = internal::MatrixFreeFunctions::TaskInfo::color
210  };
211 
215  AdditionalData(
216  const TasksParallelScheme tasks_parallel_scheme = partition_partition,
217  const unsigned int tasks_block_size = 0,
218  const UpdateFlags mapping_update_flags = update_gradients |
219  update_JxW_values,
220  const UpdateFlags mapping_update_flags_boundary_faces = update_default,
221  const UpdateFlags mapping_update_flags_inner_faces = update_default,
222  const UpdateFlags mapping_update_flags_faces_by_cells = update_default,
223  const unsigned int mg_level = numbers::invalid_unsigned_int,
224  const bool store_plain_indices = true,
225  const bool initialize_indices = true,
226  const bool initialize_mapping = true,
227  const bool overlap_communication_computation = true,
228  const bool hold_all_faces_to_owned_cells = false,
229  const bool cell_vectorization_categories_strict = false)
230  : tasks_parallel_scheme(tasks_parallel_scheme)
231  , tasks_block_size(tasks_block_size)
232  , mapping_update_flags(mapping_update_flags)
233  , mapping_update_flags_boundary_faces(mapping_update_flags_boundary_faces)
234  , mapping_update_flags_inner_faces(mapping_update_flags_inner_faces)
235  , mapping_update_flags_faces_by_cells(mapping_update_flags_faces_by_cells)
236  , mg_level(mg_level)
237  , level_mg_handler(this->mg_level)
238  , store_plain_indices(store_plain_indices)
239  , initialize_indices(initialize_indices)
240  , initialize_mapping(initialize_mapping)
241  , overlap_communication_computation(overlap_communication_computation)
242  , hold_all_faces_to_owned_cells(hold_all_faces_to_owned_cells)
243  , cell_vectorization_categories_strict(
244  cell_vectorization_categories_strict)
245  {}
246 
250  AdditionalData(const AdditionalData &other)
251  : tasks_parallel_scheme(other.tasks_parallel_scheme)
252  , tasks_block_size(other.tasks_block_size)
253  , mapping_update_flags(other.mapping_update_flags)
254  , mapping_update_flags_boundary_faces(
255  other.mapping_update_flags_boundary_faces)
256  , mapping_update_flags_inner_faces(other.mapping_update_flags_inner_faces)
257  , mapping_update_flags_faces_by_cells(
258  other.mapping_update_flags_faces_by_cells)
259  , mg_level(other.mg_level)
260  , level_mg_handler(this->mg_level)
261  , store_plain_indices(other.store_plain_indices)
262  , initialize_indices(other.initialize_indices)
263  , initialize_mapping(other.initialize_mapping)
264  , overlap_communication_computation(
265  other.overlap_communication_computation)
266  , hold_all_faces_to_owned_cells(other.hold_all_faces_to_owned_cells)
267  , cell_vectorization_category(other.cell_vectorization_category)
268  , cell_vectorization_categories_strict(
269  other.cell_vectorization_categories_strict)
270  {}
271 
275  AdditionalData &
276  operator=(const AdditionalData &other)
277  {
278  tasks_parallel_scheme = other.tasks_parallel_scheme;
279  tasks_block_size = other.tasks_block_size;
280  mapping_update_flags = other.mapping_update_flags;
281  mapping_update_flags_boundary_faces =
282  other.mapping_update_flags_boundary_faces;
283  mapping_update_flags_inner_faces = other.mapping_update_flags_inner_faces;
284  mapping_update_flags_faces_by_cells =
285  other.mapping_update_flags_faces_by_cells;
286  mg_level = other.mg_level;
287  store_plain_indices = other.store_plain_indices;
288  initialize_indices = other.initialize_indices;
289  initialize_mapping = other.initialize_mapping;
290  overlap_communication_computation =
291  other.overlap_communication_computation;
292  hold_all_faces_to_owned_cells = other.hold_all_faces_to_owned_cells;
293  cell_vectorization_category = other.cell_vectorization_category;
294  cell_vectorization_categories_strict =
295  other.cell_vectorization_categories_strict;
296 
297  return *this;
298  }
299 
336  TasksParallelScheme tasks_parallel_scheme;
337 
347  unsigned int tasks_block_size;
348 
360  UpdateFlags mapping_update_flags;
361 
381  UpdateFlags mapping_update_flags_boundary_faces;
382 
402  UpdateFlags mapping_update_flags_inner_faces;
403 
430  UpdateFlags mapping_update_flags_faces_by_cells;
431 
440  unsigned int mg_level;
441 
447  DEAL_II_DEPRECATED unsigned int &level_mg_handler;
448 
455  bool store_plain_indices;
456 
466  bool initialize_indices;
467 
475  bool initialize_mapping;
476 
489  bool overlap_communication_computation;
490 
498  bool hold_all_faces_to_owned_cells;
499 
516  std::vector<unsigned int> cell_vectorization_category;
517 
527  bool cell_vectorization_categories_strict;
528  };
529 
537  MatrixFree();
538 
542  MatrixFree(const MatrixFree<dim, Number, VectorizedArrayType> &other);
543 
547  ~MatrixFree() override = default;
548 
561  template <typename DoFHandlerType, typename QuadratureType, typename number2>
562  void
563  reinit(const Mapping<dim> & mapping,
564  const DoFHandlerType & dof_handler,
565  const AffineConstraints<number2> &constraint,
566  const QuadratureType & quad,
567  const AdditionalData & additional_data = AdditionalData());
568 
573  template <typename DoFHandlerType, typename QuadratureType, typename number2>
574  void
575  reinit(const DoFHandlerType & dof_handler,
576  const AffineConstraints<number2> &constraint,
577  const QuadratureType & quad,
578  const AdditionalData & additional_data = AdditionalData());
579 
587  template <typename DoFHandlerType, typename QuadratureType, typename number2>
588  DEAL_II_DEPRECATED void
589  reinit(const Mapping<dim> & mapping,
590  const DoFHandlerType & dof_handler,
591  const AffineConstraints<number2> &constraint,
592  const IndexSet & locally_owned_dofs,
593  const QuadratureType & quad,
594  const AdditionalData & additional_data = AdditionalData());
595 
617  template <typename DoFHandlerType, typename QuadratureType, typename number2>
618  void
619  reinit(const Mapping<dim> & mapping,
620  const std::vector<const DoFHandlerType *> & dof_handler,
621  const std::vector<const AffineConstraints<number2> *> &constraint,
622  const std::vector<QuadratureType> & quad,
623  const AdditionalData &additional_data = AdditionalData());
624 
629  template <typename DoFHandlerType, typename QuadratureType, typename number2>
630  void
631  reinit(const std::vector<const DoFHandlerType *> & dof_handler,
632  const std::vector<const AffineConstraints<number2> *> &constraint,
633  const std::vector<QuadratureType> & quad,
634  const AdditionalData &additional_data = AdditionalData());
635 
643  template <typename DoFHandlerType, typename QuadratureType, typename number2>
644  DEAL_II_DEPRECATED void
645  reinit(const Mapping<dim> & mapping,
646  const std::vector<const DoFHandlerType *> & dof_handler,
647  const std::vector<const AffineConstraints<number2> *> &constraint,
648  const std::vector<IndexSet> & locally_owned_set,
649  const std::vector<QuadratureType> &quad,
650  const AdditionalData & additional_data = AdditionalData());
651 
659  template <typename DoFHandlerType, typename QuadratureType, typename number2>
660  void
661  reinit(const Mapping<dim> & mapping,
662  const std::vector<const DoFHandlerType *> & dof_handler,
663  const std::vector<const AffineConstraints<number2> *> &constraint,
664  const QuadratureType & quad,
665  const AdditionalData &additional_data = AdditionalData());
666 
671  template <typename DoFHandlerType, typename QuadratureType, typename number2>
672  void
673  reinit(const std::vector<const DoFHandlerType *> & dof_handler,
674  const std::vector<const AffineConstraints<number2> *> &constraint,
675  const QuadratureType & quad,
676  const AdditionalData &additional_data = AdditionalData());
677 
683  void
684  copy_from(
685  const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free_base);
686 
696  void
697  update_mapping(const Mapping<dim> &mapping);
698 
703  void
704  clear();
705 
707 
719  enum class DataAccessOnFaces
720  {
727  none,
728 
738  values,
739 
748  values_all_faces,
749 
760  gradients,
761 
770  gradients_all_faces,
771 
778  unspecified
779  };
780 
825  template <typename OutVector, typename InVector>
826  void
827  cell_loop(const std::function<void(
828  const MatrixFree<dim, Number, VectorizedArrayType> &,
829  OutVector &,
830  const InVector &,
831  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
832  OutVector & dst,
833  const InVector & src,
834  const bool zero_dst_vector = false) const;
835 
882  template <typename CLASS, typename OutVector, typename InVector>
883  void
884  cell_loop(void (CLASS::*cell_operation)(
885  const MatrixFree &,
886  OutVector &,
887  const InVector &,
888  const std::pair<unsigned int, unsigned int> &) const,
889  const CLASS * owning_class,
890  OutVector & dst,
891  const InVector &src,
892  const bool zero_dst_vector = false) const;
893 
897  template <typename CLASS, typename OutVector, typename InVector>
898  void
899  cell_loop(void (CLASS::*cell_operation)(
900  const MatrixFree &,
901  OutVector &,
902  const InVector &,
903  const std::pair<unsigned int, unsigned int> &),
904  CLASS * owning_class,
905  OutVector & dst,
906  const InVector &src,
907  const bool zero_dst_vector = false) const;
908 
993  template <typename CLASS, typename OutVector, typename InVector>
994  void
995  cell_loop(void (CLASS::*cell_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 std::function<void(const unsigned int, const unsigned int)>
1004  &operation_before_loop,
1005  const std::function<void(const unsigned int, const unsigned int)>
1006  & operation_after_loop,
1007  const unsigned int dof_handler_index_pre_post = 0) const;
1008 
1012  template <typename CLASS, typename OutVector, typename InVector>
1013  void
1014  cell_loop(void (CLASS::*cell_operation)(
1015  const MatrixFree &,
1016  OutVector &,
1017  const InVector &,
1018  const std::pair<unsigned int, unsigned int> &),
1019  CLASS * owning_class,
1020  OutVector & dst,
1021  const InVector &src,
1022  const std::function<void(const unsigned int, const unsigned int)>
1023  &operation_before_loop,
1024  const std::function<void(const unsigned int, const unsigned int)>
1025  & operation_after_loop,
1026  const unsigned int dof_handler_index_pre_post = 0) const;
1027 
1032  template <typename OutVector, typename InVector>
1033  void
1034  cell_loop(const std::function<void(
1035  const MatrixFree<dim, Number, VectorizedArrayType> &,
1036  OutVector &,
1037  const InVector &,
1038  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1039  OutVector & dst,
1040  const InVector & src,
1041  const std::function<void(const unsigned int, const unsigned int)>
1042  &operation_before_loop,
1043  const std::function<void(const unsigned int, const unsigned int)>
1044  & operation_after_loop,
1045  const unsigned int dof_handler_index_pre_post = 0) const;
1046 
1122  template <typename OutVector, typename InVector>
1123  void
1124  loop(const std::function<
1125  void(const MatrixFree<dim, Number, VectorizedArrayType> &,
1126  OutVector &,
1127  const InVector &,
1128  const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1129  const std::function<
1130  void(const MatrixFree<dim, Number, VectorizedArrayType> &,
1131  OutVector &,
1132  const InVector &,
1133  const std::pair<unsigned int, unsigned int> &)> &face_operation,
1134  const std::function<void(
1135  const MatrixFree<dim, Number, VectorizedArrayType> &,
1136  OutVector &,
1137  const InVector &,
1138  const std::pair<unsigned int, unsigned int> &)> &boundary_operation,
1139  OutVector & dst,
1140  const InVector & src,
1141  const bool zero_dst_vector = false,
1142  const DataAccessOnFaces dst_vector_face_access =
1143  DataAccessOnFaces::unspecified,
1144  const DataAccessOnFaces src_vector_face_access =
1145  DataAccessOnFaces::unspecified) const;
1146 
1232  template <typename CLASS, typename OutVector, typename InVector>
1233  void
1234  loop(
1235  void (CLASS::*cell_operation)(const MatrixFree &,
1236  OutVector &,
1237  const InVector &,
1238  const std::pair<unsigned int, unsigned int> &)
1239  const,
1240  void (CLASS::*face_operation)(const MatrixFree &,
1241  OutVector &,
1242  const InVector &,
1243  const std::pair<unsigned int, unsigned int> &)
1244  const,
1245  void (CLASS::*boundary_operation)(
1246  const MatrixFree &,
1247  OutVector &,
1248  const InVector &,
1249  const std::pair<unsigned int, unsigned int> &) const,
1250  const CLASS * owning_class,
1251  OutVector & dst,
1252  const InVector & src,
1253  const bool zero_dst_vector = false,
1254  const DataAccessOnFaces dst_vector_face_access =
1255  DataAccessOnFaces::unspecified,
1256  const DataAccessOnFaces src_vector_face_access =
1257  DataAccessOnFaces::unspecified) const;
1258 
1262  template <typename CLASS, typename OutVector, typename InVector>
1263  void
1264  loop(void (CLASS::*cell_operation)(
1265  const MatrixFree &,
1266  OutVector &,
1267  const InVector &,
1268  const std::pair<unsigned int, unsigned int> &),
1269  void (CLASS::*face_operation)(
1270  const MatrixFree &,
1271  OutVector &,
1272  const InVector &,
1273  const std::pair<unsigned int, unsigned int> &),
1274  void (CLASS::*boundary_operation)(
1275  const MatrixFree &,
1276  OutVector &,
1277  const InVector &,
1278  const std::pair<unsigned int, unsigned int> &),
1279  CLASS * owning_class,
1280  OutVector & dst,
1281  const InVector & src,
1282  const bool zero_dst_vector = false,
1283  const DataAccessOnFaces dst_vector_face_access =
1284  DataAccessOnFaces::unspecified,
1285  const DataAccessOnFaces src_vector_face_access =
1286  DataAccessOnFaces::unspecified) const;
1287 
1352  template <typename CLASS, typename OutVector, typename InVector>
1353  void
1354  loop_cell_centric(void (CLASS::*cell_operation)(
1355  const MatrixFree &,
1356  OutVector &,
1357  const InVector &,
1358  const std::pair<unsigned int, unsigned int> &) const,
1359  const CLASS * owning_class,
1360  OutVector & dst,
1361  const InVector & src,
1362  const bool zero_dst_vector = false,
1363  const DataAccessOnFaces src_vector_face_access =
1364  DataAccessOnFaces::unspecified) const;
1365 
1369  template <typename CLASS, typename OutVector, typename InVector>
1370  void
1371  loop_cell_centric(void (CLASS::*cell_operation)(
1372  const MatrixFree &,
1373  OutVector &,
1374  const InVector &,
1375  const std::pair<unsigned int, unsigned int> &),
1376  CLASS * owning_class,
1377  OutVector & dst,
1378  const InVector & src,
1379  const bool zero_dst_vector = false,
1380  const DataAccessOnFaces src_vector_face_access =
1381  DataAccessOnFaces::unspecified) const;
1382 
1386  template <typename OutVector, typename InVector>
1387  void
1388  loop_cell_centric(
1389  const std::function<void(const MatrixFree &,
1390  OutVector &,
1391  const InVector &,
1392  const std::pair<unsigned int, unsigned int> &)>
1393  & cell_operation,
1394  OutVector & dst,
1395  const InVector & src,
1396  const bool zero_dst_vector = false,
1397  const DataAccessOnFaces src_vector_face_access =
1398  DataAccessOnFaces::unspecified) const;
1399 
1407  std::pair<unsigned int, unsigned int>
1408  create_cell_subrange_hp(const std::pair<unsigned int, unsigned int> &range,
1409  const unsigned int fe_degree,
1410  const unsigned int dof_handler_index = 0) const;
1411 
1418  std::pair<unsigned int, unsigned int>
1419  create_cell_subrange_hp_by_index(
1420  const std::pair<unsigned int, unsigned int> &range,
1421  const unsigned int fe_index,
1422  const unsigned int dof_handler_index = 0) const;
1423 
1425 
1450  template <typename VectorType>
1451  void
1452  initialize_dof_vector(VectorType & vec,
1453  const unsigned int dof_handler_index = 0) const;
1454 
1475  template <typename Number2>
1476  void
1477  initialize_dof_vector(LinearAlgebra::distributed::Vector<Number2> &vec,
1478  const unsigned int dof_handler_index = 0) const;
1479 
1490  const std::shared_ptr<const Utilities::MPI::Partitioner> &
1491  get_vector_partitioner(const unsigned int dof_handler_index = 0) const;
1492 
1496  const IndexSet &
1497  get_locally_owned_set(const unsigned int dof_handler_index = 0) const;
1498 
1502  const IndexSet &
1503  get_ghost_set(const unsigned int dof_handler_index = 0) const;
1504 
1514  const std::vector<unsigned int> &
1515  get_constrained_dofs(const unsigned int dof_handler_index = 0) const;
1516 
1527  void
1528  renumber_dofs(std::vector<types::global_dof_index> &renumbering,
1529  const unsigned int dof_handler_index = 0);
1530 
1532 
1540  template <int spacedim>
1541  static bool
1542  is_supported(const FiniteElement<dim, spacedim> &fe);
1543 
1547  unsigned int
1548  n_components() const;
1549 
1554  unsigned int
1555  n_base_elements(const unsigned int dof_handler_index) const;
1556 
1564  unsigned int
1565  n_physical_cells() const;
1566 
1576  unsigned int
1577  n_macro_cells() const;
1578 
1588  unsigned int
1589  n_cell_batches() const;
1590 
1597  unsigned int
1598  n_ghost_cell_batches() const;
1599 
1607  unsigned int
1608  n_inner_face_batches() const;
1609 
1618  unsigned int
1619  n_boundary_face_batches() const;
1620 
1625  unsigned int
1626  n_ghost_inner_face_batches() const;
1627 
1634  types::boundary_id
1635  get_boundary_id(const unsigned int macro_face) const;
1636 
1641  std::array<types::boundary_id, VectorizedArrayType::size()>
1642  get_faces_by_cells_boundary_id(const unsigned int macro_cell,
1643  const unsigned int face_number) const;
1644 
1653  template <typename DoFHandlerType = DoFHandler<dim>>
1654  const DoFHandlerType &
1655  get_dof_handler(const unsigned int dof_handler_index = 0) const;
1656 
1667  typename DoFHandler<dim>::cell_iterator
1668  get_cell_iterator(const unsigned int macro_cell_number,
1669  const unsigned int vector_number,
1670  const unsigned int dof_handler_index = 0) const;
1671 
1677  std::pair<int, int>
1678  get_cell_level_and_index(const unsigned int macro_cell_number,
1679  const unsigned int vector_number) const;
1680 
1693  std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>
1694  get_face_iterator(const unsigned int face_batch_number,
1695  const unsigned int vector_number,
1696  const bool interior = true,
1697  const unsigned int fe_component = 0) const;
1698 
1710  typename hp::DoFHandler<dim>::active_cell_iterator
1711  get_hp_cell_iterator(const unsigned int macro_cell_number,
1712  const unsigned int vector_number,
1713  const unsigned int dof_handler_index = 0) const;
1714 
1726  bool
1727  at_irregular_cell(const unsigned int macro_cell_number) const;
1728 
1737  unsigned int
1738  n_components_filled(const unsigned int cell_batch_number) const;
1739 
1748  unsigned int
1749  n_active_entries_per_cell_batch(const unsigned int cell_batch_number) const;
1750 
1760  unsigned int
1761  n_active_entries_per_face_batch(const unsigned int face_batch_number) const;
1762 
1766  unsigned int
1767  get_dofs_per_cell(const unsigned int dof_handler_index = 0,
1768  const unsigned int hp_active_fe_index = 0) const;
1769 
1773  unsigned int
1774  get_n_q_points(const unsigned int quad_index = 0,
1775  const unsigned int hp_active_fe_index = 0) const;
1776 
1781  unsigned int
1782  get_dofs_per_face(const unsigned int dof_handler_index = 0,
1783  const unsigned int hp_active_fe_index = 0) const;
1784 
1789  unsigned int
1790  get_n_q_points_face(const unsigned int quad_index = 0,
1791  const unsigned int hp_active_fe_index = 0) const;
1792 
1796  const Quadrature<dim> &
1797  get_quadrature(const unsigned int quad_index = 0,
1798  const unsigned int hp_active_fe_index = 0) const;
1799 
1803  const Quadrature<dim - 1> &
1804  get_face_quadrature(const unsigned int quad_index = 0,
1805  const unsigned int hp_active_fe_index = 0) const;
1806 
1813  unsigned int
1814  get_cell_category(const unsigned int macro_cell) const;
1815 
1820  std::pair<unsigned int, unsigned int>
1821  get_face_category(const unsigned int macro_face) const;
1822 
1826  bool
1827  indices_initialized() const;
1828 
1833  bool
1834  mapping_initialized() const;
1835 
1840  unsigned int
1841  get_mg_level() const;
1842 
1847  std::size_t
1848  memory_consumption() const;
1849 
1854  template <typename StreamType>
1855  void
1856  print_memory_consumption(StreamType &out) const;
1857 
1862  void
1863  print(std::ostream &out) const;
1864 
1866 
1876  const internal::MatrixFreeFunctions::TaskInfo &
1877  get_task_info() const;
1878 
1879  /*
1880  * Return geometry-dependent information on the cells.
1881  */
1882  const internal::MatrixFreeFunctions::
1883  MappingInfo<dim, Number, VectorizedArrayType> &
1884  get_mapping_info() const;
1885 
1889  const internal::MatrixFreeFunctions::DoFInfo &
1890  get_dof_info(const unsigned int dof_handler_index_component = 0) const;
1891 
1895  unsigned int
1896  n_constraint_pool_entries() const;
1897 
1902  const Number *
1903  constraint_pool_begin(const unsigned int pool_index) const;
1904 
1910  const Number *
1911  constraint_pool_end(const unsigned int pool_index) const;
1912 
1916  const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
1917  get_shape_info(const unsigned int dof_handler_index_component = 0,
1918  const unsigned int quad_index = 0,
1919  const unsigned int fe_base_element = 0,
1920  const unsigned int hp_active_fe_index = 0,
1921  const unsigned int hp_active_quad_index = 0) const;
1922 
1926  const internal::MatrixFreeFunctions::FaceToCellTopology<
1927  VectorizedArrayType::size()> &
1928  get_face_info(const unsigned int face_batch_number) const;
1929 
1930 
1936  const Table<3, unsigned int> &
1937  get_cell_and_face_to_plain_faces() const;
1938 
1952  AlignedVector<VectorizedArrayType> *
1953  acquire_scratch_data() const;
1954 
1958  void
1959  release_scratch_data(const AlignedVector<VectorizedArrayType> *memory) const;
1960 
1970  AlignedVector<Number> *
1971  acquire_scratch_data_non_threadsafe() const;
1972 
1976  void
1977  release_scratch_data_non_threadsafe(
1978  const AlignedVector<Number> *memory) const;
1979 
1981 
1982 private:
1987  template <typename number2, template <int, int> class DoFHandlerType>
1988  void
1989  internal_reinit(
1990  const Mapping<dim> & mapping,
1991  const std::vector<const DoFHandlerType<dim, dim> *> & dof_handler,
1992  const std::vector<const AffineConstraints<number2> *> &constraint,
1993  const std::vector<IndexSet> & locally_owned_set,
1994  const std::vector<hp::QCollection<1>> & quad,
1995  const AdditionalData & additional_data);
1996 
2003  template <typename number2>
2004  void
2005  initialize_indices(
2006  const std::vector<const AffineConstraints<number2> *> &constraint,
2007  const std::vector<IndexSet> & locally_owned_set,
2008  const AdditionalData & additional_data);
2009 
2013  void
2014  initialize_dof_handlers(
2015  const std::vector<const DoFHandler<dim> *> &dof_handlers,
2016  const AdditionalData & additional_data);
2017 
2021  void
2022  initialize_dof_handlers(
2023  const std::vector<const hp::DoFHandler<dim> *> &dof_handlers,
2024  const AdditionalData & additional_data);
2025 
2030  void
2031  make_connectivity_graph_faces(DynamicSparsityPattern &connectivity);
2032 
2038  struct DoFHandlers
2039  {
2040  DoFHandlers()
2041  : active_dof_handler(usual)
2042  , n_dof_handlers(0)
2043  {}
2044 
2045  std::vector<SmartPointer<const DoFHandler<dim>>> dof_handler;
2046  std::vector<SmartPointer<const hp::DoFHandler<dim>>> hp_dof_handler;
2047  enum ActiveDoFHandler
2048  {
2052  usual,
2056  hp
2057  } active_dof_handler;
2058  unsigned int n_dof_handlers;
2059  };
2060 
2064  DoFHandlers dof_handlers;
2065 
2070  std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
2071 
2078  std::vector<Number> constraint_pool_data;
2079 
2084  std::vector<unsigned int> constraint_pool_row_index;
2085 
2090  internal::MatrixFreeFunctions::MappingInfo<dim, Number, VectorizedArrayType>
2091  mapping_info;
2092 
2096  Table<4, internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>>
2097  shape_info;
2098 
2105  std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
2106 
2107 
2114  unsigned int cell_level_index_end_local;
2115 
2121  internal::MatrixFreeFunctions::TaskInfo task_info;
2122 
2127  internal::MatrixFreeFunctions::FaceInfo<VectorizedArrayType::size()>
2128  face_info;
2129 
2133  bool indices_are_initialized;
2134 
2138  bool mapping_is_initialized;
2139 
2147  mutable Threads::ThreadLocalStorage<
2148  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2149  scratch_pad;
2150 
2155  mutable std::list<std::pair<bool, AlignedVector<Number>>>
2156  scratch_pad_non_threadsafe;
2157 
2161  unsigned int mg_level;
2162 };
2163 
2164 
2165 
2166 /*----------------------- Inline functions ----------------------------------*/
2167 
2168 #ifndef DOXYGEN
2169 
2170 
2171 
2172 template <int dim, typename Number, typename VectorizedArrayType>
2173 template <typename VectorType>
2174 inline void
2175 MatrixFree<dim, Number, VectorizedArrayType>::initialize_dof_vector(
2176  VectorType & vec,
2177  const unsigned int comp) const
2178 {
2179  AssertIndexRange(comp, n_components());
2180  vec.reinit(dof_info[comp].vector_partitioner->size());
2181 }
2182 
2183 
2184 
2185 template <int dim, typename Number, typename VectorizedArrayType>
2186 template <typename Number2>
2187 inline void
2188 MatrixFree<dim, Number, VectorizedArrayType>::initialize_dof_vector(
2189  LinearAlgebra::distributed::Vector<Number2> &vec,
2190  const unsigned int comp) const
2191 {
2192  AssertIndexRange(comp, n_components());
2193  vec.reinit(dof_info[comp].vector_partitioner);
2194 }
2195 
2196 
2197 
2198 template <int dim, typename Number, typename VectorizedArrayType>
2199 inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2200 MatrixFree<dim, Number, VectorizedArrayType>::get_vector_partitioner(
2201  const unsigned int comp) const
2202 {
2203  AssertIndexRange(comp, n_components());
2204  return dof_info[comp].vector_partitioner;
2205 }
2206 
2207 
2208 
2209 template <int dim, typename Number, typename VectorizedArrayType>
2210 inline const std::vector<unsigned int> &
2211 MatrixFree<dim, Number, VectorizedArrayType>::get_constrained_dofs(
2212  const unsigned int comp) const
2213 {
2214  AssertIndexRange(comp, n_components());
2215  return dof_info[comp].constrained_dofs;
2216 }
2217 
2218 
2219 
2220 template <int dim, typename Number, typename VectorizedArrayType>
2221 inline unsigned int
2222 MatrixFree<dim, Number, VectorizedArrayType>::n_components() const
2223 {
2224  AssertDimension(dof_handlers.n_dof_handlers, dof_info.size());
2225  return dof_handlers.n_dof_handlers;
2226 }
2227 
2228 
2229 
2230 template <int dim, typename Number, typename VectorizedArrayType>
2231 inline unsigned int
2232 MatrixFree<dim, Number, VectorizedArrayType>::n_base_elements(
2233  const unsigned int dof_no) const
2234 {
2235  AssertDimension(dof_handlers.n_dof_handlers, dof_info.size());
2236  AssertIndexRange(dof_no, dof_handlers.n_dof_handlers);
2237  return dof_handlers.dof_handler[dof_no]->get_fe().n_base_elements();
2238 }
2239 
2240 
2241 
2242 template <int dim, typename Number, typename VectorizedArrayType>
2243 inline const internal::MatrixFreeFunctions::TaskInfo &
2244 MatrixFree<dim, Number, VectorizedArrayType>::get_task_info() const
2245 {
2246  return task_info;
2247 }
2248 
2249 
2250 
2251 template <int dim, typename Number, typename VectorizedArrayType>
2252 inline unsigned int
2253 MatrixFree<dim, Number, VectorizedArrayType>::n_macro_cells() const
2254 {
2255  return *(task_info.cell_partition_data.end() - 2);
2256 }
2257 
2258 
2259 
2260 template <int dim, typename Number, typename VectorizedArrayType>
2261 inline unsigned int
2262 MatrixFree<dim, Number, VectorizedArrayType>::n_physical_cells() const
2263 {
2264  return task_info.n_active_cells;
2265 }
2266 
2267 
2268 
2269 template <int dim, typename Number, typename VectorizedArrayType>
2270 inline unsigned int
2271 MatrixFree<dim, Number, VectorizedArrayType>::n_cell_batches() const
2272 {
2273  return *(task_info.cell_partition_data.end() - 2);
2274 }
2275 
2276 
2277 
2278 template <int dim, typename Number, typename VectorizedArrayType>
2279 inline unsigned int
2280 MatrixFree<dim, Number, VectorizedArrayType>::n_ghost_cell_batches() const
2281 {
2282  return *(task_info.cell_partition_data.end() - 1) -
2283  *(task_info.cell_partition_data.end() - 2);
2284 }
2285 
2286 
2287 
2288 template <int dim, typename Number, typename VectorizedArrayType>
2289 inline unsigned int
2290 MatrixFree<dim, Number, VectorizedArrayType>::n_inner_face_batches() const
2291 {
2292  if (task_info.face_partition_data.size() == 0)
2293  return 0;
2294  return task_info.face_partition_data.back();
2295 }
2296 
2297 
2298 
2299 template <int dim, typename Number, typename VectorizedArrayType>
2300 inline unsigned int
2301 MatrixFree<dim, Number, VectorizedArrayType>::n_boundary_face_batches() const
2302 {
2303  if (task_info.face_partition_data.size() == 0)
2304  return 0;
2305  return task_info.boundary_partition_data.back() -
2306  task_info.face_partition_data.back();
2307 }
2308 
2309 
2310 
2311 template <int dim, typename Number, typename VectorizedArrayType>
2312 inline unsigned int
2313 MatrixFree<dim, Number, VectorizedArrayType>::n_ghost_inner_face_batches() const
2314 {
2315  if (task_info.face_partition_data.size() == 0)
2316  return 0;
2317  return face_info.faces.size() - task_info.boundary_partition_data.back();
2318 }
2319 
2320 
2321 
2322 template <int dim, typename Number, typename VectorizedArrayType>
2323 inline types::boundary_id
2324 MatrixFree<dim, Number, VectorizedArrayType>::get_boundary_id(
2325  const unsigned int macro_face) const
2326 {
2327  Assert(macro_face >= task_info.boundary_partition_data[0] &&
2328  macro_face < task_info.boundary_partition_data.back(),
2329  ExcIndexRange(macro_face,
2330  task_info.boundary_partition_data[0],
2331  task_info.boundary_partition_data.back()));
2332  return types::boundary_id(face_info.faces[macro_face].exterior_face_no);
2333 }
2334 
2335 
2336 
2337 template <int dim, typename Number, typename VectorizedArrayType>
2338 inline std::array<types::boundary_id, VectorizedArrayType::size()>
2339 MatrixFree<dim, Number, VectorizedArrayType>::get_faces_by_cells_boundary_id(
2340  const unsigned int macro_cell,
2341  const unsigned int face_number) const
2342 {
2343  AssertIndexRange(macro_cell, n_macro_cells());
2344  AssertIndexRange(face_number, GeometryInfo<dim>::faces_per_cell);
2345  Assert(face_info.cell_and_face_boundary_id.size(0) >= n_macro_cells(),
2346  ExcNotInitialized());
2347  std::array<types::boundary_id, VectorizedArrayType::size()> result;
2348  result.fill(numbers::invalid_boundary_id);
2349  for (unsigned int v = 0; v < n_active_entries_per_cell_batch(macro_cell); ++v)
2350  result[v] = face_info.cell_and_face_boundary_id(macro_cell, face_number, v);
2351  return result;
2352 }
2353 
2354 
2355 
2356 template <int dim, typename Number, typename VectorizedArrayType>
2357 inline const internal::MatrixFreeFunctions::
2358  MappingInfo<dim, Number, VectorizedArrayType> &
2359  MatrixFree<dim, Number, VectorizedArrayType>::get_mapping_info() const
2360 {
2361  return mapping_info;
2362 }
2363 
2364 
2365 
2366 template <int dim, typename Number, typename VectorizedArrayType>
2367 inline const internal::MatrixFreeFunctions::DoFInfo &
2368 MatrixFree<dim, Number, VectorizedArrayType>::get_dof_info(
2369  const unsigned int dof_index) const
2370 {
2371  AssertIndexRange(dof_index, n_components());
2372  return dof_info[dof_index];
2373 }
2374 
2375 
2376 
2377 template <int dim, typename Number, typename VectorizedArrayType>
2378 inline unsigned int
2379 MatrixFree<dim, Number, VectorizedArrayType>::n_constraint_pool_entries() const
2380 {
2381  return constraint_pool_row_index.size() - 1;
2382 }
2383 
2384 
2385 
2386 template <int dim, typename Number, typename VectorizedArrayType>
2387 inline const Number *
2388 MatrixFree<dim, Number, VectorizedArrayType>::constraint_pool_begin(
2389  const unsigned int row) const
2390 {
2391  AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2392  return constraint_pool_data.empty() ?
2393  nullptr :
2394  constraint_pool_data.data() + constraint_pool_row_index[row];
2395 }
2396 
2397 
2398 
2399 template <int dim, typename Number, typename VectorizedArrayType>
2400 inline const Number *
2401 MatrixFree<dim, Number, VectorizedArrayType>::constraint_pool_end(
2402  const unsigned int row) const
2403 {
2404  AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2405  return constraint_pool_data.empty() ?
2406  nullptr :
2407  constraint_pool_data.data() + constraint_pool_row_index[row + 1];
2408 }
2409 
2410 
2411 
2412 template <int dim, typename Number, typename VectorizedArrayType>
2413 inline std::pair<unsigned int, unsigned int>
2414 MatrixFree<dim, Number, VectorizedArrayType>::create_cell_subrange_hp(
2415  const std::pair<unsigned int, unsigned int> &range,
2416  const unsigned int degree,
2417  const unsigned int dof_handler_component) const
2418 {
2419  if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2420  {
2421  AssertDimension(
2422  dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2423  AssertDimension(
2424  dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2425  if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2426  return range;
2427  else
2428  return {range.second, range.second};
2429  }
2430 
2431  const unsigned int fe_index =
2432  dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2433  if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2434  return {range.second, range.second};
2435  else
2436  return create_cell_subrange_hp_by_index(range,
2437  fe_index,
2438  dof_handler_component);
2439 }
2440 
2441 
2442 
2443 template <int dim, typename Number, typename VectorizedArrayType>
2444 inline bool
2445 MatrixFree<dim, Number, VectorizedArrayType>::at_irregular_cell(
2446  const unsigned int macro_cell) const
2447 {
2448  AssertIndexRange(macro_cell, task_info.cell_partition_data.back());
2449  return VectorizedArrayType::size() > 1 &&
2450  cell_level_index[(macro_cell + 1) * VectorizedArrayType::size() - 1] ==
2451  cell_level_index[(macro_cell + 1) * VectorizedArrayType::size() - 2];
2452 }
2453 
2454 
2455 
2456 template <int dim, typename Number, typename VectorizedArrayType>
2457 inline unsigned int
2458 MatrixFree<dim, Number, VectorizedArrayType>::n_components_filled(
2459  const unsigned int cell_batch_number) const
2460 {
2461  return n_active_entries_per_cell_batch(cell_batch_number);
2462 }
2463 
2464 
2465 
2466 template <int dim, typename Number, typename VectorizedArrayType>
2467 inline unsigned int
2468 MatrixFree<dim, Number, VectorizedArrayType>::n_active_entries_per_cell_batch(
2469  const unsigned int cell_batch_number) const
2470 {
2471  AssertIndexRange(cell_batch_number, task_info.cell_partition_data.back());
2472  unsigned int n_lanes = VectorizedArrayType::size();
2473  while (n_lanes > 1 &&
2474  cell_level_index[cell_batch_number * VectorizedArrayType::size() +
2475  n_lanes - 1] ==
2476  cell_level_index[cell_batch_number * VectorizedArrayType::size() +
2477  n_lanes - 2])
2478  --n_lanes;
2479  AssertIndexRange(n_lanes - 1, VectorizedArrayType::size());
2480  return n_lanes;
2481 }
2482 
2483 
2484 
2485 template <int dim, typename Number, typename VectorizedArrayType>
2486 inline unsigned int
2487 MatrixFree<dim, Number, VectorizedArrayType>::n_active_entries_per_face_batch(
2488  const unsigned int face_batch_number) const
2489 {
2490  AssertIndexRange(face_batch_number, face_info.faces.size());
2491  unsigned int n_lanes = VectorizedArrayType::size();
2492  while (n_lanes > 1 &&
2493  face_info.faces[face_batch_number].cells_interior[n_lanes - 1] ==
2494  numbers::invalid_unsigned_int)
2495  --n_lanes;
2496  AssertIndexRange(n_lanes - 1, VectorizedArrayType::size());
2497  return n_lanes;
2498 }
2499 
2500 
2501 
2502 template <int dim, typename Number, typename VectorizedArrayType>
2503 inline unsigned int
2504 MatrixFree<dim, Number, VectorizedArrayType>::get_dofs_per_cell(
2505  const unsigned int dof_handler_index,
2506  const unsigned int active_fe_index) const
2507 {
2508  return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2509 }
2510 
2511 
2512 
2513 template <int dim, typename Number, typename VectorizedArrayType>
2514 inline unsigned int
2515 MatrixFree<dim, Number, VectorizedArrayType>::get_n_q_points(
2516  const unsigned int quad_index,
2517  const unsigned int active_fe_index) const
2518 {
2519  AssertIndexRange(quad_index, mapping_info.cell_data.size());
2520  return mapping_info.cell_data[quad_index]
2521  .descriptor[active_fe_index]
2522  .n_q_points;
2523 }
2524 
2525 
2526 
2527 template <int dim, typename Number, typename VectorizedArrayType>
2528 inline unsigned int
2529 MatrixFree<dim, Number, VectorizedArrayType>::get_dofs_per_face(
2530  const unsigned int dof_handler_index,
2531  const unsigned int active_fe_index) const
2532 {
2533  return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2534 }
2535 
2536 
2537 
2538 template <int dim, typename Number, typename VectorizedArrayType>
2539 inline unsigned int
2540 MatrixFree<dim, Number, VectorizedArrayType>::get_n_q_points_face(
2541  const unsigned int quad_index,
2542  const unsigned int active_fe_index) const
2543 {
2544  AssertIndexRange(quad_index, mapping_info.face_data.size());
2545  return mapping_info.face_data[quad_index]
2546  .descriptor[active_fe_index]
2547  .n_q_points;
2548 }
2549 
2550 
2551 
2552 template <int dim, typename Number, typename VectorizedArrayType>
2553 inline const IndexSet &
2554 MatrixFree<dim, Number, VectorizedArrayType>::get_locally_owned_set(
2555  const unsigned int dof_handler_index) const
2556 {
2557  return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2558 }
2559 
2560 
2561 
2562 template <int dim, typename Number, typename VectorizedArrayType>
2563 inline const IndexSet &
2564 MatrixFree<dim, Number, VectorizedArrayType>::get_ghost_set(
2565  const unsigned int dof_handler_index) const
2566 {
2567  return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2568 }
2569 
2570 
2571 
2572 template <int dim, typename Number, typename VectorizedArrayType>
2573 inline const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
2574 MatrixFree<dim, Number, VectorizedArrayType>::get_shape_info(
2575  const unsigned int dof_handler_index,
2576  const unsigned int index_quad,
2577  const unsigned int index_fe,
2578  const unsigned int active_fe_index,
2579  const unsigned int active_quad_index) const
2580 {
2581  AssertIndexRange(dof_handler_index, dof_info.size());
2582  const unsigned int ind =
2583  dof_info[dof_handler_index].global_base_element_offset + index_fe;
2584  AssertIndexRange(ind, shape_info.size(0));
2585  AssertIndexRange(index_quad, shape_info.size(1));
2586  AssertIndexRange(active_fe_index, shape_info.size(2));
2587  AssertIndexRange(active_quad_index, shape_info.size(3));
2588  return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2589 }
2590 
2591 
2592 
2593 template <int dim, typename Number, typename VectorizedArrayType>
2594 inline const internal::MatrixFreeFunctions::FaceToCellTopology<
2595  VectorizedArrayType::size()> &
2596 MatrixFree<dim, Number, VectorizedArrayType>::get_face_info(
2597  const unsigned int macro_face) const
2598 {
2599  AssertIndexRange(macro_face, face_info.faces.size());
2600  return face_info.faces[macro_face];
2601 }
2602 
2603 
2604 
2605 template <int dim, typename Number, typename VectorizedArrayType>
2606 inline const Table<3, unsigned int> &
2607 MatrixFree<dim, Number, VectorizedArrayType>::get_cell_and_face_to_plain_faces()
2608  const
2609 {
2610  return face_info.cell_and_face_to_plain_faces;
2611 }
2612 
2613 
2614 
2615 template <int dim, typename Number, typename VectorizedArrayType>
2616 inline const Quadrature<dim> &
2617 MatrixFree<dim, Number, VectorizedArrayType>::get_quadrature(
2618  const unsigned int quad_index,
2619  const unsigned int active_fe_index) const
2620 {
2621  AssertIndexRange(quad_index, mapping_info.cell_data.size());
2622  return mapping_info.cell_data[quad_index]
2623  .descriptor[active_fe_index]
2624  .quadrature;
2625 }
2626 
2627 
2628 
2629 template <int dim, typename Number, typename VectorizedArrayType>
2630 inline const Quadrature<dim - 1> &
2631 MatrixFree<dim, Number, VectorizedArrayType>::get_face_quadrature(
2632  const unsigned int quad_index,
2633  const unsigned int active_fe_index) const
2634 {
2635  AssertIndexRange(quad_index, mapping_info.face_data.size());
2636  return mapping_info.face_data[quad_index]
2637  .descriptor[active_fe_index]
2638  .quadrature;
2639 }
2640 
2641 
2642 
2643 template <int dim, typename Number, typename VectorizedArrayType>
2644 inline unsigned int
2645 MatrixFree<dim, Number, VectorizedArrayType>::get_cell_category(
2646  const unsigned int macro_cell) const
2647 {
2648  AssertIndexRange(0, dof_info.size());
2649  AssertIndexRange(macro_cell, dof_info[0].cell_active_fe_index.size());
2650  if (dof_info[0].cell_active_fe_index.empty())
2651  return 0;
2652  else
2653  return dof_info[0].cell_active_fe_index[macro_cell];
2654 }
2655 
2656 
2657 
2658 template <int dim, typename Number, typename VectorizedArrayType>
2659 inline std::pair<unsigned int, unsigned int>
2660 MatrixFree<dim, Number, VectorizedArrayType>::get_face_category(
2661  const unsigned int macro_face) const
2662 {
2663  AssertIndexRange(macro_face, face_info.faces.size());
2664  if (dof_info[0].cell_active_fe_index.empty())
2665  return std::make_pair(0U, 0U);
2666 
2667  std::pair<unsigned int, unsigned int> result;
2668  for (unsigned int v = 0; v < VectorizedArrayType::size() &&
2669  face_info.faces[macro_face].cells_interior[v] !=
2670  numbers::invalid_unsigned_int;
2671  ++v)
2672  result.first = std::max(
2673  result.first,
2674  dof_info[0]
2675  .cell_active_fe_index[face_info.faces[macro_face].cells_interior[v]]);
2676  if (face_info.faces[macro_face].cells_exterior[0] !=
2677  numbers::invalid_unsigned_int)
2678  for (unsigned int v = 0; v < VectorizedArrayType::size() &&
2679  face_info.faces[macro_face].cells_exterior[v] !=
2680  numbers::invalid_unsigned_int;
2681  ++v)
2682  result.second = std::max(
2683  result.first,
2684  dof_info[0]
2685  .cell_active_fe_index[face_info.faces[macro_face].cells_exterior[v]]);
2686  else
2687  result.second = numbers::invalid_unsigned_int;
2688  return result;
2689 }
2690 
2691 
2692 
2693 template <int dim, typename Number, typename VectorizedArrayType>
2694 inline bool
2695 MatrixFree<dim, Number, VectorizedArrayType>::indices_initialized() const
2696 {
2697  return indices_are_initialized;
2698 }
2699 
2700 
2701 
2702 template <int dim, typename Number, typename VectorizedArrayType>
2703 inline bool
2704 MatrixFree<dim, Number, VectorizedArrayType>::mapping_initialized() const
2705 {
2706  return mapping_is_initialized;
2707 }
2708 
2709 
2710 template <int dim, typename Number, typename VectorizedArrayType>
2711 inline unsigned int
2712 MatrixFree<dim, Number, VectorizedArrayType>::get_mg_level() const
2713 {
2714  return mg_level;
2715 }
2716 
2717 
2718 
2719 template <int dim, typename Number, typename VectorizedArrayType>
2720 AlignedVector<VectorizedArrayType> *
2721 MatrixFree<dim, Number, VectorizedArrayType>::acquire_scratch_data() const
2722 {
2723  using list_type =
2724  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
2725  list_type &data = scratch_pad.get();
2726  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
2727  if (it->first == false)
2728  {
2729  it->first = true;
2730  return &it->second;
2731  }
2732  data.emplace_front(true, AlignedVector<VectorizedArrayType>());
2733  return &data.front().second;
2734 }
2735 
2736 
2737 
2738 template <int dim, typename Number, typename VectorizedArrayType>
2739 void
2740 MatrixFree<dim, Number, VectorizedArrayType>::release_scratch_data(
2741  const AlignedVector<VectorizedArrayType> *scratch) const
2742 {
2743  using list_type =
2744  std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
2745  list_type &data = scratch_pad.get();
2746  for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
2747  if (&it->second == scratch)
2748  {
2749  Assert(it->first == true, ExcInternalError());
2750  it->first = false;
2751  return;
2752  }
2753  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2754 }
2755 
2756 
2757 
2758 template <int dim, typename Number, typename VectorizedArrayType>
2759 AlignedVector<Number> *
2760 MatrixFree<dim, Number, VectorizedArrayType>::
2761  acquire_scratch_data_non_threadsafe() const
2762 {
2763  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
2764  scratch_pad_non_threadsafe.begin();
2765  it != scratch_pad_non_threadsafe.end();
2766  ++it)
2767  if (it->first == false)
2768  {
2769  it->first = true;
2770  return &it->second;
2771  }
2772  scratch_pad_non_threadsafe.push_front(
2773  std::make_pair(true, AlignedVector<Number>()));
2774  return &scratch_pad_non_threadsafe.front().second;
2775 }
2776 
2777 
2778 
2779 template <int dim, typename Number, typename VectorizedArrayType>
2780 void
2781 MatrixFree<dim, Number, VectorizedArrayType>::
2782  release_scratch_data_non_threadsafe(
2783  const AlignedVector<Number> *scratch) const
2784 {
2785  for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
2786  scratch_pad_non_threadsafe.begin();
2787  it != scratch_pad_non_threadsafe.end();
2788  ++it)
2789  if (&it->second == scratch)
2790  {
2791  Assert(it->first == true, ExcInternalError());
2792  it->first = false;
2793  return;
2794  }
2795  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2796 }
2797 
2798 
2799 
2800 // ------------------------------ reinit functions ---------------------------
2801 
2802 namespace internal
2803 {
2804  namespace MatrixFreeImplementation
2805  {
2806  template <typename DoFHandlerType>
2807  inline std::vector<IndexSet>
2808  extract_locally_owned_index_sets(
2809  const std::vector<const DoFHandlerType *> &dofh,
2810  const unsigned int level)
2811  {
2812  std::vector<IndexSet> locally_owned_set;
2813  locally_owned_set.reserve(dofh.size());
2814  for (unsigned int j = 0; j < dofh.size(); j++)
2815  if (level == numbers::invalid_unsigned_int)
2816  locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
2817  else
2818  AssertThrow(false, ExcNotImplemented());
2819  return locally_owned_set;
2820  }
2821 
2822  template <int dim, int spacedim>
2823  inline std::vector<IndexSet>
2824  extract_locally_owned_index_sets(
2825  const std::vector<const ::DoFHandler<dim, spacedim> *> &dofh,
2826  const unsigned int level)
2827  {
2828  std::vector<IndexSet> locally_owned_set;
2829  locally_owned_set.reserve(dofh.size());
2830  for (unsigned int j = 0; j < dofh.size(); j++)
2831  if (level == numbers::invalid_unsigned_int)
2832  locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
2833  else
2834  locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
2835  return locally_owned_set;
2836  }
2837  } // namespace MatrixFreeImplementation
2838 } // namespace internal
2839 
2840 
2841 
2842 template <int dim, typename Number, typename VectorizedArrayType>
2843 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2844 void
2845 MatrixFree<dim, Number, VectorizedArrayType>::reinit(
2846  const DoFHandlerType & dof_handler,
2847  const AffineConstraints<number2> &constraints_in,
2848  const QuadratureType & quad,
2849  const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
2850  &additional_data)
2851 {
2852  std::vector<const DoFHandlerType *> dof_handlers;
2853  std::vector<const AffineConstraints<number2> *> constraints;
2854  std::vector<QuadratureType> quads;
2855 
2856  dof_handlers.push_back(&dof_handler);
2857  constraints.push_back(&constraints_in);
2858  quads.push_back(quad);
2859 
2860  std::vector<IndexSet> locally_owned_sets =
2861  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2862  dof_handlers, additional_data.mg_level);
2863 
2864  std::vector<hp::QCollection<1>> quad_hp;
2865  quad_hp.emplace_back(quad);
2866 
2867  internal_reinit(StaticMappingQ1<dim>::mapping,
2868  dof_handlers,
2869  constraints,
2870  locally_owned_sets,
2871  quad_hp,
2872  additional_data);
2873 }
2874 
2875 
2876 
2877 template <int dim, typename Number, typename VectorizedArrayType>
2878 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2879 void
2880 MatrixFree<dim, Number, VectorizedArrayType>::reinit(
2881  const Mapping<dim> & mapping,
2882  const DoFHandlerType & dof_handler,
2883  const AffineConstraints<number2> &constraints_in,
2884  const QuadratureType & quad,
2885  const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
2886  &additional_data)
2887 {
2888  std::vector<const DoFHandlerType *> dof_handlers;
2889  std::vector<const AffineConstraints<number2> *> constraints;
2890 
2891  dof_handlers.push_back(&dof_handler);
2892  constraints.push_back(&constraints_in);
2893 
2894  std::vector<IndexSet> locally_owned_sets =
2895  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2896  dof_handlers, additional_data.mg_level);
2897 
2898  std::vector<hp::QCollection<1>> quad_hp;
2899  quad_hp.emplace_back(quad);
2900 
2901  internal_reinit(mapping,
2902  dof_handlers,
2903  constraints,
2904  locally_owned_sets,
2905  quad_hp,
2906  additional_data);
2907 }
2908 
2909 
2910 
2911 template <int dim, typename Number, typename VectorizedArrayType>
2912 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2913 void
2914 MatrixFree<dim, Number, VectorizedArrayType>::reinit(
2915  const std::vector<const DoFHandlerType *> & dof_handler,
2916  const std::vector<const AffineConstraints<number2> *> &constraint,
2917  const std::vector<QuadratureType> & quad,
2918  const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
2919  &additional_data)
2920 {
2921  std::vector<IndexSet> locally_owned_set =
2922  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2923  dof_handler, additional_data.mg_level);
2924  std::vector<hp::QCollection<1>> quad_hp;
2925  for (unsigned int q = 0; q < quad.size(); ++q)
2926  quad_hp.emplace_back(quad[q]);
2927  internal_reinit(StaticMappingQ1<dim>::mapping,
2928  dof_handler,
2929  constraint,
2930  locally_owned_set,
2931  quad_hp,
2932  additional_data);
2933 }
2934 
2935 
2936 
2937 template <int dim, typename Number, typename VectorizedArrayType>
2938 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2939 void
2940 MatrixFree<dim, Number, VectorizedArrayType>::reinit(
2941  const std::vector<const DoFHandlerType *> & dof_handler,
2942  const std::vector<const AffineConstraints<number2> *> &constraint,
2943  const QuadratureType & quad,
2944  const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
2945  &additional_data)
2946 {
2947  std::vector<IndexSet> locally_owned_set =
2948  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2949  dof_handler, additional_data.mg_level);
2950  std::vector<hp::QCollection<1>> quad_hp;
2951  quad_hp.emplace_back(quad);
2952  internal_reinit(StaticMappingQ1<dim>::mapping,
2953  dof_handler,
2954  constraint,
2955  locally_owned_set,
2956  quad_hp,
2957  additional_data);
2958 }
2959 
2960 
2961 
2962 template <int dim, typename Number, typename VectorizedArrayType>
2963 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2964 void
2965 MatrixFree<dim, Number, VectorizedArrayType>::reinit(
2966  const Mapping<dim> & mapping,
2967  const std::vector<const DoFHandlerType *> & dof_handler,
2968  const std::vector<const AffineConstraints<number2> *> &constraint,
2969  const QuadratureType & quad,
2970  const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
2971  &additional_data)
2972 {
2973  std::vector<IndexSet> locally_owned_set =
2974  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2975  dof_handler, additional_data.mg_level);
2976  std::vector<hp::QCollection<1>> quad_hp;
2977  quad_hp.emplace_back(quad);
2978  internal_reinit(mapping,
2979  dof_handler,
2980  constraint,
2981  locally_owned_set,
2982  quad_hp,
2983  additional_data);
2984 }
2985 
2986 
2987 
2988 template <int dim, typename Number, typename VectorizedArrayType>
2989 template <typename DoFHandlerType, typename QuadratureType, typename number2>
2990 void
2991 MatrixFree<dim, Number, VectorizedArrayType>::reinit(
2992  const Mapping<dim> & mapping,
2993  const std::vector<const DoFHandlerType *> & dof_handler,
2994  const std::vector<const AffineConstraints<number2> *> &constraint,
2995  const std::vector<QuadratureType> & quad,
2996  const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
2997  &additional_data)
2998 {
2999  std::vector<IndexSet> locally_owned_set =
3000  internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3001  dof_handler, additional_data.mg_level);
3002  std::vector<hp::QCollection<1>> quad_hp;
3003  for (unsigned int q = 0; q < quad.size(); ++q)
3004  quad_hp.emplace_back(quad[q]);
3005  internal_reinit(mapping,
3006  dof_handler,
3007  constraint,
3008  locally_owned_set,
3009  quad_hp,
3010  additional_data);
3011 }
3012 
3013 
3014 
3015 template <int dim, typename Number, typename VectorizedArrayType>
3016 template <typename DoFHandlerType, typename QuadratureType, typename number2>
3017 void
3018 MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3019  const Mapping<dim> & mapping,
3020  const std::vector<const DoFHandlerType *> & dof_handler,
3021  const std::vector<const AffineConstraints<number2> *> &constraint,
3022  const std::vector<IndexSet> & locally_owned_set,
3023  const std::vector<QuadratureType> & quad,
3024  const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
3025  &additional_data)
3026 {
3027  // find out whether we use a hp Quadrature or a standard quadrature
3028  std::vector<hp::QCollection<1>> quad_hp;
3029  for (unsigned int q = 0; q < quad.size(); ++q)
3030  quad_hp.emplace_back(quad[q]);
3031  internal_reinit(mapping,
3032  dof_handler,
3033  constraint,
3034  locally_owned_set,
3035  quad_hp,
3036  additional_data);
3037 }
3038 
3039 
3040 
3041 // ------------------------------ implementation of loops --------------------
3042 
3043 // internal helper functions that define how to call MPI data exchange
3044 // functions: for generic vectors, do nothing at all. For distributed vectors,
3045 // call update_ghost_values_start function and so on. If we have collections
3046 // of vectors, just do the individual functions of the components. In order to
3047 // keep ghost values consistent (whether we are in read or write mode), we
3048 // also reset the values at the end. the whole situation is a bit complicated
3049 // by the fact that we need to treat block vectors differently, which use some
3050 // additional helper functions to select the blocks and template magic.
3051 namespace internal
3052 {
3056  template <int dim, typename Number, typename VectorizedArrayType>
3057  struct VectorDataExchange
3058  {
3059  // A shift for the MPI messages to reduce the risk for accidental
3060  // interaction with other open communications that a user program might
3061  // set up (parallel vectors support unfinished communication). We let
3062  // the other vectors use the first 20 assigned numbers and start the
3063  // matrix-free communication.
3064  static constexpr unsigned int channel_shift = 20;
3065 
3066 
3067 
3072  VectorDataExchange(
3073  const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3074  const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3075  DataAccessOnFaces vector_face_access,
3076  const unsigned int n_components)
3077  : matrix_free(matrix_free)
3078  , vector_face_access(
3079  matrix_free.get_task_info().face_partition_data.empty() ?
3080  ::MatrixFree<dim, Number, VectorizedArrayType>::
3081  DataAccessOnFaces::unspecified :
3082  vector_face_access)
3083  , ghosts_were_set(false)
3084 # ifdef DEAL_II_WITH_MPI
3085  , tmp_data(n_components)
3086  , requests(n_components)
3087 # endif
3088  {
3089  (void)n_components;
3090  if (this->vector_face_access !=
3091  ::MatrixFree<dim, Number, VectorizedArrayType>::
3092  DataAccessOnFaces::unspecified)
3093  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3094  AssertDimension(
3095  matrix_free.get_dof_info(c).vector_partitioner_face_variants.size(),
3096  5);
3097  }
3098 
3099 
3100 
3104  ~VectorDataExchange() // NOLINT
3105  {
3106 # ifdef DEAL_II_WITH_MPI
3107  for (unsigned int i = 0; i < tmp_data.size(); ++i)
3108  if (tmp_data[i] != nullptr)
3109  matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3110 # endif
3111  }
3112 
3113 
3114 
3119  template <typename VectorType>
3120  unsigned int
3121  find_vector_in_mf(const VectorType &vec,
3122  const bool check_global_compatibility = true) const
3123  {
3124  unsigned int mf_component = numbers::invalid_unsigned_int;
3125  (void)check_global_compatibility;
3126  for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3127  if (
3128 # ifdef DEBUG
3129  check_global_compatibility ?
3130  vec.get_partitioner()->is_globally_compatible(
3131  *matrix_free.get_dof_info(c).vector_partitioner) :
3132 # endif
3133  vec.get_partitioner()->is_compatible(
3134  *matrix_free.get_dof_info(c).vector_partitioner))
3135  {
3136  mf_component = c;
3137  break;
3138  }
3139  return mf_component;
3140  }
3141 
3142 
3143 
3148  const Utilities::MPI::Partitioner &
3149  get_partitioner(const unsigned int mf_component) const
3150  {
3151  AssertDimension(matrix_free.get_dof_info(mf_component)
3152  .vector_partitioner_face_variants.size(),
3153  5);
3154  if (vector_face_access ==
3155  ::MatrixFree<dim, Number, VectorizedArrayType>::
3156  DataAccessOnFaces::none)
3157  return *matrix_free.get_dof_info(mf_component)
3158  .vector_partitioner_face_variants[0];
3159  else if (vector_face_access ==
3160  ::MatrixFree<dim, Number, VectorizedArrayType>::
3161  DataAccessOnFaces::values)
3162  return *matrix_free.get_dof_info(mf_component)
3163  .vector_partitioner_face_variants[1];
3164  else if (vector_face_access ==
3165  ::MatrixFree<dim, Number, VectorizedArrayType>::
3166  DataAccessOnFaces::gradients)
3167  return *matrix_free.get_dof_info(mf_component)
3168  .vector_partitioner_face_variants[2];
3169  else if (vector_face_access ==
3170  ::MatrixFree<dim, Number, VectorizedArrayType>::
3171  DataAccessOnFaces::values_all_faces)
3172  return *matrix_free.get_dof_info(mf_component)
3173  .vector_partitioner_face_variants[3];
3174  else /*if (vector_face_access ==
3175  ::MatrixFree<dim,
3176  Number>::DataAccessOnFaces::gradients_all_faces)*/
3177  return *matrix_free.get_dof_info(mf_component)
3178  .vector_partitioner_face_variants[4];
3179  }
3180 
3181 
3182 
3186  template <typename VectorType,
3187  typename std::enable_if<is_serial_or_dummy<VectorType>::value,
3188  VectorType>::type * = nullptr>
3189  void
3190  update_ghost_values_start(const unsigned int /*component_in_block_vector*/,
3191  const VectorType & /*vec*/)
3192  {}
3193 
3194 
3199  template <typename VectorType,
3200  typename std::enable_if<
3201  !has_update_ghost_values_start<VectorType>::value &&
3202  !is_serial_or_dummy<VectorType>::value,
3203  VectorType>::type * = nullptr>
3204  void
3205  update_ghost_values_start(const unsigned int component_in_block_vector,
3206  const VectorType & vec)
3207  {
3208  (void)component_in_block_vector;
3209  bool ghosts_set = vec.has_ghost_elements();
3210  if (ghosts_set)
3211  ghosts_were_set = true;
3212 
3213  vec.update_ghost_values();
3214  }
3215 
3216 
3217 
3223  template <typename VectorType,
3224  typename std::enable_if<
3225  has_update_ghost_values_start<VectorType>::value &&
3226  !has_exchange_on_subset<VectorType>::value,
3227  VectorType>::type * = nullptr>
3228  void
3229  update_ghost_values_start(const unsigned int component_in_block_vector,
3230  const VectorType & vec)
3231  {
3232  (void)component_in_block_vector;
3233  bool ghosts_set = vec.has_ghost_elements();
3234  if (ghosts_set)
3235  ghosts_were_set = true;
3236 
3237  vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3238  }
3239 
3240 
3241 
3248  template <typename VectorType,
3249  typename std::enable_if<
3250  has_update_ghost_values_start<VectorType>::value &&
3251  has_exchange_on_subset<VectorType>::value,
3252  VectorType>::type * = nullptr>
3253  void
3254  update_ghost_values_start(const unsigned int component_in_block_vector,
3255  const VectorType & vec)
3256  {
3257  static_assert(
3258  std::is_same<Number, typename VectorType::value_type>::value,
3259  "Type mismatch between VectorType and VectorDataExchange");
3260  (void)component_in_block_vector;
3261  bool ghosts_set = vec.has_ghost_elements();
3262  if (ghosts_set)
3263  ghosts_were_set = true;
3264  if (vector_face_access ==
3265  ::MatrixFree<dim, Number, VectorizedArrayType>::
3266  DataAccessOnFaces::unspecified ||
3267  vec.size() == 0)
3268  vec.update_ghost_values_start(component_in_block_vector +
3269  channel_shift);
3270  else
3271  {
3272 # ifdef DEAL_II_WITH_MPI
3273  const unsigned int mf_component = find_vector_in_mf(vec);
3274  if (&get_partitioner(mf_component) ==
3275  matrix_free.get_dof_info(mf_component).vector_partitioner.get())
3276  {
3277  vec.update_ghost_values_start(component_in_block_vector +
3278  channel_shift);
3279  return;
3280  }
3281 
3282  const Utilities::MPI::Partitioner &part =
3283  get_partitioner(mf_component);
3284  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0)
3285  return;
3286 
3287  tmp_data[component_in_block_vector] =
3288  matrix_free.acquire_scratch_data_non_threadsafe();
3289  tmp_data[component_in_block_vector]->resize_fast(
3290  part.n_import_indices());
3291  AssertDimension(requests.size(), tmp_data.size());
3292 
3293  part.export_to_ghosted_array_start(
3294  component_in_block_vector + channel_shift,
3295  ArrayView<const Number>(vec.begin(), part.local_size()),
3296  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3297  part.n_import_indices()),
3298  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3299  vec.get_partitioner()->local_size(),
3300  vec.get_partitioner()->n_ghost_indices()),
3301  this->requests[component_in_block_vector]);
3302 # endif
3303  }
3304  }
3305 
3306 
3307 
3312  template <
3313  typename VectorType,
3314  typename std::enable_if<!has_update_ghost_values_start<VectorType>::value,
3315  VectorType>::type * = nullptr>
3316  void
3317  update_ghost_values_finish(const unsigned int /*component_in_block_vector*/,
3318  const VectorType & /*vec*/)
3319  {}
3320 
3321 
3322 
3328  template <typename VectorType,
3329  typename std::enable_if<
3330  has_update_ghost_values_start<VectorType>::value &&
3331  !has_exchange_on_subset<VectorType>::value,
3332  VectorType>::type * = nullptr>
3333  void
3334  update_ghost_values_finish(const unsigned int component_in_block_vector,
3335  const VectorType & vec)
3336  {
3337  (void)component_in_block_vector;
3338  vec.update_ghost_values_finish();
3339  }
3340 
3341 
3342 
3349  template <typename VectorType,
3350  typename std::enable_if<
3351  has_update_ghost_values_start<VectorType>::value &&
3352  has_exchange_on_subset<VectorType>::value,
3353  VectorType>::type * = nullptr>
3354  void
3355  update_ghost_values_finish(const unsigned int component_in_block_vector,
3356  const VectorType & vec)
3357  {
3358  static_assert(
3359  std::is_same<Number, typename VectorType::value_type>::value,
3360  "Type mismatch between VectorType and VectorDataExchange");
3361  (void)component_in_block_vector;
3362  if (vector_face_access ==
3363  ::MatrixFree<dim, Number, VectorizedArrayType>::
3364  DataAccessOnFaces::unspecified ||
3365  vec.size() == 0)
3366  vec.update_ghost_values_finish();
3367  else
3368  {
3369 # ifdef DEAL_II_WITH_MPI
3370 
3371  AssertIndexRange(component_in_block_vector, tmp_data.size());
3372  AssertDimension(requests.size(), tmp_data.size());
3373 
3374  const unsigned int mf_component = find_vector_in_mf(vec);
3375  const Utilities::MPI::Partitioner &part =
3376  get_partitioner(mf_component);
3377  if (&part ==
3378  matrix_free.get_dof_info(mf_component).vector_partitioner.get())
3379  {
3380  vec.update_ghost_values_finish();
3381  return;
3382  }
3383 
3384  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0)
3385  return;
3386 
3387  part.export_to_ghosted_array_finish(
3388  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3389  vec.get_partitioner()->local_size(),
3390  vec.get_partitioner()->n_ghost_indices()),
3391  this->requests[component_in_block_vector]);
3392 
3393  matrix_free.release_scratch_data_non_threadsafe(
3394  tmp_data[component_in_block_vector]);
3395  tmp_data[component_in_block_vector] = nullptr;
3396 
3397  // let vector know that ghosts are being updated and we can read from
3398  // them
3399  vec.set_ghost_state(true);
3400 # endif
3401  }
3402  }
3403 
3404 
3405 
3409  template <typename VectorType,
3410  typename std::enable_if<is_serial_or_dummy<VectorType>::value,
3411  VectorType>::type * = nullptr>
3412  void
3413  compress_start(const unsigned int /*component_in_block_vector*/,
3414  VectorType & /*vec*/)
3415  {}
3416 
3417 
3418 
3423  template <typename VectorType,
3424  typename std::enable_if<!has_compress_start<VectorType>::value &&
3425  !is_serial_or_dummy<VectorType>::value,
3426  VectorType>::type * = nullptr>
3427  void
3428  compress_start(const unsigned int component_in_block_vector,
3429  VectorType & vec)
3430  {
3431  (void)component_in_block_vector;
3432  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3433  vec.compress(::VectorOperation::add);
3434  }
3435 
3436 
3437 
3443  template <
3444  typename VectorType,
3445  typename std::enable_if<has_compress_start<VectorType>::value &&
3446  !has_exchange_on_subset<VectorType>::value,
3447  VectorType>::type * = nullptr>
3448  void
3449  compress_start(const unsigned int component_in_block_vector,
3450  VectorType & vec)
3451  {
3452  (void)component_in_block_vector;
3453  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3454  vec.compress_start(component_in_block_vector + channel_shift);
3455  }
3456 
3457 
3458 
3465  template <
3466  typename VectorType,
3467  typename std::enable_if<has_compress_start<VectorType>::value &&
3468  has_exchange_on_subset<VectorType>::value,
3469  VectorType>::type * = nullptr>
3470  void
3471  compress_start(const unsigned int component_in_block_vector,
3472  VectorType & vec)
3473  {
3474  static_assert(
3475  std::is_same<Number, typename VectorType::value_type>::value,
3476  "Type mismatch between VectorType and VectorDataExchange");
3477  (void)component_in_block_vector;
3478  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3479  if (vector_face_access ==
3480  ::MatrixFree<dim, Number, VectorizedArrayType>::
3481  DataAccessOnFaces::unspecified ||
3482  vec.size() == 0)
3483  vec.compress_start(component_in_block_vector + channel_shift);
3484  else
3485  {
3486 # ifdef DEAL_II_WITH_MPI
3487 
3488  const unsigned int mf_component = find_vector_in_mf(vec);
3489  const Utilities::MPI::Partitioner &part =
3490  get_partitioner(mf_component);
3491  if (&part ==
3492  matrix_free.get_dof_info(mf_component).vector_partitioner.get())
3493  {
3494  vec.compress_start(component_in_block_vector + channel_shift);
3495  return;
3496  }
3497 
3498  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0)
3499  return;
3500 
3501  tmp_data[component_in_block_vector] =
3502  matrix_free.acquire_scratch_data_non_threadsafe();
3503  tmp_data[component_in_block_vector]->resize_fast(
3504  part.n_import_indices());
3505  AssertDimension(requests.size(), tmp_data.size());
3506 
3507  part.import_from_ghosted_array_start(
3508  ::VectorOperation::add,
3509  component_in_block_vector + channel_shift,
3510  ArrayView<Number>(vec.begin() + vec.get_partitioner()->local_size(),
3511  vec.get_partitioner()->n_ghost_indices()),
3512  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3513  part.n_import_indices()),
3514  this->requests[component_in_block_vector]);
3515 # endif
3516  }
3517  }
3518 
3519 
3520 
3525  template <typename VectorType,
3526  typename std::enable_if<!has_compress_start<VectorType>::value,
3527  VectorType>::type * = nullptr>
3528  void
3529  compress_finish(const unsigned int /*component_in_block_vector*/,
3530  VectorType & /*vec*/)
3531  {}
3532 
3533 
3534 
3540  template <
3541  typename VectorType,
3542  typename std::enable_if<has_compress_start<VectorType>::value &&
3543  !has_exchange_on_subset<VectorType>::value,
3544  VectorType>::type * = nullptr>
3545  void
3546  compress_finish(const unsigned int component_in_block_vector,
3547  VectorType & vec)
3548  {
3549  (void)component_in_block_vector;
3550  vec.compress_finish(::VectorOperation::add);
3551  }
3552 
3553 
3554 
3561  template <
3562  typename VectorType,
3563  typename std::enable_if<has_compress_start<VectorType>::value &&
3564  has_exchange_on_subset<VectorType>::value,
3565  VectorType>::type * = nullptr>
3566  void
3567  compress_finish(const unsigned int component_in_block_vector,
3568  VectorType & vec)
3569  {
3570  static_assert(
3571  std::is_same<Number, typename VectorType::value_type>::value,
3572  "Type mismatch between VectorType and VectorDataExchange");
3573  (void)component_in_block_vector;
3574  if (vector_face_access ==
3575  ::MatrixFree<dim, Number, VectorizedArrayType>::
3576  DataAccessOnFaces::unspecified ||
3577  vec.size() == 0)
3578  vec.compress_finish(::VectorOperation::add);
3579  else
3580  {
3581 # ifdef DEAL_II_WITH_MPI
3582  AssertIndexRange(component_in_block_vector, tmp_data.size());
3583  AssertDimension(requests.size(), tmp_data.size());
3584 
3585  const unsigned int mf_component = find_vector_in_mf(vec);
3586 
3587  const Utilities::MPI::Partitioner &part =
3588  get_partitioner(mf_component);
3589  if (&part ==
3590  matrix_free.get_dof_info(mf_component).vector_partitioner.get())
3591  {
3592  vec.compress_finish(::VectorOperation::add);
3593  return;
3594  }
3595 
3596  if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0)
3597  return;
3598 
3599  part.import_from_ghosted_array_finish(
3600  VectorOperation::add,
3601  ArrayView<const Number>(
3602  tmp_data[component_in_block_vector]->begin(),
3603  part.n_import_indices()),
3604  ArrayView<Number>(vec.begin(), part.local_size()),
3605  ArrayView<Number>(vec.begin() + vec.get_partitioner()->local_size(),
3606  vec.get_partitioner()->n_ghost_indices()),
3607  this->requests[component_in_block_vector]);
3608 
3609  matrix_free.release_scratch_data_non_threadsafe(
3610  tmp_data[component_in_block_vector]);
3611  tmp_data[component_in_block_vector] = nullptr;
3612 # endif
3613  }
3614  }
3615 
3616 
3617 
3621  template <typename VectorType,
3622  typename std::enable_if<is_serial_or_dummy<VectorType>::value,
3623  VectorType>::type * = nullptr>
3624  void
3625  reset_ghost_values(const VectorType & /*vec*/) const
3626  {}
3627 
3628 
3629 
3634  template <
3635  typename VectorType,
3636  typename std::enable_if<!has_exchange_on_subset<VectorType>::value &&
3637  !is_serial_or_dummy<VectorType>::value,
3638  VectorType>::type * = nullptr>
3639  void
3640  reset_ghost_values(const VectorType &vec) const
3641  {
3642  if (ghosts_were_set == true)
3643  return;
3644 
3645  vec.zero_out_ghosts();
3646  }
3647 
3648 
3649 
3655  template <typename VectorType,
3656  typename std::enable_if<has_exchange_on_subset<VectorType>::value,
3657  VectorType>::type * = nullptr>
3658  void
3659  reset_ghost_values(const VectorType &vec) const
3660  {
3661  static_assert(
3662  std::is_same<Number, typename VectorType::value_type>::value,
3663  "Type mismatch between VectorType and VectorDataExchange");
3664  if (ghosts_were_set == true)
3665  return;
3666 
3667  if (vector_face_access ==
3668  ::MatrixFree<dim, Number, VectorizedArrayType>::
3669  DataAccessOnFaces::unspecified ||
3670  vec.size() == 0)
3671  vec.zero_out_ghosts();
3672  else
3673  {
3674 # ifdef DEAL_II_WITH_MPI
3675  AssertDimension(requests.size(), tmp_data.size());
3676 
3677  const unsigned int mf_component = find_vector_in_mf(vec);
3678  const Utilities::MPI::Partitioner &part =
3679  get_partitioner(mf_component);
3680  if (&part ==
3681  matrix_free.get_dof_info(mf_component).vector_partitioner.get())
3682  vec.zero_out_ghosts();
3683  else if (part.n_ghost_indices() > 0)
3684  {
3685  for (std::vector<std::pair<unsigned int, unsigned int>>::
3686  const_iterator my_ghosts =
3687  part.ghost_indices_within_larger_ghost_set().begin();
3688  my_ghosts !=
3689  part.ghost_indices_within_larger_ghost_set().end();
3690  ++my_ghosts)
3691  for (unsigned int j = my_ghosts->first; j < my_ghosts->second;
3692  j++)
3693  {
3694  const_cast<LinearAlgebra::distributed::Vector<Number> &>(
3695  vec)
3696  .local_element(j + part.local_size()) = 0.;
3697  }
3698  }
3699 
3700  // let vector know that it's not ghosted anymore
3701  vec.set_ghost_state(false);
3702 # endif
3703  }
3704  }
3705 
3706 
3707 
3713  template <typename VectorType,
3714  typename std::enable_if<has_exchange_on_subset<VectorType>::value,
3715  VectorType>::type * = nullptr>
3716  void
3717  zero_vector_region(const unsigned int range_index, VectorType &vec) const
3718  {
3719  static_assert(
3720  std::is_same<Number, typename VectorType::value_type>::value,
3721  "Type mismatch between VectorType and VectorDataExchange");
3722  if (range_index == numbers::invalid_unsigned_int)
3723  vec = Number();
3724  else
3725  {
3726  const unsigned int mf_component = find_vector_in_mf(vec, false);
3727  const internal::MatrixFreeFunctions::DoFInfo &dof_info =
3728  matrix_free.get_dof_info(mf_component);
3729  Assert(dof_info.vector_zero_range_list_index.empty() == false,
3730  ExcNotInitialized());
3731 
3732  Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3733  ExcInternalError());
3734  AssertIndexRange(range_index,
3735  dof_info.vector_zero_range_list_index.size() - 1);
3736  for (unsigned int id =
3737  dof_info.vector_zero_range_list_index[range_index];
3738  id != dof_info.vector_zero_range_list_index[range_index + 1];
3739  ++id)
3740  std::memset(vec.begin() + dof_info.vector_zero_range_list[id].first,
3741  0,
3742  (dof_info.vector_zero_range_list[id].second -
3743  dof_info.vector_zero_range_list[id].first) *
3744  sizeof(Number));
3745  }
3746  }
3747 
3748 
3749 
3755  template <
3756  typename VectorType,
3757  typename std::enable_if<!has_exchange_on_subset<VectorType>::value,
3758  VectorType>::type * = nullptr,
3759  typename VectorType::value_type * = nullptr>
3760  void
3761  zero_vector_region(const unsigned int range_index, VectorType &vec) const
3762  {
3763  if (range_index == numbers::invalid_unsigned_int)
3764  vec = typename VectorType::value_type();
3765  else
3766  {
3767  Assert(false, ExcNotImplemented());
3768  }
3769  }
3770 
3771 
3772 
3777  void
3778  zero_vector_region(const unsigned int, ...) const
3779  {
3780  Assert(false,
3781  ExcNotImplemented("Zeroing is only implemented for vector types "
3782  "which provide VectorType::value_type"));
3783  }
3784 
3785 
3786 
3787  const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
3788  const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3789  DataAccessOnFaces vector_face_access;
3790  bool ghosts_were_set;
3791 # ifdef DEAL_II_WITH_MPI
3792  std::vector<AlignedVector<Number> *> tmp_data;
3793  std::vector<std::vector<MPI_Request>> requests;
3794 # endif
3795  }; // VectorDataExchange
3796 
3797  template <typename VectorStruct>
3798  unsigned int
3799  n_components(const VectorStruct &vec);
3800 
3801  template <typename VectorStruct>
3802  unsigned int
3803  n_components_block(const VectorStruct &vec,
3804  std::integral_constant<bool, true>)
3805  {
3806  unsigned int components = 0;
3807  for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
3808  components += n_components(vec.block(bl));
3809  return components;
3810  }
3811 
3812  template <typename VectorStruct>
3813  unsigned int
3814  n_components_block(const VectorStruct &, std::integral_constant<bool, false>)
3815  {
3816  return 1;
3817  }
3818 
3819  template <typename VectorStruct>
3820  unsigned int
3821  n_components(const VectorStruct &vec)
3822  {
3823  return n_components_block(
3824  vec, std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3825  }
3826 
3827  template <typename VectorStruct>
3828  inline unsigned int
3829  n_components(const std::vector<VectorStruct> &vec)
3830  {
3831  unsigned int components = 0;
3832  for (unsigned int comp = 0; comp < vec.size(); comp++)
3833  components += n_components_block(
3834  vec[comp],
3835  std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3836  return components;
3837  }
3838 
3839  template <typename VectorStruct>
3840  inline unsigned int
3841  n_components(const std::vector<VectorStruct *> &vec)
3842  {
3843  unsigned int components = 0;
3844  for (unsigned int comp = 0; comp < vec.size(); comp++)
3845  components += n_components_block(
3846  *vec[comp],
3847  std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3848  return components;
3849  }
3850 
3851 
3852 
3853  // A helper function to identify block vectors with many components where we
3854  // should not try to overlap computations and communication because there
3855  // would be too many outstanding communication requests.
3856 
3857  // default value for vectors that do not have communication_block_size
3858  template <
3859  typename VectorStruct,
3860  typename std::enable_if<!has_communication_block_size<VectorStruct>::value,
3861  VectorStruct>::type * = nullptr>
3862  constexpr unsigned int
3863  get_communication_block_size(const VectorStruct &)
3864  {
3865  return numbers::invalid_unsigned_int;
3866  }
3867 
3868 
3869 
3870  template <
3871  typename VectorStruct,
3872  typename std::enable_if<has_communication_block_size<VectorStruct>::value,
3873  VectorStruct>::type * = nullptr>
3874  constexpr unsigned int
3875  get_communication_block_size(const VectorStruct &)
3876  {
3877  return VectorStruct::communication_block_size;
3878  }
3879 
3880 
3881 
3882  // --------------------------------------------------------------------------
3883  // below we have wrappers to distinguish between block and non-block vectors.
3884  // --------------------------------------------------------------------------
3885 
3886  //
3887  // update_ghost_values_start
3888  //
3889 
3890  // update_ghost_values for block vectors
3891  template <int dim,
3892  typename VectorStruct,
3893  typename Number,
3894  typename VectorizedArrayType,
3895  typename std::enable_if<IsBlockVector<VectorStruct>::value,
3896  VectorStruct>::type * = nullptr>
3897  void
3898  update_ghost_values_start(
3899  const VectorStruct & vec,
3900  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3901  const unsigned int channel = 0)
3902  {
3903  if (get_communication_block_size(vec) < vec.n_blocks())
3904  {
3905  // don't forget to set ghosts_were_set, that otherwise happens
3906  // inside VectorDataExchange::update_ghost_values_start()
3907  exchanger.ghosts_were_set = vec.has_ghost_elements();
3908  vec.update_ghost_values();
3909  }
3910  else
3911  {
3912  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
3913  update_ghost_values_start(vec.block(i), exchanger, channel + i);
3914  }
3915  }
3916 
3917 
3918 
3919  // update_ghost_values for non-block vectors
3920  template <int dim,
3921  typename VectorStruct,
3922  typename Number,
3923  typename VectorizedArrayType,
3924  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
3925  VectorStruct>::type * = nullptr>
3926  void
3927  update_ghost_values_start(
3928  const VectorStruct & vec,
3929  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3930  const unsigned int channel = 0)
3931  {
3932  exchanger.update_ghost_values_start(channel, vec);
3933  }
3934 
3935 
3936 
3937  // update_ghost_values_start() for vector of vectors
3938  template <int dim,
3939  typename VectorStruct,
3940  typename Number,
3941  typename VectorizedArrayType>
3942  inline void
3943  update_ghost_values_start(
3944  const std::vector<VectorStruct> & vec,
3945  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3946  {
3947  unsigned int component_index = 0;
3948  for (unsigned int comp = 0; comp < vec.size(); comp++)
3949  {
3950  update_ghost_values_start(vec[comp], exchanger, component_index);
3951  component_index += n_components(vec[comp]);
3952  }
3953  }
3954 
3955 
3956 
3957  // update_ghost_values_start() for vector of pointers to vectors
3958  template <int dim,
3959  typename VectorStruct,
3960  typename Number,
3961  typename VectorizedArrayType>
3962  inline void
3963  update_ghost_values_start(
3964  const std::vector<VectorStruct *> & vec,
3965  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
3966  {
3967  unsigned int component_index = 0;
3968  for (unsigned int comp = 0; comp < vec.size(); comp++)
3969  {
3970  update_ghost_values_start(*vec[comp], exchanger, component_index);
3971  component_index += n_components(*vec[comp]);
3972  }
3973  }
3974 
3975 
3976 
3977  //
3978  // update_ghost_values_finish
3979  //
3980 
3981  // for block vectors
3982  template <int dim,
3983  typename VectorStruct,
3984  typename Number,
3985  typename VectorizedArrayType,
3986  typename std::enable_if<IsBlockVector<VectorStruct>::value,
3987  VectorStruct>::type * = nullptr>
3988  void
3989  update_ghost_values_finish(
3990  const VectorStruct & vec,
3991  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
3992  const unsigned int channel = 0)
3993  {
3994  if (get_communication_block_size(vec) < vec.n_blocks())
3995  {
3996  // do nothing, everything has already been completed in the _start()
3997  // call
3998  }
3999  else
4000  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4001  update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4002  }
4003 
4004 
4005 
4006  // for non-block vectors
4007  template <int dim,
4008  typename VectorStruct,
4009  typename Number,
4010  typename VectorizedArrayType,
4011  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4012  VectorStruct>::type * = nullptr>
4013  void
4014  update_ghost_values_finish(
4015  const VectorStruct & vec,
4016  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4017  const unsigned int channel = 0)
4018  {
4019  exchanger.update_ghost_values_finish(channel, vec);
4020  }
4021 
4022 
4023 
4024  // for vector of vectors
4025  template <int dim,
4026  typename VectorStruct,
4027  typename Number,
4028  typename VectorizedArrayType>
4029  inline void
4030  update_ghost_values_finish(
4031  const std::vector<VectorStruct> & vec,
4032  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4033  {
4034  unsigned int component_index = 0;
4035  for (unsigned int comp = 0; comp < vec.size(); comp++)
4036  {
4037  update_ghost_values_finish(vec[comp], exchanger, component_index);
4038  component_index += n_components(vec[comp]);
4039  }
4040  }
4041 
4042 
4043 
4044  // for vector of pointers to vectors
4045  template <int dim,
4046  typename VectorStruct,
4047  typename Number,
4048  typename VectorizedArrayType>
4049  inline void
4050  update_ghost_values_finish(
4051  const std::vector<VectorStruct *> & vec,
4052  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4053  {
4054  unsigned int component_index = 0;
4055  for (unsigned int comp = 0; comp < vec.size(); comp++)
4056  {
4057  update_ghost_values_finish(*vec[comp], exchanger, component_index);
4058  component_index += n_components(*vec[comp]);
4059  }
4060  }
4061 
4062 
4063 
4064  //
4065  // compress_start
4066  //
4067 
4068  // for block vectors
4069  template <int dim,
4070  typename VectorStruct,
4071  typename Number,
4072  typename VectorizedArrayType,
4073  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4074  VectorStruct>::type * = nullptr>
4075  inline void
4076  compress_start(
4077  VectorStruct & vec,
4078  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4079  const unsigned int channel = 0)
4080  {
4081  if (get_communication_block_size(vec) < vec.n_blocks())
4082  vec.compress(::VectorOperation::add);
4083  else
4084  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4085  compress_start(vec.block(i), exchanger, channel + i);
4086  }
4087 
4088 
4089 
4090  // for non-block vectors
4091  template <int dim,
4092  typename VectorStruct,
4093  typename Number,
4094  typename VectorizedArrayType,
4095  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4096  VectorStruct>::type * = nullptr>
4097  inline void
4098  compress_start(
4099  VectorStruct & vec,
4100  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4101  const unsigned int channel = 0)
4102  {
4103  exchanger.compress_start(channel, vec);
4104  }
4105 
4106 
4107 
4108  // for std::vector of vectors
4109  template <int dim,
4110  typename VectorStruct,
4111  typename Number,
4112  typename VectorizedArrayType>
4113  inline void
4114  compress_start(
4115  std::vector<VectorStruct> & vec,
4116  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4117  {
4118  unsigned int component_index = 0;
4119  for (unsigned int comp = 0; comp < vec.size(); comp++)
4120  {
4121  compress_start(vec[comp], exchanger, component_index);
4122  component_index += n_components(vec[comp]);
4123  }
4124  }
4125 
4126 
4127 
4128  // for std::vector of pointer to vectors
4129  template <int dim,
4130  typename VectorStruct,
4131  typename Number,
4132  typename VectorizedArrayType>
4133  inline void
4134  compress_start(
4135  std::vector<VectorStruct *> & vec,
4136  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4137  {
4138  unsigned int component_index = 0;
4139  for (unsigned int comp = 0; comp < vec.size(); comp++)
4140  {
4141  compress_start(*vec[comp], exchanger, component_index);
4142  component_index += n_components(*vec[comp]);
4143  }
4144  }
4145 
4146 
4147 
4148  //
4149  // compress_finish
4150  //
4151 
4152  // for block vectors
4153  template <int dim,
4154  typename VectorStruct,
4155  typename Number,
4156  typename VectorizedArrayType,
4157  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4158  VectorStruct>::type * = nullptr>
4159  inline void
4160  compress_finish(
4161  VectorStruct & vec,
4162  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4163  const unsigned int channel = 0)
4164  {
4165  if (get_communication_block_size(vec) < vec.n_blocks())
4166  {
4167  // do nothing, everything has already been completed in the _start()
4168  // call
4169  }
4170  else
4171  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4172  compress_finish(vec.block(i), exchanger, channel + i);
4173  }
4174 
4175 
4176 
4177  // for non-block vectors
4178  template <int dim,
4179  typename VectorStruct,
4180  typename Number,
4181  typename VectorizedArrayType,
4182  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4183  VectorStruct>::type * = nullptr>
4184  inline void
4185  compress_finish(
4186  VectorStruct & vec,
4187  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4188  const unsigned int channel = 0)
4189  {
4190  exchanger.compress_finish(channel, vec);
4191  }
4192 
4193 
4194 
4195  // for std::vector of vectors
4196  template <int dim,
4197  typename VectorStruct,
4198  typename Number,
4199  typename VectorizedArrayType>
4200  inline void
4201  compress_finish(
4202  std::vector<VectorStruct> & vec,
4203  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4204  {
4205  unsigned int component_index = 0;
4206  for (unsigned int comp = 0; comp < vec.size(); comp++)
4207  {
4208  compress_finish(vec[comp], exchanger, component_index);
4209  component_index += n_components(vec[comp]);
4210  }
4211  }
4212 
4213 
4214 
4215  // for std::vector of pointer to vectors
4216  template <int dim,
4217  typename VectorStruct,
4218  typename Number,
4219  typename VectorizedArrayType>
4220  inline void
4221  compress_finish(
4222  std::vector<VectorStruct *> & vec,
4223  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4224  {
4225  unsigned int component_index = 0;
4226  for (unsigned int comp = 0; comp < vec.size(); comp++)
4227  {
4228  compress_finish(*vec[comp], exchanger, component_index);
4229  component_index += n_components(*vec[comp]);
4230  }
4231  }
4232 
4233 
4234 
4235  //
4236  // reset_ghost_values:
4237  //
4238  // if the input vector did not have ghosts imported, clear them here again
4239  // in order to avoid subsequent operations e.g. in linear solvers to work
4240  // with ghosts all the time
4241  //
4242 
4243  // for block vectors
4244  template <int dim,
4245  typename VectorStruct,
4246  typename Number,
4247  typename VectorizedArrayType,
4248  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4249  VectorStruct>::type * = nullptr>
4250  inline void
4251  reset_ghost_values(
4252  const VectorStruct & vec,
4253  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4254  {
4255  // return immediately if there is nothing to do.
4256  if (exchanger.ghosts_were_set == true)
4257  return;
4258 
4259  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4260  reset_ghost_values(vec.block(i), exchanger);
4261  }
4262 
4263 
4264 
4265  // for non-block vectors
4266  template <int dim,
4267  typename VectorStruct,
4268  typename Number,
4269  typename VectorizedArrayType,
4270  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4271  VectorStruct>::type * = nullptr>
4272  inline void
4273  reset_ghost_values(
4274  const VectorStruct & vec,
4275  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4276  {
4277  exchanger.reset_ghost_values(vec);
4278  }
4279 
4280 
4281 
4282  // for std::vector of vectors
4283  template <int dim,
4284  typename VectorStruct,
4285  typename Number,
4286  typename VectorizedArrayType>
4287  inline void
4288  reset_ghost_values(
4289  const std::vector<VectorStruct> & vec,
4290  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4291  {
4292  // return immediately if there is nothing to do.
4293  if (exchanger.ghosts_were_set == true)
4294  return;
4295 
4296  for (unsigned int comp = 0; comp < vec.size(); comp++)
4297  reset_ghost_values(vec[comp], exchanger);
4298  }
4299 
4300 
4301 
4302  // for std::vector of pointer to vectors
4303  template <int dim,
4304  typename VectorStruct,
4305  typename Number,
4306  typename VectorizedArrayType>
4307  inline void
4308  reset_ghost_values(
4309  const std::vector<VectorStruct *> & vec,
4310  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4311  {
4312  // return immediately if there is nothing to do.
4313  if (exchanger.ghosts_were_set == true)
4314  return;
4315 
4316  for (unsigned int comp = 0; comp < vec.size(); comp++)
4317  reset_ghost_values(*vec[comp], exchanger);
4318  }
4319 
4320 
4321 
4322  //
4323  // zero_vector_region
4324  //
4325 
4326  // for block vectors
4327  template <int dim,
4328  typename VectorStruct,
4329  typename Number,
4330  typename VectorizedArrayType,
4331  typename std::enable_if<IsBlockVector<VectorStruct>::value,
4332  VectorStruct>::type * = nullptr>
4333  inline void
4334  zero_vector_region(
4335  const unsigned int range_index,
4336  VectorStruct & vec,
4337  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4338  {
4339  for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4340  exchanger.zero_vector_region(range_index, vec.block(i));
4341  }
4342 
4343 
4344 
4345  // for non-block vectors
4346  template <int dim,
4347  typename VectorStruct,
4348  typename Number,
4349  typename VectorizedArrayType,
4350  typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4351  VectorStruct>::type * = nullptr>
4352  inline void
4353  zero_vector_region(
4354  const unsigned int range_index,
4355  VectorStruct & vec,
4356  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4357  {
4358  exchanger.zero_vector_region(range_index, vec);
4359  }
4360 
4361 
4362 
4363  // for std::vector of vectors
4364  template <int dim,
4365  typename VectorStruct,
4366  typename Number,
4367  typename VectorizedArrayType>
4368  inline void
4369  zero_vector_region(
4370  const unsigned int range_index,
4371  std::vector<VectorStruct> & vec,
4372  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4373  {
4374  for (unsigned int comp = 0; comp < vec.size(); comp++)
4375  zero_vector_region(range_index, vec[comp], exchanger);
4376  }
4377 
4378 
4379 
4380  // for std::vector of pointers to vectors
4381  template <int dim,
4382  typename VectorStruct,
4383  typename Number,
4384  typename VectorizedArrayType>
4385  inline void
4386  zero_vector_region(
4387  const unsigned int range_index,
4388  std::vector<VectorStruct *> & vec,
4389  VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4390  {
4391  for (unsigned int comp = 0; comp < vec.size(); comp++)
4392  zero_vector_region(range_index, *vec[comp], exchanger);
4393  }
4394 
4395 
4396 
4397  namespace MatrixFreeFunctions
4398  {
4399  // struct to select between a const interface and a non-const interface
4400  // for MFWorker
4401  template <typename, typename, typename, typename, bool>
4402  struct InterfaceSelector
4403  {};
4404 
4405  // Version for constant functions
4406  template <typename MF,
4407  typename InVector,
4408  typename OutVector,
4409  typename Container>
4410  struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4411  {
4412  using function_type = void (Container::*)(
4413  const MF &,
4414  OutVector &,
4415  const InVector &,
4416  const std::pair<unsigned int, unsigned int> &) const;
4417  };
4418 
4419  // Version for non-constant functions
4420  template <typename MF,
4421  typename InVector,
4422  typename OutVector,
4423  typename Container>
4424  struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4425  {
4426  using function_type =
4427  void (Container::*)(const MF &,
4428  OutVector &,
4429  const InVector &,
4430  const std::pair<unsigned int, unsigned int> &);
4431  };
4432  } // namespace MatrixFreeFunctions
4433 
4434 
4435 
4436  // A implementation class for the worker object that runs the various
4437  // operations we want to perform during the matrix-free loop
4438  template <typename MF,
4439  typename InVector,
4440  typename OutVector,
4441  typename Container,
4442  bool is_constant>
4443  class MFWorker : public MFWorkerInterface
4444  {
4445  public:
4446  // An alias to make the arguments further down more readable
4447  using function_type = typename MatrixFreeFunctions::
4448  InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4449  function_type;
4450 
4451  // constructor, binds all the arguments to this class
4452  MFWorker(const MF & matrix_free,
4453  const InVector & src,
4454  OutVector & dst,
4455  const bool zero_dst_vector_setting,
4456  const Container & container,
4457  function_type cell_function,
4458  function_type face_function,
4459  function_type boundary_function,
4460  const typename MF::DataAccessOnFaces src_vector_face_access =
4461  MF::DataAccessOnFaces::none,
4462  const typename MF::DataAccessOnFaces dst_vector_face_access =
4463  MF::DataAccessOnFaces::none,
4464  const std::function<void(const unsigned int, const unsigned int)>
4465  &operation_before_loop = {},
4466  const std::function<void(const unsigned int, const unsigned int)>
4467  & operation_after_loop = {},
4468  const unsigned int dof_handler_index_pre_post = 0)
4469  : matrix_free(matrix_free)
4470  , container(const_cast<Container &>(container))
4471  , cell_function(cell_function)
4472  , face_function(face_function)
4473  , boundary_function(boundary_function)
4474  , src(src)
4475  , dst(dst)
4476  , src_data_exchanger(matrix_free,
4477  src_vector_face_access,
4478  n_components(src))
4479  , dst_data_exchanger(matrix_free,
4480  dst_vector_face_access,
4481  n_components(dst))
4482  , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4483  , zero_dst_vector_setting(zero_dst_vector_setting &&
4484  !src_and_dst_are_same)
4485  , operation_before_loop(operation_before_loop)
4486  , operation_after_loop(operation_after_loop)
4487  , dof_handler_index_pre_post(dof_handler_index_pre_post)
4488  {}
4489 
4490  // Runs the cell work. If no function is given, nothing is done
4491  virtual void
4492  cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4493  {
4494  if (cell_function != nullptr && cell_range.second > cell_range.first)
4495  (container.*
4496  cell_function)(matrix_free, this->dst, this->src, cell_range);
4497  }
4498 
4499  // Runs the assembler on interior faces. If no function is given, nothing
4500  // is done
4501  virtual void
4502  face(const std::pair<unsigned int, unsigned int> &face_range) override
4503  {
4504  if (face_function != nullptr && face_range.second > face_range.first)
4505  (container.*
4506  face_function)(matrix_free, this->dst, this->src, face_range);
4507  }
4508 
4509  // Runs the assembler on boundary faces. If no function is given, nothing
4510  // is done
4511  virtual void
4512  boundary(const std::pair<unsigned int, unsigned int> &face_range) override
4513  {
4514  if (boundary_function != nullptr && face_range.second > face_range.first)
4515  (container.*
4516  boundary_function)(matrix_free, this->dst, this->src, face_range);
4517  }
4518 
4519  // Starts the communication for the update ghost values operation. We
4520  // cannot call this update if ghost and destination are the same because
4521  // that would introduce spurious entries in the destination (there is also
4522  // the problem that reading from a vector that we also write to is usually
4523  // not intended in case there is overlap, but this is up to the
4524  // application code to decide and we cannot catch this case here).
4525  virtual void
4526  vector_update_ghosts_start() override
4527  {
4528  if (!src_and_dst_are_same)
4529  internal::update_ghost_values_start(src, src_data_exchanger);
4530  }
4531 
4532  // Finishes the communication for the update ghost values operation
4533  virtual void
4534  vector_update_ghosts_finish() override
4535  {
4536  if (!src_and_dst_are_same)
4537  internal::update_ghost_values_finish(src, src_data_exchanger);
4538  }
4539 
4540  // Starts the communication for the vector compress operation
4541  virtual void
4542  vector_compress_start() override
4543  {
4544  internal::compress_start(dst, dst_data_exchanger);
4545  }
4546 
4547  // Finishes the communication for the vector compress operation
4548  virtual void
4549  vector_compress_finish() override
4550  {
4551  internal::compress_finish(dst, dst_data_exchanger);
4552  if (!src_and_dst_are_same)
4553  internal::reset_ghost_values(src, src_data_exchanger);
4554  }
4555 
4556  // Zeros the given input vector
4557  virtual void
4558  zero_dst_vector_range(const unsigned int range_index) override
4559  {
4560  if (zero_dst_vector_setting)
4561  internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4562  }
4563 
4564  virtual void
4565  cell_loop_pre_range(const unsigned int range_index) override
4566  {
4567  if (operation_before_loop)
4568  {
4569  const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4570  matrix_free.get_dof_info(dof_handler_index_pre_post);
4571  if (range_index == numbers::invalid_unsigned_int)
4572  {
4573  // Case with threaded loop -> currently no overlap implemented
4574  ::parallel::apply_to_subranges(
4575  0U,
4576  dof_info.vector_partitioner->local_size(),
4577  operation_before_loop,
4578  internal::VectorImplementation::minimum_parallel_grain_size);
4579  }
4580  else
4581  {
4582  AssertIndexRange(range_index,
4583  dof_info.cell_loop_pre_list_index.size() - 1);
4584  for (unsigned int id =
4585  dof_info.cell_loop_pre_list_index[range_index];
4586  id != dof_info.cell_loop_pre_list_index[range_index + 1];
4587  ++id)
4588  operation_before_loop(dof_info.cell_loop_pre_list[id].first,
4589  dof_info.cell_loop_pre_list[id].second);
4590  }
4591  }
4592  }
4593 
4594  virtual void
4595  cell_loop_post_range(const unsigned int range_index) override
4596  {
4597  if (operation_after_loop)
4598  {
4599  const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4600  matrix_free.get_dof_info(dof_handler_index_pre_post);
4601  if (range_index == numbers::invalid_unsigned_int)
4602  {
4603  // Case with threaded loop -> currently no overlap implemented
4604  ::parallel::apply_to_subranges(
4605  0U,
4606  dof_info.vector_partitioner->local_size(),
4607  operation_after_loop,
4608  internal::VectorImplementation::minimum_parallel_grain_size);
4609  }
4610  else
4611  {
4612  AssertIndexRange(range_index,
4613  dof_info.cell_loop_post_list_index.size() - 1);
4614  for (unsigned int id =
4615  dof_info.cell_loop_post_list_index[range_index];
4616  id != dof_info.cell_loop_post_list_index[range_index + 1];
4617  ++id)
4618  operation_after_loop(dof_info.cell_loop_post_list[id].first,
4619  dof_info.cell_loop_post_list[id].second);
4620  }
4621  }
4622  }
4623 
4624  private:
4625  const MF & matrix_free;
4626  Container & container;
4627  function_type cell_function;
4628  function_type face_function;
4629  function_type boundary_function;
4630 
4631  const InVector &src;
4632  OutVector & dst;
4633  VectorDataExchange<MF::dimension,
4634  typename MF::value_type,
4635  typename MF::vectorized_value_type>
4636  src_data_exchanger;
4637  VectorDataExchange<MF::dimension,
4638  typename MF::value_type,
4639  typename MF::vectorized_value_type>
4640  dst_data_exchanger;
4641  const bool src_and_dst_are_same;
4642  const bool zero_dst_vector_setting;
4643  const std::function<void(const unsigned int, const unsigned int)>
4644  operation_before_loop;
4645  const std::function<void(const unsigned int, const unsigned int)>
4646  operation_after_loop;
4647  const unsigned int dof_handler_index_pre_post;
4648  };
4649 
4650 
4651 
4656  template <class MF, typename InVector, typename OutVector>
4657  struct MFClassWrapper
4658  {
4659  using function_type =
4660  std::function<void(const MF &,
4661  OutVector &,
4662  const InVector &,
4663  const std::pair<unsigned int, unsigned int> &)>;
4664 
4665  MFClassWrapper(const function_type cell,
4666  const function_type face,
4667  const function_type boundary)
4668  : cell(cell)
4669  , face(face)
4670  , boundary(boundary)
4671  {}
4672 
4673  void
4674  cell_integrator(const MF & mf,
4675  OutVector & dst,
4676  const InVector & src,
4677  const std::pair<unsigned int, unsigned int> &range) const
4678  {
4679  if (cell)
4680  cell(mf, dst, src, range);
4681  }
4682 
4683  void
4684  face_integrator(const MF & mf,
4685  OutVector & dst,
4686  const InVector & src,
4687  const std::pair<unsigned int, unsigned int> &range) const
4688  {
4689  if (face)
4690  face(mf, dst, src, range);
4691  }
4692 
4693  void
4694  boundary_integrator(
4695  const MF & mf,
4696  OutVector & dst,
4697  const InVector & src,
4698  const std::pair<unsigned int, unsigned int> &range) const
4699  {
4700  if (boundary)
4701  boundary(mf, dst, src, range);
4702  }
4703 
4704  const function_type cell;
4705  const function_type face;
4706  const function_type boundary;
4707  };
4708 
4709 } // end of namespace internal
4710 
4711 
4712 
4713 template <int dim, typename Number, typename VectorizedArrayType>
4714 template <typename OutVector, typename InVector>
4715 inline void
4716 MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
4717  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4718  OutVector &,
4719  const InVector &,
4720  const std::pair<unsigned int, unsigned int> &)>
4721  & cell_operation,
4722  OutVector & dst,
4723  const InVector &src,
4724  const bool zero_dst_vector) const
4725 {
4726  using Wrapper =
4727  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4728  InVector,
4729  OutVector>;
4730  Wrapper wrap(cell_operation, nullptr, nullptr);
4731  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4732  InVector,
4733  OutVector,
4734  Wrapper,
4735  true>
4736  worker(*this,
4737  src,
4738  dst,
4739  zero_dst_vector,
4740  wrap,
4741  &Wrapper::cell_integrator,
4742  &Wrapper::face_integrator,
4743  &Wrapper::boundary_integrator);
4744 
4745  task_info.loop(worker);
4746 }
4747 
4748 
4749 
4750 template <int dim, typename Number, typename VectorizedArrayType>
4751 template <typename OutVector, typename InVector>
4752 inline void
4753 MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
4754  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4755  OutVector &,
4756  const InVector &,
4757  const std::pair<unsigned int, unsigned int> &)>
4758  & cell_operation,
4759  OutVector & dst,
4760  const InVector &src,
4761  const std::function<void(const unsigned int, const unsigned int)>
4762  &operation_before_loop,
4763  const std::function<void(const unsigned int, const unsigned int)>
4764  & operation_after_loop,
4765  const unsigned int dof_handler_index_pre_post) const
4766 {
4767  using Wrapper =
4768  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4769  InVector,
4770  OutVector>;
4771  Wrapper wrap(cell_operation, nullptr, nullptr);
4772  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4773  InVector,
4774  OutVector,
4775  Wrapper,
4776  true>
4777  worker(*this,
4778  src,
4779  dst,
4780  false,
4781  wrap,
4782  &Wrapper::cell_integrator,
4783  &Wrapper::face_integrator,
4784  &Wrapper::boundary_integrator,
4785  DataAccessOnFaces::none,
4786  DataAccessOnFaces::none,
4787  operation_before_loop,
4788  operation_after_loop,
4789  dof_handler_index_pre_post);
4790 
4791  task_info.loop(worker);
4792 }
4793 
4794 
4795 
4796 template <int dim, typename Number, typename VectorizedArrayType>
4797 template <typename OutVector, typename InVector>
4798 inline void
4799 MatrixFree<dim, Number, VectorizedArrayType>::loop(
4800  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4801  OutVector &,
4802  const InVector &,
4803  const std::pair<unsigned int, unsigned int> &)>
4804  &cell_operation,
4805  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4806  OutVector &,
4807  const InVector &,
4808  const std::pair<unsigned int, unsigned int> &)>
4809  &face_operation,
4810  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4811  OutVector &,
4812  const InVector &,
4813  const std::pair<unsigned int, unsigned int> &)>
4814  & boundary_operation,
4815  OutVector & dst,
4816  const InVector & src,
4817  const bool zero_dst_vector,
4818  const DataAccessOnFaces dst_vector_face_access,
4819  const DataAccessOnFaces src_vector_face_access) const
4820 {
4821  using Wrapper =
4822  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4823  InVector,
4824  OutVector>;
4825  Wrapper wrap(cell_operation, face_operation, boundary_operation);
4826  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4827  InVector,
4828  OutVector,
4829  Wrapper,
4830  true>
4831  worker(*this,
4832  src,
4833  dst,
4834  zero_dst_vector,
4835  wrap,
4836  &Wrapper::cell_integrator,
4837  &Wrapper::face_integrator,
4838  &Wrapper::boundary_integrator,
4839  src_vector_face_access,
4840  dst_vector_face_access);
4841 
4842  task_info.loop(worker);
4843 }
4844 
4845 
4846 
4847 template <int dim, typename Number, typename VectorizedArrayType>
4848 template <typename CLASS, typename OutVector, typename InVector>
4849 inline void
4850 MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
4851  void (CLASS::*function_pointer)(
4852  const MatrixFree<dim, Number, VectorizedArrayType> &,
4853  OutVector &,
4854  const InVector &,
4855  const std::pair<unsigned int, unsigned int> &) const,
4856  const CLASS * owning_class,
4857  OutVector & dst,
4858  const InVector &src,
4859  const bool zero_dst_vector) const
4860 {
4861  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4862  InVector,
4863  OutVector,
4864  CLASS,
4865  true>
4866  worker(*this,
4867  src,
4868  dst,
4869  zero_dst_vector,
4870  *owning_class,
4871  function_pointer,
4872  nullptr,
4873  nullptr);
4874  task_info.loop(worker);
4875 }
4876 
4877 
4878 
4879 template <int dim, typename Number, typename VectorizedArrayType>
4880 template <typename CLASS, typename OutVector, typename InVector>
4881 inline void
4882 MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
4883  void (CLASS::*function_pointer)(
4884  const MatrixFree<dim, Number, VectorizedArrayType> &,
4885  OutVector &,
4886  const InVector &,
4887  const std::pair<unsigned int, unsigned int> &) const,
4888  const CLASS * owning_class,
4889  OutVector & dst,
4890  const InVector &src,
4891  const std::function<void(const unsigned int, const unsigned int)>
4892  &operation_before_loop,
4893  const std::function<void(const unsigned int, const unsigned int)>
4894  & operation_after_loop,
4895  const unsigned int dof_handler_index_pre_post) const
4896 {
4897  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4898  InVector,
4899  OutVector,
4900  CLASS,
4901  true>
4902  worker(*this,
4903  src,
4904  dst,
4905  false,
4906  *owning_class,
4907  function_pointer,
4908  nullptr,
4909  nullptr,
4910  DataAccessOnFaces::none,
4911  DataAccessOnFaces::none,
4912  operation_before_loop,
4913  operation_after_loop,
4914  dof_handler_index_pre_post);
4915  task_info.loop(worker);
4916 }
4917 
4918 
4919 
4920 template <int dim, typename Number, typename VectorizedArrayType>
4921 template <typename CLASS, typename OutVector, typename InVector>
4922 inline void
4923 MatrixFree<dim, Number, VectorizedArrayType>::loop(
4924  void (CLASS::*cell_operation)(
4925  const MatrixFree<dim, Number, VectorizedArrayType> &,
4926  OutVector &,
4927  const InVector &,
4928  const std::pair<unsigned int, unsigned int> &) const,
4929  void (CLASS::*face_operation)(
4930  const MatrixFree<dim, Number, VectorizedArrayType> &,
4931  OutVector &,
4932  const InVector &,
4933  const std::pair<unsigned int, unsigned int> &) const,
4934  void (CLASS::*boundary_operation)(
4935  const MatrixFree<dim, Number, VectorizedArrayType> &,
4936  OutVector &,
4937  const InVector &,
4938  const std::pair<unsigned int, unsigned int> &) const,
4939  const CLASS * owning_class,
4940  OutVector & dst,
4941  const InVector & src,
4942  const bool zero_dst_vector,
4943  const DataAccessOnFaces dst_vector_face_access,
4944  const DataAccessOnFaces src_vector_face_access) const
4945 {
4946  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4947  InVector,
4948  OutVector,
4949  CLASS,
4950  true>
4951  worker(*this,
4952  src,
4953  dst,
4954  zero_dst_vector,
4955  *owning_class,
4956  cell_operation,
4957  face_operation,
4958  boundary_operation,
4959  src_vector_face_access,
4960  dst_vector_face_access);
4961  task_info.loop(worker);
4962 }
4963 
4964 
4965 
4966 template <int dim, typename Number, typename VectorizedArrayType>
4967 template <typename CLASS, typename OutVector, typename InVector>
4968 inline void
4969 MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
4970  void (CLASS::*function_pointer)(
4971  const MatrixFree<dim, Number, VectorizedArrayType> &,
4972  OutVector &,
4973  const InVector &,
4974  const std::pair<unsigned int, unsigned int> &),
4975  CLASS * owning_class,
4976  OutVector & dst,
4977  const InVector &src,
4978  const bool zero_dst_vector) const
4979 {
4980  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4981  InVector,
4982  OutVector,
4983  CLASS,
4984  false>
4985  worker(*this,
4986  src,
4987  dst,
4988  zero_dst_vector,
4989  *owning_class,
4990  function_pointer,
4991  nullptr,
4992  nullptr);
4993  task_info.loop(worker);
4994 }
4995 
4996 
4997 
4998 template <int dim, typename Number, typename VectorizedArrayType>
4999 template <typename CLASS, typename OutVector, typename InVector>
5000 inline void
5001 MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5002  void (CLASS::*function_pointer)(
5003  const MatrixFree<dim, Number, VectorizedArrayType> &,
5004  OutVector &,
5005  const InVector &,
5006  const std::pair<unsigned int, unsigned int> &),
5007  CLASS * owning_class,
5008  OutVector & dst,
5009  const InVector &src,
5010  const std::function<void(const unsigned int, const unsigned int)>
5011  &operation_before_loop,
5012  const std::function<void(const unsigned int, const unsigned int)>
5013  & operation_after_loop,
5014  const unsigned int dof_handler_index_pre_post) const
5015 {
5016  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5017  InVector,
5018  OutVector,
5019  CLASS,
5020  false>
5021  worker(*this,
5022  src,
5023  dst,
5024  false,
5025  *owning_class,
5026  function_pointer,
5027  nullptr,
5028  nullptr,
5029  DataAccessOnFaces::none,
5030  DataAccessOnFaces::none,
5031  operation_before_loop,
5032  operation_after_loop,
5033  dof_handler_index_pre_post);
5034  task_info.loop(worker);
5035 }
5036 
5037 
5038 
5039 template <int dim, typename Number, typename VectorizedArrayType>
5040 template <typename CLASS, typename OutVector, typename InVector>
5041 inline void
5042 MatrixFree<dim, Number, VectorizedArrayType>::loop(
5043  void (CLASS::*cell_operation)(
5044  const MatrixFree<dim, Number, VectorizedArrayType> &,
5045  OutVector &,
5046  const InVector &,
5047  const std::pair<unsigned int, unsigned int> &),
5048  void (CLASS::*face_operation)(
5049  const MatrixFree<dim, Number, VectorizedArrayType> &,
5050  OutVector &,
5051  const InVector &,
5052  const std::pair<unsigned int, unsigned int> &),
5053  void (CLASS::*boundary_operation)(
5054  const MatrixFree<dim, Number, VectorizedArrayType> &,
5055  OutVector &,
5056  const InVector &,
5057  const std::pair<unsigned int, unsigned int> &),
5058  CLASS * owning_class,
5059  OutVector & dst,
5060  const InVector & src,
5061  const bool zero_dst_vector,
5062  const DataAccessOnFaces dst_vector_face_access,
5063  const DataAccessOnFaces src_vector_face_access) const
5064 {
5065  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5066  InVector,
5067  OutVector,
5068  CLASS,
5069  false>
5070  worker(*this,
5071  src,
5072  dst,
5073  zero_dst_vector,
5074  *owning_class,
5075  cell_operation,
5076  face_operation,
5077  boundary_operation,
5078  src_vector_face_access,
5079  dst_vector_face_access);
5080  task_info.loop(worker);
5081 }
5082 
5083 
5084 
5085 template <int dim, typename Number, typename VectorizedArrayType>
5086 template <typename CLASS, typename OutVector, typename InVector>
5087 inline void
5088 MatrixFree<dim, Number, VectorizedArrayType>::loop_cell_centric(
5089  void (CLASS::*function_pointer)(
5090  const MatrixFree<dim, Number, VectorizedArrayType> &,
5091  OutVector &,
5092  const InVector &,
5093  const std::pair<unsigned int, unsigned int> &) const,
5094  const CLASS * owning_class,
5095  OutVector & dst,
5096  const InVector & src,
5097  const bool zero_dst_vector,
5098  const DataAccessOnFaces src_vector_face_access) const
5099 {
5100  auto src_vector_face_access_temp = src_vector_face_access;
5101  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5102  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5103  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5104  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5105 
5106  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5107  InVector,
5108  OutVector,
5109  CLASS,
5110  true>
5111  worker(*this,
5112  src,
5113  dst,
5114  zero_dst_vector,
5115  *owning_class,
5116  function_pointer,
5117  nullptr,
5118  nullptr,
5119  src_vector_face_access_temp,
5120  DataAccessOnFaces::none);
5121  task_info.loop(worker);
5122 }
5123 
5124 
5125 
5126 template <int dim, typename Number, typename VectorizedArrayType>
5127 template <typename CLASS, typename OutVector, typename InVector>
5128 inline void
5129 MatrixFree<dim, Number, VectorizedArrayType>::loop_cell_centric(
5130  void (CLASS::*function_pointer)(
5131  const MatrixFree<dim, Number, VectorizedArrayType> &,
5132  OutVector &,
5133  const InVector &,
5134  const std::pair<unsigned int, unsigned int> &),
5135  CLASS * owning_class,
5136  OutVector & dst,
5137  const InVector & src,
5138  const bool zero_dst_vector,
5139  const DataAccessOnFaces src_vector_face_access) const
5140 {
5141  auto src_vector_face_access_temp = src_vector_face_access;
5142  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5143  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5144  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5145  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5146 
5147  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5148  InVector,
5149  OutVector,
5150  CLASS,
5151  false>
5152  worker(*this,
5153  src,
5154  dst,
5155  zero_dst_vector,
5156  *owning_class,
5157  function_pointer,
5158  nullptr,
5159  nullptr,
5160  src_vector_face_access_temp,
5161  DataAccessOnFaces::none);
5162  task_info.loop(worker);
5163 }
5164 
5165 
5166 
5167 template <int dim, typename Number, typename VectorizedArrayType>
5168 template <typename OutVector, typename InVector>
5169 inline void
5170 MatrixFree<dim, Number, VectorizedArrayType>::loop_cell_centric(
5171  const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5172  OutVector &,
5173  const InVector &,
5174  const std::pair<unsigned int, unsigned int> &)>
5175  & cell_operation,
5176  OutVector & dst,
5177  const InVector & src,
5178  const bool zero_dst_vector,
5179  const DataAccessOnFaces src_vector_face_access) const
5180 {
5181  auto src_vector_face_access_temp = src_vector_face_access;
5182  if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5183  src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5184  else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5185  src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5186 
5187  using Wrapper =
5188  internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5189  InVector,
5190  OutVector>;
5191  Wrapper wrap(cell_operation, nullptr, nullptr);
5192 
5193  internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5194  InVector,
5195  OutVector,
5196  Wrapper,
5197  true>
5198  worker(*this,
5199  src,
5200  dst,
5201  zero_dst_vector,
5202  wrap,
5203  &Wrapper::cell_integrator,
5204  &Wrapper::face_integrator,
5205  &Wrapper::boundary_integrator,
5206  src_vector_face_access_temp,
5207  DataAccessOnFaces::none);
5208  task_info.loop(worker);
5209 }
5210 
5211 
5212 #endif // ifndef DOXYGEN
5213 
5214 
5215 
5216 DEAL_II_NAMESPACE_CLOSE
5217 
5218 #endif
MatrixFree::n_ghost_cell_batches
unsigned int n_ghost_cell_batches() const
MatrixFree::initialize_dof_vector
void initialize_dof_vector(VectorType &vec, const unsigned int dof_handler_index=0) const
MatrixFree::indices_are_initialized
bool indices_are_initialized
Definition: matrix_free.h:2133
MatrixFree::initialize_indices
void initialize_indices(const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const AdditionalData &additional_data)
IndexSet
Definition: index_set.h:74
MatrixFree::constraint_pool_end
const Number * constraint_pool_end(const unsigned int pool_index) const
MatrixFree::internal_reinit
void internal_reinit(const Mapping< dim > &mapping, const std::vector< const DoFHandlerType< dim, 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::DoFHandlers
Definition: matrix_free.h:2038
MatrixFree::n_macro_cells
unsigned int n_macro_cells() const
MatrixFree::clear
void clear()
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 dof_handler_index=0) const
MatrixFree::n_cell_batches
unsigned int n_cell_batches() const
MatrixFree::AdditionalData::mg_level
unsigned int mg_level
Definition: matrix_free.h:440
MatrixFree::get_constrained_dofs
const std::vector< unsigned int > & get_constrained_dofs(const unsigned int dof_handler_index=0) const
MatrixFree::AdditionalData::AdditionalData
AdditionalData(const AdditionalData &other)
Definition: matrix_free.h:250
MatrixFree::scratch_pad
Threads::ThreadLocalStorage< std::list< std::pair< bool, AlignedVector< VectorizedArrayType > > > > scratch_pad
Definition: matrix_free.h:2149
Utilities::MPI::Partitioner
Definition: partitioner.h:195
MatrixFree::AdditionalData::cell_vectorization_categories_strict
bool cell_vectorization_categories_strict
Definition: matrix_free.h:527
StaticMappingQ1
Definition: mapping_q1.h:88
MatrixFree< dim, Number, VectorizedArray< Number > >::DataAccessOnFaces
DataAccessOnFaces
Definition: matrix_free.h:719
MatrixFree::AdditionalData::initialize_indices
bool initialize_indices
Definition: matrix_free.h:466
MatrixFree::is_supported
static bool is_supported(const FiniteElement< dim, spacedim > &fe)
MatrixFree::DoFHandlers::n_dof_handlers
unsigned int n_dof_handlers
Definition: matrix_free.h:2058
MatrixFree::release_scratch_data
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
MatrixFree::release_scratch_data_non_threadsafe
void release_scratch_data_non_threadsafe(const AlignedVector< Number > *memory) const
LinearAlgebra::distributed::Vector
Definition: la_parallel_vector.h:226
vectorization.h
block_vector_base.h
MatrixFree::indices_initialized
bool indices_initialized() const
MatrixFree::renumber_dofs
void renumber_dofs(std::vector< types::global_dof_index > &renumbering, const unsigned int dof_handler_index=0)
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType >
hp::DoFHandler::active_cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:312
ArrayView
Definition: array_view.h:77
update_default
@ update_default
No update.
Definition: fe_update_flags.h:69
MatrixFree::DoFHandlers::active_dof_handler
enum MatrixFree::DoFHandlers::ActiveDoFHandler active_dof_handler
MatrixFree::print_memory_consumption
void print_memory_consumption(StreamType &out) const
PointerComparison::equal
static bool equal(const T *p1, const T *p2)
Definition: template_constraints.h:301
internal::MatrixFreeFunctions::FaceToCellTopology
Definition: face_info.h:56
MatrixFree::MatrixFree
MatrixFree()
MatrixFree::get_dofs_per_face
unsigned int get_dofs_per_face(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
mapping_q1.h
internal::MatrixFreeFunctions::TaskInfo::none
@ none
Definition: task_info.h:113
LAPACKSupport::U
static const char U
Definition: lapack_support.h:167
GeometryInfo
Definition: geometry_info.h:1224
VectorType
MatrixFree::AdditionalData::initialize_mapping
bool initialize_mapping
Definition: matrix_free.h:475
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
VectorizedArray< Number >
hp::QCollection
Definition: q_collection.h:48
MatrixFree::AdditionalData::mapping_update_flags
UpdateFlags mapping_update_flags
Definition: matrix_free.h:360
DoFHandler::n_components
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
VectorOperation::add
@ add
Definition: vector_operation.h:53
MatrixFree::shape_info
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > > shape_info
Definition: matrix_free.h:2097
MatrixFree::cell_level_index
std::vector< std::pair< unsigned int, unsigned int > > cell_level_index
Definition: matrix_free.h:2105
shape_info.h
MatrixFree
Definition: matrix_free.h:117
types::boundary_id
unsigned int boundary_id
Definition: types.h:129
MatrixFree::AdditionalData::partition_color
@ partition_color
Definition: matrix_free.h:203
MatrixFree::at_irregular_cell
bool at_irregular_cell(const unsigned int macro_cell_number) const
second
Point< 2 > second
Definition: grid_out.cc:4353
MatrixFree::mapping_info
internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > mapping_info
Definition: matrix_free.h:2091
DoFTools::n_components
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
MatrixFree::get_cell_category
unsigned int get_cell_category(const unsigned int macro_cell) const
internal::MatrixFreeFunctions::DoFInfo
Definition: dof_info.h:99
hp::DoFHandler< dim >
DoFHandler
Definition: dof_handler.h:205
MatrixFree::~MatrixFree
~MatrixFree() override=default
type_traits.h
la_parallel_vector.h
thread_local_storage.h
MatrixFree::get_cell_level_and_index
std::pair< int, int > get_cell_level_and_index(const unsigned int macro_cell_number, const unsigned int vector_number) const
MatrixFree::constraint_pool_row_index
std::vector< unsigned int > constraint_pool_row_index
Definition: matrix_free.h:2084
MatrixFree::n_physical_cells
unsigned int n_physical_cells() const
MatrixFree::acquire_scratch_data
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
internal::MatrixFreeFunctions::TaskInfo
Definition: task_info.h:105
MatrixFree::DoFHandlers::dof_handler
std::vector< SmartPointer< const DoFHandler< dim > > > dof_handler
Definition: matrix_free.h:2045
Table< 3, unsigned int >
MatrixFree::get_dof_info
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
internal::MatrixFreeFunctions::DoFInfo::vector_partitioner
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition: dof_info.h:504
Subscriptor
Definition: subscriptor.h:62
MatrixFree::AdditionalData::none
@ none
Definition: matrix_free.h:194
MatrixFree::get_dof_handler
const DoFHandlerType & get_dof_handler(const unsigned int dof_handler_index=0) const
MatrixFree::dof_info
std::vector< internal::MatrixFreeFunctions::DoFInfo > dof_info
Definition: matrix_free.h:2070
MatrixFree::get_vector_partitioner
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner(const unsigned int dof_handler_index=0) const
MatrixFree::AdditionalData::mapping_update_flags_boundary_faces
UpdateFlags mapping_update_flags_boundary_faces
Definition: matrix_free.h:381
level
unsigned int level
Definition: grid_out.cc:4355
AlignedVector::begin
iterator begin()
internal::MatrixFreeFunctions::DoFInfo::cell_loop_post_list_index
std::vector< unsigned int > cell_loop_post_list_index
Definition: dof_info.h:662
MatrixFree::AdditionalData::TasksParallelScheme
TasksParallelScheme
Definition: matrix_free.h:189
internal::MatrixFreeFunctions::FaceInfo
Definition: face_info.h:127
internal::MatrixFreeFunctions::DoFInfo::vector_zero_range_list_index
std::vector< unsigned int > vector_zero_range_list_index
Definition: dof_info.h:637
DataOutBase::none
@ none
Definition: data_out_base.h:1557
Mapping< dim >
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
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
MatrixFree::get_quadrature
const Quadrature< dim > & get_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
MatrixFree::AdditionalData
Definition: matrix_free.h:182
MatrixFree::copy_from
void copy_from(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free_base)
StandardExceptions::ExcIndexRange
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
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::n_active_entries_per_face_batch
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_number) const
MatrixFree::loop_cell_centric
void loop_cell_centric(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
internal::MatrixFreeFunctions::TaskInfo::partition_partition
@ partition_partition
Definition: task_info.h:114
MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells
UpdateFlags mapping_update_flags_faces_by_cells
Definition: matrix_free.h:430
vector_operation.h
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
double
TrilinosWrappers::internal::begin
VectorType::value_type * begin(VectorType &V)
Definition: trilinos_sparse_matrix.cc:51
MatrixFree::constraint_pool_data
std::vector< Number > constraint_pool_data
Definition: matrix_free.h:2078
MatrixFree::update_mapping
void update_mapping(const Mapping< dim > &mapping)
dof_info.h
grid_tools.h
MatrixFree::cell_level_index_end_local
unsigned int cell_level_index_end_local
Definition: matrix_free.h:2114
MatrixFree::AdditionalData::level_mg_handler
unsigned int & level_mg_handler
Definition: matrix_free.h:447
MatrixFree::n_components
unsigned int n_components() const
Utilities::MPI::Partitioner::n_ghost_indices
unsigned int n_ghost_indices() const
internal::MatrixFreeFunctions::DoFInfo::cell_loop_pre_list_index
std::vector< unsigned int > cell_loop_pre_list_index
Definition: dof_info.h:649
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
update_gradients
@ update_gradients
Shape function gradients.
Definition: fe_update_flags.h:84
MatrixFree::n_active_entries_per_cell_batch
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_number) const
MatrixFree::AdditionalData::mapping_update_flags_inner_faces
UpdateFlags mapping_update_flags_inner_faces
Definition: matrix_free.h:402
MatrixFree::get_face_info
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_number) const
MatrixFree::DataAccessOnFaces::gradients_all_faces
@ gradients_all_faces
fe.h
AlignedVector< VectorizedArrayType >
DynamicSparsityPattern
Definition: dynamic_sparsity_pattern.h:323
MatrixFree::mapping_is_initialized
bool mapping_is_initialized
Definition: matrix_free.h:2138
internal::VectorImplementation::minimum_parallel_grain_size
unsigned int minimum_parallel_grain_size
Definition: parallel.cc:34
MatrixFree::task_info
internal::MatrixFreeFunctions::TaskInfo task_info
Definition: matrix_free.h:2121
DEAL_II_DEPRECATED
#define DEAL_II_DEPRECATED
Definition: config.h:98
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
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
MatrixFree::get_task_info
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
StandardExceptions::ExcNotInitialized
static ::ExceptionBase & ExcNotInitialized()
MatrixFree::n_components_filled
unsigned int n_components_filled(const unsigned int cell_batch_number) const
MatrixFree::get_mapping_info
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
Threads::ThreadLocalStorage
A class that provides a separate storage location on each thread that accesses the object.
Definition: thread_local_storage.h:72
MatrixFree::get_mg_level
unsigned int get_mg_level() const
MatrixFree::acquire_scratch_data_non_threadsafe
AlignedVector< Number > * acquire_scratch_data_non_threadsafe() const
Utilities::MPI::Partitioner::n_import_indices
unsigned int n_import_indices() const
MatrixFree::DoFHandlers::DoFHandlers
DoFHandlers()
Definition: matrix_free.h:2040
task_info.h
FiniteElement
Definition: fe.h:648
Utilities::MPI::Partitioner::local_size
unsigned int local_size() const
numbers::invalid_boundary_id
const types::boundary_id invalid_boundary_id
Definition: types.h:234
exceptions.h
MatrixFree::DataAccessOnFaces::values
@ values
UpdateFlags
UpdateFlags
Definition: fe_update_flags.h:66
mapping_info.h
MatrixFree::memory_consumption
std::size_t memory_consumption() const
MatrixFree::DataAccessOnFaces::values_all_faces
@ values_all_faces
MatrixFree::mg_level
unsigned int mg_level
Definition: matrix_free.h:2161
MatrixFree::get_face_iterator
std::pair< typename DoFHandler< dim >::cell_iterator, unsigned int > get_face_iterator(const unsigned int face_batch_number, const unsigned int vector_number, const bool interior=true, const unsigned int fe_component=0) const
MatrixFree::DoFHandlers::ActiveDoFHandler
ActiveDoFHandler
Definition: matrix_free.h:2047
internal::MatrixFreeFunctions::MappingInfo
Definition: mapping_info.h:316
MatrixFree::AdditionalData::operator=
AdditionalData & operator=(const AdditionalData &other)
Definition: matrix_free.h:276
unsigned int
value
static const bool value
Definition: dof_tools_constraints.cc:433
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
AffineConstraints
Definition: affine_constraints.h:180
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition: fe_update_flags.h:129
MatrixFree::mapping_initialized
bool mapping_initialized() const
MatrixFree::AdditionalData::hold_all_faces_to_owned_cells
bool hold_all_faces_to_owned_cells
Definition: matrix_free.h:498
MatrixFree::AdditionalData::color
@ color
Definition: matrix_free.h:209
MatrixFree::AdditionalData::cell_vectorization_category
std::vector< unsigned int > cell_vectorization_category
Definition: matrix_free.h:516
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 mg_level=numbers::invalid_unsigned_int, const bool store_plain_indices=true, const bool initialize_indices=true, const bool initialize_mapping=true, const bool overlap_communication_computation=true, const bool hold_all_faces_to_owned_cells=false, const bool cell_vectorization_categories_strict=false)
Definition: matrix_free.h:215
internal::MatrixFreeFunctions::DoFInfo::cell_loop_pre_list
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_pre_list
Definition: dof_info.h:655
mapping.h
IsBlockVector
Definition: block_vector_base.h:67
MatrixFree::AdditionalData::store_plain_indices
bool store_plain_indices
Definition: matrix_free.h:455
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
MatrixFree::n_boundary_face_batches
unsigned int n_boundary_face_batches() const
aligned_vector.h
MatrixFree::DoFHandlers::hp_dof_handler
std::vector< SmartPointer< const hp::DoFHandler< dim > > > hp_dof_handler
Definition: matrix_free.h:2046
MatrixFree::scratch_pad_non_threadsafe
std::list< std::pair< bool, AlignedVector< Number > > > scratch_pad_non_threadsafe
Definition: matrix_free.h:2156
MatrixFree::get_face_category
std::pair< unsigned int, unsigned int > get_face_category(const unsigned int macro_face) const
MatrixFree::print
void print(std::ostream &out) const
MatrixFree::dof_handlers
DoFHandlers dof_handlers
Definition: matrix_free.h:2064
MatrixFree::AdditionalData::tasks_parallel_scheme
TasksParallelScheme tasks_parallel_scheme
Definition: matrix_free.h:336
dof_handler.h
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::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)
MatrixFree::DataAccessOnFaces::none
@ none
hp
Definition: hp.h:117
MatrixFree::AdditionalData::tasks_block_size
unsigned int tasks_block_size
Definition: matrix_free.h:347
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::constraint_pool_begin
const Number * constraint_pool_begin(const unsigned int pool_index) const
affine_constraints.h
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:191
MatrixFree::n_base_elements
unsigned int n_base_elements(const unsigned int dof_handler_index) const
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
quadrature.h
template_constraints.h
MatrixFree::AdditionalData::partition_partition
@ partition_partition
Definition: matrix_free.h:198
MatrixFree::dimension
static const unsigned int dimension
Definition: matrix_free.h:134
config.h
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
internal::MatrixFreeFunctions::TaskInfo::partition_color
@ partition_color
Definition: task_info.h:115
MatrixFree::DataAccessOnFaces::gradients
@ gradients
internal
Definition: aligned_vector.h:369
MatrixFree::get_cell_and_face_to_plain_faces
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
q_collection.h
MatrixFree::get_ghost_set
const IndexSet & get_ghost_set(const unsigned int dof_handler_index=0) 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
internal::MatrixFreeFunctions::TaskInfo::color
@ color
Definition: task_info.h:116
MatrixFree::get_locally_owned_set
const IndexSet & get_locally_owned_set(const unsigned int dof_handler_index=0) const
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
MatrixFree::loop
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
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::get_faces_by_cells_boundary_id
std::array< types::boundary_id, VectorizedArrayType::size()> get_faces_by_cells_boundary_id(const unsigned int macro_cell, const unsigned int face_number) const
MatrixFree< dim, Number, VectorizedArray< Number > >::value_type
Number value_type
Definition: matrix_free.h:128
Quadrature
Definition: quadrature.h:85
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
internal::MatrixFreeFunctions::DoFInfo::cell_loop_post_list
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_post_list
Definition: dof_info.h:668
MatrixFree::initialize_dof_handlers
void initialize_dof_handlers(const std::vector< const DoFHandler< dim > * > &dof_handlers, const AdditionalData &additional_data)
MatrixFree::n_inner_face_batches
unsigned int n_inner_face_batches() const
MatrixFree::n_ghost_inner_face_batches
unsigned int n_ghost_inner_face_batches() const
MatrixFree::face_info
internal::MatrixFreeFunctions::FaceInfo< VectorizedArrayType::size()> face_info
Definition: matrix_free.h:2128
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
MatrixFree::n_constraint_pool_entries
unsigned int n_constraint_pool_entries() const
MatrixFree::AdditionalData::overlap_communication_computation
bool overlap_communication_computation
Definition: matrix_free.h:489
MatrixFree::get_boundary_id
types::boundary_id get_boundary_id(const unsigned int macro_face) const
MatrixFree::DoFHandlers::usual
@ usual
Definition: matrix_free.h:2052
internal::MatrixFreeFunctions::DoFInfo::vector_zero_range_list
std::vector< std::pair< unsigned int, unsigned int > > vector_zero_range_list
Definition: dof_info.h:642
Utilities::MPI::max
T max(const T &t, const MPI_Comm &mpi_communicator)
MatrixFree::cell_loop
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
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
MatrixFree::get_shape_info
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
parallel::internal::apply_to_subranges
void apply_to_subranges(const tbb::blocked_range< RangeType > &range, const Function &f)
Definition: parallel.h:365
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::DataAccessOnFaces::unspecified
@ unspecified
dof_handler.h
Generated by   doxygen 1.8.17