Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
matrix_free.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// 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 Subscriptor
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<SmartPointer<const DoFHandler<dim>>> dof_handlers;
2256
2263 std::vector<SmartPointer<const AffineConstraints<Number>>> affine_constraints;
2264
2269 std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
2270
2277 std::vector<Number> constraint_pool_data;
2278
2283 std::vector<unsigned int> constraint_pool_row_index;
2284
2289 internal::MatrixFreeFunctions::MappingInfo<dim, Number, VectorizedArrayType>
2290 mapping_info;
2291
2295 Table<4, internal::MatrixFreeFunctions::ShapeInfo<Number>> shape_info;
2296
2303 std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
2304
2309 std::vector<unsigned int> mf_cell_indices;
2310
2317 unsigned int cell_level_index_end_local;
2318
2324 internal::MatrixFreeFunctions::TaskInfo task_info;
2325
2330 internal::MatrixFreeFunctions::FaceInfo<VectorizedArrayType::size()>
2331 face_info;
2332
2336 bool indices_are_initialized;
2337
2341 bool mapping_is_initialized;
2342
2350 mutable Threads::ThreadLocalStorage<
2351 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2352 scratch_pad;
2353
2358 mutable std::list<std::pair<bool, AlignedVector<Number>>>
2359 scratch_pad_non_threadsafe;
2360
2364 unsigned int mg_level;
2365};
2366
2367
2368
2369/*----------------------- Inline functions ----------------------------------*/
2370
2371#ifndef DOXYGEN
2372
2373
2374
2375template <int dim, typename Number, typename VectorizedArrayType>
2376template <typename T>
2377inline void
2378MatrixFree<dim, Number, VectorizedArrayType>::initialize_cell_data_vector(
2379 AlignedVector<T> &vec) const
2380{
2381 vec.resize(this->n_cell_batches() + this->n_ghost_cell_batches());
2382}
2383
2384
2385
2386template <int dim, typename Number, typename VectorizedArrayType>
2387template <typename T>
2388inline void
2389MatrixFree<dim, Number, VectorizedArrayType>::initialize_face_data_vector(
2390 AlignedVector<T> &vec) const
2391{
2392 vec.resize(this->n_inner_face_batches() + this->n_boundary_face_batches() +
2393 this->n_ghost_inner_face_batches());
2394}
2395
2396
2397
2398template <int dim, typename Number, typename VectorizedArrayType>
2399template <typename VectorType>
2400inline void
2401MatrixFree<dim, Number, VectorizedArrayType>::initialize_dof_vector(
2402 VectorType &vec,
2403 const unsigned int comp) const
2404{
2405 static_assert(IsBlockVector<VectorType>::value == false,
2406 "This function is not supported for block vectors.");
2407
2408 Assert(task_info.n_procs == 1,
2409 ExcMessage("This function can only be used in serial."));
2410
2411 AssertIndexRange(comp, n_components());
2412 vec.reinit(dof_info[comp].vector_partitioner->size());
2413}
2414
2415
2416
2417template <int dim, typename Number, typename VectorizedArrayType>
2418template <typename Number2>
2419inline void
2420MatrixFree<dim, Number, VectorizedArrayType>::initialize_dof_vector(
2421 LinearAlgebra::distributed::Vector<Number2> &vec,
2422 const unsigned int comp) const
2423{
2424 AssertIndexRange(comp, n_components());
2425 vec.reinit(dof_info[comp].vector_partitioner, task_info.communicator_sm);
2426}
2427
2428
2429
2430template <int dim, typename Number, typename VectorizedArrayType>
2431inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2432MatrixFree<dim, Number, VectorizedArrayType>::get_vector_partitioner(
2433 const unsigned int comp) const
2434{
2435 AssertIndexRange(comp, n_components());
2436 return dof_info[comp].vector_partitioner;
2437}
2438
2439
2440
2441template <int dim, typename Number, typename VectorizedArrayType>
2442inline const std::vector<unsigned int> &
2443MatrixFree<dim, Number, VectorizedArrayType>::get_constrained_dofs(
2444 const unsigned int comp) const
2445{
2446 AssertIndexRange(comp, n_components());
2447 return dof_info[comp].constrained_dofs;
2448}
2449
2450
2451
2452template <int dim, typename Number, typename VectorizedArrayType>
2453inline unsigned int
2454MatrixFree<dim, Number, VectorizedArrayType>::n_components() const
2455{
2456 AssertDimension(dof_handlers.size(), dof_info.size());
2457 return dof_handlers.size();
2458}
2459
2460
2461
2462template <int dim, typename Number, typename VectorizedArrayType>
2463inline unsigned int
2464MatrixFree<dim, Number, VectorizedArrayType>::n_base_elements(
2465 const unsigned int dof_no) const
2466{
2467 AssertDimension(dof_handlers.size(), dof_info.size());
2468 AssertIndexRange(dof_no, dof_handlers.size());
2469 return dof_handlers[dof_no]->get_fe().n_base_elements();
2470}
2471
2472
2473
2474template <int dim, typename Number, typename VectorizedArrayType>
2475inline const internal::MatrixFreeFunctions::TaskInfo &
2476MatrixFree<dim, Number, VectorizedArrayType>::get_task_info() const
2477{
2478 return task_info;
2479}
2480
2481
2482
2483template <int dim, typename Number, typename VectorizedArrayType>
2484inline unsigned int
2485MatrixFree<dim, Number, VectorizedArrayType>::n_physical_cells() const
2486{
2487 return task_info.n_active_cells;
2488}
2489
2490
2491
2492template <int dim, typename Number, typename VectorizedArrayType>
2493inline unsigned int
2494MatrixFree<dim, Number, VectorizedArrayType>::n_cell_batches() const
2495{
2496 return *(task_info.cell_partition_data.end() - 2);
2497}
2498
2499
2500
2501template <int dim, typename Number, typename VectorizedArrayType>
2502inline unsigned int
2503MatrixFree<dim, Number, VectorizedArrayType>::n_ghost_cell_batches() const
2504{
2505 return *(task_info.cell_partition_data.end() - 1) -
2506 *(task_info.cell_partition_data.end() - 2);
2507}
2508
2509
2510
2511template <int dim, typename Number, typename VectorizedArrayType>
2512inline unsigned int
2513MatrixFree<dim, Number, VectorizedArrayType>::n_inner_face_batches() const
2514{
2515 if (task_info.face_partition_data.empty())
2516 return 0;
2517 return task_info.face_partition_data.back();
2518}
2519
2520
2521
2522template <int dim, typename Number, typename VectorizedArrayType>
2523inline unsigned int
2524MatrixFree<dim, Number, VectorizedArrayType>::n_boundary_face_batches() const
2525{
2526 if (task_info.face_partition_data.empty())
2527 return 0;
2528 return task_info.boundary_partition_data.back() -
2529 task_info.face_partition_data.back();
2530}
2531
2532
2533
2534template <int dim, typename Number, typename VectorizedArrayType>
2535inline unsigned int
2536MatrixFree<dim, Number, VectorizedArrayType>::n_ghost_inner_face_batches() const
2537{
2538 if (task_info.face_partition_data.empty())
2539 return 0;
2540 return face_info.faces.size() - task_info.boundary_partition_data.back();
2541}
2542
2543
2544
2545template <int dim, typename Number, typename VectorizedArrayType>
2546inline types::boundary_id
2547MatrixFree<dim, Number, VectorizedArrayType>::get_boundary_id(
2548 const unsigned int face_batch_index) const
2549{
2550 Assert(face_batch_index >= task_info.boundary_partition_data[0] &&
2551 face_batch_index < task_info.boundary_partition_data.back(),
2552 ExcIndexRange(face_batch_index,
2553 task_info.boundary_partition_data[0],
2554 task_info.boundary_partition_data.back()));
2555 return types::boundary_id(face_info.faces[face_batch_index].exterior_face_no);
2556}
2557
2558
2559
2560template <int dim, typename Number, typename VectorizedArrayType>
2561inline std::array<types::boundary_id, VectorizedArrayType::size()>
2562MatrixFree<dim, Number, VectorizedArrayType>::get_faces_by_cells_boundary_id(
2563 const unsigned int cell_batch_index,
2564 const unsigned int face_number) const
2565{
2566 AssertIndexRange(cell_batch_index, n_cell_batches());
2567 AssertIndexRange(face_number, GeometryInfo<dim>::faces_per_cell);
2568 Assert(face_info.cell_and_face_boundary_id.size(0) >= n_cell_batches(),
2569 ExcNotInitialized());
2570 std::array<types::boundary_id, VectorizedArrayType::size()> result;
2571 result.fill(numbers::invalid_boundary_id);
2572 for (unsigned int v = 0;
2573 v < n_active_entries_per_cell_batch(cell_batch_index);
2574 ++v)
2575 result[v] =
2576 face_info.cell_and_face_boundary_id(cell_batch_index, face_number, v);
2577 return result;
2578}
2579
2580
2581
2582template <int dim, typename Number, typename VectorizedArrayType>
2583inline const internal::MatrixFreeFunctions::
2584 MappingInfo<dim, Number, VectorizedArrayType> &
2585 MatrixFree<dim, Number, VectorizedArrayType>::get_mapping_info() const
2586{
2587 return mapping_info;
2588}
2589
2590
2591
2592template <int dim, typename Number, typename VectorizedArrayType>
2593inline const internal::MatrixFreeFunctions::DoFInfo &
2594MatrixFree<dim, Number, VectorizedArrayType>::get_dof_info(
2595 const unsigned int dof_index) const
2596{
2597 AssertIndexRange(dof_index, n_components());
2598 return dof_info[dof_index];
2599}
2600
2601
2602
2603template <int dim, typename Number, typename VectorizedArrayType>
2604inline unsigned int
2605MatrixFree<dim, Number, VectorizedArrayType>::n_constraint_pool_entries() const
2606{
2607 return constraint_pool_row_index.size() - 1;
2608}
2609
2610
2611
2612template <int dim, typename Number, typename VectorizedArrayType>
2613inline const Number *
2614MatrixFree<dim, Number, VectorizedArrayType>::constraint_pool_begin(
2615 const unsigned int row) const
2616{
2617 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2618 return constraint_pool_data.empty() ?
2619 nullptr :
2620 constraint_pool_data.data() + constraint_pool_row_index[row];
2621}
2622
2623
2624
2625template <int dim, typename Number, typename VectorizedArrayType>
2626inline const Number *
2627MatrixFree<dim, Number, VectorizedArrayType>::constraint_pool_end(
2628 const unsigned int row) const
2629{
2630 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2631 return constraint_pool_data.empty() ?
2632 nullptr :
2633 constraint_pool_data.data() + constraint_pool_row_index[row + 1];
2634}
2635
2636
2637
2638template <int dim, typename Number, typename VectorizedArrayType>
2639inline std::pair<unsigned int, unsigned int>
2640MatrixFree<dim, Number, VectorizedArrayType>::create_cell_subrange_hp(
2641 const std::pair<unsigned int, unsigned int> &range,
2642 const unsigned int degree,
2643 const unsigned int dof_handler_component) const
2644{
2645 if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2646 {
2647 AssertDimension(
2648 dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2649 AssertDimension(
2650 dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2651 if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2652 return range;
2653 else
2654 return {range.second, range.second};
2655 }
2656
2657 const unsigned int fe_index =
2658 dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2659 if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2660 return {range.second, range.second};
2661 else
2662 return create_cell_subrange_hp_by_index(range,
2663 fe_index,
2664 dof_handler_component);
2665}
2666
2667
2668
2669template <int dim, typename Number, typename VectorizedArrayType>
2670inline bool
2671MatrixFree<dim, Number, VectorizedArrayType>::at_irregular_cell(
2672 const unsigned int cell_batch_index) const
2673{
2674 AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2675 return VectorizedArrayType::size() > 1 &&
2676 cell_level_index[(cell_batch_index + 1) * VectorizedArrayType::size() -
2677 1] == cell_level_index[(cell_batch_index + 1) *
2678 VectorizedArrayType::size() -
2679 2];
2680}
2681
2682
2683
2684template <int dim, typename Number, typename VectorizedArrayType>
2685unsigned int
2686MatrixFree<dim, Number, VectorizedArrayType>::n_active_fe_indices() const
2687{
2688 return shape_info.size(2);
2689}
2690
2691
2692template <int dim, typename Number, typename VectorizedArrayType>
2693unsigned int
2694MatrixFree<dim, Number, VectorizedArrayType>::get_cell_active_fe_index(
2695 const std::pair<unsigned int, unsigned int> range) const
2696{
2697 const auto &fe_indices = dof_info[0].cell_active_fe_index;
2698
2699 if (fe_indices.empty() == true ||
2700 dof_handlers[0]->get_fe_collection().size() == 1)
2701 return 0;
2702
2703 const auto index = fe_indices[range.first];
2704
2705 for (unsigned int i = range.first; i < range.second; ++i)
2706 AssertDimension(index, fe_indices[i]);
2707
2708 return index;
2709}
2710
2711
2712
2713template <int dim, typename Number, typename VectorizedArrayType>
2714unsigned int
2715MatrixFree<dim, Number, VectorizedArrayType>::get_face_active_fe_index(
2716 const std::pair<unsigned int, unsigned int> range,
2717 const bool is_interior_face) const
2718{
2719 const auto &fe_indices = dof_info[0].cell_active_fe_index;
2720
2721 if (fe_indices.empty() == true)
2722 return 0;
2723
2724 if (is_interior_face)
2725 {
2726 const unsigned int index =
2727 fe_indices[face_info.faces[range.first].cells_interior[0] /
2728 VectorizedArrayType::size()];
2729
2730 for (unsigned int i = range.first; i < range.second; ++i)
2731 AssertDimension(index,
2732 fe_indices[face_info.faces[i].cells_interior[0] /
2733 VectorizedArrayType::size()]);
2734
2735 return index;
2736 }
2737 else
2738 {
2739 const unsigned int index =
2740 fe_indices[face_info.faces[range.first].cells_exterior[0] /
2741 VectorizedArrayType::size()];
2742
2743 for (unsigned int i = range.first; i < range.second; ++i)
2744 AssertDimension(index,
2745 fe_indices[face_info.faces[i].cells_exterior[0] /
2746 VectorizedArrayType::size()]);
2747
2748 return index;
2749 }
2750}
2751
2752
2753
2754template <int dim, typename Number, typename VectorizedArrayType>
2755inline unsigned int
2756MatrixFree<dim, Number, VectorizedArrayType>::n_active_entries_per_cell_batch(
2757 const unsigned int cell_batch_index) const
2758{
2759 Assert(!dof_info.empty(), ExcNotInitialized());
2760 AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2761 const std::vector<unsigned char> &n_lanes_filled =
2762 dof_info[0].n_vectorization_lanes_filled
2763 [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell];
2764 AssertIndexRange(cell_batch_index, n_lanes_filled.size());
2765
2766 return n_lanes_filled[cell_batch_index];
2767}
2768
2769
2770
2771template <int dim, typename Number, typename VectorizedArrayType>
2772inline unsigned int
2773MatrixFree<dim, Number, VectorizedArrayType>::n_active_entries_per_face_batch(
2774 const unsigned int face_batch_index) const
2775{
2776 AssertIndexRange(face_batch_index, face_info.faces.size());
2777 Assert(!dof_info.empty(), ExcNotInitialized());
2778 const std::vector<unsigned char> &n_lanes_filled =
2779 dof_info[0].n_vectorization_lanes_filled
2780 [internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior];
2781 AssertIndexRange(face_batch_index, n_lanes_filled.size());
2782 return n_lanes_filled[face_batch_index];
2783}
2784
2785
2786
2787template <int dim, typename Number, typename VectorizedArrayType>
2788inline unsigned int
2789MatrixFree<dim, Number, VectorizedArrayType>::get_dofs_per_cell(
2790 const unsigned int dof_handler_index,
2791 const unsigned int active_fe_index) const
2792{
2793 return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2794}
2795
2796
2797
2798template <int dim, typename Number, typename VectorizedArrayType>
2799inline unsigned int
2800MatrixFree<dim, Number, VectorizedArrayType>::get_n_q_points(
2801 const unsigned int quad_index,
2802 const unsigned int active_fe_index) const
2803{
2804 AssertIndexRange(quad_index, mapping_info.cell_data.size());
2805 return mapping_info.cell_data[quad_index]
2806 .descriptor[active_fe_index]
2807 .n_q_points;
2808}
2809
2810
2811
2812template <int dim, typename Number, typename VectorizedArrayType>
2813inline unsigned int
2814MatrixFree<dim, Number, VectorizedArrayType>::get_dofs_per_face(
2815 const unsigned int dof_handler_index,
2816 const unsigned int active_fe_index) const
2817{
2818 return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2819}
2820
2821
2822
2823template <int dim, typename Number, typename VectorizedArrayType>
2824inline unsigned int
2825MatrixFree<dim, Number, VectorizedArrayType>::get_n_q_points_face(
2826 const unsigned int quad_index,
2827 const unsigned int active_fe_index) const
2828{
2829 AssertIndexRange(quad_index, mapping_info.face_data.size());
2830 return mapping_info.face_data[quad_index]
2831 .descriptor[active_fe_index]
2832 .n_q_points;
2833}
2834
2835
2836
2837template <int dim, typename Number, typename VectorizedArrayType>
2838inline const IndexSet &
2839MatrixFree<dim, Number, VectorizedArrayType>::get_locally_owned_set(
2840 const unsigned int dof_handler_index) const
2841{
2842 return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2843}
2844
2845
2846
2847template <int dim, typename Number, typename VectorizedArrayType>
2848inline const IndexSet &
2849MatrixFree<dim, Number, VectorizedArrayType>::get_ghost_set(
2850 const unsigned int dof_handler_index) const
2851{
2852 return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2853}
2854
2855
2856
2857template <int dim, typename Number, typename VectorizedArrayType>
2858inline const internal::MatrixFreeFunctions::ShapeInfo<Number> &
2859MatrixFree<dim, Number, VectorizedArrayType>::get_shape_info(
2860 const unsigned int dof_handler_index,
2861 const unsigned int index_quad,
2862 const unsigned int index_fe,
2863 const unsigned int active_fe_index,
2864 const unsigned int active_quad_index) const
2865{
2866 AssertIndexRange(dof_handler_index, dof_info.size());
2867 const unsigned int ind =
2868 dof_info[dof_handler_index].global_base_element_offset + index_fe;
2869 AssertIndexRange(ind, shape_info.size(0));
2870 AssertIndexRange(index_quad, shape_info.size(1));
2871 AssertIndexRange(active_fe_index, shape_info.size(2));
2872 AssertIndexRange(active_quad_index, shape_info.size(3));
2873 return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2874}
2875
2876
2877
2878template <int dim, typename Number, typename VectorizedArrayType>
2879inline const internal::MatrixFreeFunctions::FaceToCellTopology<
2880 VectorizedArrayType::size()> &
2881MatrixFree<dim, Number, VectorizedArrayType>::get_face_info(
2882 const unsigned int face_batch_index) const
2883{
2884 AssertIndexRange(face_batch_index, face_info.faces.size());
2885 return face_info.faces[face_batch_index];
2886}
2887
2888
2889
2890template <int dim, typename Number, typename VectorizedArrayType>
2891inline const Table<3, unsigned int> &
2892MatrixFree<dim, Number, VectorizedArrayType>::get_cell_and_face_to_plain_faces()
2893 const
2894{
2895 return face_info.cell_and_face_to_plain_faces;
2896}
2897
2898
2899
2900template <int dim, typename Number, typename VectorizedArrayType>
2901inline const Quadrature<dim> &
2902MatrixFree<dim, Number, VectorizedArrayType>::get_quadrature(
2903 const unsigned int quad_index,
2904 const unsigned int active_fe_index) const
2905{
2906 AssertIndexRange(quad_index, mapping_info.cell_data.size());
2907 return mapping_info.cell_data[quad_index]
2908 .descriptor[active_fe_index]
2909 .quadrature;
2910}
2911
2912
2913
2914template <int dim, typename Number, typename VectorizedArrayType>
2915inline const Quadrature<dim - 1> &
2916MatrixFree<dim, Number, VectorizedArrayType>::get_face_quadrature(
2917 const unsigned int quad_index,
2918 const unsigned int active_fe_index) const
2919{
2920 AssertIndexRange(quad_index, mapping_info.face_data.size());
2921 return mapping_info.face_data[quad_index]
2922 .descriptor[active_fe_index]
2923 .quadrature;
2924}
2925
2926
2927
2928template <int dim, typename Number, typename VectorizedArrayType>
2929inline unsigned int
2930MatrixFree<dim, Number, VectorizedArrayType>::get_cell_range_category(
2931 const std::pair<unsigned int, unsigned int> range) const
2932{
2933 auto result = get_cell_category(range.first);
2934
2935 for (unsigned int i = range.first; i < range.second; ++i)
2936 result = std::max(result, get_cell_category(i));
2937
2938 return result;
2939}
2940
2941
2942
2943template <int dim, typename Number, typename VectorizedArrayType>
2944inline std::pair<unsigned int, unsigned int>
2945MatrixFree<dim, Number, VectorizedArrayType>::get_face_range_category(
2946 const std::pair<unsigned int, unsigned int> range) const
2947{
2948 auto result = get_face_category(range.first);
2949
2950 for (unsigned int i = range.first; i < range.second; ++i)
2951 {
2952 result.first = std::max(result.first, get_face_category(i).first);
2953 result.second = std::max(result.second, get_face_category(i).second);
2954 }
2955
2956 return result;
2957}
2958
2959
2960
2961template <int dim, typename Number, typename VectorizedArrayType>
2962inline unsigned int
2963MatrixFree<dim, Number, VectorizedArrayType>::get_cell_category(
2964 const unsigned int cell_batch_index) const
2965{
2966 AssertIndexRange(0, dof_info.size());
2967 AssertIndexRange(cell_batch_index, dof_info[0].cell_active_fe_index.size());
2968 if (dof_info[0].cell_active_fe_index.empty())
2969 return 0;
2970 else
2971 return dof_info[0].cell_active_fe_index[cell_batch_index];
2972}
2973
2974
2975
2976template <int dim, typename Number, typename VectorizedArrayType>
2977inline std::pair<unsigned int, unsigned int>
2978MatrixFree<dim, Number, VectorizedArrayType>::get_face_category(
2979 const unsigned int face_batch_index) const
2980{
2981 AssertIndexRange(face_batch_index, face_info.faces.size());
2982 if (dof_info[0].cell_active_fe_index.empty())
2983 return std::make_pair(0U, 0U);
2984
2985 std::pair<unsigned int, unsigned int> result = std::make_pair(0U, 0U);
2986 for (unsigned int v = 0;
2987 v < VectorizedArrayType::size() &&
2988 face_info.faces[face_batch_index].cells_interior[v] !=
2989 numbers::invalid_unsigned_int;
2990 ++v)
2991 result.first = std::max(
2992 result.first,
2993 dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
2994 .cells_interior[v] /
2995 VectorizedArrayType::size()]);
2996 if (face_info.faces[face_batch_index].cells_exterior[0] !=
2997 numbers::invalid_unsigned_int)
2998 for (unsigned int v = 0;
2999 v < VectorizedArrayType::size() &&
3000 face_info.faces[face_batch_index].cells_exterior[v] !=
3001 numbers::invalid_unsigned_int;
3002 ++v)
3003 result.second = std::max(
3004 result.second,
3005 dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
3006 .cells_exterior[v] /
3007 VectorizedArrayType::size()]);
3008 else
3009 result.second = numbers::invalid_unsigned_int;
3010 return result;
3011}
3012
3013
3014
3015template <int dim, typename Number, typename VectorizedArrayType>
3016inline bool
3017MatrixFree<dim, Number, VectorizedArrayType>::indices_initialized() const
3018{
3019 return indices_are_initialized;
3020}
3021
3022
3023
3024template <int dim, typename Number, typename VectorizedArrayType>
3025inline bool
3026MatrixFree<dim, Number, VectorizedArrayType>::mapping_initialized() const
3027{
3028 return mapping_is_initialized;
3029}
3030
3031
3032template <int dim, typename Number, typename VectorizedArrayType>
3033inline unsigned int
3034MatrixFree<dim, Number, VectorizedArrayType>::get_mg_level() const
3035{
3036 return mg_level;
3037}
3038
3039
3040
3041template <int dim, typename Number, typename VectorizedArrayType>
3042AlignedVector<VectorizedArrayType> *
3043MatrixFree<dim, Number, VectorizedArrayType>::acquire_scratch_data() const
3044{
3045 using list_type =
3046 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3047 list_type &data = scratch_pad.get();
3048 for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
3049 if (it->first == false)
3050 {
3051 it->first = true;
3052 return &it->second;
3053 }
3054 data.emplace_front(true, AlignedVector<VectorizedArrayType>());
3055 return &data.front().second;
3056}
3057
3058
3059
3060template <int dim, typename Number, typename VectorizedArrayType>
3061void
3062MatrixFree<dim, Number, VectorizedArrayType>::release_scratch_data(
3063 const AlignedVector<VectorizedArrayType> *scratch) const
3064{
3065 using list_type =
3066 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3067 list_type &data = scratch_pad.get();
3068 for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
3069 if (&it->second == scratch)
3070 {
3071 Assert(it->first == true, ExcInternalError());
3072 it->first = false;
3073 return;
3074 }
3075 AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
3076}
3077
3078
3079
3080template <int dim, typename Number, typename VectorizedArrayType>
3081AlignedVector<Number> *
3082MatrixFree<dim, Number, VectorizedArrayType>::
3083 acquire_scratch_data_non_threadsafe() const
3084{
3085 for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
3086 scratch_pad_non_threadsafe.begin();
3087 it != scratch_pad_non_threadsafe.end();
3088 ++it)
3089 if (it->first == false)
3090 {
3091 it->first = true;
3092 return &it->second;
3093 }
3094 scratch_pad_non_threadsafe.push_front(
3095 std::make_pair(true, AlignedVector<Number>()));
3096 return &scratch_pad_non_threadsafe.front().second;
3097}
3098
3099
3100
3101template <int dim, typename Number, typename VectorizedArrayType>
3102void
3103MatrixFree<dim, Number, VectorizedArrayType>::
3104 release_scratch_data_non_threadsafe(
3105 const AlignedVector<Number> *scratch) const
3106{
3107 for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
3108 scratch_pad_non_threadsafe.begin();
3109 it != scratch_pad_non_threadsafe.end();
3110 ++it)
3111 if (&it->second == scratch)
3112 {
3113 Assert(it->first == true, ExcInternalError());
3114 it->first = false;
3115 return;
3116 }
3117 AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
3118}
3119
3120
3121
3122// ------------------------------ reinit functions ---------------------------
3123
3124namespace internal
3125{
3126 namespace MatrixFreeImplementation
3127 {
3128 template <int dim, int spacedim>
3129 inline std::vector<IndexSet>
3130 extract_locally_owned_index_sets(
3131 const std::vector<const DoFHandler<dim, spacedim> *> &dofh,
3132 const unsigned int level)
3133 {
3134 std::vector<IndexSet> locally_owned_set;
3135 locally_owned_set.reserve(dofh.size());
3136 for (unsigned int j = 0; j < dofh.size(); ++j)
3137 if (level == numbers::invalid_unsigned_int)
3138 locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
3139 else
3140 locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
3141 return locally_owned_set;
3142 }
3143 } // namespace MatrixFreeImplementation
3144} // namespace internal
3145
3146
3147
3148template <int dim, typename Number, typename VectorizedArrayType>
3149template <typename QuadratureType, typename number2, typename MappingType>
3150void
3151MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3152 const MappingType &mapping,
3153 const DoFHandler<dim> &dof_handler,
3154 const AffineConstraints<number2> &constraints_in,
3155 const QuadratureType &quad,
3156 const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
3157 &additional_data)
3158{
3159 std::vector<const DoFHandler<dim, dim> *> dof_handlers;
3160 std::vector<const AffineConstraints<number2> *> constraints;
3161
3162 dof_handlers.push_back(&dof_handler);
3163 constraints.push_back(&constraints_in);
3164
3165 std::vector<IndexSet> locally_owned_sets =
3166 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3167 dof_handlers, additional_data.mg_level);
3168
3169 std::vector<hp::QCollection<dim>> quad_hp;
3170 quad_hp.emplace_back(quad);
3171
3172 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3173 dof_handlers,
3174 constraints,
3175 locally_owned_sets,
3176 quad_hp,
3177 additional_data);
3178}
3179
3180
3181
3182template <int dim, typename Number, typename VectorizedArrayType>
3183template <typename QuadratureType, typename number2, typename MappingType>
3184void
3185MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3186 const MappingType &mapping,
3187 const std::vector<const DoFHandler<dim> *> &dof_handler,
3188 const std::vector<const AffineConstraints<number2> *> &constraint,
3189 const QuadratureType &quad,
3190 const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
3191 &additional_data)
3192{
3193 std::vector<IndexSet> locally_owned_set =
3194 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3195 dof_handler, additional_data.mg_level);
3196 std::vector<hp::QCollection<dim>> quad_hp;
3197 quad_hp.emplace_back(quad);
3198
3199 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3200 dof_handler,
3201 constraint,
3202 locally_owned_set,
3203 quad_hp,
3204 additional_data);
3205}
3206
3207
3208
3209template <int dim, typename Number, typename VectorizedArrayType>
3210template <typename QuadratureType, typename number2, typename MappingType>
3211void
3212MatrixFree<dim, Number, VectorizedArrayType>::reinit(
3213 const MappingType &mapping,
3214 const std::vector<const DoFHandler<dim> *> &dof_handler,
3215 const std::vector<const AffineConstraints<number2> *> &constraint,
3216 const std::vector<QuadratureType> &quad,
3217 const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
3218 &additional_data)
3219{
3220 std::vector<IndexSet> locally_owned_set =
3221 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3222 dof_handler, additional_data.mg_level);
3223 std::vector<hp::QCollection<dim>> quad_hp;
3224 for (unsigned int q = 0; q < quad.size(); ++q)
3225 quad_hp.emplace_back(quad[q]);
3226
3227 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3228 dof_handler,
3229 constraint,
3230 locally_owned_set,
3231 quad_hp,
3232 additional_data);
3233}
3234
3235
3236
3237// ------------------------------ implementation of loops --------------------
3238
3239// internal helper functions that define how to call MPI data exchange
3240// functions: for generic vectors, do nothing at all. For distributed vectors,
3241// call update_ghost_values_start function and so on. If we have collections
3242// of vectors, just do the individual functions of the components. In order to
3243// keep ghost values consistent (whether we are in read or write mode), we
3244// also reset the values at the end. the whole situation is a bit complicated
3245// by the fact that we need to treat block vectors differently, which use some
3246// additional helper functions to select the blocks and template magic.
3247namespace internal
3248{
3252 template <int dim, typename Number, typename VectorizedArrayType>
3253 struct VectorDataExchange
3254 {
3255 // A shift for the MPI messages to reduce the risk for accidental
3256 // interaction with other open communications that a user program might
3257 // set up (parallel vectors support unfinished communication). We let
3258 // the other vectors use the first 20 assigned numbers and start the
3259 // matrix-free communication.
3260 static constexpr unsigned int channel_shift = 20;
3261
3262
3263
3268 VectorDataExchange(
3269 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3270 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3271 DataAccessOnFaces vector_face_access,
3272 const unsigned int n_components)
3273 : matrix_free(matrix_free)
3274 , vector_face_access(
3275 matrix_free.get_task_info().face_partition_data.empty() ?
3276 ::MatrixFree<dim, Number, VectorizedArrayType>::
3277 DataAccessOnFaces::unspecified :
3278 vector_face_access)
3279 , ghosts_were_set(false)
3280# ifdef DEAL_II_WITH_MPI
3281 , tmp_data(n_components)
3282 , requests(n_components)
3283# endif
3284 {
3285 (void)n_components;
3286 if (this->vector_face_access !=
3287 ::MatrixFree<dim, Number, VectorizedArrayType>::
3288 DataAccessOnFaces::unspecified)
3289 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3290 AssertDimension(
3291 matrix_free.get_dof_info(c).vector_exchanger_face_variants.size(),
3292 5);
3293 }
3294
3295
3296
3300 ~VectorDataExchange() // NOLINT
3301 {
3302# ifdef DEAL_II_WITH_MPI
3303 for (unsigned int i = 0; i < tmp_data.size(); ++i)
3304 if (tmp_data[i] != nullptr)
3305 matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3306# endif
3307 }
3308
3309
3310
3315 template <typename VectorType>
3316 unsigned int
3317 find_vector_in_mf(const VectorType &vec,
3318 const bool check_global_compatibility = true) const
3319 {
3320 // case 1: vector was set up with MatrixFree::initialize_dof_vector()
3321 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3322 if (vec.get_partitioner().get() ==
3323 matrix_free.get_dof_info(c).vector_partitioner.get())
3324 return c;
3325
3326 // case 2: user provided own partitioner (compatibility mode)
3327 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3328 if (check_global_compatibility ?
3329 vec.get_partitioner()->is_globally_compatible(
3330 *matrix_free.get_dof_info(c).vector_partitioner) :
3331 vec.get_partitioner()->is_compatible(
3332 *matrix_free.get_dof_info(c).vector_partitioner))
3333 return c;
3334
3335 Assert(false,
3336 ExcNotImplemented("Could not find partitioner that fits vector"));
3337
3338 return numbers::invalid_unsigned_int;
3339 }
3340
3341
3342
3347 const internal::MatrixFreeFunctions::VectorDataExchange::Base &
3348 get_partitioner(const unsigned int mf_component) const
3349 {
3350 AssertDimension(matrix_free.get_dof_info(mf_component)
3351 .vector_exchanger_face_variants.size(),
3352 5);
3353 if (vector_face_access ==
3354 ::MatrixFree<dim, Number, VectorizedArrayType>::
3355 DataAccessOnFaces::none)
3356 return *matrix_free.get_dof_info(mf_component)
3357 .vector_exchanger_face_variants[0];
3358 else if (vector_face_access ==
3359 ::MatrixFree<dim, Number, VectorizedArrayType>::
3360 DataAccessOnFaces::values)
3361 return *matrix_free.get_dof_info(mf_component)
3362 .vector_exchanger_face_variants[1];
3363 else if (vector_face_access ==
3364 ::MatrixFree<dim, Number, VectorizedArrayType>::
3365 DataAccessOnFaces::gradients)
3366 return *matrix_free.get_dof_info(mf_component)
3367 .vector_exchanger_face_variants[2];
3368 else if (vector_face_access ==
3369 ::MatrixFree<dim, Number, VectorizedArrayType>::
3370 DataAccessOnFaces::values_all_faces)
3371 return *matrix_free.get_dof_info(mf_component)
3372 .vector_exchanger_face_variants[3];
3373 else if (vector_face_access ==
3374 ::MatrixFree<dim, Number, VectorizedArrayType>::
3375 DataAccessOnFaces::gradients_all_faces)
3376 return *matrix_free.get_dof_info(mf_component)
3377 .vector_exchanger_face_variants[4];
3378 else
3379 return *matrix_free.get_dof_info(mf_component).vector_exchanger.get();
3380 }
3381
3382
3383
3387 template <typename VectorType,
3388 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3389 * = nullptr>
3390 void
3391 update_ghost_values_start(const unsigned int /*component_in_block_vector*/,
3392 const VectorType & /*vec*/)
3393 {}
3394
3395
3400 template <typename VectorType,
3401 std::enable_if_t<!has_update_ghost_values_start<VectorType> &&
3402 !is_not_parallel_vector<VectorType>,
3403 VectorType> * = nullptr>
3404 void
3405 update_ghost_values_start(const unsigned int component_in_block_vector,
3406 const VectorType &vec)
3407 {
3408 (void)component_in_block_vector;
3409 const bool ghosts_set = vec.has_ghost_elements();
3410
3411 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3412 ghosts_set == false,
3413 ExcNotImplemented());
3414
3415 if (ghosts_set)
3416 {
3417 ghosts_were_set = true;
3418 return;
3419 }
3420
3421 vec.update_ghost_values();
3422 }
3423
3424
3425
3431 template <typename VectorType,
3432 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3433 !has_exchange_on_subset<VectorType>,
3434 VectorType> * = nullptr>
3435 void
3436 update_ghost_values_start(const unsigned int component_in_block_vector,
3437 const VectorType &vec)
3438 {
3439 (void)component_in_block_vector;
3440 const bool ghosts_set = vec.has_ghost_elements();
3441
3442 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3443 ghosts_set == false,
3444 ExcNotImplemented());
3445
3446 if (ghosts_set)
3447 {
3448 ghosts_were_set = true;
3449 return;
3450 }
3451
3452 vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3453 }
3454
3455
3456
3463 template <typename VectorType,
3464 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3465 has_exchange_on_subset<VectorType>,
3466 VectorType> * = nullptr>
3467 void
3468 update_ghost_values_start(const unsigned int component_in_block_vector,
3469 const VectorType &vec)
3470 {
3471 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3472 "Type mismatch between VectorType and VectorDataExchange");
3473 (void)component_in_block_vector;
3474 const bool ghosts_set = vec.has_ghost_elements();
3475
3476 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3477 ghosts_set == false,
3478 ExcNotImplemented());
3479
3480 if (ghosts_set)
3481 {
3482 ghosts_were_set = true;
3483 return;
3484 }
3485
3486 if (vec.size() != 0)
3487 {
3488# ifdef DEAL_II_WITH_MPI
3489 const unsigned int mf_component = find_vector_in_mf(vec);
3490
3491 const auto &part = get_partitioner(mf_component);
3492
3493 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3494 part.n_import_sm_procs() == 0)
3495 return;
3496
3497 tmp_data[component_in_block_vector] =
3498 matrix_free.acquire_scratch_data_non_threadsafe();
3499 tmp_data[component_in_block_vector]->resize_fast(
3500 part.n_import_indices());
3501 AssertDimension(requests.size(), tmp_data.size());
3502
3503 part.export_to_ghosted_array_start(
3504 component_in_block_vector * 2 + channel_shift,
3505 ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3506 vec.shared_vector_data(),
3507 ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3508 part.locally_owned_size(),
3509 matrix_free.get_dof_info(mf_component)
3510 .vector_partitioner->n_ghost_indices()),
3511 ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3512 part.n_import_indices()),
3513 this->requests[component_in_block_vector]);
3514# endif
3515 }
3516 }
3517
3518
3519
3524 template <typename VectorType,
3525 std::enable_if_t<!has_update_ghost_values_start<VectorType>,
3526 VectorType> * = nullptr>
3527 void
3528 update_ghost_values_finish(const unsigned int /*component_in_block_vector*/,
3529 const VectorType & /*vec*/)
3530 {}
3531
3532
3533
3539 template <typename VectorType,
3540 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3541 !has_exchange_on_subset<VectorType>,
3542 VectorType> * = nullptr>
3543 void
3544 update_ghost_values_finish(const unsigned int component_in_block_vector,
3545 const VectorType &vec)
3546 {
3547 (void)component_in_block_vector;
3548
3549 if (ghosts_were_set)
3550 return;
3551
3552 vec.update_ghost_values_finish();
3553 }
3554
3555
3556
3563 template <typename VectorType,
3564 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3565 has_exchange_on_subset<VectorType>,
3566 VectorType> * = nullptr>
3567 void
3568 update_ghost_values_finish(const unsigned int component_in_block_vector,
3569 const VectorType &vec)
3570 {
3571 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3572 "Type mismatch between VectorType and VectorDataExchange");
3573 (void)component_in_block_vector;
3574
3575 if (ghosts_were_set)
3576 return;
3577
3578 if (vec.size() != 0)
3579 {
3580# ifdef DEAL_II_WITH_MPI
3581 AssertIndexRange(component_in_block_vector, tmp_data.size());
3582 AssertDimension(requests.size(), tmp_data.size());
3583
3584 const unsigned int mf_component = find_vector_in_mf(vec);
3585
3586 const auto &part = get_partitioner(mf_component);
3587
3588 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3589 part.n_import_sm_procs() != 0)
3590 {
3591 part.export_to_ghosted_array_finish(
3592 ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3593 vec.shared_vector_data(),
3594 ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3595 part.locally_owned_size(),
3596 matrix_free.get_dof_info(mf_component)
3597 .vector_partitioner->n_ghost_indices()),
3598 this->requests[component_in_block_vector]);
3599
3600 matrix_free.release_scratch_data_non_threadsafe(
3601 tmp_data[component_in_block_vector]);
3602 tmp_data[component_in_block_vector] = nullptr;
3603 }
3604# endif
3605 }
3606 // let vector know that ghosts are being updated and we can read from
3607 // them
3608 vec.set_ghost_state(true);
3609 }
3610
3611
3612
3616 template <typename VectorType,
3617 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3618 * = nullptr>
3619 void
3620 compress_start(const unsigned int /*component_in_block_vector*/,
3621 VectorType & /*vec*/)
3622 {}
3623
3624
3625
3630 template <typename VectorType,
3631 std::enable_if_t<!has_compress_start<VectorType> &&
3632 !is_not_parallel_vector<VectorType>,
3633 VectorType> * = nullptr>
3634 void
3635 compress_start(const unsigned int component_in_block_vector,
3636 VectorType &vec)
3637 {
3638 (void)component_in_block_vector;
3639 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3640 vec.compress(VectorOperation::add);
3641 }
3642
3643
3644
3650 template <typename VectorType,
3651 std::enable_if_t<has_compress_start<VectorType> &&
3652 !has_exchange_on_subset<VectorType>,
3653 VectorType> * = nullptr>
3654 void
3655 compress_start(const unsigned int component_in_block_vector,
3656 VectorType &vec)
3657 {
3658 (void)component_in_block_vector;
3659 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3660 vec.compress_start(component_in_block_vector + channel_shift);
3661 }
3662
3663
3664
3671 template <typename VectorType,
3672 std::enable_if_t<has_compress_start<VectorType> &&
3673 has_exchange_on_subset<VectorType>,
3674 VectorType> * = nullptr>
3675 void
3676 compress_start(const unsigned int component_in_block_vector,
3677 VectorType &vec)
3678 {
3679 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3680 "Type mismatch between VectorType and VectorDataExchange");
3681 (void)component_in_block_vector;
3682 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3683
3684 if (vec.size() != 0)
3685 {
3686# ifdef DEAL_II_WITH_MPI
3687 const unsigned int mf_component = find_vector_in_mf(vec);
3688
3689 const auto &part = get_partitioner(mf_component);
3690
3691 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3692 part.n_import_sm_procs() == 0)
3693 return;
3694
3695 tmp_data[component_in_block_vector] =
3696 matrix_free.acquire_scratch_data_non_threadsafe();
3697 tmp_data[component_in_block_vector]->resize_fast(
3698 part.n_import_indices());
3699 AssertDimension(requests.size(), tmp_data.size());
3700
3701 part.import_from_ghosted_array_start(
3702 VectorOperation::add,
3703 component_in_block_vector * 2 + channel_shift,
3704 ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3705 vec.shared_vector_data(),
3706 ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3707 matrix_free.get_dof_info(mf_component)
3708 .vector_partitioner->n_ghost_indices()),
3709 ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3710 part.n_import_indices()),
3711 this->requests[component_in_block_vector]);
3712# endif
3713 }
3714 }
3715
3716
3717
3722 template <
3723 typename VectorType,
3724 std::enable_if_t<!has_compress_start<VectorType>, VectorType> * = nullptr>
3725 void
3726 compress_finish(const unsigned int /*component_in_block_vector*/,
3727 VectorType & /*vec*/)
3728 {}
3729
3730
3731
3737 template <typename VectorType,
3738 std::enable_if_t<has_compress_start<VectorType> &&
3739 !has_exchange_on_subset<VectorType>,
3740 VectorType> * = nullptr>
3741 void
3742 compress_finish(const unsigned int component_in_block_vector,
3743 VectorType &vec)
3744 {
3745 (void)component_in_block_vector;
3746 vec.compress_finish(VectorOperation::add);
3747 }
3748
3749
3750
3757 template <typename VectorType,
3758 std::enable_if_t<has_compress_start<VectorType> &&
3759 has_exchange_on_subset<VectorType>,
3760 VectorType> * = nullptr>
3761 void
3762 compress_finish(const unsigned int component_in_block_vector,
3763 VectorType &vec)
3764 {
3765 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3766 "Type mismatch between VectorType and VectorDataExchange");
3767 (void)component_in_block_vector;
3768 if (vec.size() != 0)
3769 {
3770# ifdef DEAL_II_WITH_MPI
3771 AssertIndexRange(component_in_block_vector, tmp_data.size());
3772 AssertDimension(requests.size(), tmp_data.size());
3773
3774 const unsigned int mf_component = find_vector_in_mf(vec);
3775
3776 const auto &part = get_partitioner(mf_component);
3777
3778 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3779 part.n_import_sm_procs() != 0)
3780 {
3781 part.import_from_ghosted_array_finish(
3782 VectorOperation::add,
3783 ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3784 vec.shared_vector_data(),
3785 ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3786 matrix_free.get_dof_info(mf_component)
3787 .vector_partitioner->n_ghost_indices()),
3788 ArrayView<const Number>(
3789 tmp_data[component_in_block_vector]->begin(),
3790 part.n_import_indices()),
3791 this->requests[component_in_block_vector]);
3792
3793 matrix_free.release_scratch_data_non_threadsafe(
3794 tmp_data[component_in_block_vector]);
3795 tmp_data[component_in_block_vector] = nullptr;
3796 }
3797
3798 if (Utilities::MPI::job_supports_mpi())
3799 {
3800 const int ierr =
3801 MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3802 AssertThrowMPI(ierr);
3803 }
3804# endif
3805 }
3806 }
3807
3808
3809
3813 template <typename VectorType,
3814 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3815 * = nullptr>
3816 void
3817 reset_ghost_values(const VectorType & /*vec*/) const
3818 {}
3819
3820
3821
3826 template <typename VectorType,
3827 std::enable_if_t<!has_exchange_on_subset<VectorType> &&
3828 !is_not_parallel_vector<VectorType>,
3829 VectorType> * = nullptr>
3830 void
3831 reset_ghost_values(const VectorType &vec) const
3832 {
3833 if (ghosts_were_set == true)
3834 return;
3835
3836 vec.zero_out_ghost_values();
3837 }
3838
3839
3840
3846 template <typename VectorType,
3847 std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3848 * = nullptr>
3849 void
3850 reset_ghost_values(const VectorType &vec) const
3851 {
3852 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3853 "Type mismatch between VectorType and VectorDataExchange");
3854 if (ghosts_were_set == true)
3855 return;
3856
3857 if (vec.size() != 0)
3858 {
3859# ifdef DEAL_II_WITH_MPI
3860 AssertDimension(requests.size(), tmp_data.size());
3861
3862 const unsigned int mf_component = find_vector_in_mf(vec);
3863
3864 const auto &part = get_partitioner(mf_component);
3865
3866 if (part.n_ghost_indices() > 0)
3867 {
3868 part.reset_ghost_values(ArrayView<Number>(
3869 const_cast<LinearAlgebra::distributed::Vector<Number> &>(vec)
3870 .begin() +
3871 part.locally_owned_size(),
3872 matrix_free.get_dof_info(mf_component)
3873 .vector_partitioner->n_ghost_indices()));
3874 }
3875
3876# endif
3877 }
3878 // let vector know that it's not ghosted anymore
3879 vec.set_ghost_state(false);
3880 }
3881
3882
3883
3889 template <typename VectorType,
3890 std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3891 * = nullptr>
3892 void
3893 zero_vector_region(const unsigned int range_index, VectorType &vec) const
3894 {
3895 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3896 "Type mismatch between VectorType and VectorDataExchange");
3897 if (range_index == numbers::invalid_unsigned_int)
3898 vec = Number();
3899 else
3900 {
3901 const unsigned int mf_component = find_vector_in_mf(vec, false);
3902 const internal::MatrixFreeFunctions::DoFInfo &dof_info =
3903 matrix_free.get_dof_info(mf_component);
3904 Assert(dof_info.vector_zero_range_list_index.empty() == false,
3905 ExcNotInitialized());
3906
3907 Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3908 ExcInternalError());
3909 AssertIndexRange(range_index,
3910 dof_info.vector_zero_range_list_index.size() - 1);
3911 for (unsigned int id =
3912 dof_info.vector_zero_range_list_index[range_index];
3913 id != dof_info.vector_zero_range_list_index[range_index + 1];
3914 ++id)
3915 std::memset(vec.begin() + dof_info.vector_zero_range_list[id].first,
3916 0,
3917 (dof_info.vector_zero_range_list[id].second -
3918 dof_info.vector_zero_range_list[id].first) *
3919 sizeof(Number));
3920 }
3921 }
3922
3923
3924
3930 template <typename VectorType,
3931 std::enable_if_t<!has_exchange_on_subset<VectorType>, VectorType>
3932 * = nullptr,
3933 typename VectorType::value_type * = nullptr>
3934 void
3935 zero_vector_region(const unsigned int range_index, VectorType &vec) const
3936 {
3937 if (range_index == numbers::invalid_unsigned_int || range_index == 0)
3938 {
3939 if constexpr (std::is_same_v<
3940 ArrayView<typename VectorType::value_type>,
3941 VectorType>)
3942 {
3943 for (unsigned int i = 0; i < vec.size(); ++i)
3944 vec[i] = typename VectorType::value_type();
3945 }
3946 else
3947 vec = typename VectorType::value_type();
3948 }
3949 }
3950
3951
3952
3957 void
3958 zero_vector_region(const unsigned int, ...) const
3959 {
3960 Assert(false,
3961 ExcNotImplemented("Zeroing is only implemented for vector types "
3962 "which provide VectorType::value_type"));
3963 }
3964
3965
3966
3967 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
3968 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3969 DataAccessOnFaces vector_face_access;
3970 bool ghosts_were_set;
3971# ifdef DEAL_II_WITH_MPI
3972 std::vector<AlignedVector<Number> *> tmp_data;
3973 std::vector<std::vector<MPI_Request>> requests;
3974# endif
3975 }; // VectorDataExchange
3976
3977 template <typename VectorStruct>
3978 unsigned int
3979 n_components(const VectorStruct &vec);
3980
3981 template <typename VectorStruct>
3982 unsigned int
3983 n_components_block(const VectorStruct &vec, const std::bool_constant<true>)
3984 {
3985 unsigned int components = 0;
3986 for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
3987 components += n_components(vec.block(bl));
3988 return components;
3989 }
3990
3991 template <typename VectorStruct>
3992 unsigned int
3993 n_components_block(const VectorStruct &, const std::bool_constant<false>)
3994 {
3995 return 1;
3996 }
3997
3998 template <typename VectorStruct>
3999 unsigned int
4000 n_components(const VectorStruct &vec)
4001 {
4002 return n_components_block(
4003 vec, std::bool_constant<IsBlockVector<VectorStruct>::value>());
4004 }
4005
4006 template <typename VectorStruct>
4007 inline unsigned int
4008 n_components(const std::vector<VectorStruct> &vec)
4009 {
4010 unsigned int components = 0;
4011 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4012 components += n_components_block(
4013 vec[comp], std::bool_constant<IsBlockVector<VectorStruct>::value>());
4014 return components;
4015 }
4016
4017 template <typename VectorStruct>
4018 inline unsigned int
4019 n_components(const std::vector<VectorStruct *> &vec)
4020 {
4021 unsigned int components = 0;
4022 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4023 components += n_components_block(
4024 *vec[comp], std::bool_constant<IsBlockVector<VectorStruct>::value>());
4025 return components;
4026 }
4027
4028
4029
4030 // A helper function to identify block vectors with many components where we
4031 // should not try to overlap computations and communication because there
4032 // would be too many outstanding communication requests.
4033
4034 // default value for vectors that do not have communication_block_size
4035 template <typename VectorStruct,
4036 std::enable_if_t<!has_communication_block_size<VectorStruct>,
4037 VectorStruct> * = nullptr>
4038 constexpr unsigned int
4039 get_communication_block_size(const VectorStruct &)
4040 {
4041 return numbers::invalid_unsigned_int;
4042 }
4043
4044
4045
4046 template <typename VectorStruct,
4047 std::enable_if_t<has_communication_block_size<VectorStruct>,
4048 VectorStruct> * = nullptr>
4049 constexpr unsigned int
4050 get_communication_block_size(const VectorStruct &)
4051 {
4052 return VectorStruct::communication_block_size;
4053 }
4054
4055
4056
4057 template <typename VectorType,
4058 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType> * =
4059 nullptr>
4060 bool
4061 has_ghost_elements(const VectorType &vec)
4062 {
4063 (void)vec;
4064 return false;
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 return vec.has_ghost_elements();
4076 }
4077
4078
4079
4080 // --------------------------------------------------------------------------
4081 // below we have wrappers to distinguish between block and non-block vectors.
4082 // --------------------------------------------------------------------------
4083
4084 //
4085 // update_ghost_values_start
4086 //
4087
4088 // update_ghost_values for block vectors
4089 template <int dim,
4090 typename VectorStruct,
4091 typename Number,
4092 typename VectorizedArrayType,
4093 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4094 * = nullptr>
4095 void
4096 update_ghost_values_start(
4097 const VectorStruct &vec,
4098 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4099 const unsigned int channel = 0)
4100 {
4101 if (get_communication_block_size(vec) < vec.n_blocks())
4102 {
4103 const bool ghosts_set = vec.has_ghost_elements();
4104
4105 Assert(exchanger.matrix_free.get_task_info()
4106 .allow_ghosted_vectors_in_loops ||
4107 ghosts_set == false,
4108 ExcNotImplemented());
4109
4110 if (ghosts_set)
4111 {
4112 exchanger.ghosts_were_set = true;
4113 return;
4114 }
4115
4116 vec.update_ghost_values();
4117 }
4118 else
4119 {
4120 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4121 update_ghost_values_start(vec.block(i), exchanger, channel + i);
4122 }
4123 }
4124
4125
4126
4127 // update_ghost_values for non-block vectors
4128 template <int dim,
4129 typename VectorStruct,
4130 typename Number,
4131 typename VectorizedArrayType,
4132 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4133 * = nullptr>
4134 void
4135 update_ghost_values_start(
4136 const VectorStruct &vec,
4137 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4138 const unsigned int channel = 0)
4139 {
4140 exchanger.update_ghost_values_start(channel, vec);
4141 }
4142
4143
4144
4145 // update_ghost_values_start() for vector of vectors
4146 template <int dim,
4147 typename VectorStruct,
4148 typename Number,
4149 typename VectorizedArrayType>
4150 inline void
4151 update_ghost_values_start(
4152 const std::vector<VectorStruct> &vec,
4153 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4154 {
4155 unsigned int component_index = 0;
4156 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4157 {
4158 update_ghost_values_start(vec[comp], exchanger, component_index);
4159 component_index += n_components(vec[comp]);
4160 }
4161 }
4162
4163
4164
4165 // update_ghost_values_start() for vector of pointers to vectors
4166 template <int dim,
4167 typename VectorStruct,
4168 typename Number,
4169 typename VectorizedArrayType>
4170 inline void
4171 update_ghost_values_start(
4172 const std::vector<VectorStruct *> &vec,
4173 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4174 {
4175 unsigned int component_index = 0;
4176 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4177 {
4178 update_ghost_values_start(*vec[comp], exchanger, component_index);
4179 component_index += n_components(*vec[comp]);
4180 }
4181 }
4182
4183
4184
4185 //
4186 // update_ghost_values_finish
4187 //
4188
4189 // for block vectors
4190 template <int dim,
4191 typename VectorStruct,
4192 typename Number,
4193 typename VectorizedArrayType,
4194 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4195 * = nullptr>
4196 void
4197 update_ghost_values_finish(
4198 const VectorStruct &vec,
4199 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4200 const unsigned int channel = 0)
4201 {
4202 if (get_communication_block_size(vec) < vec.n_blocks())
4203 {
4204 // do nothing, everything has already been completed in the _start()
4205 // call
4206 }
4207 else
4208 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4209 update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4210 }
4211
4212
4213
4214 // for non-block vectors
4215 template <int dim,
4216 typename VectorStruct,
4217 typename Number,
4218 typename VectorizedArrayType,
4219 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4220 * = nullptr>
4221 void
4222 update_ghost_values_finish(
4223 const VectorStruct &vec,
4224 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4225 const unsigned int channel = 0)
4226 {
4227 exchanger.update_ghost_values_finish(channel, vec);
4228 }
4229
4230
4231
4232 // for vector of vectors
4233 template <int dim,
4234 typename VectorStruct,
4235 typename Number,
4236 typename VectorizedArrayType>
4237 inline void
4238 update_ghost_values_finish(
4239 const std::vector<VectorStruct> &vec,
4240 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4241 {
4242 unsigned int component_index = 0;
4243 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4244 {
4245 update_ghost_values_finish(vec[comp], exchanger, component_index);
4246 component_index += n_components(vec[comp]);
4247 }
4248 }
4249
4250
4251
4252 // for vector of pointers to vectors
4253 template <int dim,
4254 typename VectorStruct,
4255 typename Number,
4256 typename VectorizedArrayType>
4257 inline void
4258 update_ghost_values_finish(
4259 const std::vector<VectorStruct *> &vec,
4260 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4261 {
4262 unsigned int component_index = 0;
4263 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4264 {
4265 update_ghost_values_finish(*vec[comp], exchanger, component_index);
4266 component_index += n_components(*vec[comp]);
4267 }
4268 }
4269
4270
4271
4272 //
4273 // compress_start
4274 //
4275
4276 // for block vectors
4277 template <int dim,
4278 typename VectorStruct,
4279 typename Number,
4280 typename VectorizedArrayType,
4281 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4282 * = nullptr>
4283 inline void
4284 compress_start(
4285 VectorStruct &vec,
4286 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4287 const unsigned int channel = 0)
4288 {
4289 if (get_communication_block_size(vec) < vec.n_blocks())
4290 vec.compress(VectorOperation::add);
4291 else
4292 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4293 compress_start(vec.block(i), exchanger, channel + i);
4294 }
4295
4296
4297
4298 // for non-block vectors
4299 template <int dim,
4300 typename VectorStruct,
4301 typename Number,
4302 typename VectorizedArrayType,
4303 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4304 * = nullptr>
4305 inline void
4306 compress_start(
4307 VectorStruct &vec,
4308 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4309 const unsigned int channel = 0)
4310 {
4311 exchanger.compress_start(channel, vec);
4312 }
4313
4314
4315
4316 // for std::vector of vectors
4317 template <int dim,
4318 typename VectorStruct,
4319 typename Number,
4320 typename VectorizedArrayType>
4321 inline void
4322 compress_start(
4323 std::vector<VectorStruct> &vec,
4324 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4325 {
4326 unsigned int component_index = 0;
4327 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4328 {
4329 compress_start(vec[comp], exchanger, component_index);
4330 component_index += n_components(vec[comp]);
4331 }
4332 }
4333
4334
4335
4336 // for std::vector of pointer to vectors
4337 template <int dim,
4338 typename VectorStruct,
4339 typename Number,
4340 typename VectorizedArrayType>
4341 inline void
4342 compress_start(
4343 std::vector<VectorStruct *> &vec,
4344 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4345 {
4346 unsigned int component_index = 0;
4347 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4348 {
4349 compress_start(*vec[comp], exchanger, component_index);
4350 component_index += n_components(*vec[comp]);
4351 }
4352 }
4353
4354
4355
4356 //
4357 // compress_finish
4358 //
4359
4360 // for block vectors
4361 template <int dim,
4362 typename VectorStruct,
4363 typename Number,
4364 typename VectorizedArrayType,
4365 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4366 * = nullptr>
4367 inline void
4368 compress_finish(
4369 VectorStruct &vec,
4370 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4371 const unsigned int channel = 0)
4372 {
4373 if (get_communication_block_size(vec) < vec.n_blocks())
4374 {
4375 // do nothing, everything has already been completed in the _start()
4376 // call
4377 }
4378 else
4379 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4380 compress_finish(vec.block(i), exchanger, channel + i);
4381 }
4382
4383
4384
4385 // for non-block vectors
4386 template <int dim,
4387 typename VectorStruct,
4388 typename Number,
4389 typename VectorizedArrayType,
4390 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4391 * = nullptr>
4392 inline void
4393 compress_finish(
4394 VectorStruct &vec,
4395 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4396 const unsigned int channel = 0)
4397 {
4398 exchanger.compress_finish(channel, vec);
4399 }
4400
4401
4402
4403 // for std::vector of vectors
4404 template <int dim,
4405 typename VectorStruct,
4406 typename Number,
4407 typename VectorizedArrayType>
4408 inline void
4409 compress_finish(
4410 std::vector<VectorStruct> &vec,
4411 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4412 {
4413 unsigned int component_index = 0;
4414 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4415 {
4416 compress_finish(vec[comp], exchanger, component_index);
4417 component_index += n_components(vec[comp]);
4418 }
4419 }
4420
4421
4422
4423 // for std::vector of pointer to vectors
4424 template <int dim,
4425 typename VectorStruct,
4426 typename Number,
4427 typename VectorizedArrayType>
4428 inline void
4429 compress_finish(
4430 std::vector<VectorStruct *> &vec,
4431 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4432 {
4433 unsigned int component_index = 0;
4434 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4435 {
4436 compress_finish(*vec[comp], exchanger, component_index);
4437 component_index += n_components(*vec[comp]);
4438 }
4439 }
4440
4441
4442
4443 //
4444 // reset_ghost_values:
4445 //
4446 // if the input vector did not have ghosts imported, clear them here again
4447 // in order to avoid subsequent operations e.g. in linear solvers to work
4448 // with ghosts all the time
4449 //
4450
4451 // for block vectors
4452 template <int dim,
4453 typename VectorStruct,
4454 typename Number,
4455 typename VectorizedArrayType,
4456 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4457 * = nullptr>
4458 inline void
4459 reset_ghost_values(
4460 const VectorStruct &vec,
4461 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4462 {
4463 // return immediately if there is nothing to do.
4464 if (exchanger.ghosts_were_set == true)
4465 return;
4466
4467 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4468 reset_ghost_values(vec.block(i), exchanger);
4469 }
4470
4471
4472
4473 // for non-block vectors
4474 template <int dim,
4475 typename VectorStruct,
4476 typename Number,
4477 typename VectorizedArrayType,
4478 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4479 * = nullptr>
4480 inline void
4481 reset_ghost_values(
4482 const VectorStruct &vec,
4483 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4484 {
4485 // return immediately if there is nothing to do.
4486 if (exchanger.ghosts_were_set == true)
4487 return;
4488
4489 exchanger.reset_ghost_values(vec);
4490 }
4491
4492
4493
4494 // for std::vector of vectors
4495 template <int dim,
4496 typename VectorStruct,
4497 typename Number,
4498 typename VectorizedArrayType>
4499 inline void
4500 reset_ghost_values(
4501 const std::vector<VectorStruct> &vec,
4502 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4503 {
4504 // return immediately if there is nothing to do.
4505 if (exchanger.ghosts_were_set == true)
4506 return;
4507
4508 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4509 reset_ghost_values(vec[comp], exchanger);
4510 }
4511
4512
4513
4514 // for std::vector of pointer to vectors
4515 template <int dim,
4516 typename VectorStruct,
4517 typename Number,
4518 typename VectorizedArrayType>
4519 inline void
4520 reset_ghost_values(
4521 const std::vector<VectorStruct *> &vec,
4522 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4523 {
4524 // return immediately if there is nothing to do.
4525 if (exchanger.ghosts_were_set == true)
4526 return;
4527
4528 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4529 reset_ghost_values(*vec[comp], exchanger);
4530 }
4531
4532
4533
4534 //
4535 // zero_vector_region
4536 //
4537
4538 // for block vectors
4539 template <int dim,
4540 typename VectorStruct,
4541 typename Number,
4542 typename VectorizedArrayType,
4543 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4544 * = nullptr>
4545 inline void
4546 zero_vector_region(
4547 const unsigned int range_index,
4548 VectorStruct &vec,
4549 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4550 {
4551 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4552 exchanger.zero_vector_region(range_index, vec.block(i));
4553 }
4554
4555
4556
4557 // for non-block vectors
4558 template <int dim,
4559 typename VectorStruct,
4560 typename Number,
4561 typename VectorizedArrayType,
4562 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4563 * = nullptr>
4564 inline void
4565 zero_vector_region(
4566 const unsigned int range_index,
4567 VectorStruct &vec,
4568 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4569 {
4570 exchanger.zero_vector_region(range_index, vec);
4571 }
4572
4573
4574
4575 // for std::vector of vectors
4576 template <int dim,
4577 typename VectorStruct,
4578 typename Number,
4579 typename VectorizedArrayType>
4580 inline void
4581 zero_vector_region(
4582 const unsigned int range_index,
4583 std::vector<VectorStruct> &vec,
4584 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4585 {
4586 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4587 zero_vector_region(range_index, vec[comp], exchanger);
4588 }
4589
4590
4591
4592 // for std::vector of pointers to vectors
4593 template <int dim,
4594 typename VectorStruct,
4595 typename Number,
4596 typename VectorizedArrayType>
4597 inline void
4598 zero_vector_region(
4599 const unsigned int range_index,
4600 std::vector<VectorStruct *> &vec,
4601 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4602 {
4603 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4604 zero_vector_region(range_index, *vec[comp], exchanger);
4605 }
4606
4607
4608
4609 // Apply a unit matrix operation to constrained DoFs: Default cases where we
4610 // cannot detect a LinearAlgebra::distributed::Vector, we do not do
4611 // anything, else we apply the constraints as a unit operation
4612 template <typename VectorStruct1, typename VectorStruct2>
4613 inline void
4614 apply_operation_to_constrained_dofs(const std::vector<unsigned int> &,
4615 const VectorStruct1 &,
4616 VectorStruct2 &)
4617 {}
4618
4619 template <typename Number>
4620 inline void
4621 apply_operation_to_constrained_dofs(
4622 const std::vector<unsigned int> &constrained_dofs,
4623 const LinearAlgebra::distributed::Vector<Number> &src,
4624 LinearAlgebra::distributed::Vector<Number> &dst)
4625 {
4626 for (const unsigned int i : constrained_dofs)
4627 dst.local_element(i) = src.local_element(i);
4628 }
4629
4630
4631 namespace MatrixFreeFunctions
4632 {
4633 // struct to select between a const interface and a non-const interface
4634 // for MFWorker
4635 template <typename, typename, typename, typename, bool>
4636 struct InterfaceSelector
4637 {};
4638
4639 // Version for constant functions
4640 template <typename MF,
4641 typename InVector,
4642 typename OutVector,
4643 typename Container>
4644 struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4645 {
4646 using function_type = void (Container::*)(
4647 const MF &,
4648 OutVector &,
4649 const InVector &,
4650 const std::pair<unsigned int, unsigned int> &) const;
4651 };
4652
4653 // Version for non-constant functions
4654 template <typename MF,
4655 typename InVector,
4656 typename OutVector,
4657 typename Container>
4658 struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4659 {
4660 using function_type =
4661 void (Container::*)(const MF &,
4662 OutVector &,
4663 const InVector &,
4664 const std::pair<unsigned int, unsigned int> &);
4665 };
4666 } // namespace MatrixFreeFunctions
4667
4668
4669
4670 // A implementation class for the worker object that runs the various
4671 // operations we want to perform during the matrix-free loop
4672 template <typename MF,
4673 typename InVector,
4674 typename OutVector,
4675 typename Container,
4676 bool is_constant>
4677 class MFWorker : public MFWorkerInterface
4678 {
4679 public:
4680 // An alias to make the arguments further down more readable
4681 using function_type = typename MatrixFreeFunctions::
4682 InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4683 function_type;
4684
4685 // constructor, binds all the arguments to this class
4686 MFWorker(const MF &matrix_free,
4687 const InVector &src,
4688 OutVector &dst,
4689 const bool zero_dst_vector_setting,
4690 const Container &container,
4691 function_type cell_function,
4692 function_type face_function,
4693 function_type boundary_function,
4694 const typename MF::DataAccessOnFaces src_vector_face_access =
4695 MF::DataAccessOnFaces::none,
4696 const typename MF::DataAccessOnFaces dst_vector_face_access =
4697 MF::DataAccessOnFaces::none,
4698 const std::function<void(const unsigned int, const unsigned int)>
4699 &operation_before_loop = {},
4700 const std::function<void(const unsigned int, const unsigned int)>
4701 &operation_after_loop = {},
4702 const unsigned int dof_handler_index_pre_post = 0)
4703 : matrix_free(matrix_free)
4704 , container(const_cast<Container &>(container))
4705 , cell_function(cell_function)
4706 , face_function(face_function)
4707 , boundary_function(boundary_function)
4708 , src(src)
4709 , dst(dst)
4710 , src_data_exchanger(matrix_free,
4711 src_vector_face_access,
4712 n_components(src))
4713 , dst_data_exchanger(matrix_free,
4714 dst_vector_face_access,
4715 n_components(dst))
4716 , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4717 , zero_dst_vector_setting(zero_dst_vector_setting &&
4718 !src_and_dst_are_same)
4719 , operation_before_loop(operation_before_loop)
4720 , operation_after_loop(operation_after_loop)
4721 , dof_handler_index_pre_post(dof_handler_index_pre_post)
4722 {
4723 Assert(!has_ghost_elements(dst),
4724 ExcMessage("The destination vector passed to the matrix-free "
4725 "loop is ghosted. This is not allowed."));
4726 }
4727
4728 // Runs the cell work. If no function is given, nothing is done
4729 virtual void
4730 cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4731 {
4732 if (cell_function != nullptr && cell_range.second > cell_range.first)
4733 for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4734 {
4735 const auto cell_subrange =
4736 matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4737
4738 if (cell_subrange.second <= cell_subrange.first)
4739 continue;
4740
4741 (container.*
4742 cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4743 }
4744 }
4745
4746 virtual void
4747 cell(const unsigned int range_index) override
4748 {
4749 process_range(cell_function,
4750 matrix_free.get_task_info().cell_partition_data_hp_ptr,
4751 matrix_free.get_task_info().cell_partition_data_hp,
4752 range_index);
4753 }
4754
4755 virtual void
4756 face(const unsigned int range_index) override
4757 {
4758 process_range(face_function,
4759 matrix_free.get_task_info().face_partition_data_hp_ptr,
4760 matrix_free.get_task_info().face_partition_data_hp,
4761 range_index);
4762 }
4763
4764 virtual void
4765 boundary(const unsigned int range_index) override
4766 {
4767 process_range(boundary_function,
4768 matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4769 matrix_free.get_task_info().boundary_partition_data_hp,
4770 range_index);
4771 }
4772
4773 private:
4774 void
4775 process_range(const function_type &fu,
4776 const std::vector<unsigned int> &ptr,
4777 const std::vector<unsigned int> &data,
4778 const unsigned int range_index)
4779 {
4780 if (fu == nullptr)
4781 return;
4782
4783 AssertIndexRange(range_index + 1, ptr.size());
4784 for (unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4785 {
4786 AssertIndexRange(2 * i + 1, data.size());
4787 (container.*fu)(matrix_free,
4788 this->dst,
4789 this->src,
4790 std::make_pair(data[2 * i], data[2 * i + 1]));
4791 }
4792 }
4793
4794 public:
4795 // Starts the communication for the update ghost values operation. We
4796 // cannot call this update if ghost and destination are the same because
4797 // that would introduce spurious entries in the destination (there is also
4798 // the problem that reading from a vector that we also write to is usually
4799 // not intended in case there is overlap, but this is up to the
4800 // application code to decide and we cannot catch this case here).
4801 virtual void
4802 vector_update_ghosts_start() override
4803 {
4804 if (!src_and_dst_are_same)
4805 internal::update_ghost_values_start(src, src_data_exchanger);
4806 }
4807
4808 // Finishes the communication for the update ghost values operation
4809 virtual void
4810 vector_update_ghosts_finish() override
4811 {
4812 if (!src_and_dst_are_same)
4813 internal::update_ghost_values_finish(src, src_data_exchanger);
4814 }
4815
4816 // Starts the communication for the vector compress operation
4817 virtual void
4818 vector_compress_start() override
4819 {
4820 internal::compress_start(dst, dst_data_exchanger);
4821 }
4822
4823 // Finishes the communication for the vector compress operation
4824 virtual void
4825 vector_compress_finish() override
4826 {
4827 internal::compress_finish(dst, dst_data_exchanger);
4828 if (!src_and_dst_are_same)
4829 internal::reset_ghost_values(src, src_data_exchanger);
4830 }
4831
4832 // Zeros the given input vector
4833 virtual void
4834 zero_dst_vector_range(const unsigned int range_index) override
4835 {
4836 if (zero_dst_vector_setting)
4837 internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4838 }
4839
4840 virtual void
4841 cell_loop_pre_range(const unsigned int range_index) override
4842 {
4843 if (operation_before_loop)
4844 {
4845 const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4846 matrix_free.get_dof_info(dof_handler_index_pre_post);
4847 if (range_index == numbers::invalid_unsigned_int)
4848 {
4849 // Case with threaded loop -> currently no overlap implemented
4850 ::parallel::apply_to_subranges(
4851 0U,
4852 dof_info.vector_partitioner->locally_owned_size(),
4853 operation_before_loop,
4854 internal::VectorImplementation::minimum_parallel_grain_size);
4855 }
4856 else
4857 {
4858 AssertIndexRange(range_index,
4859 dof_info.cell_loop_pre_list_index.size() - 1);
4860 for (unsigned int id =
4861 dof_info.cell_loop_pre_list_index[range_index];
4862 id != dof_info.cell_loop_pre_list_index[range_index + 1];
4863 ++id)
4864 operation_before_loop(dof_info.cell_loop_pre_list[id].first,
4865 dof_info.cell_loop_pre_list[id].second);
4866 }
4867 }
4868 }
4869
4870 virtual void
4871 cell_loop_post_range(const unsigned int range_index) override
4872 {
4873 if (operation_after_loop)
4874 {
4875 // Run unit matrix operation on constrained dofs if we are at the
4876 // last range
4877 const std::vector<unsigned int> &partition_row_index =
4878 matrix_free.get_task_info().partition_row_index;
4879 if (range_index ==
4880 partition_row_index[partition_row_index.size() - 2] - 1)
4881 apply_operation_to_constrained_dofs(
4882 matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
4883 src,
4884 dst);
4885
4886 const internal::MatrixFreeFunctions::DoFInfo &dof_info =
4887 matrix_free.get_dof_info(dof_handler_index_pre_post);
4888 if (range_index == numbers::invalid_unsigned_int)
4889 {
4890 // Case with threaded loop -> currently no overlap implemented
4891 ::parallel::apply_to_subranges(
4892 0U,
4893 dof_info.vector_partitioner->locally_owned_size(),
4894 operation_after_loop,
4895 internal::VectorImplementation::minimum_parallel_grain_size);
4896 }
4897 else
4898 {
4899 AssertIndexRange(range_index,
4900 dof_info.cell_loop_post_list_index.size() - 1);
4901 for (unsigned int id =
4902 dof_info.cell_loop_post_list_index[range_index];
4903 id != dof_info.cell_loop_post_list_index[range_index + 1];
4904 ++id)
4905 operation_after_loop(dof_info.cell_loop_post_list[id].first,
4906 dof_info.cell_loop_post_list[id].second);
4907 }
4908 }
4909 }
4910
4911 private:
4912 const MF &matrix_free;
4913 Container &container;
4914 function_type cell_function;
4915 function_type face_function;
4916 function_type boundary_function;
4917
4918 const InVector &src;
4919 OutVector &dst;
4920 VectorDataExchange<MF::dimension,
4921 typename MF::value_type,
4922 typename MF::vectorized_value_type>
4923 src_data_exchanger;
4924 VectorDataExchange<MF::dimension,
4925 typename MF::value_type,
4926 typename MF::vectorized_value_type>
4927 dst_data_exchanger;
4928 const bool src_and_dst_are_same;
4929 const bool zero_dst_vector_setting;
4930 const std::function<void(const unsigned int, const unsigned int)>
4931 operation_before_loop;
4932 const std::function<void(const unsigned int, const unsigned int)>
4933 operation_after_loop;
4934 const unsigned int dof_handler_index_pre_post;
4935 };
4936
4937
4938
4943 template <class MF, typename InVector, typename OutVector>
4944 struct MFClassWrapper
4945 {
4946 using function_type =
4947 std::function<void(const MF &,
4948 OutVector &,
4949 const InVector &,
4950 const std::pair<unsigned int, unsigned int> &)>;
4951
4952 MFClassWrapper(const function_type cell,
4953 const function_type face,
4954 const function_type boundary)
4955 : cell(cell)
4956 , face(face)
4957 , boundary(boundary)
4958 {}
4959
4960 void
4961 cell_integrator(const MF &mf,
4962 OutVector &dst,
4963 const InVector &src,
4964 const std::pair<unsigned int, unsigned int> &range) const
4965 {
4966 if (cell)
4967 cell(mf, dst, src, range);
4968 }
4969
4970 void
4971 face_integrator(const MF &mf,
4972 OutVector &dst,
4973 const InVector &src,
4974 const std::pair<unsigned int, unsigned int> &range) const
4975 {
4976 if (face)
4977 face(mf, dst, src, range);
4978 }
4979
4980 void
4981 boundary_integrator(
4982 const MF &mf,
4983 OutVector &dst,
4984 const InVector &src,
4985 const std::pair<unsigned int, unsigned int> &range) const
4986 {
4987 if (boundary)
4988 boundary(mf, dst, src, range);
4989 }
4990
4991 const function_type cell;
4992 const function_type face;
4993 const function_type boundary;
4994 };
4995
4996} // end of namespace internal
4997
4998
4999
5000template <int dim, typename Number, typename VectorizedArrayType>
5001template <typename OutVector, typename InVector>
5002inline void
5003MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5004 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5005 OutVector &,
5006 const InVector &,
5007 const std::pair<unsigned int, unsigned int> &)>
5008 &cell_operation,
5009 OutVector &dst,
5010 const InVector &src,
5011 const bool zero_dst_vector) const
5012{
5013 using Wrapper =
5014 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5015 InVector,
5016 OutVector>;
5017 Wrapper wrap(cell_operation, nullptr, nullptr);
5018 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5019 InVector,
5020 OutVector,
5021 Wrapper,
5022 true>
5023 worker(*this,
5024 src,
5025 dst,
5026 zero_dst_vector,
5027 wrap,
5028 &Wrapper::cell_integrator,
5029 &Wrapper::face_integrator,
5030 &Wrapper::boundary_integrator);
5031
5032 task_info.loop(worker);
5033}
5034
5035
5036
5037template <int dim, typename Number, typename VectorizedArrayType>
5038template <typename OutVector, typename InVector>
5039inline void
5040MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5041 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5042 OutVector &,
5043 const InVector &,
5044 const std::pair<unsigned int, unsigned int> &)>
5045 &cell_operation,
5046 OutVector &dst,
5047 const InVector &src,
5048 const std::function<void(const unsigned int, const unsigned int)>
5049 &operation_before_loop,
5050 const std::function<void(const unsigned int, const unsigned int)>
5051 &operation_after_loop,
5052 const unsigned int dof_handler_index_pre_post) const
5053{
5054 using Wrapper =
5055 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5056 InVector,
5057 OutVector>;
5058 Wrapper wrap(cell_operation, nullptr, nullptr);
5059 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5060 InVector,
5061 OutVector,
5062 Wrapper,
5063 true>
5064 worker(*this,
5065 src,
5066 dst,
5067 false,
5068 wrap,
5069 &Wrapper::cell_integrator,
5070 &Wrapper::face_integrator,
5071 &Wrapper::boundary_integrator,
5072 DataAccessOnFaces::none,
5073 DataAccessOnFaces::none,
5074 operation_before_loop,
5075 operation_after_loop,
5076 dof_handler_index_pre_post);
5077
5078 task_info.loop(worker);
5079}
5080
5081
5082
5083template <int dim, typename Number, typename VectorizedArrayType>
5084template <typename OutVector, typename InVector>
5085inline void
5086MatrixFree<dim, Number, VectorizedArrayType>::loop(
5087 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5088 OutVector &,
5089 const InVector &,
5090 const std::pair<unsigned int, unsigned int> &)>
5091 &cell_operation,
5092 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5093 OutVector &,
5094 const InVector &,
5095 const std::pair<unsigned int, unsigned int> &)>
5096 &inner_face_operation,
5097 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5098 OutVector &,
5099 const InVector &,
5100 const std::pair<unsigned int, unsigned int> &)>
5101 &boundary_face_operation,
5102 OutVector &dst,
5103 const InVector &src,
5104 const bool zero_dst_vector,
5105 const DataAccessOnFaces dst_vector_face_access,
5106 const DataAccessOnFaces src_vector_face_access) const
5107{
5108 using Wrapper =
5109 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5110 InVector,
5111 OutVector>;
5112 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5113 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5114 InVector,
5115 OutVector,
5116 Wrapper,
5117 true>
5118 worker(*this,
5119 src,
5120 dst,
5121 zero_dst_vector,
5122 wrap,
5123 &Wrapper::cell_integrator,
5124 &Wrapper::face_integrator,
5125 &Wrapper::boundary_integrator,
5126 src_vector_face_access,
5127 dst_vector_face_access);
5128
5129 task_info.loop(worker);
5130}
5131
5132
5133
5134template <int dim, typename Number, typename VectorizedArrayType>
5135template <typename CLASS, typename OutVector, typename InVector>
5136inline void
5137MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5138 void (CLASS::*function_pointer)(
5139 const MatrixFree<dim, Number, VectorizedArrayType> &,
5140 OutVector &,
5141 const InVector &,
5142 const std::pair<unsigned int, unsigned int> &) const,
5143 const CLASS *owning_class,
5144 OutVector &dst,
5145 const InVector &src,
5146 const bool zero_dst_vector) const
5147{
5148 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5149 InVector,
5150 OutVector,
5151 CLASS,
5152 true>
5153 worker(*this,
5154 src,
5155 dst,
5156 zero_dst_vector,
5157 *owning_class,
5158 function_pointer,
5159 nullptr,
5160 nullptr);
5161 task_info.loop(worker);
5162}
5163
5164
5165
5166template <int dim, typename Number, typename VectorizedArrayType>
5167template <typename CLASS, typename OutVector, typename InVector>
5168inline void
5169MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5170 void (CLASS::*function_pointer)(
5171 const MatrixFree<dim, Number, VectorizedArrayType> &,
5172 OutVector &,
5173 const InVector &,
5174 const std::pair<unsigned int, unsigned int> &) const,
5175 const CLASS *owning_class,
5176 OutVector &dst,
5177 const InVector &src,
5178 const std::function<void(const unsigned int, const unsigned int)>
5179 &operation_before_loop,
5180 const std::function<void(const unsigned int, const unsigned int)>
5181 &operation_after_loop,
5182 const unsigned int dof_handler_index_pre_post) const
5183{
5184 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5185 InVector,
5186 OutVector,
5187 CLASS,
5188 true>
5189 worker(*this,
5190 src,
5191 dst,
5192 false,
5193 *owning_class,
5194 function_pointer,
5195 nullptr,
5196 nullptr,
5197 DataAccessOnFaces::none,
5198 DataAccessOnFaces::none,
5199 operation_before_loop,
5200 operation_after_loop,
5201 dof_handler_index_pre_post);
5202 task_info.loop(worker);
5203}
5204
5205
5206
5207template <int dim, typename Number, typename VectorizedArrayType>
5208template <typename CLASS, typename OutVector, typename InVector>
5209inline void
5210MatrixFree<dim, Number, VectorizedArrayType>::loop(
5211 void (CLASS::*cell_operation)(
5212 const MatrixFree<dim, Number, VectorizedArrayType> &,
5213 OutVector &,
5214 const InVector &,
5215 const std::pair<unsigned int, unsigned int> &) const,
5216 void (CLASS::*inner_face_operation)(
5217 const MatrixFree<dim, Number, VectorizedArrayType> &,
5218 OutVector &,
5219 const InVector &,
5220 const std::pair<unsigned int, unsigned int> &) const,
5221 void (CLASS::*boundary_face_operation)(
5222 const MatrixFree<dim, Number, VectorizedArrayType> &,
5223 OutVector &,
5224 const InVector &,
5225 const std::pair<unsigned int, unsigned int> &) const,
5226 const CLASS *owning_class,
5227 OutVector &dst,
5228 const InVector &src,
5229 const bool zero_dst_vector,
5230 const DataAccessOnFaces dst_vector_face_access,
5231 const DataAccessOnFaces src_vector_face_access) const
5232{
5233 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5234 InVector,
5235 OutVector,
5236 CLASS,
5237 true>
5238 worker(*this,
5239 src,
5240 dst,
5241 zero_dst_vector,
5242 *owning_class,
5243 cell_operation,
5244 inner_face_operation,
5245 boundary_face_operation,
5246 src_vector_face_access,
5247 dst_vector_face_access);
5248 task_info.loop(worker);
5249}
5250
5251
5252
5253template <int dim, typename Number, typename VectorizedArrayType>
5254template <typename CLASS, typename OutVector, typename InVector>
5255inline void
5256MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5257 void (CLASS::*function_pointer)(
5258 const MatrixFree<dim, Number, VectorizedArrayType> &,
5259 OutVector &,
5260 const InVector &,
5261 const std::pair<unsigned int, unsigned int> &),
5262 CLASS *owning_class,
5263 OutVector &dst,
5264 const InVector &src,
5265 const bool zero_dst_vector) const
5266{
5267 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5268 InVector,
5269 OutVector,
5270 CLASS,
5271 false>
5272 worker(*this,
5273 src,
5274 dst,
5275 zero_dst_vector,
5276 *owning_class,
5277 function_pointer,
5278 nullptr,
5279 nullptr);
5280 task_info.loop(worker);
5281}
5282
5283
5284
5285template <int dim, typename Number, typename VectorizedArrayType>
5286template <typename CLASS, typename OutVector, typename InVector>
5287inline void
5288MatrixFree<dim, Number, VectorizedArrayType>::cell_loop(
5289 void (CLASS::*function_pointer)(
5290 const MatrixFree<dim, Number, VectorizedArrayType> &,
5291 OutVector &,
5292 const InVector &,
5293 const std::pair<unsigned int, unsigned int> &),
5294 CLASS *owning_class,
5295 OutVector &dst,
5296 const InVector &src,
5297 const std::function<void(const unsigned int, const unsigned int)>
5298 &operation_before_loop,
5299 const std::function<void(const unsigned int, const unsigned int)>
5300 &operation_after_loop,
5301 const unsigned int dof_handler_index_pre_post) const
5302{
5303 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5304 InVector,
5305 OutVector,
5306 CLASS,
5307 false>
5308 worker(*this,
5309 src,
5310 dst,
5311 false,
5312 *owning_class,
5313 function_pointer,
5314 nullptr,
5315 nullptr,
5316 DataAccessOnFaces::none,
5317 DataAccessOnFaces::none,
5318 operation_before_loop,
5319 operation_after_loop,
5320 dof_handler_index_pre_post);
5321 task_info.loop(worker);
5322}
5323
5324
5325
5326template <int dim, typename Number, typename VectorizedArrayType>
5327template <typename CLASS, typename OutVector, typename InVector>
5328inline void
5329MatrixFree<dim, Number, VectorizedArrayType>::loop(
5330 void (CLASS::*cell_operation)(
5331 const MatrixFree<dim, Number, VectorizedArrayType> &,
5332 OutVector &,
5333 const InVector &,
5334 const std::pair<unsigned int, unsigned int> &),
5335 void (CLASS::*inner_face_operation)(
5336 const MatrixFree<dim, Number, VectorizedArrayType> &,
5337 OutVector &,
5338 const InVector &,
5339 const std::pair<unsigned int, unsigned int> &),
5340 void (CLASS::*boundary_face_operation)(
5341 const MatrixFree<dim, Number, VectorizedArrayType> &,
5342 OutVector &,
5343 const InVector &,
5344 const std::pair<unsigned int, unsigned int> &),
5345 CLASS *owning_class,
5346 OutVector &dst,
5347 const InVector &src,
5348 const bool zero_dst_vector,
5349 const DataAccessOnFaces dst_vector_face_access,
5350 const DataAccessOnFaces src_vector_face_access) const
5351{
5352 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5353 InVector,
5354 OutVector,
5355 CLASS,
5356 false>
5357 worker(*this,
5358 src,
5359 dst,
5360 zero_dst_vector,
5361 *owning_class,
5362 cell_operation,
5363 inner_face_operation,
5364 boundary_face_operation,
5365 src_vector_face_access,
5366 dst_vector_face_access);
5367 task_info.loop(worker);
5368}
5369
5370
5371
5372template <int dim, typename Number, typename VectorizedArrayType>
5373template <typename OutVector, typename InVector>
5374inline void
5375MatrixFree<dim, Number, VectorizedArrayType>::loop(
5376 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5377 OutVector &,
5378 const InVector &,
5379 const std::pair<unsigned int, unsigned int> &)>
5380 &cell_operation,
5381 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5382 OutVector &,
5383 const InVector &,
5384 const std::pair<unsigned int, unsigned int> &)>
5385 &inner_face_operation,
5386 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5387 OutVector &,
5388 const InVector &,
5389 const std::pair<unsigned int, unsigned int> &)>
5390 &boundary_face_operation,
5391 OutVector &dst,
5392 const InVector &src,
5393 const std::function<void(const unsigned int, const unsigned int)>
5394 &operation_before_loop,
5395 const std::function<void(const unsigned int, const unsigned int)>
5396 &operation_after_loop,
5397 const unsigned int dof_handler_index_pre_post,
5398 const DataAccessOnFaces dst_vector_face_access,
5399 const DataAccessOnFaces src_vector_face_access) const
5400{
5401 using Wrapper =
5402 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5403 InVector,
5404 OutVector>;
5405 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5406 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5407 InVector,
5408 OutVector,
5409 Wrapper,
5410 true>
5411 worker(*this,
5412 src,
5413 dst,
5414 false,
5415 wrap,
5416 &Wrapper::cell_integrator,
5417 &Wrapper::face_integrator,
5418 &Wrapper::boundary_integrator,
5419 src_vector_face_access,
5420 dst_vector_face_access,
5421 operation_before_loop,
5422 operation_after_loop,
5423 dof_handler_index_pre_post);
5424
5425 task_info.loop(worker);
5426}
5427
5428
5429
5430template <int dim, typename Number, typename VectorizedArrayType>
5431template <typename CLASS, typename OutVector, typename InVector>
5432inline void
5433MatrixFree<dim, Number, VectorizedArrayType>::loop(
5434 void (CLASS::*cell_operation)(const MatrixFree &,
5435 OutVector &,
5436 const InVector &,
5437 const std::pair<unsigned int, unsigned int> &)
5438 const,
5439 void (CLASS::*inner_face_operation)(
5440 const MatrixFree &,
5441 OutVector &,
5442 const InVector &,
5443 const std::pair<unsigned int, unsigned int> &) const,
5444 void (CLASS::*boundary_face_operation)(
5445 const MatrixFree &,
5446 OutVector &,
5447 const InVector &,
5448 const std::pair<unsigned int, unsigned int> &) const,
5449 const CLASS *owning_class,
5450 OutVector &dst,
5451 const InVector &src,
5452 const std::function<void(const unsigned int, const unsigned int)>
5453 &operation_before_loop,
5454 const std::function<void(const unsigned int, const unsigned int)>
5455 &operation_after_loop,
5456 const unsigned int dof_handler_index_pre_post,
5457 const DataAccessOnFaces dst_vector_face_access,
5458 const DataAccessOnFaces src_vector_face_access) const
5459{
5460 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5461 InVector,
5462 OutVector,
5463 CLASS,
5464 true>
5465 worker(*this,
5466 src,
5467 dst,
5468 false,
5469 *owning_class,
5470 cell_operation,
5471 inner_face_operation,
5472 boundary_face_operation,
5473 src_vector_face_access,
5474 dst_vector_face_access,
5475 operation_before_loop,
5476 operation_after_loop,
5477 dof_handler_index_pre_post);
5478 task_info.loop(worker);
5479}
5480
5481
5482
5483template <int dim, typename Number, typename VectorizedArrayType>
5484template <typename CLASS, typename OutVector, typename InVector>
5485inline void
5486MatrixFree<dim, Number, VectorizedArrayType>::loop(
5487 void (CLASS::*cell_operation)(const MatrixFree &,
5488 OutVector &,
5489 const InVector &,
5490 const std::pair<unsigned int, unsigned int> &),
5491 void (CLASS::*inner_face_operation)(
5492 const MatrixFree &,
5493 OutVector &,
5494 const InVector &,
5495 const std::pair<unsigned int, unsigned int> &),
5496 void (CLASS::*boundary_face_operation)(
5497 const MatrixFree &,
5498 OutVector &,
5499 const InVector &,
5500 const std::pair<unsigned int, unsigned int> &),
5501 const CLASS *owning_class,
5502 OutVector &dst,
5503 const InVector &src,
5504 const std::function<void(const unsigned int, const unsigned int)>
5505 &operation_before_loop,
5506 const std::function<void(const unsigned int, const unsigned int)>
5507 &operation_after_loop,
5508 const unsigned int dof_handler_index_pre_post,
5509 const DataAccessOnFaces dst_vector_face_access,
5510 const DataAccessOnFaces src_vector_face_access) const
5511{
5512 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5513 InVector,
5514 OutVector,
5515 CLASS,
5516 false>
5517 worker(*this,
5518 src,
5519 dst,
5520 false,
5521 *owning_class,
5522 cell_operation,
5523 inner_face_operation,
5524 boundary_face_operation,
5525 src_vector_face_access,
5526 dst_vector_face_access,
5527 operation_before_loop,
5528 operation_after_loop,
5529 dof_handler_index_pre_post);
5530 task_info.loop(worker);
5531}
5532
5533
5534
5535template <int dim, typename Number, typename VectorizedArrayType>
5536template <typename CLASS, typename OutVector, typename InVector>
5537inline void
5538MatrixFree<dim, Number, VectorizedArrayType>::loop_cell_centric(
5539 void (CLASS::*function_pointer)(
5540 const MatrixFree<dim, Number, VectorizedArrayType> &,
5541 OutVector &,
5542 const InVector &,
5543 const std::pair<unsigned int, unsigned int> &) const,
5544 const CLASS *owning_class,
5545 OutVector &dst,
5546 const InVector &src,
5547 const bool zero_dst_vector,
5548 const DataAccessOnFaces src_vector_face_access) const
5549{
5550 auto src_vector_face_access_temp = src_vector_face_access;
5551 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5552 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5553 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5554 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5555
5556 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5557 InVector,
5558 OutVector,
5559 CLASS,
5560 true>
5561 worker(*this,
5562 src,
5563 dst,
5564 zero_dst_vector,
5565 *owning_class,
5566 function_pointer,
5567 nullptr,
5568 nullptr,
5569 src_vector_face_access_temp,
5570 DataAccessOnFaces::none);
5571 task_info.loop(worker);
5572}
5573
5574
5575
5576template <int dim, typename Number, typename VectorizedArrayType>
5577template <typename CLASS, typename OutVector, typename InVector>
5578inline void
5579MatrixFree<dim, Number, VectorizedArrayType>::loop_cell_centric(
5580 void (CLASS::*function_pointer)(
5581 const MatrixFree<dim, Number, VectorizedArrayType> &,
5582 OutVector &,
5583 const InVector &,
5584 const std::pair<unsigned int, unsigned int> &),
5585 CLASS *owning_class,
5586 OutVector &dst,
5587 const InVector &src,
5588 const bool zero_dst_vector,
5589 const DataAccessOnFaces src_vector_face_access) const
5590{
5591 auto src_vector_face_access_temp = src_vector_face_access;
5592 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5593 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5594 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5595 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5596
5597 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5598 InVector,
5599 OutVector,
5600 CLASS,
5601 false>
5602 worker(*this,
5603 src,
5604 dst,
5605 zero_dst_vector,
5606 *owning_class,
5607 function_pointer,
5608 nullptr,
5609 nullptr,
5610 src_vector_face_access_temp,
5611 DataAccessOnFaces::none);
5612 task_info.loop(worker);
5613}
5614
5615
5616
5617template <int dim, typename Number, typename VectorizedArrayType>
5618template <typename OutVector, typename InVector>
5619inline void
5620MatrixFree<dim, Number, VectorizedArrayType>::loop_cell_centric(
5621 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5622 OutVector &,
5623 const InVector &,
5624 const std::pair<unsigned int, unsigned int> &)>
5625 &cell_operation,
5626 OutVector &dst,
5627 const InVector &src,
5628 const bool zero_dst_vector,
5629 const DataAccessOnFaces src_vector_face_access) const
5630{
5631 auto src_vector_face_access_temp = src_vector_face_access;
5632 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5633 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5634 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5635 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5636
5637 using Wrapper =
5638 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5639 InVector,
5640 OutVector>;
5641 Wrapper wrap(cell_operation, nullptr, nullptr);
5642
5643 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5644 InVector,
5645 OutVector,
5646 Wrapper,
5647 true>
5648 worker(*this,
5649 src,
5650 dst,
5651 zero_dst_vector,
5652 wrap,
5653 &Wrapper::cell_integrator,
5654 &Wrapper::face_integrator,
5655 &Wrapper::boundary_integrator,
5656 src_vector_face_access_temp,
5657 DataAccessOnFaces::none);
5658 task_info.loop(worker);
5659}
5660
5661
5662#endif // ifndef DOXYGEN
5663
5664
5665
5666DEAL_II_NAMESPACE_CLOSE
5667
5668#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:2290
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:2324
MatrixFree::get_dof_handler
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
MatrixFree::release_scratch_data
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
MatrixFree::loop
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*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::mapping_is_initialized
bool mapping_is_initialized
Definition matrix_free.h:2341
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:2364
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:2303
MatrixFree::get_vector_partitioner
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner(const unsigned int dof_handler_index=0) const
MatrixFree::get_mapping_info
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
MatrixFree::constraint_pool_data
std::vector< Number > constraint_pool_data
Definition matrix_free.h:2277
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:2352
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:2359
MatrixFree::constraint_pool_end
const Number * constraint_pool_end(const unsigned int pool_index) const
MatrixFree::loop_cell_centric
void loop_cell_centric(const std::function< void(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
MatrixFree::dof_handlers
std::vector< SmartPointer< const DoFHandler< dim > > > dof_handlers
Definition matrix_free.h:2255
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::affine_constraints
std::vector< SmartPointer< const AffineConstraints< Number > > > affine_constraints
Definition matrix_free.h:2263
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:2309
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:2269
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:2336
MatrixFree::face_info
internal::MatrixFreeFunctions::FaceInfo< VectorizedArrayType::size()> face_info
Definition matrix_free.h:2331
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:2283
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:2317
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:2295
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
Utilities::MPI::Partitioner::locally_owned_size
unsigned int locally_owned_size() const
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
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:1638
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:1985
AssertThrowMPI
#define AssertThrowMPI(error_code)
Definition exceptions.h:2112
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exceptions.h:2053
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:1739
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
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:697
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::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:748
internal::MatrixFreeFunctions::DoFInfo::cell_loop_post_list_index
std::vector< unsigned int > cell_loop_post_list_index
Definition dof_info.h:761
internal::MatrixFreeFunctions::DoFInfo::vector_zero_range_list
std::vector< std::pair< unsigned int, unsigned int > > vector_zero_range_list
Definition dof_info.h:741
internal::MatrixFreeFunctions::DoFInfo::vector_partitioner
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition dof_info.h:592
internal::MatrixFreeFunctions::DoFInfo::cell_loop_pre_list
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_pre_list
Definition dof_info.h:754
internal::MatrixFreeFunctions::DoFInfo::vector_zero_range_list_index
std::vector< unsigned int > vector_zero_range_list_index
Definition dof_info.h:736
internal::MatrixFreeFunctions::DoFInfo::cell_loop_post_list
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_post_list
Definition dof_info.h:767
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.11.0