Reference documentation for deal.II version 9.4.1
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
matrix_free.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2011 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
17#ifndef dealii_matrix_free_h
18#define dealii_matrix_free_h
19
20#include <deal.II/base/config.h>
21
22#include <deal.II/base/aligned_vector.h>
23#include <deal.II/base/exceptions.h>
24#include <deal.II/base/quadrature.h>
25#include <deal.II/base/template_constraints.h>
26#include <deal.II/base/thread_local_storage.h>
27#include <deal.II/base/vectorization.h>
28
29#include <deal.II/dofs/dof_handler.h>
30
31#include <deal.II/fe/fe.h>
32#include <deal.II/fe/mapping.h>
33#include <deal.II/fe/mapping_q1.h>
34
35#include <deal.II/hp/dof_handler.h>
36#include <deal.II/hp/mapping_collection.h>
37#include <deal.II/hp/q_collection.h>
38
39#include <deal.II/lac/affine_constraints.h>
40#include <deal.II/lac/block_vector_base.h>
41#include <deal.II/lac/la_parallel_vector.h>
42#include <deal.II/lac/vector_operation.h>
43
44#include <deal.II/matrix_free/dof_info.h>
45#include <deal.II/matrix_free/mapping_info.h>
46#include <deal.II/matrix_free/shape_info.h>
47#include <deal.II/matrix_free/task_info.h>
48#include <deal.II/matrix_free/type_traits.h>
49
50#include <cstdlib>
51#include <limits>
52#include <list>
53#include <memory>
54
55
56DEAL_II_NAMESPACE_OPEN
57
58
59
111template <int dim,
112 typename Number = double,
113 typename VectorizedArrayType = VectorizedArray<Number>>
114class MatrixFree : public Subscriptor
115{
116 static_assert(
117 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
118 "Type of Number and of VectorizedArrayType do not match.");
119
120public:
125 using value_type = Number;
126 using vectorized_value_type = VectorizedArrayType;
127
131 static constexpr unsigned int dimension = dim;
132
179 struct AdditionalData
180 {
184 using MatrixFreeType = MatrixFree<dim, Number, VectorizedArrayType>;
185
191 enum TasksParallelScheme
192 {
196 none = internal::MatrixFreeFunctions::TaskInfo::none,
200 partition_partition =
201 internal::MatrixFreeFunctions::TaskInfo::partition_partition,
205 partition_color =
206 internal::MatrixFreeFunctions::TaskInfo::partition_color,
211 color = internal::MatrixFreeFunctions::TaskInfo::color
212 };
213
217 AdditionalData(
218 const TasksParallelScheme tasks_parallel_scheme = partition_partition,
219 const unsigned int tasks_block_size = 0,
220 const UpdateFlags mapping_update_flags = update_gradients |
221 update_JxW_values,
222 const UpdateFlags mapping_update_flags_boundary_faces = update_default,
223 const UpdateFlags mapping_update_flags_inner_faces = update_default,
224 const UpdateFlags mapping_update_flags_faces_by_cells = update_default,
225 const unsigned int mg_level = numbers::invalid_unsigned_int,
226 const bool store_plain_indices = true,
227 const bool initialize_indices = true,
228 const bool initialize_mapping = true,
229 const bool overlap_communication_computation = true,
230 const bool hold_all_faces_to_owned_cells = false,
231 const bool cell_vectorization_categories_strict = false,
232 const bool allow_ghosted_vectors_in_loops = true)
233 : tasks_parallel_scheme(tasks_parallel_scheme)
234 , tasks_block_size(tasks_block_size)
235 , mapping_update_flags(mapping_update_flags)
236 , mapping_update_flags_boundary_faces(mapping_update_flags_boundary_faces)
237 , mapping_update_flags_inner_faces(mapping_update_flags_inner_faces)
238 , mapping_update_flags_faces_by_cells(mapping_update_flags_faces_by_cells)
239 , mg_level(mg_level)
240 , store_plain_indices(store_plain_indices)
241 , initialize_indices(initialize_indices)
242 , initialize_mapping(initialize_mapping)
243 , overlap_communication_computation(overlap_communication_computation)
244 , hold_all_faces_to_owned_cells(hold_all_faces_to_owned_cells)
245 , cell_vectorization_categories_strict(
246 cell_vectorization_categories_strict)
247 , allow_ghosted_vectors_in_loops(allow_ghosted_vectors_in_loops)
248 , communicator_sm(MPI_COMM_SELF)
249 {}
250
254 AdditionalData(const AdditionalData &other)
255 : tasks_parallel_scheme(other.tasks_parallel_scheme)
256 , tasks_block_size(other.tasks_block_size)
257 , mapping_update_flags(other.mapping_update_flags)
258 , mapping_update_flags_boundary_faces(
259 other.mapping_update_flags_boundary_faces)
260 , mapping_update_flags_inner_faces(other.mapping_update_flags_inner_faces)
261 , mapping_update_flags_faces_by_cells(
262 other.mapping_update_flags_faces_by_cells)
263 , mg_level(other.mg_level)
264 , store_plain_indices(other.store_plain_indices)
265 , initialize_indices(other.initialize_indices)
266 , initialize_mapping(other.initialize_mapping)
267 , overlap_communication_computation(
268 other.overlap_communication_computation)
269 , hold_all_faces_to_owned_cells(other.hold_all_faces_to_owned_cells)
270 , cell_vectorization_category(other.cell_vectorization_category)
271 , cell_vectorization_categories_strict(
272 other.cell_vectorization_categories_strict)
273 , allow_ghosted_vectors_in_loops(other.allow_ghosted_vectors_in_loops)
274 , communicator_sm(other.communicator_sm)
275 {}
276
280 AdditionalData &
281 operator=(const AdditionalData &other)
282 {
283 tasks_parallel_scheme = other.tasks_parallel_scheme;
284 tasks_block_size = other.tasks_block_size;
285 mapping_update_flags = other.mapping_update_flags;
286 mapping_update_flags_boundary_faces =
287 other.mapping_update_flags_boundary_faces;
288 mapping_update_flags_inner_faces = other.mapping_update_flags_inner_faces;
289 mapping_update_flags_faces_by_cells =
290 other.mapping_update_flags_faces_by_cells;
291 mg_level = other.mg_level;
292 store_plain_indices = other.store_plain_indices;
293 initialize_indices = other.initialize_indices;
294 initialize_mapping = other.initialize_mapping;
295 overlap_communication_computation =
296 other.overlap_communication_computation;
297 hold_all_faces_to_owned_cells = other.hold_all_faces_to_owned_cells;
298 cell_vectorization_category = other.cell_vectorization_category;
299 cell_vectorization_categories_strict =
300 other.cell_vectorization_categories_strict;
301 allow_ghosted_vectors_in_loops = other.allow_ghosted_vectors_in_loops;
302 communicator_sm = other.communicator_sm;
303
304 return *this;
305 }
306
343 TasksParallelScheme tasks_parallel_scheme;
344
354 unsigned int tasks_block_size;
355
367 UpdateFlags mapping_update_flags;
368
387 UpdateFlags mapping_update_flags_boundary_faces;
388
407 UpdateFlags mapping_update_flags_inner_faces;
408
435 UpdateFlags mapping_update_flags_faces_by_cells;
436
445 unsigned int mg_level;
446
453 bool store_plain_indices;
454
464 bool initialize_indices;
465
473 bool initialize_mapping;
474
487 bool overlap_communication_computation;
488
496 bool hold_all_faces_to_owned_cells;
497
514 std::vector<unsigned int> cell_vectorization_category;
515
525 bool cell_vectorization_categories_strict;
526
537 bool allow_ghosted_vectors_in_loops;
538
542 MPI_Comm communicator_sm;
543 };
544
552 MatrixFree();
553
557 MatrixFree(const MatrixFree<dim, Number, VectorizedArrayType> &other);
558
562 ~MatrixFree() override = default;
563
576 template <typename QuadratureType, typename number2, typename MappingType>
577 void
578 reinit(const MappingType & mapping,
579 const DoFHandler<dim> & dof_handler,
580 const AffineConstraints<number2> &constraint,
581 const QuadratureType & quad,
582 const AdditionalData & additional_data = AdditionalData());
583
590 template <typename QuadratureType, typename number2>
591 DEAL_II_DEPRECATED void
592 reinit(const DoFHandler<dim> & dof_handler,
593 const AffineConstraints<number2> &constraint,
594 const QuadratureType & quad,
595 const AdditionalData & additional_data = AdditionalData());
596
618 template <typename QuadratureType, typename number2, typename MappingType>
619 void
620 reinit(const MappingType & mapping,
621 const std::vector<const DoFHandler<dim> *> & dof_handler,
622 const std::vector<const AffineConstraints<number2> *> &constraint,
623 const std::vector<QuadratureType> & quad,
624 const AdditionalData &additional_data = AdditionalData());
625
631 template <typename QuadratureType,
632 typename number2,
633 typename DoFHandlerType,
634 typename MappingType>
635 DEAL_II_DEPRECATED void
636 reinit(const MappingType & mapping,
637 const std::vector<const DoFHandlerType *> & dof_handler,
638 const std::vector<const AffineConstraints<number2> *> &constraint,
639 const std::vector<QuadratureType> & quad,
640 const AdditionalData &additional_data = AdditionalData());
641
648 template <typename QuadratureType, typename number2>
649 DEAL_II_DEPRECATED void
650 reinit(const std::vector<const DoFHandler<dim> *> & dof_handler,
651 const std::vector<const AffineConstraints<number2> *> &constraint,
652 const std::vector<QuadratureType> & quad,
653 const AdditionalData &additional_data = AdditionalData());
654
660 template <typename QuadratureType, typename number2, typename DoFHandlerType>
661 DEAL_II_DEPRECATED void
662 reinit(const std::vector<const DoFHandlerType *> & dof_handler,
663 const std::vector<const AffineConstraints<number2> *> &constraint,
664 const std::vector<QuadratureType> & quad,
665 const AdditionalData &additional_data = AdditionalData());
666
674 template <typename QuadratureType, typename number2, typename MappingType>
675 void
676 reinit(const MappingType & mapping,
677 const std::vector<const DoFHandler<dim> *> & dof_handler,
678 const std::vector<const AffineConstraints<number2> *> &constraint,
679 const QuadratureType & quad,
680 const AdditionalData &additional_data = AdditionalData());
681
687 template <typename QuadratureType,
688 typename number2,
689 typename DoFHandlerType,
690 typename MappingType>
691 DEAL_II_DEPRECATED void
692 reinit(const MappingType & mapping,
693 const std::vector<const DoFHandlerType *> & dof_handler,
694 const std::vector<const AffineConstraints<number2> *> &constraint,
695 const QuadratureType & quad,
696 const AdditionalData &additional_data = AdditionalData());
697
704 template <typename QuadratureType, typename number2>
705 DEAL_II_DEPRECATED void
706 reinit(const std::vector<const DoFHandler<dim> *> & dof_handler,
707 const std::vector<const AffineConstraints<number2> *> &constraint,
708 const QuadratureType & quad,
709 const AdditionalData &additional_data = AdditionalData());
710
716 template <typename QuadratureType, typename number2, typename DoFHandlerType>
717 DEAL_II_DEPRECATED void
718 reinit(const std::vector<const DoFHandlerType *> & dof_handler,
719 const std::vector<const AffineConstraints<number2> *> &constraint,
720 const QuadratureType & quad,
721 const AdditionalData &additional_data = AdditionalData());
722
728 void
729 copy_from(
730 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free_base);
731
741 void
742 update_mapping(const Mapping<dim> &mapping);
743
747 void
748 update_mapping(const std::shared_ptr<hp::MappingCollection<dim>> &mapping);
749
754 void
755 clear();
756
758
770 enum class DataAccessOnFaces
771 {
778 none,
779
790 values,
791
800 values_all_faces,
801
812 gradients,
813
822 gradients_all_faces,
823
830 unspecified
831 };
832
877 template <typename OutVector, typename InVector>
878 void
879 cell_loop(const std::function<void(
880 const MatrixFree<dim, Number, VectorizedArrayType> &,
881 OutVector &,
882 const InVector &,
883 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
884 OutVector & dst,
885 const InVector & src,
886 const bool zero_dst_vector = false) const;
887
934 template <typename CLASS, typename OutVector, typename InVector>
935 void
936 cell_loop(void (CLASS::*cell_operation)(
937 const MatrixFree &,
938 OutVector &,
939 const InVector &,
940 const std::pair<unsigned int, unsigned int> &) const,
941 const CLASS * owning_class,
942 OutVector & dst,
943 const InVector &src,
944 const bool zero_dst_vector = false) const;
945
949 template <typename CLASS, typename OutVector, typename InVector>
950 void
951 cell_loop(void (CLASS::*cell_operation)(
952 const MatrixFree &,
953 OutVector &,
954 const InVector &,
955 const std::pair<unsigned int, unsigned int> &),
956 CLASS * owning_class,
957 OutVector & dst,
958 const InVector &src,
959 const bool zero_dst_vector = false) const;
960
1045 template <typename CLASS, typename OutVector, typename InVector>
1046 void
1047 cell_loop(void (CLASS::*cell_operation)(
1048 const MatrixFree &,
1049 OutVector &,
1050 const InVector &,
1051 const std::pair<unsigned int, unsigned int> &) const,
1052 const CLASS * owning_class,
1053 OutVector & dst,
1054 const InVector &src,
1055 const std::function<void(const unsigned int, const unsigned int)>
1056 &operation_before_loop,
1057 const std::function<void(const unsigned int, const unsigned int)>
1058 & operation_after_loop,
1059 const unsigned int dof_handler_index_pre_post = 0) const;
1060
1064 template <typename CLASS, typename OutVector, typename InVector>
1065 void
1066 cell_loop(void (CLASS::*cell_operation)(
1067 const MatrixFree &,
1068 OutVector &,
1069 const InVector &,
1070 const std::pair<unsigned int, unsigned int> &),
1071 CLASS * owning_class,
1072 OutVector & dst,
1073 const InVector &src,
1074 const std::function<void(const unsigned int, const unsigned int)>
1075 &operation_before_loop,
1076 const std::function<void(const unsigned int, const unsigned int)>
1077 & operation_after_loop,
1078 const unsigned int dof_handler_index_pre_post = 0) const;
1079
1084 template <typename OutVector, typename InVector>
1085 void
1086 cell_loop(const std::function<void(
1087 const MatrixFree<dim, Number, VectorizedArrayType> &,
1088 OutVector &,
1089 const InVector &,
1090 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1091 OutVector & dst,
1092 const InVector & src,
1093 const std::function<void(const unsigned int, const unsigned int)>
1094 &operation_before_loop,
1095 const std::function<void(const unsigned int, const unsigned int)>
1096 & operation_after_loop,
1097 const unsigned int dof_handler_index_pre_post = 0) const;
1098
1174 template <typename OutVector, typename InVector>
1175 void
1176 loop(const std::function<
1177 void(const MatrixFree<dim, Number, VectorizedArrayType> &,
1178 OutVector &,
1179 const InVector &,
1180 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1181 const std::function<
1182 void(const MatrixFree<dim, Number, VectorizedArrayType> &,
1183 OutVector &,
1184 const InVector &,
1185 const std::pair<unsigned int, unsigned int> &)> &face_operation,
1186 const std::function<void(
1187 const MatrixFree<dim, Number, VectorizedArrayType> &,
1188 OutVector &,
1189 const InVector &,
1190 const std::pair<unsigned int, unsigned int> &)> &boundary_operation,
1191 OutVector & dst,
1192 const InVector & src,
1193 const bool zero_dst_vector = false,
1194 const DataAccessOnFaces dst_vector_face_access =
1195 DataAccessOnFaces::unspecified,
1196 const DataAccessOnFaces src_vector_face_access =
1197 DataAccessOnFaces::unspecified) const;
1198
1284 template <typename CLASS, typename OutVector, typename InVector>
1285 void
1286 loop(
1287 void (CLASS::*cell_operation)(const MatrixFree &,
1288 OutVector &,
1289 const InVector &,
1290 const std::pair<unsigned int, unsigned int> &)
1291 const,
1292 void (CLASS::*face_operation)(const MatrixFree &,
1293 OutVector &,
1294 const InVector &,
1295 const std::pair<unsigned int, unsigned int> &)
1296 const,
1297 void (CLASS::*boundary_operation)(
1298 const MatrixFree &,
1299 OutVector &,
1300 const InVector &,
1301 const std::pair<unsigned int, unsigned int> &) const,
1302 const CLASS * owning_class,
1303 OutVector & dst,
1304 const InVector & src,
1305 const bool zero_dst_vector = false,
1306 const DataAccessOnFaces dst_vector_face_access =
1307 DataAccessOnFaces::unspecified,
1308 const DataAccessOnFaces src_vector_face_access =
1309 DataAccessOnFaces::unspecified) const;
1310
1314 template <typename CLASS, typename OutVector, typename InVector>
1315 void
1316 loop(void (CLASS::*cell_operation)(
1317 const MatrixFree &,
1318 OutVector &,
1319 const InVector &,
1320 const std::pair<unsigned int, unsigned int> &),
1321 void (CLASS::*face_operation)(
1322 const MatrixFree &,
1323 OutVector &,
1324 const InVector &,
1325 const std::pair<unsigned int, unsigned int> &),
1326 void (CLASS::*boundary_operation)(
1327 const MatrixFree &,
1328 OutVector &,
1329 const InVector &,
1330 const std::pair<unsigned int, unsigned int> &),
1331 CLASS * owning_class,
1332 OutVector & dst,
1333 const InVector & src,
1334 const bool zero_dst_vector = false,
1335 const DataAccessOnFaces dst_vector_face_access =
1336 DataAccessOnFaces::unspecified,
1337 const DataAccessOnFaces src_vector_face_access =
1338 DataAccessOnFaces::unspecified) const;
1339
1404 template <typename CLASS, typename OutVector, typename InVector>
1405 void
1406 loop_cell_centric(void (CLASS::*cell_operation)(
1407 const MatrixFree &,
1408 OutVector &,
1409 const InVector &,
1410 const std::pair<unsigned int, unsigned int> &) const,
1411 const CLASS * owning_class,
1412 OutVector & dst,
1413 const InVector & src,
1414 const bool zero_dst_vector = false,
1415 const DataAccessOnFaces src_vector_face_access =
1416 DataAccessOnFaces::unspecified) const;
1417
1421 template <typename CLASS, typename OutVector, typename InVector>
1422 void
1423 loop_cell_centric(void (CLASS::*cell_operation)(
1424 const MatrixFree &,
1425 OutVector &,
1426 const InVector &,
1427 const std::pair<unsigned int, unsigned int> &),
1428 CLASS * owning_class,
1429 OutVector & dst,
1430 const InVector & src,
1431 const bool zero_dst_vector = false,
1432 const DataAccessOnFaces src_vector_face_access =
1433 DataAccessOnFaces::unspecified) const;
1434
1438 template <typename OutVector, typename InVector>
1439 void
1440 loop_cell_centric(
1441 const std::function<void(const MatrixFree &,
1442 OutVector &,
1443 const InVector &,
1444 const std::pair<unsigned int, unsigned int> &)>
1445 & cell_operation,
1446 OutVector & dst,
1447 const InVector & src,
1448 const bool zero_dst_vector = false,
1449 const DataAccessOnFaces src_vector_face_access =
1450 DataAccessOnFaces::unspecified) const;
1451
1459 std::pair<unsigned int, unsigned int>
1460 create_cell_subrange_hp(const std::pair<unsigned int, unsigned int> &range,
1461 const unsigned int fe_degree,
1462 const unsigned int dof_handler_index = 0) const;
1463
1470 std::pair<unsigned int, unsigned int>
1471 create_cell_subrange_hp_by_index(
1472 const std::pair<unsigned int, unsigned int> &range,
1473 const unsigned int fe_index,
1474 const unsigned int dof_handler_index = 0) const;
1475
1479 unsigned int
1480 n_active_fe_indices() const;
1481
1485 unsigned int
1486 get_cell_active_fe_index(
1487 const std::pair<unsigned int, unsigned int> range) const;
1488
1492 unsigned int
1493 get_face_active_fe_index(const std::pair<unsigned int, unsigned int> range,
1494 const bool is_interior_face = true) const;
1495
1497
1507 template <typename T>
1508 void
1509 initialize_cell_data_vector(AlignedVector<T> &vec) const;
1510
1516 template <typename T>
1517 void
1518 initialize_face_data_vector(AlignedVector<T> &vec) const;
1519
1533 template <typename VectorType>
1534 void
1535 initialize_dof_vector(VectorType & vec,
1536 const unsigned int dof_handler_index = 0) const;
1537
1555 template <typename Number2>
1556 void
1557 initialize_dof_vector(LinearAlgebra::distributed::Vector<Number2> &vec,
1558 const unsigned int dof_handler_index = 0) const;
1559
1570 const std::shared_ptr<const Utilities::MPI::Partitioner> &
1571 get_vector_partitioner(const unsigned int dof_handler_index = 0) const;
1572
1576 const IndexSet &
1577 get_locally_owned_set(const unsigned int dof_handler_index = 0) const;
1578
1582 const IndexSet &
1583 get_ghost_set(const unsigned int dof_handler_index = 0) const;
1584
1594 const std::vector<unsigned int> &
1595 get_constrained_dofs(const unsigned int dof_handler_index = 0) const;
1596
1607 void
1608 renumber_dofs(std::vector<types::global_dof_index> &renumbering,
1609 const unsigned int dof_handler_index = 0);
1610
1612
1620 template <int spacedim>
1621 static bool
1622 is_supported(const FiniteElement<dim, spacedim> &fe);
1623
1627 unsigned int
1628 n_components() const;
1629
1634 unsigned int
1635 n_base_elements(const unsigned int dof_handler_index) const;
1636
1644 unsigned int
1645 n_physical_cells() const;
1646
1650 DEAL_II_DEPRECATED unsigned int
1651 n_macro_cells() const;
1652
1662 unsigned int
1663 n_cell_batches() const;
1664
1671 unsigned int
1672 n_ghost_cell_batches() const;
1673
1683 unsigned int
1684 n_inner_face_batches() const;
1685
1696 unsigned int
1697 n_boundary_face_batches() const;
1698
1703 unsigned int
1704 n_ghost_inner_face_batches() const;
1705
1712 types::boundary_id
1713 get_boundary_id(const unsigned int face_batch_index) const;
1714
1719 std::array<types::boundary_id, VectorizedArrayType::size()>
1720 get_faces_by_cells_boundary_id(const unsigned int cell_batch_index,
1721 const unsigned int face_number) const;
1722
1727 const DoFHandler<dim> &
1728 get_dof_handler(const unsigned int dof_handler_index = 0) const;
1729
1740 template <typename DoFHandlerType>
1741 DEAL_II_DEPRECATED const DoFHandlerType &
1742 get_dof_handler(const unsigned int dof_handler_index = 0) const;
1743
1756 typename DoFHandler<dim>::cell_iterator
1757 get_cell_iterator(const unsigned int cell_batch_index,
1758 const unsigned int lane_index,
1759 const unsigned int dof_handler_index = 0) const;
1760
1766 std::pair<int, int>
1767 get_cell_level_and_index(const unsigned int cell_batch_index,
1768 const unsigned int lane_index) const;
1769
1782 std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>
1783 get_face_iterator(const unsigned int face_batch_index,
1784 const unsigned int lane_index,
1785 const bool interior = true,
1786 const unsigned int fe_component = 0) const;
1787
1793 DEAL_II_DEPRECATED typename DoFHandler<dim>::active_cell_iterator
1794 get_hp_cell_iterator(const unsigned int cell_batch_index,
1795 const unsigned int lane_index,
1796 const unsigned int dof_handler_index = 0) const;
1797
1810 bool
1811 at_irregular_cell(const unsigned int cell_batch_index) const;
1812
1816 DEAL_II_DEPRECATED unsigned int
1817 n_components_filled(const unsigned int cell_batch_number) const;
1818
1828 unsigned int
1829 n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const;
1830
1840 unsigned int
1841 n_active_entries_per_face_batch(const unsigned int face_batch_index) const;
1842
1846 unsigned int
1847 get_dofs_per_cell(const unsigned int dof_handler_index = 0,
1848 const unsigned int hp_active_fe_index = 0) const;
1849
1853 unsigned int
1854 get_n_q_points(const unsigned int quad_index = 0,
1855 const unsigned int hp_active_fe_index = 0) const;
1856
1861 unsigned int
1862 get_dofs_per_face(const unsigned int dof_handler_index = 0,
1863 const unsigned int hp_active_fe_index = 0) const;
1864
1869 unsigned int
1870 get_n_q_points_face(const unsigned int quad_index = 0,
1871 const unsigned int hp_active_fe_index = 0) const;
1872
1876 const Quadrature<dim> &
1877 get_quadrature(const unsigned int quad_index = 0,
1878 const unsigned int hp_active_fe_index = 0) const;
1879
1883 const Quadrature<dim - 1> &
1884 get_face_quadrature(const unsigned int quad_index = 0,
1885 const unsigned int hp_active_fe_index = 0) const;
1886
1899 unsigned int
1900 get_cell_range_category(
1901 const std::pair<unsigned int, unsigned int> cell_batch_range) const;
1902
1907 std::pair<unsigned int, unsigned int>
1908 get_face_range_category(
1909 const std::pair<unsigned int, unsigned int> face_batch_range) const;
1910
1922 unsigned int
1923 get_cell_category(const unsigned int cell_batch_index) const;
1924
1929 std::pair<unsigned int, unsigned int>
1930 get_face_category(const unsigned int face_batch_index) const;
1931
1935 bool
1936 indices_initialized() const;
1937
1942 bool
1943 mapping_initialized() const;
1944
1949 unsigned int
1950 get_mg_level() const;
1951
1956 std::size_t
1957 memory_consumption() const;
1958
1963 template <typename StreamType>
1964 void
1965 print_memory_consumption(StreamType &out) const;
1966
1971 void
1972 print(std::ostream &out) const;
1973
1975
1985 const internal::MatrixFreeFunctions::TaskInfo &
1986 get_task_info() const;
1987
1988 /*
1989 * Return geometry-dependent information on the cells.
1990 */
1991 const internal::MatrixFreeFunctions::
1992 MappingInfo<dim, Number, VectorizedArrayType> &
1993 get_mapping_info() const;
1994
1998 const internal::MatrixFreeFunctions::DoFInfo &
1999 get_dof_info(const unsigned int dof_handler_index_component = 0) const;
2000
2004 unsigned int
2005 n_constraint_pool_entries() const;
2006
2011 const Number *
2012 constraint_pool_begin(const unsigned int pool_index) const;
2013
2019 const Number *
2020 constraint_pool_end(const unsigned int pool_index) const;
2021
2025 const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
2026 get_shape_info(const unsigned int dof_handler_index_component = 0,
2027 const unsigned int quad_index = 0,
2028 const unsigned int fe_base_element = 0,
2029 const unsigned int hp_active_fe_index = 0,
2030 const unsigned int hp_active_quad_index = 0) const;
2031
2035 const internal::MatrixFreeFunctions::FaceToCellTopology<
2036 VectorizedArrayType::size()> &
2037 get_face_info(const unsigned int face_batch_index) const;
2038
2039
2045 const Table<3, unsigned int> &
2046 get_cell_and_face_to_plain_faces() const;
2047
2061 AlignedVector<VectorizedArrayType> *
2062 acquire_scratch_data() const;
2063
2067 void
2068 release_scratch_data(const AlignedVector<VectorizedArrayType> *memory) const;
2069
2079 AlignedVector<Number> *
2080 acquire_scratch_data_non_threadsafe() const;
2081
2085 void
2086 release_scratch_data_non_threadsafe(
2087 const AlignedVector<Number> *memory) const;
2088
2090
2091private:
2096 template <typename number2, int q_dim>
2097 void
2098 internal_reinit(
2099 const std::shared_ptr<hp::MappingCollection<dim>> & mapping,
2100 const std::vector<const DoFHandler<dim, dim> *> & dof_handlers,
2101 const std::vector<const AffineConstraints<number2> *> &constraint,
2102 const std::vector<IndexSet> & locally_owned_set,
2103 const std::vector<hp::QCollection<q_dim>> & quad,
2104 const AdditionalData & additional_data);
2105
2112 template <typename number2>
2113 void
2114 initialize_indices(
2115 const std::vector<const AffineConstraints<number2> *> &constraint,
2116 const std::vector<IndexSet> & locally_owned_set,
2117 const AdditionalData & additional_data);
2118
2122 void
2123 initialize_dof_handlers(
2124 const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
2125 const AdditionalData & additional_data);
2126
2130 std::vector<SmartPointer<const DoFHandler<dim>>> dof_handlers;
2131
2136 std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
2137
2144 std::vector<Number> constraint_pool_data;
2145
2150 std::vector<unsigned int> constraint_pool_row_index;
2151
2156 internal::MatrixFreeFunctions::MappingInfo<dim, Number, VectorizedArrayType>
2157 mapping_info;
2158
2162 Table<4, internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>>
2163 shape_info;
2164
2171 std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
2172
2173
2180 unsigned int cell_level_index_end_local;
2181
2187 internal::MatrixFreeFunctions::TaskInfo task_info;
2188
2193 internal::MatrixFreeFunctions::FaceInfo<VectorizedArrayType::size()>
2194 face_info;
2195
2199 bool indices_are_initialized;
2200
2204 bool mapping_is_initialized;
2205
2213 mutable Threads::ThreadLocalStorage<
2214 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2215 scratch_pad;
2216
2221 mutable std::list<std::pair<bool, AlignedVector<Number>>>
2222 scratch_pad_non_threadsafe;
2223
2227 unsigned int mg_level;
2228};
2229
2230
2231
2232/*----------------------- Inline functions ----------------------------------*/
2233
2234#ifndef DOXYGEN
2235
2236
2237
2238template <int dim, typename Number, typename VectorizedArrayType>
2239template <typename T>
2240inline void
2241MatrixFree<dim, Number, VectorizedArrayType>::initialize_cell_data_vector(
2242 AlignedVector<T> &vec) const
2243{
2244 vec.resize(this->n_cell_batches() + this->n_ghost_cell_batches());
2245}
2246
2247
2248
2249template <int dim, typename Number, typename VectorizedArrayType>
2250template <typename T>
2251inline void
2252MatrixFree<dim, Number, VectorizedArrayType>::initialize_face_data_vector(
2253 AlignedVector<T> &vec) const
2254{
2255 vec.resize(this->n_inner_face_batches() + this->n_boundary_face_batches() +
2256 this->n_ghost_inner_face_batches());
2257}
2258
2259
2260
2261template <int dim, typename Number, typename VectorizedArrayType>
2262template <typename VectorType>
2263inline void
2264MatrixFree<dim, Number, VectorizedArrayType>::initialize_dof_vector(
2265 VectorType & vec,
2266 const unsigned int comp) const
2267{
2268 static_assert(IsBlockVector<VectorType>::value == false,
2269 "This function is not supported for block vectors.");
2270
2271 Assert(task_info.n_procs == 1,
2272 ExcMessage("This function can only be used in serial."));
2273
2274 AssertIndexRange(comp, n_components());
2275 vec.reinit(dof_info[comp].vector_partitioner->size());
2276}
2277
2278
2279
2280template <int dim, typename Number, typename VectorizedArrayType>
2281template <typename Number2>
2282inline void
2283MatrixFree<dim, Number, VectorizedArrayType>::initialize_dof_vector(
2284 LinearAlgebra::distributed::Vector<Number2> &vec,
2285 const unsigned int comp) const
2286{
2287 AssertIndexRange(comp, n_components());
2288 vec.reinit(dof_info[comp].vector_partitioner, task_info.communicator_sm);
2289}
2290
2291
2292
2293template <int dim, typename Number, typename VectorizedArrayType>
2294inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2295MatrixFree<dim, Number, VectorizedArrayType>::get_vector_partitioner(
2296 const unsigned int comp) const
2297{
2298 AssertIndexRange(comp, n_components());
2299 return dof_info[comp].vector_partitioner;
2300}
2301
2302
2303
2304template <int dim, typename Number, typename VectorizedArrayType>
2305inline const std::vector<unsigned int> &
2306MatrixFree<dim, Number, VectorizedArrayType>::get_constrained_dofs(
2307 const unsigned int comp) const
2308{
2309 AssertIndexRange(comp, n_components());
2310 return dof_info[comp].constrained_dofs;
2311}
2312
2313
2314
2315template <int dim, typename Number, typename VectorizedArrayType>
2316inline unsigned int
2317MatrixFree<dim, Number, VectorizedArrayType>::n_components() const
2318{
2319 AssertDimension(dof_handlers.size(), dof_info.size());
2320 return dof_handlers.size();
2321}
2322
2323
2324
2325template <int dim, typename Number, typename VectorizedArrayType>
2326inline unsigned int
2327MatrixFree<dim, Number, VectorizedArrayType>::n_base_elements(
2328 const unsigned int dof_no) const
2329{
2330 AssertDimension(dof_handlers.size(), dof_info.size());
2331 AssertIndexRange(dof_no, dof_handlers.size());
2332 return dof_handlers[dof_no]->get_fe().n_base_elements();
2333}
2334
2335
2336
2337template <int dim, typename Number, typename VectorizedArrayType>
2338inline const internal::MatrixFreeFunctions::TaskInfo &
2339MatrixFree<dim, Number, VectorizedArrayType>::get_task_info() const
2340{
2341 return task_info;
2342}
2343
2344
2345
2346template <int dim, typename Number, typename VectorizedArrayType>
2347inline unsigned int
2348MatrixFree<dim, Number, VectorizedArrayType>::n_macro_cells() const
2349{
2350 return *(task_info.cell_partition_data.end() - 2);
2351}
2352
2353
2354
2355template <int dim, typename Number, typename VectorizedArrayType>
2356inline unsigned int
2357MatrixFree<dim, Number, VectorizedArrayType>::n_physical_cells() const
2358{
2359 return task_info.n_active_cells;
2360}
2361
2362
2363
2364template <int dim, typename Number, typename VectorizedArrayType>
2365inline unsigned int
2366MatrixFree<dim, Number, VectorizedArrayType>::n_cell_batches() const
2367{
2368 return *(task_info.cell_partition_data.end() - 2);
2369}
2370
2371
2372
2373template <int dim, typename Number, typename VectorizedArrayType>
2374inline unsigned int
2375MatrixFree<dim, Number, VectorizedArrayType>::n_ghost_cell_batches() const
2376{
2377 return *(task_info.cell_partition_data.end() - 1) -
2378 *(task_info.cell_partition_data.end() - 2);
2379}
2380
2381
2382
2383template <int dim, typename Number, typename VectorizedArrayType>
2384inline unsigned int
2385MatrixFree<dim, Number, VectorizedArrayType>::n_inner_face_batches() const
2386{
2387 if (task_info.face_partition_data.size() == 0)
2388 return 0;
2389 return task_info.face_partition_data.back();
2390}
2391
2392
2393
2394template <int dim, typename Number, typename VectorizedArrayType>
2395inline unsigned int
2396MatrixFree<dim, Number, VectorizedArrayType>::n_boundary_face_batches() const
2397{
2398 if (task_info.face_partition_data.size() == 0)
2399 return 0;
2400 return task_info.boundary_partition_data.back() -
2401 task_info.face_partition_data.back();
2402}
2403
2404
2405
2406template <int dim, typename Number, typename VectorizedArrayType>
2407inline unsigned int
2408MatrixFree<dim, Number, VectorizedArrayType>::n_ghost_inner_face_batches() const
2409{
2410 if (task_info.face_partition_data.size() == 0)
2411 return 0;
2412 return face_info.faces.size() - task_info.boundary_partition_data.back();
2413}
2414
2415
2416
2417template <int dim, typename Number, typename VectorizedArrayType>
2418inline types::boundary_id
2419MatrixFree<dim, Number, VectorizedArrayType>::get_boundary_id(
2420 const unsigned int face_batch_index) const
2421{
2422 Assert(face_batch_index >= task_info.boundary_partition_data[0] &&
2423 face_batch_index < task_info.boundary_partition_data.back(),
2424 ExcIndexRange(face_batch_index,
2425 task_info.boundary_partition_data[0],
2426 task_info.boundary_partition_data.back()));
2427 return types::boundary_id(face_info.faces[face_batch_index].exterior_face_no);
2428}
2429
2430
2431
2432template <int dim, typename Number, typename VectorizedArrayType>
2433inline std::array<types::boundary_id, VectorizedArrayType::size()>
2434MatrixFree<dim, Number, VectorizedArrayType>::get_faces_by_cells_boundary_id(
2435 const unsigned int cell_batch_index,
2436 const unsigned int face_number) const
2437{
2438 AssertIndexRange(cell_batch_index, n_cell_batches());
2439 AssertIndexRange(face_number, GeometryInfo<dim>::faces_per_cell);
2440 Assert(face_info.cell_and_face_boundary_id.size(0) >= n_cell_batches(),
2441 ExcNotInitialized());
2442 std::array<types::boundary_id, VectorizedArrayType::size()> result;
2443 result.fill(numbers::invalid_boundary_id);
2444 for (unsigned int v = 0;
2445 v < n_active_entries_per_cell_batch(cell_batch_index);
2446 ++v)
2447 result[v] =
2448 face_info.cell_and_face_boundary_id(cell_batch_index, face_number, v);
2449 return result;
2450}
2451
2452
2453
2454template <int dim, typename Number, typename VectorizedArrayType>
2455inline const internal::MatrixFreeFunctions::
2456 MappingInfo<dim, Number, VectorizedArrayType> &
2457 MatrixFree<dim, Number, VectorizedArrayType>::get_mapping_info() const
2458{
2459 return mapping_info;
2460}
2461
2462
2463
2464template <int dim, typename Number, typename VectorizedArrayType>
2465inline const internal::MatrixFreeFunctions::DoFInfo &
2466MatrixFree<dim, Number, VectorizedArrayType>::get_dof_info(
2467 const unsigned int dof_index) const
2468{
2469 AssertIndexRange(dof_index, n_components());
2470 return dof_info[dof_index];
2471}
2472
2473
2474
2475template <int dim, typename Number, typename VectorizedArrayType>
2476inline unsigned int
2477MatrixFree<dim, Number, VectorizedArrayType>::n_constraint_pool_entries() const
2478{
2479 return constraint_pool_row_index.size() - 1;
2480}
2481
2482
2483
2484template <int dim, typename Number, typename VectorizedArrayType>
2485inline const Number *
2486MatrixFree<dim, Number, VectorizedArrayType>::constraint_pool_begin(
2487 const unsigned int row) const
2488{
2489 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2490 return constraint_pool_data.empty() ?
2491 nullptr :
2492 constraint_pool_data.data() + constraint_pool_row_index[row];
2493}
2494
2495
2496
2497template <int dim, typename Number, typename VectorizedArrayType>
2498inline const Number *
2499MatrixFree<dim, Number, VectorizedArrayType>::constraint_pool_end(
2500 const unsigned int row) const
2501{
2502 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2503 return constraint_pool_data.empty() ?
2504 nullptr :
2505 constraint_pool_data.data() + constraint_pool_row_index[row + 1];
2506}
2507
2508
2509
2510template <int dim, typename Number, typename VectorizedArrayType>
2511inline std::pair<unsigned int, unsigned int>
2512MatrixFree<dim, Number, VectorizedArrayType>::create_cell_subrange_hp(
2513 const std::pair<unsigned int, unsigned int> &range,
2514 const unsigned int degree,
2515 const unsigned int dof_handler_component) const
2516{
2517 if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2518 {
2519 AssertDimension(
2520 dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2521 AssertDimension(
2522 dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2523 if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2524 return range;
2525 else
2526 return {range.second, range.second};
2527 }
2528
2529 const unsigned int fe_index =
2530 dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2531 if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2532 return {range.second, range.second};
2533 else
2534 return create_cell_subrange_hp_by_index(range,
2535 fe_index,
2536 dof_handler_component);
2537}
2538
2539
2540
2541template <int dim, typename Number, typename VectorizedArrayType>
2542inline bool
2543MatrixFree<dim, Number, VectorizedArrayType>::at_irregular_cell(
2544 const unsigned int cell_batch_index) const
2545{
2546 AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2547 return VectorizedArrayType::size() > 1 &&
2548 cell_level_index[(cell_batch_index + 1) * VectorizedArrayType::size() -
2549 1] == cell_level_index[(cell_batch_index + 1) *
2550 VectorizedArrayType::size() -
2551 2];
2552}
2553
2554
2555
2556template <int dim, typename Number, typename VectorizedArrayType>
2557unsigned int
2558MatrixFree<dim, Number, VectorizedArrayType>::n_active_fe_indices() const
2559{
2560 return shape_info.size(2);
2561}
2562
2563
2564template <int dim, typename Number, typename VectorizedArrayType>
2565unsigned int
2566MatrixFree<dim, Number, VectorizedArrayType>::get_cell_active_fe_index(
2567 const std::pair<unsigned int, unsigned int> range) const
2568{
2569 const auto &fe_indices = dof_info[0].cell_active_fe_index;
2570
2571 if (fe_indices.empty() == true)
2572 return 0;
2573
2574 const auto index = fe_indices[range.first];
2575
2576 for (unsigned int i = range.first; i < range.second; ++i)
2577 AssertDimension(index, fe_indices[i]);
2578
2579 return index;
2580}
2581
2582
2583
2584template <int dim, typename Number, typename VectorizedArrayType>
2585unsigned int
2586MatrixFree<dim, Number, VectorizedArrayType>::get_face_active_fe_index(
2587 const std::pair<unsigned int, unsigned int> range,
2588 const bool is_interior_face) const
2589{
2590 const auto &fe_indices = dof_info[0].cell_active_fe_index;
2591
2592 if (fe_indices.empty() == true)
2593 return 0;
2594
2595 if (is_interior_face)
2596 {
2597 const unsigned int index =
2598 fe_indices[face_info.faces[range.first].cells_interior[0] /
2599 VectorizedArrayType::size()];
2600
2601 for (unsigned int i = range.first; i < range.second; ++i)
2602 AssertDimension(index,
2603 fe_indices[face_info.faces[i].cells_interior[0] /
2604 VectorizedArrayType::size()]);
2605
2606 return index;
2607 }
2608 else
2609 {
2610 const unsigned int index =
2611 fe_indices[face_info.faces[range.first].cells_exterior[0] /
2612 VectorizedArrayType::size()];
2613
2614 for (unsigned int i = range.first; i < range.second; ++i)
2615 AssertDimension(index,
2616 fe_indices[face_info.faces[i].cells_exterior[0] /
2617 VectorizedArrayType::size()]);
2618
2619 return index;
2620 }
2621}
2622
2623
2624
2625template <int dim, typename Number, typename VectorizedArrayType>
2626inline unsigned int
2627MatrixFree<dim, Number, VectorizedArrayType>::n_components_filled(
2628 const unsigned int cell_batch_index) const
2629{
2630 return n_active_entries_per_cell_batch(cell_batch_index);
2631}
2632
2633
2634
2635template <int dim, typename Number, typename VectorizedArrayType>
2636inline unsigned int
2637MatrixFree<dim, Number, VectorizedArrayType>::n_active_entries_per_cell_batch(
2638 const unsigned int cell_batch_index) const
2639{
2640 Assert(!dof_info.empty(), ExcNotInitialized());
2641 AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2642 const std::vector<unsigned char> &n_lanes_filled =
2643 dof_info[0].n_vectorization_lanes_filled
2644 [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell];
2645 AssertIndexRange(cell_batch_index, n_lanes_filled.size());
2646
2647 return n_lanes_filled[cell_batch_index];
2648}
2649
2650
2651
2652template <int dim, typename Number, typename VectorizedArrayType>
2653inline unsigned int
2654MatrixFree<dim, Number, VectorizedArrayType>::n_active_entries_per_face_batch(
2655 const unsigned int face_batch_index) const
2656{
2657 AssertIndexRange(face_batch_index, face_info.faces.size());
2658 Assert(!dof_info.empty(), ExcNotInitialized());
2659 const std::vector<unsigned char> &n_lanes_filled =
2660 dof_info[0].n_vectorization_lanes_filled
2661 [internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior];
2662 AssertIndexRange(face_batch_index, n_lanes_filled.size());
2663 return n_lanes_filled[face_batch_index];
2664}
2665
2666
2667
2668template <int dim, typename Number, typename VectorizedArrayType>
2669inline unsigned int
2670MatrixFree<dim, Number, VectorizedArrayType>::get_dofs_per_cell(
2671 const unsigned int dof_handler_index,
2672 const unsigned int active_fe_index) const
2673{
2674 return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2675}
2676
2677
2678
2679template <int dim, typename Number, typename VectorizedArrayType>
2680inline unsigned int
2681MatrixFree<dim, Number, VectorizedArrayType>::get_n_q_points(
2682 const unsigned int quad_index,
2683 const unsigned int active_fe_index) const
2684{
2685 AssertIndexRange(quad_index, mapping_info.cell_data.size());
2686 return mapping_info.cell_data[quad_index]
2687 .descriptor[active_fe_index]
2688 .n_q_points;
2689}
2690
2691
2692
2693template <int dim, typename Number, typename VectorizedArrayType>
2694inline unsigned int
2695MatrixFree<dim, Number, VectorizedArrayType>::get_dofs_per_face(
2696 const unsigned int dof_handler_index,
2697 const unsigned int active_fe_index) const
2698{
2699 return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2700}
2701
2702
2703
2704template <int dim, typename Number, typename VectorizedArrayType>
2705inline unsigned int
2706MatrixFree<dim, Number, VectorizedArrayType>::get_n_q_points_face(
2707 const unsigned int quad_index,
2708 const unsigned int active_fe_index) const
2709{
2710 AssertIndexRange(quad_index, mapping_info.face_data.size());
2711 return mapping_info.face_data[quad_index]
2712 .descriptor[active_fe_index]
2713 .n_q_points;
2714}
2715
2716
2717
2718template <int dim, typename Number, typename VectorizedArrayType>
2719inline const IndexSet &
2720MatrixFree<dim, Number, VectorizedArrayType>::get_locally_owned_set(
2721 const unsigned int dof_handler_index) const
2722{
2723 return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2724}
2725
2726
2727
2728template <int dim, typename Number, typename VectorizedArrayType>
2729inline const IndexSet &
2730MatrixFree<dim, Number, VectorizedArrayType>::get_ghost_set(
2731 const unsigned int dof_handler_index) const
2732{
2733 return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2734}
2735
2736
2737
2738template <int dim, typename Number, typename VectorizedArrayType>
2739inline const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
2740MatrixFree<dim, Number, VectorizedArrayType>::get_shape_info(
2741 const unsigned int dof_handler_index,
2742 const unsigned int index_quad,
2743 const unsigned int index_fe,
2744 const unsigned int active_fe_index,
2745 const unsigned int active_quad_index) const
2746{
2747 AssertIndexRange(dof_handler_index, dof_info.size());
2748 const unsigned int ind =
2749 dof_info[dof_handler_index].global_base_element_offset + index_fe;
2750 AssertIndexRange(ind, shape_info.size(0));
2751 AssertIndexRange(index_quad, shape_info.size(1));
2752 AssertIndexRange(active_fe_index, shape_info.size(2));
2753 AssertIndexRange(active_quad_index, shape_info.size(3));
2754 return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2755}
2756
2757
2758
2759template <int dim, typename Number, typename VectorizedArrayType>
2760inline const internal::MatrixFreeFunctions::FaceToCellTopology<
2761 VectorizedArrayType::size()> &
2762MatrixFree<dim, Number, VectorizedArrayType>::get_face_info(
2763 const unsigned int face_batch_index) const
2764{
2765 AssertIndexRange(face_batch_index, face_info.faces.size());
2766 return face_info.faces[face_batch_index];
2767}
2768
2769
2770
2771template <int dim, typename Number, typename VectorizedArrayType>
2772inline const Table<3, unsigned int> &
2773MatrixFree<dim, Number, VectorizedArrayType>::get_cell_and_face_to_plain_faces()
2774 const
2775{
2776 return face_info.cell_and_face_to_plain_faces;
2777}
2778
2779
2780
2781template <int dim, typename Number, typename VectorizedArrayType>
2782inline const Quadrature<dim> &
2783MatrixFree<dim, Number, VectorizedArrayType>::get_quadrature(
2784 const unsigned int quad_index,
2785 const unsigned int active_fe_index) const
2786{
2787 AssertIndexRange(quad_index, mapping_info.cell_data.size());
2788 return mapping_info.cell_data[quad_index]
2789 .descriptor[active_fe_index]
2790 .quadrature;
2791}
2792
2793
2794
2795template <int dim, typename Number, typename VectorizedArrayType>
2796inline const Quadrature<dim - 1> &
2797MatrixFree<dim, Number, VectorizedArrayType>::get_face_quadrature(
2798 const unsigned int quad_index,
2799 const unsigned int active_fe_index) const
2800{
2801 AssertIndexRange(quad_index, mapping_info.face_data.size());
2802 return mapping_info.face_data[quad_index]
2803 .descriptor[active_fe_index]
2804 .quadrature;
2805}
2806
2807
2808
2809template <int dim, typename Number, typename VectorizedArrayType>
2810inline unsigned int
2811MatrixFree<dim, Number, VectorizedArrayType>::get_cell_range_category(
2812 const std::pair<unsigned int, unsigned int> range) const
2813{
2814 auto result = get_cell_category(range.first);
2815
2816 for (unsigned int i = range.first; i < range.second; ++i)
2817 result = std::max(result, get_cell_category(i));
2818
2819 return result;
2820}
2821
2822
2823
2824template <int dim, typename Number, typename VectorizedArrayType>
2825inline std::pair<unsigned int, unsigned int>
2826MatrixFree<dim, Number, VectorizedArrayType>::get_face_range_category(
2827 const std::pair<unsigned int, unsigned int> range) const
2828{
2829 auto result = get_face_category(range.first);
2830
2831 for (unsigned int i = range.first; i < range.second; ++i)
2832 {
2833 result.first = std::max(result.first, get_face_category(i).first);
2834 result.second = std::max(result.second, get_face_category(i).second);
2835 }
2836
2837 return result;
2838}
2839
2840
2841
2842template <int dim, typename Number, typename VectorizedArrayType>
2843inline unsigned int
2844MatrixFree<dim, Number, VectorizedArrayType>::get_cell_category(
2845 const unsigned int cell_batch_index) const
2846{
2847 AssertIndexRange(0, dof_info.size());
2848 AssertIndexRange(cell_batch_index, dof_info[0].cell_active_fe_index.size());
2849 if (dof_info[0].cell_active_fe_index.empty())
2850 return 0;
2851 else
2852 return dof_info[0].cell_active_fe_index[cell_batch_index];
2853}
2854
2855
2856
2857template <int dim, typename Number, typename VectorizedArrayType>
2858inline std::pair<unsigned int, unsigned int>
2859MatrixFree<dim, Number, VectorizedArrayType>::get_face_category(
2860 const unsigned int face_batch_index) const
2861{
2862 AssertIndexRange(face_batch_index, face_info.faces.size());
2863 if (dof_info[0].cell_active_fe_index.empty())
2864 return std::make_pair(0U, 0U);
2865
2866 std::pair<unsigned int, unsigned int> result = std::make_pair(0U, 0U);
2867 for (unsigned int v = 0;
2868 v < VectorizedArrayType::size() &&
2869 face_info.faces[face_batch_index].cells_interior[v] !=
2870 numbers::invalid_unsigned_int;
2871 ++v)
2872 result.first = std::max(
2873 result.first,
2874 dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
2875 .cells_interior[v] /
2876 VectorizedArrayType::size()]);
2877 if (face_info.faces[face_batch_index].cells_exterior[0] !=
2878 numbers::invalid_unsigned_int)
2879 for (unsigned int v = 0;
2880 v < VectorizedArrayType::size() &&
2881 face_info.faces[face_batch_index].cells_exterior[v] !=
2882 numbers::invalid_unsigned_int;
2883 ++v)
2884 result.second = std::max(
2885 result.second,
2886 dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
2887 .cells_exterior[v] /
2888 VectorizedArrayType::size()]);
2889 else
2890 result.second = numbers::invalid_unsigned_int;
2891 return result;
2892}
2893
2894
2895
2896template <int dim, typename Number, typename VectorizedArrayType>
2897inline bool
2898MatrixFree<dim, Number, VectorizedArrayType>::indices_initialized() const
2899{
2900 return indices_are_initialized;
2901}
2902
2903
2904
2905template <int dim, typename Number, typename VectorizedArrayType>
2906inline bool
2907MatrixFree<dim, Number, VectorizedArrayType>::mapping_initialized() const
2908{
2909 return mapping_is_initialized;
2910}
2911
2912
2913template <int dim, typename Number, typename VectorizedArrayType>
2914inline unsigned int
2915MatrixFree<dim, Number, VectorizedArrayType>::get_mg_level() const
2916{
2917 return mg_level;
2918}
2919
2920
2921
2922template <int dim, typename Number, typename VectorizedArrayType>
2923AlignedVector<VectorizedArrayType> *
2924MatrixFree<dim, Number, VectorizedArrayType>::acquire_scratch_data() const
2925{
2926 using list_type =
2927 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
2928 list_type &data = scratch_pad.get();
2929 for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
2930 if (it->first == false)
2931 {
2932 it->first = true;
2933 return &it->second;
2934 }
2935 data.emplace_front(true, AlignedVector<VectorizedArrayType>());
2936 return &data.front().second;
2937}
2938
2939
2940
2941template <int dim, typename Number, typename VectorizedArrayType>
2942void
2943MatrixFree<dim, Number, VectorizedArrayType>::release_scratch_data(
2944 const AlignedVector<VectorizedArrayType> *scratch) const
2945{
2946 using list_type =
2947 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
2948 list_type &data = scratch_pad.get();
2949 for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
2950 if (&it->second == scratch)
2951 {
2952 Assert(it->first == true, ExcInternalError());
2953 it->first = false;
2954 return;
2955 }
2956 AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2957}
2958
2959
2960
2961template <int dim, typename Number, typename VectorizedArrayType>
2962AlignedVector<Number> *
2963MatrixFree<dim, Number, VectorizedArrayType>::
2964 acquire_scratch_data_non_threadsafe() const
2965{
2966 for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
2967 scratch_pad_non_threadsafe.begin();
2968 it != scratch_pad_non_threadsafe.end();
2969 ++it)
2970 if (it->first == false)
2971 {
2972 it->first = true;
2973 return &it->second;
2974 }
2975 scratch_pad_non_threadsafe.push_front(
2976 std::make_pair(true, AlignedVector<Number>()));
2977 return &scratch_pad_non_threadsafe.front().second;
2978}
2979
2980
2981
2982template <int dim, typename Number, typename VectorizedArrayType>
2983void
2984MatrixFree<dim, Number, VectorizedArrayType>::
2985 release_scratch_data_non_threadsafe(
2986 const AlignedVector<Number> *scratch) const
2987{
2988 for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
2989 scratch_pad_non_threadsafe.begin();
2990 it != scratch_pad_non_threadsafe.end();
2991 ++it)
2992 if (&it->second == scratch)
2993 {
2994 Assert(it->first == true, ExcInternalError());
2995 it->first = false;
2996 return;
2997 }
2998 AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2999}
3000
3001
3002
3003// ------------------------------ reinit functions ---------------------------
3004
3005namespace internal
3006{
3007 namespace MatrixFreeImplementation
3008 {
3009 template <int dim, int spacedim>
3010 inline std::vector<IndexSet>
3011 extract_locally_owned_index_sets(
3012 const std::vector<const ::DoFHandler<dim, spacedim> *> &dofh,
3013 const unsigned int level)
3014 {
3015 std::vector<IndexSet> locally_owned_set;
3016 locally_owned_set.reserve(dofh.size());
3017 for (unsigned int j = 0; j < dofh.size(); ++j)
3018 if (level == numbers::invalid_unsigned_int)
3019 locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
3020 else
3021 locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
3022 return locally_owned_set;
3023 }
3024 } // namespace MatrixFreeImplementation
3025} // namespace internal
3026
3027
3028
3029template <int dim, typename Number, typename VectorizedArrayType>
3030template <typename QuadratureType, typename number2>
3031void
3032MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3033 const DoFHandler<dim> & dof_handler,
3034 const AffineConstraints<number2> &constraints_in,
3035 const QuadratureType & quad,
3036 const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
3037 &additional_data)
3038{
3039 std::vector<const DoFHandler<dim, dim> *> dof_handlers;
3040 std::vector<const AffineConstraints<number2> *> constraints;
3041 std::vector<QuadratureType> quads;
3042
3043 dof_handlers.push_back(&dof_handler);
3044 constraints.push_back(&constraints_in);
3045 quads.push_back(quad);
3046
3047 std::vector<IndexSet> locally_owned_sets =
3048 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3049 dof_handlers, additional_data.mg_level);
3050
3051 std::vector<hp::QCollection<dim>> quad_hp;
3052 quad_hp.emplace_back(quad);
3053
3054 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(
3055 StaticMappingQ1<dim>::mapping),
3056 dof_handlers,
3057 constraints,
3058 locally_owned_sets,
3059 quad_hp,
3060 additional_data);
3061}
3062
3063
3064
3065template <int dim, typename Number, typename VectorizedArrayType>
3066template <typename QuadratureType, typename number2, typename MappingType>
3067void
3068MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3069 const MappingType & mapping,
3070 const DoFHandler<dim> & dof_handler,
3071 const AffineConstraints<number2> &constraints_in,
3072 const QuadratureType & quad,
3073 const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
3074 &additional_data)
3075{
3076 std::vector<const DoFHandler<dim, dim> *> dof_handlers;
3077 std::vector<const AffineConstraints<number2> *> constraints;
3078
3079 dof_handlers.push_back(&dof_handler);
3080 constraints.push_back(&constraints_in);
3081
3082 std::vector<IndexSet> locally_owned_sets =
3083 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3084 dof_handlers, additional_data.mg_level);
3085
3086 std::vector<hp::QCollection<dim>> quad_hp;
3087 quad_hp.emplace_back(quad);
3088
3089 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3090 dof_handlers,
3091 constraints,
3092 locally_owned_sets,
3093 quad_hp,
3094 additional_data);
3095}
3096
3097
3098
3099template <int dim, typename Number, typename VectorizedArrayType>
3100template <typename QuadratureType, typename number2>
3101void
3102MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3103 const std::vector<const DoFHandler<dim> *> & dof_handler,
3104 const std::vector<const AffineConstraints<number2> *> &constraint,
3105 const std::vector<QuadratureType> & quad,
3106 const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
3107 &additional_data)
3108{
3109 std::vector<IndexSet> locally_owned_set =
3110 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3111 dof_handler, additional_data.mg_level);
3112 std::vector<hp::QCollection<dim>> quad_hp;
3113 for (unsigned int q = 0; q < quad.size(); ++q)
3114 quad_hp.emplace_back(quad[q]);
3115
3116 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(
3117 StaticMappingQ1<dim>::mapping),
3118 dof_handler,
3119 constraint,
3120 locally_owned_set,
3121 quad_hp,
3122 additional_data);
3123}
3124
3125
3126
3127template <int dim, typename Number, typename VectorizedArrayType>
3128template <typename QuadratureType,
3129 typename number2,
3130 typename DoFHandlerType,
3131 typename MappingType>
3132void
3133MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3134 const MappingType & mapping,
3135 const std::vector<const DoFHandlerType *> & dof_handler,
3136 const std::vector<const AffineConstraints<number2> *> &constraint,
3137 const std::vector<QuadratureType> & quad,
3138 const AdditionalData & additional_data)
3139{
3140 static_assert(dim == DoFHandlerType::dimension,
3141 "Dimension dim not equal to DoFHandlerType::dimension.");
3142
3143 std::vector<const DoFHandler<dim> *> dof_handlers;
3144
3145 for (const auto dh : dof_handler)
3146 dof_handlers.push_back(dh);
3147
3148 this->reinit(mapping, dof_handlers, constraint, quad, additional_data);
3149}
3150
3151
3152
3153template <int dim, typename Number, typename VectorizedArrayType>
3154template <typename QuadratureType, typename number2>
3155void
3156MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3157 const std::vector<const DoFHandler<dim> *> & dof_handler,
3158 const std::vector<const AffineConstraints<number2> *> &constraint,
3159 const QuadratureType & quad,
3160 const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
3161 &additional_data)
3162{
3163 std::vector<IndexSet> locally_owned_set =
3164 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3165 dof_handler, additional_data.mg_level);
3166 std::vector<hp::QCollection<dim>> quad_hp;
3167 quad_hp.emplace_back(quad);
3168
3169 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(
3170 StaticMappingQ1<dim>::mapping),
3171 dof_handler,
3172 constraint,
3173 locally_owned_set,
3174 quad_hp,
3175 additional_data);
3176}
3177
3178
3179
3180template <int dim, typename Number, typename VectorizedArrayType>
3181template <typename QuadratureType, typename number2, typename DoFHandlerType>
3182void
3183MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3184 const std::vector<const DoFHandlerType *> & dof_handler,
3185 const std::vector<const AffineConstraints<number2> *> &constraint,
3186 const std::vector<QuadratureType> & quad,
3187 const AdditionalData & additional_data)
3188{
3189 static_assert(dim == DoFHandlerType::dimension,
3190 "Dimension dim not equal to DoFHandlerType::dimension.");
3191
3192 std::vector<const DoFHandler<dim> *> dof_handlers;
3193
3194 for (const auto dh : dof_handler)
3195 dof_handlers.push_back(dh);
3196
3197 this->reinit(dof_handlers, constraint, quad, additional_data);
3198}
3199
3200
3201
3202template <int dim, typename Number, typename VectorizedArrayType>
3203template <typename QuadratureType, typename number2, typename MappingType>
3204void
3205MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3206 const MappingType & mapping,
3207 const std::vector<const DoFHandler<dim> *> & dof_handler,
3208 const std::vector<const AffineConstraints<number2> *> &constraint,
3209 const QuadratureType & quad,
3210 const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
3211 &additional_data)
3212{
3213 std::vector<IndexSet> locally_owned_set =
3214 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3215 dof_handler, additional_data.mg_level);
3216 std::vector<hp::QCollection<dim>> quad_hp;
3217 quad_hp.emplace_back(quad);
3218
3219 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3220 dof_handler,
3221 constraint,
3222 locally_owned_set,
3223 quad_hp,
3224 additional_data);
3225}
3226
3227
3228
3229template <int dim, typename Number, typename VectorizedArrayType>
3230template <typename QuadratureType, typename number2, typename MappingType>
3231void
3232MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3233 const MappingType & mapping,
3234 const std::vector<const DoFHandler<dim> *> & dof_handler,
3235 const std::vector<const AffineConstraints<number2> *> &constraint,
3236 const std::vector<QuadratureType> & quad,
3237 const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
3238 &additional_data)
3239{
3240 std::vector<IndexSet> locally_owned_set =
3241 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3242 dof_handler, additional_data.mg_level);
3243 std::vector<hp::QCollection<dim>> quad_hp;
3244 for (unsigned int q = 0; q < quad.size(); ++q)
3245 quad_hp.emplace_back(quad[q]);
3246
3247 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3248 dof_handler,
3249 constraint,
3250 locally_owned_set,
3251 quad_hp,
3252 additional_data);
3253}
3254
3255
3256
3257template <int dim, typename Number, typename VectorizedArrayType>
3258template <typename QuadratureType,
3259 typename number2,
3260 typename DoFHandlerType,
3261 typename MappingType>
3262void
3263MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3264 const MappingType & mapping,
3265 const std::vector<const DoFHandlerType *> & dof_handler,
3266 const std::vector<const AffineConstraints<number2> *> &constraint,
3267 const QuadratureType & quad,
3268 const AdditionalData & additional_data)
3269{
3270 static_assert(dim == DoFHandlerType::dimension,
3271 "Dimension dim not equal to DoFHandlerType::dimension.");
3272
3273 std::vector<const DoFHandler<dim> *> dof_handlers;
3274
3275 for (const auto dh : dof_handler)
3276 dof_handlers.push_back(dh);
3277
3278 this->reinit(mapping, dof_handlers, constraint, quad, additional_data);
3279}
3280
3281
3282
3283template <int dim, typename Number, typename VectorizedArrayType>
3284template <typename QuadratureType, typename number2, typename DoFHandlerType>
3285void
3286MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3287 const std::vector<const DoFHandlerType *> & dof_handler,
3288 const std::vector<const AffineConstraints<number2> *> &constraint,
3289 const QuadratureType & quad,
3290 const AdditionalData & additional_data)
3291{
3292 static_assert(dim == DoFHandlerType::dimension,
3293 "Dimension dim not equal to DoFHandlerType::dimension.");
3294
3295 std::vector<const DoFHandler<dim> *> dof_handlers;
3296
3297 for (const auto dh : dof_handler)
3298 dof_handlers.push_back(dh);
3299
3300 this->reinit(dof_handlers, constraint, quad, additional_data);
3301}
3302
3303
3304
3305// ------------------------------ implementation of loops --------------------
3306
3307// internal helper functions that define how to call MPI data exchange
3308// functions: for generic vectors, do nothing at all. For distributed vectors,
3309// call update_ghost_values_start function and so on. If we have collections
3310// of vectors, just do the individual functions of the components. In order to
3311// keep ghost values consistent (whether we are in read or write mode), we
3312// also reset the values at the end. the whole situation is a bit complicated
3313// by the fact that we need to treat block vectors differently, which use some
3314// additional helper functions to select the blocks and template magic.
3315namespace internal
3316{
3320 template <int dim, typename Number, typename VectorizedArrayType>
3321 struct VectorDataExchange
3322 {
3323 // A shift for the MPI messages to reduce the risk for accidental
3324 // interaction with other open communications that a user program might
3325 // set up (parallel vectors support unfinished communication). We let
3326 // the other vectors use the first 20 assigned numbers and start the
3327 // matrix-free communication.
3328 static constexpr unsigned int channel_shift = 20;
3329
3330
3331
3336 VectorDataExchange(
3337 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3338 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3339 DataAccessOnFaces vector_face_access,
3340 const unsigned int n_components)
3341 : matrix_free(matrix_free)
3342 , vector_face_access(
3343 matrix_free.get_task_info().face_partition_data.empty() ?
3344 ::MatrixFree<dim, Number, VectorizedArrayType>::
3345 DataAccessOnFaces::unspecified :
3346 vector_face_access)
3347 , ghosts_were_set(false)
3348# ifdef DEAL_II_WITH_MPI
3349 , tmp_data(n_components)
3350 , requests(n_components)
3351# endif
3352 {
3353 (void)n_components;
3354 if (this->vector_face_access !=
3355 ::MatrixFree<dim, Number, VectorizedArrayType>::
3356 DataAccessOnFaces::unspecified)
3357 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3358 AssertDimension(
3359 matrix_free.get_dof_info(c).vector_exchanger_face_variants.size(),
3360 5);
3361 }
3362
3363
3364
3368 ~VectorDataExchange() // NOLINT
3369 {
3370# ifdef DEAL_II_WITH_MPI
3371 for (unsigned int i = 0; i < tmp_data.size(); ++i)
3372 if (tmp_data[i] != nullptr)
3373 matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3374# endif
3375 }
3376
3377
3378
3383 template <typename VectorType>
3384 unsigned int
3385 find_vector_in_mf(const VectorType &vec,
3386 const bool check_global_compatibility = true) const
3387 {
3388 // case 1: vector was set up with MatrixFree::initialize_dof_vector()
3389 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3390 if (vec.get_partitioner().get() ==
3391 matrix_free.get_dof_info(c).vector_partitioner.get())
3392 return c;
3393
3394 // case 2: user provided own partitioner (compatibility mode)
3395 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3396 if (check_global_compatibility ?
3397 vec.get_partitioner()->is_globally_compatible(
3398 *matrix_free.get_dof_info(c).vector_partitioner) :
3399 vec.get_partitioner()->is_compatible(
3400 *matrix_free.get_dof_info(c).vector_partitioner))
3401 return c;
3402
3403 Assert(false,
3404 ExcNotImplemented("Could not find partitioner that fits vector"));
3405
3406 return numbers::invalid_unsigned_int;
3407 }
3408
3409
3410
3415 const internal::MatrixFreeFunctions::VectorDataExchange::Base &
3416 get_partitioner(const unsigned int mf_component) const
3417 {
3418 AssertDimension(matrix_free.get_dof_info(mf_component)
3419 .vector_exchanger_face_variants.size(),
3420 5);
3421 if (vector_face_access ==
3422 ::MatrixFree<dim, Number, VectorizedArrayType>::
3423 DataAccessOnFaces::none)
3424 return *matrix_free.get_dof_info(mf_component)
3425 .vector_exchanger_face_variants[0];
3426 else if (vector_face_access ==
3427 ::MatrixFree<dim, Number, VectorizedArrayType>::
3428 DataAccessOnFaces::values)
3429 return *matrix_free.get_dof_info(mf_component)
3430 .vector_exchanger_face_variants[1];
3431 else if (vector_face_access ==
3432 ::MatrixFree<dim, Number, VectorizedArrayType>::
3433 DataAccessOnFaces::gradients)
3434 return *matrix_free.get_dof_info(mf_component)
3435 .vector_exchanger_face_variants[2];
3436 else if (vector_face_access ==
3437 ::MatrixFree<dim, Number, VectorizedArrayType>::
3438 DataAccessOnFaces::values_all_faces)
3439 return *matrix_free.get_dof_info(mf_component)
3440 .vector_exchanger_face_variants[3];
3441 else if (vector_face_access ==
3442 ::MatrixFree<dim, Number, VectorizedArrayType>::
3443 DataAccessOnFaces::gradients_all_faces)
3444 return *matrix_free.get_dof_info(mf_component)
3445 .vector_exchanger_face_variants[4];
3446 else
3447 return *matrix_free.get_dof_info(mf_component).vector_exchanger.get();
3448 }
3449
3450
3451
3455 template <typename VectorType,
3456 typename std::enable_if<is_not_parallel_vector<VectorType>,
3457 VectorType>::type * = nullptr>
3458 void
3459 update_ghost_values_start(const unsigned int /*component_in_block_vector*/,
3460 const VectorType & /*vec*/)
3461 {}
3462
3463
3468 template <
3469 typename VectorType,
3470 typename std::enable_if<!has_update_ghost_values_start<VectorType> &&
3471 !is_not_parallel_vector<VectorType>,
3472 VectorType>::type * = nullptr>
3473 void
3474 update_ghost_values_start(const unsigned int component_in_block_vector,
3475 const VectorType & vec)
3476 {
3477 (void)component_in_block_vector;
3478 bool ghosts_set = vec.has_ghost_elements();
3479
3480 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3481 ghosts_set == false,
3482 ExcNotImplemented());
3483
3484 if (ghosts_set)
3485 ghosts_were_set = true;
3486
3487 vec.update_ghost_values();
3488 }
3489
3490
3491
3497 template <
3498 typename VectorType,
3499 typename std::enable_if<has_update_ghost_values_start<VectorType> &&
3500 !has_exchange_on_subset<VectorType>,
3501 VectorType>::type * = nullptr>
3502 void
3503 update_ghost_values_start(const unsigned int component_in_block_vector,
3504 const VectorType & vec)
3505 {
3506 (void)component_in_block_vector;
3507 bool ghosts_set = vec.has_ghost_elements();
3508
3509 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3510 ghosts_set == false,
3511 ExcNotImplemented());
3512
3513 if (ghosts_set)
3514 ghosts_were_set = true;
3515
3516 vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3517 }
3518
3519
3520
3527 template <
3528 typename VectorType,
3529 typename std::enable_if<has_update_ghost_values_start<VectorType> &&
3530 has_exchange_on_subset<VectorType>,
3531 VectorType>::type * = nullptr>
3532 void
3533 update_ghost_values_start(const unsigned int component_in_block_vector,
3534 const VectorType & vec)
3535 {
3536 static_assert(
3537 std::is_same<Number, typename VectorType::value_type>::value,
3538 "Type mismatch between VectorType and VectorDataExchange");
3539 (void)component_in_block_vector;
3540 bool ghosts_set = vec.has_ghost_elements();
3541
3542 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3543 ghosts_set == false,
3544 ExcNotImplemented());
3545
3546 if (ghosts_set)
3547 ghosts_were_set = true;
3548
3549 if (vec.size() != 0)
3550 {
3551# ifdef DEAL_II_WITH_MPI
3552 const unsigned int mf_component = find_vector_in_mf(vec);
3553
3554 const auto &part = get_partitioner(mf_component);
3555
3556 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3557 part.n_import_sm_procs() == 0)
3558 return;
3559
3560 tmp_data[component_in_block_vector] =
3561 matrix_free.acquire_scratch_data_non_threadsafe();
3562 tmp_data[component_in_block_vector]->resize_fast(
3563 part.n_import_indices());
3564 AssertDimension(requests.size(), tmp_data.size());
3565
3566 part.export_to_ghosted_array_start(
3567 component_in_block_vector * 2 + channel_shift,
3568 ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3569 vec.shared_vector_data(),
3570 ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3571 part.locally_owned_size(),
3572 matrix_free.get_dof_info(mf_component)
3573 .vector_partitioner->n_ghost_indices()),
3574 ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3575 part.n_import_indices()),
3576 this->requests[component_in_block_vector]);
3577# endif
3578 }
3579 }
3580
3581
3582
3587 template <
3588 typename VectorType,
3589 typename std::enable_if<!has_update_ghost_values_start<VectorType>,
3590 VectorType>::type * = nullptr>
3591 void
3592 update_ghost_values_finish(const unsigned int /*component_in_block_vector*/,
3593 const VectorType & /*vec*/)
3594 {}
3595
3596
3597
3603 template <
3604 typename VectorType,
3605 typename std::enable_if<has_update_ghost_values_start<VectorType> &&
3606 !has_exchange_on_subset<VectorType>,
3607 VectorType>::type * = nullptr>
3608 void
3609 update_ghost_values_finish(const unsigned int component_in_block_vector,
3610 const VectorType & vec)
3611 {
3612 (void)component_in_block_vector;
3613 vec.update_ghost_values_finish();
3614 }
3615
3616
3617
3624 template <
3625 typename VectorType,
3626 typename std::enable_if<has_update_ghost_values_start<VectorType> &&
3627 has_exchange_on_subset<VectorType>,
3628 VectorType>::type * = nullptr>
3629 void
3630 update_ghost_values_finish(const unsigned int component_in_block_vector,
3631 const VectorType & vec)
3632 {
3633 static_assert(
3634 std::is_same<Number, typename VectorType::value_type>::value,
3635 "Type mismatch between VectorType and VectorDataExchange");
3636 (void)component_in_block_vector;
3637
3638 if (vec.size() != 0)
3639 {
3640# ifdef DEAL_II_WITH_MPI
3641 AssertIndexRange(component_in_block_vector, tmp_data.size());
3642 AssertDimension(requests.size(), tmp_data.size());
3643
3644 const unsigned int mf_component = find_vector_in_mf(vec);
3645
3646 const auto &part = get_partitioner(mf_component);
3647
3648 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3649 part.n_import_sm_procs() != 0)
3650 {
3651 part.export_to_ghosted_array_finish(
3652 ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3653 vec.shared_vector_data(),
3654 ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3655 part.locally_owned_size(),
3656 matrix_free.get_dof_info(mf_component)
3657 .vector_partitioner->n_ghost_indices()),
3658 this->requests[component_in_block_vector]);
3659
3660 matrix_free.release_scratch_data_non_threadsafe(
3661 tmp_data[component_in_block_vector]);
3662 tmp_data[component_in_block_vector] = nullptr;
3663 }
3664# endif
3665 }
3666 // let vector know that ghosts are being updated and we can read from
3667 // them
3668 vec.set_ghost_state(true);
3669 }
3670
3671
3672
3676 template <typename VectorType,
3677 typename std::enable_if<is_not_parallel_vector<VectorType>,
3678 VectorType>::type * = nullptr>
3679 void
3680 compress_start(const unsigned int /*component_in_block_vector*/,
3681 VectorType & /*vec*/)
3682 {}
3683
3684
3685
3690 template <typename VectorType,
3691 typename std::enable_if<!has_compress_start<VectorType> &&
3692 !is_not_parallel_vector<VectorType>,
3693 VectorType>::type * = nullptr>
3694 void
3695 compress_start(const unsigned int component_in_block_vector,
3696 VectorType & vec)
3697 {
3698 (void)component_in_block_vector;
3699 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3700 vec.compress(::VectorOperation::add);
3701 }
3702
3703
3704
3710 template <typename VectorType,
3711 typename std::enable_if<has_compress_start<VectorType> &&
3712 !has_exchange_on_subset<VectorType>,
3713 VectorType>::type * = nullptr>
3714 void
3715 compress_start(const unsigned int component_in_block_vector,
3716 VectorType & vec)
3717 {
3718 (void)component_in_block_vector;
3719 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3720 vec.compress_start(component_in_block_vector + channel_shift);
3721 }
3722
3723
3724
3731 template <typename VectorType,
3732 typename std::enable_if<has_compress_start<VectorType> &&
3733 has_exchange_on_subset<VectorType>,
3734 VectorType>::type * = nullptr>
3735 void
3736 compress_start(const unsigned int component_in_block_vector,
3737 VectorType & vec)
3738 {
3739 static_assert(
3740 std::is_same<Number, typename VectorType::value_type>::value,
3741 "Type mismatch between VectorType and VectorDataExchange");
3742 (void)component_in_block_vector;
3743 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3744
3745 if (vec.size() != 0)
3746 {
3747# ifdef DEAL_II_WITH_MPI
3748 const unsigned int mf_component = find_vector_in_mf(vec);
3749
3750 const auto &part = get_partitioner(mf_component);
3751
3752 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3753 part.n_import_sm_procs() == 0)
3754 return;
3755
3756 tmp_data[component_in_block_vector] =
3757 matrix_free.acquire_scratch_data_non_threadsafe();
3758 tmp_data[component_in_block_vector]->resize_fast(
3759 part.n_import_indices());
3760 AssertDimension(requests.size(), tmp_data.size());
3761
3762 part.import_from_ghosted_array_start(
3763 ::VectorOperation::add,
3764 component_in_block_vector * 2 + channel_shift,
3765 ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3766 vec.shared_vector_data(),
3767 ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3768 matrix_free.get_dof_info(mf_component)
3769 .vector_partitioner->n_ghost_indices()),
3770 ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3771 part.n_import_indices()),
3772 this->requests[component_in_block_vector]);
3773# endif
3774 }
3775 }
3776
3777
3778
3783 template <typename VectorType,
3784 typename std::enable_if<!has_compress_start<VectorType>,
3785 VectorType>::type * = nullptr>
3786 void
3787 compress_finish(const unsigned int /*component_in_block_vector*/,
3788 VectorType & /*vec*/)
3789 {}
3790
3791
3792
3798 template <typename VectorType,
3799 typename std::enable_if<has_compress_start<VectorType> &&
3800 !has_exchange_on_subset<VectorType>,
3801 VectorType>::type * = nullptr>
3802 void
3803 compress_finish(const unsigned int component_in_block_vector,
3804 VectorType & vec)
3805 {
3806 (void)component_in_block_vector;
3807 vec.compress_finish(::VectorOperation::add);
3808 }
3809
3810
3811
3818 template <typename VectorType,
3819 typename std::enable_if<has_compress_start<VectorType> &&
3820 has_exchange_on_subset<VectorType>,
3821 VectorType>::type * = nullptr>
3822 void
3823 compress_finish(const unsigned int component_in_block_vector,
3824 VectorType & vec)
3825 {
3826 static_assert(
3827 std::is_same<Number, typename VectorType::value_type>::value,
3828 "Type mismatch between VectorType and VectorDataExchange");
3829 (void)component_in_block_vector;
3830 if (vec.size() != 0)
3831 {
3832# ifdef DEAL_II_WITH_MPI
3833 AssertIndexRange(component_in_block_vector, tmp_data.size());
3834 AssertDimension(requests.size(), tmp_data.size());
3835
3836 const unsigned int mf_component = find_vector_in_mf(vec);
3837
3838 const auto &part = get_partitioner(mf_component);
3839
3840 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3841 part.n_import_sm_procs() != 0)
3842 {
3843 part.import_from_ghosted_array_finish(
3844 VectorOperation::add,
3845 ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3846 vec.shared_vector_data(),
3847 ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3848 matrix_free.get_dof_info(mf_component)
3849 .vector_partitioner->n_ghost_indices()),
3850 ArrayView<const Number>(
3851 tmp_data[component_in_block_vector]->begin(),
3852 part.n_import_indices()),
3853 this->requests[component_in_block_vector]);
3854
3855 matrix_free.release_scratch_data_non_threadsafe(
3856 tmp_data[component_in_block_vector]);
3857 tmp_data[component_in_block_vector] = nullptr;
3858 }
3859
3860 if (Utilities::MPI::job_supports_mpi())
3861 {
3862 const int ierr =
3863 MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3864 AssertThrowMPI(ierr);
3865 }
3866# endif
3867 }
3868 }
3869
3870
3871
3875 template <typename VectorType,
3876 typename std::enable_if<is_not_parallel_vector<VectorType>,
3877 VectorType>::type * = nullptr>
3878 void
3879 reset_ghost_values(const VectorType & /*vec*/) const
3880 {}
3881
3882
3883
3888 template <typename VectorType,
3889 typename std::enable_if<!has_exchange_on_subset<VectorType> &&
3890 !is_not_parallel_vector<VectorType>,
3891 VectorType>::type * = nullptr>
3892 void
3893 reset_ghost_values(const VectorType &vec) const
3894 {
3895 if (ghosts_were_set == true)
3896 return;
3897
3898 vec.zero_out_ghost_values();
3899 }
3900
3901
3902
3908 template <typename VectorType,
3909 typename std::enable_if<has_exchange_on_subset<VectorType>,
3910 VectorType>::type * = nullptr>
3911 void
3912 reset_ghost_values(const VectorType &vec) const
3913 {
3914 static_assert(
3915 std::is_same<Number, typename VectorType::value_type>::value,
3916 "Type mismatch between VectorType and VectorDataExchange");
3917 if (ghosts_were_set == true)
3918 return;
3919
3920 if (vec.size() != 0)
3921 {
3922# ifdef DEAL_II_WITH_MPI
3923 AssertDimension(requests.size(), tmp_data.size());
3924
3925 const unsigned int mf_component = find_vector_in_mf(vec);
3926
3927 const auto &part = get_partitioner(mf_component);
3928
3929 if (part.n_ghost_indices() > 0)
3930 {
3931 part.reset_ghost_values(ArrayView<Number>(
3932 const_cast<LinearAlgebra::distributed::Vector<Number> &>(vec)
3933 .begin() +
3934 part.locally_owned_size(),
3935 matrix_free.get_dof_info(mf_component)
3936 .vector_partitioner->n_ghost_indices()));
3937 }
3938
3939# endif
3940 }
3941 // let vector know that it's not ghosted anymore
3942 vec.set_ghost_state(false);
3943 }
3944
3945
3946
3952 template <typename VectorType,
3953 typename std::enable_if<has_exchange_on_subset<VectorType>,
3954 VectorType>::type * = nullptr>
3955 void
3956 zero_vector_region(const unsigned int range_index, VectorType &vec) const
3957 {
3958 static_assert(
3959 std::is_same<Number, typename VectorType::value_type>::value,
3960 "Type mismatch between VectorType and VectorDataExchange");
3961 if (range_index == numbers::invalid_unsigned_int)
3962 vec = Number();
3963 else
3964 {
3965 const unsigned int mf_component = find_vector_in_mf(vec, false);
3966 const internal::MatrixFreeFunctions::DoFInfo &dof_info =
3967 matrix_free.get_dof_info(mf_component);
3968 Assert(dof_info.vector_zero_range_list_index.empty() == false,
3969 ExcNotInitialized());
3970
3971 Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3972 ExcInternalError());
3973 AssertIndexRange(range_index,
3974 dof_info.vector_zero_range_list_index.size() - 1);
3975 for (unsigned int id =
3976 dof_info.vector_zero_range_list_index[range_index];
3977 id != dof_info.vector_zero_range_list_index[range_index + 1];
3978 ++id)
3979 std::memset(vec.begin() + dof_info.vector_zero_range_list[id].first,
3980 0,
3981 (dof_info.vector_zero_range_list[id].second -
3982 dof_info.vector_zero_range_list[id].first) *
3983 sizeof(Number));
3984 }
3985 }
3986
3987
3988
3994 template <typename VectorType,
3995 typename std::enable_if<!has_exchange_on_subset<VectorType>,
3996 VectorType>::type * = nullptr,
3997 typename VectorType::value_type * = nullptr>
3998 void
3999 zero_vector_region(const unsigned int range_index, VectorType &vec) const
4000 {
4001 if (range_index == numbers::invalid_unsigned_int || range_index == 0)
4002 vec = typename VectorType::value_type();
4003 }
4004
4005
4006
4011 void
4012 zero_vector_region(const unsigned int, ...) const
4013 {
4014 Assert(false,
4015 ExcNotImplemented("Zeroing is only implemented for vector types "
4016 "which provide VectorType::value_type"));
4017 }
4018
4019
4020
4021 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
4022 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
4023 DataAccessOnFaces vector_face_access;
4024 bool ghosts_were_set;
4025# ifdef DEAL_II_WITH_MPI
4026 std::vector<AlignedVector<Number> *> tmp_data;
4027 std::vector<std::vector<MPI_Request>> requests;
4028# endif
4029 }; // VectorDataExchange
4030
4031 template <typename VectorStruct>
4032 unsigned int
4033 n_components(const VectorStruct &vec);
4034
4035 template <typename VectorStruct>
4036 unsigned int
4037 n_components_block(const VectorStruct &vec,
4038 std::integral_constant<bool, true>)
4039 {
4040 unsigned int components = 0;
4041 for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
4042 components += n_components(vec.block(bl));
4043 return components;
4044 }
4045
4046 template <typename VectorStruct>
4047 unsigned int
4048 n_components_block(const VectorStruct &, std::integral_constant<bool, false>)
4049 {
4050 return 1;
4051 }
4052
4053 template <typename VectorStruct>
4054 unsigned int
4055 n_components(const VectorStruct &vec)
4056 {
4057 return n_components_block(
4058 vec, std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
4059 }
4060
4061 template <typename VectorStruct>
4062 inline unsigned int
4063 n_components(const std::vector<VectorStruct> &vec)
4064 {
4065 unsigned int components = 0;
4066 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4067 components += n_components_block(
4068 vec[comp],
4069 std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
4070 return components;
4071 }
4072
4073 template <typename VectorStruct>
4074 inline unsigned int
4075 n_components(const std::vector<VectorStruct *> &vec)
4076 {
4077 unsigned int components = 0;
4078 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4079 components += n_components_block(
4080 *vec[comp],
4081 std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
4082 return components;
4083 }
4084
4085
4086
4087 // A helper function to identify block vectors with many components where we
4088 // should not try to overlap computations and communication because there
4089 // would be too many outstanding communication requests.
4090
4091 // default value for vectors that do not have communication_block_size
4092 template <typename VectorStruct,
4093 typename std::enable_if<!has_communication_block_size<VectorStruct>,
4094 VectorStruct>::type * = nullptr>
4095 constexpr unsigned int
4096 get_communication_block_size(const VectorStruct &)
4097 {
4098 return numbers::invalid_unsigned_int;
4099 }
4100
4101
4102
4103 template <typename VectorStruct,
4104 typename std::enable_if<has_communication_block_size<VectorStruct>,
4105 VectorStruct>::type * = nullptr>
4106 constexpr unsigned int
4107 get_communication_block_size(const VectorStruct &)
4108 {
4109 return VectorStruct::communication_block_size;
4110 }
4111
4112
4113
4114 // --------------------------------------------------------------------------
4115 // below we have wrappers to distinguish between block and non-block vectors.
4116 // --------------------------------------------------------------------------
4117
4118 //
4119 // update_ghost_values_start
4120 //
4121
4122 // update_ghost_values for block vectors
4123 template <int dim,
4124 typename VectorStruct,
4125 typename Number,
4126 typename VectorizedArrayType,
4127 typename std::enable_if<IsBlockVector<VectorStruct>::value,
4128 VectorStruct>::type * = nullptr>
4129 void
4130 update_ghost_values_start(
4131 const VectorStruct & vec,
4132 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4133 const unsigned int channel = 0)
4134 {
4135 if (get_communication_block_size(vec) < vec.n_blocks())
4136 {
4137 // don't forget to set ghosts_were_set, that otherwise happens
4138 // inside VectorDataExchange::update_ghost_values_start()
4139 exchanger.ghosts_were_set = vec.has_ghost_elements();
4140 vec.update_ghost_values();
4141 }
4142 else
4143 {
4144 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4145 update_ghost_values_start(vec.block(i), exchanger, channel + i);
4146 }
4147 }
4148
4149
4150
4151 // update_ghost_values for non-block vectors
4152 template <int dim,
4153 typename VectorStruct,
4154 typename Number,
4155 typename VectorizedArrayType,
4156 typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4157 VectorStruct>::type * = nullptr>
4158 void
4159 update_ghost_values_start(
4160 const VectorStruct & vec,
4161 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4162 const unsigned int channel = 0)
4163 {
4164 exchanger.update_ghost_values_start(channel, vec);
4165 }
4166
4167
4168
4169 // update_ghost_values_start() for vector of vectors
4170 template <int dim,
4171 typename VectorStruct,
4172 typename Number,
4173 typename VectorizedArrayType>
4174 inline void
4175 update_ghost_values_start(
4176 const std::vector<VectorStruct> & vec,
4177 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4178 {
4179 unsigned int component_index = 0;
4180 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4181 {
4182 update_ghost_values_start(vec[comp], exchanger, component_index);
4183 component_index += n_components(vec[comp]);
4184 }
4185 }
4186
4187
4188
4189 // update_ghost_values_start() for vector of pointers to vectors
4190 template <int dim,
4191 typename VectorStruct,
4192 typename Number,
4193 typename VectorizedArrayType>
4194 inline void
4195 update_ghost_values_start(
4196 const std::vector<VectorStruct *> & vec,
4197 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4198 {
4199 unsigned int component_index = 0;
4200 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4201 {
4202 update_ghost_values_start(*vec[comp], exchanger, component_index);
4203 component_index += n_components(*vec[comp]);
4204 }
4205 }
4206
4207
4208
4209 //
4210 // update_ghost_values_finish
4211 //
4212
4213 // for block vectors
4214 template <int dim,
4215 typename VectorStruct,
4216 typename Number,
4217 typename VectorizedArrayType,
4218 typename std::enable_if<IsBlockVector<VectorStruct>::value,
4219 VectorStruct>::type * = nullptr>
4220 void
4221 update_ghost_values_finish(
4222 const VectorStruct & vec,
4223 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4224 const unsigned int channel = 0)
4225 {
4226 if (get_communication_block_size(vec) < vec.n_blocks())
4227 {
4228 // do nothing, everything has already been completed in the _start()
4229 // call
4230 }
4231 else
4232 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4233 update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4234 }
4235
4236
4237
4238 // for non-block vectors
4239 template <int dim,
4240 typename VectorStruct,
4241 typename Number,
4242 typename VectorizedArrayType,
4243 typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4244 VectorStruct>::type * = nullptr>
4245 void
4246 update_ghost_values_finish(
4247 const VectorStruct & vec,
4248 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4249 const unsigned int channel = 0)
4250 {
4251 exchanger.update_ghost_values_finish(channel, vec);
4252 }
4253
4254
4255
4256 // for vector of vectors
4257 template <int dim,
4258 typename VectorStruct,
4259 typename Number,
4260 typename VectorizedArrayType>
4261 inline void
4262 update_ghost_values_finish(
4263 const std::vector<VectorStruct> & vec,
4264 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4265 {
4266 unsigned int component_index = 0;
4267 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4268 {
4269 update_ghost_values_finish(vec[comp], exchanger, component_index);
4270 component_index += n_components(vec[comp]);
4271 }
4272 }
4273
4274
4275
4276 // for vector of pointers to vectors
4277 template <int dim,
4278 typename VectorStruct,
4279 typename Number,
4280 typename VectorizedArrayType>
4281 inline void
4282 update_ghost_values_finish(
4283 const std::vector<VectorStruct *> & vec,
4284 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4285 {
4286 unsigned int component_index = 0;
4287 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4288 {
4289 update_ghost_values_finish(*vec[comp], exchanger, component_index);
4290 component_index += n_components(*vec[comp]);
4291 }
4292 }
4293
4294
4295
4296 //
4297 // compress_start
4298 //
4299
4300 // for block vectors
4301 template <int dim,
4302 typename VectorStruct,
4303 typename Number,
4304 typename VectorizedArrayType,
4305 typename std::enable_if<IsBlockVector<VectorStruct>::value,
4306 VectorStruct>::type * = nullptr>
4307 inline void
4308 compress_start(
4309 VectorStruct & vec,
4310 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4311 const unsigned int channel = 0)
4312 {
4313 if (get_communication_block_size(vec) < vec.n_blocks())
4314 vec.compress(::VectorOperation::add);
4315 else
4316 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4317 compress_start(vec.block(i), exchanger, channel + i);
4318 }
4319
4320
4321
4322 // for non-block vectors
4323 template <int dim,
4324 typename VectorStruct,
4325 typename Number,
4326 typename VectorizedArrayType,
4327 typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4328 VectorStruct>::type * = nullptr>
4329 inline void
4330 compress_start(
4331 VectorStruct & vec,
4332 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4333 const unsigned int channel = 0)
4334 {
4335 exchanger.compress_start(channel, vec);
4336 }
4337
4338
4339
4340 // for std::vector of vectors
4341 template <int dim,
4342 typename VectorStruct,
4343 typename Number,
4344 typename VectorizedArrayType>
4345 inline void
4346 compress_start(
4347 std::vector<VectorStruct> & vec,
4348 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4349 {
4350 unsigned int component_index = 0;
4351 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4352 {
4353 compress_start(vec[comp], exchanger, component_index);
4354 component_index += n_components(vec[comp]);
4355 }
4356 }
4357
4358
4359
4360 // for std::vector of pointer to vectors
4361 template <int dim,
4362 typename VectorStruct,
4363 typename Number,
4364 typename VectorizedArrayType>
4365 inline void
4366 compress_start(
4367 std::vector<VectorStruct *> & vec,
4368 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4369 {
4370 unsigned int component_index = 0;
4371 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4372 {
4373 compress_start(*vec[comp], exchanger, component_index);
4374 component_index += n_components(*vec[comp]);
4375 }
4376 }
4377
4378
4379
4380 //
4381 // compress_finish
4382 //
4383
4384 // for block vectors
4385 template <int dim,
4386 typename VectorStruct,
4387 typename Number,
4388 typename VectorizedArrayType,
4389 typename std::enable_if<IsBlockVector<VectorStruct>::value,
4390 VectorStruct>::type * = nullptr>
4391 inline void
4392 compress_finish(
4393 VectorStruct & vec,
4394 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4395 const unsigned int channel = 0)
4396 {
4397 if (get_communication_block_size(vec) < vec.n_blocks())
4398 {
4399 // do nothing, everything has already been completed in the _start()
4400 // call
4401 }
4402 else
4403 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4404 compress_finish(vec.block(i), exchanger, channel + i);
4405 }
4406
4407
4408
4409 // for non-block vectors
4410 template <int dim,
4411 typename VectorStruct,
4412 typename Number,
4413 typename VectorizedArrayType,
4414 typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4415 VectorStruct>::type * = nullptr>
4416 inline void
4417 compress_finish(
4418 VectorStruct & vec,
4419 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4420 const unsigned int channel = 0)
4421 {
4422 exchanger.compress_finish(channel, vec);
4423 }
4424
4425
4426
4427 // for std::vector of vectors
4428 template <int dim,
4429 typename VectorStruct,
4430 typename Number,
4431 typename VectorizedArrayType>
4432 inline void
4433 compress_finish(
4434 std::vector<VectorStruct> & vec,
4435 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4436 {
4437 unsigned int component_index = 0;
4438 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4439 {
4440 compress_finish(vec[comp], exchanger, component_index);
4441 component_index += n_components(vec[comp]);
4442 }
4443 }
4444
4445
4446
4447 // for std::vector of pointer to vectors
4448 template <int dim,
4449 typename VectorStruct,
4450 typename Number,
4451 typename VectorizedArrayType>
4452 inline void
4453 compress_finish(
4454 std::vector<VectorStruct *> & vec,
4455 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4456 {
4457 unsigned int component_index = 0;
4458 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4459 {
4460 compress_finish(*vec[comp], exchanger, component_index);
4461 component_index += n_components(*vec[comp]);
4462 }
4463 }
4464
4465
4466
4467 //
4468 // reset_ghost_values:
4469 //
4470 // if the input vector did not have ghosts imported, clear them here again
4471 // in order to avoid subsequent operations e.g. in linear solvers to work
4472 // with ghosts all the time
4473 //
4474
4475 // for block vectors
4476 template <int dim,
4477 typename VectorStruct,
4478 typename Number,
4479 typename VectorizedArrayType,
4480 typename std::enable_if<IsBlockVector<VectorStruct>::value,
4481 VectorStruct>::type * = nullptr>
4482 inline void
4483 reset_ghost_values(
4484 const VectorStruct & vec,
4485 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4486 {
4487 // return immediately if there is nothing to do.
4488 if (exchanger.ghosts_were_set == true)
4489 return;
4490
4491 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4492 reset_ghost_values(vec.block(i), exchanger);
4493 }
4494
4495
4496
4497 // for non-block vectors
4498 template <int dim,
4499 typename VectorStruct,
4500 typename Number,
4501 typename VectorizedArrayType,
4502 typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4503 VectorStruct>::type * = nullptr>
4504 inline void
4505 reset_ghost_values(
4506 const VectorStruct & vec,
4507 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4508 {
4509 exchanger.reset_ghost_values(vec);
4510 }
4511
4512
4513
4514 // for std::vector of vectors
4515 template <int dim,
4516 typename VectorStruct,
4517 typename Number,
4518 typename VectorizedArrayType>
4519 inline void
4520 reset_ghost_values(
4521 const std::vector<VectorStruct> & vec,
4522 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4523 {
4524 // return immediately if there is nothing to do.
4525 if (exchanger.ghosts_were_set == true)
4526 return;
4527
4528 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4529 reset_ghost_values(vec[comp], exchanger);
4530 }
4531
4532
4533
4534 // for std::vector of pointer to vectors
4535 template <int dim,
4536 typename VectorStruct,
4537 typename Number,
4538 typename VectorizedArrayType>
4539 inline void
4540 reset_ghost_values(
4541 const std::vector<VectorStruct *> & vec,
4542 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4543 {
4544 // return immediately if there is nothing to do.
4545 if (exchanger.ghosts_were_set == true)
4546 return;
4547
4548 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4549 reset_ghost_values(*vec[comp], exchanger);
4550 }
4551
4552
4553
4554 //
4555 // zero_vector_region
4556 //
4557
4558 // for block vectors
4559 template <int dim,
4560 typename VectorStruct,
4561 typename Number,
4562 typename VectorizedArrayType,
4563 typename std::enable_if<IsBlockVector<VectorStruct>::value,
4564 VectorStruct>::type * = nullptr>
4565 inline void
4566 zero_vector_region(
4567 const unsigned int range_index,
4568 VectorStruct & vec,
4569 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4570 {
4571 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4572 exchanger.zero_vector_region(range_index, vec.block(i));
4573 }
4574
4575
4576
4577 // for non-block vectors
4578 template <int dim,
4579 typename VectorStruct,
4580 typename Number,
4581 typename VectorizedArrayType,
4582 typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4583 VectorStruct>::type * = nullptr>
4584 inline void
4585 zero_vector_region(
4586 const unsigned int range_index,
4587 VectorStruct & vec,
4588 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4589 {
4590 exchanger.zero_vector_region(range_index, vec);
4591 }
4592
4593
4594
4595 // for std::vector of vectors
4596 template <int dim,
4597 typename VectorStruct,
4598 typename Number,
4599 typename VectorizedArrayType>
4600 inline void
4601 zero_vector_region(
4602 const unsigned int range_index,
4603 std::vector<VectorStruct> & vec,
4604 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4605 {
4606 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4607 zero_vector_region(range_index, vec[comp], exchanger);
4608 }
4609
4610
4611
4612 // for std::vector of pointers to vectors
4613 template <int dim,
4614 typename VectorStruct,
4615 typename Number,
4616 typename VectorizedArrayType>
4617 inline void
4618 zero_vector_region(
4619 const unsigned int range_index,
4620 std::vector<VectorStruct *> & vec,
4621 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4622 {
4623 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4624 zero_vector_region(range_index, *vec[comp], exchanger);
4625 }
4626
4627
4628
4629 // Apply a unit matrix operation to constrained DoFs: Default cases where we
4630 // cannot detect a LinearAlgebra::distributed::Vector, we do not do
4631 // anything, else we apply the constraints as a unit operation
4632 template <typename VectorStruct1, typename VectorStruct2>
4633 inline void
4634 apply_operation_to_constrained_dofs(const std::vector<unsigned int> &,
4635 const VectorStruct1 &,
4636 VectorStruct2 &)
4637 {}
4638
4639 template <typename Number>
4640 inline void
4641 apply_operation_to_constrained_dofs(
4642 const std::vector<unsigned int> & constrained_dofs,
4643 const LinearAlgebra::distributed::Vector<Number> &src,
4644 LinearAlgebra::distributed::Vector<Number> & dst)
4645 {
4646 for (const unsigned int i : constrained_dofs)
4647 dst.local_element(i) = src.local_element(i);
4648 }
4649
4650
4651 namespace MatrixFreeFunctions
4652 {
4653 // struct to select between a const interface and a non-const interface
4654 // for MFWorker
4655 template <typename, typename, typename, typename, bool>
4656 struct InterfaceSelector
4657 {};
4658
4659 // Version for constant functions
4660 template <typename MF,
4661 typename InVector,
4662 typename OutVector,
4663 typename Container>
4664 struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4665 {
4666 using function_type = void (Container::*)(
4667 const MF &,
4668 OutVector &,
4669 const InVector &,
4670 const std::pair<unsigned int, unsigned int> &) const;
4671 };
4672
4673 // Version for non-constant functions
4674 template <typename MF,
4675 typename InVector,
4676 typename OutVector,
4677 typename Container>
4678 struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4679 {
4680 using function_type =
4681 void (Container::*)(const MF &,
4682 OutVector &,
4683 const InVector &,
4684 const std::pair<unsigned int, unsigned int> &);
4685 };
4686 } // namespace MatrixFreeFunctions
4687
4688
4689
4690 // A implementation class for the worker object that runs the various
4691 // operations we want to perform during the matrix-free loop
4692 template <typename MF,
4693 typename InVector,
4694 typename OutVector,
4695 typename Container,
4696 bool is_constant>
4697 class MFWorker : public MFWorkerInterface
4698 {
4699 public:
4700 // An alias to make the arguments further down more readable
4701 using function_type = typename MatrixFreeFunctions::
4702 InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4703 function_type;
4704
4705 // constructor, binds all the arguments to this class
4706 MFWorker(const MF & matrix_free,
4707 const InVector & src,
4708 OutVector & dst,
4709 const bool zero_dst_vector_setting,
4710 const Container & container,
4711 function_type cell_function,
4712 function_type face_function,
4713 function_type boundary_function,
4714 const typename MF::DataAccessOnFaces src_vector_face_access =
4715 MF::DataAccessOnFaces::none,
4716 const typename MF::DataAccessOnFaces dst_vector_face_access =
4717 MF::DataAccessOnFaces::none,
4718 const std::function<void(const unsigned int, const unsigned int)>
4719 &operation_before_loop = {},
4720 const std::function<void(const unsigned int, const unsigned int)>
4721 & operation_after_loop = {},
4722 const unsigned int dof_handler_index_pre_post = 0)
4723 : matrix_free(matrix_free)
4724 , container(const_cast<Container &>(container))
4725 , cell_function(cell_function)
4726 , face_function(face_function)
4727 , boundary_function(boundary_function)
4728 , src(src)
4729 , dst(dst)
4730 , src_data_exchanger(matrix_free,
4731 src_vector_face_access,
4732 n_components(src))
4733 , dst_data_exchanger(matrix_free,
4734 dst_vector_face_access,
4735 n_components(dst))
4736 , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4737 , zero_dst_vector_setting(zero_dst_vector_setting &&
4738 !src_and_dst_are_same)
4739 , operation_before_loop(operation_before_loop)
4740 , operation_after_loop(operation_after_loop)
4741 , dof_handler_index_pre_post(dof_handler_index_pre_post)
4742 {}
4743
4744 // Runs the cell work. If no function is given, nothing is done
4745 virtual void
4746 cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4747 {
4748 if (cell_function != nullptr && cell_range.second > cell_range.first)
4749 for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4750 {
4751 const auto cell_subrange =
4752 matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4753
4754 if (cell_subrange.second <= cell_subrange.first)
4755 continue;
4756
4757 (container.*
4758 cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4759 }
4760 }
4761
4762 virtual void
4763 cell(const unsigned int range_index) override
4764 {
4765 process_range(cell_function,
4766 matrix_free.get_task_info().cell_partition_data_hp_ptr,
4767 matrix_free.get_task_info().cell_partition_data_hp,
4768 range_index);
4769 }
4770
4771 virtual void
4772 face(const unsigned int range_index) override
4773 {
4774 process_range(face_function,
4775 matrix_free.get_task_info().face_partition_data_hp_ptr,
4776 matrix_free.get_task_info().face_partition_data_hp,
4777 range_index);
4778 }
4779
4780 virtual void
4781 boundary(const unsigned int range_index) override
4782 {
4783 process_range(boundary_function,
4784 matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4785 matrix_free.get_task_info().boundary_partition_data_hp,
4786 range_index);
4787 }
4788
4789 private:
4790 void
4791 process_range(const function_type & fu,
4792 const std::vector<unsigned int> &ptr,
4793 const std::vector<unsigned int> &data,
4794 const unsigned int range_index)
4795 {
4796 if (fu == nullptr)
4797 return;
4798
4799 AssertIndexRange(range_index + 1, ptr.size());
4800 for (unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4801 {
4802 AssertIndexRange(2 * i + 1, data.size());
4803 (container.*fu)(matrix_free,
4804 this->dst,
4805 this->src,
4806 std::make_pair(data[2 * i], data[2 * i + 1]));
4807 }
4808 }
4809
4810 public:
4811 // Starts the communication for the update ghost values operation. We
4812 // cannot call this update if ghost and destination are the same because
4813 // that would introduce spurious entries in the destination (there is also
4814 // the problem that reading from a vector that we also write to is usually
4815 // not intended in case there is overlap, but this is up to the
4816 // application code to decide and we cannot catch this case here).
4817 virtual void
4818 vector_update_ghosts_start() override
4819 {
4820 if (!src_and_dst_are_same)
4821 internal::update_ghost_values_start(src, src_data_exchanger);
4822 }
4823
4824 // Finishes the communication for the update ghost values operation
4825 virtual void
4826 vector_update_ghosts_finish() override
4827 {
4828 if (!src_and_dst_are_same)
4829 internal::update_ghost_values_finish(src, src_data_exchanger);
4830 }
4831
4832 // Starts the communication for the vector compress operation
4833 virtual void
4834 vector_compress_start() override
4835 {
4836 internal::compress_start(dst, dst_data_exchanger);
4837 }
4838
4839 // Finishes the communication for the vector compress operation
4840 virtual void
4841 vector_compress_finish() override
4842 {
4843 internal::compress_finish(dst, dst_data_exchanger);
4844 if (!src_and_dst_are_same)
4845 internal::reset_ghost_values(src, src_data_exchanger);
4846 }
4847
4848 // Zeros the given input vector
4849 virtual void
4850 zero_dst_vector_range(const unsigned int range_index) override
4851 {
4852 if (zero_dst_vector_setting)
4853 internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4854 }
4855
4856 virtual void
4857 cell_loop_pre_range(const unsigned int range_index) override
4858 {
4859 if (operation_before_loop)
4860 {
4861 const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4862 matrix_free.get_dof_info(dof_handler_index_pre_post);
4863 if (range_index == numbers::invalid_unsigned_int)
4864 {
4865 // Case with threaded loop -> currently no overlap implemented
4866 ::parallel::apply_to_subranges(
4867 0U,
4868 dof_info.vector_partitioner->locally_owned_size(),
4869 operation_before_loop,
4870 internal::VectorImplementation::minimum_parallel_grain_size);
4871 }
4872 else
4873 {
4874 AssertIndexRange(range_index,
4875 dof_info.cell_loop_pre_list_index.size() - 1);
4876 for (unsigned int id =
4877 dof_info.cell_loop_pre_list_index[range_index];
4878 id != dof_info.cell_loop_pre_list_index[range_index + 1];
4879 ++id)
4880 operation_before_loop(dof_info.cell_loop_pre_list[id].first,
4881 dof_info.cell_loop_pre_list[id].second);
4882 }
4883 }
4884 }
4885
4886 virtual void
4887 cell_loop_post_range(const unsigned int range_index) override
4888 {
4889 if (operation_after_loop)
4890 {
4891 // Run unit matrix operation on constrained dofs if we are at the
4892 // last range
4893 const std::vector<unsigned int> &partition_row_index =
4894 matrix_free.get_task_info().partition_row_index;
4895 if (range_index ==
4896 partition_row_index[partition_row_index.size() - 2] - 1)
4897 apply_operation_to_constrained_dofs(
4898 matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
4899 src,
4900 dst);
4901
4902 const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4903 matrix_free.get_dof_info(dof_handler_index_pre_post);
4904 if (range_index == numbers::invalid_unsigned_int)
4905 {
4906 // Case with threaded loop -> currently no overlap implemented
4907 ::parallel::apply_to_subranges(
4908 0U,
4909 dof_info.vector_partitioner->locally_owned_size(),
4910 operation_after_loop,
4911 internal::VectorImplementation::minimum_parallel_grain_size);
4912 }
4913 else
4914 {
4915 AssertIndexRange(range_index,
4916 dof_info.cell_loop_post_list_index.size() - 1);
4917 for (unsigned int id =
4918 dof_info.cell_loop_post_list_index[range_index];
4919 id != dof_info.cell_loop_post_list_index[range_index + 1];
4920 ++id)
4921 operation_after_loop(dof_info.cell_loop_post_list[id].first,
4922 dof_info.cell_loop_post_list[id].second);
4923 }
4924 }
4925 }
4926
4927 private:
4928 const MF & matrix_free;
4929 Container & container;
4930 function_type cell_function;
4931 function_type face_function;
4932 function_type boundary_function;
4933
4934 const InVector &src;
4935 OutVector & dst;
4936 VectorDataExchange<MF::dimension,
4937 typename MF::value_type,
4938 typename MF::vectorized_value_type>
4939 src_data_exchanger;
4940 VectorDataExchange<MF::dimension,
4941 typename MF::value_type,
4942 typename MF::vectorized_value_type>
4943 dst_data_exchanger;
4944 const bool src_and_dst_are_same;
4945 const bool zero_dst_vector_setting;
4946 const std::function<void(const unsigned int, const unsigned int)>
4947 operation_before_loop;
4948 const std::function<void(const unsigned int, const unsigned int)>
4949 operation_after_loop;
4950 const unsigned int dof_handler_index_pre_post;
4951 };
4952
4953
4954
4959 template <class MF, typename InVector, typename OutVector>
4960 struct MFClassWrapper
4961 {
4962 using function_type =
4963 std::function<void(const MF &,
4964 OutVector &,
4965 const InVector &,
4966 const std::pair<unsigned int, unsigned int> &)>;
4967
4968 MFClassWrapper(const function_type cell,
4969 const function_type face,
4970 const function_type boundary)
4971 : cell(cell)
4972 , face(face)
4973 , boundary(boundary)
4974 {}
4975
4976 void
4977 cell_integrator(const MF & mf,
4978 OutVector & dst,
4979 const InVector & src,
4980 const std::pair<unsigned int, unsigned int> &range) const
4981 {
4982 if (cell)
4983 cell(mf, dst, src, range);
4984 }
4985
4986 void
4987 face_integrator(const MF & mf,
4988 OutVector & dst,
4989 const InVector & src,
4990 const std::pair<unsigned int, unsigned int> &range) const
4991 {
4992 if (face)
4993 face(mf, dst, src, range);
4994 }
4995
4996 void
4997 boundary_integrator(
4998 const MF & mf,
4999 OutVector & dst,
5000 const InVector & src,
5001 const std::pair<unsigned int, unsigned int> &range) const
5002 {
5003 if (boundary)
5004 boundary(mf, dst, src, range);
5005 }
5006
5007 const function_type cell;
5008 const function_type face;
5009 const function_type boundary;
5010 };
5011
5012} // end of namespace internal
5013
5014
5015
5016template <int dim, typename Number, typename VectorizedArrayType>
5017template <typename OutVector, typename InVector>
5018inline void
5019MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5020 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5021 OutVector &,
5022 const InVector &,
5023 const std::pair<unsigned int, unsigned int> &)>
5024 & cell_operation,
5025 OutVector & dst,
5026 const InVector &src,
5027 const bool zero_dst_vector) const
5028{
5029 using Wrapper =
5030 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5031 InVector,
5032 OutVector>;
5033 Wrapper wrap(cell_operation, nullptr, nullptr);
5034 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5035 InVector,
5036 OutVector,
5037 Wrapper,
5038 true>
5039 worker(*this,
5040 src,
5041 dst,
5042 zero_dst_vector,
5043 wrap,
5044 &Wrapper::cell_integrator,
5045 &Wrapper::face_integrator,
5046 &Wrapper::boundary_integrator);
5047
5048 task_info.loop(worker);
5049}
5050
5051
5052
5053template <int dim, typename Number, typename VectorizedArrayType>
5054template <typename OutVector, typename InVector>
5055inline void
5056MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5057 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5058 OutVector &,
5059 const InVector &,
5060 const std::pair<unsigned int, unsigned int> &)>
5061 & cell_operation,
5062 OutVector & dst,
5063 const InVector &src,
5064 const std::function<void(const unsigned int, const unsigned int)>
5065 &operation_before_loop,
5066 const std::function<void(const unsigned int, const unsigned int)>
5067 & operation_after_loop,
5068 const unsigned int dof_handler_index_pre_post) const
5069{
5070 using Wrapper =
5071 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5072 InVector,
5073 OutVector>;
5074 Wrapper wrap(cell_operation, nullptr, nullptr);
5075 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5076 InVector,
5077 OutVector,
5078 Wrapper,
5079 true>
5080 worker(*this,
5081 src,
5082 dst,
5083 false,
5084 wrap,
5085 &Wrapper::cell_integrator,
5086 &Wrapper::face_integrator,
5087 &Wrapper::boundary_integrator,
5088 DataAccessOnFaces::none,
5089 DataAccessOnFaces::none,
5090 operation_before_loop,
5091 operation_after_loop,
5092 dof_handler_index_pre_post);
5093
5094 task_info.loop(worker);
5095}
5096
5097
5098
5099template <int dim, typename Number, typename VectorizedArrayType>
5100template <typename OutVector, typename InVector>
5101inline void
5102MatrixFree<dim, Number, VectorizedArrayType>::loop(
5103 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5104 OutVector &,
5105 const InVector &,
5106 const std::pair<unsigned int, unsigned int> &)>
5107 &cell_operation,
5108 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5109 OutVector &,
5110 const InVector &,
5111 const std::pair<unsigned int, unsigned int> &)>
5112 &face_operation,
5113 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5114 OutVector &,
5115 const InVector &,
5116 const std::pair<unsigned int, unsigned int> &)>
5117 & boundary_operation,
5118 OutVector & dst,
5119 const InVector & src,
5120 const bool zero_dst_vector,
5121 const DataAccessOnFaces dst_vector_face_access,
5122 const DataAccessOnFaces src_vector_face_access) const
5123{
5124 using Wrapper =
5125 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5126 InVector,
5127 OutVector>;
5128 Wrapper wrap(cell_operation, face_operation, boundary_operation);
5129 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5130 InVector,
5131 OutVector,
5132 Wrapper,
5133 true>
5134 worker(*this,
5135 src,
5136 dst,
5137 zero_dst_vector,
5138 wrap,
5139 &Wrapper::cell_integrator,
5140 &Wrapper::face_integrator,
5141 &Wrapper::boundary_integrator,
5142 src_vector_face_access,
5143 dst_vector_face_access);
5144
5145 task_info.loop(worker);
5146}
5147
5148
5149
5150template <int dim, typename Number, typename VectorizedArrayType>
5151template <typename CLASS, typename OutVector, typename InVector>
5152inline void
5153MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5154 void (CLASS::*function_pointer)(
5155 const MatrixFree<dim, Number, VectorizedArrayType> &,
5156 OutVector &,
5157 const InVector &,
5158 const std::pair<unsigned int, unsigned int> &) const,
5159 const CLASS * owning_class,
5160 OutVector & dst,
5161 const InVector &src,
5162 const bool zero_dst_vector) const
5163{
5164 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5165 InVector,
5166 OutVector,
5167 CLASS,
5168 true>
5169 worker(*this,
5170 src,
5171 dst,
5172 zero_dst_vector,
5173 *owning_class,
5174 function_pointer,
5175 nullptr,
5176 nullptr);
5177 task_info.loop(worker);
5178}
5179
5180
5181
5182template <int dim, typename Number, typename VectorizedArrayType>
5183template <typename CLASS, typename OutVector, typename InVector>
5184inline void
5185MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5186 void (CLASS::*function_pointer)(
5187 const MatrixFree<dim, Number, VectorizedArrayType> &,
5188 OutVector &,
5189 const InVector &,
5190 const std::pair<unsigned int, unsigned int> &) const,
5191 const CLASS * owning_class,
5192 OutVector & dst,
5193 const InVector &src,
5194 const std::function<void(const unsigned int, const unsigned int)>
5195 &operation_before_loop,
5196 const std::function<void(const unsigned int, const unsigned int)>
5197 & operation_after_loop,
5198 const unsigned int dof_handler_index_pre_post) const
5199{
5200 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5201 InVector,
5202 OutVector,
5203 CLASS,
5204 true>
5205 worker(*this,
5206 src,
5207 dst,
5208 false,
5209 *owning_class,
5210 function_pointer,
5211 nullptr,
5212 nullptr,
5213 DataAccessOnFaces::none,
5214 DataAccessOnFaces::none,
5215 operation_before_loop,
5216 operation_after_loop,
5217 dof_handler_index_pre_post);
5218 task_info.loop(worker);
5219}
5220
5221
5222
5223template <int dim, typename Number, typename VectorizedArrayType>
5224template <typename CLASS, typename OutVector, typename InVector>
5225inline void
5226MatrixFree<dim, Number, VectorizedArrayType>::loop(
5227 void (CLASS::*cell_operation)(
5228 const MatrixFree<dim, Number, VectorizedArrayType> &,
5229 OutVector &,
5230 const InVector &,
5231 const std::pair<unsigned int, unsigned int> &) const,
5232 void (CLASS::*face_operation)(
5233 const MatrixFree<dim, Number, VectorizedArrayType> &,
5234 OutVector &,
5235 const InVector &,
5236 const std::pair<unsigned int, unsigned int> &) const,
5237 void (CLASS::*boundary_operation)(
5238 const MatrixFree<dim, Number, VectorizedArrayType> &,
5239 OutVector &,
5240 const InVector &,
5241 const std::pair<unsigned int, unsigned int> &) const,
5242 const CLASS * owning_class,
5243 OutVector & dst,
5244 const InVector & src,
5245 const bool zero_dst_vector,
5246 const DataAccessOnFaces dst_vector_face_access,
5247 const DataAccessOnFaces src_vector_face_access) const
5248{
5249 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5250 InVector,
5251 OutVector,
5252 CLASS,
5253 true>
5254 worker(*this,
5255 src,
5256 dst,
5257 zero_dst_vector,
5258 *owning_class,
5259 cell_operation,
5260 face_operation,
5261 boundary_operation,
5262 src_vector_face_access,
5263 dst_vector_face_access);
5264 task_info.loop(worker);
5265}
5266
5267
5268
5269template <int dim, typename Number, typename VectorizedArrayType>
5270template <typename CLASS, typename OutVector, typename InVector>
5271inline void
5272MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5273 void (CLASS::*function_pointer)(
5274 const MatrixFree<dim, Number, VectorizedArrayType> &,
5275 OutVector &,
5276 const InVector &,
5277 const std::pair<unsigned int, unsigned int> &),
5278 CLASS * owning_class,
5279 OutVector & dst,
5280 const InVector &src,
5281 const bool zero_dst_vector) const
5282{
5283 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5284 InVector,
5285 OutVector,
5286 CLASS,
5287 false>
5288 worker(*this,
5289 src,
5290 dst,
5291 zero_dst_vector,
5292 *owning_class,
5293 function_pointer,
5294 nullptr,
5295 nullptr);
5296 task_info.loop(worker);
5297}
5298
5299
5300
5301template <int dim, typename Number, typename VectorizedArrayType>
5302template <typename CLASS, typename OutVector, typename InVector>
5303inline void
5304MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5305 void (CLASS::*function_pointer)(
5306 const MatrixFree<dim, Number, VectorizedArrayType> &,
5307 OutVector &,
5308 const InVector &,
5309 const std::pair<unsigned int, unsigned int> &),
5310 CLASS * owning_class,
5311 OutVector & dst,
5312 const InVector &src,
5313 const std::function<void(const unsigned int, const unsigned int)>
5314 &operation_before_loop,
5315 const std::function<void(const unsigned int, const unsigned int)>
5316 & operation_after_loop,
5317 const unsigned int dof_handler_index_pre_post) const
5318{
5319 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5320 InVector,
5321 OutVector,
5322 CLASS,
5323 false>
5324 worker(*this,
5325 src,
5326 dst,
5327 false,
5328 *owning_class,
5329 function_pointer,
5330 nullptr,
5331 nullptr,
5332 DataAccessOnFaces::none,
5333 DataAccessOnFaces::none,
5334 operation_before_loop,
5335 operation_after_loop,
5336 dof_handler_index_pre_post);
5337 task_info.loop(worker);
5338}
5339
5340
5341
5342template <int dim, typename Number, typename VectorizedArrayType>
5343template <typename CLASS, typename OutVector, typename InVector>
5344inline void
5345MatrixFree<dim, Number, VectorizedArrayType>::loop(
5346 void (CLASS::*cell_operation)(
5347 const MatrixFree<dim, Number, VectorizedArrayType> &,
5348 OutVector &,
5349 const InVector &,
5350 const std::pair<unsigned int, unsigned int> &),
5351 void (CLASS::*face_operation)(
5352 const MatrixFree<dim, Number, VectorizedArrayType> &,
5353 OutVector &,
5354 const InVector &,
5355 const std::pair<unsigned int, unsigned int> &),
5356 void (CLASS::*boundary_operation)(
5357 const MatrixFree<dim, Number, VectorizedArrayType> &,
5358 OutVector &,
5359 const InVector &,
5360 const std::pair<unsigned int, unsigned int> &),
5361 CLASS * owning_class,
5362 OutVector & dst,
5363 const InVector & src,
5364 const bool zero_dst_vector,
5365 const DataAccessOnFaces dst_vector_face_access,
5366 const DataAccessOnFaces src_vector_face_access) const
5367{
5368 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5369 InVector,
5370 OutVector,
5371 CLASS,
5372 false>
5373 worker(*this,
5374 src,
5375 dst,
5376 zero_dst_vector,
5377 *owning_class,
5378 cell_operation,
5379 face_operation,
5380 boundary_operation,
5381 src_vector_face_access,
5382 dst_vector_face_access);
5383 task_info.loop(worker);
5384}
5385
5386
5387
5388template <int dim, typename Number, typename VectorizedArrayType>
5389template <typename CLASS, typename OutVector, typename InVector>
5390inline void
5391MatrixFree<dim, Number, VectorizedArrayType>::loop_cell_centric(
5392 void (CLASS::*function_pointer)(
5393 const MatrixFree<dim, Number, VectorizedArrayType> &,
5394 OutVector &,
5395 const InVector &,
5396 const std::pair<unsigned int, unsigned int> &) const,
5397 const CLASS * owning_class,
5398 OutVector & dst,
5399 const InVector & src,
5400 const bool zero_dst_vector,
5401 const DataAccessOnFaces src_vector_face_access) const
5402{
5403 auto src_vector_face_access_temp = src_vector_face_access;
5404 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5405 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5406 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5407 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5408
5409 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5410 InVector,
5411 OutVector,
5412 CLASS,
5413 true>
5414 worker(*this,
5415 src,
5416 dst,
5417 zero_dst_vector,
5418 *owning_class,
5419 function_pointer,
5420 nullptr,
5421 nullptr,
5422 src_vector_face_access_temp,
5423 DataAccessOnFaces::none);
5424 task_info.loop(worker);
5425}
5426
5427
5428
5429template <int dim, typename Number, typename VectorizedArrayType>
5430template <typename CLASS, typename OutVector, typename InVector>
5431inline void
5432MatrixFree<dim, Number, VectorizedArrayType>::loop_cell_centric(
5433 void (CLASS::*function_pointer)(
5434 const MatrixFree<dim, Number, VectorizedArrayType> &,
5435 OutVector &,
5436 const InVector &,
5437 const std::pair<unsigned int, unsigned int> &),
5438 CLASS * owning_class,
5439 OutVector & dst,
5440 const InVector & src,
5441 const bool zero_dst_vector,
5442 const DataAccessOnFaces src_vector_face_access) const
5443{
5444 auto src_vector_face_access_temp = src_vector_face_access;
5445 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5446 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5447 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5448 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5449
5450 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5451 InVector,
5452 OutVector,
5453 CLASS,
5454 false>
5455 worker(*this,
5456 src,
5457 dst,
5458 zero_dst_vector,
5459 *owning_class,
5460 function_pointer,
5461 nullptr,
5462 nullptr,
5463 src_vector_face_access_temp,
5464 DataAccessOnFaces::none);
5465 task_info.loop(worker);
5466}
5467
5468
5469
5470template <int dim, typename Number, typename VectorizedArrayType>
5471template <typename OutVector, typename InVector>
5472inline void
5473MatrixFree<dim, Number, VectorizedArrayType>::loop_cell_centric(
5474 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5475 OutVector &,
5476 const InVector &,
5477 const std::pair<unsigned int, unsigned int> &)>
5478 & cell_operation,
5479 OutVector & dst,
5480 const InVector & src,
5481 const bool zero_dst_vector,
5482 const DataAccessOnFaces src_vector_face_access) const
5483{
5484 auto src_vector_face_access_temp = src_vector_face_access;
5485 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5486 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5487 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5488 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5489
5490 using Wrapper =
5491 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5492 InVector,
5493 OutVector>;
5494 Wrapper wrap(cell_operation, nullptr, nullptr);
5495
5496 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5497 InVector,
5498 OutVector,
5499 Wrapper,
5500 true>
5501 worker(*this,
5502 src,
5503 dst,
5504 zero_dst_vector,
5505 wrap,
5506 &Wrapper::cell_integrator,
5507 &Wrapper::face_integrator,
5508 &Wrapper::boundary_integrator,
5509 src_vector_face_access_temp,
5510 DataAccessOnFaces::none);
5511 task_info.loop(worker);
5512}
5513
5514
5515#endif // ifndef DOXYGEN
5516
5517
5518
5519DEAL_II_NAMESPACE_CLOSE
5520
5521#endif
affine_constraints.h
aligned_vector.h
block_vector_base.h
AlignedVector::begin
iterator begin()
AlignedVector::resize
void resize(const size_type new_size)
MPI_Comm
Mapping
Abstract base class for mapping classes.
Definition: mapping.h:311
MatrixFree::n_active_fe_indices
unsigned int n_active_fe_indices() const
MatrixFree::get_face_iterator
std::pair< typename DoFHandler< dim >::cell_iterator, unsigned int > get_face_iterator(const unsigned int face_batch_index, const unsigned int lane_index, const bool interior=true, const unsigned int fe_component=0) const
MatrixFree::loop_cell_centric
void loop_cell_centric(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
MatrixFree::mapping_info
internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > mapping_info
Definition: matrix_free.h:2157
MatrixFree::get_cell_active_fe_index
unsigned int get_cell_active_fe_index(const std::pair< unsigned int, unsigned int > range) const
MatrixFree::get_task_info
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
MatrixFree::get_hp_cell_iterator
DoFHandler< dim >::active_cell_iterator get_hp_cell_iterator(const unsigned int cell_batch_index, const unsigned int lane_index, const unsigned int dof_handler_index=0) const
MatrixFree::n_ghost_cell_batches
unsigned int n_ghost_cell_batches() const
MatrixFree::loop
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
MatrixFree::get_shape_info
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
MatrixFree::initialize_dof_handlers
void initialize_dof_handlers(const std::vector< const DoFHandler< dim, dim > * > &dof_handlers, const AdditionalData &additional_data)
MatrixFree::get_boundary_id
types::boundary_id get_boundary_id(const unsigned int face_batch_index) const
MatrixFree::shape_info
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > > shape_info
Definition: matrix_free.h:2163
MatrixFree::print
void print(std::ostream &out) const
MatrixFree::update_mapping
void update_mapping(const std::shared_ptr< hp::MappingCollection< dim > > &mapping)
MatrixFree::get_cell_and_face_to_plain_faces
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
MatrixFree::reinit
void reinit(const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
MatrixFree::n_inner_face_batches
unsigned int n_inner_face_batches() const
MatrixFree::n_active_entries_per_cell_batch
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
MatrixFree::update_mapping
void update_mapping(const Mapping< dim > &mapping)
MatrixFree::get_quadrature
const Quadrature< dim > & get_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
MatrixFree::get_face_info
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_index) const
MatrixFree::get_mg_level
unsigned int get_mg_level() const
MatrixFree::n_components_filled
unsigned int n_components_filled(const unsigned int cell_batch_number) const
MatrixFree::print_memory_consumption
void print_memory_consumption(StreamType &out) const
MatrixFree::reinit
void reinit(const std::vector< const DoFHandlerType * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< QuadratureType > &quad, const AdditionalData &additional_data=AdditionalData())
MatrixFree::mapping_initialized
bool mapping_initialized() const
MatrixFree::cell_loop
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
MatrixFree::clear
void clear()
MatrixFree::internal_reinit
void internal_reinit(const std::shared_ptr< hp::MappingCollection< dim > > &mapping, const std::vector< const DoFHandler< dim, dim > * > &dof_handlers, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const std::vector< hp::QCollection< q_dim > > &quad, const AdditionalData &additional_data)
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
MatrixFree::get_dof_info
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
MatrixFree::acquire_scratch_data
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
MatrixFree::at_irregular_cell
bool at_irregular_cell(const unsigned int cell_batch_index) const
MatrixFree::~MatrixFree
~MatrixFree() override=default
MatrixFree::get_dof_handler
const DoFHandlerType & get_dof_handler(const unsigned int dof_handler_index=0) const
MatrixFree::copy_from
void copy_from(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free_base)
MatrixFree::initialize_face_data_vector
void initialize_face_data_vector(AlignedVector< T > &vec) 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_constraint_pool_entries
unsigned int n_constraint_pool_entries() const
MatrixFree::cell_loop
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
MatrixFree::constraint_pool_begin
const Number * constraint_pool_begin(const unsigned int pool_index) const
MatrixFree::acquire_scratch_data_non_threadsafe
AlignedVector< Number > * acquire_scratch_data_non_threadsafe() const
MatrixFree::get_locally_owned_set
const IndexSet & get_locally_owned_set(const unsigned int dof_handler_index=0) const
MatrixFree::get_dofs_per_face
unsigned int get_dofs_per_face(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
MatrixFree::task_info
internal::MatrixFreeFunctions::TaskInfo task_info
Definition: matrix_free.h:2187
MatrixFree::get_dof_handler
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
MatrixFree::release_scratch_data
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
MatrixFree::loop
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*boundary_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
MatrixFree::mapping_is_initialized
bool mapping_is_initialized
Definition: matrix_free.h:2204
MatrixFree::loop_cell_centric
void loop_cell_centric(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
MatrixFree::initialize_dof_vector
void initialize_dof_vector(LinearAlgebra::distributed::Vector< Number2 > &vec, const unsigned int dof_handler_index=0) const
MatrixFree::get_cell_iterator
DoFHandler< dim >::cell_iterator get_cell_iterator(const unsigned int cell_batch_index, const unsigned int lane_index, const unsigned int dof_handler_index=0) const
MatrixFree::MatrixFree
MatrixFree(const MatrixFree< dim, Number, VectorizedArrayType > &other)
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::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::mg_level
unsigned int mg_level
Definition: matrix_free.h:2227
MatrixFree::get_face_range_category
std::pair< unsigned int, unsigned int > get_face_range_category(const std::pair< unsigned int, unsigned int > face_batch_range) const
MatrixFree::cell_loop
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
MatrixFree::cell_level_index
std::vector< std::pair< unsigned int, unsigned int > > cell_level_index
Definition: matrix_free.h:2171
MatrixFree::get_vector_partitioner
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner(const unsigned int dof_handler_index=0) const
MatrixFree::get_mapping_info
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
MatrixFree::reinit
void reinit(const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< QuadratureType > &quad, const AdditionalData &additional_data=AdditionalData())
MatrixFree::constraint_pool_data
std::vector< Number > constraint_pool_data
Definition: matrix_free.h:2144
MatrixFree::vectorized_value_type
VectorizedArrayType vectorized_value_type
Definition: matrix_free.h:126
MatrixFree::n_cell_batches
unsigned int n_cell_batches() const
MatrixFree::get_ghost_set
const IndexSet & get_ghost_set(const unsigned int dof_handler_index=0) const
MatrixFree::get_cell_level_and_index
std::pair< int, int > get_cell_level_and_index(const unsigned int cell_batch_index, const unsigned int lane_index) const
MatrixFree::indices_initialized
bool indices_initialized() const
MatrixFree::get_face_quadrature
const Quadrature< dim - 1 > & get_face_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
MatrixFree::get_cell_category
unsigned int get_cell_category(const unsigned int cell_batch_index) const
MatrixFree::reinit
void reinit(const MappingType &mapping, const std::vector< const DoFHandlerType * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
MatrixFree::scratch_pad
Threads::ThreadLocalStorage< std::list< std::pair< bool, AlignedVector< VectorizedArrayType > > > > scratch_pad
Definition: matrix_free.h:2215
MatrixFree::value_type
Number value_type
Definition: matrix_free.h:125
MatrixFree::memory_consumption
std::size_t memory_consumption() const
MatrixFree::get_face_category
std::pair< unsigned int, unsigned int > get_face_category(const unsigned int face_batch_index) const
MatrixFree::n_boundary_face_batches
unsigned int n_boundary_face_batches() const
MatrixFree::cell_loop
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
MatrixFree::scratch_pad_non_threadsafe
std::list< std::pair< bool, AlignedVector< Number > > > scratch_pad_non_threadsafe
Definition: matrix_free.h:2222
MatrixFree::constraint_pool_end
const Number * constraint_pool_end(const unsigned int pool_index) const
MatrixFree::loop_cell_centric
void loop_cell_centric(const std::function< void(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
MatrixFree::dof_handlers
std::vector< SmartPointer< const DoFHandler< dim > > > dof_handlers
Definition: matrix_free.h:2130
MatrixFree::reinit
void reinit(const std::vector< const DoFHandlerType * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
MatrixFree::initialize_indices
void initialize_indices(const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const AdditionalData &additional_data)
MatrixFree::n_components
unsigned int n_components() const
MatrixFree::n_physical_cells
unsigned int n_physical_cells() const
MatrixFree::initialize_dof_vector
void initialize_dof_vector(VectorType &vec, const unsigned int dof_handler_index=0) 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
MatrixFree::reinit
void reinit(const MappingType &mapping, const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
MatrixFree::cell_loop
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
MatrixFree::n_macro_cells
unsigned int n_macro_cells() const
MatrixFree::reinit
void reinit(const MappingType &mapping, const std::vector< const DoFHandlerType * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< QuadratureType > &quad, const AdditionalData &additional_data=AdditionalData())
MatrixFree::cell_loop
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
MatrixFree::dof_info
std::vector< internal::MatrixFreeFunctions::DoFInfo > dof_info
Definition: matrix_free.h:2136
MatrixFree::initialize_cell_data_vector
void initialize_cell_data_vector(AlignedVector< T > &vec) const
MatrixFree::n_ghost_inner_face_batches
unsigned int n_ghost_inner_face_batches() const
MatrixFree::release_scratch_data_non_threadsafe
void release_scratch_data_non_threadsafe(const AlignedVector< Number > *memory) const
MatrixFree::renumber_dofs
void renumber_dofs(std::vector< types::global_dof_index > &renumbering, const unsigned int dof_handler_index=0)
MatrixFree::get_face_active_fe_index
unsigned int get_face_active_fe_index(const std::pair< unsigned int, unsigned int > range, const bool is_interior_face=true) const
MatrixFree::n_active_entries_per_face_batch
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
MatrixFree::indices_are_initialized
bool indices_are_initialized
Definition: matrix_free.h:2199
MatrixFree::face_info
internal::MatrixFreeFunctions::FaceInfo< VectorizedArrayType::size()> face_info
Definition: matrix_free.h:2194
MatrixFree::reinit
void reinit(const MappingType &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< number2 > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
MatrixFree::get_cell_range_category
unsigned int get_cell_range_category(const std::pair< unsigned int, unsigned int > cell_batch_range) const
MatrixFree::loop
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*boundary_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, 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::reinit
void reinit(const MappingType &mapping, const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< QuadratureType > &quad, const AdditionalData &additional_data=AdditionalData())
MatrixFree::get_faces_by_cells_boundary_id
std::array< types::boundary_id, VectorizedArrayType::size()> get_faces_by_cells_boundary_id(const unsigned int cell_batch_index, const unsigned int face_number) const
MatrixFree::constraint_pool_row_index
std::vector< unsigned int > constraint_pool_row_index
Definition: matrix_free.h:2150
MatrixFree::get_constrained_dofs
const std::vector< unsigned int > & get_constrained_dofs(const unsigned int dof_handler_index=0) const
MatrixFree::reinit
void reinit(const DoFHandler< dim > &dof_handler, const AffineConstraints< number2 > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
MatrixFree::cell_level_index_end_local
unsigned int cell_level_index_end_local
Definition: matrix_free.h:2180
MatrixFree::n_base_elements
unsigned int n_base_elements(const unsigned int dof_handler_index) const
MatrixFree::is_supported
static bool is_supported(const FiniteElement< dim, spacedim > &fe)
MatrixFree::dimension
static constexpr unsigned int dimension
Definition: matrix_free.h:131
Threads::ThreadLocalStorage
A class that provides a separate storage location on each thread that accesses the object.
Definition: thread_local_storage.h:103
DEAL_II_DEPRECATED
#define DEAL_II_DEPRECATED
Definition: config.h:164
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
dof_handler.h
UpdateFlags
UpdateFlags
Definition: fe_update_flags.h:67
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition: fe_update_flags.h:129
update_gradients
@ update_gradients
Shape function gradients.
Definition: fe_update_flags.h:84
update_default
@ update_default
No update.
Definition: fe_update_flags.h:69
second
Point< 2 > second
Definition: grid_out.cc:4604
first
Point< 2 > first
Definition: grid_out.cc:4603
level
unsigned int level
Definition: grid_out.cc:4606
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1473
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
AssertThrowMPI
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1790
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcIndexRange
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
StandardExceptions::ExcNotInitialized
static ::ExceptionBase & ExcNotInitialized()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
DoFHandler::cell_iterator
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:466
DoFHandler::active_cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:438
LinearAlgebra::distributed::Vector::local_element
Number local_element(const size_type local_index) const
LinearAlgebra::distributed::Vector::reinit
void reinit(const size_type size, const bool omit_zeroing_entries=false)
VectorOperation::add
@ add
Definition: vector_operation.h:53
dof_handler.h
la_parallel_vector.h
mapping_collection.h
mapping_q1.h
dof_info.h
mapping_info.h
TrilinosWrappers::internal::begin
VectorType::value_type * begin(VectorType &V)
Definition: trilinos_sparse_matrix.cc:53
Utilities::MPI::job_supports_mpi
bool job_supports_mpi()
Definition: mpi.cc:1014
internal::VectorImplementation::minimum_parallel_grain_size
unsigned int minimum_parallel_grain_size
Definition: parallel.cc:34
internal::reinit
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
numbers::invalid_boundary_id
const types::boundary_id invalid_boundary_id
Definition: types.h:244
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:201
parallel::apply_to_subranges
void apply_to_subranges(const RangeType &begin, const typename identity< RangeType >::type &end, const Function &f, const unsigned int grainsize)
Definition: parallel.h:382
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5588
types::boundary_id
unsigned int boundary_id
Definition: types.h:129
q_collection.h
shape_info.h
MatrixFree::AdditionalData::tasks_parallel_scheme
TasksParallelScheme tasks_parallel_scheme
Definition: matrix_free.h:343
MatrixFree::AdditionalData::mapping_update_flags_inner_faces
UpdateFlags mapping_update_flags_inner_faces
Definition: matrix_free.h:407
MatrixFree::AdditionalData::AdditionalData
AdditionalData(const TasksParallelScheme tasks_parallel_scheme=partition_partition, const unsigned int tasks_block_size=0, const UpdateFlags mapping_update_flags=update_gradients|update_JxW_values, const UpdateFlags mapping_update_flags_boundary_faces=update_default, const UpdateFlags mapping_update_flags_inner_faces=update_default, const UpdateFlags mapping_update_flags_faces_by_cells=update_default, const unsigned int mg_level=numbers::invalid_unsigned_int, const bool store_plain_indices=true, const bool initialize_indices=true, const bool initialize_mapping=true, const bool overlap_communication_computation=true, const bool hold_all_faces_to_owned_cells=false, const bool cell_vectorization_categories_strict=false, const bool allow_ghosted_vectors_in_loops=true)
Definition: matrix_free.h:217
MatrixFree::AdditionalData::cell_vectorization_category
std::vector< unsigned int > cell_vectorization_category
Definition: matrix_free.h:514
MatrixFree::AdditionalData::operator=
AdditionalData & operator=(const AdditionalData &other)
Definition: matrix_free.h:281
MatrixFree::AdditionalData::mapping_update_flags_boundary_faces
UpdateFlags mapping_update_flags_boundary_faces
Definition: matrix_free.h:387
MatrixFree::AdditionalData::mapping_update_flags
UpdateFlags mapping_update_flags
Definition: matrix_free.h:367
MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells
UpdateFlags mapping_update_flags_faces_by_cells
Definition: matrix_free.h:435
MatrixFree::AdditionalData::AdditionalData
AdditionalData(const AdditionalData &other)
Definition: matrix_free.h:254
MatrixFree::AdditionalData::tasks_block_size
unsigned int tasks_block_size
Definition: matrix_free.h:354
PointerComparison::equal
static bool equal(const T *p1, const T *p2)
Definition: template_constraints.h:393
internal::MatrixFreeFunctions::DoFInfo::cell_loop_pre_list_index
std::vector< unsigned int > cell_loop_pre_list_index
Definition: dof_info.h:777
internal::MatrixFreeFunctions::DoFInfo::cell_loop_post_list_index
std::vector< unsigned int > cell_loop_post_list_index
Definition: dof_info.h:790
internal::MatrixFreeFunctions::DoFInfo::vector_zero_range_list
std::vector< std::pair< unsigned int, unsigned int > > vector_zero_range_list
Definition: dof_info.h:770
internal::MatrixFreeFunctions::DoFInfo::vector_partitioner
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition: dof_info.h:621
internal::MatrixFreeFunctions::DoFInfo::cell_loop_pre_list
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_pre_list
Definition: dof_info.h:783
internal::MatrixFreeFunctions::DoFInfo::vector_zero_range_list_index
std::vector< unsigned int > vector_zero_range_list_index
Definition: dof_info.h:765
internal::MatrixFreeFunctions::DoFInfo::cell_loop_post_list
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_post_list
Definition: dof_info.h:796
task_info.h
template_constraints.h
thread_local_storage.h
type_traits.h
vector_operation.h
Generated by doxygen 1.9.5