Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
matrix_free.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2011 - 2023 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
34#include <deal.II/hp/mapping_collection.h>
35#include <deal.II/hp/q_collection.h>
36
37#include <deal.II/lac/affine_constraints.h>
38#include <deal.II/lac/block_vector_base.h>
39#include <deal.II/lac/la_parallel_vector.h>
40#include <deal.II/lac/vector_operation.h>
41
42#include <deal.II/matrix_free/dof_info.h>
43#include <deal.II/matrix_free/mapping_info.h>
44#include <deal.II/matrix_free/shape_info.h>
45#include <deal.II/matrix_free/task_info.h>
46#include <deal.II/matrix_free/type_traits.h>
47
48#include <cstdlib>
49#include <limits>
50#include <list>
51#include <memory>
52
53
54DEAL_II_NAMESPACE_OPEN
55
56
57
109template <int dim,
110 typename Number = double,
111 typename VectorizedArrayType = VectorizedArray<Number>>
112class MatrixFree : public Subscriptor
113{
114 static_assert(
115 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
116 "Type of Number and of VectorizedArrayType do not match.");
117
118public:
123 using value_type = Number;
124 using vectorized_value_type = VectorizedArrayType;
125
129 static constexpr unsigned int dimension = dim;
130
183 struct AdditionalData
184 {
188 using MatrixFreeType = MatrixFree<dim, Number, VectorizedArrayType>;
189
195 enum TasksParallelScheme
196 {
200 none = internal::MatrixFreeFunctions::TaskInfo::none,
204 partition_partition =
205 internal::MatrixFreeFunctions::TaskInfo::partition_partition,
209 partition_color =
210 internal::MatrixFreeFunctions::TaskInfo::partition_color,
215 color = internal::MatrixFreeFunctions::TaskInfo::color
216 };
217
221 AdditionalData(
222 const TasksParallelScheme tasks_parallel_scheme = partition_partition,
223 const unsigned int tasks_block_size = 0,
224 const UpdateFlags mapping_update_flags = update_gradients |
225 update_JxW_values,
226 const UpdateFlags mapping_update_flags_boundary_faces = update_default,
227 const UpdateFlags mapping_update_flags_inner_faces = update_default,
228 const UpdateFlags mapping_update_flags_faces_by_cells = update_default,
229 const unsigned int mg_level = numbers::invalid_unsigned_int,
230 const bool store_plain_indices = true,
231 const bool initialize_indices = true,
232 const bool initialize_mapping = true,
233 const bool overlap_communication_computation = true,
234 const bool hold_all_faces_to_owned_cells = false,
235 const bool cell_vectorization_categories_strict = false,
236 const bool allow_ghosted_vectors_in_loops = true)
237 : tasks_parallel_scheme(tasks_parallel_scheme)
238 , tasks_block_size(tasks_block_size)
239 , mapping_update_flags(mapping_update_flags)
240 , mapping_update_flags_boundary_faces(mapping_update_flags_boundary_faces)
241 , mapping_update_flags_inner_faces(mapping_update_flags_inner_faces)
242 , mapping_update_flags_faces_by_cells(mapping_update_flags_faces_by_cells)
243 , mg_level(mg_level)
244 , store_plain_indices(store_plain_indices)
245 , initialize_indices(initialize_indices)
246 , initialize_mapping(initialize_mapping)
247 , overlap_communication_computation(overlap_communication_computation)
248 , hold_all_faces_to_owned_cells(hold_all_faces_to_owned_cells)
249 , cell_vectorization_categories_strict(
250 cell_vectorization_categories_strict)
251 , allow_ghosted_vectors_in_loops(allow_ghosted_vectors_in_loops)
252 , communicator_sm(MPI_COMM_SELF)
253 {}
254
258 AdditionalData(const AdditionalData &other)
259 : tasks_parallel_scheme(other.tasks_parallel_scheme)
260 , tasks_block_size(other.tasks_block_size)
261 , mapping_update_flags(other.mapping_update_flags)
262 , mapping_update_flags_boundary_faces(
263 other.mapping_update_flags_boundary_faces)
264 , mapping_update_flags_inner_faces(other.mapping_update_flags_inner_faces)
265 , mapping_update_flags_faces_by_cells(
266 other.mapping_update_flags_faces_by_cells)
267 , mg_level(other.mg_level)
268 , store_plain_indices(other.store_plain_indices)
269 , initialize_indices(other.initialize_indices)
270 , initialize_mapping(other.initialize_mapping)
271 , overlap_communication_computation(
272 other.overlap_communication_computation)
273 , hold_all_faces_to_owned_cells(other.hold_all_faces_to_owned_cells)
274 , cell_vectorization_category(other.cell_vectorization_category)
275 , cell_vectorization_categories_strict(
276 other.cell_vectorization_categories_strict)
277 , allow_ghosted_vectors_in_loops(other.allow_ghosted_vectors_in_loops)
278 , communicator_sm(other.communicator_sm)
279 {}
280
284 AdditionalData &
285 operator=(const AdditionalData &other)
286 {
287 tasks_parallel_scheme = other.tasks_parallel_scheme;
288 tasks_block_size = other.tasks_block_size;
289 mapping_update_flags = other.mapping_update_flags;
290 mapping_update_flags_boundary_faces =
291 other.mapping_update_flags_boundary_faces;
292 mapping_update_flags_inner_faces = other.mapping_update_flags_inner_faces;
293 mapping_update_flags_faces_by_cells =
294 other.mapping_update_flags_faces_by_cells;
295 mg_level = other.mg_level;
296 store_plain_indices = other.store_plain_indices;
297 initialize_indices = other.initialize_indices;
298 initialize_mapping = other.initialize_mapping;
299 overlap_communication_computation =
300 other.overlap_communication_computation;
301 hold_all_faces_to_owned_cells = other.hold_all_faces_to_owned_cells;
302 cell_vectorization_category = other.cell_vectorization_category;
303 cell_vectorization_categories_strict =
304 other.cell_vectorization_categories_strict;
305 allow_ghosted_vectors_in_loops = other.allow_ghosted_vectors_in_loops;
306 communicator_sm = other.communicator_sm;
307
308 return *this;
309 }
310
347 TasksParallelScheme tasks_parallel_scheme;
348
358 unsigned int tasks_block_size;
359
373 UpdateFlags mapping_update_flags;
374
393 UpdateFlags mapping_update_flags_boundary_faces;
394
413 UpdateFlags mapping_update_flags_inner_faces;
414
441 UpdateFlags mapping_update_flags_faces_by_cells;
442
451 unsigned int mg_level;
452
459 bool store_plain_indices;
460
470 bool initialize_indices;
471
479 bool initialize_mapping;
480
493 bool overlap_communication_computation;
494
502 bool hold_all_faces_to_owned_cells;
503
525 std::vector<unsigned int> cell_vectorization_category;
526
536 bool cell_vectorization_categories_strict;
537
548 bool allow_ghosted_vectors_in_loops;
549
553 MPI_Comm communicator_sm;
554 };
555
563 MatrixFree();
564
568 MatrixFree(const MatrixFree<dim, Number, VectorizedArrayType> &other);
569
573 ~MatrixFree() override = default;
574
587 template <typename QuadratureType, typename number2, typename MappingType>
588 void
589 reinit(const MappingType & mapping,
590 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
631 template <typename QuadratureType, typename number2, typename MappingType>
632 void
633 reinit(const MappingType & mapping,
634 const std::vector<const DoFHandler<dim> *> & dof_handler,
635 const std::vector<const AffineConstraints<number2> *> &constraint,
636 const QuadratureType & quad,
637 const AdditionalData &additional_data = AdditionalData());
638
644 void
645 copy_from(
646 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free_base);
647
657 void
658 update_mapping(const Mapping<dim> &mapping);
659
663 void
664 update_mapping(const std::shared_ptr<hp::MappingCollection<dim>> &mapping);
665
670 void
671 clear();
672
686 enum class DataAccessOnFaces
687 {
694 none,
695
706 values,
707
716 values_all_faces,
717
728 gradients,
729
738 gradients_all_faces,
739
746 unspecified
747 };
748
793 template <typename OutVector, typename InVector>
794 void
795 cell_loop(const std::function<void(
796 const MatrixFree<dim, Number, VectorizedArrayType> &,
797 OutVector &,
798 const InVector &,
799 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
800 OutVector & dst,
801 const InVector & src,
802 const bool zero_dst_vector = false) const;
803
850 template <typename CLASS, typename OutVector, typename InVector>
851 void
852 cell_loop(void (CLASS::*cell_operation)(
853 const MatrixFree &,
854 OutVector &,
855 const InVector &,
856 const std::pair<unsigned int, unsigned int> &) const,
857 const CLASS * owning_class,
858 OutVector & dst,
859 const InVector &src,
860 const bool zero_dst_vector = false) const;
861
865 template <typename CLASS, typename OutVector, typename InVector>
866 void
867 cell_loop(void (CLASS::*cell_operation)(
868 const MatrixFree &,
869 OutVector &,
870 const InVector &,
871 const std::pair<unsigned int, unsigned int> &),
872 CLASS * owning_class,
873 OutVector & dst,
874 const InVector &src,
875 const bool zero_dst_vector = false) const;
876
961 template <typename CLASS, typename OutVector, typename InVector>
962 void
963 cell_loop(void (CLASS::*cell_operation)(
964 const MatrixFree &,
965 OutVector &,
966 const InVector &,
967 const std::pair<unsigned int, unsigned int> &) const,
968 const CLASS * owning_class,
969 OutVector & dst,
970 const InVector &src,
971 const std::function<void(const unsigned int, const unsigned int)>
972 &operation_before_loop,
973 const std::function<void(const unsigned int, const unsigned int)>
974 & operation_after_loop,
975 const unsigned int dof_handler_index_pre_post = 0) const;
976
980 template <typename CLASS, typename OutVector, typename InVector>
981 void
982 cell_loop(void (CLASS::*cell_operation)(
983 const MatrixFree &,
984 OutVector &,
985 const InVector &,
986 const std::pair<unsigned int, unsigned int> &),
987 CLASS * owning_class,
988 OutVector & dst,
989 const InVector &src,
990 const std::function<void(const unsigned int, const unsigned int)>
991 &operation_before_loop,
992 const std::function<void(const unsigned int, const unsigned int)>
993 & operation_after_loop,
994 const unsigned int dof_handler_index_pre_post = 0) const;
995
1000 template <typename OutVector, typename InVector>
1001 void
1002 cell_loop(const std::function<void(
1003 const MatrixFree<dim, Number, VectorizedArrayType> &,
1004 OutVector &,
1005 const InVector &,
1006 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1007 OutVector & dst,
1008 const InVector & src,
1009 const std::function<void(const unsigned int, const unsigned int)>
1010 &operation_before_loop,
1011 const std::function<void(const unsigned int, const unsigned int)>
1012 & operation_after_loop,
1013 const unsigned int dof_handler_index_pre_post = 0) const;
1014
1090 template <typename OutVector, typename InVector>
1091 void
1092 loop(const std::function<
1093 void(const MatrixFree<dim, Number, VectorizedArrayType> &,
1094 OutVector &,
1095 const InVector &,
1096 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1097 const std::function<
1098 void(const MatrixFree<dim, Number, VectorizedArrayType> &,
1099 OutVector &,
1100 const InVector &,
1101 const std::pair<unsigned int, unsigned int> &)> &face_operation,
1102 const std::function<void(
1103 const MatrixFree<dim, Number, VectorizedArrayType> &,
1104 OutVector &,
1105 const InVector &,
1106 const std::pair<unsigned int, unsigned int> &)> &boundary_operation,
1107 OutVector & dst,
1108 const InVector & src,
1109 const bool zero_dst_vector = false,
1110 const DataAccessOnFaces dst_vector_face_access =
1111 DataAccessOnFaces::unspecified,
1112 const DataAccessOnFaces src_vector_face_access =
1113 DataAccessOnFaces::unspecified) const;
1114
1200 template <typename CLASS, typename OutVector, typename InVector>
1201 void
1202 loop(
1203 void (CLASS::*cell_operation)(const MatrixFree &,
1204 OutVector &,
1205 const InVector &,
1206 const std::pair<unsigned int, unsigned int> &)
1207 const,
1208 void (CLASS::*face_operation)(const MatrixFree &,
1209 OutVector &,
1210 const InVector &,
1211 const std::pair<unsigned int, unsigned int> &)
1212 const,
1213 void (CLASS::*boundary_operation)(
1214 const MatrixFree &,
1215 OutVector &,
1216 const InVector &,
1217 const std::pair<unsigned int, unsigned int> &) const,
1218 const CLASS * owning_class,
1219 OutVector & dst,
1220 const InVector & src,
1221 const bool zero_dst_vector = false,
1222 const DataAccessOnFaces dst_vector_face_access =
1223 DataAccessOnFaces::unspecified,
1224 const DataAccessOnFaces src_vector_face_access =
1225 DataAccessOnFaces::unspecified) const;
1226
1230 template <typename CLASS, typename OutVector, typename InVector>
1231 void
1232 loop(void (CLASS::*cell_operation)(
1233 const MatrixFree &,
1234 OutVector &,
1235 const InVector &,
1236 const std::pair<unsigned int, unsigned int> &),
1237 void (CLASS::*face_operation)(
1238 const MatrixFree &,
1239 OutVector &,
1240 const InVector &,
1241 const std::pair<unsigned int, unsigned int> &),
1242 void (CLASS::*boundary_operation)(
1243 const MatrixFree &,
1244 OutVector &,
1245 const InVector &,
1246 const std::pair<unsigned int, unsigned int> &),
1247 CLASS * owning_class,
1248 OutVector & dst,
1249 const InVector & src,
1250 const bool zero_dst_vector = false,
1251 const DataAccessOnFaces dst_vector_face_access =
1252 DataAccessOnFaces::unspecified,
1253 const DataAccessOnFaces src_vector_face_access =
1254 DataAccessOnFaces::unspecified) const;
1255
1375 template <typename CLASS, typename OutVector, typename InVector>
1376 void
1377 loop(
1378 void (CLASS::*cell_operation)(const MatrixFree &,
1379 OutVector &,
1380 const InVector &,
1381 const std::pair<unsigned int, unsigned int> &)
1382 const,
1383 void (CLASS::*face_operation)(const MatrixFree &,
1384 OutVector &,
1385 const InVector &,
1386 const std::pair<unsigned int, unsigned int> &)
1387 const,
1388 void (CLASS::*boundary_operation)(
1389 const MatrixFree &,
1390 OutVector &,
1391 const InVector &,
1392 const std::pair<unsigned int, unsigned int> &) const,
1393 const CLASS * owning_class,
1394 OutVector & dst,
1395 const InVector &src,
1396 const std::function<void(const unsigned int, const unsigned int)>
1397 &operation_before_loop,
1398 const std::function<void(const unsigned int, const unsigned int)>
1399 & operation_after_loop,
1400 const unsigned int dof_handler_index_pre_post = 0,
1401 const DataAccessOnFaces dst_vector_face_access =
1402 DataAccessOnFaces::unspecified,
1403 const DataAccessOnFaces src_vector_face_access =
1404 DataAccessOnFaces::unspecified) const;
1405
1409 template <typename CLASS, typename OutVector, typename InVector>
1410 void
1411 loop(void (CLASS::*cell_operation)(
1412 const MatrixFree &,
1413 OutVector &,
1414 const InVector &,
1415 const std::pair<unsigned int, unsigned int> &),
1416 void (CLASS::*face_operation)(
1417 const MatrixFree &,
1418 OutVector &,
1419 const InVector &,
1420 const std::pair<unsigned int, unsigned int> &),
1421 void (CLASS::*boundary_operation)(
1422 const MatrixFree &,
1423 OutVector &,
1424 const InVector &,
1425 const std::pair<unsigned int, unsigned int> &),
1426 const CLASS * owning_class,
1427 OutVector & dst,
1428 const InVector &src,
1429 const std::function<void(const unsigned int, const unsigned int)>
1430 &operation_before_loop,
1431 const std::function<void(const unsigned int, const unsigned int)>
1432 & operation_after_loop,
1433 const unsigned int dof_handler_index_pre_post = 0,
1434 const DataAccessOnFaces dst_vector_face_access =
1435 DataAccessOnFaces::unspecified,
1436 const DataAccessOnFaces src_vector_face_access =
1437 DataAccessOnFaces::unspecified) const;
1438
1444 template <typename OutVector, typename InVector>
1445 void
1446 loop(const std::function<
1447 void(const MatrixFree<dim, Number, VectorizedArrayType> &,
1448 OutVector &,
1449 const InVector &,
1450 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1451 const std::function<
1452 void(const MatrixFree<dim, Number, VectorizedArrayType> &,
1453 OutVector &,
1454 const InVector &,
1455 const std::pair<unsigned int, unsigned int> &)> &face_operation,
1456 const std::function<void(
1457 const MatrixFree<dim, Number, VectorizedArrayType> &,
1458 OutVector &,
1459 const InVector &,
1460 const std::pair<unsigned int, unsigned int> &)> &boundary_operation,
1461 OutVector & dst,
1462 const InVector & src,
1463 const std::function<void(const unsigned int, const unsigned int)>
1464 &operation_before_loop,
1465 const std::function<void(const unsigned int, const unsigned int)>
1466 & operation_after_loop,
1467 const unsigned int dof_handler_index_pre_post = 0,
1468 const DataAccessOnFaces dst_vector_face_access =
1469 DataAccessOnFaces::unspecified,
1470 const DataAccessOnFaces src_vector_face_access =
1471 DataAccessOnFaces::unspecified) const;
1472
1537 template <typename CLASS, typename OutVector, typename InVector>
1538 void
1539 loop_cell_centric(void (CLASS::*cell_operation)(
1540 const MatrixFree &,
1541 OutVector &,
1542 const InVector &,
1543 const std::pair<unsigned int, unsigned int> &) const,
1544 const CLASS * owning_class,
1545 OutVector & dst,
1546 const InVector & src,
1547 const bool zero_dst_vector = false,
1548 const DataAccessOnFaces src_vector_face_access =
1549 DataAccessOnFaces::unspecified) const;
1550
1554 template <typename CLASS, typename OutVector, typename InVector>
1555 void
1556 loop_cell_centric(void (CLASS::*cell_operation)(
1557 const MatrixFree &,
1558 OutVector &,
1559 const InVector &,
1560 const std::pair<unsigned int, unsigned int> &),
1561 CLASS * owning_class,
1562 OutVector & dst,
1563 const InVector & src,
1564 const bool zero_dst_vector = false,
1565 const DataAccessOnFaces src_vector_face_access =
1566 DataAccessOnFaces::unspecified) const;
1567
1571 template <typename OutVector, typename InVector>
1572 void
1573 loop_cell_centric(
1574 const std::function<void(const MatrixFree &,
1575 OutVector &,
1576 const InVector &,
1577 const std::pair<unsigned int, unsigned int> &)>
1578 & cell_operation,
1579 OutVector & dst,
1580 const InVector & src,
1581 const bool zero_dst_vector = false,
1582 const DataAccessOnFaces src_vector_face_access =
1583 DataAccessOnFaces::unspecified) const;
1584
1592 std::pair<unsigned int, unsigned int>
1593 create_cell_subrange_hp(const std::pair<unsigned int, unsigned int> &range,
1594 const unsigned int fe_degree,
1595 const unsigned int dof_handler_index = 0) const;
1596
1603 std::pair<unsigned int, unsigned int>
1604 create_cell_subrange_hp_by_index(
1605 const std::pair<unsigned int, unsigned int> &range,
1606 const unsigned int fe_index,
1607 const unsigned int dof_handler_index = 0) const;
1608
1612 unsigned int
1613 n_active_fe_indices() const;
1614
1618 unsigned int
1619 get_cell_active_fe_index(
1620 const std::pair<unsigned int, unsigned int> range) const;
1621
1625 unsigned int
1626 get_face_active_fe_index(const std::pair<unsigned int, unsigned int> range,
1627 const bool is_interior_face = true) const;
1628
1640 template <typename T>
1641 void
1642 initialize_cell_data_vector(AlignedVector<T> &vec) const;
1643
1649 template <typename T>
1650 void
1651 initialize_face_data_vector(AlignedVector<T> &vec) const;
1652
1666 template <typename VectorType>
1667 void
1668 initialize_dof_vector(VectorType & vec,
1669 const unsigned int dof_handler_index = 0) const;
1670
1688 template <typename Number2>
1689 void
1690 initialize_dof_vector(LinearAlgebra::distributed::Vector<Number2> &vec,
1691 const unsigned int dof_handler_index = 0) const;
1692
1703 const std::shared_ptr<const Utilities::MPI::Partitioner> &
1704 get_vector_partitioner(const unsigned int dof_handler_index = 0) const;
1705
1709 const IndexSet &
1710 get_locally_owned_set(const unsigned int dof_handler_index = 0) const;
1711
1715 const IndexSet &
1716 get_ghost_set(const unsigned int dof_handler_index = 0) const;
1717
1727 const std::vector<unsigned int> &
1728 get_constrained_dofs(const unsigned int dof_handler_index = 0) const;
1729
1740 void
1741 renumber_dofs(std::vector<types::global_dof_index> &renumbering,
1742 const unsigned int dof_handler_index = 0);
1743
1753 template <int spacedim>
1754 static bool
1755 is_supported(const FiniteElement<dim, spacedim> &fe);
1756
1760 unsigned int
1761 n_components() const;
1762
1767 unsigned int
1768 n_base_elements(const unsigned int dof_handler_index) const;
1769
1777 unsigned int
1778 n_physical_cells() const;
1779
1789 unsigned int
1790 n_cell_batches() const;
1791
1798 unsigned int
1799 n_ghost_cell_batches() const;
1800
1810 unsigned int
1811 n_inner_face_batches() const;
1812
1823 unsigned int
1824 n_boundary_face_batches() const;
1825
1830 unsigned int
1831 n_ghost_inner_face_batches() const;
1832
1843 types::boundary_id
1844 get_boundary_id(const unsigned int face_batch_index) const;
1845
1850 std::array<types::boundary_id, VectorizedArrayType::size()>
1851 get_faces_by_cells_boundary_id(const unsigned int cell_batch_index,
1852 const unsigned int face_number) const;
1853
1858 const DoFHandler<dim> &
1859 get_dof_handler(const unsigned int dof_handler_index = 0) const;
1860
1867 const AffineConstraints<Number> &
1868 get_affine_constraints(const unsigned int dof_handler_index = 0) const;
1869
1882 typename DoFHandler<dim>::cell_iterator
1883 get_cell_iterator(const unsigned int cell_batch_index,
1884 const unsigned int lane_index,
1885 const unsigned int dof_handler_index = 0) const;
1886
1892 std::pair<int, int>
1893 get_cell_level_and_index(const unsigned int cell_batch_index,
1894 const unsigned int lane_index) const;
1895
1902 unsigned int
1903 get_matrix_free_cell_index(
1904 const typename Triangulation<dim>::cell_iterator &cell) const;
1905
1918 std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>
1919 get_face_iterator(const unsigned int face_batch_index,
1920 const unsigned int lane_index,
1921 const bool interior = true,
1922 const unsigned int fe_component = 0) const;
1923
1936 bool
1937 at_irregular_cell(const unsigned int cell_batch_index) const;
1938
1948 unsigned int
1949 n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const;
1950
1960 unsigned int
1961 n_active_entries_per_face_batch(const unsigned int face_batch_index) const;
1962
1966 unsigned int
1967 get_dofs_per_cell(const unsigned int dof_handler_index = 0,
1968 const unsigned int hp_active_fe_index = 0) const;
1969
1973 unsigned int
1974 get_n_q_points(const unsigned int quad_index = 0,
1975 const unsigned int hp_active_fe_index = 0) const;
1976
1981 unsigned int
1982 get_dofs_per_face(const unsigned int dof_handler_index = 0,
1983 const unsigned int hp_active_fe_index = 0) const;
1984
1989 unsigned int
1990 get_n_q_points_face(const unsigned int quad_index = 0,
1991 const unsigned int hp_active_fe_index = 0) const;
1992
1996 const Quadrature<dim> &
1997 get_quadrature(const unsigned int quad_index = 0,
1998 const unsigned int hp_active_fe_index = 0) const;
1999
2003 const Quadrature<dim - 1> &
2004 get_face_quadrature(const unsigned int quad_index = 0,
2005 const unsigned int hp_active_fe_index = 0) const;
2006
2019 unsigned int
2020 get_cell_range_category(
2021 const std::pair<unsigned int, unsigned int> cell_batch_range) const;
2022
2027 std::pair<unsigned int, unsigned int>
2028 get_face_range_category(
2029 const std::pair<unsigned int, unsigned int> face_batch_range) const;
2030
2042 unsigned int
2043 get_cell_category(const unsigned int cell_batch_index) const;
2044
2049 std::pair<unsigned int, unsigned int>
2050 get_face_category(const unsigned int face_batch_index) const;
2051
2055 bool
2056 indices_initialized() const;
2057
2062 bool
2063 mapping_initialized() const;
2064
2069 unsigned int
2070 get_mg_level() const;
2071
2076 std::size_t
2077 memory_consumption() const;
2078
2083 template <typename StreamType>
2084 void
2085 print_memory_consumption(StreamType &out) const;
2086
2091 void
2092 print(std::ostream &out) const;
2093
2105 const internal::MatrixFreeFunctions::TaskInfo &
2106 get_task_info() const;
2107
2108 /*
2109 * Return geometry-dependent information on the cells.
2110 */
2111 const internal::MatrixFreeFunctions::
2112 MappingInfo<dim, Number, VectorizedArrayType> &
2113 get_mapping_info() const;
2114
2118 const internal::MatrixFreeFunctions::DoFInfo &
2119 get_dof_info(const unsigned int dof_handler_index_component = 0) const;
2120
2124 unsigned int
2125 n_constraint_pool_entries() const;
2126
2131 const Number *
2132 constraint_pool_begin(const unsigned int pool_index) const;
2133
2139 const Number *
2140 constraint_pool_end(const unsigned int pool_index) const;
2141
2145 const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
2146 get_shape_info(const unsigned int dof_handler_index_component = 0,
2147 const unsigned int quad_index = 0,
2148 const unsigned int fe_base_element = 0,
2149 const unsigned int hp_active_fe_index = 0,
2150 const unsigned int hp_active_quad_index = 0) const;
2151
2155 const internal::MatrixFreeFunctions::FaceToCellTopology<
2156 VectorizedArrayType::size()> &
2157 get_face_info(const unsigned int face_batch_index) const;
2158
2159
2165 const Table<3, unsigned int> &
2166 get_cell_and_face_to_plain_faces() const;
2167
2181 AlignedVector<VectorizedArrayType> *
2182 acquire_scratch_data() const;
2183
2187 void
2188 release_scratch_data(const AlignedVector<VectorizedArrayType> *memory) const;
2189
2199 AlignedVector<Number> *
2200 acquire_scratch_data_non_threadsafe() const;
2201
2205 void
2206 release_scratch_data_non_threadsafe(
2207 const AlignedVector<Number> *memory) const;
2208
2211private:
2216 template <typename number2, int q_dim>
2217 void
2218 internal_reinit(
2219 const std::shared_ptr<hp::MappingCollection<dim>> & mapping,
2220 const std::vector<const DoFHandler<dim, dim> *> & dof_handlers,
2221 const std::vector<const AffineConstraints<number2> *> &constraint,
2222 const std::vector<IndexSet> & locally_owned_set,
2223 const std::vector<hp::QCollection<q_dim>> & quad,
2224 const AdditionalData & additional_data);
2225
2232 template <typename number2>
2233 void
2234 initialize_indices(
2235 const std::vector<const AffineConstraints<number2> *> &constraint,
2236 const std::vector<IndexSet> & locally_owned_set,
2237 const AdditionalData & additional_data);
2238
2242 void
2243 initialize_dof_handlers(
2244 const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
2245 const AdditionalData & additional_data);
2246
2250 std::vector<SmartPointer<const DoFHandler<dim>>> dof_handlers;
2251
2258 std::vector<SmartPointer<const AffineConstraints<Number>>> affine_constraints;
2259
2264 std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
2265
2272 std::vector<Number> constraint_pool_data;
2273
2278 std::vector<unsigned int> constraint_pool_row_index;
2279
2284 internal::MatrixFreeFunctions::MappingInfo<dim, Number, VectorizedArrayType>
2285 mapping_info;
2286
2290 Table<4, internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>>
2291 shape_info;
2292
2299 std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
2300
2305 std::vector<unsigned int> mf_cell_indices;
2306
2313 unsigned int cell_level_index_end_local;
2314
2320 internal::MatrixFreeFunctions::TaskInfo task_info;
2321
2326 internal::MatrixFreeFunctions::FaceInfo<VectorizedArrayType::size()>
2327 face_info;
2328
2332 bool indices_are_initialized;
2333
2337 bool mapping_is_initialized;
2338
2346 mutable Threads::ThreadLocalStorage<
2347 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2348 scratch_pad;
2349
2354 mutable std::list<std::pair<bool, AlignedVector<Number>>>
2355 scratch_pad_non_threadsafe;
2356
2360 unsigned int mg_level;
2361};
2362
2363
2364
2365/*----------------------- Inline functions ----------------------------------*/
2366
2367#ifndef DOXYGEN
2368
2369
2370
2371template <int dim, typename Number, typename VectorizedArrayType>
2372template <typename T>
2373inline void
2374MatrixFree<dim, Number, VectorizedArrayType>::initialize_cell_data_vector(
2375 AlignedVector<T> &vec) const
2376{
2377 vec.resize(this->n_cell_batches() + this->n_ghost_cell_batches());
2378}
2379
2380
2381
2382template <int dim, typename Number, typename VectorizedArrayType>
2383template <typename T>
2384inline void
2385MatrixFree<dim, Number, VectorizedArrayType>::initialize_face_data_vector(
2386 AlignedVector<T> &vec) const
2387{
2388 vec.resize(this->n_inner_face_batches() + this->n_boundary_face_batches() +
2389 this->n_ghost_inner_face_batches());
2390}
2391
2392
2393
2394template <int dim, typename Number, typename VectorizedArrayType>
2395template <typename VectorType>
2396inline void
2397MatrixFree<dim, Number, VectorizedArrayType>::initialize_dof_vector(
2398 VectorType & vec,
2399 const unsigned int comp) const
2400{
2401 static_assert(IsBlockVector<VectorType>::value == false,
2402 "This function is not supported for block vectors.");
2403
2404 Assert(task_info.n_procs == 1,
2405 ExcMessage("This function can only be used in serial."));
2406
2407 AssertIndexRange(comp, n_components());
2408 vec.reinit(dof_info[comp].vector_partitioner->size());
2409}
2410
2411
2412
2413template <int dim, typename Number, typename VectorizedArrayType>
2414template <typename Number2>
2415inline void
2416MatrixFree<dim, Number, VectorizedArrayType>::initialize_dof_vector(
2417 LinearAlgebra::distributed::Vector<Number2> &vec,
2418 const unsigned int comp) const
2419{
2420 AssertIndexRange(comp, n_components());
2421 vec.reinit(dof_info[comp].vector_partitioner, task_info.communicator_sm);
2422}
2423
2424
2425
2426template <int dim, typename Number, typename VectorizedArrayType>
2427inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2428MatrixFree<dim, Number, VectorizedArrayType>::get_vector_partitioner(
2429 const unsigned int comp) const
2430{
2431 AssertIndexRange(comp, n_components());
2432 return dof_info[comp].vector_partitioner;
2433}
2434
2435
2436
2437template <int dim, typename Number, typename VectorizedArrayType>
2438inline const std::vector<unsigned int> &
2439MatrixFree<dim, Number, VectorizedArrayType>::get_constrained_dofs(
2440 const unsigned int comp) const
2441{
2442 AssertIndexRange(comp, n_components());
2443 return dof_info[comp].constrained_dofs;
2444}
2445
2446
2447
2448template <int dim, typename Number, typename VectorizedArrayType>
2449inline unsigned int
2450MatrixFree<dim, Number, VectorizedArrayType>::n_components() const
2451{
2452 AssertDimension(dof_handlers.size(), dof_info.size());
2453 return dof_handlers.size();
2454}
2455
2456
2457
2458template <int dim, typename Number, typename VectorizedArrayType>
2459inline unsigned int
2460MatrixFree<dim, Number, VectorizedArrayType>::n_base_elements(
2461 const unsigned int dof_no) const
2462{
2463 AssertDimension(dof_handlers.size(), dof_info.size());
2464 AssertIndexRange(dof_no, dof_handlers.size());
2465 return dof_handlers[dof_no]->get_fe().n_base_elements();
2466}
2467
2468
2469
2470template <int dim, typename Number, typename VectorizedArrayType>
2471inline const internal::MatrixFreeFunctions::TaskInfo &
2472MatrixFree<dim, Number, VectorizedArrayType>::get_task_info() const
2473{
2474 return task_info;
2475}
2476
2477
2478
2479template <int dim, typename Number, typename VectorizedArrayType>
2480inline unsigned int
2481MatrixFree<dim, Number, VectorizedArrayType>::n_physical_cells() const
2482{
2483 return task_info.n_active_cells;
2484}
2485
2486
2487
2488template <int dim, typename Number, typename VectorizedArrayType>
2489inline unsigned int
2490MatrixFree<dim, Number, VectorizedArrayType>::n_cell_batches() const
2491{
2492 return *(task_info.cell_partition_data.end() - 2);
2493}
2494
2495
2496
2497template <int dim, typename Number, typename VectorizedArrayType>
2498inline unsigned int
2499MatrixFree<dim, Number, VectorizedArrayType>::n_ghost_cell_batches() const
2500{
2501 return *(task_info.cell_partition_data.end() - 1) -
2502 *(task_info.cell_partition_data.end() - 2);
2503}
2504
2505
2506
2507template <int dim, typename Number, typename VectorizedArrayType>
2508inline unsigned int
2509MatrixFree<dim, Number, VectorizedArrayType>::n_inner_face_batches() const
2510{
2511 if (task_info.face_partition_data.size() == 0)
2512 return 0;
2513 return task_info.face_partition_data.back();
2514}
2515
2516
2517
2518template <int dim, typename Number, typename VectorizedArrayType>
2519inline unsigned int
2520MatrixFree<dim, Number, VectorizedArrayType>::n_boundary_face_batches() const
2521{
2522 if (task_info.face_partition_data.size() == 0)
2523 return 0;
2524 return task_info.boundary_partition_data.back() -
2525 task_info.face_partition_data.back();
2526}
2527
2528
2529
2530template <int dim, typename Number, typename VectorizedArrayType>
2531inline unsigned int
2532MatrixFree<dim, Number, VectorizedArrayType>::n_ghost_inner_face_batches() const
2533{
2534 if (task_info.face_partition_data.size() == 0)
2535 return 0;
2536 return face_info.faces.size() - task_info.boundary_partition_data.back();
2537}
2538
2539
2540
2541template <int dim, typename Number, typename VectorizedArrayType>
2542inline types::boundary_id
2543MatrixFree<dim, Number, VectorizedArrayType>::get_boundary_id(
2544 const unsigned int face_batch_index) const
2545{
2546 Assert(face_batch_index >= task_info.boundary_partition_data[0] &&
2547 face_batch_index < task_info.boundary_partition_data.back(),
2548 ExcIndexRange(face_batch_index,
2549 task_info.boundary_partition_data[0],
2550 task_info.boundary_partition_data.back()));
2551 return types::boundary_id(face_info.faces[face_batch_index].exterior_face_no);
2552}
2553
2554
2555
2556template <int dim, typename Number, typename VectorizedArrayType>
2557inline std::array<types::boundary_id, VectorizedArrayType::size()>
2558MatrixFree<dim, Number, VectorizedArrayType>::get_faces_by_cells_boundary_id(
2559 const unsigned int cell_batch_index,
2560 const unsigned int face_number) const
2561{
2562 AssertIndexRange(cell_batch_index, n_cell_batches());
2563 AssertIndexRange(face_number, GeometryInfo<dim>::faces_per_cell);
2564 Assert(face_info.cell_and_face_boundary_id.size(0) >= n_cell_batches(),
2565 ExcNotInitialized());
2566 std::array<types::boundary_id, VectorizedArrayType::size()> result;
2567 result.fill(numbers::invalid_boundary_id);
2568 for (unsigned int v = 0;
2569 v < n_active_entries_per_cell_batch(cell_batch_index);
2570 ++v)
2571 result[v] =
2572 face_info.cell_and_face_boundary_id(cell_batch_index, face_number, v);
2573 return result;
2574}
2575
2576
2577
2578template <int dim, typename Number, typename VectorizedArrayType>
2579inline const internal::MatrixFreeFunctions::
2580 MappingInfo<dim, Number, VectorizedArrayType> &
2581 MatrixFree<dim, Number, VectorizedArrayType>::get_mapping_info() const
2582{
2583 return mapping_info;
2584}
2585
2586
2587
2588template <int dim, typename Number, typename VectorizedArrayType>
2589inline const internal::MatrixFreeFunctions::DoFInfo &
2590MatrixFree<dim, Number, VectorizedArrayType>::get_dof_info(
2591 const unsigned int dof_index) const
2592{
2593 AssertIndexRange(dof_index, n_components());
2594 return dof_info[dof_index];
2595}
2596
2597
2598
2599template <int dim, typename Number, typename VectorizedArrayType>
2600inline unsigned int
2601MatrixFree<dim, Number, VectorizedArrayType>::n_constraint_pool_entries() const
2602{
2603 return constraint_pool_row_index.size() - 1;
2604}
2605
2606
2607
2608template <int dim, typename Number, typename VectorizedArrayType>
2609inline const Number *
2610MatrixFree<dim, Number, VectorizedArrayType>::constraint_pool_begin(
2611 const unsigned int row) const
2612{
2613 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2614 return constraint_pool_data.empty() ?
2615 nullptr :
2616 constraint_pool_data.data() + constraint_pool_row_index[row];
2617}
2618
2619
2620
2621template <int dim, typename Number, typename VectorizedArrayType>
2622inline const Number *
2623MatrixFree<dim, Number, VectorizedArrayType>::constraint_pool_end(
2624 const unsigned int row) const
2625{
2626 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2627 return constraint_pool_data.empty() ?
2628 nullptr :
2629 constraint_pool_data.data() + constraint_pool_row_index[row + 1];
2630}
2631
2632
2633
2634template <int dim, typename Number, typename VectorizedArrayType>
2635inline std::pair<unsigned int, unsigned int>
2636MatrixFree<dim, Number, VectorizedArrayType>::create_cell_subrange_hp(
2637 const std::pair<unsigned int, unsigned int> &range,
2638 const unsigned int degree,
2639 const unsigned int dof_handler_component) const
2640{
2641 if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2642 {
2643 AssertDimension(
2644 dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2645 AssertDimension(
2646 dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2647 if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2648 return range;
2649 else
2650 return {range.second, range.second};
2651 }
2652
2653 const unsigned int fe_index =
2654 dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2655 if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2656 return {range.second, range.second};
2657 else
2658 return create_cell_subrange_hp_by_index(range,
2659 fe_index,
2660 dof_handler_component);
2661}
2662
2663
2664
2665template <int dim, typename Number, typename VectorizedArrayType>
2666inline bool
2667MatrixFree<dim, Number, VectorizedArrayType>::at_irregular_cell(
2668 const unsigned int cell_batch_index) const
2669{
2670 AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2671 return VectorizedArrayType::size() > 1 &&
2672 cell_level_index[(cell_batch_index + 1) * VectorizedArrayType::size() -
2673 1] == cell_level_index[(cell_batch_index + 1) *
2674 VectorizedArrayType::size() -
2675 2];
2676}
2677
2678
2679
2680template <int dim, typename Number, typename VectorizedArrayType>
2681unsigned int
2682MatrixFree<dim, Number, VectorizedArrayType>::n_active_fe_indices() const
2683{
2684 return shape_info.size(2);
2685}
2686
2687
2688template <int dim, typename Number, typename VectorizedArrayType>
2689unsigned int
2690MatrixFree<dim, Number, VectorizedArrayType>::get_cell_active_fe_index(
2691 const std::pair<unsigned int, unsigned int> range) const
2692{
2693 const auto &fe_indices = dof_info[0].cell_active_fe_index;
2694
2695 if (fe_indices.empty() == true)
2696 return 0;
2697
2698 const auto index = fe_indices[range.first];
2699
2700 for (unsigned int i = range.first; i < range.second; ++i)
2701 AssertDimension(index, fe_indices[i]);
2702
2703 return index;
2704}
2705
2706
2707
2708template <int dim, typename Number, typename VectorizedArrayType>
2709unsigned int
2710MatrixFree<dim, Number, VectorizedArrayType>::get_face_active_fe_index(
2711 const std::pair<unsigned int, unsigned int> range,
2712 const bool is_interior_face) const
2713{
2714 const auto &fe_indices = dof_info[0].cell_active_fe_index;
2715
2716 if (fe_indices.empty() == true)
2717 return 0;
2718
2719 if (is_interior_face)
2720 {
2721 const unsigned int index =
2722 fe_indices[face_info.faces[range.first].cells_interior[0] /
2723 VectorizedArrayType::size()];
2724
2725 for (unsigned int i = range.first; i < range.second; ++i)
2726 AssertDimension(index,
2727 fe_indices[face_info.faces[i].cells_interior[0] /
2728 VectorizedArrayType::size()]);
2729
2730 return index;
2731 }
2732 else
2733 {
2734 const unsigned int index =
2735 fe_indices[face_info.faces[range.first].cells_exterior[0] /
2736 VectorizedArrayType::size()];
2737
2738 for (unsigned int i = range.first; i < range.second; ++i)
2739 AssertDimension(index,
2740 fe_indices[face_info.faces[i].cells_exterior[0] /
2741 VectorizedArrayType::size()]);
2742
2743 return index;
2744 }
2745}
2746
2747
2748
2749template <int dim, typename Number, typename VectorizedArrayType>
2750inline unsigned int
2751MatrixFree<dim, Number, VectorizedArrayType>::n_active_entries_per_cell_batch(
2752 const unsigned int cell_batch_index) const
2753{
2754 Assert(!dof_info.empty(), ExcNotInitialized());
2755 AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2756 const std::vector<unsigned char> &n_lanes_filled =
2757 dof_info[0].n_vectorization_lanes_filled
2758 [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell];
2759 AssertIndexRange(cell_batch_index, n_lanes_filled.size());
2760
2761 return n_lanes_filled[cell_batch_index];
2762}
2763
2764
2765
2766template <int dim, typename Number, typename VectorizedArrayType>
2767inline unsigned int
2768MatrixFree<dim, Number, VectorizedArrayType>::n_active_entries_per_face_batch(
2769 const unsigned int face_batch_index) const
2770{
2771 AssertIndexRange(face_batch_index, face_info.faces.size());
2772 Assert(!dof_info.empty(), ExcNotInitialized());
2773 const std::vector<unsigned char> &n_lanes_filled =
2774 dof_info[0].n_vectorization_lanes_filled
2775 [internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior];
2776 AssertIndexRange(face_batch_index, n_lanes_filled.size());
2777 return n_lanes_filled[face_batch_index];
2778}
2779
2780
2781
2782template <int dim, typename Number, typename VectorizedArrayType>
2783inline unsigned int
2784MatrixFree<dim, Number, VectorizedArrayType>::get_dofs_per_cell(
2785 const unsigned int dof_handler_index,
2786 const unsigned int active_fe_index) const
2787{
2788 return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2789}
2790
2791
2792
2793template <int dim, typename Number, typename VectorizedArrayType>
2794inline unsigned int
2795MatrixFree<dim, Number, VectorizedArrayType>::get_n_q_points(
2796 const unsigned int quad_index,
2797 const unsigned int active_fe_index) const
2798{
2799 AssertIndexRange(quad_index, mapping_info.cell_data.size());
2800 return mapping_info.cell_data[quad_index]
2801 .descriptor[active_fe_index]
2802 .n_q_points;
2803}
2804
2805
2806
2807template <int dim, typename Number, typename VectorizedArrayType>
2808inline unsigned int
2809MatrixFree<dim, Number, VectorizedArrayType>::get_dofs_per_face(
2810 const unsigned int dof_handler_index,
2811 const unsigned int active_fe_index) const
2812{
2813 return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2814}
2815
2816
2817
2818template <int dim, typename Number, typename VectorizedArrayType>
2819inline unsigned int
2820MatrixFree<dim, Number, VectorizedArrayType>::get_n_q_points_face(
2821 const unsigned int quad_index,
2822 const unsigned int active_fe_index) const
2823{
2824 AssertIndexRange(quad_index, mapping_info.face_data.size());
2825 return mapping_info.face_data[quad_index]
2826 .descriptor[active_fe_index]
2827 .n_q_points;
2828}
2829
2830
2831
2832template <int dim, typename Number, typename VectorizedArrayType>
2833inline const IndexSet &
2834MatrixFree<dim, Number, VectorizedArrayType>::get_locally_owned_set(
2835 const unsigned int dof_handler_index) const
2836{
2837 return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2838}
2839
2840
2841
2842template <int dim, typename Number, typename VectorizedArrayType>
2843inline const IndexSet &
2844MatrixFree<dim, Number, VectorizedArrayType>::get_ghost_set(
2845 const unsigned int dof_handler_index) const
2846{
2847 return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2848}
2849
2850
2851
2852template <int dim, typename Number, typename VectorizedArrayType>
2853inline const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
2854MatrixFree<dim, Number, VectorizedArrayType>::get_shape_info(
2855 const unsigned int dof_handler_index,
2856 const unsigned int index_quad,
2857 const unsigned int index_fe,
2858 const unsigned int active_fe_index,
2859 const unsigned int active_quad_index) const
2860{
2861 AssertIndexRange(dof_handler_index, dof_info.size());
2862 const unsigned int ind =
2863 dof_info[dof_handler_index].global_base_element_offset + index_fe;
2864 AssertIndexRange(ind, shape_info.size(0));
2865 AssertIndexRange(index_quad, shape_info.size(1));
2866 AssertIndexRange(active_fe_index, shape_info.size(2));
2867 AssertIndexRange(active_quad_index, shape_info.size(3));
2868 return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2869}
2870
2871
2872
2873template <int dim, typename Number, typename VectorizedArrayType>
2874inline const internal::MatrixFreeFunctions::FaceToCellTopology<
2875 VectorizedArrayType::size()> &
2876MatrixFree<dim, Number, VectorizedArrayType>::get_face_info(
2877 const unsigned int face_batch_index) const
2878{
2879 AssertIndexRange(face_batch_index, face_info.faces.size());
2880 return face_info.faces[face_batch_index];
2881}
2882
2883
2884
2885template <int dim, typename Number, typename VectorizedArrayType>
2886inline const Table<3, unsigned int> &
2887MatrixFree<dim, Number, VectorizedArrayType>::get_cell_and_face_to_plain_faces()
2888 const
2889{
2890 return face_info.cell_and_face_to_plain_faces;
2891}
2892
2893
2894
2895template <int dim, typename Number, typename VectorizedArrayType>
2896inline const Quadrature<dim> &
2897MatrixFree<dim, Number, VectorizedArrayType>::get_quadrature(
2898 const unsigned int quad_index,
2899 const unsigned int active_fe_index) const
2900{
2901 AssertIndexRange(quad_index, mapping_info.cell_data.size());
2902 return mapping_info.cell_data[quad_index]
2903 .descriptor[active_fe_index]
2904 .quadrature;
2905}
2906
2907
2908
2909template <int dim, typename Number, typename VectorizedArrayType>
2910inline const Quadrature<dim - 1> &
2911MatrixFree<dim, Number, VectorizedArrayType>::get_face_quadrature(
2912 const unsigned int quad_index,
2913 const unsigned int active_fe_index) const
2914{
2915 AssertIndexRange(quad_index, mapping_info.face_data.size());
2916 return mapping_info.face_data[quad_index]
2917 .descriptor[active_fe_index]
2918 .quadrature;
2919}
2920
2921
2922
2923template <int dim, typename Number, typename VectorizedArrayType>
2924inline unsigned int
2925MatrixFree<dim, Number, VectorizedArrayType>::get_cell_range_category(
2926 const std::pair<unsigned int, unsigned int> range) const
2927{
2928 auto result = get_cell_category(range.first);
2929
2930 for (unsigned int i = range.first; i < range.second; ++i)
2931 result = std::max(result, get_cell_category(i));
2932
2933 return result;
2934}
2935
2936
2937
2938template <int dim, typename Number, typename VectorizedArrayType>
2939inline std::pair<unsigned int, unsigned int>
2940MatrixFree<dim, Number, VectorizedArrayType>::get_face_range_category(
2941 const std::pair<unsigned int, unsigned int> range) const
2942{
2943 auto result = get_face_category(range.first);
2944
2945 for (unsigned int i = range.first; i < range.second; ++i)
2946 {
2947 result.first = std::max(result.first, get_face_category(i).first);
2948 result.second = std::max(result.second, get_face_category(i).second);
2949 }
2950
2951 return result;
2952}
2953
2954
2955
2956template <int dim, typename Number, typename VectorizedArrayType>
2957inline unsigned int
2958MatrixFree<dim, Number, VectorizedArrayType>::get_cell_category(
2959 const unsigned int cell_batch_index) const
2960{
2961 AssertIndexRange(0, dof_info.size());
2962 AssertIndexRange(cell_batch_index, dof_info[0].cell_active_fe_index.size());
2963 if (dof_info[0].cell_active_fe_index.empty())
2964 return 0;
2965 else
2966 return dof_info[0].cell_active_fe_index[cell_batch_index];
2967}
2968
2969
2970
2971template <int dim, typename Number, typename VectorizedArrayType>
2972inline std::pair<unsigned int, unsigned int>
2973MatrixFree<dim, Number, VectorizedArrayType>::get_face_category(
2974 const unsigned int face_batch_index) const
2975{
2976 AssertIndexRange(face_batch_index, face_info.faces.size());
2977 if (dof_info[0].cell_active_fe_index.empty())
2978 return std::make_pair(0U, 0U);
2979
2980 std::pair<unsigned int, unsigned int> result = std::make_pair(0U, 0U);
2981 for (unsigned int v = 0;
2982 v < VectorizedArrayType::size() &&
2983 face_info.faces[face_batch_index].cells_interior[v] !=
2984 numbers::invalid_unsigned_int;
2985 ++v)
2986 result.first = std::max(
2987 result.first,
2988 dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
2989 .cells_interior[v] /
2990 VectorizedArrayType::size()]);
2991 if (face_info.faces[face_batch_index].cells_exterior[0] !=
2992 numbers::invalid_unsigned_int)
2993 for (unsigned int v = 0;
2994 v < VectorizedArrayType::size() &&
2995 face_info.faces[face_batch_index].cells_exterior[v] !=
2996 numbers::invalid_unsigned_int;
2997 ++v)
2998 result.second = std::max(
2999 result.second,
3000 dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
3001 .cells_exterior[v] /
3002 VectorizedArrayType::size()]);
3003 else
3004 result.second = numbers::invalid_unsigned_int;
3005 return result;
3006}
3007
3008
3009
3010template <int dim, typename Number, typename VectorizedArrayType>
3011inline bool
3012MatrixFree<dim, Number, VectorizedArrayType>::indices_initialized() const
3013{
3014 return indices_are_initialized;
3015}
3016
3017
3018
3019template <int dim, typename Number, typename VectorizedArrayType>
3020inline bool
3021MatrixFree<dim, Number, VectorizedArrayType>::mapping_initialized() const
3022{
3023 return mapping_is_initialized;
3024}
3025
3026
3027template <int dim, typename Number, typename VectorizedArrayType>
3028inline unsigned int
3029MatrixFree<dim, Number, VectorizedArrayType>::get_mg_level() const
3030{
3031 return mg_level;
3032}
3033
3034
3035
3036template <int dim, typename Number, typename VectorizedArrayType>
3037AlignedVector<VectorizedArrayType> *
3038MatrixFree<dim, Number, VectorizedArrayType>::acquire_scratch_data() const
3039{
3040 using list_type =
3041 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3042 list_type &data = scratch_pad.get();
3043 for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
3044 if (it->first == false)
3045 {
3046 it->first = true;
3047 return &it->second;
3048 }
3049 data.emplace_front(true, AlignedVector<VectorizedArrayType>());
3050 return &data.front().second;
3051}
3052
3053
3054
3055template <int dim, typename Number, typename VectorizedArrayType>
3056void
3057MatrixFree<dim, Number, VectorizedArrayType>::release_scratch_data(
3058 const AlignedVector<VectorizedArrayType> *scratch) const
3059{
3060 using list_type =
3061 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3062 list_type &data = scratch_pad.get();
3063 for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
3064 if (&it->second == scratch)
3065 {
3066 Assert(it->first == true, ExcInternalError());
3067 it->first = false;
3068 return;
3069 }
3070 AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
3071}
3072
3073
3074
3075template <int dim, typename Number, typename VectorizedArrayType>
3076AlignedVector<Number> *
3077MatrixFree<dim, Number, VectorizedArrayType>::
3078 acquire_scratch_data_non_threadsafe() const
3079{
3080 for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
3081 scratch_pad_non_threadsafe.begin();
3082 it != scratch_pad_non_threadsafe.end();
3083 ++it)
3084 if (it->first == false)
3085 {
3086 it->first = true;
3087 return &it->second;
3088 }
3089 scratch_pad_non_threadsafe.push_front(
3090 std::make_pair(true, AlignedVector<Number>()));
3091 return &scratch_pad_non_threadsafe.front().second;
3092}
3093
3094
3095
3096template <int dim, typename Number, typename VectorizedArrayType>
3097void
3098MatrixFree<dim, Number, VectorizedArrayType>::
3099 release_scratch_data_non_threadsafe(
3100 const AlignedVector<Number> *scratch) const
3101{
3102 for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
3103 scratch_pad_non_threadsafe.begin();
3104 it != scratch_pad_non_threadsafe.end();
3105 ++it)
3106 if (&it->second == scratch)
3107 {
3108 Assert(it->first == true, ExcInternalError());
3109 it->first = false;
3110 return;
3111 }
3112 AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
3113}
3114
3115
3116
3117// ------------------------------ reinit functions ---------------------------
3118
3119namespace internal
3120{
3121 namespace MatrixFreeImplementation
3122 {
3123 template <int dim, int spacedim>
3124 inline std::vector<IndexSet>
3125 extract_locally_owned_index_sets(
3126 const std::vector<const DoFHandler<dim, spacedim> *> &dofh,
3127 const unsigned int level)
3128 {
3129 std::vector<IndexSet> locally_owned_set;
3130 locally_owned_set.reserve(dofh.size());
3131 for (unsigned int j = 0; j < dofh.size(); ++j)
3132 if (level == numbers::invalid_unsigned_int)
3133 locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
3134 else
3135 locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
3136 return locally_owned_set;
3137 }
3138 } // namespace MatrixFreeImplementation
3139} // namespace internal
3140
3141
3142
3143template <int dim, typename Number, typename VectorizedArrayType>
3144template <typename QuadratureType, typename number2, typename MappingType>
3145void
3146MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3147 const MappingType & mapping,
3148 const DoFHandler<dim> & dof_handler,
3149 const AffineConstraints<number2> &constraints_in,
3150 const QuadratureType & quad,
3151 const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
3152 &additional_data)
3153{
3154 std::vector<const DoFHandler<dim, dim> *> dof_handlers;
3155 std::vector<const AffineConstraints<number2> *> constraints;
3156
3157 dof_handlers.push_back(&dof_handler);
3158 constraints.push_back(&constraints_in);
3159
3160 std::vector<IndexSet> locally_owned_sets =
3161 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3162 dof_handlers, additional_data.mg_level);
3163
3164 std::vector<hp::QCollection<dim>> quad_hp;
3165 quad_hp.emplace_back(quad);
3166
3167 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3168 dof_handlers,
3169 constraints,
3170 locally_owned_sets,
3171 quad_hp,
3172 additional_data);
3173}
3174
3175
3176
3177template <int dim, typename Number, typename VectorizedArrayType>
3178template <typename QuadratureType, typename number2, typename MappingType>
3179void
3180MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3181 const MappingType & mapping,
3182 const std::vector<const DoFHandler<dim> *> & dof_handler,
3183 const std::vector<const AffineConstraints<number2> *> &constraint,
3184 const QuadratureType & quad,
3185 const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
3186 &additional_data)
3187{
3188 std::vector<IndexSet> locally_owned_set =
3189 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3190 dof_handler, additional_data.mg_level);
3191 std::vector<hp::QCollection<dim>> quad_hp;
3192 quad_hp.emplace_back(quad);
3193
3194 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3195 dof_handler,
3196 constraint,
3197 locally_owned_set,
3198 quad_hp,
3199 additional_data);
3200}
3201
3202
3203
3204template <int dim, typename Number, typename VectorizedArrayType>
3205template <typename QuadratureType, typename number2, typename MappingType>
3206void
3207MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3208 const MappingType & mapping,
3209 const std::vector<const DoFHandler<dim> *> & dof_handler,
3210 const std::vector<const AffineConstraints<number2> *> &constraint,
3211 const std::vector<QuadratureType> & quad,
3212 const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
3213 &additional_data)
3214{
3215 std::vector<IndexSet> locally_owned_set =
3216 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3217 dof_handler, additional_data.mg_level);
3218 std::vector<hp::QCollection<dim>> quad_hp;
3219 for (unsigned int q = 0; q < quad.size(); ++q)
3220 quad_hp.emplace_back(quad[q]);
3221
3222 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3223 dof_handler,
3224 constraint,
3225 locally_owned_set,
3226 quad_hp,
3227 additional_data);
3228}
3229
3230
3231
3232// ------------------------------ implementation of loops --------------------
3233
3234// internal helper functions that define how to call MPI data exchange
3235// functions: for generic vectors, do nothing at all. For distributed vectors,
3236// call update_ghost_values_start function and so on. If we have collections
3237// of vectors, just do the individual functions of the components. In order to
3238// keep ghost values consistent (whether we are in read or write mode), we
3239// also reset the values at the end. the whole situation is a bit complicated
3240// by the fact that we need to treat block vectors differently, which use some
3241// additional helper functions to select the blocks and template magic.
3242namespace internal
3243{
3247 template <int dim, typename Number, typename VectorizedArrayType>
3248 struct VectorDataExchange
3249 {
3250 // A shift for the MPI messages to reduce the risk for accidental
3251 // interaction with other open communications that a user program might
3252 // set up (parallel vectors support unfinished communication). We let
3253 // the other vectors use the first 20 assigned numbers and start the
3254 // matrix-free communication.
3255 static constexpr unsigned int channel_shift = 20;
3256
3257
3258
3263 VectorDataExchange(
3264 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3265 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3266 DataAccessOnFaces vector_face_access,
3267 const unsigned int n_components)
3268 : matrix_free(matrix_free)
3269 , vector_face_access(
3270 matrix_free.get_task_info().face_partition_data.empty() ?
3271 ::MatrixFree<dim, Number, VectorizedArrayType>::
3272 DataAccessOnFaces::unspecified :
3273 vector_face_access)
3274 , ghosts_were_set(false)
3275# ifdef DEAL_II_WITH_MPI
3276 , tmp_data(n_components)
3277 , requests(n_components)
3278# endif
3279 {
3280 (void)n_components;
3281 if (this->vector_face_access !=
3282 ::MatrixFree<dim, Number, VectorizedArrayType>::
3283 DataAccessOnFaces::unspecified)
3284 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3285 AssertDimension(
3286 matrix_free.get_dof_info(c).vector_exchanger_face_variants.size(),
3287 5);
3288 }
3289
3290
3291
3295 ~VectorDataExchange() // NOLINT
3296 {
3297# ifdef DEAL_II_WITH_MPI
3298 for (unsigned int i = 0; i < tmp_data.size(); ++i)
3299 if (tmp_data[i] != nullptr)
3300 matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3301# endif
3302 }
3303
3304
3305
3310 template <typename VectorType>
3311 unsigned int
3312 find_vector_in_mf(const VectorType &vec,
3313 const bool check_global_compatibility = true) const
3314 {
3315 // case 1: vector was set up with MatrixFree::initialize_dof_vector()
3316 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3317 if (vec.get_partitioner().get() ==
3318 matrix_free.get_dof_info(c).vector_partitioner.get())
3319 return c;
3320
3321 // case 2: user provided own partitioner (compatibility mode)
3322 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3323 if (check_global_compatibility ?
3324 vec.get_partitioner()->is_globally_compatible(
3325 *matrix_free.get_dof_info(c).vector_partitioner) :
3326 vec.get_partitioner()->is_compatible(
3327 *matrix_free.get_dof_info(c).vector_partitioner))
3328 return c;
3329
3330 Assert(false,
3331 ExcNotImplemented("Could not find partitioner that fits vector"));
3332
3333 return numbers::invalid_unsigned_int;
3334 }
3335
3336
3337
3342 const internal::MatrixFreeFunctions::VectorDataExchange::Base &
3343 get_partitioner(const unsigned int mf_component) const
3344 {
3345 AssertDimension(matrix_free.get_dof_info(mf_component)
3346 .vector_exchanger_face_variants.size(),
3347 5);
3348 if (vector_face_access ==
3349 ::MatrixFree<dim, Number, VectorizedArrayType>::
3350 DataAccessOnFaces::none)
3351 return *matrix_free.get_dof_info(mf_component)
3352 .vector_exchanger_face_variants[0];
3353 else if (vector_face_access ==
3354 ::MatrixFree<dim, Number, VectorizedArrayType>::
3355 DataAccessOnFaces::values)
3356 return *matrix_free.get_dof_info(mf_component)
3357 .vector_exchanger_face_variants[1];
3358 else if (vector_face_access ==
3359 ::MatrixFree<dim, Number, VectorizedArrayType>::
3360 DataAccessOnFaces::gradients)
3361 return *matrix_free.get_dof_info(mf_component)
3362 .vector_exchanger_face_variants[2];
3363 else if (vector_face_access ==
3364 ::MatrixFree<dim, Number, VectorizedArrayType>::
3365 DataAccessOnFaces::values_all_faces)
3366 return *matrix_free.get_dof_info(mf_component)
3367 .vector_exchanger_face_variants[3];
3368 else if (vector_face_access ==
3369 ::MatrixFree<dim, Number, VectorizedArrayType>::
3370 DataAccessOnFaces::gradients_all_faces)
3371 return *matrix_free.get_dof_info(mf_component)
3372 .vector_exchanger_face_variants[4];
3373 else
3374 return *matrix_free.get_dof_info(mf_component).vector_exchanger.get();
3375 }
3376
3377
3378
3382 template <typename VectorType,
3383 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3384 * = nullptr>
3385 void
3386 update_ghost_values_start(const unsigned int /*component_in_block_vector*/,
3387 const VectorType & /*vec*/)
3388 {}
3389
3390
3395 template <typename VectorType,
3396 std::enable_if_t<!has_update_ghost_values_start<VectorType> &&
3397 !is_not_parallel_vector<VectorType>,
3398 VectorType> * = nullptr>
3399 void
3400 update_ghost_values_start(const unsigned int component_in_block_vector,
3401 const VectorType & vec)
3402 {
3403 (void)component_in_block_vector;
3404 const bool ghosts_set = vec.has_ghost_elements();
3405
3406 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3407 ghosts_set == false,
3408 ExcNotImplemented());
3409
3410 if (ghosts_set)
3411 ghosts_were_set = true;
3412
3413 vec.update_ghost_values();
3414 }
3415
3416
3417
3423 template <typename VectorType,
3424 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3425 !has_exchange_on_subset<VectorType>,
3426 VectorType> * = nullptr>
3427 void
3428 update_ghost_values_start(const unsigned int component_in_block_vector,
3429 const VectorType & vec)
3430 {
3431 (void)component_in_block_vector;
3432 const bool ghosts_set = vec.has_ghost_elements();
3433
3434 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3435 ghosts_set == false,
3436 ExcNotImplemented());
3437
3438 if (ghosts_set)
3439 ghosts_were_set = true;
3440
3441 vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3442 }
3443
3444
3445
3452 template <typename VectorType,
3453 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3454 has_exchange_on_subset<VectorType>,
3455 VectorType> * = nullptr>
3456 void
3457 update_ghost_values_start(const unsigned int component_in_block_vector,
3458 const VectorType & vec)
3459 {
3460 static_assert(
3461 std::is_same<Number, typename VectorType::value_type>::value,
3462 "Type mismatch between VectorType and VectorDataExchange");
3463 (void)component_in_block_vector;
3464 const bool ghosts_set = vec.has_ghost_elements();
3465
3466 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3467 ghosts_set == false,
3468 ExcNotImplemented());
3469
3470 if (ghosts_set)
3471 ghosts_were_set = true;
3472
3473 if (vec.size() != 0)
3474 {
3475# ifdef DEAL_II_WITH_MPI
3476 const unsigned int mf_component = find_vector_in_mf(vec);
3477
3478 const auto &part = get_partitioner(mf_component);
3479
3480 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3481 part.n_import_sm_procs() == 0)
3482 return;
3483
3484 tmp_data[component_in_block_vector] =
3485 matrix_free.acquire_scratch_data_non_threadsafe();
3486 tmp_data[component_in_block_vector]->resize_fast(
3487 part.n_import_indices());
3488 AssertDimension(requests.size(), tmp_data.size());
3489
3490 part.export_to_ghosted_array_start(
3491 component_in_block_vector * 2 + channel_shift,
3492 ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3493 vec.shared_vector_data(),
3494 ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3495 part.locally_owned_size(),
3496 matrix_free.get_dof_info(mf_component)
3497 .vector_partitioner->n_ghost_indices()),
3498 ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3499 part.n_import_indices()),
3500 this->requests[component_in_block_vector]);
3501# endif
3502 }
3503 }
3504
3505
3506
3511 template <typename VectorType,
3512 std::enable_if_t<!has_update_ghost_values_start<VectorType>,
3513 VectorType> * = nullptr>
3514 void
3515 update_ghost_values_finish(const unsigned int /*component_in_block_vector*/,
3516 const VectorType & /*vec*/)
3517 {}
3518
3519
3520
3526 template <typename VectorType,
3527 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3528 !has_exchange_on_subset<VectorType>,
3529 VectorType> * = nullptr>
3530 void
3531 update_ghost_values_finish(const unsigned int component_in_block_vector,
3532 const VectorType & vec)
3533 {
3534 (void)component_in_block_vector;
3535 vec.update_ghost_values_finish();
3536 }
3537
3538
3539
3546 template <typename VectorType,
3547 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3548 has_exchange_on_subset<VectorType>,
3549 VectorType> * = nullptr>
3550 void
3551 update_ghost_values_finish(const unsigned int component_in_block_vector,
3552 const VectorType & vec)
3553 {
3554 static_assert(
3555 std::is_same<Number, typename VectorType::value_type>::value,
3556 "Type mismatch between VectorType and VectorDataExchange");
3557 (void)component_in_block_vector;
3558
3559 if (vec.size() != 0)
3560 {
3561# ifdef DEAL_II_WITH_MPI
3562 AssertIndexRange(component_in_block_vector, tmp_data.size());
3563 AssertDimension(requests.size(), tmp_data.size());
3564
3565 const unsigned int mf_component = find_vector_in_mf(vec);
3566
3567 const auto &part = get_partitioner(mf_component);
3568
3569 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3570 part.n_import_sm_procs() != 0)
3571 {
3572 part.export_to_ghosted_array_finish(
3573 ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3574 vec.shared_vector_data(),
3575 ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3576 part.locally_owned_size(),
3577 matrix_free.get_dof_info(mf_component)
3578 .vector_partitioner->n_ghost_indices()),
3579 this->requests[component_in_block_vector]);
3580
3581 matrix_free.release_scratch_data_non_threadsafe(
3582 tmp_data[component_in_block_vector]);
3583 tmp_data[component_in_block_vector] = nullptr;
3584 }
3585# endif
3586 }
3587 // let vector know that ghosts are being updated and we can read from
3588 // them
3589 vec.set_ghost_state(true);
3590 }
3591
3592
3593
3597 template <typename VectorType,
3598 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3599 * = nullptr>
3600 void
3601 compress_start(const unsigned int /*component_in_block_vector*/,
3602 VectorType & /*vec*/)
3603 {}
3604
3605
3606
3611 template <typename VectorType,
3612 std::enable_if_t<!has_compress_start<VectorType> &&
3613 !is_not_parallel_vector<VectorType>,
3614 VectorType> * = nullptr>
3615 void
3616 compress_start(const unsigned int component_in_block_vector,
3617 VectorType & vec)
3618 {
3619 (void)component_in_block_vector;
3620 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3621 vec.compress(VectorOperation::add);
3622 }
3623
3624
3625
3631 template <typename VectorType,
3632 std::enable_if_t<has_compress_start<VectorType> &&
3633 !has_exchange_on_subset<VectorType>,
3634 VectorType> * = nullptr>
3635 void
3636 compress_start(const unsigned int component_in_block_vector,
3637 VectorType & vec)
3638 {
3639 (void)component_in_block_vector;
3640 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3641 vec.compress_start(component_in_block_vector + channel_shift);
3642 }
3643
3644
3645
3652 template <typename VectorType,
3653 std::enable_if_t<has_compress_start<VectorType> &&
3654 has_exchange_on_subset<VectorType>,
3655 VectorType> * = nullptr>
3656 void
3657 compress_start(const unsigned int component_in_block_vector,
3658 VectorType & vec)
3659 {
3660 static_assert(
3661 std::is_same<Number, typename VectorType::value_type>::value,
3662 "Type mismatch between VectorType and VectorDataExchange");
3663 (void)component_in_block_vector;
3664 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3665
3666 if (vec.size() != 0)
3667 {
3668# ifdef DEAL_II_WITH_MPI
3669 const unsigned int mf_component = find_vector_in_mf(vec);
3670
3671 const auto &part = get_partitioner(mf_component);
3672
3673 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3674 part.n_import_sm_procs() == 0)
3675 return;
3676
3677 tmp_data[component_in_block_vector] =
3678 matrix_free.acquire_scratch_data_non_threadsafe();
3679 tmp_data[component_in_block_vector]->resize_fast(
3680 part.n_import_indices());
3681 AssertDimension(requests.size(), tmp_data.size());
3682
3683 part.import_from_ghosted_array_start(
3684 VectorOperation::add,
3685 component_in_block_vector * 2 + channel_shift,
3686 ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3687 vec.shared_vector_data(),
3688 ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3689 matrix_free.get_dof_info(mf_component)
3690 .vector_partitioner->n_ghost_indices()),
3691 ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3692 part.n_import_indices()),
3693 this->requests[component_in_block_vector]);
3694# endif
3695 }
3696 }
3697
3698
3699
3704 template <
3705 typename VectorType,
3706 std::enable_if_t<!has_compress_start<VectorType>, VectorType> * = nullptr>
3707 void
3708 compress_finish(const unsigned int /*component_in_block_vector*/,
3709 VectorType & /*vec*/)
3710 {}
3711
3712
3713
3719 template <typename VectorType,
3720 std::enable_if_t<has_compress_start<VectorType> &&
3721 !has_exchange_on_subset<VectorType>,
3722 VectorType> * = nullptr>
3723 void
3724 compress_finish(const unsigned int component_in_block_vector,
3725 VectorType & vec)
3726 {
3727 (void)component_in_block_vector;
3728 vec.compress_finish(VectorOperation::add);
3729 }
3730
3731
3732
3739 template <typename VectorType,
3740 std::enable_if_t<has_compress_start<VectorType> &&
3741 has_exchange_on_subset<VectorType>,
3742 VectorType> * = nullptr>
3743 void
3744 compress_finish(const unsigned int component_in_block_vector,
3745 VectorType & vec)
3746 {
3747 static_assert(
3748 std::is_same<Number, typename VectorType::value_type>::value,
3749 "Type mismatch between VectorType and VectorDataExchange");
3750 (void)component_in_block_vector;
3751 if (vec.size() != 0)
3752 {
3753# ifdef DEAL_II_WITH_MPI
3754 AssertIndexRange(component_in_block_vector, tmp_data.size());
3755 AssertDimension(requests.size(), tmp_data.size());
3756
3757 const unsigned int mf_component = find_vector_in_mf(vec);
3758
3759 const auto &part = get_partitioner(mf_component);
3760
3761 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3762 part.n_import_sm_procs() != 0)
3763 {
3764 part.import_from_ghosted_array_finish(
3765 VectorOperation::add,
3766 ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3767 vec.shared_vector_data(),
3768 ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3769 matrix_free.get_dof_info(mf_component)
3770 .vector_partitioner->n_ghost_indices()),
3771 ArrayView<const Number>(
3772 tmp_data[component_in_block_vector]->begin(),
3773 part.n_import_indices()),
3774 this->requests[component_in_block_vector]);
3775
3776 matrix_free.release_scratch_data_non_threadsafe(
3777 tmp_data[component_in_block_vector]);
3778 tmp_data[component_in_block_vector] = nullptr;
3779 }
3780
3781 if (Utilities::MPI::job_supports_mpi())
3782 {
3783 const int ierr =
3784 MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3785 AssertThrowMPI(ierr);
3786 }
3787# endif
3788 }
3789 }
3790
3791
3792
3796 template <typename VectorType,
3797 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3798 * = nullptr>
3799 void
3800 reset_ghost_values(const VectorType & /*vec*/) const
3801 {}
3802
3803
3804
3809 template <typename VectorType,
3810 std::enable_if_t<!has_exchange_on_subset<VectorType> &&
3811 !is_not_parallel_vector<VectorType>,
3812 VectorType> * = nullptr>
3813 void
3814 reset_ghost_values(const VectorType &vec) const
3815 {
3816 if (ghosts_were_set == true)
3817 return;
3818
3819 vec.zero_out_ghost_values();
3820 }
3821
3822
3823
3829 template <typename VectorType,
3830 std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3831 * = nullptr>
3832 void
3833 reset_ghost_values(const VectorType &vec) const
3834 {
3835 static_assert(
3836 std::is_same<Number, typename VectorType::value_type>::value,
3837 "Type mismatch between VectorType and VectorDataExchange");
3838 if (ghosts_were_set == true)
3839 return;
3840
3841 if (vec.size() != 0)
3842 {
3843# ifdef DEAL_II_WITH_MPI
3844 AssertDimension(requests.size(), tmp_data.size());
3845
3846 const unsigned int mf_component = find_vector_in_mf(vec);
3847
3848 const auto &part = get_partitioner(mf_component);
3849
3850 if (part.n_ghost_indices() > 0)
3851 {
3852 part.reset_ghost_values(ArrayView<Number>(
3853 const_cast<LinearAlgebra::distributed::Vector<Number> &>(vec)
3854 .begin() +
3855 part.locally_owned_size(),
3856 matrix_free.get_dof_info(mf_component)
3857 .vector_partitioner->n_ghost_indices()));
3858 }
3859
3860# endif
3861 }
3862 // let vector know that it's not ghosted anymore
3863 vec.set_ghost_state(false);
3864 }
3865
3866
3867
3873 template <typename VectorType,
3874 std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3875 * = nullptr>
3876 void
3877 zero_vector_region(const unsigned int range_index, VectorType &vec) const
3878 {
3879 static_assert(
3880 std::is_same<Number, typename VectorType::value_type>::value,
3881 "Type mismatch between VectorType and VectorDataExchange");
3882 if (range_index == numbers::invalid_unsigned_int)
3883 vec = Number();
3884 else
3885 {
3886 const unsigned int mf_component = find_vector_in_mf(vec, false);
3887 const internal::MatrixFreeFunctions::DoFInfo &dof_info =
3888 matrix_free.get_dof_info(mf_component);
3889 Assert(dof_info.vector_zero_range_list_index.empty() == false,
3890 ExcNotInitialized());
3891
3892 Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3893 ExcInternalError());
3894 AssertIndexRange(range_index,
3895 dof_info.vector_zero_range_list_index.size() - 1);
3896 for (unsigned int id =
3897 dof_info.vector_zero_range_list_index[range_index];
3898 id != dof_info.vector_zero_range_list_index[range_index + 1];
3899 ++id)
3900 std::memset(vec.begin() + dof_info.vector_zero_range_list[id].first,
3901 0,
3902 (dof_info.vector_zero_range_list[id].second -
3903 dof_info.vector_zero_range_list[id].first) *
3904 sizeof(Number));
3905 }
3906 }
3907
3908
3909
3915 template <typename VectorType,
3916 std::enable_if_t<!has_exchange_on_subset<VectorType>, VectorType>
3917 * = nullptr,
3918 typename VectorType::value_type * = nullptr>
3919 void
3920 zero_vector_region(const unsigned int range_index, VectorType &vec) const
3921 {
3922 if (range_index == numbers::invalid_unsigned_int || range_index == 0)
3923 vec = typename VectorType::value_type();
3924 }
3925
3926
3927
3932 void
3933 zero_vector_region(const unsigned int, ...) const
3934 {
3935 Assert(false,
3936 ExcNotImplemented("Zeroing is only implemented for vector types "
3937 "which provide VectorType::value_type"));
3938 }
3939
3940
3941
3942 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
3943 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3944 DataAccessOnFaces vector_face_access;
3945 bool ghosts_were_set;
3946# ifdef DEAL_II_WITH_MPI
3947 std::vector<AlignedVector<Number> *> tmp_data;
3948 std::vector<std::vector<MPI_Request>> requests;
3949# endif
3950 }; // VectorDataExchange
3951
3952 template <typename VectorStruct>
3953 unsigned int
3954 n_components(const VectorStruct &vec);
3955
3956 template <typename VectorStruct>
3957 unsigned int
3958 n_components_block(const VectorStruct &vec,
3959 std::integral_constant<bool, true>)
3960 {
3961 unsigned int components = 0;
3962 for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
3963 components += n_components(vec.block(bl));
3964 return components;
3965 }
3966
3967 template <typename VectorStruct>
3968 unsigned int
3969 n_components_block(const VectorStruct &, std::integral_constant<bool, false>)
3970 {
3971 return 1;
3972 }
3973
3974 template <typename VectorStruct>
3975 unsigned int
3976 n_components(const VectorStruct &vec)
3977 {
3978 return n_components_block(
3979 vec, std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3980 }
3981
3982 template <typename VectorStruct>
3983 inline unsigned int
3984 n_components(const std::vector<VectorStruct> &vec)
3985 {
3986 unsigned int components = 0;
3987 for (unsigned int comp = 0; comp < vec.size(); ++comp)
3988 components += n_components_block(
3989 vec[comp],
3990 std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
3991 return components;
3992 }
3993
3994 template <typename VectorStruct>
3995 inline unsigned int
3996 n_components(const std::vector<VectorStruct *> &vec)
3997 {
3998 unsigned int components = 0;
3999 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4000 components += n_components_block(
4001 *vec[comp],
4002 std::integral_constant<bool, IsBlockVector<VectorStruct>::value>());
4003 return components;
4004 }
4005
4006
4007
4008 // A helper function to identify block vectors with many components where we
4009 // should not try to overlap computations and communication because there
4010 // would be too many outstanding communication requests.
4011
4012 // default value for vectors that do not have communication_block_size
4013 template <typename VectorStruct,
4014 std::enable_if_t<!has_communication_block_size<VectorStruct>,
4015 VectorStruct> * = nullptr>
4016 constexpr unsigned int
4017 get_communication_block_size(const VectorStruct &)
4018 {
4019 return numbers::invalid_unsigned_int;
4020 }
4021
4022
4023
4024 template <typename VectorStruct,
4025 std::enable_if_t<has_communication_block_size<VectorStruct>,
4026 VectorStruct> * = nullptr>
4027 constexpr unsigned int
4028 get_communication_block_size(const VectorStruct &)
4029 {
4030 return VectorStruct::communication_block_size;
4031 }
4032
4033
4034
4035 template <typename VectorType,
4036 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType> * =
4037 nullptr>
4038 bool
4039 has_ghost_elements(const VectorType &vec)
4040 {
4041 (void)vec;
4042 return false;
4043 }
4044
4045
4046
4047 template <typename VectorType,
4048 std::enable_if_t<!is_not_parallel_vector<VectorType>, VectorType>
4049 * = nullptr>
4050 bool
4051 has_ghost_elements(const VectorType &vec)
4052 {
4053 return vec.has_ghost_elements();
4054 }
4055
4056
4057
4058 // --------------------------------------------------------------------------
4059 // below we have wrappers to distinguish between block and non-block vectors.
4060 // --------------------------------------------------------------------------
4061
4062 //
4063 // update_ghost_values_start
4064 //
4065
4066 // update_ghost_values for block vectors
4067 template <int dim,
4068 typename VectorStruct,
4069 typename Number,
4070 typename VectorizedArrayType,
4071 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4072 * = nullptr>
4073 void
4074 update_ghost_values_start(
4075 const VectorStruct & vec,
4076 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4077 const unsigned int channel = 0)
4078 {
4079 if (get_communication_block_size(vec) < vec.n_blocks())
4080 {
4081 const bool ghosts_set = vec.has_ghost_elements();
4082
4083 Assert(exchanger.matrix_free.get_task_info()
4084 .allow_ghosted_vectors_in_loops ||
4085 ghosts_set == false,
4086 ExcNotImplemented());
4087
4088 if (ghosts_set)
4089 exchanger.ghosts_were_set = true;
4090
4091 vec.update_ghost_values();
4092 }
4093 else
4094 {
4095 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4096 update_ghost_values_start(vec.block(i), exchanger, channel + i);
4097 }
4098 }
4099
4100
4101
4102 // update_ghost_values for non-block vectors
4103 template <int dim,
4104 typename VectorStruct,
4105 typename Number,
4106 typename VectorizedArrayType,
4107 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4108 * = nullptr>
4109 void
4110 update_ghost_values_start(
4111 const VectorStruct & vec,
4112 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4113 const unsigned int channel = 0)
4114 {
4115 exchanger.update_ghost_values_start(channel, vec);
4116 }
4117
4118
4119
4120 // update_ghost_values_start() for vector of vectors
4121 template <int dim,
4122 typename VectorStruct,
4123 typename Number,
4124 typename VectorizedArrayType>
4125 inline void
4126 update_ghost_values_start(
4127 const std::vector<VectorStruct> & vec,
4128 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4129 {
4130 unsigned int component_index = 0;
4131 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4132 {
4133 update_ghost_values_start(vec[comp], exchanger, component_index);
4134 component_index += n_components(vec[comp]);
4135 }
4136 }
4137
4138
4139
4140 // update_ghost_values_start() for vector of pointers to vectors
4141 template <int dim,
4142 typename VectorStruct,
4143 typename Number,
4144 typename VectorizedArrayType>
4145 inline void
4146 update_ghost_values_start(
4147 const std::vector<VectorStruct *> & vec,
4148 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4149 {
4150 unsigned int component_index = 0;
4151 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4152 {
4153 update_ghost_values_start(*vec[comp], exchanger, component_index);
4154 component_index += n_components(*vec[comp]);
4155 }
4156 }
4157
4158
4159
4160 //
4161 // update_ghost_values_finish
4162 //
4163
4164 // for block vectors
4165 template <int dim,
4166 typename VectorStruct,
4167 typename Number,
4168 typename VectorizedArrayType,
4169 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4170 * = nullptr>
4171 void
4172 update_ghost_values_finish(
4173 const VectorStruct & vec,
4174 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4175 const unsigned int channel = 0)
4176 {
4177 if (get_communication_block_size(vec) < vec.n_blocks())
4178 {
4179 // do nothing, everything has already been completed in the _start()
4180 // call
4181 }
4182 else
4183 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4184 update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4185 }
4186
4187
4188
4189 // for non-block vectors
4190 template <int dim,
4191 typename VectorStruct,
4192 typename Number,
4193 typename VectorizedArrayType,
4194 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4195 * = nullptr>
4196 void
4197 update_ghost_values_finish(
4198 const VectorStruct & vec,
4199 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4200 const unsigned int channel = 0)
4201 {
4202 exchanger.update_ghost_values_finish(channel, vec);
4203 }
4204
4205
4206
4207 // for vector of vectors
4208 template <int dim,
4209 typename VectorStruct,
4210 typename Number,
4211 typename VectorizedArrayType>
4212 inline void
4213 update_ghost_values_finish(
4214 const std::vector<VectorStruct> & vec,
4215 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4216 {
4217 unsigned int component_index = 0;
4218 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4219 {
4220 update_ghost_values_finish(vec[comp], exchanger, component_index);
4221 component_index += n_components(vec[comp]);
4222 }
4223 }
4224
4225
4226
4227 // for vector of pointers to vectors
4228 template <int dim,
4229 typename VectorStruct,
4230 typename Number,
4231 typename VectorizedArrayType>
4232 inline void
4233 update_ghost_values_finish(
4234 const std::vector<VectorStruct *> & vec,
4235 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4236 {
4237 unsigned int component_index = 0;
4238 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4239 {
4240 update_ghost_values_finish(*vec[comp], exchanger, component_index);
4241 component_index += n_components(*vec[comp]);
4242 }
4243 }
4244
4245
4246
4247 //
4248 // compress_start
4249 //
4250
4251 // for block vectors
4252 template <int dim,
4253 typename VectorStruct,
4254 typename Number,
4255 typename VectorizedArrayType,
4256 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4257 * = nullptr>
4258 inline void
4259 compress_start(
4260 VectorStruct & vec,
4261 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4262 const unsigned int channel = 0)
4263 {
4264 if (get_communication_block_size(vec) < vec.n_blocks())
4265 vec.compress(VectorOperation::add);
4266 else
4267 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4268 compress_start(vec.block(i), exchanger, channel + i);
4269 }
4270
4271
4272
4273 // for non-block vectors
4274 template <int dim,
4275 typename VectorStruct,
4276 typename Number,
4277 typename VectorizedArrayType,
4278 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4279 * = nullptr>
4280 inline void
4281 compress_start(
4282 VectorStruct & vec,
4283 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4284 const unsigned int channel = 0)
4285 {
4286 exchanger.compress_start(channel, vec);
4287 }
4288
4289
4290
4291 // for std::vector of vectors
4292 template <int dim,
4293 typename VectorStruct,
4294 typename Number,
4295 typename VectorizedArrayType>
4296 inline void
4297 compress_start(
4298 std::vector<VectorStruct> & vec,
4299 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4300 {
4301 unsigned int component_index = 0;
4302 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4303 {
4304 compress_start(vec[comp], exchanger, component_index);
4305 component_index += n_components(vec[comp]);
4306 }
4307 }
4308
4309
4310
4311 // for std::vector of pointer to vectors
4312 template <int dim,
4313 typename VectorStruct,
4314 typename Number,
4315 typename VectorizedArrayType>
4316 inline void
4317 compress_start(
4318 std::vector<VectorStruct *> & vec,
4319 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4320 {
4321 unsigned int component_index = 0;
4322 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4323 {
4324 compress_start(*vec[comp], exchanger, component_index);
4325 component_index += n_components(*vec[comp]);
4326 }
4327 }
4328
4329
4330
4331 //
4332 // compress_finish
4333 //
4334
4335 // for block vectors
4336 template <int dim,
4337 typename VectorStruct,
4338 typename Number,
4339 typename VectorizedArrayType,
4340 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4341 * = nullptr>
4342 inline void
4343 compress_finish(
4344 VectorStruct & vec,
4345 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4346 const unsigned int channel = 0)
4347 {
4348 if (get_communication_block_size(vec) < vec.n_blocks())
4349 {
4350 // do nothing, everything has already been completed in the _start()
4351 // call
4352 }
4353 else
4354 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4355 compress_finish(vec.block(i), exchanger, channel + i);
4356 }
4357
4358
4359
4360 // for non-block vectors
4361 template <int dim,
4362 typename VectorStruct,
4363 typename Number,
4364 typename VectorizedArrayType,
4365 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4366 * = nullptr>
4367 inline void
4368 compress_finish(
4369 VectorStruct & vec,
4370 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4371 const unsigned int channel = 0)
4372 {
4373 exchanger.compress_finish(channel, vec);
4374 }
4375
4376
4377
4378 // for std::vector of vectors
4379 template <int dim,
4380 typename VectorStruct,
4381 typename Number,
4382 typename VectorizedArrayType>
4383 inline void
4384 compress_finish(
4385 std::vector<VectorStruct> & vec,
4386 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4387 {
4388 unsigned int component_index = 0;
4389 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4390 {
4391 compress_finish(vec[comp], exchanger, component_index);
4392 component_index += n_components(vec[comp]);
4393 }
4394 }
4395
4396
4397
4398 // for std::vector of pointer to vectors
4399 template <int dim,
4400 typename VectorStruct,
4401 typename Number,
4402 typename VectorizedArrayType>
4403 inline void
4404 compress_finish(
4405 std::vector<VectorStruct *> & vec,
4406 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4407 {
4408 unsigned int component_index = 0;
4409 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4410 {
4411 compress_finish(*vec[comp], exchanger, component_index);
4412 component_index += n_components(*vec[comp]);
4413 }
4414 }
4415
4416
4417
4418 //
4419 // reset_ghost_values:
4420 //
4421 // if the input vector did not have ghosts imported, clear them here again
4422 // in order to avoid subsequent operations e.g. in linear solvers to work
4423 // with ghosts all the time
4424 //
4425
4426 // for block vectors
4427 template <int dim,
4428 typename VectorStruct,
4429 typename Number,
4430 typename VectorizedArrayType,
4431 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4432 * = nullptr>
4433 inline void
4434 reset_ghost_values(
4435 const VectorStruct & vec,
4436 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4437 {
4438 // return immediately if there is nothing to do.
4439 if (exchanger.ghosts_were_set == true)
4440 return;
4441
4442 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4443 reset_ghost_values(vec.block(i), exchanger);
4444 }
4445
4446
4447
4448 // for non-block vectors
4449 template <int dim,
4450 typename VectorStruct,
4451 typename Number,
4452 typename VectorizedArrayType,
4453 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4454 * = nullptr>
4455 inline void
4456 reset_ghost_values(
4457 const VectorStruct & vec,
4458 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4459 {
4460 exchanger.reset_ghost_values(vec);
4461 }
4462
4463
4464
4465 // for std::vector of vectors
4466 template <int dim,
4467 typename VectorStruct,
4468 typename Number,
4469 typename VectorizedArrayType>
4470 inline void
4471 reset_ghost_values(
4472 const std::vector<VectorStruct> & vec,
4473 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4474 {
4475 // return immediately if there is nothing to do.
4476 if (exchanger.ghosts_were_set == true)
4477 return;
4478
4479 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4480 reset_ghost_values(vec[comp], exchanger);
4481 }
4482
4483
4484
4485 // for std::vector of pointer to vectors
4486 template <int dim,
4487 typename VectorStruct,
4488 typename Number,
4489 typename VectorizedArrayType>
4490 inline void
4491 reset_ghost_values(
4492 const std::vector<VectorStruct *> & vec,
4493 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4494 {
4495 // return immediately if there is nothing to do.
4496 if (exchanger.ghosts_were_set == true)
4497 return;
4498
4499 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4500 reset_ghost_values(*vec[comp], exchanger);
4501 }
4502
4503
4504
4505 //
4506 // zero_vector_region
4507 //
4508
4509 // for block vectors
4510 template <int dim,
4511 typename VectorStruct,
4512 typename Number,
4513 typename VectorizedArrayType,
4514 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4515 * = nullptr>
4516 inline void
4517 zero_vector_region(
4518 const unsigned int range_index,
4519 VectorStruct & vec,
4520 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4521 {
4522 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4523 exchanger.zero_vector_region(range_index, vec.block(i));
4524 }
4525
4526
4527
4528 // for non-block vectors
4529 template <int dim,
4530 typename VectorStruct,
4531 typename Number,
4532 typename VectorizedArrayType,
4533 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4534 * = nullptr>
4535 inline void
4536 zero_vector_region(
4537 const unsigned int range_index,
4538 VectorStruct & vec,
4539 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4540 {
4541 exchanger.zero_vector_region(range_index, vec);
4542 }
4543
4544
4545
4546 // for std::vector of vectors
4547 template <int dim,
4548 typename VectorStruct,
4549 typename Number,
4550 typename VectorizedArrayType>
4551 inline void
4552 zero_vector_region(
4553 const unsigned int range_index,
4554 std::vector<VectorStruct> & vec,
4555 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4556 {
4557 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4558 zero_vector_region(range_index, vec[comp], exchanger);
4559 }
4560
4561
4562
4563 // for std::vector of pointers to vectors
4564 template <int dim,
4565 typename VectorStruct,
4566 typename Number,
4567 typename VectorizedArrayType>
4568 inline void
4569 zero_vector_region(
4570 const unsigned int range_index,
4571 std::vector<VectorStruct *> & vec,
4572 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4573 {
4574 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4575 zero_vector_region(range_index, *vec[comp], exchanger);
4576 }
4577
4578
4579
4580 // Apply a unit matrix operation to constrained DoFs: Default cases where we
4581 // cannot detect a LinearAlgebra::distributed::Vector, we do not do
4582 // anything, else we apply the constraints as a unit operation
4583 template <typename VectorStruct1, typename VectorStruct2>
4584 inline void
4585 apply_operation_to_constrained_dofs(const std::vector<unsigned int> &,
4586 const VectorStruct1 &,
4587 VectorStruct2 &)
4588 {}
4589
4590 template <typename Number>
4591 inline void
4592 apply_operation_to_constrained_dofs(
4593 const std::vector<unsigned int> & constrained_dofs,
4594 const LinearAlgebra::distributed::Vector<Number> &src,
4595 LinearAlgebra::distributed::Vector<Number> & dst)
4596 {
4597 for (const unsigned int i : constrained_dofs)
4598 dst.local_element(i) = src.local_element(i);
4599 }
4600
4601
4602 namespace MatrixFreeFunctions
4603 {
4604 // struct to select between a const interface and a non-const interface
4605 // for MFWorker
4606 template <typename, typename, typename, typename, bool>
4607 struct InterfaceSelector
4608 {};
4609
4610 // Version for constant functions
4611 template <typename MF,
4612 typename InVector,
4613 typename OutVector,
4614 typename Container>
4615 struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4616 {
4617 using function_type = void (Container::*)(
4618 const MF &,
4619 OutVector &,
4620 const InVector &,
4621 const std::pair<unsigned int, unsigned int> &) const;
4622 };
4623
4624 // Version for non-constant functions
4625 template <typename MF,
4626 typename InVector,
4627 typename OutVector,
4628 typename Container>
4629 struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4630 {
4631 using function_type =
4632 void (Container::*)(const MF &,
4633 OutVector &,
4634 const InVector &,
4635 const std::pair<unsigned int, unsigned int> &);
4636 };
4637 } // namespace MatrixFreeFunctions
4638
4639
4640
4641 // A implementation class for the worker object that runs the various
4642 // operations we want to perform during the matrix-free loop
4643 template <typename MF,
4644 typename InVector,
4645 typename OutVector,
4646 typename Container,
4647 bool is_constant>
4648 class MFWorker : public MFWorkerInterface
4649 {
4650 public:
4651 // An alias to make the arguments further down more readable
4652 using function_type = typename MatrixFreeFunctions::
4653 InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4654 function_type;
4655
4656 // constructor, binds all the arguments to this class
4657 MFWorker(const MF & matrix_free,
4658 const InVector & src,
4659 OutVector & dst,
4660 const bool zero_dst_vector_setting,
4661 const Container & container,
4662 function_type cell_function,
4663 function_type face_function,
4664 function_type boundary_function,
4665 const typename MF::DataAccessOnFaces src_vector_face_access =
4666 MF::DataAccessOnFaces::none,
4667 const typename MF::DataAccessOnFaces dst_vector_face_access =
4668 MF::DataAccessOnFaces::none,
4669 const std::function<void(const unsigned int, const unsigned int)>
4670 &operation_before_loop = {},
4671 const std::function<void(const unsigned int, const unsigned int)>
4672 & operation_after_loop = {},
4673 const unsigned int dof_handler_index_pre_post = 0)
4674 : matrix_free(matrix_free)
4675 , container(const_cast<Container &>(container))
4676 , cell_function(cell_function)
4677 , face_function(face_function)
4678 , boundary_function(boundary_function)
4679 , src(src)
4680 , dst(dst)
4681 , src_data_exchanger(matrix_free,
4682 src_vector_face_access,
4683 n_components(src))
4684 , dst_data_exchanger(matrix_free,
4685 dst_vector_face_access,
4686 n_components(dst))
4687 , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4688 , zero_dst_vector_setting(zero_dst_vector_setting &&
4689 !src_and_dst_are_same)
4690 , operation_before_loop(operation_before_loop)
4691 , operation_after_loop(operation_after_loop)
4692 , dof_handler_index_pre_post(dof_handler_index_pre_post)
4693 {
4694 Assert(!has_ghost_elements(dst),
4695 ExcMessage("The destination vector passed to the matrix-free "
4696 "loop is ghosted. This is not allowed."));
4697 }
4698
4699 // Runs the cell work. If no function is given, nothing is done
4700 virtual void
4701 cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4702 {
4703 if (cell_function != nullptr && cell_range.second > cell_range.first)
4704 for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4705 {
4706 const auto cell_subrange =
4707 matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4708
4709 if (cell_subrange.second <= cell_subrange.first)
4710 continue;
4711
4712 (container.*
4713 cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4714 }
4715 }
4716
4717 virtual void
4718 cell(const unsigned int range_index) override
4719 {
4720 process_range(cell_function,
4721 matrix_free.get_task_info().cell_partition_data_hp_ptr,
4722 matrix_free.get_task_info().cell_partition_data_hp,
4723 range_index);
4724 }
4725
4726 virtual void
4727 face(const unsigned int range_index) override
4728 {
4729 process_range(face_function,
4730 matrix_free.get_task_info().face_partition_data_hp_ptr,
4731 matrix_free.get_task_info().face_partition_data_hp,
4732 range_index);
4733 }
4734
4735 virtual void
4736 boundary(const unsigned int range_index) override
4737 {
4738 process_range(boundary_function,
4739 matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4740 matrix_free.get_task_info().boundary_partition_data_hp,
4741 range_index);
4742 }
4743
4744 private:
4745 void
4746 process_range(const function_type & fu,
4747 const std::vector<unsigned int> &ptr,
4748 const std::vector<unsigned int> &data,
4749 const unsigned int range_index)
4750 {
4751 if (fu == nullptr)
4752 return;
4753
4754 AssertIndexRange(range_index + 1, ptr.size());
4755 for (unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4756 {
4757 AssertIndexRange(2 * i + 1, data.size());
4758 (container.*fu)(matrix_free,
4759 this->dst,
4760 this->src,
4761 std::make_pair(data[2 * i], data[2 * i + 1]));
4762 }
4763 }
4764
4765 public:
4766 // Starts the communication for the update ghost values operation. We
4767 // cannot call this update if ghost and destination are the same because
4768 // that would introduce spurious entries in the destination (there is also
4769 // the problem that reading from a vector that we also write to is usually
4770 // not intended in case there is overlap, but this is up to the
4771 // application code to decide and we cannot catch this case here).
4772 virtual void
4773 vector_update_ghosts_start() override
4774 {
4775 if (!src_and_dst_are_same)
4776 internal::update_ghost_values_start(src, src_data_exchanger);
4777 }
4778
4779 // Finishes the communication for the update ghost values operation
4780 virtual void
4781 vector_update_ghosts_finish() override
4782 {
4783 if (!src_and_dst_are_same)
4784 internal::update_ghost_values_finish(src, src_data_exchanger);
4785 }
4786
4787 // Starts the communication for the vector compress operation
4788 virtual void
4789 vector_compress_start() override
4790 {
4791 internal::compress_start(dst, dst_data_exchanger);
4792 }
4793
4794 // Finishes the communication for the vector compress operation
4795 virtual void
4796 vector_compress_finish() override
4797 {
4798 internal::compress_finish(dst, dst_data_exchanger);
4799 if (!src_and_dst_are_same)
4800 internal::reset_ghost_values(src, src_data_exchanger);
4801 }
4802
4803 // Zeros the given input vector
4804 virtual void
4805 zero_dst_vector_range(const unsigned int range_index) override
4806 {
4807 if (zero_dst_vector_setting)
4808 internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4809 }
4810
4811 virtual void
4812 cell_loop_pre_range(const unsigned int range_index) override
4813 {
4814 if (operation_before_loop)
4815 {
4816 const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4817 matrix_free.get_dof_info(dof_handler_index_pre_post);
4818 if (range_index == numbers::invalid_unsigned_int)
4819 {
4820 // Case with threaded loop -> currently no overlap implemented
4821 ::parallel::apply_to_subranges(
4822 0U,
4823 dof_info.vector_partitioner->locally_owned_size(),
4824 operation_before_loop,
4825 internal::VectorImplementation::minimum_parallel_grain_size);
4826 }
4827 else
4828 {
4829 AssertIndexRange(range_index,
4830 dof_info.cell_loop_pre_list_index.size() - 1);
4831 for (unsigned int id =
4832 dof_info.cell_loop_pre_list_index[range_index];
4833 id != dof_info.cell_loop_pre_list_index[range_index + 1];
4834 ++id)
4835 operation_before_loop(dof_info.cell_loop_pre_list[id].first,
4836 dof_info.cell_loop_pre_list[id].second);
4837 }
4838 }
4839 }
4840
4841 virtual void
4842 cell_loop_post_range(const unsigned int range_index) override
4843 {
4844 if (operation_after_loop)
4845 {
4846 // Run unit matrix operation on constrained dofs if we are at the
4847 // last range
4848 const std::vector<unsigned int> &partition_row_index =
4849 matrix_free.get_task_info().partition_row_index;
4850 if (range_index ==
4851 partition_row_index[partition_row_index.size() - 2] - 1)
4852 apply_operation_to_constrained_dofs(
4853 matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
4854 src,
4855 dst);
4856
4857 const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4858 matrix_free.get_dof_info(dof_handler_index_pre_post);
4859 if (range_index == numbers::invalid_unsigned_int)
4860 {
4861 // Case with threaded loop -> currently no overlap implemented
4862 ::parallel::apply_to_subranges(
4863 0U,
4864 dof_info.vector_partitioner->locally_owned_size(),
4865 operation_after_loop,
4866 internal::VectorImplementation::minimum_parallel_grain_size);
4867 }
4868 else
4869 {
4870 AssertIndexRange(range_index,
4871 dof_info.cell_loop_post_list_index.size() - 1);
4872 for (unsigned int id =
4873 dof_info.cell_loop_post_list_index[range_index];
4874 id != dof_info.cell_loop_post_list_index[range_index + 1];
4875 ++id)
4876 operation_after_loop(dof_info.cell_loop_post_list[id].first,
4877 dof_info.cell_loop_post_list[id].second);
4878 }
4879 }
4880 }
4881
4882 private:
4883 const MF & matrix_free;
4884 Container & container;
4885 function_type cell_function;
4886 function_type face_function;
4887 function_type boundary_function;
4888
4889 const InVector &src;
4890 OutVector & dst;
4891 VectorDataExchange<MF::dimension,
4892 typename MF::value_type,
4893 typename MF::vectorized_value_type>
4894 src_data_exchanger;
4895 VectorDataExchange<MF::dimension,
4896 typename MF::value_type,
4897 typename MF::vectorized_value_type>
4898 dst_data_exchanger;
4899 const bool src_and_dst_are_same;
4900 const bool zero_dst_vector_setting;
4901 const std::function<void(const unsigned int, const unsigned int)>
4902 operation_before_loop;
4903 const std::function<void(const unsigned int, const unsigned int)>
4904 operation_after_loop;
4905 const unsigned int dof_handler_index_pre_post;
4906 };
4907
4908
4909
4914 template <class MF, typename InVector, typename OutVector>
4915 struct MFClassWrapper
4916 {
4917 using function_type =
4918 std::function<void(const MF &,
4919 OutVector &,
4920 const InVector &,
4921 const std::pair<unsigned int, unsigned int> &)>;
4922
4923 MFClassWrapper(const function_type cell,
4924 const function_type face,
4925 const function_type boundary)
4926 : cell(cell)
4927 , face(face)
4928 , boundary(boundary)
4929 {}
4930
4931 void
4932 cell_integrator(const MF & mf,
4933 OutVector & dst,
4934 const InVector & src,
4935 const std::pair<unsigned int, unsigned int> &range) const
4936 {
4937 if (cell)
4938 cell(mf, dst, src, range);
4939 }
4940
4941 void
4942 face_integrator(const MF & mf,
4943 OutVector & dst,
4944 const InVector & src,
4945 const std::pair<unsigned int, unsigned int> &range) const
4946 {
4947 if (face)
4948 face(mf, dst, src, range);
4949 }
4950
4951 void
4952 boundary_integrator(
4953 const MF & mf,
4954 OutVector & dst,
4955 const InVector & src,
4956 const std::pair<unsigned int, unsigned int> &range) const
4957 {
4958 if (boundary)
4959 boundary(mf, dst, src, range);
4960 }
4961
4962 const function_type cell;
4963 const function_type face;
4964 const function_type boundary;
4965 };
4966
4967} // end of namespace internal
4968
4969
4970
4971template <int dim, typename Number, typename VectorizedArrayType>
4972template <typename OutVector, typename InVector>
4973inline void
4974MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
4975 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4976 OutVector &,
4977 const InVector &,
4978 const std::pair<unsigned int, unsigned int> &)>
4979 & cell_operation,
4980 OutVector & dst,
4981 const InVector &src,
4982 const bool zero_dst_vector) const
4983{
4984 using Wrapper =
4985 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4986 InVector,
4987 OutVector>;
4988 Wrapper wrap(cell_operation, nullptr, nullptr);
4989 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4990 InVector,
4991 OutVector,
4992 Wrapper,
4993 true>
4994 worker(*this,
4995 src,
4996 dst,
4997 zero_dst_vector,
4998 wrap,
4999 &Wrapper::cell_integrator,
5000 &Wrapper::face_integrator,
5001 &Wrapper::boundary_integrator);
5002
5003 task_info.loop(worker);
5004}
5005
5006
5007
5008template <int dim, typename Number, typename VectorizedArrayType>
5009template <typename OutVector, typename InVector>
5010inline void
5011MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5012 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5013 OutVector &,
5014 const InVector &,
5015 const std::pair<unsigned int, unsigned int> &)>
5016 & cell_operation,
5017 OutVector & dst,
5018 const InVector &src,
5019 const std::function<void(const unsigned int, const unsigned int)>
5020 &operation_before_loop,
5021 const std::function<void(const unsigned int, const unsigned int)>
5022 & operation_after_loop,
5023 const unsigned int dof_handler_index_pre_post) const
5024{
5025 using Wrapper =
5026 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5027 InVector,
5028 OutVector>;
5029 Wrapper wrap(cell_operation, nullptr, nullptr);
5030 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5031 InVector,
5032 OutVector,
5033 Wrapper,
5034 true>
5035 worker(*this,
5036 src,
5037 dst,
5038 false,
5039 wrap,
5040 &Wrapper::cell_integrator,
5041 &Wrapper::face_integrator,
5042 &Wrapper::boundary_integrator,
5043 DataAccessOnFaces::none,
5044 DataAccessOnFaces::none,
5045 operation_before_loop,
5046 operation_after_loop,
5047 dof_handler_index_pre_post);
5048
5049 task_info.loop(worker);
5050}
5051
5052
5053
5054template <int dim, typename Number, typename VectorizedArrayType>
5055template <typename OutVector, typename InVector>
5056inline void
5057MatrixFree<dim, Number, VectorizedArrayType>::loop(
5058 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5059 OutVector &,
5060 const InVector &,
5061 const std::pair<unsigned int, unsigned int> &)>
5062 &cell_operation,
5063 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5064 OutVector &,
5065 const InVector &,
5066 const std::pair<unsigned int, unsigned int> &)>
5067 &face_operation,
5068 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5069 OutVector &,
5070 const InVector &,
5071 const std::pair<unsigned int, unsigned int> &)>
5072 & boundary_operation,
5073 OutVector & dst,
5074 const InVector & src,
5075 const bool zero_dst_vector,
5076 const DataAccessOnFaces dst_vector_face_access,
5077 const DataAccessOnFaces src_vector_face_access) const
5078{
5079 using Wrapper =
5080 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5081 InVector,
5082 OutVector>;
5083 Wrapper wrap(cell_operation, face_operation, boundary_operation);
5084 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5085 InVector,
5086 OutVector,
5087 Wrapper,
5088 true>
5089 worker(*this,
5090 src,
5091 dst,
5092 zero_dst_vector,
5093 wrap,
5094 &Wrapper::cell_integrator,
5095 &Wrapper::face_integrator,
5096 &Wrapper::boundary_integrator,
5097 src_vector_face_access,
5098 dst_vector_face_access);
5099
5100 task_info.loop(worker);
5101}
5102
5103
5104
5105template <int dim, typename Number, typename VectorizedArrayType>
5106template <typename CLASS, typename OutVector, typename InVector>
5107inline void
5108MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5109 void (CLASS::*function_pointer)(
5110 const MatrixFree<dim, Number, VectorizedArrayType> &,
5111 OutVector &,
5112 const InVector &,
5113 const std::pair<unsigned int, unsigned int> &) const,
5114 const CLASS * owning_class,
5115 OutVector & dst,
5116 const InVector &src,
5117 const bool zero_dst_vector) const
5118{
5119 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5120 InVector,
5121 OutVector,
5122 CLASS,
5123 true>
5124 worker(*this,
5125 src,
5126 dst,
5127 zero_dst_vector,
5128 *owning_class,
5129 function_pointer,
5130 nullptr,
5131 nullptr);
5132 task_info.loop(worker);
5133}
5134
5135
5136
5137template <int dim, typename Number, typename VectorizedArrayType>
5138template <typename CLASS, typename OutVector, typename InVector>
5139inline void
5140MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5141 void (CLASS::*function_pointer)(
5142 const MatrixFree<dim, Number, VectorizedArrayType> &,
5143 OutVector &,
5144 const InVector &,
5145 const std::pair<unsigned int, unsigned int> &) const,
5146 const CLASS * owning_class,
5147 OutVector & dst,
5148 const InVector &src,
5149 const std::function<void(const unsigned int, const unsigned int)>
5150 &operation_before_loop,
5151 const std::function<void(const unsigned int, const unsigned int)>
5152 & operation_after_loop,
5153 const unsigned int dof_handler_index_pre_post) const
5154{
5155 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5156 InVector,
5157 OutVector,
5158 CLASS,
5159 true>
5160 worker(*this,
5161 src,
5162 dst,
5163 false,
5164 *owning_class,
5165 function_pointer,
5166 nullptr,
5167 nullptr,
5168 DataAccessOnFaces::none,
5169 DataAccessOnFaces::none,
5170 operation_before_loop,
5171 operation_after_loop,
5172 dof_handler_index_pre_post);
5173 task_info.loop(worker);
5174}
5175
5176
5177
5178template <int dim, typename Number, typename VectorizedArrayType>
5179template <typename CLASS, typename OutVector, typename InVector>
5180inline void
5181MatrixFree<dim, Number, VectorizedArrayType>::loop(
5182 void (CLASS::*cell_operation)(
5183 const MatrixFree<dim, Number, VectorizedArrayType> &,
5184 OutVector &,
5185 const InVector &,
5186 const std::pair<unsigned int, unsigned int> &) const,
5187 void (CLASS::*face_operation)(
5188 const MatrixFree<dim, Number, VectorizedArrayType> &,
5189 OutVector &,
5190 const InVector &,
5191 const std::pair<unsigned int, unsigned int> &) const,
5192 void (CLASS::*boundary_operation)(
5193 const MatrixFree<dim, Number, VectorizedArrayType> &,
5194 OutVector &,
5195 const InVector &,
5196 const std::pair<unsigned int, unsigned int> &) const,
5197 const CLASS * owning_class,
5198 OutVector & dst,
5199 const InVector & src,
5200 const bool zero_dst_vector,
5201 const DataAccessOnFaces dst_vector_face_access,
5202 const DataAccessOnFaces src_vector_face_access) const
5203{
5204 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5205 InVector,
5206 OutVector,
5207 CLASS,
5208 true>
5209 worker(*this,
5210 src,
5211 dst,
5212 zero_dst_vector,
5213 *owning_class,
5214 cell_operation,
5215 face_operation,
5216 boundary_operation,
5217 src_vector_face_access,
5218 dst_vector_face_access);
5219 task_info.loop(worker);
5220}
5221
5222
5223
5224template <int dim, typename Number, typename VectorizedArrayType>
5225template <typename CLASS, typename OutVector, typename InVector>
5226inline void
5227MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5228 void (CLASS::*function_pointer)(
5229 const MatrixFree<dim, Number, VectorizedArrayType> &,
5230 OutVector &,
5231 const InVector &,
5232 const std::pair<unsigned int, unsigned int> &),
5233 CLASS * owning_class,
5234 OutVector & dst,
5235 const InVector &src,
5236 const bool zero_dst_vector) const
5237{
5238 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5239 InVector,
5240 OutVector,
5241 CLASS,
5242 false>
5243 worker(*this,
5244 src,
5245 dst,
5246 zero_dst_vector,
5247 *owning_class,
5248 function_pointer,
5249 nullptr,
5250 nullptr);
5251 task_info.loop(worker);
5252}
5253
5254
5255
5256template <int dim, typename Number, typename VectorizedArrayType>
5257template <typename CLASS, typename OutVector, typename InVector>
5258inline void
5259MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5260 void (CLASS::*function_pointer)(
5261 const MatrixFree<dim, Number, VectorizedArrayType> &,
5262 OutVector &,
5263 const InVector &,
5264 const std::pair<unsigned int, unsigned int> &),
5265 CLASS * owning_class,
5266 OutVector & dst,
5267 const InVector &src,
5268 const std::function<void(const unsigned int, const unsigned int)>
5269 &operation_before_loop,
5270 const std::function<void(const unsigned int, const unsigned int)>
5271 & operation_after_loop,
5272 const unsigned int dof_handler_index_pre_post) const
5273{
5274 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5275 InVector,
5276 OutVector,
5277 CLASS,
5278 false>
5279 worker(*this,
5280 src,
5281 dst,
5282 false,
5283 *owning_class,
5284 function_pointer,
5285 nullptr,
5286 nullptr,
5287 DataAccessOnFaces::none,
5288 DataAccessOnFaces::none,
5289 operation_before_loop,
5290 operation_after_loop,
5291 dof_handler_index_pre_post);
5292 task_info.loop(worker);
5293}
5294
5295
5296
5297template <int dim, typename Number, typename VectorizedArrayType>
5298template <typename CLASS, typename OutVector, typename InVector>
5299inline void
5300MatrixFree<dim, Number, VectorizedArrayType>::loop(
5301 void (CLASS::*cell_operation)(
5302 const MatrixFree<dim, Number, VectorizedArrayType> &,
5303 OutVector &,
5304 const InVector &,
5305 const std::pair<unsigned int, unsigned int> &),
5306 void (CLASS::*face_operation)(
5307 const MatrixFree<dim, Number, VectorizedArrayType> &,
5308 OutVector &,
5309 const InVector &,
5310 const std::pair<unsigned int, unsigned int> &),
5311 void (CLASS::*boundary_operation)(
5312 const MatrixFree<dim, Number, VectorizedArrayType> &,
5313 OutVector &,
5314 const InVector &,
5315 const std::pair<unsigned int, unsigned int> &),
5316 CLASS * owning_class,
5317 OutVector & dst,
5318 const InVector & src,
5319 const bool zero_dst_vector,
5320 const DataAccessOnFaces dst_vector_face_access,
5321 const DataAccessOnFaces src_vector_face_access) const
5322{
5323 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5324 InVector,
5325 OutVector,
5326 CLASS,
5327 false>
5328 worker(*this,
5329 src,
5330 dst,
5331 zero_dst_vector,
5332 *owning_class,
5333 cell_operation,
5334 face_operation,
5335 boundary_operation,
5336 src_vector_face_access,
5337 dst_vector_face_access);
5338 task_info.loop(worker);
5339}
5340
5341
5342
5343template <int dim, typename Number, typename VectorizedArrayType>
5344template <typename OutVector, typename InVector>
5345inline void
5346MatrixFree<dim, Number, VectorizedArrayType>::loop(
5347 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5348 OutVector &,
5349 const InVector &,
5350 const std::pair<unsigned int, unsigned int> &)>
5351 &cell_operation,
5352 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5353 OutVector &,
5354 const InVector &,
5355 const std::pair<unsigned int, unsigned int> &)>
5356 &face_operation,
5357 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5358 OutVector &,
5359 const InVector &,
5360 const std::pair<unsigned int, unsigned int> &)>
5361 & boundary_operation,
5362 OutVector & dst,
5363 const InVector &src,
5364 const std::function<void(const unsigned int, const unsigned int)>
5365 &operation_before_loop,
5366 const std::function<void(const unsigned int, const unsigned int)>
5367 & operation_after_loop,
5368 const unsigned int dof_handler_index_pre_post,
5369 const DataAccessOnFaces dst_vector_face_access,
5370 const DataAccessOnFaces src_vector_face_access) const
5371{
5372 using Wrapper =
5373 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5374 InVector,
5375 OutVector>;
5376 Wrapper wrap(cell_operation, face_operation, boundary_operation);
5377 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5378 InVector,
5379 OutVector,
5380 Wrapper,
5381 true>
5382 worker(*this,
5383 src,
5384 dst,
5385 false,
5386 wrap,
5387 &Wrapper::cell_integrator,
5388 &Wrapper::face_integrator,
5389 &Wrapper::boundary_integrator,
5390 src_vector_face_access,
5391 dst_vector_face_access,
5392 operation_before_loop,
5393 operation_after_loop,
5394 dof_handler_index_pre_post);
5395
5396 task_info.loop(worker);
5397}
5398
5399
5400
5401template <int dim, typename Number, typename VectorizedArrayType>
5402template <typename CLASS, typename OutVector, typename InVector>
5403inline void
5404MatrixFree<dim, Number, VectorizedArrayType>::loop(
5405 void (CLASS::*cell_operation)(const MatrixFree &,
5406 OutVector &,
5407 const InVector &,
5408 const std::pair<unsigned int, unsigned int> &)
5409 const,
5410 void (CLASS::*face_operation)(const MatrixFree &,
5411 OutVector &,
5412 const InVector &,
5413 const std::pair<unsigned int, unsigned int> &)
5414 const,
5415 void (CLASS::*boundary_operation)(
5416 const MatrixFree &,
5417 OutVector &,
5418 const InVector &,
5419 const std::pair<unsigned int, unsigned int> &) const,
5420 const CLASS * owning_class,
5421 OutVector & dst,
5422 const InVector &src,
5423 const std::function<void(const unsigned int, const unsigned int)>
5424 &operation_before_loop,
5425 const std::function<void(const unsigned int, const unsigned int)>
5426 & operation_after_loop,
5427 const unsigned int dof_handler_index_pre_post,
5428 const DataAccessOnFaces dst_vector_face_access,
5429 const DataAccessOnFaces src_vector_face_access) const
5430{
5431 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5432 InVector,
5433 OutVector,
5434 CLASS,
5435 true>
5436 worker(*this,
5437 src,
5438 dst,
5439 false,
5440 *owning_class,
5441 cell_operation,
5442 face_operation,
5443 boundary_operation,
5444 src_vector_face_access,
5445 dst_vector_face_access,
5446 operation_before_loop,
5447 operation_after_loop,
5448 dof_handler_index_pre_post);
5449 task_info.loop(worker);
5450}
5451
5452
5453
5454template <int dim, typename Number, typename VectorizedArrayType>
5455template <typename CLASS, typename OutVector, typename InVector>
5456inline void
5457MatrixFree<dim, Number, VectorizedArrayType>::loop(
5458 void (CLASS::*cell_operation)(const MatrixFree &,
5459 OutVector &,
5460 const InVector &,
5461 const std::pair<unsigned int, unsigned int> &),
5462 void (CLASS::*face_operation)(const MatrixFree &,
5463 OutVector &,
5464 const InVector &,
5465 const std::pair<unsigned int, unsigned int> &),
5466 void (CLASS::*boundary_operation)(
5467 const MatrixFree &,
5468 OutVector &,
5469 const InVector &,
5470 const std::pair<unsigned int, unsigned int> &),
5471 const CLASS * owning_class,
5472 OutVector & dst,
5473 const InVector &src,
5474 const std::function<void(const unsigned int, const unsigned int)>
5475 &operation_before_loop,
5476 const std::function<void(const unsigned int, const unsigned int)>
5477 & operation_after_loop,
5478 const unsigned int dof_handler_index_pre_post,
5479 const DataAccessOnFaces dst_vector_face_access,
5480 const DataAccessOnFaces src_vector_face_access) const
5481{
5482 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5483 InVector,
5484 OutVector,
5485 CLASS,
5486 false>
5487 worker(*this,
5488 src,
5489 dst,
5490 false,
5491 *owning_class,
5492 cell_operation,
5493 face_operation,
5494 boundary_operation,
5495 src_vector_face_access,
5496 dst_vector_face_access,
5497 operation_before_loop,
5498 operation_after_loop,
5499 dof_handler_index_pre_post);
5500 task_info.loop(worker);
5501}
5502
5503
5504
5505template <int dim, typename Number, typename VectorizedArrayType>
5506template <typename CLASS, typename OutVector, typename InVector>
5507inline void
5508MatrixFree<dim, Number, VectorizedArrayType>::loop_cell_centric(
5509 void (CLASS::*function_pointer)(
5510 const MatrixFree<dim, Number, VectorizedArrayType> &,
5511 OutVector &,
5512 const InVector &,
5513 const std::pair<unsigned int, unsigned int> &) const,
5514 const CLASS * owning_class,
5515 OutVector & dst,
5516 const InVector & src,
5517 const bool zero_dst_vector,
5518 const DataAccessOnFaces src_vector_face_access) const
5519{
5520 auto src_vector_face_access_temp = src_vector_face_access;
5521 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5522 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5523 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5524 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5525
5526 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5527 InVector,
5528 OutVector,
5529 CLASS,
5530 true>
5531 worker(*this,
5532 src,
5533 dst,
5534 zero_dst_vector,
5535 *owning_class,
5536 function_pointer,
5537 nullptr,
5538 nullptr,
5539 src_vector_face_access_temp,
5540 DataAccessOnFaces::none);
5541 task_info.loop(worker);
5542}
5543
5544
5545
5546template <int dim, typename Number, typename VectorizedArrayType>
5547template <typename CLASS, typename OutVector, typename InVector>
5548inline void
5549MatrixFree<dim, Number, VectorizedArrayType>::loop_cell_centric(
5550 void (CLASS::*function_pointer)(
5551 const MatrixFree<dim, Number, VectorizedArrayType> &,
5552 OutVector &,
5553 const InVector &,
5554 const std::pair<unsigned int, unsigned int> &),
5555 CLASS * owning_class,
5556 OutVector & dst,
5557 const InVector & src,
5558 const bool zero_dst_vector,
5559 const DataAccessOnFaces src_vector_face_access) const
5560{
5561 auto src_vector_face_access_temp = src_vector_face_access;
5562 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5563 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5564 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5565 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5566
5567 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5568 InVector,
5569 OutVector,
5570 CLASS,
5571 false>
5572 worker(*this,
5573 src,
5574 dst,
5575 zero_dst_vector,
5576 *owning_class,
5577 function_pointer,
5578 nullptr,
5579 nullptr,
5580 src_vector_face_access_temp,
5581 DataAccessOnFaces::none);
5582 task_info.loop(worker);
5583}
5584
5585
5586
5587template <int dim, typename Number, typename VectorizedArrayType>
5588template <typename OutVector, typename InVector>
5589inline void
5590MatrixFree<dim, Number, VectorizedArrayType>::loop_cell_centric(
5591 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5592 OutVector &,
5593 const InVector &,
5594 const std::pair<unsigned int, unsigned int> &)>
5595 & cell_operation,
5596 OutVector & dst,
5597 const InVector & src,
5598 const bool zero_dst_vector,
5599 const DataAccessOnFaces src_vector_face_access) const
5600{
5601 auto src_vector_face_access_temp = src_vector_face_access;
5602 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5603 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5604 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5605 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5606
5607 using Wrapper =
5608 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5609 InVector,
5610 OutVector>;
5611 Wrapper wrap(cell_operation, nullptr, nullptr);
5612
5613 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5614 InVector,
5615 OutVector,
5616 Wrapper,
5617 true>
5618 worker(*this,
5619 src,
5620 dst,
5621 zero_dst_vector,
5622 wrap,
5623 &Wrapper::cell_integrator,
5624 &Wrapper::face_integrator,
5625 &Wrapper::boundary_integrator,
5626 src_vector_face_access_temp,
5627 DataAccessOnFaces::none);
5628 task_info.loop(worker);
5629}
5630
5631
5632#endif // ifndef DOXYGEN
5633
5634
5635
5636DEAL_II_NAMESPACE_CLOSE
5637
5638#endif
affine_constraints.h
aligned_vector.h
block_vector_base.h
AlignedVector::begin
iterator begin()
AlignedVector::resize
void resize(const size_type new_size)
LinearAlgebra::distributed::Vector::reinit
void reinit(const size_type size, const bool omit_zeroing_entries=false)
MPI_Comm
Mapping
Abstract base class for mapping classes.
Definition mapping.h:317
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:2285
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::n_ghost_cell_batches
unsigned int n_ghost_cell_batches() const
MatrixFree::loop
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
MatrixFree::get_shape_info
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
MatrixFree::initialize_dof_handlers
void initialize_dof_handlers(const std::vector< const DoFHandler< dim, dim > * > &dof_handlers, const AdditionalData &additional_data)
MatrixFree::get_boundary_id
types::boundary_id get_boundary_id(const unsigned int face_batch_index) const
MatrixFree::shape_info
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > > shape_info
Definition matrix_free.h:2291
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::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::print_memory_consumption
void print_memory_consumption(StreamType &out) const
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_affine_constraints
const AffineConstraints< Number > & get_affine_constraints(const unsigned int dof_handler_index=0) const
MatrixFree::copy_from
void copy_from(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free_base)
MatrixFree::initialize_face_data_vector
void initialize_face_data_vector(AlignedVector< T > &vec) const
MatrixFree::get_dofs_per_cell
unsigned int get_dofs_per_cell(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
MatrixFree::n_constraint_pool_entries
unsigned int n_constraint_pool_entries() const
MatrixFree::cell_loop
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
MatrixFree::constraint_pool_begin
const Number * constraint_pool_begin(const unsigned int pool_index) const
MatrixFree::acquire_scratch_data_non_threadsafe
AlignedVector< Number > * acquire_scratch_data_non_threadsafe() const
MatrixFree::get_locally_owned_set
const IndexSet & get_locally_owned_set(const unsigned int dof_handler_index=0) const
MatrixFree::get_dofs_per_face
unsigned int get_dofs_per_face(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
MatrixFree::task_info
internal::MatrixFreeFunctions::TaskInfo task_info
Definition matrix_free.h:2320
MatrixFree::get_dof_handler
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
MatrixFree::release_scratch_data
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
MatrixFree::loop
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*boundary_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
MatrixFree::mapping_is_initialized
bool mapping_is_initialized
Definition matrix_free.h:2337
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:2360
MatrixFree::get_face_range_category
std::pair< unsigned int, unsigned int > get_face_range_category(const std::pair< unsigned int, unsigned int > face_batch_range) const
MatrixFree::cell_loop
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
MatrixFree::cell_level_index
std::vector< std::pair< unsigned int, unsigned int > > cell_level_index
Definition matrix_free.h:2299
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::constraint_pool_data
std::vector< Number > constraint_pool_data
Definition matrix_free.h:2272
MatrixFree::vectorized_value_type
VectorizedArrayType vectorized_value_type
Definition matrix_free.h:124
MatrixFree::n_cell_batches
unsigned int n_cell_batches() const
MatrixFree::get_ghost_set
const IndexSet & get_ghost_set(const unsigned int dof_handler_index=0) const
MatrixFree::get_cell_level_and_index
std::pair< int, int > get_cell_level_and_index(const unsigned int cell_batch_index, const unsigned int lane_index) const
MatrixFree::indices_initialized
bool indices_initialized() const
MatrixFree::get_face_quadrature
const Quadrature< dim - 1 > & get_face_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
MatrixFree::get_cell_category
unsigned int get_cell_category(const unsigned int cell_batch_index) const
MatrixFree::scratch_pad
Threads::ThreadLocalStorage< std::list< std::pair< bool, AlignedVector< VectorizedArrayType > > > > scratch_pad
Definition matrix_free.h:2348
MatrixFree::get_matrix_free_cell_index
unsigned int get_matrix_free_cell_index(const typename Triangulation< dim >::cell_iterator &cell) const
MatrixFree::value_type
Number value_type
Definition matrix_free.h:123
MatrixFree::memory_consumption
std::size_t memory_consumption() const
MatrixFree::get_face_category
std::pair< unsigned int, unsigned int > get_face_category(const unsigned int face_batch_index) const
MatrixFree::n_boundary_face_batches
unsigned int n_boundary_face_batches() const
MatrixFree::cell_loop
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
MatrixFree::scratch_pad_non_threadsafe
std::list< std::pair< bool, AlignedVector< Number > > > scratch_pad_non_threadsafe
Definition matrix_free.h:2355
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:2250
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::affine_constraints
std::vector< SmartPointer< const AffineConstraints< Number > > > affine_constraints
Definition matrix_free.h:2258
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::mf_cell_indices
std::vector< unsigned int > mf_cell_indices
Definition matrix_free.h:2305
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::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 > &), 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 DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) 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 bool zero_dst_vector=false) const
MatrixFree::dof_info
std::vector< internal::MatrixFreeFunctions::DoFInfo > dof_info
Definition matrix_free.h:2264
MatrixFree::initialize_cell_data_vector
void initialize_cell_data_vector(AlignedVector< T > &vec) const
MatrixFree::n_ghost_inner_face_batches
unsigned int n_ghost_inner_face_batches() const
MatrixFree::release_scratch_data_non_threadsafe
void release_scratch_data_non_threadsafe(const AlignedVector< Number > *memory) const
MatrixFree::renumber_dofs
void renumber_dofs(std::vector< types::global_dof_index > &renumbering, const unsigned int dof_handler_index=0)
MatrixFree::get_face_active_fe_index
unsigned int get_face_active_fe_index(const std::pair< unsigned int, unsigned int > range, const bool is_interior_face=true) const
MatrixFree::n_active_entries_per_face_batch
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
MatrixFree::indices_are_initialized
bool indices_are_initialized
Definition matrix_free.h:2332
MatrixFree::face_info
internal::MatrixFreeFunctions::FaceInfo< VectorizedArrayType::size()> face_info
Definition matrix_free.h:2327
MatrixFree::reinit
void reinit(const MappingType &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< number2 > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
MatrixFree::get_cell_range_category
unsigned int get_cell_range_category(const std::pair< unsigned int, unsigned int > cell_batch_range) const
MatrixFree::loop
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*boundary_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
MatrixFree::reinit
void reinit(const MappingType &mapping, const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< QuadratureType > &quad, const AdditionalData &additional_data=AdditionalData())
MatrixFree::get_faces_by_cells_boundary_id
std::array< types::boundary_id, VectorizedArrayType::size()> get_faces_by_cells_boundary_id(const unsigned int cell_batch_index, const unsigned int face_number) const
MatrixFree::constraint_pool_row_index
std::vector< unsigned int > constraint_pool_row_index
Definition matrix_free.h:2278
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 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 DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
MatrixFree::get_constrained_dofs
const std::vector< unsigned int > & get_constrained_dofs(const unsigned int dof_handler_index=0) const
MatrixFree::cell_level_index_end_local
unsigned int cell_level_index_end_local
Definition matrix_free.h:2313
MatrixFree::n_base_elements
unsigned int n_base_elements(const unsigned int dof_handler_index) 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 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 DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
MatrixFree::is_supported
static bool is_supported(const FiniteElement< dim, spacedim > &fe)
MatrixFree::dimension
static constexpr unsigned int dimension
Definition matrix_free.h:129
Threads::ThreadLocalStorage
A class that provides a separate storage location on each thread that accesses the object.
Definition thread_local_storage.h:105
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
dof_handler.h
second
Point< 2 > second
Definition grid_out.cc:4616
first
Point< 2 > first
Definition grid_out.cc:4615
level
unsigned int level
Definition grid_out.cc:4618
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exceptions.h:1614
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:1787
AssertThrowMPI
#define AssertThrowMPI(error_code)
Definition exceptions.h:1914
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exceptions.h:1855
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcIndexRange
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
StandardExceptions::ExcNotInitialized
static ::ExceptionBase & ExcNotInitialized()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition exceptions.h:1703
DoFHandler::cell_iterator
typename ActiveSelector::cell_iterator cell_iterator
Definition dof_handler.h:469
UpdateFlags
UpdateFlags
Definition fe_update_flags.h:65
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition fe_update_flags.h:135
update_gradients
@ update_gradients
Shape function gradients.
Definition fe_update_flags.h:82
update_default
@ update_default
No update.
Definition fe_update_flags.h:67
la_parallel_vector.h
mapping_collection.h
dof_info.h
mapping_info.h
type_traits.h
TrilinosWrappers::internal::begin
VectorType::value_type * begin(VectorType &V)
Definition trilinos_sparse_matrix.cc:51
Utilities::MPI::job_supports_mpi
bool job_supports_mpi()
Definition mpi.cc:1057
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:261
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition types.h:213
parallel::apply_to_subranges
void apply_to_subranges(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, const Function &f, const unsigned int grainsize)
Definition parallel.h:435
std
STL namespace.
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6024
types::fe_index
unsigned short int fe_index
Definition types.h:60
types::boundary_id
unsigned int boundary_id
Definition types.h:141
q_collection.h
shape_info.h
MatrixFree::AdditionalData::tasks_parallel_scheme
TasksParallelScheme tasks_parallel_scheme
Definition matrix_free.h:347
MatrixFree::AdditionalData::mapping_update_flags_inner_faces
UpdateFlags mapping_update_flags_inner_faces
Definition matrix_free.h:413
MatrixFree::AdditionalData::AdditionalData
AdditionalData(const TasksParallelScheme tasks_parallel_scheme=partition_partition, const unsigned int tasks_block_size=0, const UpdateFlags mapping_update_flags=update_gradients|update_JxW_values, const UpdateFlags mapping_update_flags_boundary_faces=update_default, const UpdateFlags mapping_update_flags_inner_faces=update_default, const UpdateFlags mapping_update_flags_faces_by_cells=update_default, const unsigned int mg_level=numbers::invalid_unsigned_int, const bool store_plain_indices=true, const bool initialize_indices=true, const bool initialize_mapping=true, const bool overlap_communication_computation=true, const bool hold_all_faces_to_owned_cells=false, const bool cell_vectorization_categories_strict=false, const bool allow_ghosted_vectors_in_loops=true)
Definition matrix_free.h:221
MatrixFree::AdditionalData::cell_vectorization_category
std::vector< unsigned int > cell_vectorization_category
Definition matrix_free.h:525
MatrixFree::AdditionalData::operator=
AdditionalData & operator=(const AdditionalData &other)
Definition matrix_free.h:285
MatrixFree::AdditionalData::mapping_update_flags_boundary_faces
UpdateFlags mapping_update_flags_boundary_faces
Definition matrix_free.h:393
MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells
UpdateFlags mapping_update_flags_faces_by_cells
Definition matrix_free.h:441
MatrixFree::AdditionalData::AdditionalData
AdditionalData(const AdditionalData &other)
Definition matrix_free.h:258
internal::MatrixFreeFunctions::DoFInfo::cell_loop_pre_list_index
std::vector< unsigned int > cell_loop_pre_list_index
Definition dof_info.h:745
internal::MatrixFreeFunctions::DoFInfo::cell_loop_post_list_index
std::vector< unsigned int > cell_loop_post_list_index
Definition dof_info.h:758
internal::MatrixFreeFunctions::DoFInfo::vector_zero_range_list
std::vector< std::pair< unsigned int, unsigned int > > vector_zero_range_list
Definition dof_info.h:738
internal::MatrixFreeFunctions::DoFInfo::vector_partitioner
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition dof_info.h:589
internal::MatrixFreeFunctions::DoFInfo::cell_loop_pre_list
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_pre_list
Definition dof_info.h:751
internal::MatrixFreeFunctions::DoFInfo::vector_zero_range_list_index
std::vector< unsigned int > vector_zero_range_list_index
Definition dof_info.h:733
internal::MatrixFreeFunctions::DoFInfo::cell_loop_post_list
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_post_list
Definition dof_info.h:764
internal::MatrixFreeFunctions::TaskInfo::loop
void loop(MFWorkerInterface &worker) const
Definition task_info.cc:351
task_info.h
template_constraints.h
thread_local_storage.h
vector_operation.h
Generated by doxygen 1.9.7