Reference documentation for deal.II version GIT relicensing-756-gc37f912904 2024-05-25 08:20:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
matrix_free.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2012 - 2023 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
27
29
30#include <deal.II/fe/fe.h>
31#include <deal.II/fe/mapping.h>
32
35
40
47
48#include <cstdlib>
49#include <limits>
50#include <list>
51#include <memory>
52
53
55
56
57
109template <int dim,
110 typename Number = double,
111 typename VectorizedArrayType = VectorizedArray<Number>>
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
184 {
189
217
223 const unsigned int tasks_block_size = 0,
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)
252 , communicator_sm(MPI_COMM_SELF)
253 {}
254
280
310
348
358 unsigned int tasks_block_size;
359
374
394
414
442
451 unsigned int mg_level;
452
460
471
480
494
503
525 std::vector<unsigned int> cell_vectorization_category;
526
537
549
554 };
555
564
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
646 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free_base);
647
657 void
659
663 void
664 update_mapping(const std::shared_ptr<hp::MappingCollection<dim>> &mapping);
665
670 void
672
687 {
694 none,
695
706 values,
707
717
728 gradients,
729
739
747 };
748
793 template <typename OutVector, typename InVector>
794 void
795 cell_loop(const std::function<void(
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(
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
1093 const std::function<
1095 OutVector &,
1096 const InVector &,
1097 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1098 const std::function<void(
1100 OutVector &,
1101 const InVector &,
1102 const std::pair<unsigned int, unsigned int> &)> &inner_face_operation,
1103 const std::function<void(
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 =
1113 const DataAccessOnFaces src_vector_face_access =
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 =
1227 const DataAccessOnFaces src_vector_face_access =
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 =
1256 const DataAccessOnFaces src_vector_face_access =
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 =
1407 const DataAccessOnFaces src_vector_face_access =
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 =
1440 const DataAccessOnFaces src_vector_face_access =
1442
1448 template <typename OutVector, typename InVector>
1449 void
1451 const std::function<
1453 OutVector &,
1454 const InVector &,
1455 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1456 const std::function<void(
1458 OutVector &,
1459 const InVector &,
1460 const std::pair<unsigned int, unsigned int> &)> &inner_face_operation,
1461 const std::function<void(
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 =
1475 const DataAccessOnFaces src_vector_face_access =
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 =
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 =
1572
1576 template <typename OutVector, typename InVector>
1577 void
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 =
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>
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
1619
1623 unsigned int
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
1648
1654 template <typename T>
1655 void
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
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
1761
1765 unsigned int
1767
1772 unsigned int
1773 n_base_elements(const unsigned int dof_handler_index) const;
1774
1782 unsigned int
1784
1794 unsigned int
1796
1803 unsigned int
1805
1815 unsigned int
1817
1828 unsigned int
1830
1835 unsigned int
1837
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
1873 get_affine_constraints(const unsigned int dof_handler_index = 0) const;
1874
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
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
2026 const std::pair<unsigned int, unsigned int> cell_batch_range) const;
2027
2032 std::pair<unsigned int, unsigned int>
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
2062
2067 bool
2069
2074 unsigned int
2076
2081 std::size_t
2083
2088 template <typename StreamType>
2089 void
2090 print_memory_consumption(StreamType &out) const;
2091
2096 void
2097 print(std::ostream &out) const;
2098
2112
2113 /*
2114 * Return geometry-dependent information on the cells.
2115 */
2116 const internal::MatrixFreeFunctions::
2117 MappingInfo<dim, Number, VectorizedArrayType> &
2119
2124 get_dof_info(const unsigned int dof_handler_index_component = 0) const;
2125
2129 unsigned int
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
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
2161 VectorizedArrayType::size()> &
2162 get_face_info(const unsigned int face_batch_index) const;
2163
2164
2172
2188
2192 void
2194
2206
2210 void
2212 const AlignedVector<Number> *memory) const;
2213
2216private:
2221 template <typename number2, int q_dim>
2222 void
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
2240 const std::vector<const AffineConstraints<number2> *> &constraint,
2241 const std::vector<IndexSet> &locally_owned_set,
2242 const AdditionalData &additional_data);
2243
2247 void
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
2291
2296
2303 std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
2304
2309 std::vector<unsigned int> mf_cell_indices;
2310
2318
2325
2330 internal::MatrixFreeFunctions::FaceInfo<VectorizedArrayType::size()>
2332
2337
2342
2351 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2353
2358 mutable std::list<std::pair<bool, AlignedVector<Number>>>
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
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
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
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
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> &
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> &
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
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
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>
2477{
2478 return task_info;
2479}
2480
2481
2482
2483template <int dim, typename Number, typename VectorizedArrayType>
2484inline unsigned int
2486{
2487 return task_info.n_active_cells;
2488}
2489
2490
2491
2492template <int dim, typename Number, typename VectorizedArrayType>
2493inline unsigned int
2495{
2496 return *(task_info.cell_partition_data.end() - 2);
2497}
2498
2499
2500
2501template <int dim, typename Number, typename VectorizedArrayType>
2502inline unsigned int
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
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
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
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
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()>
2563 const unsigned int cell_batch_index,
2564 const unsigned int face_number) const
2565{
2566 AssertIndexRange(cell_batch_index, n_cell_batches());
2568 Assert(face_info.cell_and_face_boundary_id.size(0) >= n_cell_batches(),
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> &
2586{
2587 return mapping_info;
2588}
2589
2590
2591
2592template <int dim, typename Number, typename VectorizedArrayType>
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
2606{
2607 return constraint_pool_row_index.size() - 1;
2608}
2609
2610
2611
2612template <int dim, typename Number, typename VectorizedArrayType>
2613inline const Number *
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 *
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>
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 {
2648 dof_info[dof_handler_component].fe_index_conversion.size(), 1);
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
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
2687{
2688 return shape_info.size(2);
2689}
2690
2691
2692template <int dim, typename Number, typename VectorizedArrayType>
2693unsigned int
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
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
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
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
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
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
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
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
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
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 &
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 &
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>
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>
2880 VectorizedArrayType::size()> &
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> &
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> &
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> &
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
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>
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
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>
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] !=
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] !=
2998 for (unsigned int v = 0;
2999 v < VectorizedArrayType::size() &&
3000 face_info.faces[face_batch_index].cells_exterior[v] !=
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
3018{
3019 return indices_are_initialized;
3020}
3021
3022
3023
3024template <int dim, typename Number, typename VectorizedArrayType>
3025inline bool
3027{
3028 return mapping_is_initialized;
3029}
3030
3031
3032template <int dim, typename Number, typename VectorizedArrayType>
3033inline unsigned int
3035{
3036 return mg_level;
3037}
3038
3039
3040
3041template <int dim, typename Number, typename VectorizedArrayType>
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
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>
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
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)
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
3152 const MappingType &mapping,
3153 const DoFHandler<dim> &dof_handler,
3154 const AffineConstraints<number2> &constraints_in,
3155 const QuadratureType &quad,
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
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,
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
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,
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 !=
3288 DataAccessOnFaces::unspecified)
3289 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
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
3339 }
3340
3341
3342
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 ==
3355 DataAccessOnFaces::none)
3356 return *matrix_free.get_dof_info(mf_component)
3357 .vector_exchanger_face_variants[0];
3358 else if (vector_face_access ==
3360 DataAccessOnFaces::values)
3361 return *matrix_free.get_dof_info(mf_component)
3362 .vector_exchanger_face_variants[1];
3363 else if (vector_face_access ==
3365 DataAccessOnFaces::gradients)
3366 return *matrix_free.get_dof_info(mf_component)
3367 .vector_exchanger_face_variants[2];
3368 else if (vector_face_access ==
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 ==
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,
3414
3415 if (ghosts_set)
3416 ghosts_were_set = true;
3417
3418 vec.update_ghost_values();
3419 }
3420
3421
3422
3428 template <typename VectorType,
3429 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3430 !has_exchange_on_subset<VectorType>,
3431 VectorType> * = nullptr>
3432 void
3433 update_ghost_values_start(const unsigned int component_in_block_vector,
3434 const VectorType &vec)
3435 {
3436 (void)component_in_block_vector;
3437 const bool ghosts_set = vec.has_ghost_elements();
3438
3439 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3440 ghosts_set == false,
3442
3443 if (ghosts_set)
3444 ghosts_were_set = true;
3445
3446 vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3447 }
3448
3449
3450
3457 template <typename VectorType,
3458 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3459 has_exchange_on_subset<VectorType>,
3460 VectorType> * = nullptr>
3461 void
3462 update_ghost_values_start(const unsigned int component_in_block_vector,
3463 const VectorType &vec)
3464 {
3465 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3466 "Type mismatch between VectorType and VectorDataExchange");
3467 (void)component_in_block_vector;
3468 const bool ghosts_set = vec.has_ghost_elements();
3469
3470 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3471 ghosts_set == false,
3473
3474 if (ghosts_set)
3475 ghosts_were_set = true;
3476
3477 if (vec.size() != 0)
3478 {
3479# ifdef DEAL_II_WITH_MPI
3480 const unsigned int mf_component = find_vector_in_mf(vec);
3481
3482 const auto &part = get_partitioner(mf_component);
3483
3484 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3485 part.n_import_sm_procs() == 0)
3486 return;
3487
3488 tmp_data[component_in_block_vector] =
3489 matrix_free.acquire_scratch_data_non_threadsafe();
3490 tmp_data[component_in_block_vector]->resize_fast(
3491 part.n_import_indices());
3492 AssertDimension(requests.size(), tmp_data.size());
3493
3494 part.export_to_ghosted_array_start(
3495 component_in_block_vector * 2 + channel_shift,
3496 ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3497 vec.shared_vector_data(),
3498 ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3499 part.locally_owned_size(),
3500 matrix_free.get_dof_info(mf_component)
3501 .vector_partitioner->n_ghost_indices()),
3502 ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3503 part.n_import_indices()),
3504 this->requests[component_in_block_vector]);
3505# endif
3506 }
3507 }
3508
3509
3510
3515 template <typename VectorType,
3516 std::enable_if_t<!has_update_ghost_values_start<VectorType>,
3517 VectorType> * = nullptr>
3518 void
3519 update_ghost_values_finish(const unsigned int /*component_in_block_vector*/,
3520 const VectorType & /*vec*/)
3521 {}
3522
3523
3524
3530 template <typename VectorType,
3531 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3532 !has_exchange_on_subset<VectorType>,
3533 VectorType> * = nullptr>
3534 void
3535 update_ghost_values_finish(const unsigned int component_in_block_vector,
3536 const VectorType &vec)
3537 {
3538 (void)component_in_block_vector;
3539 vec.update_ghost_values_finish();
3540 }
3541
3542
3543
3550 template <typename VectorType,
3551 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3552 has_exchange_on_subset<VectorType>,
3553 VectorType> * = nullptr>
3554 void
3555 update_ghost_values_finish(const unsigned int component_in_block_vector,
3556 const VectorType &vec)
3557 {
3558 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3559 "Type mismatch between VectorType and VectorDataExchange");
3560 (void)component_in_block_vector;
3561
3562 if (vec.size() != 0)
3563 {
3564# ifdef DEAL_II_WITH_MPI
3565 AssertIndexRange(component_in_block_vector, tmp_data.size());
3566 AssertDimension(requests.size(), tmp_data.size());
3567
3568 const unsigned int mf_component = find_vector_in_mf(vec);
3569
3570 const auto &part = get_partitioner(mf_component);
3571
3572 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3573 part.n_import_sm_procs() != 0)
3574 {
3575 part.export_to_ghosted_array_finish(
3576 ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3577 vec.shared_vector_data(),
3578 ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3579 part.locally_owned_size(),
3580 matrix_free.get_dof_info(mf_component)
3581 .vector_partitioner->n_ghost_indices()),
3582 this->requests[component_in_block_vector]);
3583
3584 matrix_free.release_scratch_data_non_threadsafe(
3585 tmp_data[component_in_block_vector]);
3586 tmp_data[component_in_block_vector] = nullptr;
3587 }
3588# endif
3589 }
3590 // let vector know that ghosts are being updated and we can read from
3591 // them
3592 vec.set_ghost_state(true);
3593 }
3594
3595
3596
3600 template <typename VectorType,
3601 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3602 * = nullptr>
3603 void
3604 compress_start(const unsigned int /*component_in_block_vector*/,
3605 VectorType & /*vec*/)
3606 {}
3607
3608
3609
3614 template <typename VectorType,
3615 std::enable_if_t<!has_compress_start<VectorType> &&
3616 !is_not_parallel_vector<VectorType>,
3617 VectorType> * = nullptr>
3618 void
3619 compress_start(const unsigned int component_in_block_vector,
3620 VectorType &vec)
3621 {
3622 (void)component_in_block_vector;
3623 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3624 vec.compress(VectorOperation::add);
3625 }
3626
3627
3628
3634 template <typename VectorType,
3635 std::enable_if_t<has_compress_start<VectorType> &&
3636 !has_exchange_on_subset<VectorType>,
3637 VectorType> * = nullptr>
3638 void
3639 compress_start(const unsigned int component_in_block_vector,
3640 VectorType &vec)
3641 {
3642 (void)component_in_block_vector;
3643 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3644 vec.compress_start(component_in_block_vector + channel_shift);
3645 }
3646
3647
3648
3655 template <typename VectorType,
3656 std::enable_if_t<has_compress_start<VectorType> &&
3657 has_exchange_on_subset<VectorType>,
3658 VectorType> * = nullptr>
3659 void
3660 compress_start(const unsigned int component_in_block_vector,
3661 VectorType &vec)
3662 {
3663 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3664 "Type mismatch between VectorType and VectorDataExchange");
3665 (void)component_in_block_vector;
3666 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3667
3668 if (vec.size() != 0)
3669 {
3670# ifdef DEAL_II_WITH_MPI
3671 const unsigned int mf_component = find_vector_in_mf(vec);
3672
3673 const auto &part = get_partitioner(mf_component);
3674
3675 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3676 part.n_import_sm_procs() == 0)
3677 return;
3678
3679 tmp_data[component_in_block_vector] =
3680 matrix_free.acquire_scratch_data_non_threadsafe();
3681 tmp_data[component_in_block_vector]->resize_fast(
3682 part.n_import_indices());
3683 AssertDimension(requests.size(), tmp_data.size());
3684
3685 part.import_from_ghosted_array_start(
3687 component_in_block_vector * 2 + channel_shift,
3688 ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3689 vec.shared_vector_data(),
3690 ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3691 matrix_free.get_dof_info(mf_component)
3692 .vector_partitioner->n_ghost_indices()),
3693 ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3694 part.n_import_indices()),
3695 this->requests[component_in_block_vector]);
3696# endif
3697 }
3698 }
3699
3700
3701
3706 template <
3707 typename VectorType,
3708 std::enable_if_t<!has_compress_start<VectorType>, VectorType> * = nullptr>
3709 void
3710 compress_finish(const unsigned int /*component_in_block_vector*/,
3711 VectorType & /*vec*/)
3712 {}
3713
3714
3715
3721 template <typename VectorType,
3722 std::enable_if_t<has_compress_start<VectorType> &&
3723 !has_exchange_on_subset<VectorType>,
3724 VectorType> * = nullptr>
3725 void
3726 compress_finish(const unsigned int component_in_block_vector,
3727 VectorType &vec)
3728 {
3729 (void)component_in_block_vector;
3730 vec.compress_finish(VectorOperation::add);
3731 }
3732
3733
3734
3741 template <typename VectorType,
3742 std::enable_if_t<has_compress_start<VectorType> &&
3743 has_exchange_on_subset<VectorType>,
3744 VectorType> * = nullptr>
3745 void
3746 compress_finish(const unsigned int component_in_block_vector,
3747 VectorType &vec)
3748 {
3749 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3750 "Type mismatch between VectorType and VectorDataExchange");
3751 (void)component_in_block_vector;
3752 if (vec.size() != 0)
3753 {
3754# ifdef DEAL_II_WITH_MPI
3755 AssertIndexRange(component_in_block_vector, tmp_data.size());
3756 AssertDimension(requests.size(), tmp_data.size());
3757
3758 const unsigned int mf_component = find_vector_in_mf(vec);
3759
3760 const auto &part = get_partitioner(mf_component);
3761
3762 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3763 part.n_import_sm_procs() != 0)
3764 {
3765 part.import_from_ghosted_array_finish(
3767 ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3768 vec.shared_vector_data(),
3769 ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3770 matrix_free.get_dof_info(mf_component)
3771 .vector_partitioner->n_ghost_indices()),
3773 tmp_data[component_in_block_vector]->begin(),
3774 part.n_import_indices()),
3775 this->requests[component_in_block_vector]);
3776
3777 matrix_free.release_scratch_data_non_threadsafe(
3778 tmp_data[component_in_block_vector]);
3779 tmp_data[component_in_block_vector] = nullptr;
3780 }
3781
3783 {
3784 const int ierr =
3785 MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3786 AssertThrowMPI(ierr);
3787 }
3788# endif
3789 }
3790 }
3791
3792
3793
3797 template <typename VectorType,
3798 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3799 * = nullptr>
3800 void
3801 reset_ghost_values(const VectorType & /*vec*/) const
3802 {}
3803
3804
3805
3810 template <typename VectorType,
3811 std::enable_if_t<!has_exchange_on_subset<VectorType> &&
3812 !is_not_parallel_vector<VectorType>,
3813 VectorType> * = nullptr>
3814 void
3815 reset_ghost_values(const VectorType &vec) const
3816 {
3817 if (ghosts_were_set == true)
3818 return;
3819
3820 vec.zero_out_ghost_values();
3821 }
3822
3823
3824
3830 template <typename VectorType,
3831 std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3832 * = nullptr>
3833 void
3834 reset_ghost_values(const VectorType &vec) const
3835 {
3836 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3837 "Type mismatch between VectorType and VectorDataExchange");
3838 if (ghosts_were_set == true)
3839 return;
3840
3841 if (vec.size() != 0)
3842 {
3843# ifdef DEAL_II_WITH_MPI
3844 AssertDimension(requests.size(), tmp_data.size());
3845
3846 const unsigned int mf_component = find_vector_in_mf(vec);
3847
3848 const auto &part = get_partitioner(mf_component);
3849
3850 if (part.n_ghost_indices() > 0)
3851 {
3852 part.reset_ghost_values(ArrayView<Number>(
3854 .begin() +
3855 part.locally_owned_size(),
3856 matrix_free.get_dof_info(mf_component)
3857 .vector_partitioner->n_ghost_indices()));
3858 }
3859
3860# endif
3861 }
3862 // let vector know that it's not ghosted anymore
3863 vec.set_ghost_state(false);
3864 }
3865
3866
3867
3873 template <typename VectorType,
3874 std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3875 * = nullptr>
3876 void
3877 zero_vector_region(const unsigned int range_index, VectorType &vec) const
3878 {
3879 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3880 "Type mismatch between VectorType and VectorDataExchange");
3881 if (range_index == numbers::invalid_unsigned_int)
3882 vec = Number();
3883 else
3884 {
3885 const unsigned int mf_component = find_vector_in_mf(vec, false);
3887 matrix_free.get_dof_info(mf_component);
3888 Assert(dof_info.vector_zero_range_list_index.empty() == false,
3890
3891 Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3893 AssertIndexRange(range_index,
3894 dof_info.vector_zero_range_list_index.size() - 1);
3895 for (unsigned int id =
3896 dof_info.vector_zero_range_list_index[range_index];
3897 id != dof_info.vector_zero_range_list_index[range_index + 1];
3898 ++id)
3899 std::memset(vec.begin() + dof_info.vector_zero_range_list[id].first,
3900 0,
3901 (dof_info.vector_zero_range_list[id].second -
3902 dof_info.vector_zero_range_list[id].first) *
3903 sizeof(Number));
3904 }
3905 }
3906
3907
3908
3914 template <typename VectorType,
3915 std::enable_if_t<!has_exchange_on_subset<VectorType>, VectorType>
3916 * = nullptr,
3917 typename VectorType::value_type * = nullptr>
3918 void
3919 zero_vector_region(const unsigned int range_index, VectorType &vec) const
3920 {
3921 if (range_index == numbers::invalid_unsigned_int || range_index == 0)
3922 {
3923 if constexpr (std::is_same_v<
3925 VectorType>)
3926 {
3927 for (unsigned int i = 0; i < vec.size(); ++i)
3928 vec[i] = typename VectorType::value_type();
3929 }
3930 else
3931 vec = typename VectorType::value_type();
3932 }
3933 }
3934
3935
3936
3941 void
3942 zero_vector_region(const unsigned int, ...) const
3943 {
3944 Assert(false,
3945 ExcNotImplemented("Zeroing is only implemented for vector types "
3946 "which provide VectorType::value_type"));
3947 }
3948
3949
3950
3951 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
3952 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3953 DataAccessOnFaces vector_face_access;
3954 bool ghosts_were_set;
3955# ifdef DEAL_II_WITH_MPI
3956 std::vector<AlignedVector<Number> *> tmp_data;
3957 std::vector<std::vector<MPI_Request>> requests;
3958# endif
3959 }; // VectorDataExchange
3960
3961 template <typename VectorStruct>
3962 unsigned int
3963 n_components(const VectorStruct &vec);
3964
3965 template <typename VectorStruct>
3966 unsigned int
3967 n_components_block(const VectorStruct &vec, const std::bool_constant<true>)
3968 {
3969 unsigned int components = 0;
3970 for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
3971 components += n_components(vec.block(bl));
3972 return components;
3973 }
3974
3975 template <typename VectorStruct>
3976 unsigned int
3977 n_components_block(const VectorStruct &, const std::bool_constant<false>)
3978 {
3979 return 1;
3980 }
3981
3982 template <typename VectorStruct>
3983 unsigned int
3984 n_components(const VectorStruct &vec)
3985 {
3986 return n_components_block(
3987 vec, std::bool_constant<IsBlockVector<VectorStruct>::value>());
3988 }
3989
3990 template <typename VectorStruct>
3991 inline unsigned int
3992 n_components(const std::vector<VectorStruct> &vec)
3993 {
3994 unsigned int components = 0;
3995 for (unsigned int comp = 0; comp < vec.size(); ++comp)
3996 components += n_components_block(
3997 vec[comp], std::bool_constant<IsBlockVector<VectorStruct>::value>());
3998 return components;
3999 }
4000
4001 template <typename VectorStruct>
4002 inline unsigned int
4003 n_components(const std::vector<VectorStruct *> &vec)
4004 {
4005 unsigned int components = 0;
4006 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4007 components += n_components_block(
4008 *vec[comp], std::bool_constant<IsBlockVector<VectorStruct>::value>());
4009 return components;
4010 }
4011
4012
4013
4014 // A helper function to identify block vectors with many components where we
4015 // should not try to overlap computations and communication because there
4016 // would be too many outstanding communication requests.
4017
4018 // default value for vectors that do not have communication_block_size
4019 template <typename VectorStruct,
4020 std::enable_if_t<!has_communication_block_size<VectorStruct>,
4021 VectorStruct> * = nullptr>
4022 constexpr unsigned int
4023 get_communication_block_size(const VectorStruct &)
4024 {
4026 }
4027
4028
4029
4030 template <typename VectorStruct,
4031 std::enable_if_t<has_communication_block_size<VectorStruct>,
4032 VectorStruct> * = nullptr>
4033 constexpr unsigned int
4034 get_communication_block_size(const VectorStruct &)
4035 {
4036 return VectorStruct::communication_block_size;
4037 }
4038
4039
4040
4041 template <typename VectorType,
4042 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType> * =
4043 nullptr>
4044 bool
4045 has_ghost_elements(const VectorType &vec)
4046 {
4047 (void)vec;
4048 return false;
4049 }
4050
4051
4052
4053 template <typename VectorType,
4054 std::enable_if_t<!is_not_parallel_vector<VectorType>, VectorType>
4055 * = nullptr>
4056 bool
4057 has_ghost_elements(const VectorType &vec)
4058 {
4059 return vec.has_ghost_elements();
4060 }
4061
4062
4063
4064 // --------------------------------------------------------------------------
4065 // below we have wrappers to distinguish between block and non-block vectors.
4066 // --------------------------------------------------------------------------
4067
4068 //
4069 // update_ghost_values_start
4070 //
4071
4072 // update_ghost_values for block vectors
4073 template <int dim,
4074 typename VectorStruct,
4075 typename Number,
4076 typename VectorizedArrayType,
4077 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4078 * = nullptr>
4079 void
4080 update_ghost_values_start(
4081 const VectorStruct &vec,
4082 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4083 const unsigned int channel = 0)
4084 {
4085 if (get_communication_block_size(vec) < vec.n_blocks())
4086 {
4087 const bool ghosts_set = vec.has_ghost_elements();
4088
4089 Assert(exchanger.matrix_free.get_task_info()
4090 .allow_ghosted_vectors_in_loops ||
4091 ghosts_set == false,
4093
4094 if (ghosts_set)
4095 exchanger.ghosts_were_set = true;
4096
4097 vec.update_ghost_values();
4098 }
4099 else
4100 {
4101 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4102 update_ghost_values_start(vec.block(i), exchanger, channel + i);
4103 }
4104 }
4105
4106
4107
4108 // update_ghost_values for non-block vectors
4109 template <int dim,
4110 typename VectorStruct,
4111 typename Number,
4112 typename VectorizedArrayType,
4113 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4114 * = nullptr>
4115 void
4116 update_ghost_values_start(
4117 const VectorStruct &vec,
4118 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4119 const unsigned int channel = 0)
4120 {
4121 exchanger.update_ghost_values_start(channel, vec);
4122 }
4123
4124
4125
4126 // update_ghost_values_start() for vector of vectors
4127 template <int dim,
4128 typename VectorStruct,
4129 typename Number,
4130 typename VectorizedArrayType>
4131 inline void
4132 update_ghost_values_start(
4133 const std::vector<VectorStruct> &vec,
4134 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4135 {
4136 unsigned int component_index = 0;
4137 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4138 {
4139 update_ghost_values_start(vec[comp], exchanger, component_index);
4140 component_index += n_components(vec[comp]);
4141 }
4142 }
4143
4144
4145
4146 // update_ghost_values_start() for vector of pointers to vectors
4147 template <int dim,
4148 typename VectorStruct,
4149 typename Number,
4150 typename VectorizedArrayType>
4151 inline void
4152 update_ghost_values_start(
4153 const std::vector<VectorStruct *> &vec,
4154 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4155 {
4156 unsigned int component_index = 0;
4157 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4158 {
4159 update_ghost_values_start(*vec[comp], exchanger, component_index);
4160 component_index += n_components(*vec[comp]);
4161 }
4162 }
4163
4164
4165
4166 //
4167 // update_ghost_values_finish
4168 //
4169
4170 // for block vectors
4171 template <int dim,
4172 typename VectorStruct,
4173 typename Number,
4174 typename VectorizedArrayType,
4175 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4176 * = nullptr>
4177 void
4178 update_ghost_values_finish(
4179 const VectorStruct &vec,
4180 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4181 const unsigned int channel = 0)
4182 {
4183 if (get_communication_block_size(vec) < vec.n_blocks())
4184 {
4185 // do nothing, everything has already been completed in the _start()
4186 // call
4187 }
4188 else
4189 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4190 update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4191 }
4192
4193
4194
4195 // for non-block vectors
4196 template <int dim,
4197 typename VectorStruct,
4198 typename Number,
4199 typename VectorizedArrayType,
4200 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4201 * = nullptr>
4202 void
4203 update_ghost_values_finish(
4204 const VectorStruct &vec,
4205 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4206 const unsigned int channel = 0)
4207 {
4208 exchanger.update_ghost_values_finish(channel, vec);
4209 }
4210
4211
4212
4213 // for vector of vectors
4214 template <int dim,
4215 typename VectorStruct,
4216 typename Number,
4217 typename VectorizedArrayType>
4218 inline void
4219 update_ghost_values_finish(
4220 const std::vector<VectorStruct> &vec,
4221 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4222 {
4223 unsigned int component_index = 0;
4224 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4225 {
4226 update_ghost_values_finish(vec[comp], exchanger, component_index);
4227 component_index += n_components(vec[comp]);
4228 }
4229 }
4230
4231
4232
4233 // for vector of pointers to vectors
4234 template <int dim,
4235 typename VectorStruct,
4236 typename Number,
4237 typename VectorizedArrayType>
4238 inline void
4239 update_ghost_values_finish(
4240 const std::vector<VectorStruct *> &vec,
4241 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4242 {
4243 unsigned int component_index = 0;
4244 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4245 {
4246 update_ghost_values_finish(*vec[comp], exchanger, component_index);
4247 component_index += n_components(*vec[comp]);
4248 }
4249 }
4250
4251
4252
4253 //
4254 // compress_start
4255 //
4256
4257 // for block vectors
4258 template <int dim,
4259 typename VectorStruct,
4260 typename Number,
4261 typename VectorizedArrayType,
4262 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4263 * = nullptr>
4264 inline void
4265 compress_start(
4266 VectorStruct &vec,
4267 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4268 const unsigned int channel = 0)
4269 {
4270 if (get_communication_block_size(vec) < vec.n_blocks())
4271 vec.compress(VectorOperation::add);
4272 else
4273 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4274 compress_start(vec.block(i), exchanger, channel + i);
4275 }
4276
4277
4278
4279 // for non-block vectors
4280 template <int dim,
4281 typename VectorStruct,
4282 typename Number,
4283 typename VectorizedArrayType,
4284 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4285 * = nullptr>
4286 inline void
4287 compress_start(
4288 VectorStruct &vec,
4289 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4290 const unsigned int channel = 0)
4291 {
4292 exchanger.compress_start(channel, vec);
4293 }
4294
4295
4296
4297 // for std::vector of vectors
4298 template <int dim,
4299 typename VectorStruct,
4300 typename Number,
4301 typename VectorizedArrayType>
4302 inline void
4303 compress_start(
4304 std::vector<VectorStruct> &vec,
4305 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4306 {
4307 unsigned int component_index = 0;
4308 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4309 {
4310 compress_start(vec[comp], exchanger, component_index);
4311 component_index += n_components(vec[comp]);
4312 }
4313 }
4314
4315
4316
4317 // for std::vector of pointer to vectors
4318 template <int dim,
4319 typename VectorStruct,
4320 typename Number,
4321 typename VectorizedArrayType>
4322 inline void
4323 compress_start(
4324 std::vector<VectorStruct *> &vec,
4325 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4326 {
4327 unsigned int component_index = 0;
4328 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4329 {
4330 compress_start(*vec[comp], exchanger, component_index);
4331 component_index += n_components(*vec[comp]);
4332 }
4333 }
4334
4335
4336
4337 //
4338 // compress_finish
4339 //
4340
4341 // for block vectors
4342 template <int dim,
4343 typename VectorStruct,
4344 typename Number,
4345 typename VectorizedArrayType,
4346 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4347 * = nullptr>
4348 inline void
4349 compress_finish(
4350 VectorStruct &vec,
4351 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4352 const unsigned int channel = 0)
4353 {
4354 if (get_communication_block_size(vec) < vec.n_blocks())
4355 {
4356 // do nothing, everything has already been completed in the _start()
4357 // call
4358 }
4359 else
4360 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4361 compress_finish(vec.block(i), exchanger, channel + i);
4362 }
4363
4364
4365
4366 // for non-block vectors
4367 template <int dim,
4368 typename VectorStruct,
4369 typename Number,
4370 typename VectorizedArrayType,
4371 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4372 * = nullptr>
4373 inline void
4374 compress_finish(
4375 VectorStruct &vec,
4376 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4377 const unsigned int channel = 0)
4378 {
4379 exchanger.compress_finish(channel, vec);
4380 }
4381
4382
4383
4384 // for std::vector of vectors
4385 template <int dim,
4386 typename VectorStruct,
4387 typename Number,
4388 typename VectorizedArrayType>
4389 inline void
4390 compress_finish(
4391 std::vector<VectorStruct> &vec,
4392 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4393 {
4394 unsigned int component_index = 0;
4395 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4396 {
4397 compress_finish(vec[comp], exchanger, component_index);
4398 component_index += n_components(vec[comp]);
4399 }
4400 }
4401
4402
4403
4404 // for std::vector of pointer to vectors
4405 template <int dim,
4406 typename VectorStruct,
4407 typename Number,
4408 typename VectorizedArrayType>
4409 inline void
4410 compress_finish(
4411 std::vector<VectorStruct *> &vec,
4412 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4413 {
4414 unsigned int component_index = 0;
4415 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4416 {
4417 compress_finish(*vec[comp], exchanger, component_index);
4418 component_index += n_components(*vec[comp]);
4419 }
4420 }
4421
4422
4423
4424 //
4425 // reset_ghost_values:
4426 //
4427 // if the input vector did not have ghosts imported, clear them here again
4428 // in order to avoid subsequent operations e.g. in linear solvers to work
4429 // with ghosts all the time
4430 //
4431
4432 // for block vectors
4433 template <int dim,
4434 typename VectorStruct,
4435 typename Number,
4436 typename VectorizedArrayType,
4437 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4438 * = nullptr>
4439 inline void
4440 reset_ghost_values(
4441 const VectorStruct &vec,
4442 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4443 {
4444 // return immediately if there is nothing to do.
4445 if (exchanger.ghosts_were_set == true)
4446 return;
4447
4448 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4449 reset_ghost_values(vec.block(i), exchanger);
4450 }
4451
4452
4453
4454 // for non-block vectors
4455 template <int dim,
4456 typename VectorStruct,
4457 typename Number,
4458 typename VectorizedArrayType,
4459 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4460 * = nullptr>
4461 inline void
4462 reset_ghost_values(
4463 const VectorStruct &vec,
4464 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4465 {
4466 exchanger.reset_ghost_values(vec);
4467 }
4468
4469
4470
4471 // for std::vector of vectors
4472 template <int dim,
4473 typename VectorStruct,
4474 typename Number,
4475 typename VectorizedArrayType>
4476 inline void
4477 reset_ghost_values(
4478 const std::vector<VectorStruct> &vec,
4479 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4480 {
4481 // return immediately if there is nothing to do.
4482 if (exchanger.ghosts_were_set == true)
4483 return;
4484
4485 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4486 reset_ghost_values(vec[comp], exchanger);
4487 }
4488
4489
4490
4491 // for std::vector of pointer to vectors
4492 template <int dim,
4493 typename VectorStruct,
4494 typename Number,
4495 typename VectorizedArrayType>
4496 inline void
4497 reset_ghost_values(
4498 const std::vector<VectorStruct *> &vec,
4499 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4500 {
4501 // return immediately if there is nothing to do.
4502 if (exchanger.ghosts_were_set == true)
4503 return;
4504
4505 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4506 reset_ghost_values(*vec[comp], exchanger);
4507 }
4508
4509
4510
4511 //
4512 // zero_vector_region
4513 //
4514
4515 // for block vectors
4516 template <int dim,
4517 typename VectorStruct,
4518 typename Number,
4519 typename VectorizedArrayType,
4520 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4521 * = nullptr>
4522 inline void
4523 zero_vector_region(
4524 const unsigned int range_index,
4525 VectorStruct &vec,
4526 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4527 {
4528 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4529 exchanger.zero_vector_region(range_index, vec.block(i));
4530 }
4531
4532
4533
4534 // for non-block vectors
4535 template <int dim,
4536 typename VectorStruct,
4537 typename Number,
4538 typename VectorizedArrayType,
4539 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4540 * = nullptr>
4541 inline void
4542 zero_vector_region(
4543 const unsigned int range_index,
4544 VectorStruct &vec,
4545 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4546 {
4547 exchanger.zero_vector_region(range_index, vec);
4548 }
4549
4550
4551
4552 // for std::vector of vectors
4553 template <int dim,
4554 typename VectorStruct,
4555 typename Number,
4556 typename VectorizedArrayType>
4557 inline void
4558 zero_vector_region(
4559 const unsigned int range_index,
4560 std::vector<VectorStruct> &vec,
4561 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4562 {
4563 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4564 zero_vector_region(range_index, vec[comp], exchanger);
4565 }
4566
4567
4568
4569 // for std::vector of pointers to vectors
4570 template <int dim,
4571 typename VectorStruct,
4572 typename Number,
4573 typename VectorizedArrayType>
4574 inline void
4575 zero_vector_region(
4576 const unsigned int range_index,
4577 std::vector<VectorStruct *> &vec,
4578 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4579 {
4580 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4581 zero_vector_region(range_index, *vec[comp], exchanger);
4582 }
4583
4584
4585
4586 // Apply a unit matrix operation to constrained DoFs: Default cases where we
4587 // cannot detect a LinearAlgebra::distributed::Vector, we do not do
4588 // anything, else we apply the constraints as a unit operation
4589 template <typename VectorStruct1, typename VectorStruct2>
4590 inline void
4591 apply_operation_to_constrained_dofs(const std::vector<unsigned int> &,
4592 const VectorStruct1 &,
4593 VectorStruct2 &)
4594 {}
4595
4596 template <typename Number>
4597 inline void
4598 apply_operation_to_constrained_dofs(
4599 const std::vector<unsigned int> &constrained_dofs,
4602 {
4603 for (const unsigned int i : constrained_dofs)
4604 dst.local_element(i) = src.local_element(i);
4605 }
4606
4607
4608 namespace MatrixFreeFunctions
4609 {
4610 // struct to select between a const interface and a non-const interface
4611 // for MFWorker
4612 template <typename, typename, typename, typename, bool>
4613 struct InterfaceSelector
4614 {};
4615
4616 // Version for constant functions
4617 template <typename MF,
4618 typename InVector,
4619 typename OutVector,
4620 typename Container>
4621 struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4622 {
4623 using function_type = void (Container::*)(
4624 const MF &,
4625 OutVector &,
4626 const InVector &,
4627 const std::pair<unsigned int, unsigned int> &) const;
4628 };
4629
4630 // Version for non-constant functions
4631 template <typename MF,
4632 typename InVector,
4633 typename OutVector,
4634 typename Container>
4635 struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4636 {
4637 using function_type =
4638 void (Container::*)(const MF &,
4639 OutVector &,
4640 const InVector &,
4641 const std::pair<unsigned int, unsigned int> &);
4642 };
4643 } // namespace MatrixFreeFunctions
4644
4645
4646
4647 // A implementation class for the worker object that runs the various
4648 // operations we want to perform during the matrix-free loop
4649 template <typename MF,
4650 typename InVector,
4651 typename OutVector,
4652 typename Container,
4653 bool is_constant>
4654 class MFWorker : public MFWorkerInterface
4655 {
4656 public:
4657 // An alias to make the arguments further down more readable
4658 using function_type = typename MatrixFreeFunctions::
4659 InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4660 function_type;
4661
4662 // constructor, binds all the arguments to this class
4663 MFWorker(const MF &matrix_free,
4664 const InVector &src,
4665 OutVector &dst,
4666 const bool zero_dst_vector_setting,
4667 const Container &container,
4668 function_type cell_function,
4669 function_type face_function,
4670 function_type boundary_function,
4671 const typename MF::DataAccessOnFaces src_vector_face_access =
4672 MF::DataAccessOnFaces::none,
4673 const typename MF::DataAccessOnFaces dst_vector_face_access =
4674 MF::DataAccessOnFaces::none,
4675 const std::function<void(const unsigned int, const unsigned int)>
4676 &operation_before_loop = {},
4677 const std::function<void(const unsigned int, const unsigned int)>
4678 &operation_after_loop = {},
4679 const unsigned int dof_handler_index_pre_post = 0)
4680 : matrix_free(matrix_free)
4681 , container(const_cast<Container &>(container))
4682 , cell_function(cell_function)
4683 , face_function(face_function)
4684 , boundary_function(boundary_function)
4685 , src(src)
4686 , dst(dst)
4687 , src_data_exchanger(matrix_free,
4688 src_vector_face_access,
4689 n_components(src))
4690 , dst_data_exchanger(matrix_free,
4691 dst_vector_face_access,
4692 n_components(dst))
4693 , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4694 , zero_dst_vector_setting(zero_dst_vector_setting &&
4695 !src_and_dst_are_same)
4696 , operation_before_loop(operation_before_loop)
4697 , operation_after_loop(operation_after_loop)
4698 , dof_handler_index_pre_post(dof_handler_index_pre_post)
4699 {
4700 Assert(!has_ghost_elements(dst),
4701 ExcMessage("The destination vector passed to the matrix-free "
4702 "loop is ghosted. This is not allowed."));
4703 }
4704
4705 // Runs the cell work. If no function is given, nothing is done
4706 virtual void
4707 cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4708 {
4709 if (cell_function != nullptr && cell_range.second > cell_range.first)
4710 for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4711 {
4712 const auto cell_subrange =
4713 matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4714
4715 if (cell_subrange.second <= cell_subrange.first)
4716 continue;
4717
4718 (container.*
4719 cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4720 }
4721 }
4722
4723 virtual void
4724 cell(const unsigned int range_index) override
4725 {
4726 process_range(cell_function,
4727 matrix_free.get_task_info().cell_partition_data_hp_ptr,
4728 matrix_free.get_task_info().cell_partition_data_hp,
4729 range_index);
4730 }
4731
4732 virtual void
4733 face(const unsigned int range_index) override
4734 {
4735 process_range(face_function,
4736 matrix_free.get_task_info().face_partition_data_hp_ptr,
4737 matrix_free.get_task_info().face_partition_data_hp,
4738 range_index);
4739 }
4740
4741 virtual void
4742 boundary(const unsigned int range_index) override
4743 {
4744 process_range(boundary_function,
4745 matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4746 matrix_free.get_task_info().boundary_partition_data_hp,
4747 range_index);
4748 }
4749
4750 private:
4751 void
4752 process_range(const function_type &fu,
4753 const std::vector<unsigned int> &ptr,
4754 const std::vector<unsigned int> &data,
4755 const unsigned int range_index)
4756 {
4757 if (fu == nullptr)
4758 return;
4759
4760 AssertIndexRange(range_index + 1, ptr.size());
4761 for (unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4762 {
4763 AssertIndexRange(2 * i + 1, data.size());
4764 (container.*fu)(matrix_free,
4765 this->dst,
4766 this->src,
4767 std::make_pair(data[2 * i], data[2 * i + 1]));
4768 }
4769 }
4770
4771 public:
4772 // Starts the communication for the update ghost values operation. We
4773 // cannot call this update if ghost and destination are the same because
4774 // that would introduce spurious entries in the destination (there is also
4775 // the problem that reading from a vector that we also write to is usually
4776 // not intended in case there is overlap, but this is up to the
4777 // application code to decide and we cannot catch this case here).
4778 virtual void
4779 vector_update_ghosts_start() override
4780 {
4781 if (!src_and_dst_are_same)
4782 internal::update_ghost_values_start(src, src_data_exchanger);
4783 }
4784
4785 // Finishes the communication for the update ghost values operation
4786 virtual void
4787 vector_update_ghosts_finish() override
4788 {
4789 if (!src_and_dst_are_same)
4790 internal::update_ghost_values_finish(src, src_data_exchanger);
4791 }
4792
4793 // Starts the communication for the vector compress operation
4794 virtual void
4795 vector_compress_start() override
4796 {
4797 internal::compress_start(dst, dst_data_exchanger);
4798 }
4799
4800 // Finishes the communication for the vector compress operation
4801 virtual void
4802 vector_compress_finish() override
4803 {
4804 internal::compress_finish(dst, dst_data_exchanger);
4805 if (!src_and_dst_are_same)
4806 internal::reset_ghost_values(src, src_data_exchanger);
4807 }
4808
4809 // Zeros the given input vector
4810 virtual void
4811 zero_dst_vector_range(const unsigned int range_index) override
4812 {
4813 if (zero_dst_vector_setting)
4814 internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4815 }
4816
4817 virtual void
4818 cell_loop_pre_range(const unsigned int range_index) override
4819 {
4820 if (operation_before_loop)
4821 {
4823 matrix_free.get_dof_info(dof_handler_index_pre_post);
4824 if (range_index == numbers::invalid_unsigned_int)
4825 {
4826 // Case with threaded loop -> currently no overlap implemented
4828 0U,
4829 dof_info.vector_partitioner->locally_owned_size(),
4830 operation_before_loop,
4832 }
4833 else
4834 {
4835 AssertIndexRange(range_index,
4836 dof_info.cell_loop_pre_list_index.size() - 1);
4837 for (unsigned int id =
4838 dof_info.cell_loop_pre_list_index[range_index];
4839 id != dof_info.cell_loop_pre_list_index[range_index + 1];
4840 ++id)
4841 operation_before_loop(dof_info.cell_loop_pre_list[id].first,
4842 dof_info.cell_loop_pre_list[id].second);
4843 }
4844 }
4845 }
4846
4847 virtual void
4848 cell_loop_post_range(const unsigned int range_index) override
4849 {
4850 if (operation_after_loop)
4851 {
4852 // Run unit matrix operation on constrained dofs if we are at the
4853 // last range
4854 const std::vector<unsigned int> &partition_row_index =
4855 matrix_free.get_task_info().partition_row_index;
4856 if (range_index ==
4857 partition_row_index[partition_row_index.size() - 2] - 1)
4858 apply_operation_to_constrained_dofs(
4859 matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
4860 src,
4861 dst);
4862
4864 matrix_free.get_dof_info(dof_handler_index_pre_post);
4865 if (range_index == numbers::invalid_unsigned_int)
4866 {
4867 // Case with threaded loop -> currently no overlap implemented
4869 0U,
4870 dof_info.vector_partitioner->locally_owned_size(),
4871 operation_after_loop,
4873 }
4874 else
4875 {
4876 AssertIndexRange(range_index,
4877 dof_info.cell_loop_post_list_index.size() - 1);
4878 for (unsigned int id =
4879 dof_info.cell_loop_post_list_index[range_index];
4880 id != dof_info.cell_loop_post_list_index[range_index + 1];
4881 ++id)
4882 operation_after_loop(dof_info.cell_loop_post_list[id].first,
4883 dof_info.cell_loop_post_list[id].second);
4884 }
4885 }
4886 }
4887
4888 private:
4889 const MF &matrix_free;
4890 Container &container;
4891 function_type cell_function;
4892 function_type face_function;
4893 function_type boundary_function;
4894
4895 const InVector &src;
4896 OutVector &dst;
4897 VectorDataExchange<MF::dimension,
4898 typename MF::value_type,
4899 typename MF::vectorized_value_type>
4900 src_data_exchanger;
4901 VectorDataExchange<MF::dimension,
4902 typename MF::value_type,
4903 typename MF::vectorized_value_type>
4904 dst_data_exchanger;
4905 const bool src_and_dst_are_same;
4906 const bool zero_dst_vector_setting;
4907 const std::function<void(const unsigned int, const unsigned int)>
4908 operation_before_loop;
4909 const std::function<void(const unsigned int, const unsigned int)>
4910 operation_after_loop;
4911 const unsigned int dof_handler_index_pre_post;
4912 };
4913
4914
4915
4920 template <class MF, typename InVector, typename OutVector>
4921 struct MFClassWrapper
4922 {
4923 using function_type =
4924 std::function<void(const MF &,
4925 OutVector &,
4926 const InVector &,
4927 const std::pair<unsigned int, unsigned int> &)>;
4928
4929 MFClassWrapper(const function_type cell,
4930 const function_type face,
4931 const function_type boundary)
4932 : cell(cell)
4933 , face(face)
4934 , boundary(boundary)
4935 {}
4936
4937 void
4938 cell_integrator(const MF &mf,
4939 OutVector &dst,
4940 const InVector &src,
4941 const std::pair<unsigned int, unsigned int> &range) const
4942 {
4943 if (cell)
4944 cell(mf, dst, src, range);
4945 }
4946
4947 void
4948 face_integrator(const MF &mf,
4949 OutVector &dst,
4950 const InVector &src,
4951 const std::pair<unsigned int, unsigned int> &range) const
4952 {
4953 if (face)
4954 face(mf, dst, src, range);
4955 }
4956
4957 void
4958 boundary_integrator(
4959 const MF &mf,
4960 OutVector &dst,
4961 const InVector &src,
4962 const std::pair<unsigned int, unsigned int> &range) const
4963 {
4964 if (boundary)
4965 boundary(mf, dst, src, range);
4966 }
4967
4968 const function_type cell;
4969 const function_type face;
4970 const function_type boundary;
4971 };
4972
4973} // end of namespace internal
4974
4975
4976
4977template <int dim, typename Number, typename VectorizedArrayType>
4978template <typename OutVector, typename InVector>
4979inline void
4981 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4982 OutVector &,
4983 const InVector &,
4984 const std::pair<unsigned int, unsigned int> &)>
4985 &cell_operation,
4986 OutVector &dst,
4987 const InVector &src,
4988 const bool zero_dst_vector) const
4989{
4990 using Wrapper =
4991 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4992 InVector,
4993 OutVector>;
4994 Wrapper wrap(cell_operation, nullptr, nullptr);
4995 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4996 InVector,
4997 OutVector,
4998 Wrapper,
4999 true>
5000 worker(*this,
5001 src,
5002 dst,
5003 zero_dst_vector,
5004 wrap,
5005 &Wrapper::cell_integrator,
5006 &Wrapper::face_integrator,
5007 &Wrapper::boundary_integrator);
5008
5009 task_info.loop(worker);
5010}
5011
5012
5013
5014template <int dim, typename Number, typename VectorizedArrayType>
5015template <typename OutVector, typename InVector>
5016inline void
5018 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5019 OutVector &,
5020 const InVector &,
5021 const std::pair<unsigned int, unsigned int> &)>
5022 &cell_operation,
5023 OutVector &dst,
5024 const InVector &src,
5025 const std::function<void(const unsigned int, const unsigned int)>
5026 &operation_before_loop,
5027 const std::function<void(const unsigned int, const unsigned int)>
5028 &operation_after_loop,
5029 const unsigned int dof_handler_index_pre_post) const
5030{
5031 using Wrapper =
5032 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5033 InVector,
5034 OutVector>;
5035 Wrapper wrap(cell_operation, nullptr, nullptr);
5036 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5037 InVector,
5038 OutVector,
5039 Wrapper,
5040 true>
5041 worker(*this,
5042 src,
5043 dst,
5044 false,
5045 wrap,
5046 &Wrapper::cell_integrator,
5047 &Wrapper::face_integrator,
5048 &Wrapper::boundary_integrator,
5049 DataAccessOnFaces::none,
5050 DataAccessOnFaces::none,
5051 operation_before_loop,
5052 operation_after_loop,
5053 dof_handler_index_pre_post);
5054
5055 task_info.loop(worker);
5056}
5057
5058
5059
5060template <int dim, typename Number, typename VectorizedArrayType>
5061template <typename OutVector, typename InVector>
5062inline void
5064 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5065 OutVector &,
5066 const InVector &,
5067 const std::pair<unsigned int, unsigned int> &)>
5068 &cell_operation,
5069 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5070 OutVector &,
5071 const InVector &,
5072 const std::pair<unsigned int, unsigned int> &)>
5073 &inner_face_operation,
5074 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5075 OutVector &,
5076 const InVector &,
5077 const std::pair<unsigned int, unsigned int> &)>
5078 &boundary_face_operation,
5079 OutVector &dst,
5080 const InVector &src,
5081 const bool zero_dst_vector,
5082 const DataAccessOnFaces dst_vector_face_access,
5083 const DataAccessOnFaces src_vector_face_access) const
5084{
5085 using Wrapper =
5086 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5087 InVector,
5088 OutVector>;
5089 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5090 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5091 InVector,
5092 OutVector,
5093 Wrapper,
5094 true>
5095 worker(*this,
5096 src,
5097 dst,
5098 zero_dst_vector,
5099 wrap,
5100 &Wrapper::cell_integrator,
5101 &Wrapper::face_integrator,
5102 &Wrapper::boundary_integrator,
5103 src_vector_face_access,
5104 dst_vector_face_access);
5105
5106 task_info.loop(worker);
5107}
5108
5109
5110
5111template <int dim, typename Number, typename VectorizedArrayType>
5112template <typename CLASS, typename OutVector, typename InVector>
5113inline void
5115 void (CLASS::*function_pointer)(
5116 const MatrixFree<dim, Number, VectorizedArrayType> &,
5117 OutVector &,
5118 const InVector &,
5119 const std::pair<unsigned int, unsigned int> &) const,
5120 const CLASS *owning_class,
5121 OutVector &dst,
5122 const InVector &src,
5123 const bool zero_dst_vector) const
5124{
5125 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5126 InVector,
5127 OutVector,
5128 CLASS,
5129 true>
5130 worker(*this,
5131 src,
5132 dst,
5133 zero_dst_vector,
5134 *owning_class,
5135 function_pointer,
5136 nullptr,
5137 nullptr);
5138 task_info.loop(worker);
5139}
5140
5141
5142
5143template <int dim, typename Number, typename VectorizedArrayType>
5144template <typename CLASS, typename OutVector, typename InVector>
5145inline void
5147 void (CLASS::*function_pointer)(
5148 const MatrixFree<dim, Number, VectorizedArrayType> &,
5149 OutVector &,
5150 const InVector &,
5151 const std::pair<unsigned int, unsigned int> &) const,
5152 const CLASS *owning_class,
5153 OutVector &dst,
5154 const InVector &src,
5155 const std::function<void(const unsigned int, const unsigned int)>
5156 &operation_before_loop,
5157 const std::function<void(const unsigned int, const unsigned int)>
5158 &operation_after_loop,
5159 const unsigned int dof_handler_index_pre_post) const
5160{
5161 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5162 InVector,
5163 OutVector,
5164 CLASS,
5165 true>
5166 worker(*this,
5167 src,
5168 dst,
5169 false,
5170 *owning_class,
5171 function_pointer,
5172 nullptr,
5173 nullptr,
5176 operation_before_loop,
5177 operation_after_loop,
5178 dof_handler_index_pre_post);
5179 task_info.loop(worker);
5180}
5181
5182
5183
5184template <int dim, typename Number, typename VectorizedArrayType>
5185template <typename CLASS, typename OutVector, typename InVector>
5186inline void
5188 void (CLASS::*cell_operation)(
5189 const MatrixFree<dim, Number, VectorizedArrayType> &,
5190 OutVector &,
5191 const InVector &,
5192 const std::pair<unsigned int, unsigned int> &) const,
5193 void (CLASS::*inner_face_operation)(
5194 const MatrixFree<dim, Number, VectorizedArrayType> &,
5195 OutVector &,
5196 const InVector &,
5197 const std::pair<unsigned int, unsigned int> &) const,
5198 void (CLASS::*boundary_face_operation)(
5199 const MatrixFree<dim, Number, VectorizedArrayType> &,
5200 OutVector &,
5201 const InVector &,
5202 const std::pair<unsigned int, unsigned int> &) const,
5203 const CLASS *owning_class,
5204 OutVector &dst,
5205 const InVector &src,
5206 const bool zero_dst_vector,
5207 const DataAccessOnFaces dst_vector_face_access,
5208 const DataAccessOnFaces src_vector_face_access) const
5209{
5210 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5211 InVector,
5212 OutVector,
5213 CLASS,
5214 true>
5215 worker(*this,
5216 src,
5217 dst,
5218 zero_dst_vector,
5219 *owning_class,
5220 cell_operation,
5221 inner_face_operation,
5222 boundary_face_operation,
5223 src_vector_face_access,
5224 dst_vector_face_access);
5225 task_info.loop(worker);
5226}
5227
5228
5229
5230template <int dim, typename Number, typename VectorizedArrayType>
5231template <typename CLASS, typename OutVector, typename InVector>
5232inline void
5234 void (CLASS::*function_pointer)(
5235 const MatrixFree<dim, Number, VectorizedArrayType> &,
5236 OutVector &,
5237 const InVector &,
5238 const std::pair<unsigned int, unsigned int> &),
5239 CLASS *owning_class,
5240 OutVector &dst,
5241 const InVector &src,
5242 const bool zero_dst_vector) const
5243{
5244 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5245 InVector,
5246 OutVector,
5247 CLASS,
5248 false>
5249 worker(*this,
5250 src,
5251 dst,
5252 zero_dst_vector,
5253 *owning_class,
5254 function_pointer,
5255 nullptr,
5256 nullptr);
5257 task_info.loop(worker);
5258}
5259
5260
5261
5262template <int dim, typename Number, typename VectorizedArrayType>
5263template <typename CLASS, typename OutVector, typename InVector>
5264inline void
5266 void (CLASS::*function_pointer)(
5267 const MatrixFree<dim, Number, VectorizedArrayType> &,
5268 OutVector &,
5269 const InVector &,
5270 const std::pair<unsigned int, unsigned int> &),
5271 CLASS *owning_class,
5272 OutVector &dst,
5273 const InVector &src,
5274 const std::function<void(const unsigned int, const unsigned int)>
5275 &operation_before_loop,
5276 const std::function<void(const unsigned int, const unsigned int)>
5277 &operation_after_loop,
5278 const unsigned int dof_handler_index_pre_post) const
5279{
5280 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5281 InVector,
5282 OutVector,
5283 CLASS,
5284 false>
5285 worker(*this,
5286 src,
5287 dst,
5288 false,
5289 *owning_class,
5290 function_pointer,
5291 nullptr,
5292 nullptr,
5295 operation_before_loop,
5296 operation_after_loop,
5297 dof_handler_index_pre_post);
5298 task_info.loop(worker);
5299}
5300
5301
5302
5303template <int dim, typename Number, typename VectorizedArrayType>
5304template <typename CLASS, typename OutVector, typename InVector>
5305inline void
5307 void (CLASS::*cell_operation)(
5308 const MatrixFree<dim, Number, VectorizedArrayType> &,
5309 OutVector &,
5310 const InVector &,
5311 const std::pair<unsigned int, unsigned int> &),
5312 void (CLASS::*inner_face_operation)(
5313 const MatrixFree<dim, Number, VectorizedArrayType> &,
5314 OutVector &,
5315 const InVector &,
5316 const std::pair<unsigned int, unsigned int> &),
5317 void (CLASS::*boundary_face_operation)(
5318 const MatrixFree<dim, Number, VectorizedArrayType> &,
5319 OutVector &,
5320 const InVector &,
5321 const std::pair<unsigned int, unsigned int> &),
5322 CLASS *owning_class,
5323 OutVector &dst,
5324 const InVector &src,
5325 const bool zero_dst_vector,
5326 const DataAccessOnFaces dst_vector_face_access,
5327 const DataAccessOnFaces src_vector_face_access) const
5328{
5329 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5330 InVector,
5331 OutVector,
5332 CLASS,
5333 false>
5334 worker(*this,
5335 src,
5336 dst,
5337 zero_dst_vector,
5338 *owning_class,
5339 cell_operation,
5340 inner_face_operation,
5341 boundary_face_operation,
5342 src_vector_face_access,
5343 dst_vector_face_access);
5344 task_info.loop(worker);
5345}
5346
5347
5348
5349template <int dim, typename Number, typename VectorizedArrayType>
5350template <typename OutVector, typename InVector>
5351inline void
5353 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5354 OutVector &,
5355 const InVector &,
5356 const std::pair<unsigned int, unsigned int> &)>
5357 &cell_operation,
5358 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5359 OutVector &,
5360 const InVector &,
5361 const std::pair<unsigned int, unsigned int> &)>
5362 &inner_face_operation,
5363 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5364 OutVector &,
5365 const InVector &,
5366 const std::pair<unsigned int, unsigned int> &)>
5367 &boundary_face_operation,
5368 OutVector &dst,
5369 const InVector &src,
5370 const std::function<void(const unsigned int, const unsigned int)>
5371 &operation_before_loop,
5372 const std::function<void(const unsigned int, const unsigned int)>
5373 &operation_after_loop,
5374 const unsigned int dof_handler_index_pre_post,
5375 const DataAccessOnFaces dst_vector_face_access,
5376 const DataAccessOnFaces src_vector_face_access) const
5377{
5378 using Wrapper =
5379 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5380 InVector,
5381 OutVector>;
5382 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5383 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5384 InVector,
5385 OutVector,
5386 Wrapper,
5387 true>
5388 worker(*this,
5389 src,
5390 dst,
5391 false,
5392 wrap,
5393 &Wrapper::cell_integrator,
5394 &Wrapper::face_integrator,
5395 &Wrapper::boundary_integrator,
5396 src_vector_face_access,
5397 dst_vector_face_access,
5398 operation_before_loop,
5399 operation_after_loop,
5400 dof_handler_index_pre_post);
5401
5402 task_info.loop(worker);
5403}
5404
5405
5406
5407template <int dim, typename Number, typename VectorizedArrayType>
5408template <typename CLASS, typename OutVector, typename InVector>
5409inline void
5411 void (CLASS::*cell_operation)(const MatrixFree &,
5412 OutVector &,
5413 const InVector &,
5414 const std::pair<unsigned int, unsigned int> &)
5415 const,
5416 void (CLASS::*inner_face_operation)(
5417 const MatrixFree &,
5418 OutVector &,
5419 const InVector &,
5420 const std::pair<unsigned int, unsigned int> &) const,
5421 void (CLASS::*boundary_face_operation)(
5422 const MatrixFree &,
5423 OutVector &,
5424 const InVector &,
5425 const std::pair<unsigned int, unsigned int> &) const,
5426 const CLASS *owning_class,
5427 OutVector &dst,
5428 const InVector &src,
5429 const std::function<void(const unsigned int, const unsigned int)>
5430 &operation_before_loop,
5431 const std::function<void(const unsigned int, const unsigned int)>
5432 &operation_after_loop,
5433 const unsigned int dof_handler_index_pre_post,
5434 const DataAccessOnFaces dst_vector_face_access,
5435 const DataAccessOnFaces src_vector_face_access) const
5436{
5437 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5438 InVector,
5439 OutVector,
5440 CLASS,
5441 true>
5442 worker(*this,
5443 src,
5444 dst,
5445 false,
5446 *owning_class,
5447 cell_operation,
5448 inner_face_operation,
5449 boundary_face_operation,
5450 src_vector_face_access,
5451 dst_vector_face_access,
5452 operation_before_loop,
5453 operation_after_loop,
5454 dof_handler_index_pre_post);
5455 task_info.loop(worker);
5456}
5457
5458
5459
5460template <int dim, typename Number, typename VectorizedArrayType>
5461template <typename CLASS, typename OutVector, typename InVector>
5462inline void
5464 void (CLASS::*cell_operation)(const MatrixFree &,
5465 OutVector &,
5466 const InVector &,
5467 const std::pair<unsigned int, unsigned int> &),
5468 void (CLASS::*inner_face_operation)(
5469 const MatrixFree &,
5470 OutVector &,
5471 const InVector &,
5472 const std::pair<unsigned int, unsigned int> &),
5473 void (CLASS::*boundary_face_operation)(
5474 const MatrixFree &,
5475 OutVector &,
5476 const InVector &,
5477 const std::pair<unsigned int, unsigned int> &),
5478 const CLASS *owning_class,
5479 OutVector &dst,
5480 const InVector &src,
5481 const std::function<void(const unsigned int, const unsigned int)>
5482 &operation_before_loop,
5483 const std::function<void(const unsigned int, const unsigned int)>
5484 &operation_after_loop,
5485 const unsigned int dof_handler_index_pre_post,
5486 const DataAccessOnFaces dst_vector_face_access,
5487 const DataAccessOnFaces src_vector_face_access) const
5488{
5489 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5490 InVector,
5491 OutVector,
5492 CLASS,
5493 false>
5494 worker(*this,
5495 src,
5496 dst,
5497 false,
5498 *owning_class,
5499 cell_operation,
5500 inner_face_operation,
5501 boundary_face_operation,
5502 src_vector_face_access,
5503 dst_vector_face_access,
5504 operation_before_loop,
5505 operation_after_loop,
5506 dof_handler_index_pre_post);
5507 task_info.loop(worker);
5508}
5509
5510
5511
5512template <int dim, typename Number, typename VectorizedArrayType>
5513template <typename CLASS, typename OutVector, typename InVector>
5514inline void
5516 void (CLASS::*function_pointer)(
5517 const MatrixFree<dim, Number, VectorizedArrayType> &,
5518 OutVector &,
5519 const InVector &,
5520 const std::pair<unsigned int, unsigned int> &) const,
5521 const CLASS *owning_class,
5522 OutVector &dst,
5523 const InVector &src,
5524 const bool zero_dst_vector,
5525 const DataAccessOnFaces src_vector_face_access) const
5526{
5527 auto src_vector_face_access_temp = src_vector_face_access;
5528 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5529 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5530 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5531 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5532
5533 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5534 InVector,
5535 OutVector,
5536 CLASS,
5537 true>
5538 worker(*this,
5539 src,
5540 dst,
5541 zero_dst_vector,
5542 *owning_class,
5543 function_pointer,
5544 nullptr,
5545 nullptr,
5546 src_vector_face_access_temp,
5548 task_info.loop(worker);
5549}
5550
5551
5552
5553template <int dim, typename Number, typename VectorizedArrayType>
5554template <typename CLASS, typename OutVector, typename InVector>
5555inline void
5557 void (CLASS::*function_pointer)(
5558 const MatrixFree<dim, Number, VectorizedArrayType> &,
5559 OutVector &,
5560 const InVector &,
5561 const std::pair<unsigned int, unsigned int> &),
5562 CLASS *owning_class,
5563 OutVector &dst,
5564 const InVector &src,
5565 const bool zero_dst_vector,
5566 const DataAccessOnFaces src_vector_face_access) const
5567{
5568 auto src_vector_face_access_temp = src_vector_face_access;
5569 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5570 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5571 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5572 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5573
5574 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5575 InVector,
5576 OutVector,
5577 CLASS,
5578 false>
5579 worker(*this,
5580 src,
5581 dst,
5582 zero_dst_vector,
5583 *owning_class,
5584 function_pointer,
5585 nullptr,
5586 nullptr,
5587 src_vector_face_access_temp,
5589 task_info.loop(worker);
5590}
5591
5592
5593
5594template <int dim, typename Number, typename VectorizedArrayType>
5595template <typename OutVector, typename InVector>
5596inline void
5598 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5599 OutVector &,
5600 const InVector &,
5601 const std::pair<unsigned int, unsigned int> &)>
5602 &cell_operation,
5603 OutVector &dst,
5604 const InVector &src,
5605 const bool zero_dst_vector,
5606 const DataAccessOnFaces src_vector_face_access) const
5607{
5608 auto src_vector_face_access_temp = src_vector_face_access;
5609 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5610 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5611 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5612 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5613
5614 using Wrapper =
5615 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5616 InVector,
5617 OutVector>;
5618 Wrapper wrap(cell_operation, nullptr, nullptr);
5619
5620 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5621 InVector,
5622 OutVector,
5623 Wrapper,
5624 true>
5625 worker(*this,
5626 src,
5627 dst,
5628 zero_dst_vector,
5629 wrap,
5630 &Wrapper::cell_integrator,
5631 &Wrapper::face_integrator,
5632 &Wrapper::boundary_integrator,
5633 src_vector_face_access_temp,
5634 DataAccessOnFaces::none);
5635 task_info.loop(worker);
5636}
5637
5638
5639#endif // ifndef DOXYGEN
5640
5641
5642
5644
5645#endif
iterator begin()
void resize(const size_type new_size)
void reinit(const size_type size, const bool omit_zeroing_entries=false)
Abstract base class for mapping classes.
Definition mapping.h:318
unsigned int n_active_fe_indices() const
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
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
internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > mapping_info
unsigned int get_cell_active_fe_index(const std::pair< unsigned int, unsigned int > range) const
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
unsigned int n_ghost_cell_batches() const
void initialize_dof_handlers(const std::vector< const DoFHandler< dim, dim > * > &dof_handlers, const AdditionalData &additional_data)
types::boundary_id get_boundary_id(const unsigned int face_batch_index) const
void print(std::ostream &out) const
void update_mapping(const std::shared_ptr< hp::MappingCollection< dim > > &mapping)
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
unsigned int n_inner_face_batches() const
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
void update_mapping(const Mapping< dim > &mapping)
const Quadrature< dim > & get_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_index) const
unsigned int get_mg_level() const
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
void print_memory_consumption(StreamType &out) const
bool mapping_initialized() const
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
void clear()
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)
unsigned int get_n_q_points_face(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
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
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
bool at_irregular_cell(const unsigned int cell_batch_index) const
~MatrixFree() override=default
const AffineConstraints< Number > & get_affine_constraints(const unsigned int dof_handler_index=0) const
void copy_from(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free_base)
void initialize_face_data_vector(AlignedVector< T > &vec) const
unsigned int get_dofs_per_cell(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int n_constraint_pool_entries() const
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
const Number * constraint_pool_begin(const unsigned int pool_index) const
AlignedVector< Number > * acquire_scratch_data_non_threadsafe() const
const IndexSet & get_locally_owned_set(const unsigned int dof_handler_index=0) const
unsigned int get_dofs_per_face(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
internal::MatrixFreeFunctions::TaskInfo task_info
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
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
bool mapping_is_initialized
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
void initialize_dof_vector(LinearAlgebra::distributed::Vector< Number2 > &vec, const unsigned int dof_handler_index=0) const
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(const MatrixFree< dim, Number, VectorizedArrayType > &other)
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
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
unsigned int mg_level
std::pair< unsigned int, unsigned int > get_face_range_category(const std::pair< unsigned int, unsigned int > face_batch_range) const
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
std::vector< std::pair< unsigned int, unsigned int > > cell_level_index
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner(const unsigned int dof_handler_index=0) const
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
std::vector< Number > constraint_pool_data
VectorizedArrayType vectorized_value_type
unsigned int n_cell_batches() const
const IndexSet & get_ghost_set(const unsigned int dof_handler_index=0) const
std::pair< int, int > get_cell_level_and_index(const unsigned int cell_batch_index, const unsigned int lane_index) const
bool indices_initialized() const
const Quadrature< dim - 1 > & get_face_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int get_cell_category(const unsigned int cell_batch_index) const
Threads::ThreadLocalStorage< std::list< std::pair< bool, AlignedVector< VectorizedArrayType > > > > scratch_pad
unsigned int get_matrix_free_cell_index(const typename Triangulation< dim >::cell_iterator &cell) const
Number value_type
std::size_t memory_consumption() const
std::pair< unsigned int, unsigned int > get_face_category(const unsigned int face_batch_index) const
unsigned int n_boundary_face_batches() const
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
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
std::list< std::pair< bool, AlignedVector< Number > > > scratch_pad_non_threadsafe
const Number * constraint_pool_end(const unsigned int pool_index) const
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
std::vector< SmartPointer< const DoFHandler< dim > > > dof_handlers
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
void initialize_indices(const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const AdditionalData &additional_data)
std::vector< SmartPointer< const AffineConstraints< Number > > > affine_constraints
unsigned int n_components() const
unsigned int n_physical_cells() const
void initialize_dof_vector(VectorType &vec, const unsigned int dof_handler_index=0) const
unsigned int get_n_q_points(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
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())
std::vector< unsigned int > mf_cell_indices
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
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
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
std::vector< internal::MatrixFreeFunctions::DoFInfo > dof_info
void initialize_cell_data_vector(AlignedVector< T > &vec) const
unsigned int n_ghost_inner_face_batches() const
void release_scratch_data_non_threadsafe(const AlignedVector< Number > *memory) const
void renumber_dofs(std::vector< types::global_dof_index > &renumbering, const unsigned int dof_handler_index=0)
unsigned int get_face_active_fe_index(const std::pair< unsigned int, unsigned int > range, const bool is_interior_face=true) const
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
bool indices_are_initialized
internal::MatrixFreeFunctions::FaceInfo< VectorizedArrayType::size()> face_info
void reinit(const MappingType &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< number2 > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
unsigned int get_cell_range_category(const std::pair< unsigned int, unsigned int > cell_batch_range) const
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())
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
std::vector< unsigned int > constraint_pool_row_index
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
const std::vector< unsigned int > & get_constrained_dofs(const unsigned int dof_handler_index=0) const
unsigned int cell_level_index_end_local
unsigned int n_base_elements(const unsigned int dof_handler_index) const
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< Number > > shape_info
static bool is_supported(const FiniteElement< dim, spacedim > &fe)
static constexpr unsigned int dimension
A class that provides a separate storage location on each thread that accesses the object.
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
unsigned int level
Definition grid_out.cc:4626
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ActiveSelector::cell_iterator cell_iterator
UpdateFlags
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_default
No update.
for(unsigned int j=best_vertex+1;j< vertices.size();++j) if(vertices_to_use[j])
bool job_supports_mpi()
Definition mpi.cc:697
unsigned int minimum_parallel_grain_size
Definition parallel.cc:33
const types::boundary_id invalid_boundary_id
Definition types.h:292
static const unsigned int invalid_unsigned_int
Definition types.h:220
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
const InputIterator OutputIterator const Function & function
Definition parallel.h:168
const Iterator & begin
Definition parallel.h:609
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
unsigned short int fe_index
Definition types.h:59
unsigned int boundary_id
Definition types.h:144
TasksParallelScheme tasks_parallel_scheme
UpdateFlags mapping_update_flags_inner_faces
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)
std::vector< unsigned int > cell_vectorization_category
AdditionalData & operator=(const AdditionalData &other)
UpdateFlags mapping_update_flags_boundary_faces
UpdateFlags mapping_update_flags_faces_by_cells
AdditionalData(const AdditionalData &other)
std::vector< unsigned int > cell_loop_pre_list_index
Definition dof_info.h:748
std::vector< unsigned int > cell_loop_post_list_index
Definition dof_info.h:761
std::vector< std::pair< unsigned int, unsigned int > > vector_zero_range_list
Definition dof_info.h:741
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition dof_info.h:592
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_pre_list
Definition dof_info.h:754
std::vector< unsigned int > vector_zero_range_list_index
Definition dof_info.h:736
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_post_list
Definition dof_info.h:767
void loop(MFWorkerInterface &worker) const
Definition task_info.cc:350