deal.II version GIT relicensing-2165-gc91f007519 2024-11-20 01:40: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(ArrayView<Number>(
3881 const_cast<LinearAlgebra::distributed::Vector<Number> &>(vec)
3882 .begin() +
3883 part.locally_owned_size(),
3884 matrix_free.get_dof_info(mf_component)
3885 .vector_partitioner->n_ghost_indices()));
3886 }
3887
3888# endif
3889 }
3890 // let vector know that it's not ghosted anymore
3891 vec.set_ghost_state(false);
3892 }
3893
3894
3895
3901 template <typename VectorType,
3902 std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3903 * = nullptr>
3904 void
3905 zero_vector_region(const unsigned int range_index, VectorType &vec) const
3906 {
3907 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3908 "Type mismatch between VectorType and VectorDataExchange");
3909 if (range_index == numbers::invalid_unsigned_int)
3910 vec = Number();
3911 else
3912 {
3913 const unsigned int mf_component = find_vector_in_mf(vec, false);
3914 const internal::MatrixFreeFunctions::DoFInfo &dof_info =
3915 matrix_free.get_dof_info(mf_component);
3916 Assert(dof_info.vector_zero_range_list_index.empty() == false,
3917 ExcNotInitialized());
3918
3919 Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3920 ExcInternalError());
3921 AssertIndexRange(range_index,
3922 dof_info.vector_zero_range_list_index.size() - 1);
3923 for (unsigned int id =
3924 dof_info.vector_zero_range_list_index[range_index];
3925 id != dof_info.vector_zero_range_list_index[range_index + 1];
3926 ++id)
3927 std::memset(vec.begin() + dof_info.vector_zero_range_list[id].first,
3928 0,
3929 (dof_info.vector_zero_range_list[id].second -
3930 dof_info.vector_zero_range_list[id].first) *
3931 sizeof(Number));
3932 }
3933 }
3934
3935
3936
3942 template <typename VectorType,
3943 std::enable_if_t<!has_exchange_on_subset<VectorType>, VectorType>
3944 * = nullptr,
3945 typename VectorType::value_type * = nullptr>
3946 void
3947 zero_vector_region(const unsigned int range_index, VectorType &vec) const
3948 {
3949 if (range_index == numbers::invalid_unsigned_int || range_index == 0)
3950 {
3951 if constexpr (std::is_same_v<
3952 ArrayView<typename VectorType::value_type>,
3953 VectorType>)
3954 {
3955 for (unsigned int i = 0; i < vec.size(); ++i)
3956 vec[i] = typename VectorType::value_type();
3957 }
3958 else
3959 vec = typename VectorType::value_type();
3960 }
3961 }
3962
3963
3964
3969 void
3970 zero_vector_region(const unsigned int, ...) const
3971 {
3972 Assert(false,
3973 ExcNotImplemented("Zeroing is only implemented for vector types "
3974 "which provide VectorType::value_type"));
3975 }
3976
3977
3978
3979 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
3980 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3981 DataAccessOnFaces vector_face_access;
3982 bool ghosts_were_set;
3983# ifdef DEAL_II_WITH_MPI
3984 std::vector<AlignedVector<Number> *> tmp_data;
3985 std::vector<std::vector<MPI_Request>> requests;
3986# endif
3987 }; // VectorDataExchange
3988
3989 template <typename VectorStruct>
3990 unsigned int
3991 n_components(const VectorStruct &vec);
3992
3993 template <typename VectorStruct>
3994 unsigned int
3995 n_components_block(const VectorStruct &vec, const std::bool_constant<true>)
3996 {
3997 unsigned int components = 0;
3998 for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
3999 components += n_components(vec.block(bl));
4000 return components;
4001 }
4002
4003 template <typename VectorStruct>
4004 unsigned int
4005 n_components_block(const VectorStruct &, const std::bool_constant<false>)
4006 {
4007 return 1;
4008 }
4009
4010 template <typename VectorStruct>
4011 unsigned int
4012 n_components(const VectorStruct &vec)
4013 {
4014 return n_components_block(
4015 vec, std::bool_constant<IsBlockVector<VectorStruct>::value>());
4016 }
4017
4018 template <typename VectorStruct>
4019 inline unsigned int
4020 n_components(const std::vector<VectorStruct> &vec)
4021 {
4022 unsigned int components = 0;
4023 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4024 components += n_components_block(
4025 vec[comp], std::bool_constant<IsBlockVector<VectorStruct>::value>());
4026 return components;
4027 }
4028
4029 template <typename VectorStruct>
4030 inline unsigned int
4031 n_components(const std::vector<VectorStruct *> &vec)
4032 {
4033 unsigned int components = 0;
4034 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4035 components += n_components_block(
4036 *vec[comp], std::bool_constant<IsBlockVector<VectorStruct>::value>());
4037 return components;
4038 }
4039
4040
4041
4042 // A helper function to identify block vectors with many components where we
4043 // should not try to overlap computations and communication because there
4044 // would be too many outstanding communication requests.
4045
4046 // default value for vectors that do not have communication_block_size
4047 template <typename VectorStruct,
4048 std::enable_if_t<!has_communication_block_size<VectorStruct>,
4049 VectorStruct> * = nullptr>
4050 constexpr unsigned int
4051 get_communication_block_size(const VectorStruct &)
4052 {
4053 return numbers::invalid_unsigned_int;
4054 }
4055
4056
4057
4058 template <typename VectorStruct,
4059 std::enable_if_t<has_communication_block_size<VectorStruct>,
4060 VectorStruct> * = nullptr>
4061 constexpr unsigned int
4062 get_communication_block_size(const VectorStruct &)
4063 {
4064 return VectorStruct::communication_block_size;
4065 }
4066
4067
4068
4069 template <typename VectorType,
4070 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType> * =
4071 nullptr>
4072 bool
4073 has_ghost_elements(const VectorType &vec)
4074 {
4075 (void)vec;
4076 return false;
4077 }
4078
4079
4080
4081 template <typename VectorType,
4082 std::enable_if_t<!is_not_parallel_vector<VectorType>, VectorType>
4083 * = nullptr>
4084 bool
4085 has_ghost_elements(const VectorType &vec)
4086 {
4087 return vec.has_ghost_elements();
4088 }
4089
4090
4091
4092 // --------------------------------------------------------------------------
4093 // below we have wrappers to distinguish between block and non-block vectors.
4094 // --------------------------------------------------------------------------
4095
4096 //
4097 // update_ghost_values_start
4098 //
4099
4100 // update_ghost_values for block vectors
4101 template <int dim,
4102 typename VectorStruct,
4103 typename Number,
4104 typename VectorizedArrayType,
4105 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4106 * = nullptr>
4107 void
4108 update_ghost_values_start(
4109 const VectorStruct &vec,
4110 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4111 const unsigned int channel = 0)
4112 {
4113 if (get_communication_block_size(vec) < vec.n_blocks())
4114 {
4115 const bool ghosts_set = vec.has_ghost_elements();
4116
4117 Assert(exchanger.matrix_free.get_task_info()
4118 .allow_ghosted_vectors_in_loops ||
4119 ghosts_set == false,
4120 ExcNotImplemented());
4121
4122 if (ghosts_set)
4123 {
4124 exchanger.ghosts_were_set = true;
4125 return;
4126 }
4127
4128 vec.update_ghost_values();
4129 }
4130 else
4131 {
4132 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4133 update_ghost_values_start(vec.block(i), exchanger, channel + i);
4134 }
4135 }
4136
4137
4138
4139 // update_ghost_values for non-block vectors
4140 template <int dim,
4141 typename VectorStruct,
4142 typename Number,
4143 typename VectorizedArrayType,
4144 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4145 * = nullptr>
4146 void
4147 update_ghost_values_start(
4148 const VectorStruct &vec,
4149 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4150 const unsigned int channel = 0)
4151 {
4152 exchanger.update_ghost_values_start(channel, vec);
4153 }
4154
4155
4156
4157 // update_ghost_values_start() for vector of vectors
4158 template <int dim,
4159 typename VectorStruct,
4160 typename Number,
4161 typename VectorizedArrayType>
4162 inline void
4163 update_ghost_values_start(
4164 const std::vector<VectorStruct> &vec,
4165 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4166 {
4167 unsigned int component_index = 0;
4168 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4169 {
4170 update_ghost_values_start(vec[comp], exchanger, component_index);
4171 component_index += n_components(vec[comp]);
4172 }
4173 }
4174
4175
4176
4177 // update_ghost_values_start() for vector of pointers to vectors
4178 template <int dim,
4179 typename VectorStruct,
4180 typename Number,
4181 typename VectorizedArrayType>
4182 inline void
4183 update_ghost_values_start(
4184 const std::vector<VectorStruct *> &vec,
4185 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4186 {
4187 unsigned int component_index = 0;
4188 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4189 {
4190 update_ghost_values_start(*vec[comp], exchanger, component_index);
4191 component_index += n_components(*vec[comp]);
4192 }
4193 }
4194
4195
4196
4197 //
4198 // update_ghost_values_finish
4199 //
4200
4201 // for block vectors
4202 template <int dim,
4203 typename VectorStruct,
4204 typename Number,
4205 typename VectorizedArrayType,
4206 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4207 * = nullptr>
4208 void
4209 update_ghost_values_finish(
4210 const VectorStruct &vec,
4211 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4212 const unsigned int channel = 0)
4213 {
4214 if (get_communication_block_size(vec) < vec.n_blocks())
4215 {
4216 // do nothing, everything has already been completed in the _start()
4217 // call
4218 }
4219 else
4220 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4221 update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4222 }
4223
4224
4225
4226 // for non-block vectors
4227 template <int dim,
4228 typename VectorStruct,
4229 typename Number,
4230 typename VectorizedArrayType,
4231 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4232 * = nullptr>
4233 void
4234 update_ghost_values_finish(
4235 const VectorStruct &vec,
4236 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4237 const unsigned int channel = 0)
4238 {
4239 exchanger.update_ghost_values_finish(channel, vec);
4240 }
4241
4242
4243
4244 // for vector of vectors
4245 template <int dim,
4246 typename VectorStruct,
4247 typename Number,
4248 typename VectorizedArrayType>
4249 inline void
4250 update_ghost_values_finish(
4251 const std::vector<VectorStruct> &vec,
4252 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4253 {
4254 unsigned int component_index = 0;
4255 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4256 {
4257 update_ghost_values_finish(vec[comp], exchanger, component_index);
4258 component_index += n_components(vec[comp]);
4259 }
4260 }
4261
4262
4263
4264 // for vector of pointers to vectors
4265 template <int dim,
4266 typename VectorStruct,
4267 typename Number,
4268 typename VectorizedArrayType>
4269 inline void
4270 update_ghost_values_finish(
4271 const std::vector<VectorStruct *> &vec,
4272 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4273 {
4274 unsigned int component_index = 0;
4275 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4276 {
4277 update_ghost_values_finish(*vec[comp], exchanger, component_index);
4278 component_index += n_components(*vec[comp]);
4279 }
4280 }
4281
4282
4283
4284 //
4285 // compress_start
4286 //
4287
4288 // for block vectors
4289 template <int dim,
4290 typename VectorStruct,
4291 typename Number,
4292 typename VectorizedArrayType,
4293 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4294 * = nullptr>
4295 inline void
4296 compress_start(
4297 VectorStruct &vec,
4298 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4299 const unsigned int channel = 0)
4300 {
4301 if (get_communication_block_size(vec) < vec.n_blocks())
4302 vec.compress(VectorOperation::add);
4303 else
4304 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4305 compress_start(vec.block(i), exchanger, channel + i);
4306 }
4307
4308
4309
4310 // for non-block vectors
4311 template <int dim,
4312 typename VectorStruct,
4313 typename Number,
4314 typename VectorizedArrayType,
4315 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4316 * = nullptr>
4317 inline void
4318 compress_start(
4319 VectorStruct &vec,
4320 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4321 const unsigned int channel = 0)
4322 {
4323 exchanger.compress_start(channel, vec);
4324 }
4325
4326
4327
4328 // for std::vector of vectors
4329 template <int dim,
4330 typename VectorStruct,
4331 typename Number,
4332 typename VectorizedArrayType>
4333 inline void
4334 compress_start(
4335 std::vector<VectorStruct> &vec,
4336 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4337 {
4338 unsigned int component_index = 0;
4339 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4340 {
4341 compress_start(vec[comp], exchanger, component_index);
4342 component_index += n_components(vec[comp]);
4343 }
4344 }
4345
4346
4347
4348 // for std::vector of pointer to vectors
4349 template <int dim,
4350 typename VectorStruct,
4351 typename Number,
4352 typename VectorizedArrayType>
4353 inline void
4354 compress_start(
4355 std::vector<VectorStruct *> &vec,
4356 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4357 {
4358 unsigned int component_index = 0;
4359 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4360 {
4361 compress_start(*vec[comp], exchanger, component_index);
4362 component_index += n_components(*vec[comp]);
4363 }
4364 }
4365
4366
4367
4368 //
4369 // compress_finish
4370 //
4371
4372 // for block vectors
4373 template <int dim,
4374 typename VectorStruct,
4375 typename Number,
4376 typename VectorizedArrayType,
4377 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4378 * = nullptr>
4379 inline void
4380 compress_finish(
4381 VectorStruct &vec,
4382 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4383 const unsigned int channel = 0)
4384 {
4385 if (get_communication_block_size(vec) < vec.n_blocks())
4386 {
4387 // do nothing, everything has already been completed in the _start()
4388 // call
4389 }
4390 else
4391 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4392 compress_finish(vec.block(i), exchanger, channel + i);
4393 }
4394
4395
4396
4397 // for non-block vectors
4398 template <int dim,
4399 typename VectorStruct,
4400 typename Number,
4401 typename VectorizedArrayType,
4402 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4403 * = nullptr>
4404 inline void
4405 compress_finish(
4406 VectorStruct &vec,
4407 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4408 const unsigned int channel = 0)
4409 {
4410 exchanger.compress_finish(channel, vec);
4411 }
4412
4413
4414
4415 // for std::vector of vectors
4416 template <int dim,
4417 typename VectorStruct,
4418 typename Number,
4419 typename VectorizedArrayType>
4420 inline void
4421 compress_finish(
4422 std::vector<VectorStruct> &vec,
4423 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4424 {
4425 unsigned int component_index = 0;
4426 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4427 {
4428 compress_finish(vec[comp], exchanger, component_index);
4429 component_index += n_components(vec[comp]);
4430 }
4431 }
4432
4433
4434
4435 // for std::vector of pointer to vectors
4436 template <int dim,
4437 typename VectorStruct,
4438 typename Number,
4439 typename VectorizedArrayType>
4440 inline void
4441 compress_finish(
4442 std::vector<VectorStruct *> &vec,
4443 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4444 {
4445 unsigned int component_index = 0;
4446 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4447 {
4448 compress_finish(*vec[comp], exchanger, component_index);
4449 component_index += n_components(*vec[comp]);
4450 }
4451 }
4452
4453
4454
4455 //
4456 // reset_ghost_values:
4457 //
4458 // if the input vector did not have ghosts imported, clear them here again
4459 // in order to avoid subsequent operations e.g. in linear solvers to work
4460 // with ghosts all the time
4461 //
4462
4463 // for block vectors
4464 template <int dim,
4465 typename VectorStruct,
4466 typename Number,
4467 typename VectorizedArrayType,
4468 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4469 * = nullptr>
4470 inline void
4471 reset_ghost_values(
4472 const VectorStruct &vec,
4473 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4474 {
4475 // return immediately if there is nothing to do.
4476 if (exchanger.ghosts_were_set == true)
4477 return;
4478
4479 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4480 reset_ghost_values(vec.block(i), exchanger);
4481 }
4482
4483
4484
4485 // for non-block vectors
4486 template <int dim,
4487 typename VectorStruct,
4488 typename Number,
4489 typename VectorizedArrayType,
4490 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4491 * = nullptr>
4492 inline void
4493 reset_ghost_values(
4494 const VectorStruct &vec,
4495 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4496 {
4497 // return immediately if there is nothing to do.
4498 if (exchanger.ghosts_were_set == true)
4499 return;
4500
4501 exchanger.reset_ghost_values(vec);
4502 }
4503
4504
4505
4506 // for std::vector of vectors
4507 template <int dim,
4508 typename VectorStruct,
4509 typename Number,
4510 typename VectorizedArrayType>
4511 inline void
4512 reset_ghost_values(
4513 const std::vector<VectorStruct> &vec,
4514 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4515 {
4516 // return immediately if there is nothing to do.
4517 if (exchanger.ghosts_were_set == true)
4518 return;
4519
4520 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4521 reset_ghost_values(vec[comp], exchanger);
4522 }
4523
4524
4525
4526 // for std::vector of pointer to vectors
4527 template <int dim,
4528 typename VectorStruct,
4529 typename Number,
4530 typename VectorizedArrayType>
4531 inline void
4532 reset_ghost_values(
4533 const std::vector<VectorStruct *> &vec,
4534 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4535 {
4536 // return immediately if there is nothing to do.
4537 if (exchanger.ghosts_were_set == true)
4538 return;
4539
4540 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4541 reset_ghost_values(*vec[comp], exchanger);
4542 }
4543
4544
4545
4546 //
4547 // zero_vector_region
4548 //
4549
4550 // for block vectors
4551 template <int dim,
4552 typename VectorStruct,
4553 typename Number,
4554 typename VectorizedArrayType,
4555 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4556 * = nullptr>
4557 inline void
4558 zero_vector_region(
4559 const unsigned int range_index,
4560 VectorStruct &vec,
4561 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4562 {
4563 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4564 exchanger.zero_vector_region(range_index, vec.block(i));
4565 }
4566
4567
4568
4569 // for non-block vectors
4570 template <int dim,
4571 typename VectorStruct,
4572 typename Number,
4573 typename VectorizedArrayType,
4574 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4575 * = nullptr>
4576 inline void
4577 zero_vector_region(
4578 const unsigned int range_index,
4579 VectorStruct &vec,
4580 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4581 {
4582 exchanger.zero_vector_region(range_index, vec);
4583 }
4584
4585
4586
4587 // for std::vector of vectors
4588 template <int dim,
4589 typename VectorStruct,
4590 typename Number,
4591 typename VectorizedArrayType>
4592 inline void
4593 zero_vector_region(
4594 const unsigned int range_index,
4595 std::vector<VectorStruct> &vec,
4596 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4597 {
4598 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4599 zero_vector_region(range_index, vec[comp], exchanger);
4600 }
4601
4602
4603
4604 // for std::vector of pointers to vectors
4605 template <int dim,
4606 typename VectorStruct,
4607 typename Number,
4608 typename VectorizedArrayType>
4609 inline void
4610 zero_vector_region(
4611 const unsigned int range_index,
4612 std::vector<VectorStruct *> &vec,
4613 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4614 {
4615 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4616 zero_vector_region(range_index, *vec[comp], exchanger);
4617 }
4618
4619
4620
4621 // Apply a unit matrix operation to constrained DoFs: Default cases where we
4622 // cannot detect a LinearAlgebra::distributed::Vector, we do not do
4623 // anything, else we apply the constraints as a unit operation
4624 template <typename VectorStruct1, typename VectorStruct2>
4625 inline void
4626 apply_operation_to_constrained_dofs(const std::vector<unsigned int> &,
4627 const VectorStruct1 &,
4628 VectorStruct2 &)
4629 {}
4630
4631 template <typename Number>
4632 inline void
4633 apply_operation_to_constrained_dofs(
4634 const std::vector<unsigned int> &constrained_dofs,
4635 const LinearAlgebra::distributed::Vector<Number> &src,
4636 LinearAlgebra::distributed::Vector<Number> &dst)
4637 {
4638 for (const unsigned int i : constrained_dofs)
4639 dst.local_element(i) = src.local_element(i);
4640 }
4641
4642
4643 namespace MatrixFreeFunctions
4644 {
4645 // struct to select between a const interface and a non-const interface
4646 // for MFWorker
4647 template <typename, typename, typename, typename, bool>
4648 struct InterfaceSelector
4649 {};
4650
4651 // Version for constant functions
4652 template <typename MF,
4653 typename InVector,
4654 typename OutVector,
4655 typename Container>
4656 struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4657 {
4658 using function_type = void (Container::*)(
4659 const MF &,
4660 OutVector &,
4661 const InVector &,
4662 const std::pair<unsigned int, unsigned int> &) const;
4663 };
4664
4665 // Version for non-constant functions
4666 template <typename MF,
4667 typename InVector,
4668 typename OutVector,
4669 typename Container>
4670 struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4671 {
4672 using function_type =
4673 void (Container::*)(const MF &,
4674 OutVector &,
4675 const InVector &,
4676 const std::pair<unsigned int, unsigned int> &);
4677 };
4678 } // namespace MatrixFreeFunctions
4679
4680
4681
4682 // A implementation class for the worker object that runs the various
4683 // operations we want to perform during the matrix-free loop
4684 template <typename MF,
4685 typename InVector,
4686 typename OutVector,
4687 typename Container,
4688 bool is_constant>
4689 class MFWorker : public MFWorkerInterface
4690 {
4691 public:
4692 // An alias to make the arguments further down more readable
4693 using function_type = typename MatrixFreeFunctions::
4694 InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4695 function_type;
4696
4697 // constructor, binds all the arguments to this class
4698 MFWorker(const MF &matrix_free,
4699 const InVector &src,
4700 OutVector &dst,
4701 const bool zero_dst_vector_setting,
4702 const Container &container,
4703 function_type cell_function,
4704 function_type face_function,
4705 function_type boundary_function,
4706 const typename MF::DataAccessOnFaces src_vector_face_access =
4707 MF::DataAccessOnFaces::none,
4708 const typename MF::DataAccessOnFaces dst_vector_face_access =
4709 MF::DataAccessOnFaces::none,
4710 const std::function<void(const unsigned int, const unsigned int)>
4711 &operation_before_loop = {},
4712 const std::function<void(const unsigned int, const unsigned int)>
4713 &operation_after_loop = {},
4714 const unsigned int dof_handler_index_pre_post = 0)
4715 : matrix_free(matrix_free)
4716 , container(const_cast<Container &>(container))
4717 , cell_function(cell_function)
4718 , face_function(face_function)
4719 , boundary_function(boundary_function)
4720 , src(src)
4721 , dst(dst)
4722 , src_data_exchanger(matrix_free,
4723 src_vector_face_access,
4724 n_components(src))
4725 , dst_data_exchanger(matrix_free,
4726 dst_vector_face_access,
4727 n_components(dst))
4728 , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4729 , zero_dst_vector_setting(zero_dst_vector_setting &&
4730 !src_and_dst_are_same)
4731 , operation_before_loop(operation_before_loop)
4732 , operation_after_loop(operation_after_loop)
4733 , dof_handler_index_pre_post(dof_handler_index_pre_post)
4734 {
4735 Assert(!has_ghost_elements(dst),
4736 ExcMessage("The destination vector passed to the matrix-free "
4737 "loop is ghosted. This is not allowed."));
4738 }
4739
4740 // Runs the cell work. If no function is given, nothing is done
4741 virtual void
4742 cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4743 {
4744 if (cell_function != nullptr && cell_range.second > cell_range.first)
4745 for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4746 {
4747 const auto cell_subrange =
4748 matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4749
4750 if (cell_subrange.second <= cell_subrange.first)
4751 continue;
4752
4753 (container.*
4754 cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4755 }
4756 }
4757
4758 virtual void
4759 cell(const unsigned int range_index) override
4760 {
4761 process_range(cell_function,
4762 matrix_free.get_task_info().cell_partition_data_hp_ptr,
4763 matrix_free.get_task_info().cell_partition_data_hp,
4764 range_index);
4765 }
4766
4767 virtual void
4768 face(const unsigned int range_index) override
4769 {
4770 process_range(face_function,
4771 matrix_free.get_task_info().face_partition_data_hp_ptr,
4772 matrix_free.get_task_info().face_partition_data_hp,
4773 range_index);
4774 }
4775
4776 virtual void
4777 boundary(const unsigned int range_index) override
4778 {
4779 process_range(boundary_function,
4780 matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4781 matrix_free.get_task_info().boundary_partition_data_hp,
4782 range_index);
4783 }
4784
4785 private:
4786 void
4787 process_range(const function_type &fu,
4788 const std::vector<unsigned int> &ptr,
4789 const std::vector<unsigned int> &data,
4790 const unsigned int range_index)
4791 {
4792 if (fu == nullptr)
4793 return;
4794
4795 AssertIndexRange(range_index + 1, ptr.size());
4796 for (unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4797 {
4798 AssertIndexRange(2 * i + 1, data.size());
4799 (container.*fu)(matrix_free,
4800 this->dst,
4801 this->src,
4802 std::make_pair(data[2 * i], data[2 * i + 1]));
4803 }
4804 }
4805
4806 public:
4807 // Starts the communication for the update ghost values operation. We
4808 // cannot call this update if ghost and destination are the same because
4809 // that would introduce spurious entries in the destination (there is also
4810 // the problem that reading from a vector that we also write to is usually
4811 // not intended in case there is overlap, but this is up to the
4812 // application code to decide and we cannot catch this case here).
4813 virtual void
4814 vector_update_ghosts_start() override
4815 {
4816 if (!src_and_dst_are_same)
4817 internal::update_ghost_values_start(src, src_data_exchanger);
4818 }
4819
4820 // Finishes the communication for the update ghost values operation
4821 virtual void
4822 vector_update_ghosts_finish() override
4823 {
4824 if (!src_and_dst_are_same)
4825 internal::update_ghost_values_finish(src, src_data_exchanger);
4826 }
4827
4828 // Starts the communication for the vector compress operation
4829 virtual void
4830 vector_compress_start() override
4831 {
4832 internal::compress_start(dst, dst_data_exchanger);
4833 }
4834
4835 // Finishes the communication for the vector compress operation
4836 virtual void
4837 vector_compress_finish() override
4838 {
4839 internal::compress_finish(dst, dst_data_exchanger);
4840 if (!src_and_dst_are_same)
4841 internal::reset_ghost_values(src, src_data_exchanger);
4842 }
4843
4844 // Zeros the given input vector
4845 virtual void
4846 zero_dst_vector_range(const unsigned int range_index) override
4847 {
4848 if (zero_dst_vector_setting)
4849 internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4850 }
4851
4852 virtual void
4853 cell_loop_pre_range(const unsigned int range_index) override
4854 {
4855 if (operation_before_loop)
4856 {
4857 const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4858 matrix_free.get_dof_info(dof_handler_index_pre_post);
4859 if (range_index == numbers::invalid_unsigned_int)
4860 {
4861 // Case with threaded loop -> currently no overlap implemented
4862 ::parallel::apply_to_subranges(
4863 0U,
4864 dof_info.vector_partitioner->locally_owned_size(),
4865 operation_before_loop,
4866 internal::VectorImplementation::minimum_parallel_grain_size);
4867 }
4868 else
4869 {
4870 AssertIndexRange(range_index,
4871 dof_info.cell_loop_pre_list_index.size() - 1);
4872 for (unsigned int id =
4873 dof_info.cell_loop_pre_list_index[range_index];
4874 id != dof_info.cell_loop_pre_list_index[range_index + 1];
4875 ++id)
4876 operation_before_loop(dof_info.cell_loop_pre_list[id].first,
4877 dof_info.cell_loop_pre_list[id].second);
4878 }
4879 }
4880 }
4881
4882 virtual void
4883 cell_loop_post_range(const unsigned int range_index) override
4884 {
4885 if (operation_after_loop)
4886 {
4887 // Run unit matrix operation on constrained dofs if we are at the
4888 // last range
4889 const std::vector<unsigned int> &partition_row_index =
4890 matrix_free.get_task_info().partition_row_index;
4891 if (range_index ==
4892 partition_row_index[partition_row_index.size() - 2] - 1)
4893 apply_operation_to_constrained_dofs(
4894 matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
4895 src,
4896 dst);
4897
4898 const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4899 matrix_free.get_dof_info(dof_handler_index_pre_post);
4900 if (range_index == numbers::invalid_unsigned_int)
4901 {
4902 // Case with threaded loop -> currently no overlap implemented
4903 ::parallel::apply_to_subranges(
4904 0U,
4905 dof_info.vector_partitioner->locally_owned_size(),
4906 operation_after_loop,
4907 internal::VectorImplementation::minimum_parallel_grain_size);
4908 }
4909 else
4910 {
4911 AssertIndexRange(range_index,
4912 dof_info.cell_loop_post_list_index.size() - 1);
4913 for (unsigned int id =
4914 dof_info.cell_loop_post_list_index[range_index];
4915 id != dof_info.cell_loop_post_list_index[range_index + 1];
4916 ++id)
4917 operation_after_loop(dof_info.cell_loop_post_list[id].first,
4918 dof_info.cell_loop_post_list[id].second);
4919 }
4920 }
4921 }
4922
4923 private:
4924 const MF &matrix_free;
4925 Container &container;
4926 function_type cell_function;
4927 function_type face_function;
4928 function_type boundary_function;
4929
4930 const InVector &src;
4931 OutVector &dst;
4932 VectorDataExchange<MF::dimension,
4933 typename MF::value_type,
4934 typename MF::vectorized_value_type>
4935 src_data_exchanger;
4936 VectorDataExchange<MF::dimension,
4937 typename MF::value_type,
4938 typename MF::vectorized_value_type>
4939 dst_data_exchanger;
4940 const bool src_and_dst_are_same;
4941 const bool zero_dst_vector_setting;
4942 const std::function<void(const unsigned int, const unsigned int)>
4943 operation_before_loop;
4944 const std::function<void(const unsigned int, const unsigned int)>
4945 operation_after_loop;
4946 const unsigned int dof_handler_index_pre_post;
4947 };
4948
4949
4950
4955 template <class MF, typename InVector, typename OutVector>
4956 struct MFClassWrapper
4957 {
4958 using function_type =
4959 std::function<void(const MF &,
4960 OutVector &,
4961 const InVector &,
4962 const std::pair<unsigned int, unsigned int> &)>;
4963
4964 MFClassWrapper(const function_type cell,
4965 const function_type face,
4966 const function_type boundary)
4967 : cell(cell)
4968 , face(face)
4969 , boundary(boundary)
4970 {}
4971
4972 void
4973 cell_integrator(const MF &mf,
4974 OutVector &dst,
4975 const InVector &src,
4976 const std::pair<unsigned int, unsigned int> &range) const
4977 {
4978 if (cell)
4979 cell(mf, dst, src, range);
4980 }
4981
4982 void
4983 face_integrator(const MF &mf,
4984 OutVector &dst,
4985 const InVector &src,
4986 const std::pair<unsigned int, unsigned int> &range) const
4987 {
4988 if (face)
4989 face(mf, dst, src, range);
4990 }
4991
4992 void
4993 boundary_integrator(
4994 const MF &mf,
4995 OutVector &dst,
4996 const InVector &src,
4997 const std::pair<unsigned int, unsigned int> &range) const
4998 {
4999 if (boundary)
5000 boundary(mf, dst, src, range);
5001 }
5002
5003 const function_type cell;
5004 const function_type face;
5005 const function_type boundary;
5006 };
5007
5008} // end of namespace internal
5009
5010
5011
5012template <int dim, typename Number, typename VectorizedArrayType>
5013template <typename OutVector, typename InVector>
5014inline void
5015MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5016 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5017 OutVector &,
5018 const InVector &,
5019 const std::pair<unsigned int, unsigned int> &)>
5020 &cell_operation,
5021 OutVector &dst,
5022 const InVector &src,
5023 const bool zero_dst_vector) const
5024{
5025 using Wrapper =
5026 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5027 InVector,
5028 OutVector>;
5029 Wrapper wrap(cell_operation, nullptr, nullptr);
5030 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5031 InVector,
5032 OutVector,
5033 Wrapper,
5034 true>
5035 worker(*this,
5036 src,
5037 dst,
5038 zero_dst_vector,
5039 wrap,
5040 &Wrapper::cell_integrator,
5041 &Wrapper::face_integrator,
5042 &Wrapper::boundary_integrator);
5043
5044 task_info.loop(worker);
5045}
5046
5047
5048
5049template <int dim, typename Number, typename VectorizedArrayType>
5050template <typename OutVector, typename InVector>
5051inline void
5052MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5053 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5054 OutVector &,
5055 const InVector &,
5056 const std::pair<unsigned int, unsigned int> &)>
5057 &cell_operation,
5058 OutVector &dst,
5059 const InVector &src,
5060 const std::function<void(const unsigned int, const unsigned int)>
5061 &operation_before_loop,
5062 const std::function<void(const unsigned int, const unsigned int)>
5063 &operation_after_loop,
5064 const unsigned int dof_handler_index_pre_post) const
5065{
5066 using Wrapper =
5067 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5068 InVector,
5069 OutVector>;
5070 Wrapper wrap(cell_operation, nullptr, nullptr);
5071 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5072 InVector,
5073 OutVector,
5074 Wrapper,
5075 true>
5076 worker(*this,
5077 src,
5078 dst,
5079 false,
5080 wrap,
5081 &Wrapper::cell_integrator,
5082 &Wrapper::face_integrator,
5083 &Wrapper::boundary_integrator,
5084 DataAccessOnFaces::none,
5085 DataAccessOnFaces::none,
5086 operation_before_loop,
5087 operation_after_loop,
5088 dof_handler_index_pre_post);
5089
5090 task_info.loop(worker);
5091}
5092
5093
5094
5095template <int dim, typename Number, typename VectorizedArrayType>
5096template <typename OutVector, typename InVector>
5097inline void
5098MatrixFree<dim, Number, VectorizedArrayType>::loop(
5099 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5100 OutVector &,
5101 const InVector &,
5102 const std::pair<unsigned int, unsigned int> &)>
5103 &cell_operation,
5104 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5105 OutVector &,
5106 const InVector &,
5107 const std::pair<unsigned int, unsigned int> &)>
5108 &inner_face_operation,
5109 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5110 OutVector &,
5111 const InVector &,
5112 const std::pair<unsigned int, unsigned int> &)>
5113 &boundary_face_operation,
5114 OutVector &dst,
5115 const InVector &src,
5116 const bool zero_dst_vector,
5117 const DataAccessOnFaces dst_vector_face_access,
5118 const DataAccessOnFaces src_vector_face_access) const
5119{
5120 using Wrapper =
5121 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5122 InVector,
5123 OutVector>;
5124 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5125 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5126 InVector,
5127 OutVector,
5128 Wrapper,
5129 true>
5130 worker(*this,
5131 src,
5132 dst,
5133 zero_dst_vector,
5134 wrap,
5135 &Wrapper::cell_integrator,
5136 &Wrapper::face_integrator,
5137 &Wrapper::boundary_integrator,
5138 src_vector_face_access,
5139 dst_vector_face_access);
5140
5141 task_info.loop(worker);
5142}
5143
5144
5145
5146template <int dim, typename Number, typename VectorizedArrayType>
5147template <typename CLASS, typename OutVector, typename InVector>
5148inline void
5149MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5150 void (CLASS::*function_pointer)(
5151 const MatrixFree<dim, Number, VectorizedArrayType> &,
5152 OutVector &,
5153 const InVector &,
5154 const std::pair<unsigned int, unsigned int> &) const,
5155 const CLASS *owning_class,
5156 OutVector &dst,
5157 const InVector &src,
5158 const bool zero_dst_vector) const
5159{
5160 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5161 InVector,
5162 OutVector,
5163 CLASS,
5164 true>
5165 worker(*this,
5166 src,
5167 dst,
5168 zero_dst_vector,
5169 *owning_class,
5170 function_pointer,
5171 nullptr,
5172 nullptr);
5173 task_info.loop(worker);
5174}
5175
5176
5177
5178template <int dim, typename Number, typename VectorizedArrayType>
5179template <typename CLASS, typename OutVector, typename InVector>
5180inline void
5181MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5182 void (CLASS::*function_pointer)(
5183 const MatrixFree<dim, Number, VectorizedArrayType> &,
5184 OutVector &,
5185 const InVector &,
5186 const std::pair<unsigned int, unsigned int> &) const,
5187 const CLASS *owning_class,
5188 OutVector &dst,
5189 const InVector &src,
5190 const std::function<void(const unsigned int, const unsigned int)>
5191 &operation_before_loop,
5192 const std::function<void(const unsigned int, const unsigned int)>
5193 &operation_after_loop,
5194 const unsigned int dof_handler_index_pre_post) const
5195{
5196 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5197 InVector,
5198 OutVector,
5199 CLASS,
5200 true>
5201 worker(*this,
5202 src,
5203 dst,
5204 false,
5205 *owning_class,
5206 function_pointer,
5207 nullptr,
5208 nullptr,
5209 DataAccessOnFaces::none,
5210 DataAccessOnFaces::none,
5211 operation_before_loop,
5212 operation_after_loop,
5213 dof_handler_index_pre_post);
5214 task_info.loop(worker);
5215}
5216
5217
5218
5219template <int dim, typename Number, typename VectorizedArrayType>
5220template <typename CLASS, typename OutVector, typename InVector>
5221inline void
5222MatrixFree<dim, Number, VectorizedArrayType>::loop(
5223 void (CLASS::*cell_operation)(
5224 const MatrixFree<dim, Number, VectorizedArrayType> &,
5225 OutVector &,
5226 const InVector &,
5227 const std::pair<unsigned int, unsigned int> &) const,
5228 void (CLASS::*inner_face_operation)(
5229 const MatrixFree<dim, Number, VectorizedArrayType> &,
5230 OutVector &,
5231 const InVector &,
5232 const std::pair<unsigned int, unsigned int> &) const,
5233 void (CLASS::*boundary_face_operation)(
5234 const MatrixFree<dim, Number, VectorizedArrayType> &,
5235 OutVector &,
5236 const InVector &,
5237 const std::pair<unsigned int, unsigned int> &) const,
5238 const CLASS *owning_class,
5239 OutVector &dst,
5240 const InVector &src,
5241 const bool zero_dst_vector,
5242 const DataAccessOnFaces dst_vector_face_access,
5243 const DataAccessOnFaces src_vector_face_access) const
5244{
5245 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5246 InVector,
5247 OutVector,
5248 CLASS,
5249 true>
5250 worker(*this,
5251 src,
5252 dst,
5253 zero_dst_vector,
5254 *owning_class,
5255 cell_operation,
5256 inner_face_operation,
5257 boundary_face_operation,
5258 src_vector_face_access,
5259 dst_vector_face_access);
5260 task_info.loop(worker);
5261}
5262
5263
5264
5265template <int dim, typename Number, typename VectorizedArrayType>
5266template <typename CLASS, typename OutVector, typename InVector>
5267inline void
5268MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5269 void (CLASS::*function_pointer)(
5270 const MatrixFree<dim, Number, VectorizedArrayType> &,
5271 OutVector &,
5272 const InVector &,
5273 const std::pair<unsigned int, unsigned int> &),
5274 CLASS *owning_class,
5275 OutVector &dst,
5276 const InVector &src,
5277 const bool zero_dst_vector) const
5278{
5279 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5280 InVector,
5281 OutVector,
5282 CLASS,
5283 false>
5284 worker(*this,
5285 src,
5286 dst,
5287 zero_dst_vector,
5288 *owning_class,
5289 function_pointer,
5290 nullptr,
5291 nullptr);
5292 task_info.loop(worker);
5293}
5294
5295
5296
5297template <int dim, typename Number, typename VectorizedArrayType>
5298template <typename CLASS, typename OutVector, typename InVector>
5299inline void
5300MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5301 void (CLASS::*function_pointer)(
5302 const MatrixFree<dim, Number, VectorizedArrayType> &,
5303 OutVector &,
5304 const InVector &,
5305 const std::pair<unsigned int, unsigned int> &),
5306 CLASS *owning_class,
5307 OutVector &dst,
5308 const InVector &src,
5309 const std::function<void(const unsigned int, const unsigned int)>
5310 &operation_before_loop,
5311 const std::function<void(const unsigned int, const unsigned int)>
5312 &operation_after_loop,
5313 const unsigned int dof_handler_index_pre_post) const
5314{
5315 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5316 InVector,
5317 OutVector,
5318 CLASS,
5319 false>
5320 worker(*this,
5321 src,
5322 dst,
5323 false,
5324 *owning_class,
5325 function_pointer,
5326 nullptr,
5327 nullptr,
5328 DataAccessOnFaces::none,
5329 DataAccessOnFaces::none,
5330 operation_before_loop,
5331 operation_after_loop,
5332 dof_handler_index_pre_post);
5333 task_info.loop(worker);
5334}
5335
5336
5337
5338template <int dim, typename Number, typename VectorizedArrayType>
5339template <typename CLASS, typename OutVector, typename InVector>
5340inline void
5341MatrixFree<dim, Number, VectorizedArrayType>::loop(
5342 void (CLASS::*cell_operation)(
5343 const MatrixFree<dim, Number, VectorizedArrayType> &,
5344 OutVector &,
5345 const InVector &,
5346 const std::pair<unsigned int, unsigned int> &),
5347 void (CLASS::*inner_face_operation)(
5348 const MatrixFree<dim, Number, VectorizedArrayType> &,
5349 OutVector &,
5350 const InVector &,
5351 const std::pair<unsigned int, unsigned int> &),
5352 void (CLASS::*boundary_face_operation)(
5353 const MatrixFree<dim, Number, VectorizedArrayType> &,
5354 OutVector &,
5355 const InVector &,
5356 const std::pair<unsigned int, unsigned int> &),
5357 CLASS *owning_class,
5358 OutVector &dst,
5359 const InVector &src,
5360 const bool zero_dst_vector,
5361 const DataAccessOnFaces dst_vector_face_access,
5362 const DataAccessOnFaces src_vector_face_access) const
5363{
5364 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5365 InVector,
5366 OutVector,
5367 CLASS,
5368 false>
5369 worker(*this,
5370 src,
5371 dst,
5372 zero_dst_vector,
5373 *owning_class,
5374 cell_operation,
5375 inner_face_operation,
5376 boundary_face_operation,
5377 src_vector_face_access,
5378 dst_vector_face_access);
5379 task_info.loop(worker);
5380}
5381
5382
5383
5384template <int dim, typename Number, typename VectorizedArrayType>
5385template <typename OutVector, typename InVector>
5386inline void
5387MatrixFree<dim, Number, VectorizedArrayType>::loop(
5388 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5389 OutVector &,
5390 const InVector &,
5391 const std::pair<unsigned int, unsigned int> &)>
5392 &cell_operation,
5393 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5394 OutVector &,
5395 const InVector &,
5396 const std::pair<unsigned int, unsigned int> &)>
5397 &inner_face_operation,
5398 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5399 OutVector &,
5400 const InVector &,
5401 const std::pair<unsigned int, unsigned int> &)>
5402 &boundary_face_operation,
5403 OutVector &dst,
5404 const InVector &src,
5405 const std::function<void(const unsigned int, const unsigned int)>
5406 &operation_before_loop,
5407 const std::function<void(const unsigned int, const unsigned int)>
5408 &operation_after_loop,
5409 const unsigned int dof_handler_index_pre_post,
5410 const DataAccessOnFaces dst_vector_face_access,
5411 const DataAccessOnFaces src_vector_face_access) const
5412{
5413 using Wrapper =
5414 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5415 InVector,
5416 OutVector>;
5417 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5418 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5419 InVector,
5420 OutVector,
5421 Wrapper,
5422 true>
5423 worker(*this,
5424 src,
5425 dst,
5426 false,
5427 wrap,
5428 &Wrapper::cell_integrator,
5429 &Wrapper::face_integrator,
5430 &Wrapper::boundary_integrator,
5431 src_vector_face_access,
5432 dst_vector_face_access,
5433 operation_before_loop,
5434 operation_after_loop,
5435 dof_handler_index_pre_post);
5436
5437 task_info.loop(worker);
5438}
5439
5440
5441
5442template <int dim, typename Number, typename VectorizedArrayType>
5443template <typename CLASS, typename OutVector, typename InVector>
5444inline void
5445MatrixFree<dim, Number, VectorizedArrayType>::loop(
5446 void (CLASS::*cell_operation)(const MatrixFree &,
5447 OutVector &,
5448 const InVector &,
5449 const std::pair<unsigned int, unsigned int> &)
5450 const,
5451 void (CLASS::*inner_face_operation)(
5452 const MatrixFree &,
5453 OutVector &,
5454 const InVector &,
5455 const std::pair<unsigned int, unsigned int> &) const,
5456 void (CLASS::*boundary_face_operation)(
5457 const MatrixFree &,
5458 OutVector &,
5459 const InVector &,
5460 const std::pair<unsigned int, unsigned int> &) const,
5461 const CLASS *owning_class,
5462 OutVector &dst,
5463 const InVector &src,
5464 const std::function<void(const unsigned int, const unsigned int)>
5465 &operation_before_loop,
5466 const std::function<void(const unsigned int, const unsigned int)>
5467 &operation_after_loop,
5468 const unsigned int dof_handler_index_pre_post,
5469 const DataAccessOnFaces dst_vector_face_access,
5470 const DataAccessOnFaces src_vector_face_access) const
5471{
5472 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5473 InVector,
5474 OutVector,
5475 CLASS,
5476 true>
5477 worker(*this,
5478 src,
5479 dst,
5480 false,
5481 *owning_class,
5482 cell_operation,
5483 inner_face_operation,
5484 boundary_face_operation,
5485 src_vector_face_access,
5486 dst_vector_face_access,
5487 operation_before_loop,
5488 operation_after_loop,
5489 dof_handler_index_pre_post);
5490 task_info.loop(worker);
5491}
5492
5493
5494
5495template <int dim, typename Number, typename VectorizedArrayType>
5496template <typename CLASS, typename OutVector, typename InVector>
5497inline void
5498MatrixFree<dim, Number, VectorizedArrayType>::loop(
5499 void (CLASS::*cell_operation)(const MatrixFree &,
5500 OutVector &,
5501 const InVector &,
5502 const std::pair<unsigned int, unsigned int> &),
5503 void (CLASS::*inner_face_operation)(
5504 const MatrixFree &,
5505 OutVector &,
5506 const InVector &,
5507 const std::pair<unsigned int, unsigned int> &),
5508 void (CLASS::*boundary_face_operation)(
5509 const MatrixFree &,
5510 OutVector &,
5511 const InVector &,
5512 const std::pair<unsigned int, unsigned int> &),
5513 const CLASS *owning_class,
5514 OutVector &dst,
5515 const InVector &src,
5516 const std::function<void(const unsigned int, const unsigned int)>
5517 &operation_before_loop,
5518 const std::function<void(const unsigned int, const unsigned int)>
5519 &operation_after_loop,
5520 const unsigned int dof_handler_index_pre_post,
5521 const DataAccessOnFaces dst_vector_face_access,
5522 const DataAccessOnFaces src_vector_face_access) const
5523{
5524 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5525 InVector,
5526 OutVector,
5527 CLASS,
5528 false>
5529 worker(*this,
5530 src,
5531 dst,
5532 false,
5533 *owning_class,
5534 cell_operation,
5535 inner_face_operation,
5536 boundary_face_operation,
5537 src_vector_face_access,
5538 dst_vector_face_access,
5539 operation_before_loop,
5540 operation_after_loop,
5541 dof_handler_index_pre_post);
5542 task_info.loop(worker);
5543}
5544
5545
5546
5547template <int dim, typename Number, typename VectorizedArrayType>
5548template <typename CLASS, typename OutVector, typename InVector>
5549inline void
5550MatrixFree<dim, Number, VectorizedArrayType>::loop_cell_centric(
5551 void (CLASS::*function_pointer)(
5552 const MatrixFree<dim, Number, VectorizedArrayType> &,
5553 OutVector &,
5554 const InVector &,
5555 const std::pair<unsigned int, unsigned int> &) const,
5556 const CLASS *owning_class,
5557 OutVector &dst,
5558 const InVector &src,
5559 const bool zero_dst_vector,
5560 const DataAccessOnFaces src_vector_face_access) const
5561{
5562 auto src_vector_face_access_temp = src_vector_face_access;
5563 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5564 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5565 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5566 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5567
5568 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5569 InVector,
5570 OutVector,
5571 CLASS,
5572 true>
5573 worker(*this,
5574 src,
5575 dst,
5576 zero_dst_vector,
5577 *owning_class,
5578 function_pointer,
5579 nullptr,
5580 nullptr,
5581 src_vector_face_access_temp,
5582 DataAccessOnFaces::none);
5583 task_info.loop(worker);
5584}
5585
5586
5587
5588template <int dim, typename Number, typename VectorizedArrayType>
5589template <typename CLASS, typename OutVector, typename InVector>
5590inline void
5591MatrixFree<dim, Number, VectorizedArrayType>::loop_cell_centric(
5592 void (CLASS::*function_pointer)(
5593 const MatrixFree<dim, Number, VectorizedArrayType> &,
5594 OutVector &,
5595 const InVector &,
5596 const std::pair<unsigned int, unsigned int> &),
5597 CLASS *owning_class,
5598 OutVector &dst,
5599 const InVector &src,
5600 const bool zero_dst_vector,
5601 const DataAccessOnFaces src_vector_face_access) const
5602{
5603 auto src_vector_face_access_temp = src_vector_face_access;
5604 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5605 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5606 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5607 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5608
5609 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5610 InVector,
5611 OutVector,
5612 CLASS,
5613 false>
5614 worker(*this,
5615 src,
5616 dst,
5617 zero_dst_vector,
5618 *owning_class,
5619 function_pointer,
5620 nullptr,
5621 nullptr,
5622 src_vector_face_access_temp,
5623 DataAccessOnFaces::none);
5624 task_info.loop(worker);
5625}
5626
5627
5628
5629template <int dim, typename Number, typename VectorizedArrayType>
5630template <typename OutVector, typename InVector>
5631inline void
5632MatrixFree<dim, Number, VectorizedArrayType>::loop_cell_centric(
5633 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5634 OutVector &,
5635 const InVector &,
5636 const std::pair<unsigned int, unsigned int> &)>
5637 &cell_operation,
5638 OutVector &dst,
5639 const InVector &src,
5640 const bool zero_dst_vector,
5641 const DataAccessOnFaces src_vector_face_access) const
5642{
5643 auto src_vector_face_access_temp = src_vector_face_access;
5644 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5645 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5646 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5647 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5648
5649 using Wrapper =
5650 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5651 InVector,
5652 OutVector>;
5653 Wrapper wrap(cell_operation, nullptr, nullptr);
5654
5655 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5656 InVector,
5657 OutVector,
5658 Wrapper,
5659 true>
5660 worker(*this,
5661 src,
5662 dst,
5663 zero_dst_vector,
5664 wrap,
5665 &Wrapper::cell_integrator,
5666 &Wrapper::face_integrator,
5667 &Wrapper::boundary_integrator,
5668 src_vector_face_access_temp,
5669 DataAccessOnFaces::none);
5670 task_info.loop(worker);
5671}
5672
5673
5674#endif // ifndef DOXYGEN
5675
5676
5677
5678DEAL_II_NAMESPACE_CLOSE
5679
5680#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:4624
first
Point< 2 > first
Definition grid_out.cc:4623
level
unsigned int level
Definition grid_out.cc:4626
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