deal.II version GIT relicensing-2289-g1e5549a87a 2024-12-21 21:30:00+00:00
\(\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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2012 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
16#ifndef dealii_matrix_free_h
17#define dealii_matrix_free_h
18
19#include <deal.II/base/config.h>
20
21#include <deal.II/base/aligned_vector.h>
22#include <deal.II/base/exceptions.h>
23#include <deal.II/base/quadrature.h>
24#include <deal.II/base/template_constraints.h>
25#include <deal.II/base/thread_local_storage.h>
26#include <deal.II/base/vectorization.h>
27
28#include <deal.II/dofs/dof_handler.h>
29
30#include <deal.II/fe/fe.h>
31#include <deal.II/fe/mapping.h>
32
33#include <deal.II/hp/mapping_collection.h>
34#include <deal.II/hp/q_collection.h>
35
36#include <deal.II/lac/affine_constraints.h>
37#include <deal.II/lac/block_vector_base.h>
38#include <deal.II/lac/la_parallel_vector.h>
39#include <deal.II/lac/vector_operation.h>
40
41#include <deal.II/matrix_free/dof_info.h>
42#include <deal.II/matrix_free/mapping_info.h>
43#include <deal.II/matrix_free/shape_info.h>
44#include <deal.II/matrix_free/task_info.h>
45#include <deal.II/matrix_free/type_traits.h>
46#include <deal.II/matrix_free/vector_data_exchange.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 EnableObserverPointer
113{
114 static_assert(
115 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
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
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
280
284 AdditionalData &
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(
1093 const std::function<
1094 void(const MatrixFree<dim, Number, VectorizedArrayType> &,
1095 OutVector &,
1096 const InVector &,
1097 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1098 const std::function<void(
1099 const MatrixFree<dim, Number, VectorizedArrayType> &,
1100 OutVector &,
1101 const InVector &,
1102 const std::pair<unsigned int, unsigned int> &)> &inner_face_operation,
1103 const std::function<void(
1104 const MatrixFree<dim, Number, VectorizedArrayType> &,
1105 OutVector &,
1106 const InVector &,
1107 const std::pair<unsigned int, unsigned int> &)> &boundary_face_operation,
1108 OutVector &dst,
1109 const InVector &src,
1110 const bool zero_dst_vector = false,
1111 const DataAccessOnFaces dst_vector_face_access =
1112 DataAccessOnFaces::unspecified,
1113 const DataAccessOnFaces src_vector_face_access =
1114 DataAccessOnFaces::unspecified) const;
1115
1204 template <typename CLASS, typename OutVector, typename InVector>
1205 void
1206 loop(void (CLASS::*cell_operation)(
1207 const MatrixFree &,
1208 OutVector &,
1209 const InVector &,
1210 const std::pair<unsigned int, unsigned int> &) const,
1211 void (CLASS::*inner_face_operation)(
1212 const MatrixFree &,
1213 OutVector &,
1214 const InVector &,
1215 const std::pair<unsigned int, unsigned int> &) const,
1216 void (CLASS::*boundary_face_operation)(
1217 const MatrixFree &,
1218 OutVector &,
1219 const InVector &,
1220 const std::pair<unsigned int, unsigned int> &) const,
1221 const CLASS *owning_class,
1222 OutVector &dst,
1223 const InVector &src,
1224 const bool zero_dst_vector = false,
1225 const DataAccessOnFaces dst_vector_face_access =
1226 DataAccessOnFaces::unspecified,
1227 const DataAccessOnFaces src_vector_face_access =
1228 DataAccessOnFaces::unspecified) const;
1229
1233 template <typename CLASS, typename OutVector, typename InVector>
1234 void
1235 loop(void (CLASS::*cell_operation)(
1236 const MatrixFree &,
1237 OutVector &,
1238 const InVector &,
1239 const std::pair<unsigned int, unsigned int> &),
1240 void (CLASS::*inner_face_operation)(
1241 const MatrixFree &,
1242 OutVector &,
1243 const InVector &,
1244 const std::pair<unsigned int, unsigned int> &),
1245 void (CLASS::*boundary_face_operation)(
1246 const MatrixFree &,
1247 OutVector &,
1248 const InVector &,
1249 const std::pair<unsigned int, unsigned int> &),
1250 CLASS *owning_class,
1251 OutVector &dst,
1252 const InVector &src,
1253 const bool zero_dst_vector = false,
1254 const DataAccessOnFaces dst_vector_face_access =
1255 DataAccessOnFaces::unspecified,
1256 const DataAccessOnFaces src_vector_face_access =
1257 DataAccessOnFaces::unspecified) const;
1258
1380 template <typename CLASS, typename OutVector, typename InVector>
1381 void
1382 loop(void (CLASS::*cell_operation)(
1383 const MatrixFree &,
1384 OutVector &,
1385 const InVector &,
1386 const std::pair<unsigned int, unsigned int> &) const,
1387 void (CLASS::*inner_face_operation)(
1388 const MatrixFree &,
1389 OutVector &,
1390 const InVector &,
1391 const std::pair<unsigned int, unsigned int> &) const,
1392 void (CLASS::*boundary_face_operation)(
1393 const MatrixFree &,
1394 OutVector &,
1395 const InVector &,
1396 const std::pair<unsigned int, unsigned int> &) const,
1397 const CLASS *owning_class,
1398 OutVector &dst,
1399 const InVector &src,
1400 const std::function<void(const unsigned int, const unsigned int)>
1401 &operation_before_loop,
1402 const std::function<void(const unsigned int, const unsigned int)>
1403 &operation_after_loop,
1404 const unsigned int dof_handler_index_pre_post = 0,
1405 const DataAccessOnFaces dst_vector_face_access =
1406 DataAccessOnFaces::unspecified,
1407 const DataAccessOnFaces src_vector_face_access =
1408 DataAccessOnFaces::unspecified) const;
1409
1413 template <typename CLASS, typename OutVector, typename InVector>
1414 void
1415 loop(void (CLASS::*cell_operation)(
1416 const MatrixFree &,
1417 OutVector &,
1418 const InVector &,
1419 const std::pair<unsigned int, unsigned int> &),
1420 void (CLASS::*inner_face_operation)(
1421 const MatrixFree &,
1422 OutVector &,
1423 const InVector &,
1424 const std::pair<unsigned int, unsigned int> &),
1425 void (CLASS::*boundary_face_operation)(
1426 const MatrixFree &,
1427 OutVector &,
1428 const InVector &,
1429 const std::pair<unsigned int, unsigned int> &),
1430 const CLASS *owning_class,
1431 OutVector &dst,
1432 const InVector &src,
1433 const std::function<void(const unsigned int, const unsigned int)>
1434 &operation_before_loop,
1435 const std::function<void(const unsigned int, const unsigned int)>
1436 &operation_after_loop,
1437 const unsigned int dof_handler_index_pre_post = 0,
1438 const DataAccessOnFaces dst_vector_face_access =
1439 DataAccessOnFaces::unspecified,
1440 const DataAccessOnFaces src_vector_face_access =
1441 DataAccessOnFaces::unspecified) const;
1442
1448 template <typename OutVector, typename InVector>
1449 void
1450 loop(
1451 const std::function<
1452 void(const MatrixFree<dim, Number, VectorizedArrayType> &,
1453 OutVector &,
1454 const InVector &,
1455 const std::pair<unsigned int, unsigned int> &)> &cell_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> &)> &inner_face_operation,
1461 const std::function<void(
1462 const MatrixFree<dim, Number, VectorizedArrayType> &,
1463 OutVector &,
1464 const InVector &,
1465 const std::pair<unsigned int, unsigned int> &)> &boundary_face_operation,
1466 OutVector &dst,
1467 const InVector &src,
1468 const std::function<void(const unsigned int, const unsigned int)>
1469 &operation_before_loop,
1470 const std::function<void(const unsigned int, const unsigned int)>
1471 &operation_after_loop,
1472 const unsigned int dof_handler_index_pre_post = 0,
1473 const DataAccessOnFaces dst_vector_face_access =
1474 DataAccessOnFaces::unspecified,
1475 const DataAccessOnFaces src_vector_face_access =
1476 DataAccessOnFaces::unspecified) const;
1477
1542 template <typename CLASS, typename OutVector, typename InVector>
1543 void
1544 loop_cell_centric(void (CLASS::*cell_operation)(
1545 const MatrixFree &,
1546 OutVector &,
1547 const InVector &,
1548 const std::pair<unsigned int, unsigned int> &) const,
1549 const CLASS *owning_class,
1550 OutVector &dst,
1551 const InVector &src,
1552 const bool zero_dst_vector = false,
1553 const DataAccessOnFaces src_vector_face_access =
1554 DataAccessOnFaces::unspecified) const;
1555
1559 template <typename CLASS, typename OutVector, typename InVector>
1560 void
1561 loop_cell_centric(void (CLASS::*cell_operation)(
1562 const MatrixFree &,
1563 OutVector &,
1564 const InVector &,
1565 const std::pair<unsigned int, unsigned int> &),
1566 CLASS *owning_class,
1567 OutVector &dst,
1568 const InVector &src,
1569 const bool zero_dst_vector = false,
1570 const DataAccessOnFaces src_vector_face_access =
1571 DataAccessOnFaces::unspecified) const;
1572
1576 template <typename OutVector, typename InVector>
1577 void
1578 loop_cell_centric(
1579 const std::function<void(const MatrixFree &,
1580 OutVector &,
1581 const InVector &,
1582 const std::pair<unsigned int, unsigned int> &)>
1583 &cell_operation,
1584 OutVector &dst,
1585 const InVector &src,
1586 const bool zero_dst_vector = false,
1587 const DataAccessOnFaces src_vector_face_access =
1588 DataAccessOnFaces::unspecified) const;
1589
1597 std::pair<unsigned int, unsigned int>
1598 create_cell_subrange_hp(const std::pair<unsigned int, unsigned int> &range,
1599 const unsigned int fe_degree,
1600 const unsigned int dof_handler_index = 0) const;
1601
1608 std::pair<unsigned int, unsigned int>
1609 create_cell_subrange_hp_by_index(
1610 const std::pair<unsigned int, unsigned int> &range,
1611 const unsigned int fe_index,
1612 const unsigned int dof_handler_index = 0) const;
1613
1617 unsigned int
1618 n_active_fe_indices() const;
1619
1623 unsigned int
1624 get_cell_active_fe_index(
1625 const std::pair<unsigned int, unsigned int> range) const;
1626
1630 unsigned int
1631 get_face_active_fe_index(const std::pair<unsigned int, unsigned int> range,
1632 const bool is_interior_face = true) const;
1633
1645 template <typename T>
1646 void
1647 initialize_cell_data_vector(AlignedVector<T> &vec) const;
1648
1654 template <typename T>
1655 void
1656 initialize_face_data_vector(AlignedVector<T> &vec) const;
1657
1671 template <typename VectorType>
1672 void
1673 initialize_dof_vector(VectorType &vec,
1674 const unsigned int dof_handler_index = 0) const;
1675
1693 template <typename Number2>
1694 void
1695 initialize_dof_vector(LinearAlgebra::distributed::Vector<Number2> &vec,
1696 const unsigned int dof_handler_index = 0) const;
1697
1708 const std::shared_ptr<const Utilities::MPI::Partitioner> &
1709 get_vector_partitioner(const unsigned int dof_handler_index = 0) const;
1710
1714 const IndexSet &
1715 get_locally_owned_set(const unsigned int dof_handler_index = 0) const;
1716
1720 const IndexSet &
1721 get_ghost_set(const unsigned int dof_handler_index = 0) const;
1722
1732 const std::vector<unsigned int> &
1733 get_constrained_dofs(const unsigned int dof_handler_index = 0) const;
1734
1745 void
1746 renumber_dofs(std::vector<types::global_dof_index> &renumbering,
1747 const unsigned int dof_handler_index = 0);
1748
1758 template <int spacedim>
1759 static bool
1760 is_supported(const FiniteElement<dim, spacedim> &fe);
1761
1765 unsigned int
1766 n_components() const;
1767
1772 unsigned int
1773 n_base_elements(const unsigned int dof_handler_index) const;
1774
1782 unsigned int
1783 n_physical_cells() const;
1784
1794 unsigned int
1795 n_cell_batches() const;
1796
1803 unsigned int
1804 n_ghost_cell_batches() const;
1805
1815 unsigned int
1816 n_inner_face_batches() const;
1817
1828 unsigned int
1829 n_boundary_face_batches() const;
1830
1835 unsigned int
1836 n_ghost_inner_face_batches() const;
1837
1848 types::boundary_id
1849 get_boundary_id(const unsigned int face_batch_index) const;
1850
1855 std::array<types::boundary_id, VectorizedArrayType::size()>
1856 get_faces_by_cells_boundary_id(const unsigned int cell_batch_index,
1857 const unsigned int face_number) const;
1858
1863 const DoFHandler<dim> &
1864 get_dof_handler(const unsigned int dof_handler_index = 0) const;
1865
1872 const AffineConstraints<Number> &
1873 get_affine_constraints(const unsigned int dof_handler_index = 0) const;
1874
1887 typename DoFHandler<dim>::cell_iterator
1888 get_cell_iterator(const unsigned int cell_batch_index,
1889 const unsigned int lane_index,
1890 const unsigned int dof_handler_index = 0) const;
1891
1897 std::pair<int, int>
1898 get_cell_level_and_index(const unsigned int cell_batch_index,
1899 const unsigned int lane_index) const;
1900
1907 unsigned int
1908 get_matrix_free_cell_index(
1909 const typename Triangulation<dim>::cell_iterator &cell) const;
1910
1923 std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>
1924 get_face_iterator(const unsigned int face_batch_index,
1925 const unsigned int lane_index,
1926 const bool interior = true,
1927 const unsigned int fe_component = 0) const;
1928
1941 bool
1942 at_irregular_cell(const unsigned int cell_batch_index) const;
1943
1953 unsigned int
1954 n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const;
1955
1965 unsigned int
1966 n_active_entries_per_face_batch(const unsigned int face_batch_index) const;
1967
1971 unsigned int
1972 get_dofs_per_cell(const unsigned int dof_handler_index = 0,
1973 const unsigned int hp_active_fe_index = 0) const;
1974
1978 unsigned int
1979 get_n_q_points(const unsigned int quad_index = 0,
1980 const unsigned int hp_active_fe_index = 0) const;
1981
1986 unsigned int
1987 get_dofs_per_face(const unsigned int dof_handler_index = 0,
1988 const unsigned int hp_active_fe_index = 0) const;
1989
1994 unsigned int
1995 get_n_q_points_face(const unsigned int quad_index = 0,
1996 const unsigned int hp_active_fe_index = 0) const;
1997
2001 const Quadrature<dim> &
2002 get_quadrature(const unsigned int quad_index = 0,
2003 const unsigned int hp_active_fe_index = 0) const;
2004
2008 const Quadrature<dim - 1> &
2009 get_face_quadrature(const unsigned int quad_index = 0,
2010 const unsigned int hp_active_fe_index = 0) const;
2011
2024 unsigned int
2025 get_cell_range_category(
2026 const std::pair<unsigned int, unsigned int> cell_batch_range) const;
2027
2032 std::pair<unsigned int, unsigned int>
2033 get_face_range_category(
2034 const std::pair<unsigned int, unsigned int> face_batch_range) const;
2035
2047 unsigned int
2048 get_cell_category(const unsigned int cell_batch_index) const;
2049
2054 std::pair<unsigned int, unsigned int>
2055 get_face_category(const unsigned int face_batch_index) const;
2056
2060 bool
2061 indices_initialized() const;
2062
2067 bool
2068 mapping_initialized() const;
2069
2074 unsigned int
2075 get_mg_level() const;
2076
2081 std::size_t
2082 memory_consumption() const;
2083
2088 template <typename StreamType>
2089 void
2090 print_memory_consumption(StreamType &out) const;
2091
2096 void
2097 print(std::ostream &out) const;
2098
2110 const internal::MatrixFreeFunctions::TaskInfo &
2111 get_task_info() const;
2112
2113 /*
2114 * Return geometry-dependent information on the cells.
2115 */
2116 const internal::MatrixFreeFunctions::
2117 MappingInfo<dim, Number, VectorizedArrayType> &
2118 get_mapping_info() const;
2119
2123 const internal::MatrixFreeFunctions::DoFInfo &
2124 get_dof_info(const unsigned int dof_handler_index_component = 0) const;
2125
2129 unsigned int
2130 n_constraint_pool_entries() const;
2131
2136 const Number *
2137 constraint_pool_begin(const unsigned int pool_index) const;
2138
2144 const Number *
2145 constraint_pool_end(const unsigned int pool_index) const;
2146
2150 const internal::MatrixFreeFunctions::ShapeInfo<Number> &
2151 get_shape_info(const unsigned int dof_handler_index_component = 0,
2152 const unsigned int quad_index = 0,
2153 const unsigned int fe_base_element = 0,
2154 const unsigned int hp_active_fe_index = 0,
2155 const unsigned int hp_active_quad_index = 0) const;
2156
2160 const internal::MatrixFreeFunctions::FaceToCellTopology<
2161 VectorizedArrayType::size()> &
2162 get_face_info(const unsigned int face_batch_index) const;
2163
2164
2170 const Table<3, unsigned int> &
2171 get_cell_and_face_to_plain_faces() const;
2172
2186 AlignedVector<VectorizedArrayType> *
2187 acquire_scratch_data() const;
2188
2192 void
2193 release_scratch_data(const AlignedVector<VectorizedArrayType> *memory) const;
2194
2204 AlignedVector<Number> *
2205 acquire_scratch_data_non_threadsafe() const;
2206
2210 void
2211 release_scratch_data_non_threadsafe(
2212 const AlignedVector<Number> *memory) const;
2213
2216private:
2221 template <typename number2, int q_dim>
2222 void
2223 internal_reinit(
2224 const std::shared_ptr<hp::MappingCollection<dim>> &mapping,
2225 const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
2226 const std::vector<const AffineConstraints<number2> *> &constraint,
2227 const std::vector<IndexSet> &locally_owned_set,
2228 const std::vector<hp::QCollection<q_dim>> &quad,
2229 const AdditionalData &additional_data);
2230
2237 template <typename number2>
2238 void
2239 initialize_indices(
2240 const std::vector<const AffineConstraints<number2> *> &constraint,
2241 const std::vector<IndexSet> &locally_owned_set,
2242 const AdditionalData &additional_data);
2243
2247 void
2248 initialize_dof_handlers(
2249 const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
2250 const AdditionalData &additional_data);
2251
2255 std::vector<ObserverPointer<const DoFHandler<dim>>> dof_handlers;
2256
2263 std::vector<ObserverPointer<const AffineConstraints<Number>>>
2264 affine_constraints;
2265
2270 std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
2271
2278 std::vector<Number> constraint_pool_data;
2279
2284 std::vector<unsigned int> constraint_pool_row_index;
2285
2290 internal::MatrixFreeFunctions::MappingInfo<dim, Number, VectorizedArrayType>
2291 mapping_info;
2292
2296 Table<4, internal::MatrixFreeFunctions::ShapeInfo<Number>> shape_info;
2297
2304 std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
2305
2310 std::vector<unsigned int> mf_cell_indices;
2311
2318 unsigned int cell_level_index_end_local;
2319
2325 internal::MatrixFreeFunctions::TaskInfo task_info;
2326
2331 internal::MatrixFreeFunctions::FaceInfo<VectorizedArrayType::size()>
2332 face_info;
2333
2337 bool indices_are_initialized;
2338
2342 bool mapping_is_initialized;
2343
2351 mutable Threads::ThreadLocalStorage<
2352 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2353 scratch_pad;
2354
2359 mutable std::list<std::pair<bool, AlignedVector<Number>>>
2360 scratch_pad_non_threadsafe;
2361
2365 unsigned int mg_level;
2366
2371 unsigned int first_hp_dof_handler_index;
2372};
2373
2374
2375
2376/*----------------------- Inline functions ----------------------------------*/
2377
2378#ifndef DOXYGEN
2379
2380
2381
2382template <int dim, typename Number, typename VectorizedArrayType>
2383template <typename T>
2384inline void
2385MatrixFree<dim, Number, VectorizedArrayType>::initialize_cell_data_vector(
2386 AlignedVector<T> &vec) const
2387{
2388 vec.resize(this->n_cell_batches() + this->n_ghost_cell_batches());
2389}
2390
2391
2392
2393template <int dim, typename Number, typename VectorizedArrayType>
2394template <typename T>
2395inline void
2396MatrixFree<dim, Number, VectorizedArrayType>::initialize_face_data_vector(
2397 AlignedVector<T> &vec) const
2398{
2399 vec.resize(this->n_inner_face_batches() + this->n_boundary_face_batches() +
2400 this->n_ghost_inner_face_batches());
2401}
2402
2403
2404
2405template <int dim, typename Number, typename VectorizedArrayType>
2406template <typename VectorType>
2407inline void
2408MatrixFree<dim, Number, VectorizedArrayType>::initialize_dof_vector(
2409 VectorType &vec,
2410 const unsigned int comp) const
2411{
2412 static_assert(IsBlockVector<VectorType>::value == false,
2413 "This function is not supported for block vectors.");
2414
2415 Assert(task_info.n_procs == 1,
2416 ExcMessage("This function can only be used in serial."));
2417
2418 AssertIndexRange(comp, n_components());
2419 vec.reinit(dof_info[comp].vector_partitioner->size());
2420}
2421
2422
2423
2424template <int dim, typename Number, typename VectorizedArrayType>
2425template <typename Number2>
2426inline void
2427MatrixFree<dim, Number, VectorizedArrayType>::initialize_dof_vector(
2428 LinearAlgebra::distributed::Vector<Number2> &vec,
2429 const unsigned int comp) const
2430{
2431 AssertIndexRange(comp, n_components());
2432 vec.reinit(dof_info[comp].vector_partitioner, task_info.communicator_sm);
2433}
2434
2435
2436
2437template <int dim, typename Number, typename VectorizedArrayType>
2438inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2439MatrixFree<dim, Number, VectorizedArrayType>::get_vector_partitioner(
2440 const unsigned int comp) const
2441{
2442 AssertIndexRange(comp, n_components());
2443 return dof_info[comp].vector_partitioner;
2444}
2445
2446
2447
2448template <int dim, typename Number, typename VectorizedArrayType>
2449inline const std::vector<unsigned int> &
2450MatrixFree<dim, Number, VectorizedArrayType>::get_constrained_dofs(
2451 const unsigned int comp) const
2452{
2453 AssertIndexRange(comp, n_components());
2454 return dof_info[comp].constrained_dofs;
2455}
2456
2457
2458
2459template <int dim, typename Number, typename VectorizedArrayType>
2460inline unsigned int
2461MatrixFree<dim, Number, VectorizedArrayType>::n_components() const
2462{
2463 AssertDimension(dof_handlers.size(), dof_info.size());
2464 return dof_handlers.size();
2465}
2466
2467
2468
2469template <int dim, typename Number, typename VectorizedArrayType>
2470inline unsigned int
2471MatrixFree<dim, Number, VectorizedArrayType>::n_base_elements(
2472 const unsigned int dof_no) const
2473{
2474 AssertDimension(dof_handlers.size(), dof_info.size());
2475 AssertIndexRange(dof_no, dof_handlers.size());
2476 return dof_handlers[dof_no]->get_fe().n_base_elements();
2477}
2478
2479
2480
2481template <int dim, typename Number, typename VectorizedArrayType>
2482inline const internal::MatrixFreeFunctions::TaskInfo &
2483MatrixFree<dim, Number, VectorizedArrayType>::get_task_info() const
2484{
2485 return task_info;
2486}
2487
2488
2489
2490template <int dim, typename Number, typename VectorizedArrayType>
2491inline unsigned int
2492MatrixFree<dim, Number, VectorizedArrayType>::n_physical_cells() const
2493{
2494 return task_info.n_active_cells;
2495}
2496
2497
2498
2499template <int dim, typename Number, typename VectorizedArrayType>
2500inline unsigned int
2501MatrixFree<dim, Number, VectorizedArrayType>::n_cell_batches() const
2502{
2503 return *(task_info.cell_partition_data.end() - 2);
2504}
2505
2506
2507
2508template <int dim, typename Number, typename VectorizedArrayType>
2509inline unsigned int
2510MatrixFree<dim, Number, VectorizedArrayType>::n_ghost_cell_batches() const
2511{
2512 return *(task_info.cell_partition_data.end() - 1) -
2513 *(task_info.cell_partition_data.end() - 2);
2514}
2515
2516
2517
2518template <int dim, typename Number, typename VectorizedArrayType>
2519inline unsigned int
2520MatrixFree<dim, Number, VectorizedArrayType>::n_inner_face_batches() const
2521{
2522 if (task_info.face_partition_data.empty())
2523 return 0;
2524 return task_info.face_partition_data.back();
2525}
2526
2527
2528
2529template <int dim, typename Number, typename VectorizedArrayType>
2530inline unsigned int
2531MatrixFree<dim, Number, VectorizedArrayType>::n_boundary_face_batches() const
2532{
2533 if (task_info.face_partition_data.empty())
2534 return 0;
2535 return task_info.boundary_partition_data.back() -
2536 task_info.face_partition_data.back();
2537}
2538
2539
2540
2541template <int dim, typename Number, typename VectorizedArrayType>
2542inline unsigned int
2543MatrixFree<dim, Number, VectorizedArrayType>::n_ghost_inner_face_batches() const
2544{
2545 if (task_info.face_partition_data.empty())
2546 return 0;
2547 return face_info.faces.size() - task_info.boundary_partition_data.back();
2548}
2549
2550
2551
2552template <int dim, typename Number, typename VectorizedArrayType>
2553inline types::boundary_id
2554MatrixFree<dim, Number, VectorizedArrayType>::get_boundary_id(
2555 const unsigned int face_batch_index) const
2556{
2557 Assert(face_batch_index >= task_info.boundary_partition_data[0] &&
2558 face_batch_index < task_info.boundary_partition_data.back(),
2559 ExcIndexRange(face_batch_index,
2560 task_info.boundary_partition_data[0],
2561 task_info.boundary_partition_data.back()));
2562 return types::boundary_id(face_info.faces[face_batch_index].exterior_face_no);
2563}
2564
2565
2566
2567template <int dim, typename Number, typename VectorizedArrayType>
2568inline std::array<types::boundary_id, VectorizedArrayType::size()>
2569MatrixFree<dim, Number, VectorizedArrayType>::get_faces_by_cells_boundary_id(
2570 const unsigned int cell_batch_index,
2571 const unsigned int face_number) const
2572{
2573 AssertIndexRange(cell_batch_index, n_cell_batches());
2574 AssertIndexRange(face_number, GeometryInfo<dim>::faces_per_cell);
2575 Assert(face_info.cell_and_face_boundary_id.size(0) >= n_cell_batches(),
2576 ExcNotInitialized());
2577 std::array<types::boundary_id, VectorizedArrayType::size()> result;
2578 result.fill(numbers::invalid_boundary_id);
2579 for (unsigned int v = 0;
2580 v < n_active_entries_per_cell_batch(cell_batch_index);
2581 ++v)
2582 result[v] =
2583 face_info.cell_and_face_boundary_id(cell_batch_index, face_number, v);
2584 return result;
2585}
2586
2587
2588
2589template <int dim, typename Number, typename VectorizedArrayType>
2590inline const internal::MatrixFreeFunctions::
2591 MappingInfo<dim, Number, VectorizedArrayType> &
2592 MatrixFree<dim, Number, VectorizedArrayType>::get_mapping_info() const
2593{
2594 return mapping_info;
2595}
2596
2597
2598
2599template <int dim, typename Number, typename VectorizedArrayType>
2600inline const internal::MatrixFreeFunctions::DoFInfo &
2601MatrixFree<dim, Number, VectorizedArrayType>::get_dof_info(
2602 const unsigned int dof_index) const
2603{
2604 AssertIndexRange(dof_index, n_components());
2605 return dof_info[dof_index];
2606}
2607
2608
2609
2610template <int dim, typename Number, typename VectorizedArrayType>
2611inline unsigned int
2612MatrixFree<dim, Number, VectorizedArrayType>::n_constraint_pool_entries() const
2613{
2614 return constraint_pool_row_index.size() - 1;
2615}
2616
2617
2618
2619template <int dim, typename Number, typename VectorizedArrayType>
2620inline const Number *
2621MatrixFree<dim, Number, VectorizedArrayType>::constraint_pool_begin(
2622 const unsigned int row) const
2623{
2624 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2625 return constraint_pool_data.empty() ?
2626 nullptr :
2627 constraint_pool_data.data() + constraint_pool_row_index[row];
2628}
2629
2630
2631
2632template <int dim, typename Number, typename VectorizedArrayType>
2633inline const Number *
2634MatrixFree<dim, Number, VectorizedArrayType>::constraint_pool_end(
2635 const unsigned int row) const
2636{
2637 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2638 return constraint_pool_data.empty() ?
2639 nullptr :
2640 constraint_pool_data.data() + constraint_pool_row_index[row + 1];
2641}
2642
2643
2644
2645template <int dim, typename Number, typename VectorizedArrayType>
2646inline std::pair<unsigned int, unsigned int>
2647MatrixFree<dim, Number, VectorizedArrayType>::create_cell_subrange_hp(
2648 const std::pair<unsigned int, unsigned int> &range,
2649 const unsigned int degree,
2650 const unsigned int dof_handler_component) const
2651{
2652 if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2653 {
2654 AssertDimension(
2655 dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2656 AssertDimension(
2657 dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2658 if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2659 return range;
2660 else
2661 return {range.second, range.second};
2662 }
2663
2664 const unsigned int fe_index =
2665 dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2666 if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2667 return {range.second, range.second};
2668 else
2669 return create_cell_subrange_hp_by_index(range,
2670 fe_index,
2671 dof_handler_component);
2672}
2673
2674
2675
2676template <int dim, typename Number, typename VectorizedArrayType>
2677inline bool
2678MatrixFree<dim, Number, VectorizedArrayType>::at_irregular_cell(
2679 const unsigned int cell_batch_index) const
2680{
2681 AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2682 return VectorizedArrayType::size() > 1 &&
2683 cell_level_index[(cell_batch_index + 1) * VectorizedArrayType::size() -
2684 1] == cell_level_index[(cell_batch_index + 1) *
2685 VectorizedArrayType::size() -
2686 2];
2687}
2688
2689
2690
2691template <int dim, typename Number, typename VectorizedArrayType>
2692unsigned int
2693MatrixFree<dim, Number, VectorizedArrayType>::n_active_fe_indices() const
2694{
2695 return shape_info.size(2);
2696}
2697
2698
2699template <int dim, typename Number, typename VectorizedArrayType>
2700unsigned int
2701MatrixFree<dim, Number, VectorizedArrayType>::get_cell_active_fe_index(
2702 const std::pair<unsigned int, unsigned int> range) const
2703{
2704 const auto &fe_indices =
2705 dof_info[first_hp_dof_handler_index].cell_active_fe_index;
2706
2707 if (fe_indices.empty() == true ||
2708 dof_handlers[0]->get_fe_collection().size() == 1)
2709 return 0;
2710
2711 const auto index = fe_indices[range.first];
2712
2713 for (unsigned int i = range.first; i < range.second; ++i)
2714 AssertDimension(index, fe_indices[i]);
2715
2716 return index;
2717}
2718
2719
2720
2721template <int dim, typename Number, typename VectorizedArrayType>
2722unsigned int
2723MatrixFree<dim, Number, VectorizedArrayType>::get_face_active_fe_index(
2724 const std::pair<unsigned int, unsigned int> range,
2725 const bool is_interior_face) const
2726{
2727 const auto &fe_indices =
2728 dof_info[first_hp_dof_handler_index].cell_active_fe_index;
2729
2730 if (fe_indices.empty() == true)
2731 return 0;
2732
2733 if (is_interior_face)
2734 {
2735 const unsigned int index =
2736 fe_indices[face_info.faces[range.first].cells_interior[0] /
2737 VectorizedArrayType::size()];
2738
2739 for (unsigned int i = range.first; i < range.second; ++i)
2740 AssertDimension(index,
2741 fe_indices[face_info.faces[i].cells_interior[0] /
2742 VectorizedArrayType::size()]);
2743
2744 return index;
2745 }
2746 else
2747 {
2748 const unsigned int index =
2749 fe_indices[face_info.faces[range.first].cells_exterior[0] /
2750 VectorizedArrayType::size()];
2751
2752 for (unsigned int i = range.first; i < range.second; ++i)
2753 AssertDimension(index,
2754 fe_indices[face_info.faces[i].cells_exterior[0] /
2755 VectorizedArrayType::size()]);
2756
2757 return index;
2758 }
2759}
2760
2761
2762
2763template <int dim, typename Number, typename VectorizedArrayType>
2764inline unsigned int
2765MatrixFree<dim, Number, VectorizedArrayType>::n_active_entries_per_cell_batch(
2766 const unsigned int cell_batch_index) const
2767{
2768 Assert(!dof_info.empty(), ExcNotInitialized());
2769 AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2770 const std::vector<unsigned char> &n_lanes_filled =
2771 dof_info[0].n_vectorization_lanes_filled
2772 [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell];
2773 AssertIndexRange(cell_batch_index, n_lanes_filled.size());
2774
2775 return n_lanes_filled[cell_batch_index];
2776}
2777
2778
2779
2780template <int dim, typename Number, typename VectorizedArrayType>
2781inline unsigned int
2782MatrixFree<dim, Number, VectorizedArrayType>::n_active_entries_per_face_batch(
2783 const unsigned int face_batch_index) const
2784{
2785 AssertIndexRange(face_batch_index, face_info.faces.size());
2786 Assert(!dof_info.empty(), ExcNotInitialized());
2787 const std::vector<unsigned char> &n_lanes_filled =
2788 dof_info[0].n_vectorization_lanes_filled
2789 [internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior];
2790 AssertIndexRange(face_batch_index, n_lanes_filled.size());
2791 return n_lanes_filled[face_batch_index];
2792}
2793
2794
2795
2796template <int dim, typename Number, typename VectorizedArrayType>
2797inline unsigned int
2798MatrixFree<dim, Number, VectorizedArrayType>::get_dofs_per_cell(
2799 const unsigned int dof_handler_index,
2800 const unsigned int active_fe_index) const
2801{
2802 return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2803}
2804
2805
2806
2807template <int dim, typename Number, typename VectorizedArrayType>
2808inline unsigned int
2809MatrixFree<dim, Number, VectorizedArrayType>::get_n_q_points(
2810 const unsigned int quad_index,
2811 const unsigned int active_fe_index) const
2812{
2813 AssertIndexRange(quad_index, mapping_info.cell_data.size());
2814 return mapping_info.cell_data[quad_index]
2815 .descriptor[active_fe_index]
2816 .n_q_points;
2817}
2818
2819
2820
2821template <int dim, typename Number, typename VectorizedArrayType>
2822inline unsigned int
2823MatrixFree<dim, Number, VectorizedArrayType>::get_dofs_per_face(
2824 const unsigned int dof_handler_index,
2825 const unsigned int active_fe_index) const
2826{
2827 return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2828}
2829
2830
2831
2832template <int dim, typename Number, typename VectorizedArrayType>
2833inline unsigned int
2834MatrixFree<dim, Number, VectorizedArrayType>::get_n_q_points_face(
2835 const unsigned int quad_index,
2836 const unsigned int active_fe_index) const
2837{
2838 AssertIndexRange(quad_index, mapping_info.face_data.size());
2839 return mapping_info.face_data[quad_index]
2840 .descriptor[active_fe_index]
2841 .n_q_points;
2842}
2843
2844
2845
2846template <int dim, typename Number, typename VectorizedArrayType>
2847inline const IndexSet &
2848MatrixFree<dim, Number, VectorizedArrayType>::get_locally_owned_set(
2849 const unsigned int dof_handler_index) const
2850{
2851 return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2852}
2853
2854
2855
2856template <int dim, typename Number, typename VectorizedArrayType>
2857inline const IndexSet &
2858MatrixFree<dim, Number, VectorizedArrayType>::get_ghost_set(
2859 const unsigned int dof_handler_index) const
2860{
2861 return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2862}
2863
2864
2865
2866template <int dim, typename Number, typename VectorizedArrayType>
2867inline const internal::MatrixFreeFunctions::ShapeInfo<Number> &
2868MatrixFree<dim, Number, VectorizedArrayType>::get_shape_info(
2869 const unsigned int dof_handler_index,
2870 const unsigned int index_quad,
2871 const unsigned int index_fe,
2872 const unsigned int active_fe_index,
2873 const unsigned int active_quad_index) const
2874{
2875 AssertIndexRange(dof_handler_index, dof_info.size());
2876 const unsigned int ind =
2877 dof_info[dof_handler_index].global_base_element_offset + index_fe;
2878 AssertIndexRange(ind, shape_info.size(0));
2879 AssertIndexRange(index_quad, shape_info.size(1));
2880 AssertIndexRange(active_fe_index, shape_info.size(2));
2881 AssertIndexRange(active_quad_index, shape_info.size(3));
2882 return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2883}
2884
2885
2886
2887template <int dim, typename Number, typename VectorizedArrayType>
2888inline const internal::MatrixFreeFunctions::FaceToCellTopology<
2889 VectorizedArrayType::size()> &
2890MatrixFree<dim, Number, VectorizedArrayType>::get_face_info(
2891 const unsigned int face_batch_index) const
2892{
2893 AssertIndexRange(face_batch_index, face_info.faces.size());
2894 return face_info.faces[face_batch_index];
2895}
2896
2897
2898
2899template <int dim, typename Number, typename VectorizedArrayType>
2900inline const Table<3, unsigned int> &
2901MatrixFree<dim, Number, VectorizedArrayType>::get_cell_and_face_to_plain_faces()
2902 const
2903{
2904 return face_info.cell_and_face_to_plain_faces;
2905}
2906
2907
2908
2909template <int dim, typename Number, typename VectorizedArrayType>
2910inline const Quadrature<dim> &
2911MatrixFree<dim, Number, VectorizedArrayType>::get_quadrature(
2912 const unsigned int quad_index,
2913 const unsigned int active_fe_index) const
2914{
2915 AssertIndexRange(quad_index, mapping_info.cell_data.size());
2916 return mapping_info.cell_data[quad_index]
2917 .descriptor[active_fe_index]
2918 .quadrature;
2919}
2920
2921
2922
2923template <int dim, typename Number, typename VectorizedArrayType>
2924inline const Quadrature<dim - 1> &
2925MatrixFree<dim, Number, VectorizedArrayType>::get_face_quadrature(
2926 const unsigned int quad_index,
2927 const unsigned int active_fe_index) const
2928{
2929 AssertIndexRange(quad_index, mapping_info.face_data.size());
2930 return mapping_info.face_data[quad_index]
2931 .descriptor[active_fe_index]
2932 .quadrature;
2933}
2934
2935
2936
2937template <int dim, typename Number, typename VectorizedArrayType>
2938inline unsigned int
2939MatrixFree<dim, Number, VectorizedArrayType>::get_cell_range_category(
2940 const std::pair<unsigned int, unsigned int> range) const
2941{
2942 auto result = get_cell_category(range.first);
2943
2944 for (unsigned int i = range.first; i < range.second; ++i)
2945 result = std::max(result, get_cell_category(i));
2946
2947 return result;
2948}
2949
2950
2951
2952template <int dim, typename Number, typename VectorizedArrayType>
2953inline std::pair<unsigned int, unsigned int>
2954MatrixFree<dim, Number, VectorizedArrayType>::get_face_range_category(
2955 const std::pair<unsigned int, unsigned int> range) const
2956{
2957 auto result = get_face_category(range.first);
2958
2959 for (unsigned int i = range.first; i < range.second; ++i)
2960 {
2961 result.first = std::max(result.first, get_face_category(i).first);
2962 result.second = std::max(result.second, get_face_category(i).second);
2963 }
2964
2965 return result;
2966}
2967
2968
2969
2970template <int dim, typename Number, typename VectorizedArrayType>
2971inline unsigned int
2972MatrixFree<dim, Number, VectorizedArrayType>::get_cell_category(
2973 const unsigned int cell_batch_index) const
2974{
2975 AssertIndexRange(0, dof_info.size());
2976 AssertIndexRange(
2977 cell_batch_index,
2978 dof_info[first_hp_dof_handler_index].cell_active_fe_index.size());
2979 if (dof_info[first_hp_dof_handler_index].cell_active_fe_index.empty())
2980 return 0;
2981 else
2982 return dof_info[first_hp_dof_handler_index]
2983 .cell_active_fe_index[cell_batch_index];
2984}
2985
2986
2987
2988template <int dim, typename Number, typename VectorizedArrayType>
2989inline std::pair<unsigned int, unsigned int>
2990MatrixFree<dim, Number, VectorizedArrayType>::get_face_category(
2991 const unsigned int face_batch_index) const
2992{
2993 AssertIndexRange(face_batch_index, face_info.faces.size());
2994 if (dof_info[first_hp_dof_handler_index].cell_active_fe_index.empty())
2995 return std::make_pair(0U, 0U);
2996
2997 std::pair<unsigned int, unsigned int> result = std::make_pair(0U, 0U);
2998 for (unsigned int v = 0;
2999 v < VectorizedArrayType::size() &&
3000 face_info.faces[face_batch_index].cells_interior[v] !=
3001 numbers::invalid_unsigned_int;
3002 ++v)
3003 result.first =
3004 std::max(result.first,
3005 dof_info[first_hp_dof_handler_index].cell_active_fe_index
3006 [face_info.faces[face_batch_index].cells_interior[v] /
3007 VectorizedArrayType::size()]);
3008 if (face_info.faces[face_batch_index].cells_exterior[0] !=
3009 numbers::invalid_unsigned_int)
3010 for (unsigned int v = 0;
3011 v < VectorizedArrayType::size() &&
3012 face_info.faces[face_batch_index].cells_exterior[v] !=
3013 numbers::invalid_unsigned_int;
3014 ++v)
3015 result.second =
3016 std::max(result.second,
3017 dof_info[first_hp_dof_handler_index].cell_active_fe_index
3018 [face_info.faces[face_batch_index].cells_exterior[v] /
3019 VectorizedArrayType::size()]);
3020 else
3021 result.second = numbers::invalid_unsigned_int;
3022 return result;
3023}
3024
3025
3026
3027template <int dim, typename Number, typename VectorizedArrayType>
3028inline bool
3029MatrixFree<dim, Number, VectorizedArrayType>::indices_initialized() const
3030{
3031 return indices_are_initialized;
3032}
3033
3034
3035
3036template <int dim, typename Number, typename VectorizedArrayType>
3037inline bool
3038MatrixFree<dim, Number, VectorizedArrayType>::mapping_initialized() const
3039{
3040 return mapping_is_initialized;
3041}
3042
3043
3044template <int dim, typename Number, typename VectorizedArrayType>
3045inline unsigned int
3046MatrixFree<dim, Number, VectorizedArrayType>::get_mg_level() const
3047{
3048 return mg_level;
3049}
3050
3051
3052
3053template <int dim, typename Number, typename VectorizedArrayType>
3054AlignedVector<VectorizedArrayType> *
3055MatrixFree<dim, Number, VectorizedArrayType>::acquire_scratch_data() const
3056{
3057 using list_type =
3058 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3059 list_type &data = scratch_pad.get();
3060 for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
3061 if (it->first == false)
3062 {
3063 it->first = true;
3064 return &it->second;
3065 }
3066 data.emplace_front(true, AlignedVector<VectorizedArrayType>());
3067 return &data.front().second;
3068}
3069
3070
3071
3072template <int dim, typename Number, typename VectorizedArrayType>
3073void
3074MatrixFree<dim, Number, VectorizedArrayType>::release_scratch_data(
3075 const AlignedVector<VectorizedArrayType> *scratch) const
3076{
3077 using list_type =
3078 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3079 list_type &data = scratch_pad.get();
3080 for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
3081 if (&it->second == scratch)
3082 {
3083 Assert(it->first == true, ExcInternalError());
3084 it->first = false;
3085 return;
3086 }
3087 AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
3088}
3089
3090
3091
3092template <int dim, typename Number, typename VectorizedArrayType>
3093AlignedVector<Number> *
3094MatrixFree<dim, Number, VectorizedArrayType>::
3095 acquire_scratch_data_non_threadsafe() const
3096{
3097 for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
3098 scratch_pad_non_threadsafe.begin();
3099 it != scratch_pad_non_threadsafe.end();
3100 ++it)
3101 if (it->first == false)
3102 {
3103 it->first = true;
3104 return &it->second;
3105 }
3106 scratch_pad_non_threadsafe.push_front(
3107 std::make_pair(true, AlignedVector<Number>()));
3108 return &scratch_pad_non_threadsafe.front().second;
3109}
3110
3111
3112
3113template <int dim, typename Number, typename VectorizedArrayType>
3114void
3115MatrixFree<dim, Number, VectorizedArrayType>::
3116 release_scratch_data_non_threadsafe(
3117 const AlignedVector<Number> *scratch) const
3118{
3119 for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
3120 scratch_pad_non_threadsafe.begin();
3121 it != scratch_pad_non_threadsafe.end();
3122 ++it)
3123 if (&it->second == scratch)
3124 {
3125 Assert(it->first == true, ExcInternalError());
3126 it->first = false;
3127 return;
3128 }
3129 AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
3130}
3131
3132
3133
3134// ------------------------------ reinit functions ---------------------------
3135
3136namespace internal
3137{
3138 namespace MatrixFreeImplementation
3139 {
3140 template <int dim, int spacedim>
3141 inline std::vector<IndexSet>
3142 extract_locally_owned_index_sets(
3143 const std::vector<const DoFHandler<dim, spacedim> *> &dofh,
3144 const unsigned int level)
3145 {
3146 std::vector<IndexSet> locally_owned_set;
3147 locally_owned_set.reserve(dofh.size());
3148 for (unsigned int j = 0; j < dofh.size(); ++j)
3149 if (level == numbers::invalid_unsigned_int)
3150 locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
3151 else
3152 locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
3153 return locally_owned_set;
3154 }
3155 } // namespace MatrixFreeImplementation
3156} // namespace internal
3157
3158
3159
3160template <int dim, typename Number, typename VectorizedArrayType>
3161template <typename QuadratureType, typename number2, typename MappingType>
3162void
3163MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3164 const MappingType &mapping,
3165 const DoFHandler<dim> &dof_handler,
3166 const AffineConstraints<number2> &constraints_in,
3167 const QuadratureType &quad,
3168 const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
3169 &additional_data)
3170{
3171 std::vector<const DoFHandler<dim, dim> *> dof_handlers;
3172 std::vector<const AffineConstraints<number2> *> constraints;
3173
3174 dof_handlers.push_back(&dof_handler);
3175 constraints.push_back(&constraints_in);
3176
3177 std::vector<IndexSet> locally_owned_sets =
3178 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3179 dof_handlers, additional_data.mg_level);
3180
3181 std::vector<hp::QCollection<dim>> quad_hp;
3182 quad_hp.emplace_back(quad);
3183
3184 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3185 dof_handlers,
3186 constraints,
3187 locally_owned_sets,
3188 quad_hp,
3189 additional_data);
3190}
3191
3192
3193
3194template <int dim, typename Number, typename VectorizedArrayType>
3195template <typename QuadratureType, typename number2, typename MappingType>
3196void
3197MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3198 const MappingType &mapping,
3199 const std::vector<const DoFHandler<dim> *> &dof_handler,
3200 const std::vector<const AffineConstraints<number2> *> &constraint,
3201 const QuadratureType &quad,
3202 const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
3203 &additional_data)
3204{
3205 std::vector<IndexSet> locally_owned_set =
3206 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3207 dof_handler, additional_data.mg_level);
3208 std::vector<hp::QCollection<dim>> quad_hp;
3209 quad_hp.emplace_back(quad);
3210
3211 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3212 dof_handler,
3213 constraint,
3214 locally_owned_set,
3215 quad_hp,
3216 additional_data);
3217}
3218
3219
3220
3221template <int dim, typename Number, typename VectorizedArrayType>
3222template <typename QuadratureType, typename number2, typename MappingType>
3223void
3224MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3225 const MappingType &mapping,
3226 const std::vector<const DoFHandler<dim> *> &dof_handler,
3227 const std::vector<const AffineConstraints<number2> *> &constraint,
3228 const std::vector<QuadratureType> &quad,
3229 const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
3230 &additional_data)
3231{
3232 std::vector<IndexSet> locally_owned_set =
3233 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3234 dof_handler, additional_data.mg_level);
3235 std::vector<hp::QCollection<dim>> quad_hp;
3236 for (unsigned int q = 0; q < quad.size(); ++q)
3237 quad_hp.emplace_back(quad[q]);
3238
3239 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3240 dof_handler,
3241 constraint,
3242 locally_owned_set,
3243 quad_hp,
3244 additional_data);
3245}
3246
3247
3248
3249// ------------------------------ implementation of loops --------------------
3250
3251// internal helper functions that define how to call MPI data exchange
3252// functions: for generic vectors, do nothing at all. For distributed vectors,
3253// call update_ghost_values_start function and so on. If we have collections
3254// of vectors, just do the individual functions of the components. In order to
3255// keep ghost values consistent (whether we are in read or write mode), we
3256// also reset the values at the end. the whole situation is a bit complicated
3257// by the fact that we need to treat block vectors differently, which use some
3258// additional helper functions to select the blocks and template magic.
3259namespace internal
3260{
3264 template <int dim, typename Number, typename VectorizedArrayType>
3265 struct VectorDataExchange
3266 {
3267 // A shift for the MPI messages to reduce the risk for accidental
3268 // interaction with other open communications that a user program might
3269 // set up (parallel vectors support unfinished communication). We let
3270 // the other vectors use the first 20 assigned numbers and start the
3271 // matrix-free communication.
3272 static constexpr unsigned int channel_shift = 20;
3273
3274
3275
3280 VectorDataExchange(
3281 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3282 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3283 DataAccessOnFaces vector_face_access,
3284 const unsigned int n_components)
3285 : matrix_free(matrix_free)
3286 , vector_face_access(
3287 matrix_free.get_task_info().face_partition_data.empty() ?
3288 ::MatrixFree<dim, Number, VectorizedArrayType>::
3289 DataAccessOnFaces::unspecified :
3290 vector_face_access)
3291 , ghosts_were_set(false)
3292# ifdef DEAL_II_WITH_MPI
3293 , tmp_data(n_components)
3294 , requests(n_components)
3295# endif
3296 {
3297 (void)n_components;
3298 if (this->vector_face_access !=
3299 ::MatrixFree<dim, Number, VectorizedArrayType>::
3300 DataAccessOnFaces::unspecified)
3301 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3302 AssertDimension(
3303 matrix_free.get_dof_info(c).vector_exchanger_face_variants.size(),
3304 5);
3305 }
3306
3307
3308
3312 ~VectorDataExchange() // NOLINT
3313 {
3314# ifdef DEAL_II_WITH_MPI
3315 for (unsigned int i = 0; i < tmp_data.size(); ++i)
3316 if (tmp_data[i] != nullptr)
3317 matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3318# endif
3319 }
3320
3321
3322
3327 template <typename VectorType>
3328 unsigned int
3329 find_vector_in_mf(const VectorType &vec,
3330 const bool check_global_compatibility = true) const
3331 {
3332 // case 1: vector was set up with MatrixFree::initialize_dof_vector()
3333 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3334 if (vec.get_partitioner().get() ==
3335 matrix_free.get_dof_info(c).vector_partitioner.get())
3336 return c;
3337
3338 // case 2: user provided own partitioner (compatibility mode)
3339 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3340 if (check_global_compatibility ?
3341 vec.get_partitioner()->is_globally_compatible(
3342 *matrix_free.get_dof_info(c).vector_partitioner) :
3343 vec.get_partitioner()->is_compatible(
3344 *matrix_free.get_dof_info(c).vector_partitioner))
3345 return c;
3346
3347 Assert(false,
3348 ExcNotImplemented("Could not find partitioner that fits vector"));
3349
3350 return numbers::invalid_unsigned_int;
3351 }
3352
3353
3354
3359 const internal::MatrixFreeFunctions::VectorDataExchange::Base &
3360 get_partitioner(const unsigned int mf_component) const
3361 {
3362 AssertDimension(matrix_free.get_dof_info(mf_component)
3363 .vector_exchanger_face_variants.size(),
3364 5);
3365 if (vector_face_access ==
3366 ::MatrixFree<dim, Number, VectorizedArrayType>::
3367 DataAccessOnFaces::none)
3368 return *matrix_free.get_dof_info(mf_component)
3369 .vector_exchanger_face_variants[0];
3370 else if (vector_face_access ==
3371 ::MatrixFree<dim, Number, VectorizedArrayType>::
3372 DataAccessOnFaces::values)
3373 return *matrix_free.get_dof_info(mf_component)
3374 .vector_exchanger_face_variants[1];
3375 else if (vector_face_access ==
3376 ::MatrixFree<dim, Number, VectorizedArrayType>::
3377 DataAccessOnFaces::gradients)
3378 return *matrix_free.get_dof_info(mf_component)
3379 .vector_exchanger_face_variants[2];
3380 else if (vector_face_access ==
3381 ::MatrixFree<dim, Number, VectorizedArrayType>::
3382 DataAccessOnFaces::values_all_faces)
3383 return *matrix_free.get_dof_info(mf_component)
3384 .vector_exchanger_face_variants[3];
3385 else if (vector_face_access ==
3386 ::MatrixFree<dim, Number, VectorizedArrayType>::
3387 DataAccessOnFaces::gradients_all_faces)
3388 return *matrix_free.get_dof_info(mf_component)
3389 .vector_exchanger_face_variants[4];
3390 else
3391 return *matrix_free.get_dof_info(mf_component).vector_exchanger.get();
3392 }
3393
3394
3395
3399 template <typename VectorType,
3400 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3401 * = nullptr>
3402 void
3403 update_ghost_values_start(const unsigned int /*component_in_block_vector*/,
3404 const VectorType & /*vec*/)
3405 {}
3406
3407
3412 template <typename VectorType,
3413 std::enable_if_t<!has_update_ghost_values_start<VectorType> &&
3414 !is_not_parallel_vector<VectorType>,
3415 VectorType> * = nullptr>
3416 void
3417 update_ghost_values_start(const unsigned int component_in_block_vector,
3418 const VectorType &vec)
3419 {
3420 (void)component_in_block_vector;
3421 const bool ghosts_set = vec.has_ghost_elements();
3422
3423 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3424 ghosts_set == false,
3425 ExcNotImplemented());
3426
3427 if (ghosts_set)
3428 {
3429 ghosts_were_set = true;
3430 return;
3431 }
3432
3433 vec.update_ghost_values();
3434 }
3435
3436
3437
3443 template <typename VectorType,
3444 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3445 !has_exchange_on_subset<VectorType>,
3446 VectorType> * = nullptr>
3447 void
3448 update_ghost_values_start(const unsigned int component_in_block_vector,
3449 const VectorType &vec)
3450 {
3451 (void)component_in_block_vector;
3452 const bool ghosts_set = vec.has_ghost_elements();
3453
3454 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3455 ghosts_set == false,
3456 ExcNotImplemented());
3457
3458 if (ghosts_set)
3459 {
3460 ghosts_were_set = true;
3461 return;
3462 }
3463
3464 vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3465 }
3466
3467
3468
3475 template <typename VectorType,
3476 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3477 has_exchange_on_subset<VectorType>,
3478 VectorType> * = nullptr>
3479 void
3480 update_ghost_values_start(const unsigned int component_in_block_vector,
3481 const VectorType &vec)
3482 {
3483 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3484 "Type mismatch between VectorType and VectorDataExchange");
3485 (void)component_in_block_vector;
3486 const bool ghosts_set = vec.has_ghost_elements();
3487
3488 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3489 ghosts_set == false,
3490 ExcNotImplemented());
3491
3492 if (ghosts_set)
3493 {
3494 ghosts_were_set = true;
3495 return;
3496 }
3497
3498 if (vec.size() != 0)
3499 {
3500# ifdef DEAL_II_WITH_MPI
3501 const unsigned int mf_component = find_vector_in_mf(vec);
3502
3503 const auto &part = get_partitioner(mf_component);
3504
3505 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3506 part.n_import_sm_procs() == 0)
3507 return;
3508
3509 tmp_data[component_in_block_vector] =
3510 matrix_free.acquire_scratch_data_non_threadsafe();
3511 tmp_data[component_in_block_vector]->resize_fast(
3512 part.n_import_indices());
3513 AssertDimension(requests.size(), tmp_data.size());
3514
3515 part.export_to_ghosted_array_start(
3516 component_in_block_vector * 2 + channel_shift,
3517 ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3518 vec.shared_vector_data(),
3519 ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3520 part.locally_owned_size(),
3521 matrix_free.get_dof_info(mf_component)
3522 .vector_partitioner->n_ghost_indices()),
3523 ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3524 part.n_import_indices()),
3525 this->requests[component_in_block_vector]);
3526# endif
3527 }
3528 }
3529
3530
3531
3536 template <typename VectorType,
3537 std::enable_if_t<!has_update_ghost_values_start<VectorType>,
3538 VectorType> * = nullptr>
3539 void
3540 update_ghost_values_finish(const unsigned int /*component_in_block_vector*/,
3541 const VectorType & /*vec*/)
3542 {}
3543
3544
3545
3551 template <typename VectorType,
3552 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3553 !has_exchange_on_subset<VectorType>,
3554 VectorType> * = nullptr>
3555 void
3556 update_ghost_values_finish(const unsigned int component_in_block_vector,
3557 const VectorType &vec)
3558 {
3559 (void)component_in_block_vector;
3560
3561 if (ghosts_were_set)
3562 return;
3563
3564 vec.update_ghost_values_finish();
3565 }
3566
3567
3568
3575 template <typename VectorType,
3576 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3577 has_exchange_on_subset<VectorType>,
3578 VectorType> * = nullptr>
3579 void
3580 update_ghost_values_finish(const unsigned int component_in_block_vector,
3581 const VectorType &vec)
3582 {
3583 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3584 "Type mismatch between VectorType and VectorDataExchange");
3585 (void)component_in_block_vector;
3586
3587 if (ghosts_were_set)
3588 return;
3589
3590 if (vec.size() != 0)
3591 {
3592# ifdef DEAL_II_WITH_MPI
3593 AssertIndexRange(component_in_block_vector, tmp_data.size());
3594 AssertDimension(requests.size(), tmp_data.size());
3595
3596 const unsigned int mf_component = find_vector_in_mf(vec);
3597
3598 const auto &part = get_partitioner(mf_component);
3599
3600 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3601 part.n_import_sm_procs() != 0)
3602 {
3603 part.export_to_ghosted_array_finish(
3604 ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3605 vec.shared_vector_data(),
3606 ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3607 part.locally_owned_size(),
3608 matrix_free.get_dof_info(mf_component)
3609 .vector_partitioner->n_ghost_indices()),
3610 this->requests[component_in_block_vector]);
3611
3612 matrix_free.release_scratch_data_non_threadsafe(
3613 tmp_data[component_in_block_vector]);
3614 tmp_data[component_in_block_vector] = nullptr;
3615 }
3616# endif
3617 }
3618 // let vector know that ghosts are being updated and we can read from
3619 // them
3620 vec.set_ghost_state(true);
3621 }
3622
3623
3624
3628 template <typename VectorType,
3629 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3630 * = nullptr>
3631 void
3632 compress_start(const unsigned int /*component_in_block_vector*/,
3633 VectorType & /*vec*/)
3634 {}
3635
3636
3637
3642 template <typename VectorType,
3643 std::enable_if_t<!has_compress_start<VectorType> &&
3644 !is_not_parallel_vector<VectorType>,
3645 VectorType> * = nullptr>
3646 void
3647 compress_start(const unsigned int component_in_block_vector,
3648 VectorType &vec)
3649 {
3650 (void)component_in_block_vector;
3651 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3652 vec.compress(VectorOperation::add);
3653 }
3654
3655
3656
3662 template <typename VectorType,
3663 std::enable_if_t<has_compress_start<VectorType> &&
3664 !has_exchange_on_subset<VectorType>,
3665 VectorType> * = nullptr>
3666 void
3667 compress_start(const unsigned int component_in_block_vector,
3668 VectorType &vec)
3669 {
3670 (void)component_in_block_vector;
3671 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3672 vec.compress_start(component_in_block_vector + channel_shift);
3673 }
3674
3675
3676
3683 template <typename VectorType,
3684 std::enable_if_t<has_compress_start<VectorType> &&
3685 has_exchange_on_subset<VectorType>,
3686 VectorType> * = nullptr>
3687 void
3688 compress_start(const unsigned int component_in_block_vector,
3689 VectorType &vec)
3690 {
3691 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3692 "Type mismatch between VectorType and VectorDataExchange");
3693 (void)component_in_block_vector;
3694 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3695
3696 if (vec.size() != 0)
3697 {
3698# ifdef DEAL_II_WITH_MPI
3699 const unsigned int mf_component = find_vector_in_mf(vec);
3700
3701 const auto &part = get_partitioner(mf_component);
3702
3703 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3704 part.n_import_sm_procs() == 0)
3705 return;
3706
3707 tmp_data[component_in_block_vector] =
3708 matrix_free.acquire_scratch_data_non_threadsafe();
3709 tmp_data[component_in_block_vector]->resize_fast(
3710 part.n_import_indices());
3711 AssertDimension(requests.size(), tmp_data.size());
3712
3713 part.import_from_ghosted_array_start(
3714 VectorOperation::add,
3715 component_in_block_vector * 2 + channel_shift,
3716 ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3717 vec.shared_vector_data(),
3718 ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3719 matrix_free.get_dof_info(mf_component)
3720 .vector_partitioner->n_ghost_indices()),
3721 ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3722 part.n_import_indices()),
3723 this->requests[component_in_block_vector]);
3724# endif
3725 }
3726 }
3727
3728
3729
3734 template <
3735 typename VectorType,
3736 std::enable_if_t<!has_compress_start<VectorType>, VectorType> * = nullptr>
3737 void
3738 compress_finish(const unsigned int /*component_in_block_vector*/,
3739 VectorType & /*vec*/)
3740 {}
3741
3742
3743
3749 template <typename VectorType,
3750 std::enable_if_t<has_compress_start<VectorType> &&
3751 !has_exchange_on_subset<VectorType>,
3752 VectorType> * = nullptr>
3753 void
3754 compress_finish(const unsigned int component_in_block_vector,
3755 VectorType &vec)
3756 {
3757 (void)component_in_block_vector;
3758 vec.compress_finish(VectorOperation::add);
3759 }
3760
3761
3762
3769 template <typename VectorType,
3770 std::enable_if_t<has_compress_start<VectorType> &&
3771 has_exchange_on_subset<VectorType>,
3772 VectorType> * = nullptr>
3773 void
3774 compress_finish(const unsigned int component_in_block_vector,
3775 VectorType &vec)
3776 {
3777 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3778 "Type mismatch between VectorType and VectorDataExchange");
3779 (void)component_in_block_vector;
3780 if (vec.size() != 0)
3781 {
3782# ifdef DEAL_II_WITH_MPI
3783 AssertIndexRange(component_in_block_vector, tmp_data.size());
3784 AssertDimension(requests.size(), tmp_data.size());
3785
3786 const unsigned int mf_component = find_vector_in_mf(vec);
3787
3788 const auto &part = get_partitioner(mf_component);
3789
3790 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3791 part.n_import_sm_procs() != 0)
3792 {
3793 part.import_from_ghosted_array_finish(
3794 VectorOperation::add,
3795 ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3796 vec.shared_vector_data(),
3797 ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3798 matrix_free.get_dof_info(mf_component)
3799 .vector_partitioner->n_ghost_indices()),
3800 ArrayView<const Number>(
3801 tmp_data[component_in_block_vector]->begin(),
3802 part.n_import_indices()),
3803 this->requests[component_in_block_vector]);
3804
3805 matrix_free.release_scratch_data_non_threadsafe(
3806 tmp_data[component_in_block_vector]);
3807 tmp_data[component_in_block_vector] = nullptr;
3808 }
3809
3810 if (Utilities::MPI::job_supports_mpi())
3811 {
3812 const int ierr =
3813 MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3814 AssertThrowMPI(ierr);
3815 }
3816# endif
3817 }
3818 }
3819
3820
3821
3825 template <typename VectorType,
3826 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3827 * = nullptr>
3828 void
3829 reset_ghost_values(const VectorType & /*vec*/) const
3830 {}
3831
3832
3833
3838 template <typename VectorType,
3839 std::enable_if_t<!has_exchange_on_subset<VectorType> &&
3840 !is_not_parallel_vector<VectorType>,
3841 VectorType> * = nullptr>
3842 void
3843 reset_ghost_values(const VectorType &vec) const
3844 {
3845 if (ghosts_were_set == true)
3846 return;
3847
3848 vec.zero_out_ghost_values();
3849 }
3850
3851
3852
3858 template <typename VectorType,
3859 std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3860 * = nullptr>
3861 void
3862 reset_ghost_values(const VectorType &vec) const
3863 {
3864 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3865 "Type mismatch between VectorType and VectorDataExchange");
3866 if (ghosts_were_set == true)
3867 return;
3868
3869 if (vec.size() != 0)
3870 {
3871# ifdef DEAL_II_WITH_MPI
3872 AssertDimension(requests.size(), tmp_data.size());
3873
3874 const unsigned int mf_component = find_vector_in_mf(vec);
3875
3876 const auto &part = get_partitioner(mf_component);
3877
3878 if (part.n_ghost_indices() > 0)
3879 {
3880 part.reset_ghost_values(
3881 ArrayView<Number>(const_cast<VectorType &>(vec).begin() +
3882 part.locally_owned_size(),
3883 matrix_free.get_dof_info(mf_component)
3884 .vector_partitioner->n_ghost_indices()));
3885 }
3886
3887# endif
3888 }
3889 // let vector know that it's not ghosted anymore
3890 vec.set_ghost_state(false);
3891 }
3892
3893
3894
3900 template <typename VectorType,
3901 std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3902 * = nullptr>
3903 void
3904 zero_vector_region(const unsigned int range_index, VectorType &vec) const
3905 {
3906 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3907 "Type mismatch between VectorType and VectorDataExchange");
3908 if (range_index == numbers::invalid_unsigned_int)
3909 vec = Number();
3910 else
3911 {
3912 const unsigned int mf_component = find_vector_in_mf(vec, false);
3913 const internal::MatrixFreeFunctions::DoFInfo &dof_info =
3914 matrix_free.get_dof_info(mf_component);
3915 Assert(dof_info.vector_zero_range_list_index.empty() == false,
3916 ExcNotInitialized());
3917
3918 Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3919 ExcInternalError());
3920 AssertIndexRange(range_index,
3921 dof_info.vector_zero_range_list_index.size() - 1);
3922 for (unsigned int id =
3923 dof_info.vector_zero_range_list_index[range_index];
3924 id != dof_info.vector_zero_range_list_index[range_index + 1];
3925 ++id)
3926 std::memset(vec.begin() + dof_info.vector_zero_range_list[id].first,
3927 0,
3928 (dof_info.vector_zero_range_list[id].second -
3929 dof_info.vector_zero_range_list[id].first) *
3930 sizeof(Number));
3931 }
3932 }
3933
3934
3935
3941 template <typename VectorType,
3942 std::enable_if_t<!has_exchange_on_subset<VectorType>, VectorType>
3943 * = nullptr,
3944 typename VectorType::value_type * = nullptr>
3945 void
3946 zero_vector_region(const unsigned int range_index, VectorType &vec) const
3947 {
3948 if (range_index == numbers::invalid_unsigned_int || range_index == 0)
3949 {
3950 if constexpr (std::is_same_v<
3951 ArrayView<typename VectorType::value_type>,
3952 VectorType>)
3953 {
3954 for (unsigned int i = 0; i < vec.size(); ++i)
3955 vec[i] = typename VectorType::value_type();
3956 }
3957 else
3958 vec = typename VectorType::value_type();
3959 }
3960 }
3961
3962
3963
3968 void
3969 zero_vector_region(const unsigned int, ...) const
3970 {
3971 Assert(false,
3972 ExcNotImplemented("Zeroing is only implemented for vector types "
3973 "which provide VectorType::value_type"));
3974 }
3975
3976
3977
3978 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
3979 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3980 DataAccessOnFaces vector_face_access;
3981 bool ghosts_were_set;
3982# ifdef DEAL_II_WITH_MPI
3983 std::vector<AlignedVector<Number> *> tmp_data;
3984 std::vector<std::vector<MPI_Request>> requests;
3985# endif
3986 }; // VectorDataExchange
3987
3988 template <typename VectorStruct>
3989 unsigned int
3990 n_components(const VectorStruct &vec);
3991
3992 template <typename VectorStruct>
3993 unsigned int
3994 n_components_block(const VectorStruct &vec, const std::bool_constant<true>)
3995 {
3996 unsigned int components = 0;
3997 for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
3998 components += n_components(vec.block(bl));
3999 return components;
4000 }
4001
4002 template <typename VectorStruct>
4003 unsigned int
4004 n_components_block(const VectorStruct &, const std::bool_constant<false>)
4005 {
4006 return 1;
4007 }
4008
4009 template <typename VectorStruct>
4010 unsigned int
4011 n_components(const VectorStruct &vec)
4012 {
4013 return n_components_block(
4014 vec, std::bool_constant<IsBlockVector<VectorStruct>::value>());
4015 }
4016
4017 template <typename VectorStruct>
4018 inline unsigned int
4019 n_components(const std::vector<VectorStruct> &vec)
4020 {
4021 unsigned int components = 0;
4022 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4023 components += n_components_block(
4024 vec[comp], std::bool_constant<IsBlockVector<VectorStruct>::value>());
4025 return components;
4026 }
4027
4028 template <typename VectorStruct>
4029 inline unsigned int
4030 n_components(const std::vector<VectorStruct *> &vec)
4031 {
4032 unsigned int components = 0;
4033 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4034 components += n_components_block(
4035 *vec[comp], std::bool_constant<IsBlockVector<VectorStruct>::value>());
4036 return components;
4037 }
4038
4039
4040
4041 // A helper function to identify block vectors with many components where we
4042 // should not try to overlap computations and communication because there
4043 // would be too many outstanding communication requests.
4044
4045 // default value for vectors that do not have communication_block_size
4046 template <typename VectorStruct,
4047 std::enable_if_t<!has_communication_block_size<VectorStruct>,
4048 VectorStruct> * = nullptr>
4049 constexpr unsigned int
4050 get_communication_block_size(const VectorStruct &)
4051 {
4052 return numbers::invalid_unsigned_int;
4053 }
4054
4055
4056
4057 template <typename VectorStruct,
4058 std::enable_if_t<has_communication_block_size<VectorStruct>,
4059 VectorStruct> * = nullptr>
4060 constexpr unsigned int
4061 get_communication_block_size(const VectorStruct &)
4062 {
4063 return VectorStruct::communication_block_size;
4064 }
4065
4066
4067
4068 template <typename VectorType,
4069 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType> * =
4070 nullptr>
4071 bool
4072 has_ghost_elements(const VectorType &vec)
4073 {
4074 (void)vec;
4075 return false;
4076 }
4077
4078
4079
4080 template <typename VectorType,
4081 std::enable_if_t<!is_not_parallel_vector<VectorType>, VectorType>
4082 * = nullptr>
4083 bool
4084 has_ghost_elements(const VectorType &vec)
4085 {
4086 return vec.has_ghost_elements();
4087 }
4088
4089
4090
4091 // --------------------------------------------------------------------------
4092 // below we have wrappers to distinguish between block and non-block vectors.
4093 // --------------------------------------------------------------------------
4094
4095 //
4096 // update_ghost_values_start
4097 //
4098
4099 // update_ghost_values for block vectors
4100 template <int dim,
4101 typename VectorStruct,
4102 typename Number,
4103 typename VectorizedArrayType,
4104 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4105 * = nullptr>
4106 void
4107 update_ghost_values_start(
4108 const VectorStruct &vec,
4109 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4110 const unsigned int channel = 0)
4111 {
4112 if (get_communication_block_size(vec) < vec.n_blocks())
4113 {
4114 const bool ghosts_set = vec.has_ghost_elements();
4115
4116 Assert(exchanger.matrix_free.get_task_info()
4117 .allow_ghosted_vectors_in_loops ||
4118 ghosts_set == false,
4119 ExcNotImplemented());
4120
4121 if (ghosts_set)
4122 {
4123 exchanger.ghosts_were_set = true;
4124 return;
4125 }
4126
4127 vec.update_ghost_values();
4128 }
4129 else
4130 {
4131 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4132 update_ghost_values_start(vec.block(i), exchanger, channel + i);
4133 }
4134 }
4135
4136
4137
4138 // update_ghost_values for non-block vectors
4139 template <int dim,
4140 typename VectorStruct,
4141 typename Number,
4142 typename VectorizedArrayType,
4143 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4144 * = nullptr>
4145 void
4146 update_ghost_values_start(
4147 const VectorStruct &vec,
4148 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4149 const unsigned int channel = 0)
4150 {
4151 exchanger.update_ghost_values_start(channel, vec);
4152 }
4153
4154
4155
4156 // update_ghost_values_start() for vector of vectors
4157 template <int dim,
4158 typename VectorStruct,
4159 typename Number,
4160 typename VectorizedArrayType>
4161 inline void
4162 update_ghost_values_start(
4163 const std::vector<VectorStruct> &vec,
4164 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4165 {
4166 unsigned int component_index = 0;
4167 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4168 {
4169 update_ghost_values_start(vec[comp], exchanger, component_index);
4170 component_index += n_components(vec[comp]);
4171 }
4172 }
4173
4174
4175
4176 // update_ghost_values_start() for vector of pointers to vectors
4177 template <int dim,
4178 typename VectorStruct,
4179 typename Number,
4180 typename VectorizedArrayType>
4181 inline void
4182 update_ghost_values_start(
4183 const std::vector<VectorStruct *> &vec,
4184 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4185 {
4186 unsigned int component_index = 0;
4187 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4188 {
4189 update_ghost_values_start(*vec[comp], exchanger, component_index);
4190 component_index += n_components(*vec[comp]);
4191 }
4192 }
4193
4194
4195
4196 //
4197 // update_ghost_values_finish
4198 //
4199
4200 // for block vectors
4201 template <int dim,
4202 typename VectorStruct,
4203 typename Number,
4204 typename VectorizedArrayType,
4205 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4206 * = nullptr>
4207 void
4208 update_ghost_values_finish(
4209 const VectorStruct &vec,
4210 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4211 const unsigned int channel = 0)
4212 {
4213 if (get_communication_block_size(vec) < vec.n_blocks())
4214 {
4215 // do nothing, everything has already been completed in the _start()
4216 // call
4217 }
4218 else
4219 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4220 update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4221 }
4222
4223
4224
4225 // for non-block vectors
4226 template <int dim,
4227 typename VectorStruct,
4228 typename Number,
4229 typename VectorizedArrayType,
4230 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4231 * = nullptr>
4232 void
4233 update_ghost_values_finish(
4234 const VectorStruct &vec,
4235 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4236 const unsigned int channel = 0)
4237 {
4238 exchanger.update_ghost_values_finish(channel, vec);
4239 }
4240
4241
4242
4243 // for vector of vectors
4244 template <int dim,
4245 typename VectorStruct,
4246 typename Number,
4247 typename VectorizedArrayType>
4248 inline void
4249 update_ghost_values_finish(
4250 const std::vector<VectorStruct> &vec,
4251 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4252 {
4253 unsigned int component_index = 0;
4254 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4255 {
4256 update_ghost_values_finish(vec[comp], exchanger, component_index);
4257 component_index += n_components(vec[comp]);
4258 }
4259 }
4260
4261
4262
4263 // for vector of pointers to vectors
4264 template <int dim,
4265 typename VectorStruct,
4266 typename Number,
4267 typename VectorizedArrayType>
4268 inline void
4269 update_ghost_values_finish(
4270 const std::vector<VectorStruct *> &vec,
4271 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4272 {
4273 unsigned int component_index = 0;
4274 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4275 {
4276 update_ghost_values_finish(*vec[comp], exchanger, component_index);
4277 component_index += n_components(*vec[comp]);
4278 }
4279 }
4280
4281
4282
4283 //
4284 // compress_start
4285 //
4286
4287 // for block vectors
4288 template <int dim,
4289 typename VectorStruct,
4290 typename Number,
4291 typename VectorizedArrayType,
4292 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4293 * = nullptr>
4294 inline void
4295 compress_start(
4296 VectorStruct &vec,
4297 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4298 const unsigned int channel = 0)
4299 {
4300 if (get_communication_block_size(vec) < vec.n_blocks())
4301 vec.compress(VectorOperation::add);
4302 else
4303 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4304 compress_start(vec.block(i), exchanger, channel + i);
4305 }
4306
4307
4308
4309 // for non-block vectors
4310 template <int dim,
4311 typename VectorStruct,
4312 typename Number,
4313 typename VectorizedArrayType,
4314 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4315 * = nullptr>
4316 inline void
4317 compress_start(
4318 VectorStruct &vec,
4319 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4320 const unsigned int channel = 0)
4321 {
4322 exchanger.compress_start(channel, vec);
4323 }
4324
4325
4326
4327 // for std::vector of vectors
4328 template <int dim,
4329 typename VectorStruct,
4330 typename Number,
4331 typename VectorizedArrayType>
4332 inline void
4333 compress_start(
4334 std::vector<VectorStruct> &vec,
4335 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4336 {
4337 unsigned int component_index = 0;
4338 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4339 {
4340 compress_start(vec[comp], exchanger, component_index);
4341 component_index += n_components(vec[comp]);
4342 }
4343 }
4344
4345
4346
4347 // for std::vector of pointer to vectors
4348 template <int dim,
4349 typename VectorStruct,
4350 typename Number,
4351 typename VectorizedArrayType>
4352 inline void
4353 compress_start(
4354 std::vector<VectorStruct *> &vec,
4355 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4356 {
4357 unsigned int component_index = 0;
4358 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4359 {
4360 compress_start(*vec[comp], exchanger, component_index);
4361 component_index += n_components(*vec[comp]);
4362 }
4363 }
4364
4365
4366
4367 //
4368 // compress_finish
4369 //
4370
4371 // for block vectors
4372 template <int dim,
4373 typename VectorStruct,
4374 typename Number,
4375 typename VectorizedArrayType,
4376 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4377 * = nullptr>
4378 inline void
4379 compress_finish(
4380 VectorStruct &vec,
4381 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4382 const unsigned int channel = 0)
4383 {
4384 if (get_communication_block_size(vec) < vec.n_blocks())
4385 {
4386 // do nothing, everything has already been completed in the _start()
4387 // call
4388 }
4389 else
4390 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4391 compress_finish(vec.block(i), exchanger, channel + i);
4392 }
4393
4394
4395
4396 // for non-block vectors
4397 template <int dim,
4398 typename VectorStruct,
4399 typename Number,
4400 typename VectorizedArrayType,
4401 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4402 * = nullptr>
4403 inline void
4404 compress_finish(
4405 VectorStruct &vec,
4406 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4407 const unsigned int channel = 0)
4408 {
4409 exchanger.compress_finish(channel, vec);
4410 }
4411
4412
4413
4414 // for std::vector of vectors
4415 template <int dim,
4416 typename VectorStruct,
4417 typename Number,
4418 typename VectorizedArrayType>
4419 inline void
4420 compress_finish(
4421 std::vector<VectorStruct> &vec,
4422 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4423 {
4424 unsigned int component_index = 0;
4425 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4426 {
4427 compress_finish(vec[comp], exchanger, component_index);
4428 component_index += n_components(vec[comp]);
4429 }
4430 }
4431
4432
4433
4434 // for std::vector of pointer to vectors
4435 template <int dim,
4436 typename VectorStruct,
4437 typename Number,
4438 typename VectorizedArrayType>
4439 inline void
4440 compress_finish(
4441 std::vector<VectorStruct *> &vec,
4442 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4443 {
4444 unsigned int component_index = 0;
4445 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4446 {
4447 compress_finish(*vec[comp], exchanger, component_index);
4448 component_index += n_components(*vec[comp]);
4449 }
4450 }
4451
4452
4453
4454 //
4455 // reset_ghost_values:
4456 //
4457 // if the input vector did not have ghosts imported, clear them here again
4458 // in order to avoid subsequent operations e.g. in linear solvers to work
4459 // with ghosts all the time
4460 //
4461
4462 // for block vectors
4463 template <int dim,
4464 typename VectorStruct,
4465 typename Number,
4466 typename VectorizedArrayType,
4467 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4468 * = nullptr>
4469 inline void
4470 reset_ghost_values(
4471 const VectorStruct &vec,
4472 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4473 {
4474 // return immediately if there is nothing to do.
4475 if (exchanger.ghosts_were_set == true)
4476 return;
4477
4478 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4479 reset_ghost_values(vec.block(i), exchanger);
4480 }
4481
4482
4483
4484 // for non-block vectors
4485 template <int dim,
4486 typename VectorStruct,
4487 typename Number,
4488 typename VectorizedArrayType,
4489 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4490 * = nullptr>
4491 inline void
4492 reset_ghost_values(
4493 const VectorStruct &vec,
4494 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4495 {
4496 // return immediately if there is nothing to do.
4497 if (exchanger.ghosts_were_set == true)
4498 return;
4499
4500 exchanger.reset_ghost_values(vec);
4501 }
4502
4503
4504
4505 // for std::vector of vectors
4506 template <int dim,
4507 typename VectorStruct,
4508 typename Number,
4509 typename VectorizedArrayType>
4510 inline void
4511 reset_ghost_values(
4512 const std::vector<VectorStruct> &vec,
4513 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4514 {
4515 // return immediately if there is nothing to do.
4516 if (exchanger.ghosts_were_set == true)
4517 return;
4518
4519 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4520 reset_ghost_values(vec[comp], exchanger);
4521 }
4522
4523
4524
4525 // for std::vector of pointer to vectors
4526 template <int dim,
4527 typename VectorStruct,
4528 typename Number,
4529 typename VectorizedArrayType>
4530 inline void
4531 reset_ghost_values(
4532 const std::vector<VectorStruct *> &vec,
4533 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4534 {
4535 // return immediately if there is nothing to do.
4536 if (exchanger.ghosts_were_set == true)
4537 return;
4538
4539 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4540 reset_ghost_values(*vec[comp], exchanger);
4541 }
4542
4543
4544
4545 //
4546 // zero_vector_region
4547 //
4548
4549 // for block vectors
4550 template <int dim,
4551 typename VectorStruct,
4552 typename Number,
4553 typename VectorizedArrayType,
4554 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4555 * = nullptr>
4556 inline void
4557 zero_vector_region(
4558 const unsigned int range_index,
4559 VectorStruct &vec,
4560 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4561 {
4562 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4563 exchanger.zero_vector_region(range_index, vec.block(i));
4564 }
4565
4566
4567
4568 // for non-block vectors
4569 template <int dim,
4570 typename VectorStruct,
4571 typename Number,
4572 typename VectorizedArrayType,
4573 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4574 * = nullptr>
4575 inline void
4576 zero_vector_region(
4577 const unsigned int range_index,
4578 VectorStruct &vec,
4579 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4580 {
4581 exchanger.zero_vector_region(range_index, vec);
4582 }
4583
4584
4585
4586 // for std::vector of vectors
4587 template <int dim,
4588 typename VectorStruct,
4589 typename Number,
4590 typename VectorizedArrayType>
4591 inline void
4592 zero_vector_region(
4593 const unsigned int range_index,
4594 std::vector<VectorStruct> &vec,
4595 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4596 {
4597 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4598 zero_vector_region(range_index, vec[comp], exchanger);
4599 }
4600
4601
4602
4603 // for std::vector of pointers to vectors
4604 template <int dim,
4605 typename VectorStruct,
4606 typename Number,
4607 typename VectorizedArrayType>
4608 inline void
4609 zero_vector_region(
4610 const unsigned int range_index,
4611 std::vector<VectorStruct *> &vec,
4612 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4613 {
4614 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4615 zero_vector_region(range_index, *vec[comp], exchanger);
4616 }
4617
4618
4619
4620 // Apply a unit matrix operation to constrained DoFs: Default cases where we
4621 // cannot detect a LinearAlgebra::distributed::Vector, we do not do
4622 // anything, else we apply the constraints as a unit operation
4623 template <typename VectorStruct1, typename VectorStruct2>
4624 inline void
4625 apply_operation_to_constrained_dofs(const std::vector<unsigned int> &,
4626 const VectorStruct1 &,
4627 VectorStruct2 &)
4628 {}
4629
4630 template <typename Number>
4631 inline void
4632 apply_operation_to_constrained_dofs(
4633 const std::vector<unsigned int> &constrained_dofs,
4634 const LinearAlgebra::distributed::Vector<Number> &src,
4635 LinearAlgebra::distributed::Vector<Number> &dst)
4636 {
4637 for (const unsigned int i : constrained_dofs)
4638 dst.local_element(i) = src.local_element(i);
4639 }
4640
4641
4642 namespace MatrixFreeFunctions
4643 {
4644 // struct to select between a const interface and a non-const interface
4645 // for MFWorker
4646 template <typename, typename, typename, typename, bool>
4647 struct InterfaceSelector
4648 {};
4649
4650 // Version for constant functions
4651 template <typename MF,
4652 typename InVector,
4653 typename OutVector,
4654 typename Container>
4655 struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4656 {
4657 using function_type = void (Container::*)(
4658 const MF &,
4659 OutVector &,
4660 const InVector &,
4661 const std::pair<unsigned int, unsigned int> &) const;
4662 };
4663
4664 // Version for non-constant functions
4665 template <typename MF,
4666 typename InVector,
4667 typename OutVector,
4668 typename Container>
4669 struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4670 {
4671 using function_type =
4672 void (Container::*)(const MF &,
4673 OutVector &,
4674 const InVector &,
4675 const std::pair<unsigned int, unsigned int> &);
4676 };
4677 } // namespace MatrixFreeFunctions
4678
4679
4680
4681 // A implementation class for the worker object that runs the various
4682 // operations we want to perform during the matrix-free loop
4683 template <typename MF,
4684 typename InVector,
4685 typename OutVector,
4686 typename Container,
4687 bool is_constant>
4688 class MFWorker : public MFWorkerInterface
4689 {
4690 public:
4691 // An alias to make the arguments further down more readable
4692 using function_type = typename MatrixFreeFunctions::
4693 InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4694 function_type;
4695
4696 // constructor, binds all the arguments to this class
4697 MFWorker(const MF &matrix_free,
4698 const InVector &src,
4699 OutVector &dst,
4700 const bool zero_dst_vector_setting,
4701 const Container &container,
4702 function_type cell_function,
4703 function_type face_function,
4704 function_type boundary_function,
4705 const typename MF::DataAccessOnFaces src_vector_face_access =
4706 MF::DataAccessOnFaces::none,
4707 const typename MF::DataAccessOnFaces dst_vector_face_access =
4708 MF::DataAccessOnFaces::none,
4709 const std::function<void(const unsigned int, const unsigned int)>
4710 &operation_before_loop = {},
4711 const std::function<void(const unsigned int, const unsigned int)>
4712 &operation_after_loop = {},
4713 const unsigned int dof_handler_index_pre_post = 0)
4714 : matrix_free(matrix_free)
4715 , container(const_cast<Container &>(container))
4716 , cell_function(cell_function)
4717 , face_function(face_function)
4718 , boundary_function(boundary_function)
4719 , src(src)
4720 , dst(dst)
4721 , src_data_exchanger(matrix_free,
4722 src_vector_face_access,
4723 n_components(src))
4724 , dst_data_exchanger(matrix_free,
4725 dst_vector_face_access,
4726 n_components(dst))
4727 , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4728 , zero_dst_vector_setting(zero_dst_vector_setting &&
4729 !src_and_dst_are_same)
4730 , operation_before_loop(operation_before_loop)
4731 , operation_after_loop(operation_after_loop)
4732 , dof_handler_index_pre_post(dof_handler_index_pre_post)
4733 {
4734 Assert(!has_ghost_elements(dst),
4735 ExcMessage("The destination vector passed to the matrix-free "
4736 "loop is ghosted. This is not allowed."));
4737 }
4738
4739 // Runs the cell work. If no function is given, nothing is done
4740 virtual void
4741 cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4742 {
4743 if (cell_function != nullptr && cell_range.second > cell_range.first)
4744 for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4745 {
4746 const auto cell_subrange =
4747 matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4748
4749 if (cell_subrange.second <= cell_subrange.first)
4750 continue;
4751
4752 (container.*
4753 cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4754 }
4755 }
4756
4757 virtual void
4758 cell(const unsigned int range_index) override
4759 {
4760 process_range(cell_function,
4761 matrix_free.get_task_info().cell_partition_data_hp_ptr,
4762 matrix_free.get_task_info().cell_partition_data_hp,
4763 range_index);
4764 }
4765
4766 virtual void
4767 face(const unsigned int range_index) override
4768 {
4769 process_range(face_function,
4770 matrix_free.get_task_info().face_partition_data_hp_ptr,
4771 matrix_free.get_task_info().face_partition_data_hp,
4772 range_index);
4773 }
4774
4775 virtual void
4776 boundary(const unsigned int range_index) override
4777 {
4778 process_range(boundary_function,
4779 matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4780 matrix_free.get_task_info().boundary_partition_data_hp,
4781 range_index);
4782 }
4783
4784 private:
4785 void
4786 process_range(const function_type &fu,
4787 const std::vector<unsigned int> &ptr,
4788 const std::vector<unsigned int> &data,
4789 const unsigned int range_index)
4790 {
4791 if (fu == nullptr)
4792 return;
4793
4794 AssertIndexRange(range_index + 1, ptr.size());
4795 for (unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4796 {
4797 AssertIndexRange(2 * i + 1, data.size());
4798 (container.*fu)(matrix_free,
4799 this->dst,
4800 this->src,
4801 std::make_pair(data[2 * i], data[2 * i + 1]));
4802 }
4803 }
4804
4805 public:
4806 // Starts the communication for the update ghost values operation. We
4807 // cannot call this update if ghost and destination are the same because
4808 // that would introduce spurious entries in the destination (there is also
4809 // the problem that reading from a vector that we also write to is usually
4810 // not intended in case there is overlap, but this is up to the
4811 // application code to decide and we cannot catch this case here).
4812 virtual void
4813 vector_update_ghosts_start() override
4814 {
4815 if (!src_and_dst_are_same)
4816 internal::update_ghost_values_start(src, src_data_exchanger);
4817 }
4818
4819 // Finishes the communication for the update ghost values operation
4820 virtual void
4821 vector_update_ghosts_finish() override
4822 {
4823 if (!src_and_dst_are_same)
4824 internal::update_ghost_values_finish(src, src_data_exchanger);
4825 }
4826
4827 // Starts the communication for the vector compress operation
4828 virtual void
4829 vector_compress_start() override
4830 {
4831 internal::compress_start(dst, dst_data_exchanger);
4832 }
4833
4834 // Finishes the communication for the vector compress operation
4835 virtual void
4836 vector_compress_finish() override
4837 {
4838 internal::compress_finish(dst, dst_data_exchanger);
4839 if (!src_and_dst_are_same)
4840 internal::reset_ghost_values(src, src_data_exchanger);
4841 }
4842
4843 // Zeros the given input vector
4844 virtual void
4845 zero_dst_vector_range(const unsigned int range_index) override
4846 {
4847 if (zero_dst_vector_setting)
4848 internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4849 }
4850
4851 virtual void
4852 cell_loop_pre_range(const unsigned int range_index) override
4853 {
4854 if (operation_before_loop)
4855 {
4856 const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4857 matrix_free.get_dof_info(dof_handler_index_pre_post);
4858 if (range_index == numbers::invalid_unsigned_int)
4859 {
4860 // Case with threaded loop -> currently no overlap implemented
4861 ::parallel::apply_to_subranges(
4862 0U,
4863 dof_info.vector_partitioner->locally_owned_size(),
4864 operation_before_loop,
4865 internal::VectorImplementation::minimum_parallel_grain_size);
4866 }
4867 else
4868 {
4869 AssertIndexRange(range_index,
4870 dof_info.cell_loop_pre_list_index.size() - 1);
4871 for (unsigned int id =
4872 dof_info.cell_loop_pre_list_index[range_index];
4873 id != dof_info.cell_loop_pre_list_index[range_index + 1];
4874 ++id)
4875 operation_before_loop(dof_info.cell_loop_pre_list[id].first,
4876 dof_info.cell_loop_pre_list[id].second);
4877 }
4878 }
4879 }
4880
4881 virtual void
4882 cell_loop_post_range(const unsigned int range_index) override
4883 {
4884 if (operation_after_loop)
4885 {
4886 // Run unit matrix operation on constrained dofs if we are at the
4887 // last range
4888 const std::vector<unsigned int> &partition_row_index =
4889 matrix_free.get_task_info().partition_row_index;
4890 if (range_index ==
4891 partition_row_index[partition_row_index.size() - 2] - 1)
4892 apply_operation_to_constrained_dofs(
4893 matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
4894 src,
4895 dst);
4896
4897 const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4898 matrix_free.get_dof_info(dof_handler_index_pre_post);
4899 if (range_index == numbers::invalid_unsigned_int)
4900 {
4901 // Case with threaded loop -> currently no overlap implemented
4902 ::parallel::apply_to_subranges(
4903 0U,
4904 dof_info.vector_partitioner->locally_owned_size(),
4905 operation_after_loop,
4906 internal::VectorImplementation::minimum_parallel_grain_size);
4907 }
4908 else
4909 {
4910 AssertIndexRange(range_index,
4911 dof_info.cell_loop_post_list_index.size() - 1);
4912 for (unsigned int id =
4913 dof_info.cell_loop_post_list_index[range_index];
4914 id != dof_info.cell_loop_post_list_index[range_index + 1];
4915 ++id)
4916 operation_after_loop(dof_info.cell_loop_post_list[id].first,
4917 dof_info.cell_loop_post_list[id].second);
4918 }
4919 }
4920 }
4921
4922 private:
4923 const MF &matrix_free;
4924 Container &container;
4925 function_type cell_function;
4926 function_type face_function;
4927 function_type boundary_function;
4928
4929 const InVector &src;
4930 OutVector &dst;
4931 VectorDataExchange<MF::dimension,
4932 typename MF::value_type,
4933 typename MF::vectorized_value_type>
4934 src_data_exchanger;
4935 VectorDataExchange<MF::dimension,
4936 typename MF::value_type,
4937 typename MF::vectorized_value_type>
4938 dst_data_exchanger;
4939 const bool src_and_dst_are_same;
4940 const bool zero_dst_vector_setting;
4941 const std::function<void(const unsigned int, const unsigned int)>
4942 operation_before_loop;
4943 const std::function<void(const unsigned int, const unsigned int)>
4944 operation_after_loop;
4945 const unsigned int dof_handler_index_pre_post;
4946 };
4947
4948
4949
4954 template <class MF, typename InVector, typename OutVector>
4955 struct MFClassWrapper
4956 {
4957 using function_type =
4958 std::function<void(const MF &,
4959 OutVector &,
4960 const InVector &,
4961 const std::pair<unsigned int, unsigned int> &)>;
4962
4963 MFClassWrapper(const function_type cell,
4964 const function_type face,
4965 const function_type boundary)
4966 : cell(cell)
4967 , face(face)
4968 , boundary(boundary)
4969 {}
4970
4971 void
4972 cell_integrator(const MF &mf,
4973 OutVector &dst,
4974 const InVector &src,
4975 const std::pair<unsigned int, unsigned int> &range) const
4976 {
4977 if (cell)
4978 cell(mf, dst, src, range);
4979 }
4980
4981 void
4982 face_integrator(const MF &mf,
4983 OutVector &dst,
4984 const InVector &src,
4985 const std::pair<unsigned int, unsigned int> &range) const
4986 {
4987 if (face)
4988 face(mf, dst, src, range);
4989 }
4990
4991 void
4992 boundary_integrator(
4993 const MF &mf,
4994 OutVector &dst,
4995 const InVector &src,
4996 const std::pair<unsigned int, unsigned int> &range) const
4997 {
4998 if (boundary)
4999 boundary(mf, dst, src, range);
5000 }
5001
5002 const function_type cell;
5003 const function_type face;
5004 const function_type boundary;
5005 };
5006
5007} // end of namespace internal
5008
5009
5010
5011template <int dim, typename Number, typename VectorizedArrayType>
5012template <typename OutVector, typename InVector>
5013inline void
5014MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5015 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5016 OutVector &,
5017 const InVector &,
5018 const std::pair<unsigned int, unsigned int> &)>
5019 &cell_operation,
5020 OutVector &dst,
5021 const InVector &src,
5022 const bool zero_dst_vector) const
5023{
5024 using Wrapper =
5025 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5026 InVector,
5027 OutVector>;
5028 Wrapper wrap(cell_operation, nullptr, nullptr);
5029 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5030 InVector,
5031 OutVector,
5032 Wrapper,
5033 true>
5034 worker(*this,
5035 src,
5036 dst,
5037 zero_dst_vector,
5038 wrap,
5039 &Wrapper::cell_integrator,
5040 &Wrapper::face_integrator,
5041 &Wrapper::boundary_integrator);
5042
5043 task_info.loop(worker);
5044}
5045
5046
5047
5048template <int dim, typename Number, typename VectorizedArrayType>
5049template <typename OutVector, typename InVector>
5050inline void
5051MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5052 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5053 OutVector &,
5054 const InVector &,
5055 const std::pair<unsigned int, unsigned int> &)>
5056 &cell_operation,
5057 OutVector &dst,
5058 const InVector &src,
5059 const std::function<void(const unsigned int, const unsigned int)>
5060 &operation_before_loop,
5061 const std::function<void(const unsigned int, const unsigned int)>
5062 &operation_after_loop,
5063 const unsigned int dof_handler_index_pre_post) const
5064{
5065 using Wrapper =
5066 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5067 InVector,
5068 OutVector>;
5069 Wrapper wrap(cell_operation, nullptr, nullptr);
5070 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5071 InVector,
5072 OutVector,
5073 Wrapper,
5074 true>
5075 worker(*this,
5076 src,
5077 dst,
5078 false,
5079 wrap,
5080 &Wrapper::cell_integrator,
5081 &Wrapper::face_integrator,
5082 &Wrapper::boundary_integrator,
5083 DataAccessOnFaces::none,
5084 DataAccessOnFaces::none,
5085 operation_before_loop,
5086 operation_after_loop,
5087 dof_handler_index_pre_post);
5088
5089 task_info.loop(worker);
5090}
5091
5092
5093
5094template <int dim, typename Number, typename VectorizedArrayType>
5095template <typename OutVector, typename InVector>
5096inline void
5097MatrixFree<dim, Number, VectorizedArrayType>::loop(
5098 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5099 OutVector &,
5100 const InVector &,
5101 const std::pair<unsigned int, unsigned int> &)>
5102 &cell_operation,
5103 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5104 OutVector &,
5105 const InVector &,
5106 const std::pair<unsigned int, unsigned int> &)>
5107 &inner_face_operation,
5108 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5109 OutVector &,
5110 const InVector &,
5111 const std::pair<unsigned int, unsigned int> &)>
5112 &boundary_face_operation,
5113 OutVector &dst,
5114 const InVector &src,
5115 const bool zero_dst_vector,
5116 const DataAccessOnFaces dst_vector_face_access,
5117 const DataAccessOnFaces src_vector_face_access) const
5118{
5119 using Wrapper =
5120 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5121 InVector,
5122 OutVector>;
5123 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5124 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5125 InVector,
5126 OutVector,
5127 Wrapper,
5128 true>
5129 worker(*this,
5130 src,
5131 dst,
5132 zero_dst_vector,
5133 wrap,
5134 &Wrapper::cell_integrator,
5135 &Wrapper::face_integrator,
5136 &Wrapper::boundary_integrator,
5137 src_vector_face_access,
5138 dst_vector_face_access);
5139
5140 task_info.loop(worker);
5141}
5142
5143
5144
5145template <int dim, typename Number, typename VectorizedArrayType>
5146template <typename CLASS, typename OutVector, typename InVector>
5147inline void
5148MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5149 void (CLASS::*function_pointer)(
5150 const MatrixFree<dim, Number, VectorizedArrayType> &,
5151 OutVector &,
5152 const InVector &,
5153 const std::pair<unsigned int, unsigned int> &) const,
5154 const CLASS *owning_class,
5155 OutVector &dst,
5156 const InVector &src,
5157 const bool zero_dst_vector) const
5158{
5159 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5160 InVector,
5161 OutVector,
5162 CLASS,
5163 true>
5164 worker(*this,
5165 src,
5166 dst,
5167 zero_dst_vector,
5168 *owning_class,
5169 function_pointer,
5170 nullptr,
5171 nullptr);
5172 task_info.loop(worker);
5173}
5174
5175
5176
5177template <int dim, typename Number, typename VectorizedArrayType>
5178template <typename CLASS, typename OutVector, typename InVector>
5179inline void
5180MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5181 void (CLASS::*function_pointer)(
5182 const MatrixFree<dim, Number, VectorizedArrayType> &,
5183 OutVector &,
5184 const InVector &,
5185 const std::pair<unsigned int, unsigned int> &) const,
5186 const CLASS *owning_class,
5187 OutVector &dst,
5188 const InVector &src,
5189 const std::function<void(const unsigned int, const unsigned int)>
5190 &operation_before_loop,
5191 const std::function<void(const unsigned int, const unsigned int)>
5192 &operation_after_loop,
5193 const unsigned int dof_handler_index_pre_post) const
5194{
5195 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5196 InVector,
5197 OutVector,
5198 CLASS,
5199 true>
5200 worker(*this,
5201 src,
5202 dst,
5203 false,
5204 *owning_class,
5205 function_pointer,
5206 nullptr,
5207 nullptr,
5208 DataAccessOnFaces::none,
5209 DataAccessOnFaces::none,
5210 operation_before_loop,
5211 operation_after_loop,
5212 dof_handler_index_pre_post);
5213 task_info.loop(worker);
5214}
5215
5216
5217
5218template <int dim, typename Number, typename VectorizedArrayType>
5219template <typename CLASS, typename OutVector, typename InVector>
5220inline void
5221MatrixFree<dim, Number, VectorizedArrayType>::loop(
5222 void (CLASS::*cell_operation)(
5223 const MatrixFree<dim, Number, VectorizedArrayType> &,
5224 OutVector &,
5225 const InVector &,
5226 const std::pair<unsigned int, unsigned int> &) const,
5227 void (CLASS::*inner_face_operation)(
5228 const MatrixFree<dim, Number, VectorizedArrayType> &,
5229 OutVector &,
5230 const InVector &,
5231 const std::pair<unsigned int, unsigned int> &) const,
5232 void (CLASS::*boundary_face_operation)(
5233 const MatrixFree<dim, Number, VectorizedArrayType> &,
5234 OutVector &,
5235 const InVector &,
5236 const std::pair<unsigned int, unsigned int> &) const,
5237 const CLASS *owning_class,
5238 OutVector &dst,
5239 const InVector &src,
5240 const bool zero_dst_vector,
5241 const DataAccessOnFaces dst_vector_face_access,
5242 const DataAccessOnFaces src_vector_face_access) const
5243{
5244 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5245 InVector,
5246 OutVector,
5247 CLASS,
5248 true>
5249 worker(*this,
5250 src,
5251 dst,
5252 zero_dst_vector,
5253 *owning_class,
5254 cell_operation,
5255 inner_face_operation,
5256 boundary_face_operation,
5257 src_vector_face_access,
5258 dst_vector_face_access);
5259 task_info.loop(worker);
5260}
5261
5262
5263
5264template <int dim, typename Number, typename VectorizedArrayType>
5265template <typename CLASS, typename OutVector, typename InVector>
5266inline void
5267MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5268 void (CLASS::*function_pointer)(
5269 const MatrixFree<dim, Number, VectorizedArrayType> &,
5270 OutVector &,
5271 const InVector &,
5272 const std::pair<unsigned int, unsigned int> &),
5273 CLASS *owning_class,
5274 OutVector &dst,
5275 const InVector &src,
5276 const bool zero_dst_vector) const
5277{
5278 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5279 InVector,
5280 OutVector,
5281 CLASS,
5282 false>
5283 worker(*this,
5284 src,
5285 dst,
5286 zero_dst_vector,
5287 *owning_class,
5288 function_pointer,
5289 nullptr,
5290 nullptr);
5291 task_info.loop(worker);
5292}
5293
5294
5295
5296template <int dim, typename Number, typename VectorizedArrayType>
5297template <typename CLASS, typename OutVector, typename InVector>
5298inline void
5299MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5300 void (CLASS::*function_pointer)(
5301 const MatrixFree<dim, Number, VectorizedArrayType> &,
5302 OutVector &,
5303 const InVector &,
5304 const std::pair<unsigned int, unsigned int> &),
5305 CLASS *owning_class,
5306 OutVector &dst,
5307 const InVector &src,
5308 const std::function<void(const unsigned int, const unsigned int)>
5309 &operation_before_loop,
5310 const std::function<void(const unsigned int, const unsigned int)>
5311 &operation_after_loop,
5312 const unsigned int dof_handler_index_pre_post) const
5313{
5314 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5315 InVector,
5316 OutVector,
5317 CLASS,
5318 false>
5319 worker(*this,
5320 src,
5321 dst,
5322 false,
5323 *owning_class,
5324 function_pointer,
5325 nullptr,
5326 nullptr,
5327 DataAccessOnFaces::none,
5328 DataAccessOnFaces::none,
5329 operation_before_loop,
5330 operation_after_loop,
5331 dof_handler_index_pre_post);
5332 task_info.loop(worker);
5333}
5334
5335
5336
5337template <int dim, typename Number, typename VectorizedArrayType>
5338template <typename CLASS, typename OutVector, typename InVector>
5339inline void
5340MatrixFree<dim, Number, VectorizedArrayType>::loop(
5341 void (CLASS::*cell_operation)(
5342 const MatrixFree<dim, Number, VectorizedArrayType> &,
5343 OutVector &,
5344 const InVector &,
5345 const std::pair<unsigned int, unsigned int> &),
5346 void (CLASS::*inner_face_operation)(
5347 const MatrixFree<dim, Number, VectorizedArrayType> &,
5348 OutVector &,
5349 const InVector &,
5350 const std::pair<unsigned int, unsigned int> &),
5351 void (CLASS::*boundary_face_operation)(
5352 const MatrixFree<dim, Number, VectorizedArrayType> &,
5353 OutVector &,
5354 const InVector &,
5355 const std::pair<unsigned int, unsigned int> &),
5356 CLASS *owning_class,
5357 OutVector &dst,
5358 const InVector &src,
5359 const bool zero_dst_vector,
5360 const DataAccessOnFaces dst_vector_face_access,
5361 const DataAccessOnFaces src_vector_face_access) const
5362{
5363 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5364 InVector,
5365 OutVector,
5366 CLASS,
5367 false>
5368 worker(*this,
5369 src,
5370 dst,
5371 zero_dst_vector,
5372 *owning_class,
5373 cell_operation,
5374 inner_face_operation,
5375 boundary_face_operation,
5376 src_vector_face_access,
5377 dst_vector_face_access);
5378 task_info.loop(worker);
5379}
5380
5381
5382
5383template <int dim, typename Number, typename VectorizedArrayType>
5384template <typename OutVector, typename InVector>
5385inline void
5386MatrixFree<dim, Number, VectorizedArrayType>::loop(
5387 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5388 OutVector &,
5389 const InVector &,
5390 const std::pair<unsigned int, unsigned int> &)>
5391 &cell_operation,
5392 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5393 OutVector &,
5394 const InVector &,
5395 const std::pair<unsigned int, unsigned int> &)>
5396 &inner_face_operation,
5397 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5398 OutVector &,
5399 const InVector &,
5400 const std::pair<unsigned int, unsigned int> &)>
5401 &boundary_face_operation,
5402 OutVector &dst,
5403 const InVector &src,
5404 const std::function<void(const unsigned int, const unsigned int)>
5405 &operation_before_loop,
5406 const std::function<void(const unsigned int, const unsigned int)>
5407 &operation_after_loop,
5408 const unsigned int dof_handler_index_pre_post,
5409 const DataAccessOnFaces dst_vector_face_access,
5410 const DataAccessOnFaces src_vector_face_access) const
5411{
5412 using Wrapper =
5413 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5414 InVector,
5415 OutVector>;
5416 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5417 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5418 InVector,
5419 OutVector,
5420 Wrapper,
5421 true>
5422 worker(*this,
5423 src,
5424 dst,
5425 false,
5426 wrap,
5427 &Wrapper::cell_integrator,
5428 &Wrapper::face_integrator,
5429 &Wrapper::boundary_integrator,
5430 src_vector_face_access,
5431 dst_vector_face_access,
5432 operation_before_loop,
5433 operation_after_loop,
5434 dof_handler_index_pre_post);
5435
5436 task_info.loop(worker);
5437}
5438
5439
5440
5441template <int dim, typename Number, typename VectorizedArrayType>
5442template <typename CLASS, typename OutVector, typename InVector>
5443inline void
5444MatrixFree<dim, Number, VectorizedArrayType>::loop(
5445 void (CLASS::*cell_operation)(const MatrixFree &,
5446 OutVector &,
5447 const InVector &,
5448 const std::pair<unsigned int, unsigned int> &)
5449 const,
5450 void (CLASS::*inner_face_operation)(
5451 const MatrixFree &,
5452 OutVector &,
5453 const InVector &,
5454 const std::pair<unsigned int, unsigned int> &) const,
5455 void (CLASS::*boundary_face_operation)(
5456 const MatrixFree &,
5457 OutVector &,
5458 const InVector &,
5459 const std::pair<unsigned int, unsigned int> &) const,
5460 const CLASS *owning_class,
5461 OutVector &dst,
5462 const InVector &src,
5463 const std::function<void(const unsigned int, const unsigned int)>
5464 &operation_before_loop,
5465 const std::function<void(const unsigned int, const unsigned int)>
5466 &operation_after_loop,
5467 const unsigned int dof_handler_index_pre_post,
5468 const DataAccessOnFaces dst_vector_face_access,
5469 const DataAccessOnFaces src_vector_face_access) const
5470{
5471 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5472 InVector,
5473 OutVector,
5474 CLASS,
5475 true>
5476 worker(*this,
5477 src,
5478 dst,
5479 false,
5480 *owning_class,
5481 cell_operation,
5482 inner_face_operation,
5483 boundary_face_operation,
5484 src_vector_face_access,
5485 dst_vector_face_access,
5486 operation_before_loop,
5487 operation_after_loop,
5488 dof_handler_index_pre_post);
5489 task_info.loop(worker);
5490}
5491
5492
5493
5494template <int dim, typename Number, typename VectorizedArrayType>
5495template <typename CLASS, typename OutVector, typename InVector>
5496inline void
5497MatrixFree<dim, Number, VectorizedArrayType>::loop(
5498 void (CLASS::*cell_operation)(const MatrixFree &,
5499 OutVector &,
5500 const InVector &,
5501 const std::pair<unsigned int, unsigned int> &),
5502 void (CLASS::*inner_face_operation)(
5503 const MatrixFree &,
5504 OutVector &,
5505 const InVector &,
5506 const std::pair<unsigned int, unsigned int> &),
5507 void (CLASS::*boundary_face_operation)(
5508 const MatrixFree &,
5509 OutVector &,
5510 const InVector &,
5511 const std::pair<unsigned int, unsigned int> &),
5512 const CLASS *owning_class,
5513 OutVector &dst,
5514 const InVector &src,
5515 const std::function<void(const unsigned int, const unsigned int)>
5516 &operation_before_loop,
5517 const std::function<void(const unsigned int, const unsigned int)>
5518 &operation_after_loop,
5519 const unsigned int dof_handler_index_pre_post,
5520 const DataAccessOnFaces dst_vector_face_access,
5521 const DataAccessOnFaces src_vector_face_access) const
5522{
5523 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5524 InVector,
5525 OutVector,
5526 CLASS,
5527 false>
5528 worker(*this,
5529 src,
5530 dst,
5531 false,
5532 *owning_class,
5533 cell_operation,
5534 inner_face_operation,
5535 boundary_face_operation,
5536 src_vector_face_access,
5537 dst_vector_face_access,
5538 operation_before_loop,
5539 operation_after_loop,
5540 dof_handler_index_pre_post);
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> &) const,
5555 const 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 true>
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 CLASS, typename OutVector, typename InVector>
5589inline void
5590MatrixFree<dim, Number, VectorizedArrayType>::loop_cell_centric(
5591 void (CLASS::*function_pointer)(
5592 const MatrixFree<dim, Number, VectorizedArrayType> &,
5593 OutVector &,
5594 const InVector &,
5595 const std::pair<unsigned int, unsigned int> &),
5596 CLASS *owning_class,
5597 OutVector &dst,
5598 const InVector &src,
5599 const bool zero_dst_vector,
5600 const DataAccessOnFaces src_vector_face_access) const
5601{
5602 auto src_vector_face_access_temp = src_vector_face_access;
5603 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5604 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5605 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5606 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5607
5608 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5609 InVector,
5610 OutVector,
5611 CLASS,
5612 false>
5613 worker(*this,
5614 src,
5615 dst,
5616 zero_dst_vector,
5617 *owning_class,
5618 function_pointer,
5619 nullptr,
5620 nullptr,
5621 src_vector_face_access_temp,
5622 DataAccessOnFaces::none);
5623 task_info.loop(worker);
5624}
5625
5626
5627
5628template <int dim, typename Number, typename VectorizedArrayType>
5629template <typename OutVector, typename InVector>
5630inline void
5631MatrixFree<dim, Number, VectorizedArrayType>::loop_cell_centric(
5632 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5633 OutVector &,
5634 const InVector &,
5635 const std::pair<unsigned int, unsigned int> &)>
5636 &cell_operation,
5637 OutVector &dst,
5638 const InVector &src,
5639 const bool zero_dst_vector,
5640 const DataAccessOnFaces src_vector_face_access) const
5641{
5642 auto src_vector_face_access_temp = src_vector_face_access;
5643 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5644 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5645 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5646 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5647
5648 using Wrapper =
5649 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5650 InVector,
5651 OutVector>;
5652 Wrapper wrap(cell_operation, nullptr, nullptr);
5653
5654 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5655 InVector,
5656 OutVector,
5657 Wrapper,
5658 true>
5659 worker(*this,
5660 src,
5661 dst,
5662 zero_dst_vector,
5663 wrap,
5664 &Wrapper::cell_integrator,
5665 &Wrapper::face_integrator,
5666 &Wrapper::boundary_integrator,
5667 src_vector_face_access_temp,
5668 DataAccessOnFaces::none);
5669 task_info.loop(worker);
5670}
5671
5672
5673#endif // ifndef DOXYGEN
5674
5675
5676
5677DEAL_II_NAMESPACE_CLOSE
5678
5679#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:318
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:2291
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::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::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::loop
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*boundary_face_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::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::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 > &)> &inner_face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_face_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_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:2325
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::dof_handlers
std::vector< ObserverPointer< const DoFHandler< dim > > > dof_handlers
Definition matrix_free.h:2255
MatrixFree::loop
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*boundary_face_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::first_hp_dof_handler_index
unsigned int first_hp_dof_handler_index
Definition matrix_free.h:2371
MatrixFree::mapping_is_initialized
bool mapping_is_initialized
Definition matrix_free.h:2342
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:2365
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:2304
MatrixFree::get_vector_partitioner
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner(const unsigned int dof_handler_index=0) const
MatrixFree::affine_constraints
std::vector< ObserverPointer< const AffineConstraints< Number > > > affine_constraints
Definition matrix_free.h:2264
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:2278
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:2353
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::loop
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*boundary_face_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::scratch_pad_non_threadsafe
std::list< std::pair< bool, AlignedVector< Number > > > scratch_pad_non_threadsafe
Definition matrix_free.h:2360
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::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 > &)> &inner_face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_face_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::initialize_indices
void initialize_indices(const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const AdditionalData &additional_data)
MatrixFree::n_components
unsigned int n_components() const
MatrixFree::n_physical_cells
unsigned int n_physical_cells() const
MatrixFree::initialize_dof_vector
void initialize_dof_vector(VectorType &vec, const unsigned int dof_handler_index=0) const
MatrixFree::get_n_q_points
unsigned int get_n_q_points(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
MatrixFree::reinit
void reinit(const MappingType &mapping, const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
MatrixFree::mf_cell_indices
std::vector< unsigned int > mf_cell_indices
Definition matrix_free.h:2310
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::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*boundary_face_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::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:2270
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:2337
MatrixFree::face_info
internal::MatrixFreeFunctions::FaceInfo< VectorizedArrayType::size()> face_info
Definition matrix_free.h:2332
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::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:2284
MatrixFree::get_shape_info
const internal::MatrixFreeFunctions::ShapeInfo< Number > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
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:2318
MatrixFree::n_base_elements
unsigned int n_base_elements(const unsigned int dof_handler_index) const
MatrixFree::shape_info
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< Number > > shape_info
Definition matrix_free.h:2296
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:498
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
dof_handler.h
second
Point< 2 > second
Definition grid_out.cc:4630
first
Point< 2 > first
Definition grid_out.cc:4629
level
unsigned int level
Definition grid_out.cc:4632
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exceptions.h:1659
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:2026
AssertThrowMPI
#define AssertThrowMPI(error_code)
Definition exceptions.h:2174
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exceptions.h:2115
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:1780
DoFHandler::cell_iterator
typename ActiveSelector::cell_iterator cell_iterator
Definition dof_handler.h:468
UpdateFlags
UpdateFlags
Definition fe_update_flags.h:64
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition fe_update_flags.h:134
update_gradients
@ update_gradients
Shape function gradients.
Definition fe_update_flags.h:81
update_default
@ update_default
No update.
Definition fe_update_flags.h:66
la_parallel_vector.h
mapping_collection.h
dof_info.h
mapping_info.h
type_traits.h
data
std::vector< index_type > data
Definition mpi.cc:735
size
std::size_t size
Definition mpi.cc:734
LinearAlgebra::TpetraWrappers::TpetraTypes::VectorType
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
Definition trilinos_tpetra_types.h:112
TrilinosWrappers::internal::begin
VectorType::value_type * begin(VectorType &V)
Definition trilinos_sparse_matrix.cc:50
Utilities::MPI::job_supports_mpi
bool job_supports_mpi()
Definition mpi.cc:681
internal::VectorImplementation::minimum_parallel_grain_size
unsigned int minimum_parallel_grain_size
Definition parallel.cc:33
numbers::invalid_boundary_id
const types::boundary_id invalid_boundary_id
Definition types.h:292
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition types.h:220
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:452
std
STL namespace.
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6943
types::fe_index
unsigned short int fe_index
Definition types.h:59
types::boundary_id
unsigned int boundary_id
Definition types.h:144
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:749
internal::MatrixFreeFunctions::DoFInfo::cell_loop_post_list_index
std::vector< unsigned int > cell_loop_post_list_index
Definition dof_info.h:762
internal::MatrixFreeFunctions::DoFInfo::vector_zero_range_list
std::vector< std::pair< unsigned int, unsigned int > > vector_zero_range_list
Definition dof_info.h:742
internal::MatrixFreeFunctions::DoFInfo::vector_partitioner
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition dof_info.h:593
internal::MatrixFreeFunctions::DoFInfo::cell_loop_pre_list
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_pre_list
Definition dof_info.h:755
internal::MatrixFreeFunctions::DoFInfo::vector_zero_range_list_index
std::vector< unsigned int > vector_zero_range_list_index
Definition dof_info.h:737
internal::MatrixFreeFunctions::DoFInfo::cell_loop_post_list
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_post_list
Definition dof_info.h:768
internal::MatrixFreeFunctions::TaskInfo::loop
void loop(MFWorkerInterface &worker) const
Definition task_info.cc:350
task_info.h
template_constraints.h
thread_local_storage.h
vector_data_exchange.h
vector_operation.h
Generated by doxygen 1.9.8