Reference documentation for deal.II version 9.3.3
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
matrix_free.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2011 - 2021 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
17#ifndef dealii_matrix_free_h
18#define dealii_matrix_free_h
19
20#include <deal.II/base/config.h>
21
22#include <deal.II/base/aligned_vector.h>
23#include <deal.II/base/exceptions.h>
24#include <deal.II/base/quadrature.h>
25#include <deal.II/base/template_constraints.h>
26#include <deal.II/base/thread_local_storage.h>
27#include <deal.II/base/vectorization.h>
28
29#include <deal.II/dofs/dof_handler.h>
30
31#include <deal.II/fe/fe.h>
32#include <deal.II/fe/mapping.h>
33#include <deal.II/fe/mapping_q1.h>
34
35#include <deal.II/grid/grid_tools.h>
36
37#include <deal.II/hp/dof_handler.h>
38#include <deal.II/hp/mapping_collection.h>
39#include <deal.II/hp/q_collection.h>
40
41#include <deal.II/lac/affine_constraints.h>
42#include <deal.II/lac/block_vector_base.h>
43#include <deal.II/lac/la_parallel_vector.h>
44#include <deal.II/lac/vector_operation.h>
45
46#include <deal.II/matrix_free/dof_info.h>
47#include <deal.II/matrix_free/mapping_info.h>
48#include <deal.II/matrix_free/shape_info.h>
49#include <deal.II/matrix_free/task_info.h>
50#include <deal.II/matrix_free/type_traits.h>
51
52#include <cstdlib>
53#include <limits>
54#include <list>
55#include <memory>
56
57
58DEAL_II_NAMESPACE_OPEN
59
60
61
113template <int dim,
114 typename Number = double,
115 typename VectorizedArrayType = VectorizedArray<Number>>
116class MatrixFree : public Subscriptor
117{
118 static_assert(
119 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
120 "Type of Number and of VectorizedArrayType do not match.");
121
122public:
127 using value_type = Number;
128 using vectorized_value_type = VectorizedArrayType;
129
133 static const unsigned int dimension = dim;
134
181 struct AdditionalData
182 {
188 enum TasksParallelScheme
189 {
193 none = internal::MatrixFreeFunctions::TaskInfo::none,
197 partition_partition =
198 internal::MatrixFreeFunctions::TaskInfo::partition_partition,
202 partition_color =
203 internal::MatrixFreeFunctions::TaskInfo::partition_color,
208 color = internal::MatrixFreeFunctions::TaskInfo::color
209 };
210
211 // remove with level_mg_handler
212 DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
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 : tasks_parallel_scheme(tasks_parallel_scheme)
233 , tasks_block_size(tasks_block_size)
234 , mapping_update_flags(mapping_update_flags)
235 , mapping_update_flags_boundary_faces(mapping_update_flags_boundary_faces)
236 , mapping_update_flags_inner_faces(mapping_update_flags_inner_faces)
237 , mapping_update_flags_faces_by_cells(mapping_update_flags_faces_by_cells)
238 , mg_level(mg_level)
239 , level_mg_handler(this->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 , communicator_sm(MPI_COMM_SELF)
248 {}
249
253 AdditionalData(const AdditionalData &other)
254 : tasks_parallel_scheme(other.tasks_parallel_scheme)
255 , tasks_block_size(other.tasks_block_size)
256 , mapping_update_flags(other.mapping_update_flags)
257 , mapping_update_flags_boundary_faces(
258 other.mapping_update_flags_boundary_faces)
259 , mapping_update_flags_inner_faces(other.mapping_update_flags_inner_faces)
260 , mapping_update_flags_faces_by_cells(
261 other.mapping_update_flags_faces_by_cells)
262 , mg_level(other.mg_level)
263 , level_mg_handler(this->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 , communicator_sm(other.communicator_sm)
274 {}
275
276 // remove with level_mg_handler
277 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
278
282 AdditionalData &
283 operator=(const AdditionalData &other)
284 {
285 tasks_parallel_scheme = other.tasks_parallel_scheme;
286 tasks_block_size = other.tasks_block_size;
287 mapping_update_flags = other.mapping_update_flags;
288 mapping_update_flags_boundary_faces =
289 other.mapping_update_flags_boundary_faces;
290 mapping_update_flags_inner_faces = other.mapping_update_flags_inner_faces;
291 mapping_update_flags_faces_by_cells =
292 other.mapping_update_flags_faces_by_cells;
293 mg_level = other.mg_level;
294 store_plain_indices = other.store_plain_indices;
295 initialize_indices = other.initialize_indices;
296 initialize_mapping = other.initialize_mapping;
297 overlap_communication_computation =
298 other.overlap_communication_computation;
299 hold_all_faces_to_owned_cells = other.hold_all_faces_to_owned_cells;
300 cell_vectorization_category = other.cell_vectorization_category;
301 cell_vectorization_categories_strict =
302 other.cell_vectorization_categories_strict;
303 communicator_sm = other.communicator_sm;
304
305 return *this;
306 }
307
344 TasksParallelScheme tasks_parallel_scheme;
345
355 unsigned int tasks_block_size;
356
368 UpdateFlags mapping_update_flags;
369
389 UpdateFlags mapping_update_flags_boundary_faces;
390
410 UpdateFlags mapping_update_flags_inner_faces;
411
438 UpdateFlags mapping_update_flags_faces_by_cells;
439
448 unsigned int mg_level;
449
455 DEAL_II_DEPRECATED unsigned int &level_mg_handler;
456
463 bool store_plain_indices;
464
474 bool initialize_indices;
475
483 bool initialize_mapping;
484
497 bool overlap_communication_computation;
498
506 bool hold_all_faces_to_owned_cells;
507
524 std::vector<unsigned int> cell_vectorization_category;
525
535 bool cell_vectorization_categories_strict;
536
540 MPI_Comm communicator_sm;
541 };
542
550 MatrixFree();
551
555 MatrixFree(const MatrixFree<dim, Number, VectorizedArrayType> &other);
556
560 ~MatrixFree() override = default;
561
574 template <typename QuadratureType, typename number2, typename MappingType>
575 void
576 reinit(const MappingType & mapping,
577 const DoFHandler<dim> & dof_handler,
578 const AffineConstraints<number2> &constraint,
579 const QuadratureType & quad,
580 const AdditionalData & additional_data = AdditionalData());
581
588 template <typename QuadratureType, typename number2>
589 DEAL_II_DEPRECATED void
590 reinit(const DoFHandler<dim> & dof_handler,
591 const AffineConstraints<number2> &constraint,
592 const QuadratureType & quad,
593 const AdditionalData & additional_data = AdditionalData());
594
616 template <typename QuadratureType, typename number2, typename MappingType>
617 void
618 reinit(const MappingType & mapping,
619 const std::vector<const DoFHandler<dim> *> & dof_handler,
620 const std::vector<const AffineConstraints<number2> *> &constraint,
621 const std::vector<QuadratureType> & quad,
622 const AdditionalData &additional_data = AdditionalData());
623
629 template <typename QuadratureType,
630 typename number2,
631 typename DoFHandlerType,
632 typename MappingType>
633 DEAL_II_DEPRECATED void
634 reinit(const MappingType & mapping,
635 const std::vector<const DoFHandlerType *> & dof_handler,
636 const std::vector<const AffineConstraints<number2> *> &constraint,
637 const std::vector<QuadratureType> & quad,
638 const AdditionalData &additional_data = AdditionalData());
639
646 template <typename QuadratureType, typename number2>
647 DEAL_II_DEPRECATED void
648 reinit(const std::vector<const DoFHandler<dim> *> & dof_handler,
649 const std::vector<const AffineConstraints<number2> *> &constraint,
650 const std::vector<QuadratureType> & quad,
651 const AdditionalData &additional_data = AdditionalData());
652
658 template <typename QuadratureType, typename number2, typename DoFHandlerType>
659 DEAL_II_DEPRECATED void
660 reinit(const std::vector<const DoFHandlerType *> & dof_handler,
661 const std::vector<const AffineConstraints<number2> *> &constraint,
662 const std::vector<QuadratureType> & quad,
663 const AdditionalData &additional_data = AdditionalData());
664
672 template <typename QuadratureType, typename number2, typename MappingType>
673 void
674 reinit(const MappingType & mapping,
675 const std::vector<const DoFHandler<dim> *> & dof_handler,
676 const std::vector<const AffineConstraints<number2> *> &constraint,
677 const QuadratureType & quad,
678 const AdditionalData &additional_data = AdditionalData());
679
685 template <typename QuadratureType,
686 typename number2,
687 typename DoFHandlerType,
688 typename MappingType>
689 DEAL_II_DEPRECATED void
690 reinit(const MappingType & mapping,
691 const std::vector<const DoFHandlerType *> & dof_handler,
692 const std::vector<const AffineConstraints<number2> *> &constraint,
693 const QuadratureType & quad,
694 const AdditionalData &additional_data = AdditionalData());
695
702 template <typename QuadratureType, typename number2>
703 DEAL_II_DEPRECATED void
704 reinit(const std::vector<const DoFHandler<dim> *> & dof_handler,
705 const std::vector<const AffineConstraints<number2> *> &constraint,
706 const QuadratureType & quad,
707 const AdditionalData &additional_data = AdditionalData());
708
714 template <typename QuadratureType, typename number2, typename DoFHandlerType>
715 DEAL_II_DEPRECATED void
716 reinit(const std::vector<const DoFHandlerType *> & dof_handler,
717 const std::vector<const AffineConstraints<number2> *> &constraint,
718 const QuadratureType & quad,
719 const AdditionalData &additional_data = AdditionalData());
720
726 void
727 copy_from(
728 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free_base);
729
739 void
740 update_mapping(const Mapping<dim> &mapping);
741
745 void
746 update_mapping(const std::shared_ptr<hp::MappingCollection<dim>> &mapping);
747
752 void
753 clear();
754
756
768 enum class DataAccessOnFaces
769 {
776 none,
777
788 values,
789
798 values_all_faces,
799
810 gradients,
811
820 gradients_all_faces,
821
828 unspecified
829 };
830
875 template <typename OutVector, typename InVector>
876 void
877 cell_loop(const std::function<void(
878 const MatrixFree<dim, Number, VectorizedArrayType> &,
879 OutVector &,
880 const InVector &,
881 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
882 OutVector & dst,
883 const InVector & src,
884 const bool zero_dst_vector = false) const;
885
932 template <typename CLASS, typename OutVector, typename InVector>
933 void
934 cell_loop(void (CLASS::*cell_operation)(
935 const MatrixFree &,
936 OutVector &,
937 const InVector &,
938 const std::pair<unsigned int, unsigned int> &) const,
939 const CLASS * owning_class,
940 OutVector & dst,
941 const InVector &src,
942 const bool zero_dst_vector = false) const;
943
947 template <typename CLASS, typename OutVector, typename InVector>
948 void
949 cell_loop(void (CLASS::*cell_operation)(
950 const MatrixFree &,
951 OutVector &,
952 const InVector &,
953 const std::pair<unsigned int, unsigned int> &),
954 CLASS * owning_class,
955 OutVector & dst,
956 const InVector &src,
957 const bool zero_dst_vector = false) const;
958
1043 template <typename CLASS, typename OutVector, typename InVector>
1044 void
1045 cell_loop(void (CLASS::*cell_operation)(
1046 const MatrixFree &,
1047 OutVector &,
1048 const InVector &,
1049 const std::pair<unsigned int, unsigned int> &) const,
1050 const CLASS * owning_class,
1051 OutVector & dst,
1052 const InVector &src,
1053 const std::function<void(const unsigned int, const unsigned int)>
1054 &operation_before_loop,
1055 const std::function<void(const unsigned int, const unsigned int)>
1056 & operation_after_loop,
1057 const unsigned int dof_handler_index_pre_post = 0) const;
1058
1062 template <typename CLASS, typename OutVector, typename InVector>
1063 void
1064 cell_loop(void (CLASS::*cell_operation)(
1065 const MatrixFree &,
1066 OutVector &,
1067 const InVector &,
1068 const std::pair<unsigned int, unsigned int> &),
1069 CLASS * owning_class,
1070 OutVector & dst,
1071 const InVector &src,
1072 const std::function<void(const unsigned int, const unsigned int)>
1073 &operation_before_loop,
1074 const std::function<void(const unsigned int, const unsigned int)>
1075 & operation_after_loop,
1076 const unsigned int dof_handler_index_pre_post = 0) const;
1077
1082 template <typename OutVector, typename InVector>
1083 void
1084 cell_loop(const std::function<void(
1085 const MatrixFree<dim, Number, VectorizedArrayType> &,
1086 OutVector &,
1087 const InVector &,
1088 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1089 OutVector & dst,
1090 const InVector & src,
1091 const std::function<void(const unsigned int, const unsigned int)>
1092 &operation_before_loop,
1093 const std::function<void(const unsigned int, const unsigned int)>
1094 & operation_after_loop,
1095 const unsigned int dof_handler_index_pre_post = 0) const;
1096
1172 template <typename OutVector, typename InVector>
1173 void
1174 loop(const std::function<
1175 void(const MatrixFree<dim, Number, VectorizedArrayType> &,
1176 OutVector &,
1177 const InVector &,
1178 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1179 const std::function<
1180 void(const MatrixFree<dim, Number, VectorizedArrayType> &,
1181 OutVector &,
1182 const InVector &,
1183 const std::pair<unsigned int, unsigned int> &)> &face_operation,
1184 const std::function<void(
1185 const MatrixFree<dim, Number, VectorizedArrayType> &,
1186 OutVector &,
1187 const InVector &,
1188 const std::pair<unsigned int, unsigned int> &)> &boundary_operation,
1189 OutVector & dst,
1190 const InVector & src,
1191 const bool zero_dst_vector = false,
1192 const DataAccessOnFaces dst_vector_face_access =
1193 DataAccessOnFaces::unspecified,
1194 const DataAccessOnFaces src_vector_face_access =
1195 DataAccessOnFaces::unspecified) const;
1196
1282 template <typename CLASS, typename OutVector, typename InVector>
1283 void
1284 loop(
1285 void (CLASS::*cell_operation)(const MatrixFree &,
1286 OutVector &,
1287 const InVector &,
1288 const std::pair<unsigned int, unsigned int> &)
1289 const,
1290 void (CLASS::*face_operation)(const MatrixFree &,
1291 OutVector &,
1292 const InVector &,
1293 const std::pair<unsigned int, unsigned int> &)
1294 const,
1295 void (CLASS::*boundary_operation)(
1296 const MatrixFree &,
1297 OutVector &,
1298 const InVector &,
1299 const std::pair<unsigned int, unsigned int> &) const,
1300 const CLASS * owning_class,
1301 OutVector & dst,
1302 const InVector & src,
1303 const bool zero_dst_vector = false,
1304 const DataAccessOnFaces dst_vector_face_access =
1305 DataAccessOnFaces::unspecified,
1306 const DataAccessOnFaces src_vector_face_access =
1307 DataAccessOnFaces::unspecified) const;
1308
1312 template <typename CLASS, typename OutVector, typename InVector>
1313 void
1314 loop(void (CLASS::*cell_operation)(
1315 const MatrixFree &,
1316 OutVector &,
1317 const InVector &,
1318 const std::pair<unsigned int, unsigned int> &),
1319 void (CLASS::*face_operation)(
1320 const MatrixFree &,
1321 OutVector &,
1322 const InVector &,
1323 const std::pair<unsigned int, unsigned int> &),
1324 void (CLASS::*boundary_operation)(
1325 const MatrixFree &,
1326 OutVector &,
1327 const InVector &,
1328 const std::pair<unsigned int, unsigned int> &),
1329 CLASS * owning_class,
1330 OutVector & dst,
1331 const InVector & src,
1332 const bool zero_dst_vector = false,
1333 const DataAccessOnFaces dst_vector_face_access =
1334 DataAccessOnFaces::unspecified,
1335 const DataAccessOnFaces src_vector_face_access =
1336 DataAccessOnFaces::unspecified) const;
1337
1402 template <typename CLASS, typename OutVector, typename InVector>
1403 void
1404 loop_cell_centric(void (CLASS::*cell_operation)(
1405 const MatrixFree &,
1406 OutVector &,
1407 const InVector &,
1408 const std::pair<unsigned int, unsigned int> &) const,
1409 const CLASS * owning_class,
1410 OutVector & dst,
1411 const InVector & src,
1412 const bool zero_dst_vector = false,
1413 const DataAccessOnFaces src_vector_face_access =
1414 DataAccessOnFaces::unspecified) const;
1415
1419 template <typename CLASS, typename OutVector, typename InVector>
1420 void
1421 loop_cell_centric(void (CLASS::*cell_operation)(
1422 const MatrixFree &,
1423 OutVector &,
1424 const InVector &,
1425 const std::pair<unsigned int, unsigned int> &),
1426 CLASS * owning_class,
1427 OutVector & dst,
1428 const InVector & src,
1429 const bool zero_dst_vector = false,
1430 const DataAccessOnFaces src_vector_face_access =
1431 DataAccessOnFaces::unspecified) const;
1432
1436 template <typename OutVector, typename InVector>
1437 void
1438 loop_cell_centric(
1439 const std::function<void(const MatrixFree &,
1440 OutVector &,
1441 const InVector &,
1442 const std::pair<unsigned int, unsigned int> &)>
1443 & cell_operation,
1444 OutVector & dst,
1445 const InVector & src,
1446 const bool zero_dst_vector = false,
1447 const DataAccessOnFaces src_vector_face_access =
1448 DataAccessOnFaces::unspecified) const;
1449
1457 std::pair<unsigned int, unsigned int>
1458 create_cell_subrange_hp(const std::pair<unsigned int, unsigned int> &range,
1459 const unsigned int fe_degree,
1460 const unsigned int dof_handler_index = 0) const;
1461
1468 std::pair<unsigned int, unsigned int>
1469 create_cell_subrange_hp_by_index(
1470 const std::pair<unsigned int, unsigned int> &range,
1471 const unsigned int fe_index,
1472 const unsigned int dof_handler_index = 0) const;
1473
1477 unsigned int
1478 n_active_fe_indices() const;
1479
1483 unsigned int
1484 get_cell_active_fe_index(
1485 const std::pair<unsigned int, unsigned int> range) const;
1486
1490 unsigned int
1491 get_face_active_fe_index(const std::pair<unsigned int, unsigned int> range,
1492 const bool is_interior_face = true) const;
1493
1495
1520 template <typename VectorType>
1521 void
1522 initialize_dof_vector(VectorType & vec,
1523 const unsigned int dof_handler_index = 0) const;
1524
1545 template <typename Number2>
1546 void
1547 initialize_dof_vector(LinearAlgebra::distributed::Vector<Number2> &vec,
1548 const unsigned int dof_handler_index = 0) const;
1549
1560 const std::shared_ptr<const Utilities::MPI::Partitioner> &
1561 get_vector_partitioner(const unsigned int dof_handler_index = 0) const;
1562
1566 const IndexSet &
1567 get_locally_owned_set(const unsigned int dof_handler_index = 0) const;
1568
1572 const IndexSet &
1573 get_ghost_set(const unsigned int dof_handler_index = 0) const;
1574
1584 const std::vector<unsigned int> &
1585 get_constrained_dofs(const unsigned int dof_handler_index = 0) const;
1586
1597 void
1598 renumber_dofs(std::vector<types::global_dof_index> &renumbering,
1599 const unsigned int dof_handler_index = 0);
1600
1602
1610 template <int spacedim>
1611 static bool
1612 is_supported(const FiniteElement<dim, spacedim> &fe);
1613
1617 unsigned int
1618 n_components() const;
1619
1624 unsigned int
1625 n_base_elements(const unsigned int dof_handler_index) const;
1626
1634 unsigned int
1635 n_physical_cells() const;
1636
1640 DEAL_II_DEPRECATED unsigned int
1641 n_macro_cells() const;
1642
1652 unsigned int
1653 n_cell_batches() const;
1654
1661 unsigned int
1662 n_ghost_cell_batches() const;
1663
1671 unsigned int
1672 n_inner_face_batches() const;
1673
1682 unsigned int
1683 n_boundary_face_batches() const;
1684
1689 unsigned int
1690 n_ghost_inner_face_batches() const;
1691
1698 types::boundary_id
1699 get_boundary_id(const unsigned int macro_face) const;
1700
1705 std::array<types::boundary_id, VectorizedArrayType::size()>
1706 get_faces_by_cells_boundary_id(const unsigned int cell_batch_index,
1707 const unsigned int face_number) const;
1708
1713 const DoFHandler<dim> &
1714 get_dof_handler(const unsigned int dof_handler_index = 0) const;
1715
1726 template <typename DoFHandlerType>
1727 DEAL_II_DEPRECATED const DoFHandlerType &
1728 get_dof_handler(const unsigned int dof_handler_index = 0) const;
1729
1742 typename DoFHandler<dim>::cell_iterator
1743 get_cell_iterator(const unsigned int cell_batch_index,
1744 const unsigned int lane_index,
1745 const unsigned int dof_handler_index = 0) const;
1746
1752 std::pair<int, int>
1753 get_cell_level_and_index(const unsigned int cell_batch_index,
1754 const unsigned int lane_index) const;
1755
1768 std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>
1769 get_face_iterator(const unsigned int face_batch_index,
1770 const unsigned int lane_index,
1771 const bool interior = true,
1772 const unsigned int fe_component = 0) const;
1773
1779 DEAL_II_DEPRECATED typename DoFHandler<dim>::active_cell_iterator
1780 get_hp_cell_iterator(const unsigned int cell_batch_index,
1781 const unsigned int lane_index,
1782 const unsigned int dof_handler_index = 0) const;
1783
1796 bool
1797 at_irregular_cell(const unsigned int cell_batch_index) const;
1798
1802 DEAL_II_DEPRECATED unsigned int
1803 n_components_filled(const unsigned int cell_batch_number) const;
1804
1814 unsigned int
1815 n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const;
1816
1826 unsigned int
1827 n_active_entries_per_face_batch(const unsigned int face_batch_index) const;
1828
1832 unsigned int
1833 get_dofs_per_cell(const unsigned int dof_handler_index = 0,
1834 const unsigned int hp_active_fe_index = 0) const;
1835
1839 unsigned int
1840 get_n_q_points(const unsigned int quad_index = 0,
1841 const unsigned int hp_active_fe_index = 0) const;
1842
1847 unsigned int
1848 get_dofs_per_face(const unsigned int dof_handler_index = 0,
1849 const unsigned int hp_active_fe_index = 0) const;
1850
1855 unsigned int
1856 get_n_q_points_face(const unsigned int quad_index = 0,
1857 const unsigned int hp_active_fe_index = 0) const;
1858
1862 const Quadrature<dim> &
1863 get_quadrature(const unsigned int quad_index = 0,
1864 const unsigned int hp_active_fe_index = 0) const;
1865
1869 const Quadrature<dim - 1> &
1870 get_face_quadrature(const unsigned int quad_index = 0,
1871 const unsigned int hp_active_fe_index = 0) const;
1872
1879 unsigned int
1880 get_cell_category(const unsigned int cell_batch_index) const;
1881
1886 std::pair<unsigned int, unsigned int>
1887 get_face_category(const unsigned int macro_face) const;
1888
1892 bool
1893 indices_initialized() const;
1894
1899 bool
1900 mapping_initialized() const;
1901
1906 unsigned int
1907 get_mg_level() const;
1908
1913 std::size_t
1914 memory_consumption() const;
1915
1920 template <typename StreamType>
1921 void
1922 print_memory_consumption(StreamType &out) const;
1923
1928 void
1929 print(std::ostream &out) const;
1930
1932
1942 const internal::MatrixFreeFunctions::TaskInfo &
1943 get_task_info() const;
1944
1945 /*
1946 * Return geometry-dependent information on the cells.
1947 */
1948 const internal::MatrixFreeFunctions::
1949 MappingInfo<dim, Number, VectorizedArrayType> &
1950 get_mapping_info() const;
1951
1955 const internal::MatrixFreeFunctions::DoFInfo &
1956 get_dof_info(const unsigned int dof_handler_index_component = 0) const;
1957
1961 unsigned int
1962 n_constraint_pool_entries() const;
1963
1968 const Number *
1969 constraint_pool_begin(const unsigned int pool_index) const;
1970
1976 const Number *
1977 constraint_pool_end(const unsigned int pool_index) const;
1978
1982 const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
1983 get_shape_info(const unsigned int dof_handler_index_component = 0,
1984 const unsigned int quad_index = 0,
1985 const unsigned int fe_base_element = 0,
1986 const unsigned int hp_active_fe_index = 0,
1987 const unsigned int hp_active_quad_index = 0) const;
1988
1992 const internal::MatrixFreeFunctions::FaceToCellTopology<
1993 VectorizedArrayType::size()> &
1994 get_face_info(const unsigned int face_batch_index) const;
1995
1996
2002 const Table<3, unsigned int> &
2003 get_cell_and_face_to_plain_faces() const;
2004
2018 AlignedVector<VectorizedArrayType> *
2019 acquire_scratch_data() const;
2020
2024 void
2025 release_scratch_data(const AlignedVector<VectorizedArrayType> *memory) const;
2026
2036 AlignedVector<Number> *
2037 acquire_scratch_data_non_threadsafe() const;
2038
2042 void
2043 release_scratch_data_non_threadsafe(
2044 const AlignedVector<Number> *memory) const;
2045
2047
2048private:
2053 template <typename number2, int q_dim>
2054 void
2055 internal_reinit(
2056 const std::shared_ptr<hp::MappingCollection<dim>> & mapping,
2057 const std::vector<const DoFHandler<dim, dim> *> & dof_handlers,
2058 const std::vector<const AffineConstraints<number2> *> &constraint,
2059 const std::vector<IndexSet> & locally_owned_set,
2060 const std::vector<hp::QCollection<q_dim>> & quad,
2061 const AdditionalData & additional_data);
2062
2069 template <typename number2>
2070 void
2071 initialize_indices(
2072 const std::vector<const AffineConstraints<number2> *> &constraint,
2073 const std::vector<IndexSet> & locally_owned_set,
2074 const AdditionalData & additional_data);
2075
2079 void
2080 initialize_dof_handlers(
2081 const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
2082 const AdditionalData & additional_data);
2083
2087 std::vector<SmartPointer<const DoFHandler<dim>>> dof_handlers;
2088
2093 std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
2094
2101 std::vector<Number> constraint_pool_data;
2102
2107 std::vector<unsigned int> constraint_pool_row_index;
2108
2113 internal::MatrixFreeFunctions::MappingInfo<dim, Number, VectorizedArrayType>
2114 mapping_info;
2115
2119 Table<4, internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>>
2120 shape_info;
2121
2128 std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
2129
2130
2137 unsigned int cell_level_index_end_local;
2138
2144 internal::MatrixFreeFunctions::TaskInfo task_info;
2145
2150 internal::MatrixFreeFunctions::FaceInfo<VectorizedArrayType::size()>
2151 face_info;
2152
2156 bool indices_are_initialized;
2157
2161 bool mapping_is_initialized;
2162
2170 mutable Threads::ThreadLocalStorage<
2171 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2172 scratch_pad;
2173
2178 mutable std::list<std::pair<bool, AlignedVector<Number>>>
2179 scratch_pad_non_threadsafe;
2180
2184 unsigned int mg_level;
2185};
2186
2187
2188
2189/*----------------------- Inline functions ----------------------------------*/
2190
2191#ifndef DOXYGEN
2192
2193
2194
2195template <int dim, typename Number, typename VectorizedArrayType>
2196template <typename VectorType>
2197inline void
2198MatrixFree<dim, Number, VectorizedArrayType>::initialize_dof_vector(
2199 VectorType & vec,
2200 const unsigned int comp) const
2201{
2202 AssertIndexRange(comp, n_components());
2203 vec.reinit(dof_info[comp].vector_partitioner->size());
2204}
2205
2206
2207
2208template <int dim, typename Number, typename VectorizedArrayType>
2209template <typename Number2>
2210inline void
2211MatrixFree<dim, Number, VectorizedArrayType>::initialize_dof_vector(
2212 LinearAlgebra::distributed::Vector<Number2> &vec,
2213 const unsigned int comp) const
2214{
2215 AssertIndexRange(comp, n_components());
2216 vec.reinit(dof_info[comp].vector_partitioner, task_info.communicator_sm);
2217}
2218
2219
2220
2221template <int dim, typename Number, typename VectorizedArrayType>
2222inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2223MatrixFree<dim, Number, VectorizedArrayType>::get_vector_partitioner(
2224 const unsigned int comp) const
2225{
2226 AssertIndexRange(comp, n_components());
2227 return dof_info[comp].vector_partitioner;
2228}
2229
2230
2231
2232template <int dim, typename Number, typename VectorizedArrayType>
2233inline const std::vector<unsigned int> &
2234MatrixFree<dim, Number, VectorizedArrayType>::get_constrained_dofs(
2235 const unsigned int comp) const
2236{
2237 AssertIndexRange(comp, n_components());
2238 return dof_info[comp].constrained_dofs;
2239}
2240
2241
2242
2243template <int dim, typename Number, typename VectorizedArrayType>
2244inline unsigned int
2245MatrixFree<dim, Number, VectorizedArrayType>::n_components() const
2246{
2247 AssertDimension(dof_handlers.size(), dof_info.size());
2248 return dof_handlers.size();
2249}
2250
2251
2252
2253template <int dim, typename Number, typename VectorizedArrayType>
2254inline unsigned int
2255MatrixFree<dim, Number, VectorizedArrayType>::n_base_elements(
2256 const unsigned int dof_no) const
2257{
2258 AssertDimension(dof_handlers.size(), dof_info.size());
2259 AssertIndexRange(dof_no, dof_handlers.size());
2260 return dof_handlers[dof_no]->get_fe().n_base_elements();
2261}
2262
2263
2264
2265template <int dim, typename Number, typename VectorizedArrayType>
2266inline const internal::MatrixFreeFunctions::TaskInfo &
2267MatrixFree<dim, Number, VectorizedArrayType>::get_task_info() const
2268{
2269 return task_info;
2270}
2271
2272
2273
2274template <int dim, typename Number, typename VectorizedArrayType>
2275inline unsigned int
2276MatrixFree<dim, Number, VectorizedArrayType>::n_macro_cells() const
2277{
2278 return *(task_info.cell_partition_data.end() - 2);
2279}
2280
2281
2282
2283template <int dim, typename Number, typename VectorizedArrayType>
2284inline unsigned int
2285MatrixFree<dim, Number, VectorizedArrayType>::n_physical_cells() const
2286{
2287 return task_info.n_active_cells;
2288}
2289
2290
2291
2292template <int dim, typename Number, typename VectorizedArrayType>
2293inline unsigned int
2294MatrixFree<dim, Number, VectorizedArrayType>::n_cell_batches() const
2295{
2296 return *(task_info.cell_partition_data.end() - 2);
2297}
2298
2299
2300
2301template <int dim, typename Number, typename VectorizedArrayType>
2302inline unsigned int
2303MatrixFree<dim, Number, VectorizedArrayType>::n_ghost_cell_batches() const
2304{
2305 return *(task_info.cell_partition_data.end() - 1) -
2306 *(task_info.cell_partition_data.end() - 2);
2307}
2308
2309
2310
2311template <int dim, typename Number, typename VectorizedArrayType>
2312inline unsigned int
2313MatrixFree<dim, Number, VectorizedArrayType>::n_inner_face_batches() const
2314{
2315 if (task_info.face_partition_data.size() == 0)
2316 return 0;
2317 return task_info.face_partition_data.back();
2318}
2319
2320
2321
2322template <int dim, typename Number, typename VectorizedArrayType>
2323inline unsigned int
2324MatrixFree<dim, Number, VectorizedArrayType>::n_boundary_face_batches() const
2325{
2326 if (task_info.face_partition_data.size() == 0)
2327 return 0;
2328 return task_info.boundary_partition_data.back() -
2329 task_info.face_partition_data.back();
2330}
2331
2332
2333
2334template <int dim, typename Number, typename VectorizedArrayType>
2335inline unsigned int
2336MatrixFree<dim, Number, VectorizedArrayType>::n_ghost_inner_face_batches() const
2337{
2338 if (task_info.face_partition_data.size() == 0)
2339 return 0;
2340 return face_info.faces.size() - task_info.boundary_partition_data.back();
2341}
2342
2343
2344
2345template <int dim, typename Number, typename VectorizedArrayType>
2346inline types::boundary_id
2347MatrixFree<dim, Number, VectorizedArrayType>::get_boundary_id(
2348 const unsigned int macro_face) const
2349{
2350 Assert(macro_face >= task_info.boundary_partition_data[0] &&
2351 macro_face < task_info.boundary_partition_data.back(),
2352 ExcIndexRange(macro_face,
2353 task_info.boundary_partition_data[0],
2354 task_info.boundary_partition_data.back()));
2355 return types::boundary_id(face_info.faces[macro_face].exterior_face_no);
2356}
2357
2358
2359
2360template <int dim, typename Number, typename VectorizedArrayType>
2361inline std::array<types::boundary_id, VectorizedArrayType::size()>
2362MatrixFree<dim, Number, VectorizedArrayType>::get_faces_by_cells_boundary_id(
2363 const unsigned int cell_batch_index,
2364 const unsigned int face_number) const
2365{
2366 AssertIndexRange(cell_batch_index, n_cell_batches());
2367 AssertIndexRange(face_number, GeometryInfo<dim>::faces_per_cell);
2368 Assert(face_info.cell_and_face_boundary_id.size(0) >= n_cell_batches(),
2369 ExcNotInitialized());
2370 std::array<types::boundary_id, VectorizedArrayType::size()> result;
2371 result.fill(numbers::invalid_boundary_id);
2372 for (unsigned int v = 0;
2373 v < n_active_entries_per_cell_batch(cell_batch_index);
2374 ++v)
2375 result[v] =
2376 face_info.cell_and_face_boundary_id(cell_batch_index, face_number, v);
2377 return result;
2378}
2379
2380
2381
2382template <int dim, typename Number, typename VectorizedArrayType>
2383inline const internal::MatrixFreeFunctions::
2384 MappingInfo<dim, Number, VectorizedArrayType> &
2385 MatrixFree<dim, Number, VectorizedArrayType>::get_mapping_info() const
2386{
2387 return mapping_info;
2388}
2389
2390
2391
2392template <int dim, typename Number, typename VectorizedArrayType>
2393inline const internal::MatrixFreeFunctions::DoFInfo &
2394MatrixFree<dim, Number, VectorizedArrayType>::get_dof_info(
2395 const unsigned int dof_index) const
2396{
2397 AssertIndexRange(dof_index, n_components());
2398 return dof_info[dof_index];
2399}
2400
2401
2402
2403template <int dim, typename Number, typename VectorizedArrayType>
2404inline unsigned int
2405MatrixFree<dim, Number, VectorizedArrayType>::n_constraint_pool_entries() const
2406{
2407 return constraint_pool_row_index.size() - 1;
2408}
2409
2410
2411
2412template <int dim, typename Number, typename VectorizedArrayType>
2413inline const Number *
2414MatrixFree<dim, Number, VectorizedArrayType>::constraint_pool_begin(
2415 const unsigned int row) const
2416{
2417 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2418 return constraint_pool_data.empty() ?
2419 nullptr :
2420 constraint_pool_data.data() + constraint_pool_row_index[row];
2421}
2422
2423
2424
2425template <int dim, typename Number, typename VectorizedArrayType>
2426inline const Number *
2427MatrixFree<dim, Number, VectorizedArrayType>::constraint_pool_end(
2428 const unsigned int row) const
2429{
2430 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2431 return constraint_pool_data.empty() ?
2432 nullptr :
2433 constraint_pool_data.data() + constraint_pool_row_index[row + 1];
2434}
2435
2436
2437
2438template <int dim, typename Number, typename VectorizedArrayType>
2439inline std::pair<unsigned int, unsigned int>
2440MatrixFree<dim, Number, VectorizedArrayType>::create_cell_subrange_hp(
2441 const std::pair<unsigned int, unsigned int> &range,
2442 const unsigned int degree,
2443 const unsigned int dof_handler_component) const
2444{
2445 if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2446 {
2447 AssertDimension(
2448 dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2449 AssertDimension(
2450 dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2451 if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2452 return range;
2453 else
2454 return {range.second, range.second};
2455 }
2456
2457 const unsigned int fe_index =
2458 dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2459 if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2460 return {range.second, range.second};
2461 else
2462 return create_cell_subrange_hp_by_index(range,
2463 fe_index,
2464 dof_handler_component);
2465}
2466
2467
2468
2469template <int dim, typename Number, typename VectorizedArrayType>
2470inline bool
2471MatrixFree<dim, Number, VectorizedArrayType>::at_irregular_cell(
2472 const unsigned int cell_batch_index) const
2473{
2474 AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2475 return VectorizedArrayType::size() > 1 &&
2476 cell_level_index[(cell_batch_index + 1) * VectorizedArrayType::size() -
2477 1] == cell_level_index[(cell_batch_index + 1) *
2478 VectorizedArrayType::size() -
2479 2];
2480}
2481
2482
2483
2484template <int dim, typename Number, typename VectorizedArrayType>
2485unsigned int
2486MatrixFree<dim, Number, VectorizedArrayType>::n_active_fe_indices() const
2487{
2488 return shape_info.size(2);
2489}
2490
2491
2492template <int dim, typename Number, typename VectorizedArrayType>
2493unsigned int
2494MatrixFree<dim, Number, VectorizedArrayType>::get_cell_active_fe_index(
2495 const std::pair<unsigned int, unsigned int> range) const
2496{
2497 const auto &fe_indices = dof_info[0].cell_active_fe_index;
2498
2499 if (fe_indices.empty() == true)
2500 return 0;
2501
2502 const auto index = fe_indices[range.first];
2503
2504 for (unsigned int i = range.first; i < range.second; ++i)
2505 AssertDimension(index, fe_indices[i]);
2506
2507 return index;
2508}
2509
2510
2511
2512template <int dim, typename Number, typename VectorizedArrayType>
2513unsigned int
2514MatrixFree<dim, Number, VectorizedArrayType>::get_face_active_fe_index(
2515 const std::pair<unsigned int, unsigned int> range,
2516 const bool is_interior_face) const
2517{
2518 const auto &fe_indices = dof_info[0].cell_active_fe_index;
2519
2520 if (fe_indices.empty() == true)
2521 return 0;
2522
2523 if (is_interior_face)
2524 {
2525 const unsigned int index =
2526 fe_indices[face_info.faces[range.first].cells_interior[0] /
2527 VectorizedArrayType::size()];
2528
2529 for (unsigned int i = range.first; i < range.second; ++i)
2530 AssertDimension(index,
2531 fe_indices[face_info.faces[i].cells_interior[0] /
2532 VectorizedArrayType::size()]);
2533
2534 return index;
2535 }
2536 else
2537 {
2538 const unsigned int index =
2539 fe_indices[face_info.faces[range.first].cells_exterior[0] /
2540 VectorizedArrayType::size()];
2541
2542 for (unsigned int i = range.first; i < range.second; ++i)
2543 AssertDimension(index,
2544 fe_indices[face_info.faces[i].cells_exterior[0] /
2545 VectorizedArrayType::size()]);
2546
2547 return index;
2548 }
2549}
2550
2551
2552
2553template <int dim, typename Number, typename VectorizedArrayType>
2554inline unsigned int
2555MatrixFree<dim, Number, VectorizedArrayType>::n_components_filled(
2556 const unsigned int cell_batch_index) const
2557{
2558 return n_active_entries_per_cell_batch(cell_batch_index);
2559}
2560
2561
2562
2563template <int dim, typename Number, typename VectorizedArrayType>
2564inline unsigned int
2565MatrixFree<dim, Number, VectorizedArrayType>::n_active_entries_per_cell_batch(
2566 const unsigned int cell_batch_index) const
2567{
2568 AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2569 unsigned int n_lanes = VectorizedArrayType::size();
2570 while (n_lanes > 1 &&
2571 cell_level_index[cell_batch_index * VectorizedArrayType::size() +
2572 n_lanes - 1] ==
2573 cell_level_index[cell_batch_index * VectorizedArrayType::size() +
2574 n_lanes - 2])
2575 --n_lanes;
2576 AssertIndexRange(n_lanes - 1, VectorizedArrayType::size());
2577 return n_lanes;
2578}
2579
2580
2581
2582template <int dim, typename Number, typename VectorizedArrayType>
2583inline unsigned int
2584MatrixFree<dim, Number, VectorizedArrayType>::n_active_entries_per_face_batch(
2585 const unsigned int face_batch_index) const
2586{
2587 AssertIndexRange(face_batch_index, face_info.faces.size());
2588 unsigned int n_lanes = VectorizedArrayType::size();
2589 while (n_lanes > 1 &&
2590 face_info.faces[face_batch_index].cells_interior[n_lanes - 1] ==
2591 numbers::invalid_unsigned_int)
2592 --n_lanes;
2593 AssertIndexRange(n_lanes - 1, VectorizedArrayType::size());
2594 return n_lanes;
2595}
2596
2597
2598
2599template <int dim, typename Number, typename VectorizedArrayType>
2600inline unsigned int
2601MatrixFree<dim, Number, VectorizedArrayType>::get_dofs_per_cell(
2602 const unsigned int dof_handler_index,
2603 const unsigned int active_fe_index) const
2604{
2605 return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2606}
2607
2608
2609
2610template <int dim, typename Number, typename VectorizedArrayType>
2611inline unsigned int
2612MatrixFree<dim, Number, VectorizedArrayType>::get_n_q_points(
2613 const unsigned int quad_index,
2614 const unsigned int active_fe_index) const
2615{
2616 AssertIndexRange(quad_index, mapping_info.cell_data.size());
2617 return mapping_info.cell_data[quad_index]
2618 .descriptor[active_fe_index]
2619 .n_q_points;
2620}
2621
2622
2623
2624template <int dim, typename Number, typename VectorizedArrayType>
2625inline unsigned int
2626MatrixFree<dim, Number, VectorizedArrayType>::get_dofs_per_face(
2627 const unsigned int dof_handler_index,
2628 const unsigned int active_fe_index) const
2629{
2630 return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2631}
2632
2633
2634
2635template <int dim, typename Number, typename VectorizedArrayType>
2636inline unsigned int
2637MatrixFree<dim, Number, VectorizedArrayType>::get_n_q_points_face(
2638 const unsigned int quad_index,
2639 const unsigned int active_fe_index) const
2640{
2641 AssertIndexRange(quad_index, mapping_info.face_data.size());
2642 return mapping_info.face_data[quad_index]
2643 .descriptor[active_fe_index]
2644 .n_q_points;
2645}
2646
2647
2648
2649template <int dim, typename Number, typename VectorizedArrayType>
2650inline const IndexSet &
2651MatrixFree<dim, Number, VectorizedArrayType>::get_locally_owned_set(
2652 const unsigned int dof_handler_index) const
2653{
2654 return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2655}
2656
2657
2658
2659template <int dim, typename Number, typename VectorizedArrayType>
2660inline const IndexSet &
2661MatrixFree<dim, Number, VectorizedArrayType>::get_ghost_set(
2662 const unsigned int dof_handler_index) const
2663{
2664 return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2665}
2666
2667
2668
2669template <int dim, typename Number, typename VectorizedArrayType>
2670inline const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
2671MatrixFree<dim, Number, VectorizedArrayType>::get_shape_info(
2672 const unsigned int dof_handler_index,
2673 const unsigned int index_quad,
2674 const unsigned int index_fe,
2675 const unsigned int active_fe_index,
2676 const unsigned int active_quad_index) const
2677{
2678 AssertIndexRange(dof_handler_index, dof_info.size());
2679 const unsigned int ind =
2680 dof_info[dof_handler_index].global_base_element_offset + index_fe;
2681 AssertIndexRange(ind, shape_info.size(0));
2682 AssertIndexRange(index_quad, shape_info.size(1));
2683 AssertIndexRange(active_fe_index, shape_info.size(2));
2684 AssertIndexRange(active_quad_index, shape_info.size(3));
2685 return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2686}
2687
2688
2689
2690template <int dim, typename Number, typename VectorizedArrayType>
2691inline const internal::MatrixFreeFunctions::FaceToCellTopology<
2692 VectorizedArrayType::size()> &
2693MatrixFree<dim, Number, VectorizedArrayType>::get_face_info(
2694 const unsigned int macro_face) const
2695{
2696 AssertIndexRange(macro_face, face_info.faces.size());
2697 return face_info.faces[macro_face];
2698}
2699
2700
2701
2702template <int dim, typename Number, typename VectorizedArrayType>
2703inline const Table<3, unsigned int> &
2704MatrixFree<dim, Number, VectorizedArrayType>::get_cell_and_face_to_plain_faces()
2705 const
2706{
2707 return face_info.cell_and_face_to_plain_faces;
2708}
2709
2710
2711
2712template <int dim, typename Number, typename VectorizedArrayType>
2713inline const Quadrature<dim> &
2714MatrixFree<dim, Number, VectorizedArrayType>::get_quadrature(
2715 const unsigned int quad_index,
2716 const unsigned int active_fe_index) const
2717{
2718 AssertIndexRange(quad_index, mapping_info.cell_data.size());
2719 return mapping_info.cell_data[quad_index]
2720 .descriptor[active_fe_index]
2721 .quadrature;
2722}
2723
2724
2725
2726template <int dim, typename Number, typename VectorizedArrayType>
2727inline const Quadrature<dim - 1> &
2728MatrixFree<dim, Number, VectorizedArrayType>::get_face_quadrature(
2729 const unsigned int quad_index,
2730 const unsigned int active_fe_index) const
2731{
2732 AssertIndexRange(quad_index, mapping_info.face_data.size());
2733 return mapping_info.face_data[quad_index]
2734 .descriptor[active_fe_index]
2735 .quadrature;
2736}
2737
2738
2739
2740template <int dim, typename Number, typename VectorizedArrayType>
2741inline unsigned int
2742MatrixFree<dim, Number, VectorizedArrayType>::get_cell_category(
2743 const unsigned int cell_batch_index) const
2744{
2745 AssertIndexRange(0, dof_info.size());
2746 AssertIndexRange(cell_batch_index, dof_info[0].cell_active_fe_index.size());
2747 if (dof_info[0].cell_active_fe_index.empty())
2748 return 0;
2749 else
2750 return dof_info[0].cell_active_fe_index[cell_batch_index];
2751}
2752
2753
2754
2755template <int dim, typename Number, typename VectorizedArrayType>
2756inline std::pair<unsigned int, unsigned int>
2757MatrixFree<dim, Number, VectorizedArrayType>::get_face_category(
2758 const unsigned int macro_face) const
2759{
2760 AssertIndexRange(macro_face, face_info.faces.size());
2761 if (dof_info[0].cell_active_fe_index.empty())
2762 return std::make_pair(0U, 0U);
2763
2764 std::pair<unsigned int, unsigned int> result;
2765 for (unsigned int v = 0; v < VectorizedArrayType::size() &&
2766 face_info.faces[macro_face].cells_interior[v] !=
2767 numbers::invalid_unsigned_int;
2768 ++v)
2769 result.first = std::max(
2770 result.first,
2771 dof_info[0]
2772 .cell_active_fe_index[face_info.faces[macro_face].cells_interior[v]]);
2773 if (face_info.faces[macro_face].cells_exterior[0] !=
2774 numbers::invalid_unsigned_int)
2775 for (unsigned int v = 0; v < VectorizedArrayType::size() &&
2776 face_info.faces[macro_face].cells_exterior[v] !=
2777 numbers::invalid_unsigned_int;
2778 ++v)
2779 result.second = std::max(
2780 result.first,
2781 dof_info[0]
2782 .cell_active_fe_index[face_info.faces[macro_face].cells_exterior[v]]);
2783 else
2784 result.second = numbers::invalid_unsigned_int;
2785 return result;
2786}
2787
2788
2789
2790template <int dim, typename Number, typename VectorizedArrayType>
2791inline bool
2792MatrixFree<dim, Number, VectorizedArrayType>::indices_initialized() const
2793{
2794 return indices_are_initialized;
2795}
2796
2797
2798
2799template <int dim, typename Number, typename VectorizedArrayType>
2800inline bool
2801MatrixFree<dim, Number, VectorizedArrayType>::mapping_initialized() const
2802{
2803 return mapping_is_initialized;
2804}
2805
2806
2807template <int dim, typename Number, typename VectorizedArrayType>
2808inline unsigned int
2809MatrixFree<dim, Number, VectorizedArrayType>::get_mg_level() const
2810{
2811 return mg_level;
2812}
2813
2814
2815
2816template <int dim, typename Number, typename VectorizedArrayType>
2817AlignedVector<VectorizedArrayType> *
2818MatrixFree<dim, Number, VectorizedArrayType>::acquire_scratch_data() const
2819{
2820 using list_type =
2821 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
2822 list_type &data = scratch_pad.get();
2823 for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
2824 if (it->first == false)
2825 {
2826 it->first = true;
2827 return &it->second;
2828 }
2829 data.emplace_front(true, AlignedVector<VectorizedArrayType>());
2830 return &data.front().second;
2831}
2832
2833
2834
2835template <int dim, typename Number, typename VectorizedArrayType>
2836void
2837MatrixFree<dim, Number, VectorizedArrayType>::release_scratch_data(
2838 const AlignedVector<VectorizedArrayType> *scratch) const
2839{
2840 using list_type =
2841 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
2842 list_type &data = scratch_pad.get();
2843 for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
2844 if (&it->second == scratch)
2845 {
2846 Assert(it->first == true, ExcInternalError());
2847 it->first = false;
2848 return;
2849 }
2850 AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2851}
2852
2853
2854
2855template <int dim, typename Number, typename VectorizedArrayType>
2856AlignedVector<Number> *
2857MatrixFree<dim, Number, VectorizedArrayType>::
2858 acquire_scratch_data_non_threadsafe() const
2859{
2860 for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
2861 scratch_pad_non_threadsafe.begin();
2862 it != scratch_pad_non_threadsafe.end();
2863 ++it)
2864 if (it->first == false)
2865 {
2866 it->first = true;
2867 return &it->second;
2868 }
2869 scratch_pad_non_threadsafe.push_front(
2870 std::make_pair(true, AlignedVector<Number>()));
2871 return &scratch_pad_non_threadsafe.front().second;
2872}
2873
2874
2875
2876template <int dim, typename Number, typename VectorizedArrayType>
2877void
2878MatrixFree<dim, Number, VectorizedArrayType>::
2879 release_scratch_data_non_threadsafe(
2880 const AlignedVector<Number> *scratch) const
2881{
2882 for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
2883 scratch_pad_non_threadsafe.begin();
2884 it != scratch_pad_non_threadsafe.end();
2885 ++it)
2886 if (&it->second == scratch)
2887 {
2888 Assert(it->first == true, ExcInternalError());
2889 it->first = false;
2890 return;
2891 }
2892 AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
2893}
2894
2895
2896
2897// ------------------------------ reinit functions ---------------------------
2898
2899namespace internal
2900{
2901 namespace MatrixFreeImplementation
2902 {
2903 template <int dim, int spacedim>
2904 inline std::vector<IndexSet>
2905 extract_locally_owned_index_sets(
2906 const std::vector<const ::DoFHandler<dim, spacedim> *> &dofh,
2907 const unsigned int level)
2908 {
2909 std::vector<IndexSet> locally_owned_set;
2910 locally_owned_set.reserve(dofh.size());
2911 for (unsigned int j = 0; j < dofh.size(); j++)
2912 if (level == numbers::invalid_unsigned_int)
2913 locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
2914 else
2915 locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
2916 return locally_owned_set;
2917 }
2918 } // namespace MatrixFreeImplementation
2919} // namespace internal
2920
2921
2922
2923template <int dim, typename Number, typename VectorizedArrayType>
2924template <typename QuadratureType, typename number2>
2925void
2926MatrixFree<dim, Number, VectorizedArrayType>::reinit(
2927 const DoFHandler<dim> & dof_handler,
2928 const AffineConstraints<number2> &constraints_in,
2929 const QuadratureType & quad,
2930 const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
2931 &additional_data)
2932{
2933 std::vector<const DoFHandler<dim, dim> *> dof_handlers;
2934 std::vector<const AffineConstraints<number2> *> constraints;
2935 std::vector<QuadratureType> quads;
2936
2937 dof_handlers.push_back(&dof_handler);
2938 constraints.push_back(&constraints_in);
2939 quads.push_back(quad);
2940
2941 std::vector<IndexSet> locally_owned_sets =
2942 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2943 dof_handlers, additional_data.mg_level);
2944
2945 std::vector<hp::QCollection<dim>> quad_hp;
2946 quad_hp.emplace_back(quad);
2947
2948 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(
2949 StaticMappingQ1<dim>::mapping),
2950 dof_handlers,
2951 constraints,
2952 locally_owned_sets,
2953 quad_hp,
2954 additional_data);
2955}
2956
2957
2958
2959template <int dim, typename Number, typename VectorizedArrayType>
2960template <typename QuadratureType, typename number2, typename MappingType>
2961void
2962MatrixFree<dim, Number, VectorizedArrayType>::reinit(
2963 const MappingType & mapping,
2964 const DoFHandler<dim> & dof_handler,
2965 const AffineConstraints<number2> &constraints_in,
2966 const QuadratureType & quad,
2967 const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
2968 &additional_data)
2969{
2970 std::vector<const DoFHandler<dim, dim> *> dof_handlers;
2971 std::vector<const AffineConstraints<number2> *> constraints;
2972
2973 dof_handlers.push_back(&dof_handler);
2974 constraints.push_back(&constraints_in);
2975
2976 std::vector<IndexSet> locally_owned_sets =
2977 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
2978 dof_handlers, additional_data.mg_level);
2979
2980 std::vector<hp::QCollection<dim>> quad_hp;
2981 quad_hp.emplace_back(quad);
2982
2983 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
2984 dof_handlers,
2985 constraints,
2986 locally_owned_sets,
2987 quad_hp,
2988 additional_data);
2989}
2990
2991
2992
2993template <int dim, typename Number, typename VectorizedArrayType>
2994template <typename QuadratureType, typename number2>
2995void
2996MatrixFree<dim, Number, VectorizedArrayType>::reinit(
2997 const std::vector<const DoFHandler<dim> *> & dof_handler,
2998 const std::vector<const AffineConstraints<number2> *> &constraint,
2999 const std::vector<QuadratureType> & quad,
3000 const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
3001 &additional_data)
3002{
3003 std::vector<IndexSet> locally_owned_set =
3004 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3005 dof_handler, additional_data.mg_level);
3006 std::vector<hp::QCollection<dim>> quad_hp;
3007 for (unsigned int q = 0; q < quad.size(); ++q)
3008 quad_hp.emplace_back(quad[q]);
3009
3010 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(
3011 StaticMappingQ1<dim>::mapping),
3012 dof_handler,
3013 constraint,
3014 locally_owned_set,
3015 quad_hp,
3016 additional_data);
3017}
3018
3019
3020
3021template <int dim, typename Number, typename VectorizedArrayType>
3022template <typename QuadratureType,
3023 typename number2,
3024 typename DoFHandlerType,
3025 typename MappingType>
3026void
3027MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3028 const MappingType & mapping,
3029 const std::vector<const DoFHandlerType *> & dof_handler,
3030 const std::vector<const AffineConstraints<number2> *> &constraint,
3031 const std::vector<QuadratureType> & quad,
3032 const AdditionalData & additional_data)
3033{
3034 static_assert(dim == DoFHandlerType::dimension,
3035 "Dimension dim not equal to DoFHandlerType::dimension.");
3036
3037 std::vector<const DoFHandler<dim> *> dof_handlers;
3038
3039 for (const auto dh : dof_handler)
3040 dof_handlers.push_back(dh);
3041
3042 this->reinit(mapping, dof_handlers, constraint, quad, additional_data);
3043}
3044
3045
3046
3047template <int dim, typename Number, typename VectorizedArrayType>
3048template <typename QuadratureType, typename number2>
3049void
3050MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3051 const std::vector<const DoFHandler<dim> *> & dof_handler,
3052 const std::vector<const AffineConstraints<number2> *> &constraint,
3053 const QuadratureType & quad,
3054 const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
3055 &additional_data)
3056{
3057 std::vector<IndexSet> locally_owned_set =
3058 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3059 dof_handler, additional_data.mg_level);
3060 std::vector<hp::QCollection<dim>> quad_hp;
3061 quad_hp.emplace_back(quad);
3062
3063 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(
3064 StaticMappingQ1<dim>::mapping),
3065 dof_handler,
3066 constraint,
3067 locally_owned_set,
3068 quad_hp,
3069 additional_data);
3070}
3071
3072
3073
3074template <int dim, typename Number, typename VectorizedArrayType>
3075template <typename QuadratureType, typename number2, typename DoFHandlerType>
3076void
3077MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3078 const std::vector<const DoFHandlerType *> & dof_handler,
3079 const std::vector<const AffineConstraints<number2> *> &constraint,
3080 const std::vector<QuadratureType> & quad,
3081 const AdditionalData & additional_data)
3082{
3083 static_assert(dim == DoFHandlerType::dimension,
3084 "Dimension dim not equal to DoFHandlerType::dimension.");
3085
3086 std::vector<const DoFHandler<dim> *> dof_handlers;
3087
3088 for (const auto dh : dof_handler)
3089 dof_handlers.push_back(dof_handler);
3090
3091 this->reinit(dof_handlers, constraint, quad, additional_data);
3092}
3093
3094
3095
3096template <int dim, typename Number, typename VectorizedArrayType>
3097template <typename QuadratureType, typename number2, typename MappingType>
3098void
3099MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3100 const MappingType & mapping,
3101 const std::vector<const DoFHandler<dim> *> & dof_handler,
3102 const std::vector<const AffineConstraints<number2> *> &constraint,
3103 const QuadratureType & quad,
3104 const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
3105 &additional_data)
3106{
3107 std::vector<IndexSet> locally_owned_set =
3108 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3109 dof_handler, additional_data.mg_level);
3110 std::vector<hp::QCollection<dim>> quad_hp;
3111 quad_hp.emplace_back(quad);
3112
3113 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3114 dof_handler,
3115 constraint,
3116 locally_owned_set,
3117 quad_hp,
3118 additional_data);
3119}
3120
3121
3122
3123template <int dim, typename Number, typename VectorizedArrayType>
3124template <typename QuadratureType, typename number2, typename MappingType>
3125void
3126MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3127 const MappingType & mapping,
3128 const std::vector<const DoFHandler<dim> *> & dof_handler,
3129 const std::vector<const AffineConstraints<number2> *> &constraint,
3130 const std::vector<QuadratureType> & quad,
3131 const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
3132 &additional_data)
3133{
3134 std::vector<IndexSet> locally_owned_set =
3135 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3136 dof_handler, additional_data.mg_level);
3137 std::vector<hp::QCollection<dim>> quad_hp;
3138 for (unsigned int q = 0; q < quad.size(); ++q)
3139 quad_hp.emplace_back(quad[q]);
3140
3141 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3142 dof_handler,
3143 constraint,
3144 locally_owned_set,
3145 quad_hp,
3146 additional_data);
3147}
3148
3149
3150
3151template <int dim, typename Number, typename VectorizedArrayType>
3152template <typename QuadratureType,
3153 typename number2,
3154 typename DoFHandlerType,
3155 typename MappingType>
3156void
3157MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3158 const MappingType & mapping,
3159 const std::vector<const DoFHandlerType *> & dof_handler,
3160 const std::vector<const AffineConstraints<number2> *> &constraint,
3161 const QuadratureType & quad,
3162 const AdditionalData & additional_data)
3163{
3164 static_assert(dim == DoFHandlerType::dimension,
3165 "Dimension dim not equal to DoFHandlerType::dimension.");
3166
3167 std::vector<const DoFHandler<dim> *> dof_handlers;
3168
3169 for (const auto dh : dof_handler)
3170 dof_handlers.push_back(dof_handler);
3171
3172 this->reinit(mapping, dof_handlers, constraint, quad, additional_data);
3173}
3174
3175
3176
3177template <int dim, typename Number, typename VectorizedArrayType>
3178template <typename QuadratureType, typename number2, typename DoFHandlerType>
3179void
3180MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3181 const std::vector<const DoFHandlerType *> & dof_handler,
3182 const std::vector<const AffineConstraints<number2> *> &constraint,
3183 const QuadratureType & quad,
3184 const AdditionalData & additional_data)
3185{
3186 static_assert(dim == DoFHandlerType::dimension,
3187 "Dimension dim not equal to DoFHandlerType::dimension.");
3188
3189 std::vector<const DoFHandler<dim> *> dof_handlers;
3190
3191 for (const auto dh : dof_handler)
3192 dof_handlers.push_back(dof_handler);
3193
3194 this->reinit(dof_handlers, constraint, quad, additional_data);
3195}
3196
3197
3198
3199// ------------------------------ implementation of loops --------------------
3200
3201// internal helper functions that define how to call MPI data exchange
3202// functions: for generic vectors, do nothing at all. For distributed vectors,
3203// call update_ghost_values_start function and so on. If we have collections
3204// of vectors, just do the individual functions of the components. In order to
3205// keep ghost values consistent (whether we are in read or write mode), we
3206// also reset the values at the end. the whole situation is a bit complicated
3207// by the fact that we need to treat block vectors differently, which use some
3208// additional helper functions to select the blocks and template magic.
3209namespace internal
3210{
3214 template <int dim, typename Number, typename VectorizedArrayType>
3215 struct VectorDataExchange
3216 {
3217 // A shift for the MPI messages to reduce the risk for accidental
3218 // interaction with other open communications that a user program might
3219 // set up (parallel vectors support unfinished communication). We let
3220 // the other vectors use the first 20 assigned numbers and start the
3221 // matrix-free communication.
3222 static constexpr unsigned int channel_shift = 20;
3223
3224
3225
3230 VectorDataExchange(
3231 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3232 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3233 DataAccessOnFaces vector_face_access,
3234 const unsigned int n_components)
3235 : matrix_free(matrix_free)
3236 , vector_face_access(
3237 matrix_free.get_task_info().face_partition_data.empty() ?
3238 ::MatrixFree<dim, Number, VectorizedArrayType>::
3239 DataAccessOnFaces::unspecified :
3240 vector_face_access)
3241 , ghosts_were_set(false)
3242# ifdef DEAL_II_WITH_MPI
3243 , tmp_data(n_components)
3244 , requests(n_components)
3245# endif
3246 {
3247 (void)n_components;
3248 if (this->vector_face_access !=
3249 ::MatrixFree<dim, Number, VectorizedArrayType>::
3250 DataAccessOnFaces::unspecified)
3251 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3252 AssertDimension(
3253 matrix_free.get_dof_info(c).vector_exchanger_face_variants.size(),
3254 5);
3255 }
3256
3257
3258
3262 ~VectorDataExchange() // NOLINT
3263 {
3264# ifdef DEAL_II_WITH_MPI
3265 for (unsigned int i = 0; i < tmp_data.size(); ++i)
3266 if (tmp_data[i] != nullptr)
3267 matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3268# endif
3269 }
3270
3271
3272
3277 template <typename VectorType>
3278 unsigned int
3279 find_vector_in_mf(const VectorType &vec,
3280 const bool check_global_compatibility = true) const
3281 {
3282 // case 1: vector was set up with MatrixFree::initialize_dof_vector()
3283 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3284 if (vec.get_partitioner().get() ==
3285 matrix_free.get_dof_info(c).vector_partitioner.get())
3286 return c;
3287
3288 // case 2: user provided own partitioner (compatibility mode)
3289 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3290 if (check_global_compatibility ?
3291 vec.get_partitioner()->is_globally_compatible(
3292 *matrix_free.get_dof_info(c).vector_partitioner) :
3293 vec.get_partitioner()->is_compatible(
3294 *matrix_free.get_dof_info(c).vector_partitioner))
3295 return c;
3296
3297 Assert(false,
3298 ExcNotImplemented("Could not find partitioner that fits vector"));
3299
3300 return numbers::invalid_unsigned_int;
3301 }
3302
3303
3304
3309 const internal::MatrixFreeFunctions::VectorDataExchange::Base &
3310 get_partitioner(const unsigned int mf_component) const
3311 {
3312 AssertDimension(matrix_free.get_dof_info(mf_component)
3313 .vector_exchanger_face_variants.size(),
3314 5);
3315 if (vector_face_access ==
3316 ::MatrixFree<dim, Number, VectorizedArrayType>::
3317 DataAccessOnFaces::none)
3318 return *matrix_free.get_dof_info(mf_component)
3319 .vector_exchanger_face_variants[0];
3320 else if (vector_face_access ==
3321 ::MatrixFree<dim, Number, VectorizedArrayType>::
3322 DataAccessOnFaces::values)
3323 return *matrix_free.get_dof_info(mf_component)
3324 .vector_exchanger_face_variants[1];
3325 else if (vector_face_access ==
3326 ::MatrixFree<dim, Number, VectorizedArrayType>::
3327 DataAccessOnFaces::gradients)
3328 return *matrix_free.get_dof_info(mf_component)
3329 .vector_exchanger_face_variants[2];
3330 else if (vector_face_access ==
3331 ::MatrixFree<dim, Number, VectorizedArrayType>::
3332 DataAccessOnFaces::values_all_faces)
3333 return *matrix_free.get_dof_info(mf_component)
3334 .vector_exchanger_face_variants[3];
3335 else if (vector_face_access ==
3336 ::MatrixFree<dim, Number, VectorizedArrayType>::
3337 DataAccessOnFaces::gradients_all_faces)
3338 return *matrix_free.get_dof_info(mf_component)
3339 .vector_exchanger_face_variants[4];
3340 else
3341 return *matrix_free.get_dof_info(mf_component).vector_exchanger.get();
3342 }
3343
3344
3345
3349 template <typename VectorType,
3350 typename std::enable_if<is_serial_or_dummy<VectorType>::value,
3351 VectorType>::type * = nullptr>
3352 void
3353 update_ghost_values_start(const unsigned int /*component_in_block_vector*/,
3354 const VectorType & /*vec*/)
3355 {}
3356
3357
3362 template <typename VectorType,
3363 typename std::enable_if<
3364 !has_update_ghost_values_start<VectorType>::value &&
3365 !is_serial_or_dummy<VectorType>::value,
3366 VectorType>::type * = nullptr>
3367 void
3368 update_ghost_values_start(const unsigned int component_in_block_vector,
3369 const VectorType & vec)
3370 {
3371 (void)component_in_block_vector;
3372 bool ghosts_set = vec.has_ghost_elements();
3373 if (ghosts_set)
3374 ghosts_were_set = true;
3375
3376 vec.update_ghost_values();
3377 }
3378
3379
3380
3386 template <typename VectorType,
3387 typename std::enable_if<
3388 has_update_ghost_values_start<VectorType>::value &&
3389 !has_exchange_on_subset<VectorType>::value,
3390 VectorType>::type * = nullptr>
3391 void
3392 update_ghost_values_start(const unsigned int component_in_block_vector,
3393 const VectorType & vec)
3394 {
3395 (void)component_in_block_vector;
3396 bool ghosts_set = vec.has_ghost_elements();
3397 if (ghosts_set)
3398 ghosts_were_set = true;
3399
3400 vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3401 }
3402
3403
3404
3411 template <typename VectorType,
3412 typename std::enable_if<
3413 has_update_ghost_values_start<VectorType>::value &&
3414 has_exchange_on_subset<VectorType>::value,
3415 VectorType>::type * = nullptr>
3416 void
3417 update_ghost_values_start(const unsigned int component_in_block_vector,
3418 const VectorType & vec)
3419 {
3420 static_assert(
3421 std::is_same<Number, typename VectorType::value_type>::value,
3422 "Type mismatch between VectorType and VectorDataExchange");
3423 (void)component_in_block_vector;
3424 bool ghosts_set = vec.has_ghost_elements();
3425 if (ghosts_set)
3426 ghosts_were_set = true;
3427
3428 if (vec.size() != 0)
3429 {
3430# ifdef DEAL_II_WITH_MPI
3431 const unsigned int mf_component = find_vector_in_mf(vec);
3432
3433 const auto &part = get_partitioner(mf_component);
3434
3435 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3436 part.n_import_sm_procs() == 0)
3437 return;
3438
3439 tmp_data[component_in_block_vector] =
3440 matrix_free.acquire_scratch_data_non_threadsafe();
3441 tmp_data[component_in_block_vector]->resize_fast(
3442 part.n_import_indices());
3443 AssertDimension(requests.size(), tmp_data.size());
3444
3445 part.export_to_ghosted_array_start(
3446 component_in_block_vector * 2 + channel_shift,
3447 ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3448 vec.shared_vector_data(),
3449 ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3450 part.locally_owned_size(),
3451 matrix_free.get_dof_info(mf_component)
3452 .vector_partitioner->n_ghost_indices()),
3453 ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3454 part.n_import_indices()),
3455 this->requests[component_in_block_vector]);
3456# endif
3457 }
3458 }
3459
3460
3461
3466 template <
3467 typename VectorType,
3468 typename std::enable_if<!has_update_ghost_values_start<VectorType>::value,
3469 VectorType>::type * = nullptr>
3470 void
3471 update_ghost_values_finish(const unsigned int /*component_in_block_vector*/,
3472 const VectorType & /*vec*/)
3473 {}
3474
3475
3476
3482 template <typename VectorType,
3483 typename std::enable_if<
3484 has_update_ghost_values_start<VectorType>::value &&
3485 !has_exchange_on_subset<VectorType>::value,
3486 VectorType>::type * = nullptr>
3487 void
3488 update_ghost_values_finish(const unsigned int component_in_block_vector,
3489 const VectorType & vec)
3490 {
3491 (void)component_in_block_vector;
3492 vec.update_ghost_values_finish();
3493 }
3494
3495
3496
3503 template <typename VectorType,
3504 typename std::enable_if<
3505 has_update_ghost_values_start<VectorType>::value &&
3506 has_exchange_on_subset<VectorType>::value,
3507 VectorType>::type * = nullptr>
3508 void
3509 update_ghost_values_finish(const unsigned int component_in_block_vector,
3510 const VectorType & vec)
3511 {
3512 static_assert(
3513 std::is_same<Number, typename VectorType::value_type>::value,
3514 "Type mismatch between VectorType and VectorDataExchange");
3515 (void)component_in_block_vector;
3516
3517 if (vec.size() != 0)
3518 {
3519# ifdef DEAL_II_WITH_MPI
3520 AssertIndexRange(component_in_block_vector, tmp_data.size());
3521 AssertDimension(requests.size(), tmp_data.size());
3522
3523 const unsigned int mf_component = find_vector_in_mf(vec);
3524
3525 const auto &part = get_partitioner(mf_component);
3526
3527 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3528 part.n_import_sm_procs() != 0)
3529 {
3530 part.export_to_ghosted_array_finish(
3531 ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3532 vec.shared_vector_data(),
3533 ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3534 part.locally_owned_size(),
3535 matrix_free.get_dof_info(mf_component)
3536 .vector_partitioner->n_ghost_indices()),
3537 this->requests[component_in_block_vector]);
3538
3539 matrix_free.release_scratch_data_non_threadsafe(
3540 tmp_data[component_in_block_vector]);
3541 tmp_data[component_in_block_vector] = nullptr;
3542 }
3543# endif
3544 }
3545 // let vector know that ghosts are being updated and we can read from
3546 // them
3547 vec.set_ghost_state(true);
3548 }
3549
3550
3551
3555 template <typename VectorType,
3556 typename std::enable_if<is_serial_or_dummy<VectorType>::value,
3557 VectorType>::type * = nullptr>
3558 void
3559 compress_start(const unsigned int /*component_in_block_vector*/,
3560 VectorType & /*vec*/)
3561 {}
3562
3563
3564
3569 template <typename VectorType,
3570 typename std::enable_if<!has_compress_start<VectorType>::value &&
3571 !is_serial_or_dummy<VectorType>::value,
3572 VectorType>::type * = nullptr>
3573 void
3574 compress_start(const unsigned int component_in_block_vector,
3575 VectorType & vec)
3576 {
3577 (void)component_in_block_vector;
3578 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3579 vec.compress(::VectorOperation::add);
3580 }
3581
3582
3583
3589 template <
3590 typename VectorType,
3591 typename std::enable_if<has_compress_start<VectorType>::value &&
3592 !has_exchange_on_subset<VectorType>::value,
3593 VectorType>::type * = nullptr>
3594 void
3595 compress_start(const unsigned int component_in_block_vector,
3596 VectorType & vec)
3597 {
3598 (void)component_in_block_vector;
3599 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3600 vec.compress_start(component_in_block_vector + channel_shift);
3601 }
3602
3603
3604
3611 template <
3612 typename VectorType,
3613 typename std::enable_if<has_compress_start<VectorType>::value &&
3614 has_exchange_on_subset<VectorType>::value,
3615 VectorType>::type * = nullptr>
3616 void
3617 compress_start(const unsigned int component_in_block_vector,
3618 VectorType & vec)
3619 {
3620 static_assert(
3621 std::is_same<Number, typename VectorType::value_type>::value,
3622 "Type mismatch between VectorType and VectorDataExchange");
3623 (void)component_in_block_vector;
3624 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3625
3626 if (vec.size() != 0)
3627 {
3628# ifdef DEAL_II_WITH_MPI
3629 const unsigned int mf_component = find_vector_in_mf(vec);
3630
3631 const auto &part = get_partitioner(mf_component);
3632
3633 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3634 part.n_import_sm_procs() == 0)
3635 return;
3636
3637 tmp_data[component_in_block_vector] =
3638 matrix_free.acquire_scratch_data_non_threadsafe();
3639 tmp_data[component_in_block_vector]->resize_fast(
3640 part.n_import_indices());
3641 AssertDimension(requests.size(), tmp_data.size());
3642
3643 part.import_from_ghosted_array_start(
3644 ::VectorOperation::add,
3645 component_in_block_vector * 2 + channel_shift,
3646 ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3647 vec.shared_vector_data(),
3648 ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3649 matrix_free.get_dof_info(mf_component)
3650 .vector_partitioner->n_ghost_indices()),
3651 ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3652 part.n_import_indices()),
3653 this->requests[component_in_block_vector]);
3654# endif
3655 }
3656 }
3657
3658
3659
3664 template <typename VectorType,
3665 typename std::enable_if<!has_compress_start<VectorType>::value,
3666 VectorType>::type * = nullptr>
3667 void
3668 compress_finish(const unsigned int /*component_in_block_vector*/,
3669 VectorType & /*vec*/)
3670 {}
3671
3672
3673
3679 template <
3680 typename VectorType,
3681 typename std::enable_if<has_compress_start<VectorType>::value &&
3682 !has_exchange_on_subset<VectorType>::value,
3683 VectorType>::type * = nullptr>
3684 void
3685 compress_finish(const unsigned int component_in_block_vector,
3686 VectorType & vec)
3687 {
3688 (void)component_in_block_vector;
3689 vec.compress_finish(::VectorOperation::add);
3690 }
3691
3692
3693
3700 template <
3701 typename VectorType,
3702 typename std::enable_if<has_compress_start<VectorType>::value &&
3703 has_exchange_on_subset<VectorType>::value,
3704 VectorType>::type * = nullptr>
3705 void
3706 compress_finish(const unsigned int component_in_block_vector,
3707 VectorType & vec)
3708 {
3709 static_assert(
3710 std::is_same<Number, typename VectorType::value_type>::value,
3711 "Type mismatch between VectorType and VectorDataExchange");
3712 (void)component_in_block_vector;
3713 if (vec.size() != 0)
3714 {
3715# ifdef DEAL_II_WITH_MPI
3716 AssertIndexRange(component_in_block_vector, tmp_data.size());
3717 AssertDimension(requests.size(), tmp_data.size());
3718
3719 const unsigned int mf_component = find_vector_in_mf(vec);
3720
3721 const auto &part = get_partitioner(mf_component);
3722
3723 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3724 part.n_import_sm_procs() != 0)
3725 {
3726 part.import_from_ghosted_array_finish(
3727 VectorOperation::add,
3728 ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3729 vec.shared_vector_data(),
3730 ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3731 matrix_free.get_dof_info(mf_component)
3732 .vector_partitioner->n_ghost_indices()),
3733 ArrayView<const Number>(
3734 tmp_data[component_in_block_vector]->begin(),
3735 part.n_import_indices()),
3736 this->requests[component_in_block_vector]);
3737
3738 matrix_free.release_scratch_data_non_threadsafe(
3739 tmp_data[component_in_block_vector]);
3740 tmp_data[component_in_block_vector] = nullptr;
3741 }
3742
3743 if (Utilities::MPI::job_supports_mpi())
3744 MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3745# endif
3746 }
3747 }
3748
3749
3750
3754 template <typename VectorType,
3755 typename std::enable_if<is_serial_or_dummy<VectorType>::value,
3756 VectorType>::type * = nullptr>
3757 void
3758 reset_ghost_values(const VectorType & /*vec*/) const
3759 {}
3760
3761
3762
3767 template <
3768 typename VectorType,
3769 typename std::enable_if<!has_exchange_on_subset<VectorType>::value &&
3770 !is_serial_or_dummy<VectorType>::value,
3771 VectorType>::type * = nullptr>
3772 void
3773 reset_ghost_values(const VectorType &vec) const
3774 {
3775 if (ghosts_were_set == true)
3776 return;
3777
3778 vec.zero_out_ghost_values();
3779 }
3780
3781
3782
3788 template <typename VectorType,
3789 typename std::enable_if<has_exchange_on_subset<VectorType>::value,
3790 VectorType>::type * = nullptr>
3791 void
3792 reset_ghost_values(const VectorType &vec) const
3793 {
3794 static_assert(
3795 std::is_same<Number, typename VectorType::value_type>::value,
3796 "Type mismatch between VectorType and VectorDataExchange");
3797 if (ghosts_were_set == true)
3798 return;
3799
3800 if (vec.size() != 0)
3801 {
3802# ifdef DEAL_II_WITH_MPI
3803 AssertDimension(requests.size(), tmp_data.size());
3804
3805 const unsigned int mf_component = find_vector_in_mf(vec);
3806
3807 const auto &part = get_partitioner(mf_component);
3808
3809 if (part.n_ghost_indices() > 0)
3810 {
3811 part.reset_ghost_values(ArrayView<Number>(
3812 const_cast<LinearAlgebra::distributed::Vector<Number> &>(vec)
3813 .begin() +
3814 part.locally_owned_size(),
3815 matrix_free.get_dof_info(mf_component)
3816 .vector_partitioner->n_ghost_indices()));
3817 }
3818
3819# endif
3820 }
3821 // let vector know that it's not ghosted anymore
3822 vec.set_ghost_state(false);
3823 }
3824
3825
3826
3832 template <typename VectorType,
3833 typename std::enable_if<has_exchange_on_subset<VectorType>::value,
3834 VectorType>::type * = nullptr>
3835 void
3836 zero_vector_region(const unsigned int range_index, VectorType &vec) const
3837 {
3838 static_assert(
3839 std::is_same<Number, typename VectorType::value_type>::value,
3840 "Type mismatch between VectorType and VectorDataExchange");
3841 if (range_index == numbers::invalid_unsigned_int)
3842 vec = Number();
3843 else
3844 {
3845 const unsigned int mf_component = find_vector_in_mf(vec, false);
3846 const internal::MatrixFreeFunctions::DoFInfo &dof_info =
3847 matrix_free.get_dof_info(mf_component);
3848 Assert(dof_info.vector_zero_range_list_index.empty() == false,
3849 ExcNotInitialized());
3850
3851 Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3852 ExcInternalError());
3853 AssertIndexRange(range_index,
3854 dof_info.vector_zero_range_list_index.size() - 1);
3855 for (unsigned int id =
3856 dof_info.vector_zero_range_list_index[range_index];
3857 id != dof_info.vector_zero_range_list_index[range_index + 1];
3858 ++id)
3859 std::memset(vec.begin() + dof_info.vector_zero_range_list[id].first,
3860 0,
3861 (dof_info.vector_zero_range_list[id].second -
3862 dof_info.vector_zero_range_list[id].first) *
3863 sizeof(Number));
3864 }
3865 }
3866
3867
3868
3874 template <
3875 typename VectorType,
3876 typename std::enable_if<!has_exchange_on_subset<VectorType>::value,
3877 VectorType>::type * = nullptr,
3878 typename VectorType::value_type * = nullptr>
3879 void
3880 zero_vector_region(const unsigned int range_index, VectorType &vec) const
3881 {
3882 if (range_index == numbers::invalid_unsigned_int || range_index == 0)
3883 vec = typename VectorType::value_type();
3884 }
3885
3886
3887
3892 void
3893 zero_vector_region(const unsigned int, ...) const
3894 {
3895 Assert(false,
3896 ExcNotImplemented("Zeroing is only implemented for vector types "
3897 "which provide VectorType::value_type"));
3898 }
3899
3900
3901
3902 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
3903 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3904 DataAccessOnFaces vector_face_access;
3905 bool ghosts_were_set;
3906# ifdef DEAL_II_WITH_MPI
3907 std::vector<AlignedVector<Number> *> tmp_data;
3908 std::vector<std::vector<MPI_Request>> requests;
3909# endif
3910 }; // VectorDataExchange
3911
3912 template <typename VectorStruct>
3913 unsigned int
3914 n_components(const VectorStruct &vec);
3915
3916 template <typename VectorStruct>
3917 unsigned int
3918 n_components_block(const VectorStruct &vec,
3919 std::integral_constant<bool, true>)
3920 {
3921 unsigned int components = 0;
3922 for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
3923 components += n_components(vec.block(bl));
3924 return components;
3925 }
3926
3927 template <typename VectorStruct>
3928 unsigned int
3929 n_components_block(const VectorStruct &, std::integral_constant<bool, false>)
3930 {
3931 return 1;
3932 }
3933
3934 template <typename VectorStruct>
3935 unsigned int
3936 n_components(const VectorStruct &vec)
3937 {
3938 return n_components_block(
3939 vec, std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3940 }
3941
3942 template <typename VectorStruct>
3943 inline unsigned int
3944 n_components(const std::vector<VectorStruct> &vec)
3945 {
3946 unsigned int components = 0;
3947 for (unsigned int comp = 0; comp < vec.size(); comp++)
3948 components += n_components_block(
3949 vec[comp],
3950 std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3951 return components;
3952 }
3953
3954 template <typename VectorStruct>
3955 inline unsigned int
3956 n_components(const std::vector<VectorStruct *> &vec)
3957 {
3958 unsigned int components = 0;
3959 for (unsigned int comp = 0; comp < vec.size(); comp++)
3960 components += n_components_block(
3961 *vec[comp],
3962 std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3963 return components;
3964 }
3965
3966
3967
3968 // A helper function to identify block vectors with many components where we
3969 // should not try to overlap computations and communication because there
3970 // would be too many outstanding communication requests.
3971
3972 // default value for vectors that do not have communication_block_size
3973 template <
3974 typename VectorStruct,
3975 typename std::enable_if<!has_communication_block_size<VectorStruct>::value,
3976 VectorStruct>::type * = nullptr>
3977 constexpr unsigned int
3978 get_communication_block_size(const VectorStruct &)
3979 {
3980 return numbers::invalid_unsigned_int;
3981 }
3982
3983
3984
3985 template <
3986 typename VectorStruct,
3987 typename std::enable_if<has_communication_block_size<VectorStruct>::value,
3988 VectorStruct>::type * = nullptr>
3989 constexpr unsigned int
3990 get_communication_block_size(const VectorStruct &)
3991 {
3992 return VectorStruct::communication_block_size;
3993 }
3994
3995
3996
3997 // --------------------------------------------------------------------------
3998 // below we have wrappers to distinguish between block and non-block vectors.
3999 // --------------------------------------------------------------------------
4000
4001 //
4002 // update_ghost_values_start
4003 //
4004
4005 // update_ghost_values for block vectors
4006 template <int dim,
4007 typename VectorStruct,
4008 typename Number,
4009 typename VectorizedArrayType,
4010 typename std::enable_if<IsBlockVector<VectorStruct>::value,
4011 VectorStruct>::type * = nullptr>
4012 void
4013 update_ghost_values_start(
4014 const VectorStruct & vec,
4015 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4016 const unsigned int channel = 0)
4017 {
4018 if (get_communication_block_size(vec) < vec.n_blocks())
4019 {
4020 // don't forget to set ghosts_were_set, that otherwise happens
4021 // inside VectorDataExchange::update_ghost_values_start()
4022 exchanger.ghosts_were_set = vec.has_ghost_elements();
4023 vec.update_ghost_values();
4024 }
4025 else
4026 {
4027 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4028 update_ghost_values_start(vec.block(i), exchanger, channel + i);
4029 }
4030 }
4031
4032
4033
4034 // update_ghost_values for non-block vectors
4035 template <int dim,
4036 typename VectorStruct,
4037 typename Number,
4038 typename VectorizedArrayType,
4039 typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4040 VectorStruct>::type * = nullptr>
4041 void
4042 update_ghost_values_start(
4043 const VectorStruct & vec,
4044 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4045 const unsigned int channel = 0)
4046 {
4047 exchanger.update_ghost_values_start(channel, vec);
4048 }
4049
4050
4051
4052 // update_ghost_values_start() for vector of vectors
4053 template <int dim,
4054 typename VectorStruct,
4055 typename Number,
4056 typename VectorizedArrayType>
4057 inline void
4058 update_ghost_values_start(
4059 const std::vector<VectorStruct> & vec,
4060 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4061 {
4062 unsigned int component_index = 0;
4063 for (unsigned int comp = 0; comp < vec.size(); comp++)
4064 {
4065 update_ghost_values_start(vec[comp], exchanger, component_index);
4066 component_index += n_components(vec[comp]);
4067 }
4068 }
4069
4070
4071
4072 // update_ghost_values_start() for vector of pointers to vectors
4073 template <int dim,
4074 typename VectorStruct,
4075 typename Number,
4076 typename VectorizedArrayType>
4077 inline void
4078 update_ghost_values_start(
4079 const std::vector<VectorStruct *> & vec,
4080 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4081 {
4082 unsigned int component_index = 0;
4083 for (unsigned int comp = 0; comp < vec.size(); comp++)
4084 {
4085 update_ghost_values_start(*vec[comp], exchanger, component_index);
4086 component_index += n_components(*vec[comp]);
4087 }
4088 }
4089
4090
4091
4092 //
4093 // update_ghost_values_finish
4094 //
4095
4096 // for block vectors
4097 template <int dim,
4098 typename VectorStruct,
4099 typename Number,
4100 typename VectorizedArrayType,
4101 typename std::enable_if<IsBlockVector<VectorStruct>::value,
4102 VectorStruct>::type * = nullptr>
4103 void
4104 update_ghost_values_finish(
4105 const VectorStruct & vec,
4106 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4107 const unsigned int channel = 0)
4108 {
4109 if (get_communication_block_size(vec) < vec.n_blocks())
4110 {
4111 // do nothing, everything has already been completed in the _start()
4112 // call
4113 }
4114 else
4115 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4116 update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4117 }
4118
4119
4120
4121 // for non-block vectors
4122 template <int dim,
4123 typename VectorStruct,
4124 typename Number,
4125 typename VectorizedArrayType,
4126 typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4127 VectorStruct>::type * = nullptr>
4128 void
4129 update_ghost_values_finish(
4130 const VectorStruct & vec,
4131 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4132 const unsigned int channel = 0)
4133 {
4134 exchanger.update_ghost_values_finish(channel, vec);
4135 }
4136
4137
4138
4139 // for vector of vectors
4140 template <int dim,
4141 typename VectorStruct,
4142 typename Number,
4143 typename VectorizedArrayType>
4144 inline void
4145 update_ghost_values_finish(
4146 const std::vector<VectorStruct> & vec,
4147 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4148 {
4149 unsigned int component_index = 0;
4150 for (unsigned int comp = 0; comp < vec.size(); comp++)
4151 {
4152 update_ghost_values_finish(vec[comp], exchanger, component_index);
4153 component_index += n_components(vec[comp]);
4154 }
4155 }
4156
4157
4158
4159 // for vector of pointers to vectors
4160 template <int dim,
4161 typename VectorStruct,
4162 typename Number,
4163 typename VectorizedArrayType>
4164 inline void
4165 update_ghost_values_finish(
4166 const std::vector<VectorStruct *> & vec,
4167 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4168 {
4169 unsigned int component_index = 0;
4170 for (unsigned int comp = 0; comp < vec.size(); comp++)
4171 {
4172 update_ghost_values_finish(*vec[comp], exchanger, component_index);
4173 component_index += n_components(*vec[comp]);
4174 }
4175 }
4176
4177
4178
4179 //
4180 // compress_start
4181 //
4182
4183 // for block vectors
4184 template <int dim,
4185 typename VectorStruct,
4186 typename Number,
4187 typename VectorizedArrayType,
4188 typename std::enable_if<IsBlockVector<VectorStruct>::value,
4189 VectorStruct>::type * = nullptr>
4190 inline void
4191 compress_start(
4192 VectorStruct & vec,
4193 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4194 const unsigned int channel = 0)
4195 {
4196 if (get_communication_block_size(vec) < vec.n_blocks())
4197 vec.compress(::VectorOperation::add);
4198 else
4199 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4200 compress_start(vec.block(i), exchanger, channel + i);
4201 }
4202
4203
4204
4205 // for non-block vectors
4206 template <int dim,
4207 typename VectorStruct,
4208 typename Number,
4209 typename VectorizedArrayType,
4210 typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4211 VectorStruct>::type * = nullptr>
4212 inline void
4213 compress_start(
4214 VectorStruct & vec,
4215 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4216 const unsigned int channel = 0)
4217 {
4218 exchanger.compress_start(channel, vec);
4219 }
4220
4221
4222
4223 // for std::vector of vectors
4224 template <int dim,
4225 typename VectorStruct,
4226 typename Number,
4227 typename VectorizedArrayType>
4228 inline void
4229 compress_start(
4230 std::vector<VectorStruct> & vec,
4231 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4232 {
4233 unsigned int component_index = 0;
4234 for (unsigned int comp = 0; comp < vec.size(); comp++)
4235 {
4236 compress_start(vec[comp], exchanger, component_index);
4237 component_index += n_components(vec[comp]);
4238 }
4239 }
4240
4241
4242
4243 // for std::vector of pointer to vectors
4244 template <int dim,
4245 typename VectorStruct,
4246 typename Number,
4247 typename VectorizedArrayType>
4248 inline void
4249 compress_start(
4250 std::vector<VectorStruct *> & vec,
4251 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4252 {
4253 unsigned int component_index = 0;
4254 for (unsigned int comp = 0; comp < vec.size(); comp++)
4255 {
4256 compress_start(*vec[comp], exchanger, component_index);
4257 component_index += n_components(*vec[comp]);
4258 }
4259 }
4260
4261
4262
4263 //
4264 // compress_finish
4265 //
4266
4267 // for block vectors
4268 template <int dim,
4269 typename VectorStruct,
4270 typename Number,
4271 typename VectorizedArrayType,
4272 typename std::enable_if<IsBlockVector<VectorStruct>::value,
4273 VectorStruct>::type * = nullptr>
4274 inline void
4275 compress_finish(
4276 VectorStruct & vec,
4277 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4278 const unsigned int channel = 0)
4279 {
4280 if (get_communication_block_size(vec) < vec.n_blocks())
4281 {
4282 // do nothing, everything has already been completed in the _start()
4283 // call
4284 }
4285 else
4286 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4287 compress_finish(vec.block(i), exchanger, channel + i);
4288 }
4289
4290
4291
4292 // for non-block vectors
4293 template <int dim,
4294 typename VectorStruct,
4295 typename Number,
4296 typename VectorizedArrayType,
4297 typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4298 VectorStruct>::type * = nullptr>
4299 inline void
4300 compress_finish(
4301 VectorStruct & vec,
4302 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4303 const unsigned int channel = 0)
4304 {
4305 exchanger.compress_finish(channel, vec);
4306 }
4307
4308
4309
4310 // for std::vector of vectors
4311 template <int dim,
4312 typename VectorStruct,
4313 typename Number,
4314 typename VectorizedArrayType>
4315 inline void
4316 compress_finish(
4317 std::vector<VectorStruct> & vec,
4318 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4319 {
4320 unsigned int component_index = 0;
4321 for (unsigned int comp = 0; comp < vec.size(); comp++)
4322 {
4323 compress_finish(vec[comp], exchanger, component_index);
4324 component_index += n_components(vec[comp]);
4325 }
4326 }
4327
4328
4329
4330 // for std::vector of pointer to vectors
4331 template <int dim,
4332 typename VectorStruct,
4333 typename Number,
4334 typename VectorizedArrayType>
4335 inline void
4336 compress_finish(
4337 std::vector<VectorStruct *> & vec,
4338 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4339 {
4340 unsigned int component_index = 0;
4341 for (unsigned int comp = 0; comp < vec.size(); comp++)
4342 {
4343 compress_finish(*vec[comp], exchanger, component_index);
4344 component_index += n_components(*vec[comp]);
4345 }
4346 }
4347
4348
4349
4350 //
4351 // reset_ghost_values:
4352 //
4353 // if the input vector did not have ghosts imported, clear them here again
4354 // in order to avoid subsequent operations e.g. in linear solvers to work
4355 // with ghosts all the time
4356 //
4357
4358 // for block vectors
4359 template <int dim,
4360 typename VectorStruct,
4361 typename Number,
4362 typename VectorizedArrayType,
4363 typename std::enable_if<IsBlockVector<VectorStruct>::value,
4364 VectorStruct>::type * = nullptr>
4365 inline void
4366 reset_ghost_values(
4367 const VectorStruct & vec,
4368 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4369 {
4370 // return immediately if there is nothing to do.
4371 if (exchanger.ghosts_were_set == true)
4372 return;
4373
4374 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4375 reset_ghost_values(vec.block(i), exchanger);
4376 }
4377
4378
4379
4380 // for non-block vectors
4381 template <int dim,
4382 typename VectorStruct,
4383 typename Number,
4384 typename VectorizedArrayType,
4385 typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4386 VectorStruct>::type * = nullptr>
4387 inline void
4388 reset_ghost_values(
4389 const VectorStruct & vec,
4390 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4391 {
4392 exchanger.reset_ghost_values(vec);
4393 }
4394
4395
4396
4397 // for std::vector of vectors
4398 template <int dim,
4399 typename VectorStruct,
4400 typename Number,
4401 typename VectorizedArrayType>
4402 inline void
4403 reset_ghost_values(
4404 const std::vector<VectorStruct> & vec,
4405 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4406 {
4407 // return immediately if there is nothing to do.
4408 if (exchanger.ghosts_were_set == true)
4409 return;
4410
4411 for (unsigned int comp = 0; comp < vec.size(); comp++)
4412 reset_ghost_values(vec[comp], exchanger);
4413 }
4414
4415
4416
4417 // for std::vector of pointer to vectors
4418 template <int dim,
4419 typename VectorStruct,
4420 typename Number,
4421 typename VectorizedArrayType>
4422 inline void
4423 reset_ghost_values(
4424 const std::vector<VectorStruct *> & vec,
4425 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4426 {
4427 // return immediately if there is nothing to do.
4428 if (exchanger.ghosts_were_set == true)
4429 return;
4430
4431 for (unsigned int comp = 0; comp < vec.size(); comp++)
4432 reset_ghost_values(*vec[comp], exchanger);
4433 }
4434
4435
4436
4437 //
4438 // zero_vector_region
4439 //
4440
4441 // for block vectors
4442 template <int dim,
4443 typename VectorStruct,
4444 typename Number,
4445 typename VectorizedArrayType,
4446 typename std::enable_if<IsBlockVector<VectorStruct>::value,
4447 VectorStruct>::type * = nullptr>
4448 inline void
4449 zero_vector_region(
4450 const unsigned int range_index,
4451 VectorStruct & vec,
4452 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4453 {
4454 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4455 exchanger.zero_vector_region(range_index, vec.block(i));
4456 }
4457
4458
4459
4460 // for non-block vectors
4461 template <int dim,
4462 typename VectorStruct,
4463 typename Number,
4464 typename VectorizedArrayType,
4465 typename std::enable_if<!IsBlockVector<VectorStruct>::value,
4466 VectorStruct>::type * = nullptr>
4467 inline void
4468 zero_vector_region(
4469 const unsigned int range_index,
4470 VectorStruct & vec,
4471 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4472 {
4473 exchanger.zero_vector_region(range_index, vec);
4474 }
4475
4476
4477
4478 // for std::vector of vectors
4479 template <int dim,
4480 typename VectorStruct,
4481 typename Number,
4482 typename VectorizedArrayType>
4483 inline void
4484 zero_vector_region(
4485 const unsigned int range_index,
4486 std::vector<VectorStruct> & vec,
4487 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4488 {
4489 for (unsigned int comp = 0; comp < vec.size(); comp++)
4490 zero_vector_region(range_index, vec[comp], exchanger);
4491 }
4492
4493
4494
4495 // for std::vector of pointers to vectors
4496 template <int dim,
4497 typename VectorStruct,
4498 typename Number,
4499 typename VectorizedArrayType>
4500 inline void
4501 zero_vector_region(
4502 const unsigned int range_index,
4503 std::vector<VectorStruct *> & vec,
4504 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4505 {
4506 for (unsigned int comp = 0; comp < vec.size(); comp++)
4507 zero_vector_region(range_index, *vec[comp], exchanger);
4508 }
4509
4510
4511
4512 namespace MatrixFreeFunctions
4513 {
4514 // struct to select between a const interface and a non-const interface
4515 // for MFWorker
4516 template <typename, typename, typename, typename, bool>
4517 struct InterfaceSelector
4518 {};
4519
4520 // Version for constant functions
4521 template <typename MF,
4522 typename InVector,
4523 typename OutVector,
4524 typename Container>
4525 struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4526 {
4527 using function_type = void (Container::*)(
4528 const MF &,
4529 OutVector &,
4530 const InVector &,
4531 const std::pair<unsigned int, unsigned int> &) const;
4532 };
4533
4534 // Version for non-constant functions
4535 template <typename MF,
4536 typename InVector,
4537 typename OutVector,
4538 typename Container>
4539 struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4540 {
4541 using function_type =
4542 void (Container::*)(const MF &,
4543 OutVector &,
4544 const InVector &,
4545 const std::pair<unsigned int, unsigned int> &);
4546 };
4547 } // namespace MatrixFreeFunctions
4548
4549
4550
4551 // A implementation class for the worker object that runs the various
4552 // operations we want to perform during the matrix-free loop
4553 template <typename MF,
4554 typename InVector,
4555 typename OutVector,
4556 typename Container,
4557 bool is_constant>
4558 class MFWorker : public MFWorkerInterface
4559 {
4560 public:
4561 // An alias to make the arguments further down more readable
4562 using function_type = typename MatrixFreeFunctions::
4563 InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4564 function_type;
4565
4566 // constructor, binds all the arguments to this class
4567 MFWorker(const MF & matrix_free,
4568 const InVector & src,
4569 OutVector & dst,
4570 const bool zero_dst_vector_setting,
4571 const Container & container,
4572 function_type cell_function,
4573 function_type face_function,
4574 function_type boundary_function,
4575 const typename MF::DataAccessOnFaces src_vector_face_access =
4576 MF::DataAccessOnFaces::none,
4577 const typename MF::DataAccessOnFaces dst_vector_face_access =
4578 MF::DataAccessOnFaces::none,
4579 const std::function<void(const unsigned int, const unsigned int)>
4580 &operation_before_loop = {},
4581 const std::function<void(const unsigned int, const unsigned int)>
4582 & operation_after_loop = {},
4583 const unsigned int dof_handler_index_pre_post = 0)
4584 : matrix_free(matrix_free)
4585 , container(const_cast<Container &>(container))
4586 , cell_function(cell_function)
4587 , face_function(face_function)
4588 , boundary_function(boundary_function)
4589 , src(src)
4590 , dst(dst)
4591 , src_data_exchanger(matrix_free,
4592 src_vector_face_access,
4593 n_components(src))
4594 , dst_data_exchanger(matrix_free,
4595 dst_vector_face_access,
4596 n_components(dst))
4597 , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4598 , zero_dst_vector_setting(zero_dst_vector_setting &&
4599 !src_and_dst_are_same)
4600 , operation_before_loop(operation_before_loop)
4601 , operation_after_loop(operation_after_loop)
4602 , dof_handler_index_pre_post(dof_handler_index_pre_post)
4603 {}
4604
4605 // Runs the cell work. If no function is given, nothing is done
4606 virtual void
4607 cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4608 {
4609 if (cell_function != nullptr && cell_range.second > cell_range.first)
4610 for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4611 {
4612 const auto cell_subrange =
4613 matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4614
4615 if (cell_subrange.second <= cell_subrange.first)
4616 continue;
4617
4618 (container.*
4619 cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4620 }
4621 }
4622
4623 virtual void
4624 cell(const unsigned int range_index) override
4625 {
4626 process_range(cell_function,
4627 matrix_free.get_task_info().cell_partition_data_hp_ptr,
4628 matrix_free.get_task_info().cell_partition_data_hp,
4629 range_index);
4630 }
4631
4632 virtual void
4633 face(const unsigned int range_index) override
4634 {
4635 process_range(face_function,
4636 matrix_free.get_task_info().face_partition_data_hp_ptr,
4637 matrix_free.get_task_info().face_partition_data_hp,
4638 range_index);
4639 }
4640
4641 virtual void
4642 boundary(const unsigned int range_index) override
4643 {
4644 process_range(boundary_function,
4645 matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4646 matrix_free.get_task_info().boundary_partition_data_hp,
4647 range_index);
4648 }
4649
4650 private:
4651 void
4652 process_range(const function_type & fu,
4653 const std::vector<unsigned int> &ptr,
4654 const std::vector<unsigned int> &data,
4655 const unsigned int range_index)
4656 {
4657 if (fu == nullptr)
4658 return;
4659
4660 for (unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4661 (container.*fu)(matrix_free,
4662 this->dst,
4663 this->src,
4664 std::make_pair(data[2 * i], data[2 * i + 1]));
4665 }
4666
4667 public:
4668 // Starts the communication for the update ghost values operation. We
4669 // cannot call this update if ghost and destination are the same because
4670 // that would introduce spurious entries in the destination (there is also
4671 // the problem that reading from a vector that we also write to is usually
4672 // not intended in case there is overlap, but this is up to the
4673 // application code to decide and we cannot catch this case here).
4674 virtual void
4675 vector_update_ghosts_start() override
4676 {
4677 if (!src_and_dst_are_same)
4678 internal::update_ghost_values_start(src, src_data_exchanger);
4679 }
4680
4681 // Finishes the communication for the update ghost values operation
4682 virtual void
4683 vector_update_ghosts_finish() override
4684 {
4685 if (!src_and_dst_are_same)
4686 internal::update_ghost_values_finish(src, src_data_exchanger);
4687 }
4688
4689 // Starts the communication for the vector compress operation
4690 virtual void
4691 vector_compress_start() override
4692 {
4693 internal::compress_start(dst, dst_data_exchanger);
4694 }
4695
4696 // Finishes the communication for the vector compress operation
4697 virtual void
4698 vector_compress_finish() override
4699 {
4700 internal::compress_finish(dst, dst_data_exchanger);
4701 if (!src_and_dst_are_same)
4702 internal::reset_ghost_values(src, src_data_exchanger);
4703 }
4704
4705 // Zeros the given input vector
4706 virtual void
4707 zero_dst_vector_range(const unsigned int range_index) override
4708 {
4709 if (zero_dst_vector_setting)
4710 internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4711 }
4712
4713 virtual void
4714 cell_loop_pre_range(const unsigned int range_index) override
4715 {
4716 if (operation_before_loop)
4717 {
4718 const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4719 matrix_free.get_dof_info(dof_handler_index_pre_post);
4720 if (range_index == numbers::invalid_unsigned_int)
4721 {
4722 // Case with threaded loop -> currently no overlap implemented
4723 ::parallel::apply_to_subranges(
4724 0U,
4725 dof_info.vector_partitioner->locally_owned_size(),
4726 operation_before_loop,
4727 internal::VectorImplementation::minimum_parallel_grain_size);
4728 }
4729 else
4730 {
4731 AssertIndexRange(range_index,
4732 dof_info.cell_loop_pre_list_index.size() - 1);
4733 for (unsigned int id =
4734 dof_info.cell_loop_pre_list_index[range_index];
4735 id != dof_info.cell_loop_pre_list_index[range_index + 1];
4736 ++id)
4737 operation_before_loop(dof_info.cell_loop_pre_list[id].first,
4738 dof_info.cell_loop_pre_list[id].second);
4739 }
4740 }
4741 }
4742
4743 virtual void
4744 cell_loop_post_range(const unsigned int range_index) override
4745 {
4746 if (operation_after_loop)
4747 {
4748 const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4749 matrix_free.get_dof_info(dof_handler_index_pre_post);
4750 if (range_index == numbers::invalid_unsigned_int)
4751 {
4752 // Case with threaded loop -> currently no overlap implemented
4753 ::parallel::apply_to_subranges(
4754 0U,
4755 dof_info.vector_partitioner->locally_owned_size(),
4756 operation_after_loop,
4757 internal::VectorImplementation::minimum_parallel_grain_size);
4758 }
4759 else
4760 {
4761 AssertIndexRange(range_index,
4762 dof_info.cell_loop_post_list_index.size() - 1);
4763 for (unsigned int id =
4764 dof_info.cell_loop_post_list_index[range_index];
4765 id != dof_info.cell_loop_post_list_index[range_index + 1];
4766 ++id)
4767 operation_after_loop(dof_info.cell_loop_post_list[id].first,
4768 dof_info.cell_loop_post_list[id].second);
4769 }
4770 }
4771 }
4772
4773 private:
4774 const MF & matrix_free;
4775 Container & container;
4776 function_type cell_function;
4777 function_type face_function;
4778 function_type boundary_function;
4779
4780 const InVector &src;
4781 OutVector & dst;
4782 VectorDataExchange<MF::dimension,
4783 typename MF::value_type,
4784 typename MF::vectorized_value_type>
4785 src_data_exchanger;
4786 VectorDataExchange<MF::dimension,
4787 typename MF::value_type,
4788 typename MF::vectorized_value_type>
4789 dst_data_exchanger;
4790 const bool src_and_dst_are_same;
4791 const bool zero_dst_vector_setting;
4792 const std::function<void(const unsigned int, const unsigned int)>
4793 operation_before_loop;
4794 const std::function<void(const unsigned int, const unsigned int)>
4795 operation_after_loop;
4796 const unsigned int dof_handler_index_pre_post;
4797 };
4798
4799
4800
4805 template <class MF, typename InVector, typename OutVector>
4806 struct MFClassWrapper
4807 {
4808 using function_type =
4809 std::function<void(const MF &,
4810 OutVector &,
4811 const InVector &,
4812 const std::pair<unsigned int, unsigned int> &)>;
4813
4814 MFClassWrapper(const function_type cell,
4815 const function_type face,
4816 const function_type boundary)
4817 : cell(cell)
4818 , face(face)
4819 , boundary(boundary)
4820 {}
4821
4822 void
4823 cell_integrator(const MF & mf,
4824 OutVector & dst,
4825 const InVector & src,
4826 const std::pair<unsigned int, unsigned int> &range) const
4827 {
4828 if (cell)
4829 cell(mf, dst, src, range);
4830 }
4831
4832 void
4833 face_integrator(const MF & mf,
4834 OutVector & dst,
4835 const InVector & src,
4836 const std::pair<unsigned int, unsigned int> &range) const
4837 {
4838 if (face)
4839 face(mf, dst, src, range);
4840 }
4841
4842 void
4843 boundary_integrator(
4844 const MF & mf,
4845 OutVector & dst,
4846 const InVector & src,
4847 const std::pair<unsigned int, unsigned int> &range) const
4848 {
4849 if (boundary)
4850 boundary(mf, dst, src, range);
4851 }
4852
4853 const function_type cell;
4854 const function_type face;
4855 const function_type boundary;
4856 };
4857
4858} // end of namespace internal
4859
4860
4861
4862template <int dim, typename Number, typename VectorizedArrayType>
4863template <typename OutVector, typename InVector>
4864inline void
4865MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
4866 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4867 OutVector &,
4868 const InVector &,
4869 const std::pair<unsigned int, unsigned int> &)>
4870 & cell_operation,
4871 OutVector & dst,
4872 const InVector &src,
4873 const bool zero_dst_vector) const
4874{
4875 using Wrapper =
4876 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4877 InVector,
4878 OutVector>;
4879 Wrapper wrap(cell_operation, nullptr, nullptr);
4880 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4881 InVector,
4882 OutVector,
4883 Wrapper,
4884 true>
4885 worker(*this,
4886 src,
4887 dst,
4888 zero_dst_vector,
4889 wrap,
4890 &Wrapper::cell_integrator,
4891 &Wrapper::face_integrator,
4892 &Wrapper::boundary_integrator);
4893
4894 task_info.loop(worker);
4895}
4896
4897
4898
4899template <int dim, typename Number, typename VectorizedArrayType>
4900template <typename OutVector, typename InVector>
4901inline void
4902MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
4903 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4904 OutVector &,
4905 const InVector &,
4906 const std::pair<unsigned int, unsigned int> &)>
4907 & cell_operation,
4908 OutVector & dst,
4909 const InVector &src,
4910 const std::function<void(const unsigned int, const unsigned int)>
4911 &operation_before_loop,
4912 const std::function<void(const unsigned int, const unsigned int)>
4913 & operation_after_loop,
4914 const unsigned int dof_handler_index_pre_post) const
4915{
4916 using Wrapper =
4917 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4918 InVector,
4919 OutVector>;
4920 Wrapper wrap(cell_operation, nullptr, nullptr);
4921 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4922 InVector,
4923 OutVector,
4924 Wrapper,
4925 true>
4926 worker(*this,
4927 src,
4928 dst,
4929 false,
4930 wrap,
4931 &Wrapper::cell_integrator,
4932 &Wrapper::face_integrator,
4933 &Wrapper::boundary_integrator,
4934 DataAccessOnFaces::none,
4935 DataAccessOnFaces::none,
4936 operation_before_loop,
4937 operation_after_loop,
4938 dof_handler_index_pre_post);
4939
4940 task_info.loop(worker);
4941}
4942
4943
4944
4945template <int dim, typename Number, typename VectorizedArrayType>
4946template <typename OutVector, typename InVector>
4947inline void
4948MatrixFree<dim, Number, VectorizedArrayType>::loop(
4949 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4950 OutVector &,
4951 const InVector &,
4952 const std::pair<unsigned int, unsigned int> &)>
4953 &cell_operation,
4954 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4955 OutVector &,
4956 const InVector &,
4957 const std::pair<unsigned int, unsigned int> &)>
4958 &face_operation,
4959 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4960 OutVector &,
4961 const InVector &,
4962 const std::pair<unsigned int, unsigned int> &)>
4963 & boundary_operation,
4964 OutVector & dst,
4965 const InVector & src,
4966 const bool zero_dst_vector,
4967 const DataAccessOnFaces dst_vector_face_access,
4968 const DataAccessOnFaces src_vector_face_access) const
4969{
4970 using Wrapper =
4971 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4972 InVector,
4973 OutVector>;
4974 Wrapper wrap(cell_operation, face_operation, boundary_operation);
4975 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4976 InVector,
4977 OutVector,
4978 Wrapper,
4979 true>
4980 worker(*this,
4981 src,
4982 dst,
4983 zero_dst_vector,
4984 wrap,
4985 &Wrapper::cell_integrator,
4986 &Wrapper::face_integrator,
4987 &Wrapper::boundary_integrator,
4988 src_vector_face_access,
4989 dst_vector_face_access);
4990
4991 task_info.loop(worker);
4992}
4993
4994
4995
4996template <int dim, typename Number, typename VectorizedArrayType>
4997template <typename CLASS, typename OutVector, typename InVector>
4998inline void
4999MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5000 void (CLASS::*function_pointer)(
5001 const MatrixFree<dim, Number, VectorizedArrayType> &,
5002 OutVector &,
5003 const InVector &,
5004 const std::pair<unsigned int, unsigned int> &) const,
5005 const CLASS * owning_class,
5006 OutVector & dst,
5007 const InVector &src,
5008 const bool zero_dst_vector) const
5009{
5010 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5011 InVector,
5012 OutVector,
5013 CLASS,
5014 true>
5015 worker(*this,
5016 src,
5017 dst,
5018 zero_dst_vector,
5019 *owning_class,
5020 function_pointer,
5021 nullptr,
5022 nullptr);
5023 task_info.loop(worker);
5024}
5025
5026
5027
5028template <int dim, typename Number, typename VectorizedArrayType>
5029template <typename CLASS, typename OutVector, typename InVector>
5030inline void
5031MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5032 void (CLASS::*function_pointer)(
5033 const MatrixFree<dim, Number, VectorizedArrayType> &,
5034 OutVector &,
5035 const InVector &,
5036 const std::pair<unsigned int, unsigned int> &) const,
5037 const CLASS * owning_class,
5038 OutVector & dst,
5039 const InVector &src,
5040 const std::function<void(const unsigned int, const unsigned int)>
5041 &operation_before_loop,
5042 const std::function<void(const unsigned int, const unsigned int)>
5043 & operation_after_loop,
5044 const unsigned int dof_handler_index_pre_post) const
5045{
5046 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5047 InVector,
5048 OutVector,
5049 CLASS,
5050 true>
5051 worker(*this,
5052 src,
5053 dst,
5054 false,
5055 *owning_class,
5056 function_pointer,
5057 nullptr,
5058 nullptr,
5059 DataAccessOnFaces::none,
5060 DataAccessOnFaces::none,
5061 operation_before_loop,
5062 operation_after_loop,
5063 dof_handler_index_pre_post);
5064 task_info.loop(worker);
5065}
5066
5067
5068
5069template <int dim, typename Number, typename VectorizedArrayType>
5070template <typename CLASS, typename OutVector, typename InVector>
5071inline void
5072MatrixFree<dim, Number, VectorizedArrayType>::loop(
5073 void (CLASS::*cell_operation)(
5074 const MatrixFree<dim, Number, VectorizedArrayType> &,
5075 OutVector &,
5076 const InVector &,
5077 const std::pair<unsigned int, unsigned int> &) const,
5078 void (CLASS::*face_operation)(
5079 const MatrixFree<dim, Number, VectorizedArrayType> &,
5080 OutVector &,
5081 const InVector &,
5082 const std::pair<unsigned int, unsigned int> &) const,
5083 void (CLASS::*boundary_operation)(
5084 const MatrixFree<dim, Number, VectorizedArrayType> &,
5085 OutVector &,
5086 const InVector &,
5087 const std::pair<unsigned int, unsigned int> &) const,
5088 const CLASS * owning_class,
5089 OutVector & dst,
5090 const InVector & src,
5091 const bool zero_dst_vector,
5092 const DataAccessOnFaces dst_vector_face_access,
5093 const DataAccessOnFaces src_vector_face_access) const
5094{
5095 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5096 InVector,
5097 OutVector,
5098 CLASS,
5099 true>
5100 worker(*this,
5101 src,
5102 dst,
5103 zero_dst_vector,
5104 *owning_class,
5105 cell_operation,
5106 face_operation,
5107 boundary_operation,
5108 src_vector_face_access,
5109 dst_vector_face_access);
5110 task_info.loop(worker);
5111}
5112
5113
5114
5115template <int dim, typename Number, typename VectorizedArrayType>
5116template <typename CLASS, typename OutVector, typename InVector>
5117inline void
5118MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5119 void (CLASS::*function_pointer)(
5120 const MatrixFree<dim, Number, VectorizedArrayType> &,
5121 OutVector &,
5122 const InVector &,
5123 const std::pair<unsigned int, unsigned int> &),
5124 CLASS * owning_class,
5125 OutVector & dst,
5126 const InVector &src,
5127 const bool zero_dst_vector) const
5128{
5129 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5130 InVector,
5131 OutVector,
5132 CLASS,
5133 false>
5134 worker(*this,
5135 src,
5136 dst,
5137 zero_dst_vector,
5138 *owning_class,
5139 function_pointer,
5140 nullptr,
5141 nullptr);
5142 task_info.loop(worker);
5143}
5144
5145
5146
5147template <int dim, typename Number, typename VectorizedArrayType>
5148template <typename CLASS, typename OutVector, typename InVector>
5149inline void
5150MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5151 void (CLASS::*function_pointer)(
5152 const MatrixFree<dim, Number, VectorizedArrayType> &,
5153 OutVector &,
5154 const InVector &,
5155 const std::pair<unsigned int, unsigned int> &),
5156 CLASS * owning_class,
5157 OutVector & dst,
5158 const InVector &src,
5159 const std::function<void(const unsigned int, const unsigned int)>
5160 &operation_before_loop,
5161 const std::function<void(const unsigned int, const unsigned int)>
5162 & operation_after_loop,
5163 const unsigned int dof_handler_index_pre_post) const
5164{
5165 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5166 InVector,
5167 OutVector,
5168 CLASS,
5169 false>
5170 worker(*this,
5171 src,
5172 dst,
5173 false,
5174 *owning_class,
5175 function_pointer,
5176 nullptr,
5177 nullptr,
5178 DataAccessOnFaces::none,
5179 DataAccessOnFaces::none,
5180 operation_before_loop,
5181 operation_after_loop,
5182 dof_handler_index_pre_post);
5183 task_info.loop(worker);
5184}
5185
5186
5187
5188template <int dim, typename Number, typename VectorizedArrayType>
5189template <typename CLASS, typename OutVector, typename InVector>
5190inline void
5191MatrixFree<dim, Number, VectorizedArrayType>::loop(
5192 void (CLASS::*cell_operation)(
5193 const MatrixFree<dim, Number, VectorizedArrayType> &,
5194 OutVector &,
5195 const InVector &,
5196 const std::pair<unsigned int, unsigned int> &),
5197 void (CLASS::*face_operation)(
5198 const MatrixFree<dim, Number, VectorizedArrayType> &,
5199 OutVector &,
5200 const InVector &,
5201 const std::pair<unsigned int, unsigned int> &),
5202 void (CLASS::*boundary_operation)(
5203 const MatrixFree<dim, Number, VectorizedArrayType> &,
5204 OutVector &,
5205 const InVector &,
5206 const std::pair<unsigned int, unsigned int> &),
5207 CLASS * owning_class,
5208 OutVector & dst,
5209 const InVector & src,
5210 const bool zero_dst_vector,
5211 const DataAccessOnFaces dst_vector_face_access,
5212 const DataAccessOnFaces src_vector_face_access) const
5213{
5214 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5215 InVector,
5216 OutVector,
5217 CLASS,
5218 false>
5219 worker(*this,
5220 src,
5221 dst,
5222 zero_dst_vector,
5223 *owning_class,
5224 cell_operation,
5225 face_operation,
5226 boundary_operation,
5227 src_vector_face_access,
5228 dst_vector_face_access);
5229 task_info.loop(worker);
5230}
5231
5232
5233
5234template <int dim, typename Number, typename VectorizedArrayType>
5235template <typename CLASS, typename OutVector, typename InVector>
5236inline void
5237MatrixFree<dim, Number, VectorizedArrayType>::loop_cell_centric(
5238 void (CLASS::*function_pointer)(
5239 const MatrixFree<dim, Number, VectorizedArrayType> &,
5240 OutVector &,
5241 const InVector &,
5242 const std::pair<unsigned int, unsigned int> &) const,
5243 const CLASS * owning_class,
5244 OutVector & dst,
5245 const InVector & src,
5246 const bool zero_dst_vector,
5247 const DataAccessOnFaces src_vector_face_access) const
5248{
5249 auto src_vector_face_access_temp = src_vector_face_access;
5250 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5251 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5252 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5253 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5254
5255 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5256 InVector,
5257 OutVector,
5258 CLASS,
5259 true>
5260 worker(*this,
5261 src,
5262 dst,
5263 zero_dst_vector,
5264 *owning_class,
5265 function_pointer,
5266 nullptr,
5267 nullptr,
5268 src_vector_face_access_temp,
5269 DataAccessOnFaces::none);
5270 task_info.loop(worker);
5271}
5272
5273
5274
5275template <int dim, typename Number, typename VectorizedArrayType>
5276template <typename CLASS, typename OutVector, typename InVector>
5277inline void
5278MatrixFree<dim, Number, VectorizedArrayType>::loop_cell_centric(
5279 void (CLASS::*function_pointer)(
5280 const MatrixFree<dim, Number, VectorizedArrayType> &,
5281 OutVector &,
5282 const InVector &,
5283 const std::pair<unsigned int, unsigned int> &),
5284 CLASS * owning_class,
5285 OutVector & dst,
5286 const InVector & src,
5287 const bool zero_dst_vector,
5288 const DataAccessOnFaces src_vector_face_access) const
5289{
5290 auto src_vector_face_access_temp = src_vector_face_access;
5291 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5292 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5293 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5294 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5295
5296 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5297 InVector,
5298 OutVector,
5299 CLASS,
5300 false>
5301 worker(*this,
5302 src,
5303 dst,
5304 zero_dst_vector,
5305 *owning_class,
5306 function_pointer,
5307 nullptr,
5308 nullptr,
5309 src_vector_face_access_temp,
5310 DataAccessOnFaces::none);
5311 task_info.loop(worker);
5312}
5313
5314
5315
5316template <int dim, typename Number, typename VectorizedArrayType>
5317template <typename OutVector, typename InVector>
5318inline void
5319MatrixFree<dim, Number, VectorizedArrayType>::loop_cell_centric(
5320 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5321 OutVector &,
5322 const InVector &,
5323 const std::pair<unsigned int, unsigned int> &)>
5324 & cell_operation,
5325 OutVector & dst,
5326 const InVector & src,
5327 const bool zero_dst_vector,
5328 const DataAccessOnFaces src_vector_face_access) const
5329{
5330 auto src_vector_face_access_temp = src_vector_face_access;
5331 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5332 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5333 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5334 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5335
5336 using Wrapper =
5337 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5338 InVector,
5339 OutVector>;
5340 Wrapper wrap(cell_operation, nullptr, nullptr);
5341
5342 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5343 InVector,
5344 OutVector,
5345 Wrapper,
5346 true>
5347 worker(*this,
5348 src,
5349 dst,
5350 zero_dst_vector,
5351 wrap,
5352 &Wrapper::cell_integrator,
5353 &Wrapper::face_integrator,
5354 &Wrapper::boundary_integrator,
5355 src_vector_face_access_temp,
5356 DataAccessOnFaces::none);
5357 task_info.loop(worker);
5358}
5359
5360
5361#endif // ifndef DOXYGEN
5362
5363
5364
5365DEAL_II_NAMESPACE_CLOSE
5366
5367#endif
affine_constraints.h
aligned_vector.h
block_vector_base.h
MPI_Comm
Mapping
Abstract base class for mapping classes.
Definition: mapping.h:304
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:2114
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::shape_info
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > > shape_info
Definition: matrix_free.h:2120
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::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:2144
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::get_face_category
std::pair< unsigned int, unsigned int > get_face_category(const unsigned int macro_face) const
MatrixFree::mapping_is_initialized
bool mapping_is_initialized
Definition: matrix_free.h:2161
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:2184
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:2128
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:2101
MatrixFree::vectorized_value_type
VectorizedArrayType vectorized_value_type
Definition: matrix_free.h:128
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_boundary_id
types::boundary_id get_boundary_id(const unsigned int macro_face) 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:2172
MatrixFree::value_type
Number value_type
Definition: matrix_free.h:127
MatrixFree::memory_consumption
std::size_t memory_consumption() 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:2179
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:2087
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:2093
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:2156
MatrixFree::face_info
internal::MatrixFreeFunctions::FaceInfo< VectorizedArrayType::size()> face_info
Definition: matrix_free.h:2151
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::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:2107
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:2137
MatrixFree::n_base_elements
unsigned int n_base_elements(const unsigned int dof_handler_index) const
MatrixFree::dimension
static const unsigned int dimension
Definition: matrix_free.h:133
MatrixFree::is_supported
static bool is_supported(const FiniteElement< dim, spacedim > &fe)
Threads::ThreadLocalStorage
A class that provides a separate storage location on each thread that accesses the object.
Definition: thread_local_storage.h:102
DEAL_II_DEPRECATED
#define DEAL_II_DEPRECATED
Definition: config.h:162
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:416
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:454
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
level
unsigned int level
Definition: grid_out.cc:4590
grid_tools.h
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
StandardExceptions::ExcIndexRange
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1465
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
DataOutBase::none
@ none
Definition: data_out_base.h:1555
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcNotInitialized
static ::ExceptionBase & ExcNotInitialized()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
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::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_info.h
mapping_q1.h
dof_info.h
LAPACKSupport::U
static const char U
Definition: lapack_support.h:169
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:1097
internal::VectorImplementation::minimum_parallel_grain_size
unsigned int minimum_parallel_grain_size
Definition: parallel.cc:34
numbers::invalid_boundary_id
const types::boundary_id invalid_boundary_id
Definition: types.h:239
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:196
parallel::internal::apply_to_subranges
void apply_to_subranges(const tbb::blocked_range< RangeType > &range, const Function &f)
Definition: parallel.h:367
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5484
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:344
MatrixFree::AdditionalData::mapping_update_flags_inner_faces
UpdateFlags mapping_update_flags_inner_faces
Definition: matrix_free.h:410
MatrixFree::AdditionalData::AdditionalData
AdditionalData(const TasksParallelScheme tasks_parallel_scheme=partition_partition, const unsigned int tasks_block_size=0, const UpdateFlags mapping_update_flags=update_gradients|update_JxW_values, const UpdateFlags mapping_update_flags_boundary_faces=update_default, const UpdateFlags mapping_update_flags_inner_faces=update_default, const UpdateFlags mapping_update_flags_faces_by_cells=update_default, const unsigned int mg_level=numbers::invalid_unsigned_int, const bool store_plain_indices=true, const bool initialize_indices=true, const bool initialize_mapping=true, const bool overlap_communication_computation=true, const bool hold_all_faces_to_owned_cells=false, const bool cell_vectorization_categories_strict=false)
Definition: matrix_free.h:217
MatrixFree::AdditionalData::cell_vectorization_category
std::vector< unsigned int > cell_vectorization_category
Definition: matrix_free.h:524
MatrixFree::AdditionalData::operator=
AdditionalData & operator=(const AdditionalData &other)
Definition: matrix_free.h:283
MatrixFree::AdditionalData::level_mg_handler
unsigned int & level_mg_handler
Definition: matrix_free.h:455
MatrixFree::AdditionalData::mapping_update_flags_boundary_faces
UpdateFlags mapping_update_flags_boundary_faces
Definition: matrix_free.h:389
MatrixFree::AdditionalData::mapping_update_flags
UpdateFlags mapping_update_flags
Definition: matrix_free.h:368
MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells
UpdateFlags mapping_update_flags_faces_by_cells
Definition: matrix_free.h:438
MatrixFree::AdditionalData::AdditionalData
AdditionalData(const AdditionalData &other)
Definition: matrix_free.h:253
MatrixFree::AdditionalData::tasks_block_size
unsigned int tasks_block_size
Definition: matrix_free.h:355
PointerComparison::equal
static bool equal(const T *p1, const T *p2)
Definition: template_constraints.h:228
internal::MatrixFreeFunctions::DoFInfo::cell_loop_pre_list_index
std::vector< unsigned int > cell_loop_pre_list_index
Definition: dof_info.h:721
internal::MatrixFreeFunctions::DoFInfo::cell_loop_post_list_index
std::vector< unsigned int > cell_loop_post_list_index
Definition: dof_info.h:734
internal::MatrixFreeFunctions::DoFInfo::vector_zero_range_list
std::vector< std::pair< unsigned int, unsigned int > > vector_zero_range_list
Definition: dof_info.h:714
internal::MatrixFreeFunctions::DoFInfo::vector_partitioner
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition: dof_info.h:565
internal::MatrixFreeFunctions::DoFInfo::cell_loop_pre_list
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_pre_list
Definition: dof_info.h:727
internal::MatrixFreeFunctions::DoFInfo::vector_zero_range_list_index
std::vector< unsigned int > vector_zero_range_list_index
Definition: dof_info.h:709
internal::MatrixFreeFunctions::DoFInfo::cell_loop_post_list
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_post_list
Definition: dof_info.h:740
internal::MatrixFreeFunctions::FaceInfo::cell_and_face_to_plain_faces
::Table< 3, unsigned int > cell_and_face_to_plain_faces
Definition: face_info.h:172
internal::MatrixFreeFunctions::FaceInfo::faces
std::vector< FaceToCellTopology< vectorization_width > > faces
Definition: face_info.h:165
internal::MatrixFreeFunctions::FaceInfo::cell_and_face_boundary_id
::Table< 3, types::boundary_id > cell_and_face_boundary_id
Definition: face_info.h:178
internal::MatrixFreeFunctions::TaskInfo::boundary_partition_data
std::vector< unsigned int > boundary_partition_data
Definition: task_info.h:518
internal::MatrixFreeFunctions::TaskInfo::loop
void loop(MFWorkerInterface &worker) const
Definition: task_info.cc:348
internal::MatrixFreeFunctions::TaskInfo::cell_partition_data
std::vector< unsigned int > cell_partition_data
Definition: task_info.h:474
internal::MatrixFreeFunctions::TaskInfo::face_partition_data
std::vector< unsigned int > face_partition_data
Definition: task_info.h:496
task_info.h
template_constraints.h
thread_local_storage.h
type_traits.h
vector_operation.h
Generated by doxygen 1.9.4