Reference documentation for deal.II version 9.0.0
matrix_free.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2011 - 2018 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 at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii_matrix_free_h
18 #define dealii_matrix_free_h
19 
20 #include <deal.II/base/aligned_vector.h>
21 #include <deal.II/base/exceptions.h>
22 #include <deal.II/base/quadrature.h>
23 #include <deal.II/base/vectorization.h>
24 #include <deal.II/base/thread_local_storage.h>
25 #include <deal.II/base/template_constraints.h>
26 #include <deal.II/fe/fe.h>
27 #include <deal.II/fe/mapping.h>
28 #include <deal.II/fe/mapping_q1.h>
29 #include <deal.II/lac/vector_operation.h>
30 #include <deal.II/lac/la_parallel_vector.h>
31 #include <deal.II/lac/block_vector_base.h>
32 #include <deal.II/lac/constraint_matrix.h>
33 #include <deal.II/dofs/dof_handler.h>
34 #include <deal.II/grid/grid_tools.h>
35 #include <deal.II/hp/dof_handler.h>
36 #include <deal.II/hp/q_collection.h>
37 #include <deal.II/matrix_free/task_info.h>
38 #include <deal.II/matrix_free/shape_info.h>
39 #include <deal.II/matrix_free/dof_info.h>
40 #include <deal.II/matrix_free/mapping_info.h>
41 
42 #include <stdlib.h>
43 #include <memory>
44 #include <limits>
45 #include <list>
46 
47 
48 DEAL_II_NAMESPACE_OPEN
49 
50 
51 
104 template <int dim, typename Number=double>
105 class MatrixFree : public Subscriptor
106 {
107 public:
112  typedef Number value_type;
113 
117  static const unsigned int dimension = dim;
118 
165  struct AdditionalData
166  {
172  enum TasksParallelScheme
173  {
177  none = internal::MatrixFreeFunctions::TaskInfo::none,
181  partition_partition = internal::MatrixFreeFunctions::TaskInfo::partition_partition,
185  partition_color = internal::MatrixFreeFunctions::TaskInfo::partition_color,
190  color = internal::MatrixFreeFunctions::TaskInfo::color
191  };
192 
196  AdditionalData (const TasksParallelScheme tasks_parallel_scheme = partition_partition,
197  const unsigned int tasks_block_size = 0,
198  const UpdateFlags mapping_update_flags = update_gradients | update_JxW_values,
199  const UpdateFlags mapping_update_flags_boundary_faces = update_default,
200  const UpdateFlags mapping_update_flags_inner_faces = update_default,
201  const UpdateFlags mapping_update_flags_faces_by_cells = update_default,
202  const unsigned int level_mg_handler = numbers::invalid_unsigned_int,
203  const bool store_plain_indices = true,
204  const bool initialize_indices = true,
205  const bool initialize_mapping = true,
206  const bool overlap_communication_computation = true,
207  const bool hold_all_faces_to_owned_cells = false,
208  const bool cell_vectorization_categories_strict = false)
209  :
210  tasks_parallel_scheme (tasks_parallel_scheme),
211  tasks_block_size (tasks_block_size),
212  mapping_update_flags (mapping_update_flags),
213  mapping_update_flags_boundary_faces (mapping_update_flags_boundary_faces),
214  mapping_update_flags_inner_faces (mapping_update_flags_inner_faces),
215  mapping_update_flags_faces_by_cells (mapping_update_flags_faces_by_cells),
216  level_mg_handler (level_mg_handler),
217  store_plain_indices (store_plain_indices),
218  initialize_indices (initialize_indices),
219  initialize_mapping (initialize_mapping),
220  overlap_communication_computation(overlap_communication_computation),
221  hold_all_faces_to_owned_cells(hold_all_faces_to_owned_cells),
222  cell_vectorization_categories_strict(cell_vectorization_categories_strict)
223  {}
224 
261  TasksParallelScheme tasks_parallel_scheme;
262 
272  unsigned int tasks_block_size;
273 
285  UpdateFlags mapping_update_flags;
286 
306  UpdateFlags mapping_update_flags_boundary_faces;
307 
327  UpdateFlags mapping_update_flags_inner_faces;
328 
355  UpdateFlags mapping_update_flags_faces_by_cells;
356 
365  unsigned int level_mg_handler;
366 
373  bool store_plain_indices;
374 
384  bool initialize_indices;
385 
393  bool initialize_mapping;
394 
407  bool overlap_communication_computation;
408 
416  bool hold_all_faces_to_owned_cells;
417 
434  std::vector<unsigned int> cell_vectorization_category;
435 
445  bool cell_vectorization_categories_strict;
446  };
447 
455  MatrixFree ();
456 
460  MatrixFree (const MatrixFree<dim,Number> &other);
461 
465  ~MatrixFree() = default;
466 
478  template <typename DoFHandlerType, typename QuadratureType>
479  void reinit (const Mapping<dim> &mapping,
480  const DoFHandlerType &dof_handler,
481  const ConstraintMatrix &constraint,
482  const QuadratureType &quad,
483  const AdditionalData additional_data = AdditionalData());
484 
489  template <typename DoFHandlerType, typename QuadratureType>
490  void reinit (const DoFHandlerType &dof_handler,
491  const ConstraintMatrix &constraint,
492  const QuadratureType &quad,
493  const AdditionalData additional_data = AdditionalData());
494 
502  template <typename DoFHandlerType, typename QuadratureType>
503  DEAL_II_DEPRECATED
504  void reinit (const Mapping<dim> &mapping,
505  const DoFHandlerType &dof_handler,
506  const ConstraintMatrix &constraint,
507  const IndexSet &locally_owned_dofs,
508  const QuadratureType &quad,
509  const AdditionalData additional_data = AdditionalData());
510 
531  template <typename DoFHandlerType, typename QuadratureType>
532  void reinit (const Mapping<dim> &mapping,
533  const std::vector<const DoFHandlerType *> &dof_handler,
534  const std::vector<const ConstraintMatrix *> &constraint,
535  const std::vector<QuadratureType> &quad,
536  const AdditionalData additional_data = AdditionalData());
537 
542  template <typename DoFHandlerType, typename QuadratureType>
543  void reinit (const std::vector<const DoFHandlerType *> &dof_handler,
544  const std::vector<const ConstraintMatrix *> &constraint,
545  const std::vector<QuadratureType> &quad,
546  const AdditionalData additional_data = AdditionalData());
547 
555  template <typename DoFHandlerType, typename QuadratureType>
556  DEAL_II_DEPRECATED
557  void reinit (const Mapping<dim> &mapping,
558  const std::vector<const DoFHandlerType *> &dof_handler,
559  const std::vector<const ConstraintMatrix *> &constraint,
560  const std::vector<IndexSet> &locally_owned_set,
561  const std::vector<QuadratureType> &quad,
562  const AdditionalData additional_data = AdditionalData());
563 
571  template <typename DoFHandlerType, typename QuadratureType>
572  void reinit (const Mapping<dim> &mapping,
573  const std::vector<const DoFHandlerType *> &dof_handler,
574  const std::vector<const ConstraintMatrix *> &constraint,
575  const QuadratureType &quad,
576  const AdditionalData additional_data = AdditionalData());
577 
582  template <typename DoFHandlerType, typename QuadratureType>
583  void reinit (const std::vector<const DoFHandlerType *> &dof_handler,
584  const std::vector<const ConstraintMatrix *> &constraint,
585  const QuadratureType &quad,
586  const AdditionalData additional_data = AdditionalData());
587 
593  void copy_from (const MatrixFree<dim,Number> &matrix_free_base);
594 
599  void clear();
600 
602 
614  enum class DataAccessOnFaces
615  {
622  none,
623 
633  values,
634 
645  gradients,
646 
653  unspecified
654  };
655 
697  template <typename OutVector, typename InVector>
698  void cell_loop (const std::function<void (const MatrixFree<dim,Number> &,
699  OutVector &,
700  const InVector &,
701  const std::pair<unsigned int,
702  unsigned int> &)> &cell_operation,
703  OutVector &dst,
704  const InVector &src,
705  const bool zero_dst_vector = false) const;
706 
748  template <typename CLASS, typename OutVector, typename InVector>
749  void cell_loop (void (CLASS::*cell_operation)(const MatrixFree &,
750  OutVector &,
751  const InVector &,
752  const std::pair<unsigned int,
753  unsigned int> &)const,
754  const CLASS *owning_class,
755  OutVector &dst,
756  const InVector &src,
757  const bool zero_dst_vector = false) const;
758 
762  template <typename CLASS, typename OutVector, typename InVector>
763  void cell_loop (void (CLASS::*cell_operation)(const MatrixFree &,
764  OutVector &,
765  const InVector &,
766  const std::pair<unsigned int,
767  unsigned int> &),
768  CLASS *owning_class,
769  OutVector &dst,
770  const InVector &src,
771  const bool zero_dst_vector = false) const;
772 
848  template <typename OutVector, typename InVector>
849  void loop (const std::function<void (const MatrixFree<dim,Number> &,
850  OutVector &,
851  const InVector &,
852  const std::pair<unsigned int,
853  unsigned int> &)> &cell_operation,
854  const std::function<void (const MatrixFree<dim,Number> &,
855  OutVector &,
856  const InVector &,
857  const std::pair<unsigned int,
858  unsigned int> &)> &face_operation,
859  const std::function<void (const MatrixFree<dim,Number> &,
860  OutVector &,
861  const InVector &,
862  const std::pair<unsigned int,
863  unsigned int> &)> &boundary_operation,
864  OutVector &dst,
865  const InVector &src,
866  const bool zero_dst_vector = false,
867  const DataAccessOnFaces dst_vector_face_access = DataAccessOnFaces::unspecified,
868  const DataAccessOnFaces src_vector_face_access = DataAccessOnFaces::unspecified) const;
869 
954  template <typename CLASS, typename OutVector, typename InVector>
955  void loop (void (CLASS::*cell_operation)(const MatrixFree &,
956  OutVector &,
957  const InVector &,
958  const std::pair<unsigned int,
959  unsigned int> &)const,
960  void (CLASS::*face_operation)(const MatrixFree &,
961  OutVector &,
962  const InVector &,
963  const std::pair<unsigned int,
964  unsigned int> &)const,
965  void (CLASS::*boundary_operation)(const MatrixFree &,
966  OutVector &,
967  const InVector &,
968  const std::pair<unsigned int,
969  unsigned int> &)const,
970  const CLASS *owning_class,
971  OutVector &dst,
972  const InVector &src,
973  const bool zero_dst_vector = false,
974  const DataAccessOnFaces dst_vector_face_access = DataAccessOnFaces::unspecified,
975  const DataAccessOnFaces src_vector_face_access = DataAccessOnFaces::unspecified) const;
976 
980  template <typename CLASS, typename OutVector, typename InVector>
981  void loop (void (CLASS::*cell_operation)(const MatrixFree &,
982  OutVector &,
983  const InVector &,
984  const std::pair<unsigned int,
985  unsigned int> &),
986  void (CLASS::*face_operation)(const MatrixFree &,
987  OutVector &,
988  const InVector &,
989  const std::pair<unsigned int,
990  unsigned int> &),
991  void (CLASS::*boundary_operation)(const MatrixFree &,
992  OutVector &,
993  const InVector &,
994  const std::pair<unsigned int,
995  unsigned int> &),
996  CLASS *owning_class,
997  OutVector &dst,
998  const InVector &src,
999  const bool zero_dst_vector = false,
1000  const DataAccessOnFaces dst_vector_face_access = DataAccessOnFaces::unspecified,
1001  const DataAccessOnFaces src_vector_face_access = DataAccessOnFaces::unspecified) const;
1002 
1010  std::pair<unsigned int,unsigned int>
1011  create_cell_subrange_hp (const std::pair<unsigned int,unsigned int> &range,
1012  const unsigned int fe_degree,
1013  const unsigned int dof_handler_index = 0) const;
1014 
1021  std::pair<unsigned int,unsigned int>
1022  create_cell_subrange_hp_by_index (const std::pair<unsigned int,unsigned int> &range,
1023  const unsigned int fe_index,
1024  const unsigned int dof_handler_index = 0) const;
1025 
1027 
1052  template <typename VectorType>
1053  void initialize_dof_vector(VectorType &vec,
1054  const unsigned int dof_handler_index=0) const;
1055 
1076  template <typename Number2>
1077  void initialize_dof_vector(LinearAlgebra::distributed::Vector<Number2> &vec,
1078  const unsigned int dof_handler_index=0) const;
1079 
1090  const std::shared_ptr<const Utilities::MPI::Partitioner> &
1091  get_vector_partitioner (const unsigned int dof_handler_index=0) const;
1092 
1096  const IndexSet &
1097  get_locally_owned_set (const unsigned int dof_handler_index=0) const;
1098 
1102  const IndexSet &
1103  get_ghost_set (const unsigned int dof_handler_index=0) const;
1104 
1114  const std::vector<unsigned int> &
1115  get_constrained_dofs (const unsigned int dof_handler_index=0) const;
1116 
1127  void renumber_dofs (std::vector<types::global_dof_index> &renumbering,
1128  const unsigned int dof_handler_index=0);
1129 
1131 
1139  template <int spacedim>
1140  static
1141  bool is_supported (const FiniteElement<dim, spacedim> &fe);
1142 
1146  unsigned int n_components () const;
1147 
1152  unsigned int n_base_elements (const unsigned int dof_handler_index) const;
1153 
1161  unsigned int n_physical_cells () const;
1162 
1172  unsigned int n_macro_cells () const;
1173 
1183  unsigned int n_cell_batches () const;
1184 
1191  unsigned int n_ghost_cell_batches () const;
1192 
1200  unsigned int n_inner_face_batches () const;
1201 
1210  unsigned int n_boundary_face_batches () const;
1211 
1216  unsigned int n_ghost_inner_face_batches() const;
1217 
1224  types::boundary_id get_boundary_id (const unsigned int macro_face) const;
1225 
1230  std::array<types::boundary_id, VectorizedArray<Number>::n_array_elements>
1231  get_faces_by_cells_boundary_id (const unsigned int macro_cell,
1232  const unsigned int face_number) const;
1233 
1238  const DoFHandler<dim> &
1239  get_dof_handler (const unsigned int dof_handler_index = 0) const;
1240 
1251  typename DoFHandler<dim>::cell_iterator
1252  get_cell_iterator (const unsigned int macro_cell_number,
1253  const unsigned int vector_number,
1254  const unsigned int fe_component = 0) const;
1255 
1267  typename hp::DoFHandler<dim>::active_cell_iterator
1268  get_hp_cell_iterator (const unsigned int macro_cell_number,
1269  const unsigned int vector_number,
1270  const unsigned int dof_handler_index = 0) const;
1271 
1283  bool
1284  at_irregular_cell (const unsigned int macro_cell_number) const;
1285 
1294  unsigned int
1295  n_components_filled (const unsigned int cell_batch_number) const;
1296 
1305  unsigned int
1306  n_active_entries_per_cell_batch (const unsigned int cell_batch_number) const;
1307 
1317  unsigned int
1318  n_active_entries_per_face_batch (const unsigned int face_batch_number) const;
1319 
1323  unsigned int
1324  get_dofs_per_cell (const unsigned int dof_handler_index = 0,
1325  const unsigned int hp_active_fe_index = 0) const;
1326 
1330  unsigned int
1331  get_n_q_points (const unsigned int quad_index = 0,
1332  const unsigned int hp_active_fe_index = 0) const;
1333 
1338  unsigned int
1339  get_dofs_per_face (const unsigned int fe_component = 0,
1340  const unsigned int hp_active_fe_index = 0) const;
1341 
1346  unsigned int
1347  get_n_q_points_face (const unsigned int quad_index = 0,
1348  const unsigned int hp_active_fe_index = 0) const;
1349 
1353  const Quadrature<dim> &
1354  get_quadrature (const unsigned int quad_index = 0,
1355  const unsigned int hp_active_fe_index = 0) const;
1356 
1360  const Quadrature<dim-1> &
1361  get_face_quadrature (const unsigned int quad_index = 0,
1362  const unsigned int hp_active_fe_index = 0) const;
1363 
1370  unsigned int get_cell_category (const unsigned int macro_cell) const;
1371 
1376  std::pair<unsigned int,unsigned int>
1377  get_face_category (const unsigned int macro_face) const;
1378 
1382  bool indices_initialized () const;
1383 
1389  bool mapping_initialized () const;
1390 
1395  std::size_t memory_consumption() const;
1396 
1401  template <typename StreamType>
1402  void print_memory_consumption(StreamType &out) const;
1403 
1408  void print (std::ostream &out) const;
1409 
1411 
1419  const internal::MatrixFreeFunctions::TaskInfo &
1420  get_task_info () const;
1421 
1425  DEAL_II_DEPRECATED
1426  const internal::MatrixFreeFunctions::TaskInfo &
1427  get_size_info () const;
1428 
1429  /*
1430  * Return geometry-dependent information on the cells.
1431  */
1432  const internal::MatrixFreeFunctions::MappingInfo<dim,Number> &
1433  get_mapping_info () const;
1434 
1438  const internal::MatrixFreeFunctions::DoFInfo &
1439  get_dof_info (const unsigned int dof_handler_index_component = 0) const;
1440 
1444  unsigned int n_constraint_pool_entries() const;
1445 
1450  const Number *
1451  constraint_pool_begin (const unsigned int pool_index) const;
1452 
1458  const Number *
1459  constraint_pool_end (const unsigned int pool_index) const;
1460 
1464  const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &
1465  get_shape_info (const unsigned int dof_handler_index_component = 0,
1466  const unsigned int quad_index = 0,
1467  const unsigned int fe_base_element = 0,
1468  const unsigned int hp_active_fe_index = 0,
1469  const unsigned int hp_active_quad_index = 0) const;
1470 
1474  const internal::MatrixFreeFunctions::FaceToCellTopology<VectorizedArray<Number>::n_array_elements> &
1475  get_face_info (const unsigned int face_batch_number) const;
1476 
1490  AlignedVector<VectorizedArray<Number> > *acquire_scratch_data() const;
1491 
1495  void release_scratch_data(const AlignedVector<VectorizedArray<Number> > *memory) const;
1496 
1506  AlignedVector<Number> *acquire_scratch_data_non_threadsafe() const;
1507 
1511  void release_scratch_data_non_threadsafe(const AlignedVector<Number> *memory) const;
1512 
1514 
1515 private:
1516 
1521  void internal_reinit (const Mapping<dim> &mapping,
1522  const std::vector<const DoFHandler<dim> *> &dof_handler,
1523  const std::vector<const ConstraintMatrix *> &constraint,
1524  const std::vector<IndexSet> &locally_owned_set,
1525  const std::vector<hp::QCollection<1> > &quad,
1526  const AdditionalData &additional_data);
1527 
1531  void internal_reinit (const Mapping<dim> &mapping,
1532  const std::vector<const hp::DoFHandler<dim>*> &dof_handler,
1533  const std::vector<const ConstraintMatrix *> &constraint,
1534  const std::vector<IndexSet> &locally_owned_set,
1535  const std::vector<hp::QCollection<1> > &quad,
1536  const AdditionalData &additional_data);
1537 
1544  void
1545  initialize_indices (const std::vector<const ConstraintMatrix *> &constraint,
1546  const std::vector<IndexSet> &locally_owned_set,
1547  const AdditionalData &additional_data);
1548 
1552  void initialize_dof_handlers (const std::vector<const DoFHandler<dim>*> &dof_handlers,
1553  const AdditionalData &additional_data);
1554 
1558  void initialize_dof_handlers (const std::vector<const hp::DoFHandler<dim>*> &dof_handlers,
1559  const AdditionalData &additional_data);
1560 
1565  void make_connectivity_graph_faces (DynamicSparsityPattern &connectivity);
1566 
1572  struct DoFHandlers
1573  {
1574  DoFHandlers ()
1575  :
1576  active_dof_handler(usual),
1577  n_dof_handlers (0),
1578  level (numbers::invalid_unsigned_int)
1579  {}
1580 
1581  std::vector<SmartPointer<const DoFHandler<dim> > > dof_handler;
1582  std::vector<SmartPointer<const hp::DoFHandler<dim> > > hp_dof_handler;
1583  enum ActiveDoFHandler
1584  {
1588  usual,
1592  hp
1593  } active_dof_handler;
1594  unsigned int n_dof_handlers;
1595  unsigned int level;
1596  };
1597 
1601  DoFHandlers dof_handlers;
1602 
1607  std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
1608 
1615  std::vector<Number> constraint_pool_data;
1616 
1621  std::vector<unsigned int> constraint_pool_row_index;
1622 
1627  internal::MatrixFreeFunctions::MappingInfo<dim,Number> mapping_info;
1628 
1632  Table<4,internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>>> shape_info;
1633 
1640  std::vector<std::pair<unsigned int,unsigned int> > cell_level_index;
1641 
1642 
1649  unsigned int cell_level_index_end_local;
1650 
1656  internal::MatrixFreeFunctions::TaskInfo task_info;
1657 
1662  internal::MatrixFreeFunctions::FaceInfo<VectorizedArray<Number>::n_array_elements> face_info;
1663 
1667  bool indices_are_initialized;
1668 
1672  bool mapping_is_initialized;
1673 
1681  mutable Threads::ThreadLocalStorage<std::list<std::pair<bool, AlignedVector<VectorizedArray<Number> > > > > scratch_pad;
1682 
1687  mutable std::list<std::pair<bool, AlignedVector<Number> > > scratch_pad_non_threadsafe;
1688 };
1689 
1690 
1691 
1692 /*----------------------- Inline functions ----------------------------------*/
1693 
1694 #ifndef DOXYGEN
1695 
1696 
1697 
1698 template <int dim, typename Number>
1699 template <typename VectorType>
1700 inline
1701 void
1702 MatrixFree<dim,Number>::initialize_dof_vector(VectorType &vec,
1703  const unsigned int comp) const
1704 {
1705  AssertIndexRange(comp, n_components());
1706  vec.reinit(dof_info[comp].vector_partitioner->size());
1707 }
1708 
1709 
1710 
1711 template <int dim, typename Number>
1712 template <typename Number2>
1713 inline
1714 void
1715 MatrixFree<dim,Number>::initialize_dof_vector(LinearAlgebra::distributed::Vector<Number2> &vec,
1716  const unsigned int comp) const
1717 {
1718  AssertIndexRange(comp, n_components());
1719  vec.reinit(dof_info[comp].vector_partitioner);
1720 }
1721 
1722 
1723 
1724 template <int dim, typename Number>
1725 inline
1726 const std::shared_ptr<const Utilities::MPI::Partitioner> &
1727 MatrixFree<dim,Number>::get_vector_partitioner (const unsigned int comp) const
1728 {
1729  AssertIndexRange(comp, n_components());
1730  return dof_info[comp].vector_partitioner;
1731 }
1732 
1733 
1734 
1735 template <int dim, typename Number>
1736 inline
1737 const std::vector<unsigned int> &
1738 MatrixFree<dim,Number>::get_constrained_dofs (const unsigned int comp) const
1739 {
1740  AssertIndexRange(comp, n_components());
1741  return dof_info[comp].constrained_dofs;
1742 }
1743 
1744 
1745 
1746 template <int dim, typename Number>
1747 inline
1748 unsigned int
1749 MatrixFree<dim,Number>::n_components () const
1750 {
1751  AssertDimension (dof_handlers.n_dof_handlers, dof_info.size());
1752  return dof_handlers.n_dof_handlers;
1753 }
1754 
1755 
1756 
1757 template <int dim, typename Number>
1758 inline
1759 unsigned int
1760 MatrixFree<dim,Number>::n_base_elements (const unsigned int dof_no) const
1761 {
1762  AssertDimension (dof_handlers.n_dof_handlers, dof_info.size());
1763  AssertIndexRange(dof_no, dof_handlers.n_dof_handlers);
1764  return dof_handlers.dof_handler[dof_no]->get_fe().n_base_elements();
1765 }
1766 
1767 
1768 
1769 template <int dim, typename Number>
1770 inline
1771 const internal::MatrixFreeFunctions::TaskInfo &
1772 MatrixFree<dim,Number>::get_task_info () const
1773 {
1774  return task_info;
1775 }
1776 
1777 
1778 
1779 template <int dim, typename Number>
1780 inline
1781 const internal::MatrixFreeFunctions::TaskInfo &
1782 MatrixFree<dim,Number>::get_size_info () const
1783 {
1784  return task_info;
1785 }
1786 
1787 
1788 
1789 template <int dim, typename Number>
1790 inline
1791 unsigned int
1792 MatrixFree<dim,Number>::n_macro_cells () const
1793 {
1794  return *(task_info.cell_partition_data.end()-2);
1795 }
1796 
1797 
1798 
1799 template <int dim, typename Number>
1800 inline
1801 unsigned int
1802 MatrixFree<dim,Number>::n_physical_cells () const
1803 {
1804  return task_info.n_active_cells;
1805 }
1806 
1807 
1808 
1809 template <int dim, typename Number>
1810 inline
1811 unsigned int
1812 MatrixFree<dim,Number>::n_cell_batches () const
1813 {
1814  return *(task_info.cell_partition_data.end()-2);
1815 }
1816 
1817 
1818 
1819 template <int dim, typename Number>
1820 inline
1821 unsigned int
1822 MatrixFree<dim,Number>::n_ghost_cell_batches () const
1823 {
1824  return *(task_info.cell_partition_data.end()-1)-
1825  *(task_info.cell_partition_data.end()-2);
1826 }
1827 
1828 
1829 
1830 template <int dim, typename Number>
1831 inline
1832 unsigned int
1833 MatrixFree<dim,Number>::n_inner_face_batches () const
1834 {
1835  if (task_info.face_partition_data.size() == 0)
1836  return 0;
1837  return task_info.face_partition_data.back();
1838 }
1839 
1840 
1841 
1842 template <int dim, typename Number>
1843 inline
1844 unsigned int
1845 MatrixFree<dim,Number>::n_boundary_face_batches () const
1846 {
1847  if (task_info.face_partition_data.size() == 0)
1848  return 0;
1849  return task_info.boundary_partition_data.back()-task_info.face_partition_data.back();
1850 }
1851 
1852 
1853 
1854 template <int dim, typename Number>
1855 inline
1856 unsigned int
1857 MatrixFree<dim,Number>::n_ghost_inner_face_batches () const
1858 {
1859  if (task_info.face_partition_data.size() == 0)
1860  return 0;
1861  return face_info.faces.size() - task_info.boundary_partition_data.back();
1862 }
1863 
1864 
1865 
1866 template <int dim, typename Number>
1867 inline
1868 types::boundary_id
1869 MatrixFree<dim,Number>::get_boundary_id(const unsigned int macro_face) const
1870 {
1871  Assert(macro_face >= task_info.boundary_partition_data[0] &&
1872  macro_face < task_info.boundary_partition_data.back(),
1873  ExcIndexRange(macro_face,
1874  task_info.boundary_partition_data[0],
1875  task_info.boundary_partition_data.back()));
1876  return types::boundary_id(face_info.faces[macro_face].exterior_face_no);
1877 }
1878 
1879 
1880 
1881 template <int dim, typename Number>
1882 inline
1883 std::array<types::boundary_id, VectorizedArray<Number>::n_array_elements>
1884 MatrixFree<dim,Number>::get_faces_by_cells_boundary_id (const unsigned int macro_cell,
1885  const unsigned int face_number) const
1886 {
1887  AssertIndexRange(macro_cell, n_macro_cells());
1888  AssertIndexRange(face_number, GeometryInfo<dim>::faces_per_cell);
1889  Assert(face_info.cell_and_face_boundary_id.size(0)>=n_macro_cells(),
1890  ExcNotInitialized());
1891  std::array<types::boundary_id, VectorizedArray<Number>::n_array_elements> result;
1892  result.fill(numbers::invalid_boundary_id);
1893  for (unsigned int v=0; v<n_active_entries_per_cell_batch(macro_cell); ++v)
1894  result[v] = face_info.cell_and_face_boundary_id(macro_cell, face_number, v);
1895  return result;
1896 }
1897 
1898 
1899 
1900 template <int dim, typename Number>
1901 inline
1902 const internal::MatrixFreeFunctions::MappingInfo<dim,Number> &
1903 MatrixFree<dim,Number>::get_mapping_info () const
1904 {
1905  return mapping_info;
1906 }
1907 
1908 
1909 
1910 template <int dim, typename Number>
1911 inline
1912 const internal::MatrixFreeFunctions::DoFInfo &
1913 MatrixFree<dim,Number>::get_dof_info (const unsigned int dof_index) const
1914 {
1915  AssertIndexRange (dof_index, n_components());
1916  return dof_info[dof_index];
1917 }
1918 
1919 
1920 
1921 template <int dim, typename Number>
1922 inline
1923 unsigned int
1924 MatrixFree<dim,Number>::n_constraint_pool_entries() const
1925 {
1926  return constraint_pool_row_index.size()-1;
1927 }
1928 
1929 
1930 
1931 template <int dim, typename Number>
1932 inline
1933 const Number *
1934 MatrixFree<dim,Number>::constraint_pool_begin (const unsigned int row) const
1935 {
1936  AssertIndexRange (row, constraint_pool_row_index.size()-1);
1937  return constraint_pool_data.empty() ? nullptr :
1938  constraint_pool_data.data() + constraint_pool_row_index[row];
1939 }
1940 
1941 
1942 
1943 template <int dim, typename Number>
1944 inline
1945 const Number *
1946 MatrixFree<dim,Number>::constraint_pool_end (const unsigned int row) const
1947 {
1948  AssertIndexRange (row, constraint_pool_row_index.size()-1);
1949  return constraint_pool_data.empty() ? nullptr :
1950  constraint_pool_data.data() + constraint_pool_row_index[row+1];
1951 }
1952 
1953 
1954 
1955 template <int dim, typename Number>
1956 inline
1957 std::pair<unsigned int,unsigned int>
1958 MatrixFree<dim,Number>::create_cell_subrange_hp
1959 (const std::pair<unsigned int,unsigned int> &range,
1960  const unsigned int degree,
1961  const unsigned int dof_handler_component) const
1962 {
1963  if (dof_info[dof_handler_component].cell_active_fe_index.empty())
1964  {
1965  AssertDimension (dof_info[dof_handler_component].fe_index_conversion.size(),1);
1966  AssertDimension (dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
1967  if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
1968  return range;
1969  else
1970  return std::pair<unsigned int,unsigned int> (range.second,range.second);
1971  }
1972 
1973  const unsigned int fe_index =
1974  dof_info[dof_handler_component].fe_index_from_degree(0, degree);
1975  if (fe_index >= dof_info[dof_handler_component].max_fe_index)
1976  return std::pair<unsigned int,unsigned int>(range.second, range.second);
1977  else
1978  return create_cell_subrange_hp_by_index (range, fe_index, dof_handler_component);
1979 }
1980 
1981 
1982 
1983 template <int dim, typename Number>
1984 inline
1985 bool
1986 MatrixFree<dim,Number>::at_irregular_cell (const unsigned int macro_cell) const
1987 {
1988  AssertIndexRange (macro_cell, task_info.cell_partition_data.back());
1989  return VectorizedArray<Number>::n_array_elements > 1 &&
1990  cell_level_index[(macro_cell+1)*VectorizedArray<Number>::n_array_elements-1] ==
1991  cell_level_index[(macro_cell+1)*VectorizedArray<Number>::n_array_elements-2];
1992 }
1993 
1994 
1995 
1996 template <int dim, typename Number>
1997 inline
1998 unsigned int
1999 MatrixFree<dim,Number>::n_components_filled (const unsigned int cell_batch_number) const
2000 {
2001  return n_active_entries_per_cell_batch(cell_batch_number);
2002 }
2003 
2004 
2005 
2006 template <int dim, typename Number>
2007 inline
2008 unsigned int
2009 MatrixFree<dim,Number>::n_active_entries_per_cell_batch(const unsigned int cell_batch_number) const
2010 {
2011  AssertIndexRange (cell_batch_number, task_info.cell_partition_data.back());
2012  unsigned int n_components = VectorizedArray<Number>::n_array_elements;
2013  while (n_components > 1 &&
2014  cell_level_index[cell_batch_number*VectorizedArray<Number>::n_array_elements+n_components-1] ==
2015  cell_level_index[cell_batch_number*VectorizedArray<Number>::n_array_elements+n_components-2])
2016  --n_components;
2017  AssertIndexRange(n_components-1, VectorizedArray<Number>::n_array_elements);
2018  return n_components;
2019 }
2020 
2021 
2022 
2023 template <int dim, typename Number>
2024 inline
2025 unsigned int
2026 MatrixFree<dim,Number>::n_active_entries_per_face_batch(const unsigned int face_batch_number) const
2027 {
2028  AssertIndexRange (face_batch_number, face_info.faces.size());
2029  unsigned int n_components = VectorizedArray<Number>::n_array_elements;
2030  while (n_components > 1 &&
2031  face_info.faces[face_batch_number].cells_interior[n_components-1] ==
2032  numbers::invalid_unsigned_int)
2033  --n_components;
2034  AssertIndexRange(n_components-1, VectorizedArray<Number>::n_array_elements);
2035  return n_components;
2036 }
2037 
2038 
2039 
2040 template <int dim, typename Number>
2041 inline
2042 unsigned int
2043 MatrixFree<dim,Number>::get_dofs_per_cell(const unsigned int dof_handler_index,
2044  const unsigned int active_fe_index) const
2045 {
2046  return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2047 }
2048 
2049 
2050 
2051 template <int dim, typename Number>
2052 inline
2053 unsigned int
2054 MatrixFree<dim,Number>::get_n_q_points(const unsigned int quad_index,
2055  const unsigned int active_fe_index) const
2056 {
2057  AssertIndexRange (quad_index, mapping_info.cell_data.size());
2058  return mapping_info.cell_data[quad_index].descriptor[active_fe_index].n_q_points;
2059 }
2060 
2061 
2062 
2063 template <int dim, typename Number>
2064 inline
2065 unsigned int
2066 MatrixFree<dim,Number>::get_dofs_per_face(const unsigned int dof_handler_index,
2067  const unsigned int active_fe_index) const
2068 {
2069  return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2070 }
2071 
2072 
2073 
2074 template <int dim, typename Number>
2075 inline
2076 unsigned int
2077 MatrixFree<dim,Number>::get_n_q_points_face(const unsigned int quad_index,
2078  const unsigned int active_fe_index) const
2079 {
2080  AssertIndexRange (quad_index, mapping_info.face_data.size());
2081  return mapping_info.face_data[quad_index].descriptor[active_fe_index].n_q_points;
2082 }
2083 
2084 
2085 
2086 template <int dim, typename Number>
2087 inline
2088 const IndexSet &
2089 MatrixFree<dim,Number>::get_locally_owned_set(const unsigned int dof_handler_index) const
2090 {
2091  return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2092 }
2093 
2094 
2095 
2096 template <int dim, typename Number>
2097 inline
2098 const IndexSet &
2099 MatrixFree<dim,Number>::get_ghost_set(const unsigned int dof_handler_index) const
2100 {
2101  return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2102 }
2103 
2104 
2105 
2106 template <int dim, typename Number>
2107 inline
2108 const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &
2109 MatrixFree<dim,Number>::get_shape_info (const unsigned int dof_handler_index,
2110  const unsigned int index_quad,
2111  const unsigned int index_fe,
2112  const unsigned int active_fe_index,
2113  const unsigned int active_quad_index) const
2114 {
2115  AssertIndexRange(dof_handler_index, dof_info.size());
2116  const unsigned int ind = dof_info[dof_handler_index].global_base_element_offset+index_fe;
2117  AssertIndexRange (ind, shape_info.size(0));
2118  AssertIndexRange (index_quad, shape_info.size(1));
2119  AssertIndexRange (active_fe_index, shape_info.size(2));
2120  AssertIndexRange (active_quad_index, shape_info.size(3));
2121  return shape_info(ind, index_quad,
2122  active_fe_index, active_quad_index);
2123 }
2124 
2125 
2126 
2127 template <int dim, typename Number>
2128 inline
2129 const internal::MatrixFreeFunctions::FaceToCellTopology<VectorizedArray<Number>::n_array_elements> &
2130 MatrixFree<dim,Number>::get_face_info (const unsigned int macro_face) const
2131 {
2132  AssertIndexRange(macro_face, face_info.faces.size());
2133  return face_info.faces[macro_face];
2134 }
2135 
2136 
2137 
2138 template <int dim, typename Number>
2139 inline
2140 const Quadrature<dim> &
2141 MatrixFree<dim,Number>::get_quadrature (const unsigned int quad_index,
2142  const unsigned int active_fe_index) const
2143 {
2144  AssertIndexRange (quad_index, mapping_info.cell_data.size());
2145  return mapping_info.cell_data[quad_index].descriptor[active_fe_index].quadrature;
2146 }
2147 
2148 
2149 
2150 template <int dim, typename Number>
2151 inline
2152 const Quadrature<dim-1> &
2153 MatrixFree<dim,Number>::get_face_quadrature (const unsigned int quad_index,
2154  const unsigned int active_fe_index) const
2155 {
2156  AssertIndexRange (quad_index, mapping_info.face_data.size());
2157  return mapping_info.face_data[quad_index].descriptor[active_fe_index].quadrature;
2158 }
2159 
2160 
2161 
2162 template <int dim, typename Number>
2163 inline
2164 unsigned int
2165 MatrixFree<dim,Number>::get_cell_category (const unsigned int macro_cell) const
2166 {
2167  AssertIndexRange(0, dof_info.size());
2168  AssertIndexRange(macro_cell, dof_info[0].cell_active_fe_index.size());
2169  if (dof_info[0].cell_active_fe_index.empty())
2170  return 0;
2171  else
2172  return dof_info[0].cell_active_fe_index[macro_cell];
2173 }
2174 
2175 
2176 
2177 template <int dim, typename Number>
2178 inline
2179 std::pair<unsigned int,unsigned int>
2180 MatrixFree<dim,Number>::get_face_category (const unsigned int macro_face) const
2181 {
2182  AssertIndexRange(macro_face, face_info.faces.size());
2183  if (dof_info[0].cell_active_fe_index.empty())
2184  return std::make_pair(0U, 0U);
2185 
2186  std::pair<unsigned int,unsigned int> result;
2187  for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements &&
2188  face_info.faces[macro_face].cells_interior[v] != numbers::invalid_unsigned_int; ++v)
2189  result.first = std::max(result.first,
2190  dof_info[0].cell_active_fe_index[face_info.faces[macro_face].cells_interior[v]]);
2191  if (face_info.faces[macro_face].cells_exterior[0] != numbers::invalid_unsigned_int)
2192  for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements &&
2193  face_info.faces[macro_face].cells_exterior[v] != numbers::invalid_unsigned_int; ++v)
2194  result.second = std::max(result.first,
2195  dof_info[0].cell_active_fe_index[face_info.faces[macro_face].cells_exterior[v]]);
2196  else
2197  result.second = numbers::invalid_unsigned_int;
2198  return result;
2199 }
2200 
2201 
2202 
2203 template <int dim, typename Number>
2204 inline
2205 bool
2206 MatrixFree<dim,Number>::indices_initialized () const
2207 {
2208  return indices_are_initialized;
2209 }
2210 
2211 
2212 
2213 template <int dim, typename Number>
2214 inline
2215 bool
2216 MatrixFree<dim,Number>::mapping_initialized () const
2217 {
2218  return mapping_is_initialized;
2219 }
2220 
2221 
2222 
2223 template <int dim,typename Number>
2224 AlignedVector<VectorizedArray<Number> > *
2225 MatrixFree<dim,Number>::acquire_scratch_data() const
2226 {
2227  typedef std::list<std::pair<bool, AlignedVector<VectorizedArray<Number> > > > list_type;
2228  list_type &data = scratch_pad.get();
2229  for (typename list_type::iterator it=data.begin(); it!=data.end(); ++it)
2230  if (it->first == false)
2231  {
2232  it->first = true;
2233  return &it->second;
2234  }
2235  data.push_front(std::make_pair(true,AlignedVector<VectorizedArray<Number> >()));
2236  return &data.front().second;
2237 }
2238 
2239 
2240 
2241 template <int dim, typename Number>
2242 void
2243 MatrixFree<dim,Number>::release_scratch_data(const AlignedVector<VectorizedArray<Number> > *scratch) const
2244 {
2245  typedef std::list<std::pair<bool, AlignedVector<VectorizedArray<Number> > > > list_type;
2246  list_type &data = scratch_pad.get();
2247  for (typename list_type::iterator it=data.begin(); it!=data.end(); ++it)
2248  if (&it->second == scratch)
2249  {
2250  Assert(it->first == true, ExcInternalError());
2251  it->first = false;
2252  return;
2253  }
2254  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2255 }
2256 
2257 
2258 
2259 template <int dim,typename Number>
2260 AlignedVector<Number> *
2261 MatrixFree<dim,Number>::acquire_scratch_data_non_threadsafe() const
2262 {
2263  for (typename std::list<std::pair<bool, AlignedVector<Number> > >::iterator
2264  it=scratch_pad_non_threadsafe.begin(); it!=scratch_pad_non_threadsafe.end(); ++it)
2265  if (it->first == false)
2266  {
2267  it->first = true;
2268  return &it->second;
2269  }
2270  scratch_pad_non_threadsafe.push_front(std::make_pair(true,AlignedVector<Number>()));
2271  return &scratch_pad_non_threadsafe.front().second;
2272 }
2273 
2274 
2275 
2276 template <int dim, typename Number>
2277 void
2278 MatrixFree<dim,Number>::release_scratch_data_non_threadsafe(const AlignedVector<Number> *scratch) const
2279 {
2280  for (typename std::list<std::pair<bool, AlignedVector<Number> > >::iterator
2281  it=scratch_pad_non_threadsafe.begin(); it!=scratch_pad_non_threadsafe.end(); ++it)
2282  if (&it->second == scratch)
2283  {
2284  Assert(it->first == true, ExcInternalError());
2285  it->first = false;
2286  return;
2287  }
2288  AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2289 }
2290 
2291 
2292 
2293 // ------------------------------ reinit functions ---------------------------
2294 
2295 namespace internal
2296 {
2297  namespace MatrixFreeImplementation
2298  {
2299  template <typename DoFHandlerType>
2300  inline
2301  std::vector<IndexSet>
2302  extract_locally_owned_index_sets (const std::vector<const DoFHandlerType *> &dofh,
2303  const unsigned int level)
2304  {
2305  std::vector<IndexSet> locally_owned_set;
2306  locally_owned_set.reserve (dofh.size());
2307  for (unsigned int j=0; j<dofh.size(); j++)
2308  if (level == numbers::invalid_unsigned_int)
2309  locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
2310  else
2311  AssertThrow(false, ExcNotImplemented());
2312  return locally_owned_set;
2313  }
2314 
2315  template <int dim, int spacedim>
2316  inline
2317  std::vector<IndexSet>
2318  extract_locally_owned_index_sets (const std::vector<const ::DoFHandler<dim,spacedim> *> &dofh,
2319  const unsigned int level)
2320  {
2321  std::vector<IndexSet> locally_owned_set;
2322  locally_owned_set.reserve (dofh.size());
2323  for (unsigned int j=0; j<dofh.size(); j++)
2324  if (level == numbers::invalid_unsigned_int)
2325  locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
2326  else
2327  locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
2328  return locally_owned_set;
2329  }
2330  }
2331 }
2332 
2333 
2334 
2335 template <int dim, typename Number>
2336 template <typename DoFHandlerType, typename QuadratureType>
2337 void MatrixFree<dim,Number>::
2338 reinit(const DoFHandlerType &dof_handler,
2339  const ConstraintMatrix &constraints_in,
2340  const QuadratureType &quad,
2341  const typename MatrixFree<dim,Number>::AdditionalData additional_data)
2342 {
2343  std::vector<const DoFHandlerType *> dof_handlers;
2344  std::vector<const ConstraintMatrix *> constraints;
2345  std::vector<QuadratureType> quads;
2346 
2347  dof_handlers.push_back(&dof_handler);
2348  constraints.push_back (&constraints_in);
2349  quads.push_back (quad);
2350 
2351  std::vector<IndexSet> locally_owned_sets =
2352  internal::MatrixFreeImplementation::extract_locally_owned_index_sets
2353  (dof_handlers, additional_data.level_mg_handler);
2354 
2355  std::vector<hp::QCollection<1> > quad_hp;
2356  quad_hp.emplace_back (quad);
2357 
2358  internal_reinit(StaticMappingQ1<dim>::mapping, dof_handlers,constraints,
2359  locally_owned_sets, quad_hp, additional_data);
2360 }
2361 
2362 
2363 
2364 template <int dim, typename Number>
2365 template <typename DoFHandlerType, typename QuadratureType>
2366 void MatrixFree<dim,Number>::
2367 reinit(const Mapping<dim> &mapping,
2368  const DoFHandlerType &dof_handler,
2369  const ConstraintMatrix &constraints_in,
2370  const QuadratureType &quad,
2371  const typename MatrixFree<dim,Number>::AdditionalData additional_data)
2372 {
2373  std::vector<const DoFHandlerType *> dof_handlers;
2374  std::vector<const ConstraintMatrix *> constraints;
2375 
2376  dof_handlers.push_back(&dof_handler);
2377  constraints.push_back (&constraints_in);
2378 
2379  std::vector<IndexSet> locally_owned_sets =
2380  internal::MatrixFreeImplementation::extract_locally_owned_index_sets
2381  (dof_handlers, additional_data.level_mg_handler);
2382 
2383  std::vector<hp::QCollection<1> > quad_hp;
2384  quad_hp.emplace_back (quad);
2385 
2386  internal_reinit(mapping, dof_handlers,constraints,locally_owned_sets,
2387  quad_hp, additional_data);
2388 }
2389 
2390 
2391 
2392 template <int dim, typename Number>
2393 template <typename DoFHandlerType, typename QuadratureType>
2394 void MatrixFree<dim,Number>::
2395 reinit(const std::vector<const DoFHandlerType *> &dof_handler,
2396  const std::vector<const ConstraintMatrix *> &constraint,
2397  const std::vector<QuadratureType> &quad,
2398  const typename MatrixFree<dim,Number>::AdditionalData additional_data)
2399 {
2400  std::vector<IndexSet> locally_owned_set =
2401  internal::MatrixFreeImplementation::extract_locally_owned_index_sets
2402  (dof_handler, additional_data.level_mg_handler);
2403  std::vector<hp::QCollection<1> > quad_hp;
2404  for (unsigned int q=0; q<quad.size(); ++q)
2405  quad_hp.emplace_back (quad[q]);
2406  internal_reinit(StaticMappingQ1<dim>::mapping, dof_handler,constraint,
2407  locally_owned_set, quad_hp, additional_data);
2408 }
2409 
2410 
2411 
2412 template <int dim, typename Number>
2413 template <typename DoFHandlerType, typename QuadratureType>
2414 void MatrixFree<dim,Number>::
2415 reinit(const std::vector<const DoFHandlerType *> &dof_handler,
2416  const std::vector<const ConstraintMatrix *> &constraint,
2417  const QuadratureType &quad,
2418  const typename MatrixFree<dim,Number>::AdditionalData additional_data)
2419 {
2420  std::vector<IndexSet> locally_owned_set =
2421  internal::MatrixFreeImplementation::extract_locally_owned_index_sets
2422  (dof_handler, additional_data.level_mg_handler);
2423  std::vector<hp::QCollection<1> > quad_hp;
2424  quad_hp.emplace_back (quad);
2425  internal_reinit(StaticMappingQ1<dim>::mapping, dof_handler,constraint,
2426  locally_owned_set, quad_hp, additional_data);
2427 }
2428 
2429 
2430 
2431 template <int dim, typename Number>
2432 template <typename DoFHandlerType, typename QuadratureType>
2433 void MatrixFree<dim,Number>::
2434 reinit(const Mapping<dim> &mapping,
2435  const std::vector<const DoFHandlerType *> &dof_handler,
2436  const std::vector<const ConstraintMatrix *> &constraint,
2437  const QuadratureType &quad,
2438  const typename MatrixFree<dim,Number>::AdditionalData additional_data)
2439 {
2440  std::vector<IndexSet> locally_owned_set =
2441  internal::MatrixFreeImplementation::extract_locally_owned_index_sets
2442  (dof_handler, additional_data.level_mg_handler);
2443  std::vector<hp::QCollection<1> > quad_hp;
2444  quad_hp.emplace_back (quad);
2445  internal_reinit(mapping, dof_handler,constraint,
2446  locally_owned_set, quad_hp, additional_data);
2447 }
2448 
2449 
2450 
2451 template <int dim, typename Number>
2452 template <typename DoFHandlerType, typename QuadratureType>
2453 void MatrixFree<dim,Number>::
2454 reinit(const Mapping<dim> &mapping,
2455  const std::vector<const DoFHandlerType *> &dof_handler,
2456  const std::vector<const ConstraintMatrix *> &constraint,
2457  const std::vector<QuadratureType> &quad,
2458  const typename MatrixFree<dim,Number>::AdditionalData additional_data)
2459 {
2460  std::vector<IndexSet> locally_owned_set =
2461  internal::MatrixFreeImplementation::extract_locally_owned_index_sets
2462  (dof_handler, additional_data.level_mg_handler);
2463  std::vector<hp::QCollection<1> > quad_hp;
2464  for (unsigned int q=0; q<quad.size(); ++q)
2465  quad_hp.emplace_back (quad[q]);
2466  internal_reinit(mapping, dof_handler,constraint,locally_owned_set,
2467  quad_hp, additional_data);
2468 }
2469 
2470 
2471 
2472 template <int dim, typename Number>
2473 template <typename DoFHandlerType, typename QuadratureType>
2474 void MatrixFree<dim,Number>::
2475 reinit(const Mapping<dim> &mapping,
2476  const std::vector<const DoFHandlerType *> &dof_handler,
2477  const std::vector<const ConstraintMatrix *> &constraint,
2478  const std::vector<IndexSet> &locally_owned_set,
2479  const std::vector<QuadratureType> &quad,
2480  const typename MatrixFree<dim,Number>::AdditionalData additional_data)
2481 {
2482  // find out whether we use a hp Quadrature or a standard quadrature
2483  std::vector<hp::QCollection<1> > quad_hp;
2484  for (unsigned int q=0; q<quad.size(); ++q)
2485  quad_hp.emplace_back (quad[q]);
2486  internal_reinit (mapping,
2487  dof_handler,
2488  constraint, locally_owned_set, quad_hp, additional_data);
2489 }
2490 
2491 
2492 
2493 // ------------------------------ implementation of loops --------------------
2494 
2495 // internal helper functions that define how to call MPI data exchange
2496 // functions: for generic vectors, do nothing at all. For distributed vectors,
2497 // call update_ghost_values_start function and so on. If we have collections
2498 // of vectors, just do the individual functions of the components. In order to
2499 // keep ghost values consistent (whether we are in read or write mode), we
2500 // also reset the values at the end. the whole situation is a bit complicated
2501 // by the fact that we need to treat block vectors differently, which use some
2502 // additional helper functions to select the blocks and template magic.
2503 namespace internal
2504 {
2505  template <int dim, typename Number>
2506  struct VectorDataExchange
2507  {
2508  // An arbitrary shift for communication to reduce the risk for accidental
2509  // interaction with other open communications that a user program might
2510  // set up
2511  static constexpr unsigned int channel_shift = 103;
2512 
2513  VectorDataExchange (const ::MatrixFree<dim,Number> &matrix_free,
2514  const typename ::MatrixFree<dim,Number>::DataAccessOnFaces vector_face_access,
2515  const unsigned int n_components)
2516  :
2517  matrix_free (matrix_free),
2518  vector_face_access (matrix_free.get_task_info().face_partition_data.empty() ?
2519  ::MatrixFree<dim,Number>::DataAccessOnFaces::unspecified :
2520  vector_face_access),
2521  ghosts_were_set (false)
2522 #ifdef DEAL_II_WITH_MPI
2523  , tmp_data(n_components),
2524  requests(n_components)
2525 #endif
2526  {
2527  (void)n_components;
2528  if (this->vector_face_access != ::MatrixFree<dim,Number>::DataAccessOnFaces::unspecified)
2529  for (unsigned int c=0; c<matrix_free.n_components(); ++c)
2530  AssertDimension(matrix_free.get_dof_info(c).vector_partitioner_face_variants.size(), 3);
2531  }
2532 
2533  ~VectorDataExchange ()
2534  {
2535 #ifdef DEAL_II_WITH_MPI
2536  for (unsigned int i=0; i<tmp_data.size(); ++i)
2537  if (tmp_data[i] != nullptr)
2538  matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
2539 #endif
2540  }
2541 
2542  unsigned int find_vector_in_mf (const LinearAlgebra::distributed::Vector<Number> &vec,
2543  const bool check_global_compatibility = true) const
2544  {
2545  unsigned int mf_component = numbers::invalid_unsigned_int;
2546  (void)check_global_compatibility;
2547  for (unsigned int c=0; c<matrix_free.n_components(); ++c)
2548  if (
2549 #ifdef DEBUG
2550  check_global_compatibility
2551  ?
2552  vec.get_partitioner()->is_globally_compatible(*matrix_free.get_dof_info(c).vector_partitioner)
2553  :
2554 #endif
2555  vec.get_partitioner()->is_compatible(*matrix_free.get_dof_info(c).vector_partitioner))
2556  {
2557  mf_component = c;
2558  break;
2559  }
2560  return mf_component;
2561  }
2562 
2563  const Utilities::MPI::Partitioner &
2564  get_partitioner(const unsigned int mf_component) const
2565  {
2566  AssertDimension(matrix_free.get_dof_info(mf_component).vector_partitioner_face_variants.size(),3);
2567  if (vector_face_access == ::MatrixFree<dim,Number>::DataAccessOnFaces::none)
2568  return *matrix_free.get_dof_info(mf_component).vector_partitioner_face_variants[0];
2569  else if (vector_face_access == ::MatrixFree<dim,Number>::DataAccessOnFaces::values)
2570  return *matrix_free.get_dof_info(mf_component).vector_partitioner_face_variants[1];
2571  else
2572  return *matrix_free.get_dof_info(mf_component).vector_partitioner_face_variants[2];
2573  }
2574 
2575  void update_ghost_values_start(const unsigned int component_in_block_vector,
2576  const LinearAlgebra::distributed::Vector<Number> &vec)
2577  {
2578  (void)component_in_block_vector;
2579  bool ghosts_set = vec.has_ghost_elements();
2580  if (ghosts_set)
2581  ghosts_were_set = true;
2582  if (vector_face_access == ::MatrixFree<dim,Number>::DataAccessOnFaces::unspecified ||
2583  vec.size() == 0)
2584  vec.update_ghost_values_start(component_in_block_vector + channel_shift);
2585  else
2586  {
2587 #ifdef DEAL_II_WITH_MPI
2588  const unsigned int mf_component = find_vector_in_mf(vec);
2589  if (&get_partitioner(mf_component) == matrix_free.get_dof_info(mf_component)
2590  .vector_partitioner.get())
2591  {
2592  vec.update_ghost_values_start(component_in_block_vector + channel_shift);
2593  return;
2594  }
2595 
2596  const Utilities::MPI::Partitioner &part = get_partitioner(mf_component);
2597  if (part.n_ghost_indices()==0 && part.n_import_indices()==0)
2598  return;
2599 
2600  tmp_data[component_in_block_vector] = matrix_free.acquire_scratch_data_non_threadsafe();
2601  tmp_data[component_in_block_vector]->resize_fast(part.n_import_indices());
2602  AssertDimension(requests.size(), tmp_data.size());
2603 
2604  part.export_to_ghosted_array_start
2605  (component_in_block_vector+channel_shift,
2606  ArrayView<const Number>(vec.begin(), part.local_size()),
2607  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
2608  part.n_import_indices()),
2609  ArrayView<Number>(const_cast<Number *>(vec.begin()) +
2610  vec.get_partitioner()->local_size(),
2611  vec.get_partitioner()->n_ghost_indices()),
2612  this->requests[component_in_block_vector]);
2613 #endif
2614  }
2615  }
2616 
2617  void update_ghost_values_finish (const unsigned int component_in_block_vector,
2618  const LinearAlgebra::distributed::Vector<Number> &vec)
2619  {
2620  (void)component_in_block_vector;
2621  if (vector_face_access == ::MatrixFree<dim,Number>::DataAccessOnFaces::unspecified ||
2622  vec.size() == 0)
2623  vec.update_ghost_values_finish();
2624  else
2625  {
2626 #ifdef DEAL_II_WITH_MPI
2627 
2628  AssertIndexRange(component_in_block_vector, tmp_data.size());
2629  AssertDimension(requests.size(), tmp_data.size());
2630 
2631  const unsigned int mf_component = find_vector_in_mf(vec);
2632  const Utilities::MPI::Partitioner &part = get_partitioner(mf_component);
2633  if (&part == matrix_free.get_dof_info(mf_component).vector_partitioner.get())
2634  {
2635  vec.update_ghost_values_finish();
2636  return;
2637  }
2638 
2639  if (part.n_ghost_indices()==0 && part.n_import_indices()==0)
2640  return;
2641 
2642  part.export_to_ghosted_array_finish
2643  (ArrayView<Number>(const_cast<Number *>(vec.begin()) +
2644  vec.get_partitioner()->local_size(),
2645  vec.get_partitioner()->n_ghost_indices()),
2646  this->requests[component_in_block_vector]);
2647 
2648  matrix_free.release_scratch_data_non_threadsafe(tmp_data[component_in_block_vector]);
2649  tmp_data[component_in_block_vector] = nullptr;
2650 #endif
2651  }
2652  }
2653 
2654  void compress_start(const unsigned int component_in_block_vector,
2655  LinearAlgebra::distributed::Vector<Number> &vec)
2656  {
2657  (void)component_in_block_vector;
2658  Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
2659  if (vector_face_access == ::MatrixFree<dim,Number>::DataAccessOnFaces::unspecified ||
2660  vec.size() == 0)
2661  vec.compress_start(component_in_block_vector + channel_shift);
2662  else
2663  {
2664 #ifdef DEAL_II_WITH_MPI
2665 
2666  const unsigned int mf_component = find_vector_in_mf(vec);
2667  const Utilities::MPI::Partitioner &part = get_partitioner(mf_component);
2668  if (&part == matrix_free.get_dof_info(mf_component).vector_partitioner.get())
2669  {
2670  vec.compress_start(component_in_block_vector + channel_shift);
2671  return;
2672  }
2673 
2674  if (part.n_ghost_indices()==0 && part.n_import_indices()==0)
2675  return;
2676 
2677  tmp_data[component_in_block_vector] = matrix_free.acquire_scratch_data_non_threadsafe();
2678  tmp_data[component_in_block_vector]->resize_fast(part.n_import_indices());
2679  AssertDimension(requests.size(), tmp_data.size());
2680 
2681  part.import_from_ghosted_array_start
2682  (::VectorOperation::add,
2683  component_in_block_vector+channel_shift,
2684  ArrayView<Number>(vec.begin()+vec.get_partitioner()->local_size(),
2685  vec.get_partitioner()->n_ghost_indices()),
2686  ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
2687  part.n_import_indices()),
2688  this->requests[component_in_block_vector]);
2689 #endif
2690  }
2691  }
2692 
2693  void compress_finish (const unsigned int component_in_block_vector,
2694  LinearAlgebra::distributed::Vector<Number> &vec)
2695  {
2696  (void)component_in_block_vector;
2697  if (vector_face_access == ::MatrixFree<dim,Number>::DataAccessOnFaces::unspecified ||
2698  vec.size() == 0)
2699  vec.compress_finish(::VectorOperation::add);
2700  else
2701  {
2702 #ifdef DEAL_II_WITH_MPI
2703  AssertIndexRange(component_in_block_vector, tmp_data.size());
2704  AssertDimension(requests.size(), tmp_data.size());
2705 
2706  const unsigned int mf_component = find_vector_in_mf(vec);
2707 
2708  const Utilities::MPI::Partitioner &part = get_partitioner(mf_component);
2709  if (&part == matrix_free.get_dof_info(mf_component).vector_partitioner.get())
2710  {
2711  vec.compress_finish(::VectorOperation::add);
2712  return;
2713  }
2714 
2715  if (part.n_ghost_indices()==0 && part.n_import_indices()==0)
2716  return;
2717 
2718  part.import_from_ghosted_array_finish
2719  (VectorOperation::add,
2720  ArrayView<const Number>(tmp_data[component_in_block_vector]->begin(),
2721  part.n_import_indices()),
2722  ArrayView<Number>(vec.begin(), part.local_size()),
2723  ArrayView<Number>(vec.begin()+vec.get_partitioner()->local_size(),
2724  vec.get_partitioner()->n_ghost_indices()),
2725  this->requests[component_in_block_vector]);
2726 
2727  matrix_free.release_scratch_data_non_threadsafe(tmp_data[component_in_block_vector]);
2728  tmp_data[component_in_block_vector] = nullptr;
2729 #endif
2730  }
2731  }
2732 
2733  void reset_ghost_values(const LinearAlgebra::distributed::Vector<Number> &vec) const
2734  {
2735  if (ghosts_were_set == true)
2736  return;
2737 
2738  if (vector_face_access == ::MatrixFree<dim,Number>::DataAccessOnFaces::unspecified ||
2739  vec.size() == 0)
2740  vec.zero_out_ghosts();
2741  else
2742  {
2743 #ifdef DEAL_II_WITH_MPI
2744  AssertDimension(requests.size(), tmp_data.size());
2745 
2746  const unsigned int mf_component = find_vector_in_mf(vec);
2747  const Utilities::MPI::Partitioner &part = get_partitioner(mf_component);
2748  if (&part == matrix_free.get_dof_info(mf_component).vector_partitioner.get())
2749  vec.zero_out_ghosts();
2750  else if (part.n_ghost_indices() > 0)
2751  {
2752  for (std::vector<std::pair<unsigned int, unsigned int> >::const_iterator
2753  my_ghosts = part.ghost_indices_within_larger_ghost_set().begin();
2754  my_ghosts != part.ghost_indices_within_larger_ghost_set().end();
2755  ++my_ghosts)
2756  for (unsigned int j=my_ghosts->first; j<my_ghosts->second; j++)
2757  {
2758  const_cast<LinearAlgebra::distributed::Vector<Number> &>(vec)
2759  .local_element(j+part.local_size()) = 0.;
2760  }
2761  }
2762 #endif
2763  }
2764  }
2765 
2766  void zero_vector_region(const unsigned int range_index,
2767  LinearAlgebra::distributed::Vector<Number> &vec) const
2768  {
2769  if (range_index == numbers::invalid_unsigned_int)
2770  vec = 0;
2771  else
2772  {
2773  const unsigned int mf_component = find_vector_in_mf(vec, false);
2774  const internal::MatrixFreeFunctions::DoFInfo &dof_info =
2775  matrix_free.get_dof_info(mf_component);
2776  Assert(dof_info.vector_zero_range_list_index.empty() == false,
2777  ExcNotInitialized());
2778 
2779  Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
2780  ExcInternalError());
2781  AssertIndexRange(range_index, dof_info.vector_zero_range_list_index.size()-1);
2782  for (unsigned int id=dof_info.vector_zero_range_list_index[range_index];
2783  id != dof_info.vector_zero_range_list_index[range_index+1]; ++id)
2784  {
2785  const unsigned int start_pos = dof_info.vector_zero_range_list[id]*
2786  internal::MatrixFreeFunctions::DoFInfo::chunk_size_zero_vector;
2787  const unsigned int end_pos = std::min((dof_info.vector_zero_range_list[id]+1)*
2788  internal::MatrixFreeFunctions::DoFInfo::chunk_size_zero_vector,
2789  dof_info.vector_partitioner->local_size()+
2790  dof_info.vector_partitioner->n_ghost_indices());
2791  std::memset(vec.begin()+start_pos, 0, (end_pos-start_pos)*sizeof(Number));
2792  }
2793  }
2794  }
2795 
2796  const ::MatrixFree<dim,Number> &matrix_free;
2797  const typename ::MatrixFree<dim,Number>::DataAccessOnFaces vector_face_access;
2798  bool ghosts_were_set;
2799 #ifdef DEAL_II_WITH_MPI
2800  std::vector<AlignedVector<Number> *> tmp_data;
2801  std::vector<std::vector<MPI_Request> > requests;
2802 #endif
2803  };
2804 
2805  template <typename VectorStruct>
2806  unsigned int n_components (const VectorStruct &vec);
2807 
2808  template <typename VectorStruct>
2809  unsigned int n_components_block (const VectorStruct &vec,
2810  std::integral_constant<bool,true>)
2811  {
2812  unsigned int components = 0;
2813  for (unsigned int bl=0; bl<vec.n_blocks(); ++bl)
2814  components += n_components(vec.block(bl));
2815  return components;
2816  }
2817 
2818  template <typename VectorStruct>
2819  unsigned int n_components_block (const VectorStruct &,
2820  std::integral_constant<bool,false>)
2821  {
2822  return 1;
2823  }
2824 
2825  template <typename VectorStruct>
2826  unsigned int n_components (const VectorStruct &vec)
2827  {
2828  return n_components_block(vec, std::integral_constant<bool,IsBlockVector<VectorStruct>::value>());
2829  }
2830 
2831  template <typename VectorStruct>
2832  inline
2833  unsigned int n_components (const std::vector<VectorStruct> &vec)
2834  {
2835  unsigned int components = 0;
2836  for (unsigned int comp=0; comp<vec.size(); comp++)
2837  components += n_components_block(vec[comp], std::integral_constant<bool,IsBlockVector<VectorStruct>::value>());
2838  return components;
2839  }
2840 
2841  template <typename VectorStruct>
2842  inline
2843  unsigned int n_components (const std::vector<VectorStruct *> &vec)
2844  {
2845  unsigned int components = 0;
2846  for (unsigned int comp=0; comp<vec.size(); comp++)
2847  components += n_components_block(*vec[comp], std::integral_constant<bool,IsBlockVector<VectorStruct>::value>());
2848  return components;
2849  }
2850 
2851  template <int dim, typename VectorStruct, typename Number>
2852  void update_ghost_values_start_block (const VectorStruct &vec,
2853  const unsigned int channel,
2854  std::integral_constant<bool,true>,
2855  VectorDataExchange<dim,Number> &exchanger);
2856  template <int dim, typename VectorStruct, typename Number>
2857  void reset_ghost_values_block (const VectorStruct &vec,
2858  std::integral_constant<bool,true>,
2859  VectorDataExchange<dim,Number> &exchanger);
2860  template <int dim, typename VectorStruct, typename Number>
2861  void update_ghost_values_finish_block (const VectorStruct &vec,
2862  const unsigned int channel,
2863  std::integral_constant<bool,true>,
2864  VectorDataExchange<dim,Number> &exchanger);
2865  template <int dim, typename VectorStruct, typename Number>
2866  void compress_start_block (const VectorStruct &vec,
2867  const unsigned int channel,
2868  std::integral_constant<bool,true>,
2869  VectorDataExchange<dim,Number> &exchanger);
2870  template <int dim, typename VectorStruct, typename Number>
2871  void compress_finish_block (const VectorStruct &vec,
2872  const unsigned int channel,
2873  std::integral_constant<bool,true>,
2874  VectorDataExchange<dim,Number> &exchanger);
2875  template <int dim, typename VectorStruct, typename Number>
2876  void zero_vector_region_block (const unsigned int range_index,
2877  VectorStruct &,
2878  std::integral_constant<bool,true>,
2879  VectorDataExchange<dim,Number> &);
2880 
2881  template <int dim, typename VectorStruct, typename Number>
2882  void update_ghost_values_start_block (const VectorStruct &,
2883  const unsigned int ,
2884  std::integral_constant<bool,false>,
2885  VectorDataExchange<dim,Number> &)
2886  {}
2887  template <int dim, typename VectorStruct, typename Number>
2888  void reset_ghost_values_block (const VectorStruct &,
2889  std::integral_constant<bool,false>,
2890  VectorDataExchange<dim,Number> &)
2891  {}
2892  template <int dim, typename VectorStruct, typename Number>
2893  void update_ghost_values_finish_block (const VectorStruct &,
2894  const unsigned int ,
2895  std::integral_constant<bool,false>,
2896  VectorDataExchange<dim,Number> &)
2897  {}
2898  template <int dim, typename VectorStruct, typename Number>
2899  void compress_start_block (const VectorStruct &,
2900  const unsigned int ,
2901  std::integral_constant<bool,false>,
2902  VectorDataExchange<dim,Number> &)
2903  {}
2904  template <int dim, typename VectorStruct, typename Number>
2905  void compress_finish_block (const VectorStruct &,
2906  const unsigned int ,
2907  std::integral_constant<bool,false>,
2908  VectorDataExchange<dim,Number> &)
2909  {}
2910  template <int dim, typename VectorStruct, typename Number>
2911  void zero_vector_region_block (const unsigned int range_index,
2912  VectorStruct &vec,
2913  std::integral_constant<bool,false>,
2914  VectorDataExchange<dim,Number> &)
2915  {
2916  if (range_index == 0 || range_index == numbers::invalid_unsigned_int)
2917  vec = 0;
2918  }
2919 
2920 
2921 
2922  template <int dim, typename VectorStruct, typename Number>
2923  inline
2924  void update_ghost_values_start (const VectorStruct &vec,
2925  VectorDataExchange<dim,Number> &exchanger,
2926  const unsigned int channel = 0)
2927  {
2928  update_ghost_values_start_block(vec, channel,
2929  std::integral_constant<bool,
2930  IsBlockVector<VectorStruct>::value>(),
2931  exchanger);
2932  }
2933 
2934 
2935 
2936  template <int dim, typename Number, typename Number2>
2937  inline
2938  void update_ghost_values_start (const LinearAlgebra::distributed::Vector<Number> &vec,
2939  VectorDataExchange<dim,Number2> &exchanger,
2940  const unsigned int channel = 0)
2941  {
2942  exchanger.update_ghost_values_start(channel, vec);
2943  }
2944 
2945 
2946 
2947  template <int dim, typename VectorStruct, typename Number>
2948  inline
2949  void update_ghost_values_start (const std::vector<VectorStruct> &vec,
2950  VectorDataExchange<dim,Number> &exchanger)
2951  {
2952  unsigned int component_index = 0;
2953  for (unsigned int comp=0; comp<vec.size(); comp++)
2954  {
2955  update_ghost_values_start(vec[comp], exchanger, component_index);
2956  component_index += n_components(vec[comp]);
2957  }
2958  }
2959 
2960 
2961 
2962  template <int dim, typename VectorStruct, typename Number>
2963  inline
2964  void update_ghost_values_start (const std::vector<VectorStruct *> &vec,
2965  VectorDataExchange<dim,Number> &exchanger)
2966  {
2967  unsigned int component_index = 0;
2968  for (unsigned int comp=0; comp<vec.size(); comp++)
2969  {
2970  update_ghost_values_start(*vec[comp], exchanger, component_index);
2971  component_index += n_components(*vec[comp]);
2972  }
2973  }
2974 
2975 
2976 
2977  template <int dim, typename VectorStruct, typename Number>
2978  inline
2979  void update_ghost_values_start_block (const VectorStruct &vec,
2980  const unsigned int channel,
2981  std::integral_constant<bool,true>,
2982  VectorDataExchange<dim,Number> &exchanger)
2983  {
2984  for (unsigned int i=0; i<vec.n_blocks(); ++i)
2985  update_ghost_values_start(vec.block(i), exchanger, channel+i);
2986  }
2987 
2988 
2989 
2990  // if the input vector did not have ghosts imported, clear them here again
2991  // in order to avoid subsequent operations e.g. in linear solvers to work
2992  // with ghosts all the time
2993  template <int dim, typename VectorStruct, typename Number>
2994  inline
2995  void reset_ghost_values (const VectorStruct &vec,
2996  VectorDataExchange<dim,Number> &exchanger)
2997  {
2998  reset_ghost_values_block(vec,
2999  std::integral_constant<bool,
3000  IsBlockVector<VectorStruct>::value>(),
3001  exchanger);
3002  }
3003 
3004 
3005 
3006  template <int dim, typename Number, typename Number2>
3007  inline
3008  void reset_ghost_values (const LinearAlgebra::distributed::Vector<Number> &vec,
3009  VectorDataExchange<dim,Number2> &exchanger)
3010  {
3011  exchanger.reset_ghost_values(vec);
3012  }
3013 
3014 
3015 
3016  template <int dim, typename VectorStruct, typename Number>
3017  inline
3018  void reset_ghost_values (const std::vector<VectorStruct> &vec,
3019  VectorDataExchange<dim,Number> &exchanger)
3020  {
3021  for (unsigned int comp=0; comp<vec.size(); comp++)
3022  reset_ghost_values(vec[comp], exchanger);
3023  }
3024 
3025 
3026 
3027  template <int dim, typename VectorStruct, typename Number>
3028  inline
3029  void reset_ghost_values (const std::vector<VectorStruct *> &vec,
3030  VectorDataExchange<dim,Number> &exchanger)
3031  {
3032  for (unsigned int comp=0; comp<vec.size(); comp++)
3033  reset_ghost_values(*vec[comp], exchanger);
3034  }
3035 
3036 
3037 
3038  template <int dim, typename VectorStruct, typename Number>
3039  inline
3040  void reset_ghost_values_block (const VectorStruct &vec,
3041  std::integral_constant<bool,true>,
3042  VectorDataExchange<dim,Number> &exchanger)
3043  {
3044  for (unsigned int i=0; i<vec.n_blocks(); ++i)
3045  reset_ghost_values(vec.block(i), exchanger);
3046  }
3047 
3048 
3049 
3050  template <int dim, typename VectorStruct, typename Number>
3051  inline
3052  void update_ghost_values_finish (const VectorStruct &vec,
3053  VectorDataExchange<dim,Number> &exchanger,
3054  const unsigned int channel = 0)
3055  {
3056  update_ghost_values_finish_block(vec, channel,
3057  std::integral_constant<bool,
3058  IsBlockVector<VectorStruct>::value>(),
3059  exchanger);
3060  }
3061 
3062 
3063 
3064  template <int dim, typename Number, typename Number2>
3065  inline
3066  void update_ghost_values_finish (const LinearAlgebra::distributed::Vector<Number> &vec,
3067  VectorDataExchange<dim,Number2> &exchanger,
3068  const unsigned int channel = 0)
3069  {
3070  exchanger.update_ghost_values_finish(channel, vec);
3071  }
3072 
3073 
3074 
3075  template <int dim, typename VectorStruct, typename Number>
3076  inline
3077  void update_ghost_values_finish (const std::vector<VectorStruct> &vec,
3078  VectorDataExchange<dim,Number> &exchanger)
3079  {
3080  unsigned int component_index = 0;
3081  for (unsigned int comp=0; comp<vec.size(); comp++)
3082  {
3083  update_ghost_values_finish(vec[comp], exchanger, component_index);
3084  component_index += n_components(vec[comp]);
3085  }
3086  }
3087 
3088 
3089 
3090  template <int dim, typename VectorStruct, typename Number>
3091  inline
3092  void update_ghost_values_finish (const std::vector<VectorStruct *> &vec,
3093  VectorDataExchange<dim,Number> &exchanger)
3094  {
3095  unsigned int component_index = 0;
3096  for (unsigned int comp=0; comp<vec.size(); comp++)
3097  {
3098  update_ghost_values_finish(*vec[comp], exchanger, component_index);
3099  component_index += n_components(*vec[comp]);
3100  }
3101  }
3102 
3103 
3104 
3105  template <int dim, typename VectorStruct, typename Number>
3106  inline
3107  void update_ghost_values_finish_block (const VectorStruct &vec,
3108  const unsigned int channel,
3109  std::integral_constant<bool,true>,
3110  VectorDataExchange<dim,Number> &exchanger)
3111  {
3112  for (unsigned int i=0; i<vec.n_blocks(); ++i)
3113  update_ghost_values_finish(vec.block(i), exchanger, channel+i);
3114  }
3115 
3116 
3117 
3118  template <int dim, typename VectorStruct, typename Number>
3119  inline
3120  void compress_start (VectorStruct &vec,
3121  VectorDataExchange<dim, Number> &exchanger,
3122  const unsigned int channel = 0)
3123  {
3124  compress_start_block (vec, channel,
3125  std::integral_constant<bool,
3126  IsBlockVector<VectorStruct>::value>(),
3127  exchanger);
3128  }
3129 
3130 
3131 
3132  template <int dim, typename Number, typename Number2>
3133  inline
3134  void compress_start (LinearAlgebra::distributed::Vector<Number> &vec,
3135  VectorDataExchange<dim,Number2> &exchanger,
3136  const unsigned int channel = 0)
3137  {
3138  exchanger.compress_start(channel, vec);
3139  }
3140 
3141 
3142 
3143  template <int dim, typename VectorStruct, typename Number>
3144  inline
3145  void compress_start (std::vector<VectorStruct> &vec,
3146  VectorDataExchange<dim, Number> &exchanger)
3147  {
3148  unsigned int component_index = 0;
3149  for (unsigned int comp=0; comp<vec.size(); comp++)
3150  {
3151  compress_start(vec[comp], exchanger, component_index);
3152  component_index += n_components(vec[comp]);
3153  }
3154  }
3155 
3156 
3157 
3158  template <int dim, typename VectorStruct, typename Number>
3159  inline
3160  void compress_start (std::vector<VectorStruct *> &vec,
3161  VectorDataExchange<dim, Number> &exchanger)
3162  {
3163  unsigned int component_index = 0;
3164  for (unsigned int comp=0; comp<vec.size(); comp++)
3165  {
3166  compress_start(*vec[comp], exchanger, component_index);
3167  component_index += n_components(*vec[comp]);
3168  }
3169  }
3170 
3171 
3172 
3173  template <int dim, typename VectorStruct, typename Number>
3174  inline
3175  void compress_start_block (VectorStruct &vec,
3176  const unsigned int channel,
3177  std::integral_constant<bool,true>,
3178  VectorDataExchange<dim, Number> &exchanger)
3179  {
3180  for (unsigned int i=0; i<vec.n_blocks(); ++i)
3181  compress_start(vec.block(i), exchanger, channel+i);
3182  }
3183 
3184 
3185 
3186  template <int dim, typename VectorStruct, typename Number>
3187  inline
3188  void compress_finish (VectorStruct &vec,
3189  VectorDataExchange<dim, Number> &exchanger,
3190  const unsigned int channel = 0)
3191  {
3192  compress_finish_block(vec, channel,
3193  std::integral_constant<bool,
3194  IsBlockVector<VectorStruct>::value>(),
3195  exchanger);
3196  }
3197 
3198 
3199 
3200  template <int dim, typename Number, typename Number2>
3201  inline
3202  void compress_finish (LinearAlgebra::distributed::Vector<Number> &vec,
3203  VectorDataExchange<dim, Number2> &exchanger,
3204  const unsigned int channel = 0)
3205  {
3206  exchanger.compress_finish(channel, vec);
3207  }
3208 
3209 
3210 
3211  template <int dim, typename VectorStruct, typename Number>
3212  inline
3213  void compress_finish (std::vector<VectorStruct> &vec,
3214  VectorDataExchange<dim, Number> &exchanger)
3215  {
3216  unsigned int component_index = 0;
3217  for (unsigned int comp=0; comp<vec.size(); comp++)
3218  {
3219  compress_finish(vec[comp], exchanger, component_index);
3220  component_index += n_components(vec[comp]);
3221  }
3222  }
3223 
3224 
3225 
3226  template <int dim, typename VectorStruct, typename Number>
3227  inline
3228  void compress_finish (std::vector<VectorStruct *> &vec,
3229  VectorDataExchange<dim, Number> &exchanger)
3230  {
3231  unsigned int component_index = 0;
3232  for (unsigned int comp=0; comp<vec.size(); comp++)
3233  {
3234  compress_finish(*vec[comp], exchanger, component_index);
3235  component_index += n_components(*vec[comp]);
3236  }
3237  }
3238 
3239 
3240 
3241  template <int dim, typename VectorStruct, typename Number>
3242  inline
3243  void compress_finish_block (VectorStruct &vec,
3244  const unsigned int channel,
3245  std::integral_constant<bool,true>,
3246  VectorDataExchange<dim, Number> &exchanger)
3247  {
3248  for (unsigned int i=0; i<vec.n_blocks(); ++i)
3249  compress_finish(vec.block(i), exchanger, channel+i);
3250  }
3251 
3252 
3253 
3254  template <int dim, typename VectorStruct, typename Number>
3255  inline
3256  void zero_vector_region (const unsigned int range_index,
3257  VectorStruct &vec,
3258  VectorDataExchange<dim, Number> &exchanger)
3259  {
3260  zero_vector_region_block(range_index, vec,
3261  std::integral_constant<bool,
3262  IsBlockVector<VectorStruct>::value>(),
3263  exchanger);
3264  }
3265 
3266 
3267 
3268  template <int dim, typename Number, typename Number2>
3269  inline
3270  void zero_vector_region (const unsigned int range_index,
3271  LinearAlgebra::distributed::Vector<Number> &vec,
3272  VectorDataExchange<dim, Number2> &exchanger)
3273  {
3274  exchanger.zero_vector_region(range_index, vec);
3275  }
3276 
3277 
3278 
3279  template <int dim, typename VectorStruct, typename Number>
3280  inline
3281  void zero_vector_region (const unsigned int range_index,
3282  std::vector<VectorStruct> &vec,
3283  VectorDataExchange<dim, Number> &exchanger)
3284  {
3285  for (unsigned int comp=0; comp<vec.size(); comp++)
3286  zero_vector_region(range_index, vec[comp], exchanger);
3287  }
3288 
3289 
3290 
3291  template <int dim, typename VectorStruct, typename Number>
3292  inline
3293  void zero_vector_region (const unsigned int range_index,
3294  std::vector<VectorStruct *> &vec,
3295  VectorDataExchange<dim, Number> &exchanger)
3296  {
3297  for (unsigned int comp=0; comp<vec.size(); comp++)
3298  zero_vector_region(range_index, *vec[comp], exchanger);
3299  }
3300 
3301 
3302 
3303  template <int dim, typename VectorStruct, typename Number>
3304  inline
3305  void zero_vector_region_block (const unsigned int range_index,
3306  VectorStruct &vec,
3307  std::integral_constant<bool,true>,
3308  VectorDataExchange<dim, Number> &exchanger)
3309  {
3310  for (unsigned int i=0; i<vec.n_blocks(); ++i)
3311  zero_vector_region(range_index, vec.block(i), exchanger);
3312  }
3313 
3314 
3315 
3316  namespace MatrixFreeFunctions
3317  {
3318  // struct to select between a const interface and a non-const interface
3319  // for MFWorker
3320  template <typename, typename, typename, typename, bool>
3321  struct InterfaceSelector
3322  {};
3323 
3324  // Version of constant functions
3325  template <typename MF, typename InVector, typename OutVector, typename Container>
3326  struct InterfaceSelector<MF, InVector, OutVector, Container, true>
3327  {
3328  typedef void (Container::*function_type)
3329  (const MF &, OutVector &, const InVector &,
3330  const std::pair<unsigned int, unsigned int> &)const;
3331  };
3332 
3333  // Version for non-constant functions
3334  template <typename MF, typename InVector, typename OutVector, typename Container>
3335  struct InterfaceSelector<MF, InVector, OutVector, Container, false>
3336  {
3337  typedef void (Container::*function_type)
3338  (const MF &, OutVector &, const InVector &,
3339  const std::pair<unsigned int, unsigned int> &);
3340  };
3341  }
3342 
3343 
3344 
3345  // A implementation class for the worker object that runs the various
3346  // operations we want to perform during the matrix-free loop
3347  template <typename MF, typename InVector, typename OutVector,
3348  typename Container, bool is_constant>
3349  class MFWorker : public MFWorkerInterface
3350  {
3351  public:
3352  // A typedef to make the arguments further down more readable
3353  typedef typename MatrixFreeFunctions::InterfaceSelector
3354  <MF,InVector,OutVector,Container,is_constant>::function_type function_type;
3355 
3356  // constructor, binds all the arguments to this class
3357  MFWorker (const MF &matrix_free,
3358  const InVector &src,
3359  OutVector &dst,
3360  const bool zero_dst_vector_setting,
3361  const Container &container,
3362  function_type cell_function,
3363  function_type face_function,
3364  function_type boundary_function,
3365  const typename MF::DataAccessOnFaces src_vector_face_access =
3366  MF::DataAccessOnFaces::none,
3367  const typename MF::DataAccessOnFaces dst_vector_face_access =
3368  MF::DataAccessOnFaces::none)
3369  :
3370  matrix_free (matrix_free),
3371  container (const_cast<Container &>(container)),
3372  cell_function (cell_function),
3373  face_function (face_function),
3374  boundary_function (boundary_function),
3375  src (src),
3376  dst (dst),
3377  src_data_exchanger (matrix_free, src_vector_face_access,
3378  n_components(src)),
3379  dst_data_exchanger (matrix_free, dst_vector_face_access,
3380  n_components(dst)),
3381  src_and_dst_are_same (PointerComparison::equal(&src, &dst)),
3382  zero_dst_vector_setting(zero_dst_vector_setting && !src_and_dst_are_same)
3383  {}
3384 
3385  // Runs the cell work. If no function is given, nothing is done
3386  virtual void cell(const std::pair<unsigned int,unsigned int> &cell_range) override
3387  {
3388  if (cell_function != nullptr && cell_range.second > cell_range.first)
3389  (container.*cell_function)(matrix_free, this->dst, this->src, cell_range);
3390  }
3391 
3392  // Runs the assembler on interior faces. If no function is given, nothing
3393  // is done
3394  virtual void face(const std::pair<unsigned int,unsigned int> &face_range) override
3395  {
3396  if (face_function != nullptr && face_range.second > face_range.first)
3397  (container.*face_function)(matrix_free, this->dst, this->src, face_range);
3398  }
3399 
3400  // Runs the assembler on boundary faces. If no function is given, nothing
3401  // is done
3402  virtual void boundary(const std::pair<unsigned int,unsigned int> &face_range) override
3403  {
3404  if (boundary_function != nullptr && face_range.second > face_range.first)
3405  (container.*boundary_function)(matrix_free, this->dst, this->src, face_range);
3406  }
3407 
3408  // Starts the communication for the update ghost values operation. We
3409  // cannot call this update if ghost and destination are the same because
3410  // that would introduce spurious entries in the destination (there is also
3411  // the problem that reading from a vector that we also write to is usually
3412  // not intended in case there is overlap, but this is up to the
3413  // application code to decide and we cannot catch this case here).
3414  virtual void vector_update_ghosts_start() override
3415  {
3416  if (!src_and_dst_are_same)
3417  internal::update_ghost_values_start(src, src_data_exchanger);
3418  }
3419 
3420  // Finishes the communication for the update ghost values operation
3421  virtual void vector_update_ghosts_finish() override
3422  {
3423  if (!src_and_dst_are_same)
3424  internal::update_ghost_values_finish(src, src_data_exchanger);
3425  }
3426 
3427  // Starts the communication for the vector compress operation
3428  virtual void vector_compress_start() override
3429  {
3430  internal::compress_start(dst, dst_data_exchanger);
3431  }
3432 
3433  // Finishes the communication for the vector compress operation
3434  virtual void vector_compress_finish() override
3435  {
3436  internal::compress_finish(dst, dst_data_exchanger);
3437  if (!src_and_dst_are_same)
3438  internal::reset_ghost_values(src, src_data_exchanger);
3439  }
3440 
3441  // Zeros the given input vector
3442  virtual void zero_dst_vector_range(const unsigned int range_index) override
3443  {
3444  if (zero_dst_vector_setting)
3445  internal::zero_vector_region(range_index, dst, dst_data_exchanger);
3446  }
3447 
3448  private:
3449  const MF &matrix_free;
3450  Container &container;
3451  function_type cell_function;
3452  function_type face_function;
3453  function_type boundary_function;
3454 
3455  const InVector &src;
3456  OutVector &dst;
3457  VectorDataExchange<MF::dimension,typename MF::value_type> src_data_exchanger;
3458  VectorDataExchange<MF::dimension,typename MF::value_type> dst_data_exchanger;
3459  const bool src_and_dst_are_same;
3460  const bool zero_dst_vector_setting;
3461  };
3462 
3463 
3464 
3469  template <class MF, typename InVector, typename OutVector>
3470  struct MFClassWrapper
3471  {
3472  typedef std::function<void (const MF &, OutVector &, const InVector &,
3473  const std::pair<unsigned int, unsigned int> &)> function_type;
3474 
3475  MFClassWrapper (const function_type cell,
3476  const function_type face,
3477  const function_type boundary)
3478  :
3479  cell (cell),
3480  face (face),
3481  boundary (boundary)
3482  {}
3483 
3484  void cell_integrator (const MF &mf, OutVector &dst, const InVector &src,
3485  const std::pair<unsigned int, unsigned int> &range) const
3486  {
3487  if (cell)
3488  cell(mf, dst, src, range);
3489  }
3490 
3491  void face_integrator (const MF &mf, OutVector &dst, const InVector &src,
3492  const std::pair<unsigned int, unsigned int> &range) const
3493  {
3494  if (face)
3495  face(mf, dst, src, range);
3496  }
3497 
3498  void boundary_integrator (const MF &mf, OutVector &dst, const InVector &src,
3499  const std::pair<unsigned int, unsigned int> &range) const
3500  {
3501  if (boundary)
3502  boundary(mf, dst, src, range);
3503  }
3504 
3505  const function_type cell;
3506  const function_type face;
3507  const function_type boundary;
3508  };
3509 
3510 } // end of namespace internal
3511 
3512 
3513 
3514 template <int dim, typename Number>
3515 template <typename OutVector, typename InVector>
3516 inline
3517 void
3518 MatrixFree<dim, Number>::cell_loop
3519 (const std::function<void (const MatrixFree<dim,Number> &,
3520  OutVector &,
3521  const InVector &,
3522  const std::pair<unsigned int,
3523  unsigned int> &)> &cell_operation,
3524  OutVector &dst,
3525  const InVector &src,
3526  const bool zero_dst_vector) const
3527 {
3528  typedef internal::MFClassWrapper<MatrixFree<dim, Number>, InVector, OutVector> Wrapper;
3529  Wrapper wrap (cell_operation, nullptr, nullptr);
3530  internal::MFWorker<MatrixFree<dim, Number>, InVector, OutVector, Wrapper, true>
3531  worker(*this, src, dst, zero_dst_vector, wrap, &Wrapper::cell_integrator,
3532  &Wrapper::face_integrator, &Wrapper::boundary_integrator);
3533 
3534  task_info.loop (worker);
3535 }
3536 
3537 
3538 
3539 template <int dim, typename Number>
3540 template <typename OutVector, typename InVector>
3541 inline
3542 void
3543 MatrixFree<dim, Number>::loop
3544 (const std::function<void (const MatrixFree<dim,Number> &,
3545  OutVector &,
3546  const InVector &,
3547  const std::pair<unsigned int,
3548  unsigned int> &)> &cell_operation,
3549  const std::function<void (const MatrixFree<dim,Number> &,
3550  OutVector &,
3551  const InVector &,
3552  const std::pair<unsigned int,
3553  unsigned int> &)> &face_operation,
3554  const std::function<void (const MatrixFree<dim,Number> &,
3555  OutVector &,
3556  const InVector &,
3557  const std::pair<unsigned int,
3558  unsigned int> &)> &boundary_operation,
3559  OutVector &dst,
3560  const InVector &src,
3561  const bool zero_dst_vector,
3562  const DataAccessOnFaces dst_vector_face_access,
3563  const DataAccessOnFaces src_vector_face_access) const
3564 {
3565  typedef internal::MFClassWrapper<MatrixFree<dim, Number>, InVector, OutVector> Wrapper;
3566  Wrapper wrap (cell_operation, face_operation, boundary_operation);
3567  internal::MFWorker<MatrixFree<dim, Number>, InVector, OutVector, Wrapper, true>
3568  worker(*this, src, dst, zero_dst_vector, wrap, &Wrapper::cell_integrator,
3569  &Wrapper::face_integrator, &Wrapper::boundary_integrator,
3570  src_vector_face_access, dst_vector_face_access);
3571 
3572  task_info.loop(worker);
3573 }
3574 
3575 
3576 
3577 template <int dim, typename Number>
3578 template <typename CLASS, typename OutVector, typename InVector>
3579 inline
3580 void
3581 MatrixFree<dim,Number>::cell_loop
3582 (void (CLASS::*function_pointer)(const MatrixFree<dim,Number> &,
3583  OutVector &,
3584  const InVector &,
3585  const std::pair<unsigned int, unsigned int> &)const,
3586  const CLASS *owning_class,
3587  OutVector &dst,
3588  const InVector &src,
3589  const bool zero_dst_vector) const
3590 {
3591  internal::MFWorker<MatrixFree<dim, Number>, InVector, OutVector, CLASS, true>
3592  worker(*this, src, dst, zero_dst_vector, *owning_class, function_pointer, nullptr, nullptr);
3593  task_info.loop(worker);
3594 }
3595 
3596 
3597 
3598 template <int dim, typename Number>
3599 template <typename CLASS, typename OutVector, typename InVector>
3600 inline
3601 void
3602 MatrixFree<dim,Number>::loop
3603 (void (CLASS::*cell_operation)(const MatrixFree<dim,Number> &,
3604  OutVector &,
3605  const InVector &,
3606  const std::pair<unsigned int, unsigned int> &)const,
3607  void (CLASS::*face_operation)(const MatrixFree<dim,Number> &,
3608  OutVector &,
3609  const InVector &,
3610  const std::pair<unsigned int, unsigned int> &)const,
3611  void (CLASS::*boundary_operation)(const MatrixFree<dim,Number> &,
3612  OutVector &,
3613  const InVector &,
3614  const std::pair<unsigned int, unsigned int> &)const,
3615  const CLASS *owning_class,
3616  OutVector &dst,
3617  const InVector &src,
3618  const bool zero_dst_vector,
3619  const DataAccessOnFaces dst_vector_face_access,
3620  const DataAccessOnFaces src_vector_face_access) const
3621 {
3622  internal::MFWorker<MatrixFree<dim, Number>, InVector, OutVector, CLASS, true>
3623  worker(*this, src, dst, zero_dst_vector, *owning_class, cell_operation, face_operation,
3624  boundary_operation, src_vector_face_access, dst_vector_face_access);
3625  task_info.loop(worker);
3626 }
3627 
3628 
3629 
3630 template <int dim, typename Number>
3631 template <typename CLASS, typename OutVector, typename InVector>
3632 inline
3633 void
3634 MatrixFree<dim,Number>::cell_loop
3635 (void(CLASS::*function_pointer)(const MatrixFree<dim,Number> &,
3636  OutVector &,
3637  const InVector &,
3638  const std::pair<unsigned int, unsigned int> &),
3639  CLASS *owning_class,
3640  OutVector &dst,
3641  const InVector &src,
3642  const bool zero_dst_vector) const
3643 {
3644  internal::MFWorker<MatrixFree<dim, Number>, InVector, OutVector, CLASS, false>
3645  worker(*this, src, dst, zero_dst_vector, *owning_class, function_pointer, nullptr, nullptr);
3646  task_info.loop(worker);
3647 }
3648 
3649 
3650 
3651 template <int dim, typename Number>
3652 template <typename CLASS, typename OutVector, typename InVector>
3653 inline
3654 void
3655 MatrixFree<dim,Number>::loop
3656 (void(CLASS::*cell_operation)(const MatrixFree<dim,Number> &,
3657  OutVector &,
3658  const InVector &,
3659  const std::pair<unsigned int, unsigned int> &),
3660  void(CLASS::*face_operation)(const MatrixFree<dim,Number> &,
3661  OutVector &,
3662  const InVector &,
3663  const std::pair<unsigned int, unsigned int> &),
3664  void(CLASS::*boundary_operation)(const MatrixFree<dim,Number> &,
3665  OutVector &,
3666  const InVector &,
3667  const std::pair<unsigned int, unsigned int> &),
3668  CLASS *owning_class,
3669  OutVector &dst,
3670  const InVector &src,
3671  const bool zero_dst_vector,
3672  const DataAccessOnFaces dst_vector_face_access,
3673  const DataAccessOnFaces src_vector_face_access) const
3674 {
3675  internal::MFWorker<MatrixFree<dim, Number>, InVector, OutVector, CLASS, false>
3676  worker(*this, src, dst, zero_dst_vector, *owning_class, cell_operation,
3677  face_operation, boundary_operation,
3678  src_vector_face_access, dst_vector_face_access);
3679  task_info.loop(worker);
3680 }
3681 
3682 
3683 #endif // ifndef DOXYGEN
3684 
3685 
3686 
3687 DEAL_II_NAMESPACE_CLOSE
3688 
3689 #endif
MatrixFree::get_ghost_set
const IndexSet & get_ghost_set(const unsigned int dof_handler_index=0) const
update_JxW_values
Transformed quadrature weights.
Definition: fe_update_flags.h:122
MatrixFree::get_face_quadrature
const Quadrature< dim-1 > & get_face_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
MatrixFree::AdditionalData::mapping_update_flags_inner_faces
UpdateFlags mapping_update_flags_inner_faces
Definition: matrix_free.h:327
MatrixFree::n_physical_cells
unsigned int n_physical_cells() const
MatrixFree::get_cell_iterator
DoFHandler< dim >::cell_iterator get_cell_iterator(const unsigned int macro_cell_number, const unsigned int vector_number, const unsigned int fe_component=0) const
LinearAlgebra::distributed::Vector::reinit
void reinit(const size_type size, const bool omit_zeroing_entries=false)
MatrixFree::loop
void loop(const std::function< void(const MatrixFree< dim, Number > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &face_operation, const std::function< void(const MatrixFree< dim, Number > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
MatrixFree::task_info
internal::MatrixFreeFunctions::TaskInfo task_info
Definition: matrix_free.h:1656
MatrixFree::indices_are_initialized
bool indices_are_initialized
Definition: matrix_free.h:1667
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:173
MatrixFree::get_dofs_per_face
unsigned int get_dofs_per_face(const unsigned int fe_component=0, const unsigned int hp_active_fe_index=0) const
MatrixFree::AdditionalData::mapping_update_flags_boundary_faces
UpdateFlags mapping_update_flags_boundary_faces
Definition: matrix_free.h:306
MatrixFree::constraint_pool_begin
const Number * constraint_pool_begin(const unsigned int pool_index) const
MatrixFree::constraint_pool_row_index
std::vector< unsigned int > constraint_pool_row_index
Definition: matrix_free.h:1621
MatrixFree::get_n_q_points_face
unsigned int get_n_q_points_face(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1248
MatrixFree::mapping_info
internal::MatrixFreeFunctions::MappingInfo< dim, Number > mapping_info
Definition: matrix_free.h:1627
Threads::ThreadLocalStorage
A class that provides a separate storage location on each thread that accesses the object...
Definition: thread_local_storage.h:70
MatrixFree::n_ghost_cell_batches
unsigned int n_ghost_cell_batches() const
MatrixFree::dimension
static const unsigned int dimension
Definition: matrix_free.h:117
MatrixFree::scratch_pad_non_threadsafe
std::list< std::pair< bool, AlignedVector< Number > > > scratch_pad_non_threadsafe
Definition: matrix_free.h:1687
MatrixFree::get_face_info
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArray< Number >::n_array_elements > & get_face_info(const unsigned int face_batch_number) const
Utilities::MPI::Partitioner::local_size
unsigned int local_size() const
MatrixFree::dof_handlers
DoFHandlers dof_handlers
Definition: matrix_free.h:1601
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::initialize_dof_handlers
void initialize_dof_handlers(const std::vector< const DoFHandler< dim > *> &dof_handlers, const AdditionalData &additional_data)
MatrixFree::get_size_info
const internal::MatrixFreeFunctions::TaskInfo & get_size_info() const
MatrixFree::n_components
unsigned int n_components() const
MatrixFree::scratch_pad
Threads::ThreadLocalStorage< std::list< std::pair< bool, AlignedVector< VectorizedArray< Number > > > > > scratch_pad
Definition: matrix_free.h:1681
MatrixFree::get_face_category
std::pair< unsigned int, unsigned int > get_face_category(const unsigned int macro_face) const
MatrixFree::get_locally_owned_set
const IndexSet & get_locally_owned_set(const unsigned int dof_handler_index=0) const
MatrixFree::n_active_entries_per_face_batch
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_number) const
types::boundary_id
unsigned int boundary_id
Definition: types.h:110
MatrixFree::constraint_pool_data
std::vector< Number > constraint_pool_data
Definition: matrix_free.h:1615
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1284
internal::MatrixFreeFunctions::DoFInfo::vector_partitioner
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition: dof_info.h:366
internal::MatrixFreeFunctions::DoFInfo::vector_zero_range_list_index
std::vector< unsigned int > vector_zero_range_list_index
Definition: dof_info.h:491
MatrixFree::reinit
void reinit(const Mapping< dim > &mapping, const DoFHandlerType &dof_handler, const ConstraintMatrix &constraint, const QuadratureType &quad, const AdditionalData additional_data=AdditionalData())
MatrixFree::copy_from
void copy_from(const MatrixFree< dim, Number > &matrix_free_base)
MatrixFree::get_faces_by_cells_boundary_id
std::array< types::boundary_id, VectorizedArray< Number >::n_array_elements > get_faces_by_cells_boundary_id(const unsigned int macro_cell, const unsigned int face_number) const
MatrixFree::cell_level_index
std::vector< std::pair< unsigned int, unsigned int > > cell_level_index
Definition: matrix_free.h:1640
Utilities::MPI::Partitioner::ghost_indices_within_larger_ghost_set
const std::vector< std::pair< unsigned int, unsigned int > > & ghost_indices_within_larger_ghost_set() const
MatrixFree::make_connectivity_graph_faces
void make_connectivity_graph_faces(DynamicSparsityPattern &connectivity)
StandardExceptions::ExcNotInitialized
static ::ExceptionBase & ExcNotInitialized()
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1221
MatrixFree::get_hp_cell_iterator
hp::DoFHandler< dim >::active_cell_iterator get_hp_cell_iterator(const unsigned int macro_cell_number, const unsigned int vector_number, const unsigned int dof_handler_index=0) const
StandardExceptions::ExcIndexRange
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
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 > &ghost_array, const ArrayView< Number > &temporary_storage, std::vector< MPI_Request > &requests) 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::n_boundary_face_batches
unsigned int n_boundary_face_batches() const
MatrixFree::n_inner_face_batches
unsigned int n_inner_face_batches() const
hp::DoFHandler::active_cell_iterator
ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:186
MatrixFree::mapping_is_initialized
bool mapping_is_initialized
Definition: matrix_free.h:1672
MatrixFree::acquire_scratch_data
AlignedVector< VectorizedArray< Number > > * acquire_scratch_data() const
Utilities::MPI::Partitioner::export_to_ghosted_array_finish
void export_to_ghosted_array_finish(const ArrayView< Number > &ghost_array, std::vector< MPI_Request > &requests) const
MatrixFree::internal_reinit
void internal_reinit(const Mapping< dim > &mapping, const std::vector< const DoFHandler< dim > *> &dof_handler, const std::vector< const ConstraintMatrix *> &constraint, const std::vector< IndexSet > &locally_owned_set, const std::vector< hp::QCollection< 1 > > &quad, const AdditionalData &additional_data)
MatrixFree::initialize_dof_vector
void initialize_dof_vector(VectorType &vec, const unsigned int dof_handler_index=0) const
MatrixFree::n_components_filled
unsigned int n_components_filled(const unsigned int cell_batch_number) const
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
MatrixFree::clear
void clear()
MatrixFree::renumber_dofs
void renumber_dofs(std::vector< types::global_dof_index > &renumbering, const unsigned int dof_handler_index=0)
update_default
No update.
Definition: fe_update_flags.h:65
MatrixFree::get_vector_partitioner
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner(const unsigned int dof_handler_index=0) const
Utilities::MPI::Partitioner::n_ghost_indices
unsigned int n_ghost_indices() const
LinearAlgebra::distributed::Vector::size
virtual size_type size() const override
MatrixFree::at_irregular_cell
bool at_irregular_cell(const unsigned int macro_cell_number) const
MatrixFree::value_type
Number value_type
Definition: matrix_free.h:112
MatrixFree::shape_info
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< VectorizedArray< Number > > > shape_info
Definition: matrix_free.h:1632
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1142
MatrixFree::is_supported
static bool is_supported(const FiniteElement< dim, spacedim > &fe)
UpdateFlags
UpdateFlags
Definition: fe_update_flags.h:62
MatrixFree::print_memory_consumption
void print_memory_consumption(StreamType &out) const
numbers::invalid_boundary_id
const types::boundary_id invalid_boundary_id
Definition: types.h:204
MatrixFree::constraint_pool_end
const Number * constraint_pool_end(const unsigned int pool_index) const
MatrixFree::get_shape_info
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArray< Number > > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
Utilities::MPI::Partitioner::export_to_ghosted_array_start
void export_to_ghosted_array_start(const unsigned int communication_channel, const ArrayView< const Number > &locally_owned_array, const ArrayView< Number > &temporary_storage, const ArrayView< Number > &ghost_array, std::vector< MPI_Request > &requests) const
MatrixFree::n_constraint_pool_entries
unsigned int n_constraint_pool_entries() const
MatrixFree::AdditionalData::AdditionalData
AdditionalData(const TasksParallelScheme tasks_parallel_scheme=partition_partition, const unsigned int tasks_block_size=0, const UpdateFlags mapping_update_flags=update_gradients|update_JxW_values, const UpdateFlags mapping_update_flags_boundary_faces=update_default, const UpdateFlags mapping_update_flags_inner_faces=update_default, const UpdateFlags mapping_update_flags_faces_by_cells=update_default, const unsigned int level_mg_handler=numbers::invalid_unsigned_int, const bool store_plain_indices=true, const bool initialize_indices=true, const bool initialize_mapping=true, const bool overlap_communication_computation=true, const bool hold_all_faces_to_owned_cells=false, const bool cell_vectorization_categories_strict=false)
Definition: matrix_free.h:196
MatrixFree::release_scratch_data_non_threadsafe
void release_scratch_data_non_threadsafe(const AlignedVector< Number > *memory) const
MatrixFree::n_macro_cells
unsigned int n_macro_cells() const
LinearAlgebra::distributed::Vector::get_partitioner
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_partitioner() const
MatrixFree::AdditionalData::tasks_parallel_scheme
TasksParallelScheme tasks_parallel_scheme
Definition: matrix_free.h:261
MatrixFree::acquire_scratch_data_non_threadsafe
AlignedVector< Number > * acquire_scratch_data_non_threadsafe() const
MatrixFree::get_boundary_id
types::boundary_id get_boundary_id(const unsigned int macro_face) const
MatrixFree::mapping_initialized
bool mapping_initialized() const
hp
Definition: hp.h:102
MatrixFree::n_ghost_inner_face_batches
unsigned int n_ghost_inner_face_batches() const
LinearAlgebra::distributed::Vector::compress_start
void compress_start(const unsigned int communication_channel=0, ::VectorOperation::values operation=VectorOperation::add)
MatrixFree::AdditionalData::tasks_block_size
unsigned int tasks_block_size
Definition: matrix_free.h:272
MatrixFree::release_scratch_data
void release_scratch_data(const AlignedVector< VectorizedArray< Number > > *memory) const
LinearAlgebra::distributed::Vector::compress_finish
void compress_finish(::VectorOperation::values operation)
MatrixFree::~MatrixFree
~MatrixFree()=default
MatrixFree::cell_level_index_end_local
unsigned int cell_level_index_end_local
Definition: matrix_free.h:1649
MatrixFree::face_info
internal::MatrixFreeFunctions::FaceInfo< VectorizedArray< Number >::n_array_elements > face_info
Definition: matrix_free.h:1662
MatrixFree::AdditionalData::level_mg_handler
unsigned int level_mg_handler
Definition: matrix_free.h:365
MatrixFree::get_dof_handler
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
update_gradients
Shape function gradients.
Definition: fe_update_flags.h:80
Utilities::MPI::Partitioner::import_from_ghosted_array_finish
void import_from_ghosted_array_finish(const VectorOperation::values vector_operation, const ArrayView< const Number > &temporary_storage, const ArrayView< Number > &locally_owned_storage, const ArrayView< Number > &ghost_array, std::vector< MPI_Request > &requests) const
MatrixFree::get_n_q_points
unsigned int get_n_q_points(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
AlignedVector::begin
iterator begin()
MatrixFree::n_active_entries_per_cell_batch
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_number) const
LinearAlgebra::distributed::Vector::partitioners_are_compatible
bool partitioners_are_compatible(const Utilities::MPI::Partitioner &part) const
MatrixFree::n_base_elements
unsigned int n_base_elements(const unsigned int dof_handler_index) const
MatrixFree::AdditionalData::mapping_update_flags
UpdateFlags mapping_update_flags
Definition: matrix_free.h:285
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
MatrixFree::memory_consumption
std::size_t memory_consumption() const
internal::MatrixFreeFunctions::DoFInfo::vector_zero_range_list
std::vector< unsigned int > vector_zero_range_list
Definition: dof_info.h:496
Table
Definition: table.h:34
MatrixFree::initialize_indices
void initialize_indices(const std::vector< const ConstraintMatrix *> &constraint, const std::vector< IndexSet > &locally_owned_set, const AdditionalData &additional_data)
MatrixFree::indices_initialized
bool indices_initialized() const
MatrixFree::get_cell_category
unsigned int get_cell_category(const unsigned int macro_cell) const
MatrixFree::get_dof_info
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
MatrixFree::get_quadrature
const Quadrature< dim > & get_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
DoFHandler::n_components
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
MatrixFree::dof_info
std::vector< internal::MatrixFreeFunctions::DoFInfo > dof_info
Definition: matrix_free.h:1607
internal::MatrixFreeFunctions::DoFInfo::chunk_size_zero_vector
static const unsigned int chunk_size_zero_vector
Definition: dof_info.h:74
DoFTools::n_components
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
Utilities::MPI::Partitioner::n_import_indices
unsigned int n_import_indices() const
MatrixFree::n_cell_batches
unsigned int n_cell_batches() const
MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells
UpdateFlags mapping_update_flags_faces_by_cells
Definition: matrix_free.h:355
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
MatrixFree::get_task_info
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
MatrixFree::AdditionalData::cell_vectorization_category
std::vector< unsigned int > cell_vectorization_category
Definition: matrix_free.h:434
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
MatrixFree::cell_loop
void cell_loop(const std::function< void(const MatrixFree< dim, Number > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
MatrixFree::print
void print(std::ostream &out) const
LinearAlgebra::distributed::Vector::update_ghost_values_start
void update_ghost_values_start(const unsigned int communication_channel=0) const
MatrixFree::get_constrained_dofs
const std::vector< unsigned int > & get_constrained_dofs(const unsigned int dof_handler_index=0) const
Generated by   doxygen 1.8.14