deal.II version GIT relicensing-3307-g81a1c05a67 2025-05-12 20:50:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
matrix_free.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2012 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
16#ifndef dealii_matrix_free_h
17#define dealii_matrix_free_h
18
19#include <deal.II/base/config.h>
20
28
30
31#include <deal.II/fe/fe.h>
32#include <deal.II/fe/mapping.h>
33
36
41
48
49#include <cstdlib>
50#include <limits>
51#include <list>
52#include <memory>
53
54
56
57
58
110template <int dim,
111 typename Number = double,
112 typename VectorizedArrayType = VectorizedArray<Number>>
114{
115 static_assert(
116 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
117 "Type of Number and of VectorizedArrayType do not match.");
118
119public:
124 using value_type = Number;
125 using vectorized_value_type = VectorizedArrayType;
126
130 static constexpr unsigned int dimension = dim;
131
185 {
190
218
224 const unsigned int tasks_block_size = 0,
230 const unsigned int mg_level = numbers::invalid_unsigned_int,
231 const bool store_plain_indices = true,
232 const bool initialize_indices = true,
233 const bool initialize_mapping = true,
234 const bool overlap_communication_computation = true,
235 const bool hold_all_faces_to_owned_cells = false,
236 const bool cell_vectorization_categories_strict = false,
237 const bool allow_ghosted_vectors_in_loops = true)
253 , store_ghost_cells(false)
254 , communicator_sm(MPI_COMM_SELF)
255 {}
256
283
288 operator=(const AdditionalData &other) = default;
289
327
337 unsigned int tasks_block_size;
338
353
373
393
421
430 unsigned int mg_level;
431
439
450
459
473
482
504 std::vector<unsigned int> cell_vectorization_category;
505
516
528
535
540 };
541
550
555
559 ~MatrixFree() override = default;
560
573 template <typename QuadratureType, typename number2, typename MappingType>
574 void
575 reinit(const MappingType &mapping,
576 const DoFHandler<dim> &dof_handler,
577 const AffineConstraints<number2> &constraint,
578 const QuadratureType &quad,
579 const AdditionalData &additional_data = AdditionalData());
580
602 template <typename QuadratureType, typename number2, typename MappingType>
603 void
604 reinit(const MappingType &mapping,
605 const std::vector<const DoFHandler<dim> *> &dof_handler,
606 const std::vector<const AffineConstraints<number2> *> &constraint,
607 const std::vector<QuadratureType> &quad,
608 const AdditionalData &additional_data = AdditionalData());
609
617 template <typename QuadratureType, typename number2, typename MappingType>
618 void
619 reinit(const MappingType &mapping,
620 const std::vector<const DoFHandler<dim> *> &dof_handler,
621 const std::vector<const AffineConstraints<number2> *> &constraint,
622 const QuadratureType &quad,
623 const AdditionalData &additional_data = AdditionalData());
624
630 void
632 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free_base);
633
643 void
645
649 void
650 update_mapping(const std::shared_ptr<hp::MappingCollection<dim>> &mapping);
651
656 void
658
673 {
680 none,
681
692 values,
693
703
714 gradients,
715
725
733 };
734
779 template <typename OutVector, typename InVector>
780 void
781 cell_loop(const std::function<void(
783 OutVector &,
784 const InVector &,
785 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
786 OutVector &dst,
787 const InVector &src,
788 const bool zero_dst_vector = false) const;
789
836 template <typename CLASS, typename OutVector, typename InVector>
837 void
838 cell_loop(void (CLASS::*cell_operation)(
839 const MatrixFree &,
840 OutVector &,
841 const InVector &,
842 const std::pair<unsigned int, unsigned int> &) const,
843 const CLASS *owning_class,
844 OutVector &dst,
845 const InVector &src,
846 const bool zero_dst_vector = false) const;
847
851 template <typename CLASS, typename OutVector, typename InVector>
852 void
853 cell_loop(void (CLASS::*cell_operation)(
854 const MatrixFree &,
855 OutVector &,
856 const InVector &,
857 const std::pair<unsigned int, unsigned int> &),
858 CLASS *owning_class,
859 OutVector &dst,
860 const InVector &src,
861 const bool zero_dst_vector = false) const;
862
947 template <typename CLASS, typename OutVector, typename InVector>
948 void
949 cell_loop(void (CLASS::*cell_operation)(
950 const MatrixFree &,
951 OutVector &,
952 const InVector &,
953 const std::pair<unsigned int, unsigned int> &) const,
954 const CLASS *owning_class,
955 OutVector &dst,
956 const InVector &src,
957 const std::function<void(const unsigned int, const unsigned int)>
958 &operation_before_loop,
959 const std::function<void(const unsigned int, const unsigned int)>
960 &operation_after_loop,
961 const unsigned int dof_handler_index_pre_post = 0) const;
962
966 template <typename CLASS, typename OutVector, typename InVector>
967 void
968 cell_loop(void (CLASS::*cell_operation)(
969 const MatrixFree &,
970 OutVector &,
971 const InVector &,
972 const std::pair<unsigned int, unsigned int> &),
973 CLASS *owning_class,
974 OutVector &dst,
975 const InVector &src,
976 const std::function<void(const unsigned int, const unsigned int)>
977 &operation_before_loop,
978 const std::function<void(const unsigned int, const unsigned int)>
979 &operation_after_loop,
980 const unsigned int dof_handler_index_pre_post = 0) const;
981
986 template <typename OutVector, typename InVector>
987 void
988 cell_loop(const std::function<void(
990 OutVector &,
991 const InVector &,
992 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
993 OutVector &dst,
994 const InVector &src,
995 const std::function<void(const unsigned int, const unsigned int)>
996 &operation_before_loop,
997 const std::function<void(const unsigned int, const unsigned int)>
998 &operation_after_loop,
999 const unsigned int dof_handler_index_pre_post = 0) const;
1000
1076 template <typename OutVector, typename InVector>
1077 void
1079 const std::function<
1081 OutVector &,
1082 const InVector &,
1083 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1084 const std::function<void(
1086 OutVector &,
1087 const InVector &,
1088 const std::pair<unsigned int, unsigned int> &)> &inner_face_operation,
1089 const std::function<void(
1091 OutVector &,
1092 const InVector &,
1093 const std::pair<unsigned int, unsigned int> &)> &boundary_face_operation,
1094 OutVector &dst,
1095 const InVector &src,
1096 const bool zero_dst_vector = false,
1097 const DataAccessOnFaces dst_vector_face_access =
1099 const DataAccessOnFaces src_vector_face_access =
1101
1190 template <typename CLASS, typename OutVector, typename InVector>
1191 void
1192 loop(void (CLASS::*cell_operation)(
1193 const MatrixFree &,
1194 OutVector &,
1195 const InVector &,
1196 const std::pair<unsigned int, unsigned int> &) const,
1197 void (CLASS::*inner_face_operation)(
1198 const MatrixFree &,
1199 OutVector &,
1200 const InVector &,
1201 const std::pair<unsigned int, unsigned int> &) const,
1202 void (CLASS::*boundary_face_operation)(
1203 const MatrixFree &,
1204 OutVector &,
1205 const InVector &,
1206 const std::pair<unsigned int, unsigned int> &) const,
1207 const CLASS *owning_class,
1208 OutVector &dst,
1209 const InVector &src,
1210 const bool zero_dst_vector = false,
1211 const DataAccessOnFaces dst_vector_face_access =
1213 const DataAccessOnFaces src_vector_face_access =
1215
1219 template <typename CLASS, typename OutVector, typename InVector>
1220 void
1221 loop(void (CLASS::*cell_operation)(
1222 const MatrixFree &,
1223 OutVector &,
1224 const InVector &,
1225 const std::pair<unsigned int, unsigned int> &),
1226 void (CLASS::*inner_face_operation)(
1227 const MatrixFree &,
1228 OutVector &,
1229 const InVector &,
1230 const std::pair<unsigned int, unsigned int> &),
1231 void (CLASS::*boundary_face_operation)(
1232 const MatrixFree &,
1233 OutVector &,
1234 const InVector &,
1235 const std::pair<unsigned int, unsigned int> &),
1236 CLASS *owning_class,
1237 OutVector &dst,
1238 const InVector &src,
1239 const bool zero_dst_vector = false,
1240 const DataAccessOnFaces dst_vector_face_access =
1242 const DataAccessOnFaces src_vector_face_access =
1244
1366 template <typename CLASS, typename OutVector, typename InVector>
1367 void
1368 loop(void (CLASS::*cell_operation)(
1369 const MatrixFree &,
1370 OutVector &,
1371 const InVector &,
1372 const std::pair<unsigned int, unsigned int> &) const,
1373 void (CLASS::*inner_face_operation)(
1374 const MatrixFree &,
1375 OutVector &,
1376 const InVector &,
1377 const std::pair<unsigned int, unsigned int> &) const,
1378 void (CLASS::*boundary_face_operation)(
1379 const MatrixFree &,
1380 OutVector &,
1381 const InVector &,
1382 const std::pair<unsigned int, unsigned int> &) const,
1383 const CLASS *owning_class,
1384 OutVector &dst,
1385 const InVector &src,
1386 const std::function<void(const unsigned int, const unsigned int)>
1387 &operation_before_loop,
1388 const std::function<void(const unsigned int, const unsigned int)>
1389 &operation_after_loop,
1390 const unsigned int dof_handler_index_pre_post = 0,
1391 const DataAccessOnFaces dst_vector_face_access =
1393 const DataAccessOnFaces src_vector_face_access =
1395
1399 template <typename CLASS, typename OutVector, typename InVector>
1400 void
1401 loop(void (CLASS::*cell_operation)(
1402 const MatrixFree &,
1403 OutVector &,
1404 const InVector &,
1405 const std::pair<unsigned int, unsigned int> &),
1406 void (CLASS::*inner_face_operation)(
1407 const MatrixFree &,
1408 OutVector &,
1409 const InVector &,
1410 const std::pair<unsigned int, unsigned int> &),
1411 void (CLASS::*boundary_face_operation)(
1412 const MatrixFree &,
1413 OutVector &,
1414 const InVector &,
1415 const std::pair<unsigned int, unsigned int> &),
1416 const CLASS *owning_class,
1417 OutVector &dst,
1418 const InVector &src,
1419 const std::function<void(const unsigned int, const unsigned int)>
1420 &operation_before_loop,
1421 const std::function<void(const unsigned int, const unsigned int)>
1422 &operation_after_loop,
1423 const unsigned int dof_handler_index_pre_post = 0,
1424 const DataAccessOnFaces dst_vector_face_access =
1426 const DataAccessOnFaces src_vector_face_access =
1428
1434 template <typename OutVector, typename InVector>
1435 void
1437 const std::function<
1439 OutVector &,
1440 const InVector &,
1441 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1442 const std::function<void(
1444 OutVector &,
1445 const InVector &,
1446 const std::pair<unsigned int, unsigned int> &)> &inner_face_operation,
1447 const std::function<void(
1449 OutVector &,
1450 const InVector &,
1451 const std::pair<unsigned int, unsigned int> &)> &boundary_face_operation,
1452 OutVector &dst,
1453 const InVector &src,
1454 const std::function<void(const unsigned int, const unsigned int)>
1455 &operation_before_loop,
1456 const std::function<void(const unsigned int, const unsigned int)>
1457 &operation_after_loop,
1458 const unsigned int dof_handler_index_pre_post = 0,
1459 const DataAccessOnFaces dst_vector_face_access =
1461 const DataAccessOnFaces src_vector_face_access =
1463
1528 template <typename CLASS, typename OutVector, typename InVector>
1529 void
1530 loop_cell_centric(void (CLASS::*cell_operation)(
1531 const MatrixFree &,
1532 OutVector &,
1533 const InVector &,
1534 const std::pair<unsigned int, unsigned int> &) const,
1535 const CLASS *owning_class,
1536 OutVector &dst,
1537 const InVector &src,
1538 const bool zero_dst_vector = false,
1539 const DataAccessOnFaces src_vector_face_access =
1541
1545 template <typename CLASS, typename OutVector, typename InVector>
1546 void
1547 loop_cell_centric(void (CLASS::*cell_operation)(
1548 const MatrixFree &,
1549 OutVector &,
1550 const InVector &,
1551 const std::pair<unsigned int, unsigned int> &),
1552 CLASS *owning_class,
1553 OutVector &dst,
1554 const InVector &src,
1555 const bool zero_dst_vector = false,
1556 const DataAccessOnFaces src_vector_face_access =
1558
1562 template <typename OutVector, typename InVector>
1563 void
1565 const std::function<void(const MatrixFree &,
1566 OutVector &,
1567 const InVector &,
1568 const std::pair<unsigned int, unsigned int> &)>
1569 &cell_operation,
1570 OutVector &dst,
1571 const InVector &src,
1572 const bool zero_dst_vector = false,
1573 const DataAccessOnFaces src_vector_face_access =
1575
1583 std::pair<unsigned int, unsigned int>
1584 create_cell_subrange_hp(const std::pair<unsigned int, unsigned int> &range,
1585 const unsigned int fe_degree,
1586 const unsigned int dof_handler_index = 0) const;
1587
1594 std::pair<unsigned int, unsigned int>
1596 const std::pair<unsigned int, unsigned int> &range,
1597 const unsigned int fe_index,
1598 const unsigned int dof_handler_index = 0) const;
1599
1603 unsigned int
1605
1609 unsigned int
1611 const std::pair<unsigned int, unsigned int> range,
1612 const unsigned int dof_handler_index = numbers::invalid_unsigned_int) const;
1613
1617 unsigned int
1619 const std::pair<unsigned int, unsigned int> range,
1620 const bool is_interior_face = true,
1621 const unsigned int dof_handler_index = numbers::invalid_unsigned_int) const;
1622
1634 template <typename T>
1635 void
1637
1643 template <typename T>
1644 void
1646
1660 template <typename VectorType>
1661 void
1662 initialize_dof_vector(VectorType &vec,
1663 const unsigned int dof_handler_index = 0) const;
1664
1682 template <typename Number2>
1683 void
1685 const unsigned int dof_handler_index = 0) const;
1686
1697 const std::shared_ptr<const Utilities::MPI::Partitioner> &
1698 get_vector_partitioner(const unsigned int dof_handler_index = 0) const;
1699
1703 const IndexSet &
1704 get_locally_owned_set(const unsigned int dof_handler_index = 0) const;
1705
1709 const IndexSet &
1710 get_ghost_set(const unsigned int dof_handler_index = 0) const;
1711
1721 const std::vector<unsigned int> &
1722 get_constrained_dofs(const unsigned int dof_handler_index = 0) const;
1723
1734 void
1735 renumber_dofs(std::vector<types::global_dof_index> &renumbering,
1736 const unsigned int dof_handler_index = 0);
1737
1747 template <int spacedim>
1748 static bool
1750
1754 unsigned int
1756
1761 unsigned int
1762 n_base_elements(const unsigned int dof_handler_index) const;
1763
1771 unsigned int
1773
1783 unsigned int
1785
1792 unsigned int
1794
1804 unsigned int
1806
1817 unsigned int
1819
1824 unsigned int
1826
1838 get_boundary_id(const unsigned int face_batch_index) const;
1839
1844 std::array<types::boundary_id, VectorizedArrayType::size()>
1845 get_faces_by_cells_boundary_id(const unsigned int cell_batch_index,
1846 const unsigned int face_number) const;
1847
1852 const DoFHandler<dim> &
1853 get_dof_handler(const unsigned int dof_handler_index = 0) const;
1854
1862 get_affine_constraints(const unsigned int dof_handler_index = 0) const;
1863
1877 get_cell_iterator(const unsigned int cell_batch_index,
1878 const unsigned int lane_index,
1879 const unsigned int dof_handler_index = 0) const;
1880
1886 std::pair<int, int>
1887 get_cell_level_and_index(const unsigned int cell_batch_index,
1888 const unsigned int lane_index) const;
1889
1896 unsigned int
1898 const typename Triangulation<dim>::cell_iterator &cell) const;
1899
1912 std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>
1913 get_face_iterator(const unsigned int face_batch_index,
1914 const unsigned int lane_index,
1915 const bool interior = true,
1916 const unsigned int fe_component = 0) const;
1917
1930 bool
1931 at_irregular_cell(const unsigned int cell_batch_index) const;
1932
1942 unsigned int
1943 n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const;
1944
1954 unsigned int
1955 n_active_entries_per_face_batch(const unsigned int face_batch_index) const;
1956
1960 unsigned int
1961 get_dofs_per_cell(const unsigned int dof_handler_index = 0,
1962 const unsigned int hp_active_fe_index = 0) const;
1963
1967 unsigned int
1968 get_n_q_points(const unsigned int quad_index = 0,
1969 const unsigned int hp_active_fe_index = 0) const;
1970
1975 unsigned int
1976 get_dofs_per_face(const unsigned int dof_handler_index = 0,
1977 const unsigned int hp_active_fe_index = 0) const;
1978
1983 unsigned int
1984 get_n_q_points_face(const unsigned int quad_index = 0,
1985 const unsigned int hp_active_fe_index = 0) const;
1986
1990 const Quadrature<dim> &
1991 get_quadrature(const unsigned int quad_index = 0,
1992 const unsigned int hp_active_fe_index = 0) const;
1993
1997 const Quadrature<dim - 1> &
1998 get_face_quadrature(const unsigned int quad_index = 0,
1999 const unsigned int hp_active_fe_index = 0) const;
2000
2013 unsigned int
2015 const std::pair<unsigned int, unsigned int> cell_batch_range,
2016 const unsigned int dof_handler_index = numbers::invalid_unsigned_int) const;
2017
2022 std::pair<unsigned int, unsigned int>
2024 const std::pair<unsigned int, unsigned int> face_batch_range,
2025 const unsigned int dof_handler_index = numbers::invalid_unsigned_int) const;
2026
2038 unsigned int
2040 const unsigned int cell_batch_index,
2041 const unsigned int dof_handler_index = numbers::invalid_unsigned_int) const;
2042
2047 std::pair<unsigned int, unsigned int>
2049 const unsigned int face_batch_index,
2050 const unsigned int dof_handler_index = numbers::invalid_unsigned_int) const;
2051
2055 bool
2057
2062 bool
2064
2069 unsigned int
2071
2076 std::size_t
2078
2083 template <typename StreamType>
2084 void
2085 print_memory_consumption(StreamType &out) const;
2086
2091 void
2092 print(std::ostream &out) const;
2093
2107
2108 /*
2109 * Return geometry-dependent information on the cells.
2110 */
2111 const internal::MatrixFreeFunctions::
2112 MappingInfo<dim, Number, VectorizedArrayType> &
2114
2119 get_dof_info(const unsigned int dof_handler_index_component = 0) const;
2120
2124 unsigned int
2126
2131 const Number *
2132 constraint_pool_begin(const unsigned int pool_index) const;
2133
2139 const Number *
2140 constraint_pool_end(const unsigned int pool_index) const;
2141
2146 get_shape_info(const unsigned int dof_handler_index_component = 0,
2147 const unsigned int quad_index = 0,
2148 const unsigned int fe_base_element = 0,
2149 const unsigned int hp_active_fe_index = 0,
2150 const unsigned int hp_active_quad_index = 0) const;
2151
2156 VectorizedArrayType::size()> &
2157 get_face_info(const unsigned int face_batch_index) const;
2158
2159
2167
2183
2187 void
2189
2201
2205 void
2207 const AlignedVector<Number> *memory) const;
2208
2211private:
2216 template <typename number2, int q_dim>
2217 void
2219 const std::shared_ptr<hp::MappingCollection<dim>> &mapping,
2220 const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
2221 const std::vector<const AffineConstraints<number2> *> &constraint,
2222 const std::vector<IndexSet> &locally_owned_set,
2223 const std::vector<hp::QCollection<q_dim>> &quad,
2224 const AdditionalData &additional_data);
2225
2232 template <typename number2>
2233 void
2235 const std::vector<const AffineConstraints<number2> *> &constraint,
2236 const std::vector<IndexSet> &locally_owned_set,
2237 const AdditionalData &additional_data);
2238
2242 void
2244 const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
2245 const AdditionalData &additional_data);
2246
2250 std::vector<ObserverPointer<const DoFHandler<dim>>> dof_handlers;
2251
2258 std::vector<ObserverPointer<const AffineConstraints<Number>>>
2260
2265 std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
2266
2273 std::vector<Number> constraint_pool_data;
2274
2279 std::vector<unsigned int> constraint_pool_row_index;
2280
2287
2292
2299 std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
2300
2305 std::vector<unsigned int> mf_cell_indices;
2306
2314
2321
2326 internal::MatrixFreeFunctions::FaceInfo<VectorizedArrayType::size()>
2328
2333
2338
2347 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2349
2354 mutable std::list<std::pair<bool, AlignedVector<Number>>>
2356
2360 unsigned int mg_level;
2361
2367};
2368
2369
2370
2371/*----------------------- Inline functions ----------------------------------*/
2372
2373#ifndef DOXYGEN
2374
2375
2376
2377template <int dim, typename Number, typename VectorizedArrayType>
2378template <typename T>
2379inline void
2381 AlignedVector<T> &vec) const
2382{
2383 vec.resize(this->n_cell_batches() + this->n_ghost_cell_batches());
2384}
2385
2386
2387
2388template <int dim, typename Number, typename VectorizedArrayType>
2389template <typename T>
2390inline void
2392 AlignedVector<T> &vec) const
2393{
2394 vec.resize(this->n_inner_face_batches() + this->n_boundary_face_batches() +
2395 this->n_ghost_inner_face_batches());
2396}
2397
2398
2399
2400template <int dim, typename Number, typename VectorizedArrayType>
2401template <typename VectorType>
2402inline void
2404 VectorType &vec,
2405 const unsigned int comp) const
2406{
2407 static_assert(IsBlockVector<VectorType>::value == false,
2408 "This function is not supported for block vectors.");
2409
2410 Assert(task_info.n_procs == 1,
2411 ExcMessage("This function can only be used in serial."));
2412
2413 AssertIndexRange(comp, n_components());
2414 vec.reinit(dof_info[comp].vector_partitioner->size());
2415}
2416
2417
2418
2419template <int dim, typename Number, typename VectorizedArrayType>
2420template <typename Number2>
2421inline void
2424 const unsigned int comp) const
2425{
2426 AssertIndexRange(comp, n_components());
2427 vec.reinit(dof_info[comp].vector_partitioner, task_info.communicator_sm);
2428}
2429
2430
2431
2432template <int dim, typename Number, typename VectorizedArrayType>
2433inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2435 const unsigned int comp) const
2436{
2437 AssertIndexRange(comp, n_components());
2438 return dof_info[comp].vector_partitioner;
2439}
2440
2441
2442
2443template <int dim, typename Number, typename VectorizedArrayType>
2444inline const std::vector<unsigned int> &
2446 const unsigned int comp) const
2447{
2448 AssertIndexRange(comp, n_components());
2449 return dof_info[comp].constrained_dofs;
2450}
2451
2452
2453
2454template <int dim, typename Number, typename VectorizedArrayType>
2455inline unsigned int
2457{
2458 AssertDimension(dof_handlers.size(), dof_info.size());
2459 return dof_handlers.size();
2460}
2461
2462
2463
2464template <int dim, typename Number, typename VectorizedArrayType>
2465inline unsigned int
2467 const unsigned int dof_no) const
2468{
2469 AssertDimension(dof_handlers.size(), dof_info.size());
2470 AssertIndexRange(dof_no, dof_handlers.size());
2471 return dof_handlers[dof_no]->get_fe().n_base_elements();
2472}
2473
2474
2475
2476template <int dim, typename Number, typename VectorizedArrayType>
2479{
2480 return task_info;
2481}
2482
2483
2484
2485template <int dim, typename Number, typename VectorizedArrayType>
2486inline unsigned int
2488{
2489 return task_info.n_active_cells;
2490}
2491
2492
2493
2494template <int dim, typename Number, typename VectorizedArrayType>
2495inline unsigned int
2497{
2498 return *(task_info.cell_partition_data.end() - 2);
2499}
2500
2501
2502
2503template <int dim, typename Number, typename VectorizedArrayType>
2504inline unsigned int
2506{
2507 return *(task_info.cell_partition_data.end() - 1) -
2508 *(task_info.cell_partition_data.end() - 2);
2509}
2510
2511
2512
2513template <int dim, typename Number, typename VectorizedArrayType>
2514inline unsigned int
2516{
2517 if (task_info.face_partition_data.empty())
2518 return 0;
2519 return task_info.face_partition_data.back();
2520}
2521
2522
2523
2524template <int dim, typename Number, typename VectorizedArrayType>
2525inline unsigned int
2527{
2528 if (task_info.face_partition_data.empty())
2529 return 0;
2530 return task_info.boundary_partition_data.back() -
2531 task_info.face_partition_data.back();
2532}
2533
2534
2535
2536template <int dim, typename Number, typename VectorizedArrayType>
2537inline unsigned int
2539{
2540 if (task_info.face_partition_data.empty())
2541 return 0;
2542 return face_info.faces.size() - task_info.boundary_partition_data.back();
2543}
2544
2545
2546
2547template <int dim, typename Number, typename VectorizedArrayType>
2548inline types::boundary_id
2550 const unsigned int face_batch_index) const
2551{
2552 Assert(face_batch_index >= task_info.boundary_partition_data[0] &&
2553 face_batch_index < task_info.boundary_partition_data.back(),
2554 ExcIndexRange(face_batch_index,
2555 task_info.boundary_partition_data[0],
2556 task_info.boundary_partition_data.back()));
2557 return types::boundary_id(face_info.faces[face_batch_index].exterior_face_no);
2558}
2559
2560
2561
2562template <int dim, typename Number, typename VectorizedArrayType>
2563inline std::array<types::boundary_id, VectorizedArrayType::size()>
2565 const unsigned int cell_batch_index,
2566 const unsigned int face_number) const
2567{
2568 AssertIndexRange(cell_batch_index, n_cell_batches());
2569 AssertIndexRange(face_number, ReferenceCells::max_n_faces<dim>());
2570 Assert(face_info.cell_and_face_boundary_id.size(0) >= n_cell_batches(),
2572 std::array<types::boundary_id, VectorizedArrayType::size()> result;
2573 result.fill(numbers::invalid_boundary_id);
2574 for (unsigned int v = 0;
2575 v < n_active_entries_per_cell_batch(cell_batch_index);
2576 ++v)
2577 result[v] =
2578 face_info.cell_and_face_boundary_id(cell_batch_index, face_number, v);
2579 return result;
2580}
2581
2582
2583
2584template <int dim, typename Number, typename VectorizedArrayType>
2585inline const internal::MatrixFreeFunctions::
2586 MappingInfo<dim, Number, VectorizedArrayType> &
2588{
2589 return mapping_info;
2590}
2591
2592
2593
2594template <int dim, typename Number, typename VectorizedArrayType>
2597 const unsigned int dof_index) const
2598{
2599 AssertIndexRange(dof_index, n_components());
2600 return dof_info[dof_index];
2601}
2602
2603
2604
2605template <int dim, typename Number, typename VectorizedArrayType>
2606inline unsigned int
2608{
2609 return constraint_pool_row_index.size() - 1;
2610}
2611
2612
2613
2614template <int dim, typename Number, typename VectorizedArrayType>
2615inline const Number *
2617 const unsigned int row) const
2618{
2619 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2620 return constraint_pool_data.empty() ?
2621 nullptr :
2622 constraint_pool_data.data() + constraint_pool_row_index[row];
2623}
2624
2625
2626
2627template <int dim, typename Number, typename VectorizedArrayType>
2628inline const Number *
2630 const unsigned int row) const
2631{
2632 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2633 return constraint_pool_data.empty() ?
2634 nullptr :
2635 constraint_pool_data.data() + constraint_pool_row_index[row + 1];
2636}
2637
2638
2639
2640template <int dim, typename Number, typename VectorizedArrayType>
2641inline std::pair<unsigned int, unsigned int>
2643 const std::pair<unsigned int, unsigned int> &range,
2644 const unsigned int degree,
2645 const unsigned int dof_handler_component) const
2646{
2647 if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2648 {
2650 dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2652 dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2653 if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2654 return range;
2655 else
2656 return {range.second, range.second};
2657 }
2658
2659 const unsigned int fe_index =
2660 dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2661 if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2662 return {range.second, range.second};
2663 else
2664 return create_cell_subrange_hp_by_index(range,
2665 fe_index,
2666 dof_handler_component);
2667}
2668
2669
2670
2671template <int dim, typename Number, typename VectorizedArrayType>
2672inline bool
2674 const unsigned int cell_batch_index) const
2675{
2676 AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2677 return VectorizedArrayType::size() > 1 &&
2678 cell_level_index[(cell_batch_index + 1) * VectorizedArrayType::size() -
2679 1] == cell_level_index[(cell_batch_index + 1) *
2680 VectorizedArrayType::size() -
2681 2];
2682}
2683
2684
2685
2686template <int dim, typename Number, typename VectorizedArrayType>
2687unsigned int
2689{
2690 return shape_info.size(2);
2691}
2692
2693
2694template <int dim, typename Number, typename VectorizedArrayType>
2695unsigned int
2697 const std::pair<unsigned int, unsigned int> range,
2698 const unsigned int dof_handler_index) const
2699{
2700 const unsigned int dof_no =
2701 dof_handler_index == numbers::invalid_unsigned_int ?
2702 first_hp_dof_handler_index :
2703 dof_handler_index;
2704
2705 const auto &fe_indices = dof_info[dof_no].cell_active_fe_index;
2706
2707 if (fe_indices.empty() == true ||
2708 dof_handlers[dof_no]->get_fe_collection().size() == 1)
2709 return 0;
2710
2711 const auto index = fe_indices[range.first];
2712
2713 for (unsigned int i = range.first; i < range.second; ++i)
2714 AssertDimension(index, fe_indices[i]);
2715
2716 return index;
2717}
2718
2719
2720
2721template <int dim, typename Number, typename VectorizedArrayType>
2722unsigned int
2724 const std::pair<unsigned int, unsigned int> range,
2725 const bool is_interior_face,
2726 const unsigned int dof_handler_index) const
2727{
2728 const unsigned int dof_no =
2729 dof_handler_index == numbers::invalid_unsigned_int ?
2730 first_hp_dof_handler_index :
2731 dof_handler_index;
2732
2733 const auto &fe_indices = dof_info[dof_no].cell_active_fe_index;
2734
2735 if (fe_indices.empty() == true)
2736 return 0;
2737
2738 if (is_interior_face)
2739 {
2740 const unsigned int index =
2741 fe_indices[face_info.faces[range.first].cells_interior[0] /
2742 VectorizedArrayType::size()];
2743
2744 for (unsigned int i = range.first; i < range.second; ++i)
2745 AssertDimension(index,
2746 fe_indices[face_info.faces[i].cells_interior[0] /
2747 VectorizedArrayType::size()]);
2748
2749 return index;
2750 }
2751 else
2752 {
2753 const unsigned int index =
2754 fe_indices[face_info.faces[range.first].cells_exterior[0] /
2755 VectorizedArrayType::size()];
2756
2757 for (unsigned int i = range.first; i < range.second; ++i)
2758 AssertDimension(index,
2759 fe_indices[face_info.faces[i].cells_exterior[0] /
2760 VectorizedArrayType::size()]);
2761
2762 return index;
2763 }
2764}
2765
2766
2767
2768template <int dim, typename Number, typename VectorizedArrayType>
2769inline unsigned int
2771 const unsigned int cell_batch_index) const
2772{
2773 Assert(!dof_info.empty(), ExcNotInitialized());
2774 AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2775 const std::vector<unsigned char> &n_lanes_filled =
2776 dof_info[0].n_vectorization_lanes_filled
2778 AssertIndexRange(cell_batch_index, n_lanes_filled.size());
2779
2780 return n_lanes_filled[cell_batch_index];
2781}
2782
2783
2784
2785template <int dim, typename Number, typename VectorizedArrayType>
2786inline unsigned int
2788 const unsigned int face_batch_index) const
2789{
2790 AssertIndexRange(face_batch_index, face_info.faces.size());
2791 Assert(!dof_info.empty(), ExcNotInitialized());
2792 const std::vector<unsigned char> &n_lanes_filled =
2793 dof_info[0].n_vectorization_lanes_filled
2795 AssertIndexRange(face_batch_index, n_lanes_filled.size());
2796 return n_lanes_filled[face_batch_index];
2797}
2798
2799
2800
2801template <int dim, typename Number, typename VectorizedArrayType>
2802inline unsigned int
2804 const unsigned int dof_handler_index,
2805 const unsigned int active_fe_index) const
2806{
2807 return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2808}
2809
2810
2811
2812template <int dim, typename Number, typename VectorizedArrayType>
2813inline unsigned int
2815 const unsigned int quad_index,
2816 const unsigned int active_fe_index) const
2817{
2818 AssertIndexRange(quad_index, mapping_info.cell_data.size());
2819 return mapping_info.cell_data[quad_index]
2820 .descriptor[active_fe_index]
2821 .n_q_points;
2822}
2823
2824
2825
2826template <int dim, typename Number, typename VectorizedArrayType>
2827inline unsigned int
2829 const unsigned int dof_handler_index,
2830 const unsigned int active_fe_index) const
2831{
2832 return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2833}
2834
2835
2836
2837template <int dim, typename Number, typename VectorizedArrayType>
2838inline unsigned int
2840 const unsigned int quad_index,
2841 const unsigned int active_fe_index) const
2842{
2843 AssertIndexRange(quad_index, mapping_info.face_data.size());
2844 return mapping_info.face_data[quad_index]
2845 .descriptor[active_fe_index]
2846 .n_q_points;
2847}
2848
2849
2850
2851template <int dim, typename Number, typename VectorizedArrayType>
2852inline const IndexSet &
2854 const unsigned int dof_handler_index) const
2855{
2856 return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2857}
2858
2859
2860
2861template <int dim, typename Number, typename VectorizedArrayType>
2862inline const IndexSet &
2864 const unsigned int dof_handler_index) const
2865{
2866 return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2867}
2868
2869
2870
2871template <int dim, typename Number, typename VectorizedArrayType>
2874 const unsigned int dof_handler_index,
2875 const unsigned int index_quad,
2876 const unsigned int index_fe,
2877 const unsigned int active_fe_index,
2878 const unsigned int active_quad_index) const
2879{
2880 AssertIndexRange(dof_handler_index, dof_info.size());
2881 const unsigned int ind =
2882 dof_info[dof_handler_index].global_base_element_offset + index_fe;
2883 AssertIndexRange(ind, shape_info.size(0));
2884 AssertIndexRange(index_quad, shape_info.size(1));
2885 AssertIndexRange(active_fe_index, shape_info.size(2));
2886 AssertIndexRange(active_quad_index, shape_info.size(3));
2887 return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2888}
2889
2890
2891
2892template <int dim, typename Number, typename VectorizedArrayType>
2894 VectorizedArrayType::size()> &
2896 const unsigned int face_batch_index) const
2897{
2898 AssertIndexRange(face_batch_index, face_info.faces.size());
2899 return face_info.faces[face_batch_index];
2900}
2901
2902
2903
2904template <int dim, typename Number, typename VectorizedArrayType>
2905inline const Table<3, unsigned int> &
2907 const
2908{
2909 return face_info.cell_and_face_to_plain_faces;
2910}
2911
2912
2913
2914template <int dim, typename Number, typename VectorizedArrayType>
2915inline const Quadrature<dim> &
2917 const unsigned int quad_index,
2918 const unsigned int active_fe_index) const
2919{
2920 AssertIndexRange(quad_index, mapping_info.cell_data.size());
2921 return mapping_info.cell_data[quad_index]
2922 .descriptor[active_fe_index]
2923 .quadrature;
2924}
2925
2926
2927
2928template <int dim, typename Number, typename VectorizedArrayType>
2929inline const Quadrature<dim - 1> &
2931 const unsigned int quad_index,
2932 const unsigned int active_fe_index) const
2933{
2934 AssertIndexRange(quad_index, mapping_info.face_data.size());
2935 return mapping_info.face_data[quad_index]
2936 .descriptor[active_fe_index]
2937 .quadrature;
2938}
2939
2940
2941
2942template <int dim, typename Number, typename VectorizedArrayType>
2943inline unsigned int
2945 const std::pair<unsigned int, unsigned int> range,
2946 const unsigned int dof_handler_index) const
2947{
2948 auto result = get_cell_category(range.first, dof_handler_index);
2949
2950 for (unsigned int i = range.first; i < range.second; ++i)
2951 result = std::max(result, get_cell_category(i, dof_handler_index));
2952
2953 return result;
2954}
2955
2956
2957
2958template <int dim, typename Number, typename VectorizedArrayType>
2959inline std::pair<unsigned int, unsigned int>
2961 const std::pair<unsigned int, unsigned int> range,
2962 const unsigned int dof_handler_index) const
2963{
2964 auto result = get_face_category(range.first, dof_handler_index);
2965
2966 for (unsigned int i = range.first; i < range.second; ++i)
2967 {
2968 result.first =
2969 std::max(result.first, get_face_category(i, dof_handler_index).first);
2970 result.second =
2971 std::max(result.second, get_face_category(i, dof_handler_index).second);
2972 }
2973
2974 return result;
2975}
2976
2977
2978
2979template <int dim, typename Number, typename VectorizedArrayType>
2980inline unsigned int
2982 const unsigned int cell_batch_index,
2983 const unsigned int dof_handler_index) const
2984{
2985 AssertIndexRange(0, dof_info.size());
2986
2987 const unsigned int dof_no =
2988 dof_handler_index == numbers::invalid_unsigned_int ?
2989 first_hp_dof_handler_index :
2990 dof_handler_index;
2991
2992 AssertIndexRange(cell_batch_index,
2993 dof_info[dof_no].cell_active_fe_index.size());
2994 if (dof_info[dof_no].cell_active_fe_index.empty())
2995 return 0;
2996 else
2997 return dof_info[dof_no].cell_active_fe_index[cell_batch_index];
2998}
2999
3000
3001
3002template <int dim, typename Number, typename VectorizedArrayType>
3003inline std::pair<unsigned int, unsigned int>
3005 const unsigned int face_batch_index,
3006 const unsigned int dof_handler_index) const
3007{
3008 const unsigned int dof_no =
3009 dof_handler_index == numbers::invalid_unsigned_int ?
3010 first_hp_dof_handler_index :
3011 dof_handler_index;
3012
3013 AssertIndexRange(face_batch_index, face_info.faces.size());
3014 if (dof_info[dof_no].cell_active_fe_index.empty())
3015 return std::make_pair(0U, 0U);
3016
3017 std::pair<unsigned int, unsigned int> result = std::make_pair(0U, 0U);
3018 for (unsigned int v = 0;
3019 v < VectorizedArrayType::size() &&
3020 face_info.faces[face_batch_index].cells_interior[v] !=
3022 ++v)
3023 result.first = std::max(
3024 result.first,
3025 dof_info[dof_no].cell_active_fe_index[face_info.faces[face_batch_index]
3026 .cells_interior[v] /
3027 VectorizedArrayType::size()]);
3028 if (face_info.faces[face_batch_index].cells_exterior[0] !=
3030 for (unsigned int v = 0;
3031 v < VectorizedArrayType::size() &&
3032 face_info.faces[face_batch_index].cells_exterior[v] !=
3034 ++v)
3035 result.second = std::max(
3036 result.second,
3037 dof_info[dof_no].cell_active_fe_index[face_info.faces[face_batch_index]
3038 .cells_exterior[v] /
3039 VectorizedArrayType::size()]);
3040 else
3041 result.second = numbers::invalid_unsigned_int;
3042 return result;
3043}
3044
3045
3046
3047template <int dim, typename Number, typename VectorizedArrayType>
3048inline bool
3050{
3051 return indices_are_initialized;
3052}
3053
3054
3055
3056template <int dim, typename Number, typename VectorizedArrayType>
3057inline bool
3059{
3060 return mapping_is_initialized;
3061}
3062
3063
3064template <int dim, typename Number, typename VectorizedArrayType>
3065inline unsigned int
3067{
3068 return mg_level;
3069}
3070
3071
3072
3073template <int dim, typename Number, typename VectorizedArrayType>
3076{
3077 using list_type =
3078 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3079 list_type &data = scratch_pad.get();
3080 for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
3081 if (it->first == false)
3082 {
3083 it->first = true;
3084 return &it->second;
3085 }
3086 data.emplace_front(true, AlignedVector<VectorizedArrayType>());
3087 return &data.front().second;
3088}
3089
3090
3091
3092template <int dim, typename Number, typename VectorizedArrayType>
3093void
3095 const AlignedVector<VectorizedArrayType> *scratch) const
3096{
3097 using list_type =
3098 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3099 list_type &data = scratch_pad.get();
3100 for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
3101 if (&it->second == scratch)
3102 {
3103 Assert(it->first == true, ExcInternalError());
3104 it->first = false;
3105 return;
3106 }
3107 AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
3108}
3109
3110
3111
3112template <int dim, typename Number, typename VectorizedArrayType>
3116{
3117 for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
3118 scratch_pad_non_threadsafe.begin();
3119 it != scratch_pad_non_threadsafe.end();
3120 ++it)
3121 if (it->first == false)
3122 {
3123 it->first = true;
3124 return &it->second;
3125 }
3126 scratch_pad_non_threadsafe.push_front(
3127 std::make_pair(true, AlignedVector<Number>()));
3128 return &scratch_pad_non_threadsafe.front().second;
3129}
3130
3131
3132
3133template <int dim, typename Number, typename VectorizedArrayType>
3134void
3137 const AlignedVector<Number> *scratch) const
3138{
3139 for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
3140 scratch_pad_non_threadsafe.begin();
3141 it != scratch_pad_non_threadsafe.end();
3142 ++it)
3143 if (&it->second == scratch)
3144 {
3145 Assert(it->first == true, ExcInternalError());
3146 it->first = false;
3147 return;
3148 }
3149 AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
3150}
3151
3152
3153
3154// ------------------------------ reinit functions ---------------------------
3155
3156namespace internal
3157{
3158 namespace MatrixFreeImplementation
3159 {
3160 template <int dim, int spacedim>
3161 inline std::vector<IndexSet>
3162 extract_locally_owned_index_sets(
3163 const std::vector<const DoFHandler<dim, spacedim> *> &dofh,
3164 const unsigned int level)
3165 {
3166 std::vector<IndexSet> locally_owned_set;
3167 locally_owned_set.reserve(dofh.size());
3168 for (unsigned int j = 0; j < dofh.size(); ++j)
3170 locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
3171 else
3172 locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
3173 return locally_owned_set;
3174 }
3175 } // namespace MatrixFreeImplementation
3176} // namespace internal
3177
3178
3179
3180template <int dim, typename Number, typename VectorizedArrayType>
3181template <typename QuadratureType, typename number2, typename MappingType>
3182void
3184 const MappingType &mapping,
3185 const DoFHandler<dim> &dof_handler,
3186 const AffineConstraints<number2> &constraints_in,
3187 const QuadratureType &quad,
3189 &additional_data)
3190{
3191 std::vector<const DoFHandler<dim, dim> *> dof_handlers;
3192 std::vector<const AffineConstraints<number2> *> constraints;
3193
3194 dof_handlers.push_back(&dof_handler);
3195 constraints.push_back(&constraints_in);
3196
3197 std::vector<IndexSet> locally_owned_sets =
3198 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3199 dof_handlers, additional_data.mg_level);
3200
3201 std::vector<hp::QCollection<dim>> quad_hp;
3202 quad_hp.emplace_back(quad);
3203
3204 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3205 dof_handlers,
3206 constraints,
3207 locally_owned_sets,
3208 quad_hp,
3209 additional_data);
3210}
3211
3212
3213
3214template <int dim, typename Number, typename VectorizedArrayType>
3215template <typename QuadratureType, typename number2, typename MappingType>
3216void
3218 const MappingType &mapping,
3219 const std::vector<const DoFHandler<dim> *> &dof_handler,
3220 const std::vector<const AffineConstraints<number2> *> &constraint,
3221 const QuadratureType &quad,
3223 &additional_data)
3224{
3225 std::vector<IndexSet> locally_owned_set =
3226 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3227 dof_handler, additional_data.mg_level);
3228 std::vector<hp::QCollection<dim>> quad_hp;
3229 quad_hp.emplace_back(quad);
3230
3231 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3232 dof_handler,
3233 constraint,
3234 locally_owned_set,
3235 quad_hp,
3236 additional_data);
3237}
3238
3239
3240
3241template <int dim, typename Number, typename VectorizedArrayType>
3242template <typename QuadratureType, typename number2, typename MappingType>
3243void
3245 const MappingType &mapping,
3246 const std::vector<const DoFHandler<dim> *> &dof_handler,
3247 const std::vector<const AffineConstraints<number2> *> &constraint,
3248 const std::vector<QuadratureType> &quad,
3250 &additional_data)
3251{
3252 std::vector<IndexSet> locally_owned_set =
3253 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3254 dof_handler, additional_data.mg_level);
3255 std::vector<hp::QCollection<dim>> quad_hp;
3256 for (unsigned int q = 0; q < quad.size(); ++q)
3257 quad_hp.emplace_back(quad[q]);
3258
3259 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3260 dof_handler,
3261 constraint,
3262 locally_owned_set,
3263 quad_hp,
3264 additional_data);
3265}
3266
3267
3268
3269// ------------------------------ implementation of loops --------------------
3270
3271// internal helper functions that define how to call MPI data exchange
3272// functions: for generic vectors, do nothing at all. For distributed vectors,
3273// call update_ghost_values_start function and so on. If we have collections
3274// of vectors, just do the individual functions of the components. In order to
3275// keep ghost values consistent (whether we are in read or write mode), we
3276// also reset the values at the end. the whole situation is a bit complicated
3277// by the fact that we need to treat block vectors differently, which use some
3278// additional helper functions to select the blocks and template magic.
3279namespace internal
3280{
3284 template <int dim, typename Number, typename VectorizedArrayType>
3285 struct VectorDataExchange
3286 {
3287 // A shift for the MPI messages to reduce the risk for accidental
3288 // interaction with other open communications that a user program might
3289 // set up (parallel vectors support unfinished communication). We let
3290 // the other vectors use the first 20 assigned numbers and start the
3291 // matrix-free communication.
3292 static constexpr unsigned int channel_shift = 20;
3293
3294
3295
3300 VectorDataExchange(
3301 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3302 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3303 DataAccessOnFaces vector_face_access,
3304 const unsigned int n_components)
3305 : matrix_free(matrix_free)
3306 , vector_face_access(
3307 matrix_free.get_task_info().face_partition_data.empty() ?
3308 ::MatrixFree<dim, Number, VectorizedArrayType>::
3309 DataAccessOnFaces::unspecified :
3310 vector_face_access)
3311 , ghosts_were_set(false)
3312# ifdef DEAL_II_WITH_MPI
3313 , tmp_data(n_components)
3314 , requests(n_components)
3315# endif
3316 {
3317 (void)n_components;
3318 if (this->vector_face_access !=
3320 DataAccessOnFaces::unspecified)
3321 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3323 matrix_free.get_dof_info(c).vector_exchanger_face_variants.size(),
3324 5);
3325 }
3326
3327
3328
3332 ~VectorDataExchange() // NOLINT
3333 {
3334# ifdef DEAL_II_WITH_MPI
3335 for (unsigned int i = 0; i < tmp_data.size(); ++i)
3336 if (tmp_data[i] != nullptr)
3337 matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3338# endif
3339 }
3340
3341
3342
3347 template <typename VectorType>
3348 unsigned int
3349 find_vector_in_mf(const VectorType &vec,
3350 const bool check_global_compatibility = true) const
3351 {
3352 // case 1: vector was set up with MatrixFree::initialize_dof_vector()
3353 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3354 if (vec.get_partitioner().get() ==
3355 matrix_free.get_dof_info(c).vector_partitioner.get())
3356 return c;
3357
3358 // case 2: user provided own partitioner (compatibility mode)
3359 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3360 if (check_global_compatibility ?
3361 vec.get_partitioner()->is_globally_compatible(
3362 *matrix_free.get_dof_info(c).vector_partitioner) :
3363 vec.get_partitioner()->is_compatible(
3364 *matrix_free.get_dof_info(c).vector_partitioner))
3365 return c;
3366
3367 Assert(false,
3368 ExcNotImplemented("Could not find partitioner that fits vector"));
3369
3371 }
3372
3373
3374
3380 get_partitioner(const unsigned int mf_component) const
3381 {
3382 AssertDimension(matrix_free.get_dof_info(mf_component)
3383 .vector_exchanger_face_variants.size(),
3384 5);
3385 if (vector_face_access ==
3387 DataAccessOnFaces::none)
3388 return *matrix_free.get_dof_info(mf_component)
3389 .vector_exchanger_face_variants[0];
3390 else if (vector_face_access ==
3392 DataAccessOnFaces::values)
3393 return *matrix_free.get_dof_info(mf_component)
3394 .vector_exchanger_face_variants[1];
3395 else if (vector_face_access ==
3397 DataAccessOnFaces::gradients)
3398 return *matrix_free.get_dof_info(mf_component)
3399 .vector_exchanger_face_variants[2];
3400 else if (vector_face_access ==
3402 DataAccessOnFaces::values_all_faces)
3403 return *matrix_free.get_dof_info(mf_component)
3404 .vector_exchanger_face_variants[3];
3405 else if (vector_face_access ==
3407 DataAccessOnFaces::gradients_all_faces)
3408 return *matrix_free.get_dof_info(mf_component)
3409 .vector_exchanger_face_variants[4];
3410 else
3411 return *matrix_free.get_dof_info(mf_component).vector_exchanger.get();
3412 }
3413
3414
3415
3419 template <typename VectorType,
3420 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3421 * = nullptr>
3422 void
3423 update_ghost_values_start(const unsigned int /*component_in_block_vector*/,
3424 const VectorType & /*vec*/)
3425 {}
3426
3427
3432 template <typename VectorType,
3433 std::enable_if_t<!has_update_ghost_values_start<VectorType> &&
3434 !is_not_parallel_vector<VectorType>,
3435 VectorType> * = nullptr>
3436 void
3437 update_ghost_values_start(const unsigned int component_in_block_vector,
3438 const VectorType &vec)
3439 {
3440 (void)component_in_block_vector;
3441 const bool ghosts_set = vec.has_ghost_elements();
3442
3443 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3444 ghosts_set == false,
3446
3447 if (ghosts_set)
3448 {
3449 ghosts_were_set = true;
3450 return;
3451 }
3452
3453 vec.update_ghost_values();
3454 }
3455
3456
3457
3463 template <typename VectorType,
3464 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3465 !has_exchange_on_subset<VectorType>,
3466 VectorType> * = nullptr>
3467 void
3468 update_ghost_values_start(const unsigned int component_in_block_vector,
3469 const VectorType &vec)
3470 {
3471 (void)component_in_block_vector;
3472 const bool ghosts_set = vec.has_ghost_elements();
3473
3474 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3475 ghosts_set == false,
3477
3478 if (ghosts_set)
3479 {
3480 ghosts_were_set = true;
3481 return;
3482 }
3483
3484 vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3485 }
3486
3487
3488
3495 template <typename VectorType,
3496 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3497 has_exchange_on_subset<VectorType>,
3498 VectorType> * = nullptr>
3499 void
3500 update_ghost_values_start(const unsigned int component_in_block_vector,
3501 const VectorType &vec)
3502 {
3503 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3504 "Type mismatch between VectorType and VectorDataExchange");
3505 (void)component_in_block_vector;
3506 const bool ghosts_set = vec.has_ghost_elements();
3507
3508 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3509 ghosts_set == false,
3511
3512 if (ghosts_set)
3513 {
3514 ghosts_were_set = true;
3515 return;
3516 }
3517
3518 if (vec.size() != 0)
3519 {
3520# ifdef DEAL_II_WITH_MPI
3521 const unsigned int mf_component = find_vector_in_mf(vec);
3522
3523 const auto &part = get_partitioner(mf_component);
3524
3525 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3526 part.n_import_sm_procs() == 0)
3527 return;
3528
3529 tmp_data[component_in_block_vector] =
3530 matrix_free.acquire_scratch_data_non_threadsafe();
3531 tmp_data[component_in_block_vector]->resize_fast(
3532 part.n_import_indices());
3533 AssertDimension(requests.size(), tmp_data.size());
3534
3535 part.export_to_ghosted_array_start(
3536 component_in_block_vector * 2 + channel_shift,
3537 ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3538 vec.shared_vector_data(),
3539 ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3540 part.locally_owned_size(),
3541 matrix_free.get_dof_info(mf_component)
3542 .vector_partitioner->n_ghost_indices()),
3543 ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3544 part.n_import_indices()),
3545 this->requests[component_in_block_vector]);
3546# endif
3547 }
3548 }
3549
3550
3551
3556 template <typename VectorType,
3557 std::enable_if_t<!has_update_ghost_values_start<VectorType>,
3558 VectorType> * = nullptr>
3559 void
3560 update_ghost_values_finish(const unsigned int /*component_in_block_vector*/,
3561 const VectorType & /*vec*/)
3562 {}
3563
3564
3565
3571 template <typename VectorType,
3572 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3573 !has_exchange_on_subset<VectorType>,
3574 VectorType> * = nullptr>
3575 void
3576 update_ghost_values_finish(const unsigned int component_in_block_vector,
3577 const VectorType &vec)
3578 {
3579 (void)component_in_block_vector;
3580
3581 if (ghosts_were_set)
3582 return;
3583
3584 vec.update_ghost_values_finish();
3585 }
3586
3587
3588
3595 template <typename VectorType,
3596 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3597 has_exchange_on_subset<VectorType>,
3598 VectorType> * = nullptr>
3599 void
3600 update_ghost_values_finish(const unsigned int component_in_block_vector,
3601 const VectorType &vec)
3602 {
3603 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3604 "Type mismatch between VectorType and VectorDataExchange");
3605 (void)component_in_block_vector;
3606
3607 if (ghosts_were_set)
3608 return;
3609
3610 if (vec.size() != 0)
3611 {
3612# ifdef DEAL_II_WITH_MPI
3613 AssertIndexRange(component_in_block_vector, tmp_data.size());
3614 AssertDimension(requests.size(), tmp_data.size());
3615
3616 const unsigned int mf_component = find_vector_in_mf(vec);
3617
3618 const auto &part = get_partitioner(mf_component);
3619
3620 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3621 part.n_import_sm_procs() != 0)
3622 {
3623 part.export_to_ghosted_array_finish(
3624 ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3625 vec.shared_vector_data(),
3626 ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3627 part.locally_owned_size(),
3628 matrix_free.get_dof_info(mf_component)
3629 .vector_partitioner->n_ghost_indices()),
3630 this->requests[component_in_block_vector]);
3631
3632 matrix_free.release_scratch_data_non_threadsafe(
3633 tmp_data[component_in_block_vector]);
3634 tmp_data[component_in_block_vector] = nullptr;
3635 }
3636# endif
3637 }
3638 // let vector know that ghosts are being updated and we can read from
3639 // them
3640 vec.set_ghost_state(true);
3641 }
3642
3643
3644
3648 template <typename VectorType,
3649 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3650 * = nullptr>
3651 void
3652 compress_start(const unsigned int /*component_in_block_vector*/,
3653 VectorType & /*vec*/)
3654 {}
3655
3656
3657
3662 template <typename VectorType,
3663 std::enable_if_t<!has_compress_start<VectorType> &&
3664 !is_not_parallel_vector<VectorType>,
3665 VectorType> * = nullptr>
3666 void
3667 compress_start(const unsigned int component_in_block_vector,
3668 VectorType &vec)
3669 {
3670 (void)component_in_block_vector;
3671 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3672 vec.compress(VectorOperation::add);
3673 }
3674
3675
3676
3682 template <typename VectorType,
3683 std::enable_if_t<has_compress_start<VectorType> &&
3684 !has_exchange_on_subset<VectorType>,
3685 VectorType> * = nullptr>
3686 void
3687 compress_start(const unsigned int component_in_block_vector,
3688 VectorType &vec)
3689 {
3690 (void)component_in_block_vector;
3691 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3692 vec.compress_start(component_in_block_vector + channel_shift);
3693 }
3694
3695
3696
3703 template <typename VectorType,
3704 std::enable_if_t<has_compress_start<VectorType> &&
3705 has_exchange_on_subset<VectorType>,
3706 VectorType> * = nullptr>
3707 void
3708 compress_start(const unsigned int component_in_block_vector,
3709 VectorType &vec)
3710 {
3711 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3712 "Type mismatch between VectorType and VectorDataExchange");
3713 (void)component_in_block_vector;
3714 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3715
3716 if (vec.size() != 0)
3717 {
3718# ifdef DEAL_II_WITH_MPI
3719 const unsigned int mf_component = find_vector_in_mf(vec);
3720
3721 const auto &part = get_partitioner(mf_component);
3722
3723 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3724 part.n_import_sm_procs() == 0)
3725 return;
3726
3727 tmp_data[component_in_block_vector] =
3728 matrix_free.acquire_scratch_data_non_threadsafe();
3729 tmp_data[component_in_block_vector]->resize_fast(
3730 part.n_import_indices());
3731 AssertDimension(requests.size(), tmp_data.size());
3732
3733 part.import_from_ghosted_array_start(
3735 component_in_block_vector * 2 + channel_shift,
3736 ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3737 vec.shared_vector_data(),
3738 ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3739 matrix_free.get_dof_info(mf_component)
3740 .vector_partitioner->n_ghost_indices()),
3741 ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3742 part.n_import_indices()),
3743 this->requests[component_in_block_vector]);
3744# endif
3745 }
3746 }
3747
3748
3749
3754 template <
3755 typename VectorType,
3756 std::enable_if_t<!has_compress_start<VectorType>, VectorType> * = nullptr>
3757 void
3758 compress_finish(const unsigned int /*component_in_block_vector*/,
3759 VectorType & /*vec*/)
3760 {}
3761
3762
3763
3769 template <typename VectorType,
3770 std::enable_if_t<has_compress_start<VectorType> &&
3771 !has_exchange_on_subset<VectorType>,
3772 VectorType> * = nullptr>
3773 void
3774 compress_finish(const unsigned int component_in_block_vector,
3775 VectorType &vec)
3776 {
3777 (void)component_in_block_vector;
3778 vec.compress_finish(VectorOperation::add);
3779 }
3780
3781
3782
3789 template <typename VectorType,
3790 std::enable_if_t<has_compress_start<VectorType> &&
3791 has_exchange_on_subset<VectorType>,
3792 VectorType> * = nullptr>
3793 void
3794 compress_finish(const unsigned int component_in_block_vector,
3795 VectorType &vec)
3796 {
3797 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3798 "Type mismatch between VectorType and VectorDataExchange");
3799 (void)component_in_block_vector;
3800 if (vec.size() != 0)
3801 {
3802# ifdef DEAL_II_WITH_MPI
3803 AssertIndexRange(component_in_block_vector, tmp_data.size());
3804 AssertDimension(requests.size(), tmp_data.size());
3805
3806 const unsigned int mf_component = find_vector_in_mf(vec);
3807
3808 const auto &part = get_partitioner(mf_component);
3809
3810 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3811 part.n_import_sm_procs() != 0)
3812 {
3813 part.import_from_ghosted_array_finish(
3815 ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3816 vec.shared_vector_data(),
3817 ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3818 matrix_free.get_dof_info(mf_component)
3819 .vector_partitioner->n_ghost_indices()),
3821 tmp_data[component_in_block_vector]->begin(),
3822 part.n_import_indices()),
3823 this->requests[component_in_block_vector]);
3824
3825 matrix_free.release_scratch_data_non_threadsafe(
3826 tmp_data[component_in_block_vector]);
3827 tmp_data[component_in_block_vector] = nullptr;
3828 }
3829
3831 {
3832 const int ierr =
3833 MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3834 AssertThrowMPI(ierr);
3835 }
3836# endif
3837 }
3838 }
3839
3840
3841
3845 template <typename VectorType,
3846 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3847 * = nullptr>
3848 void
3849 reset_ghost_values(const VectorType & /*vec*/) const
3850 {}
3851
3852
3853
3858 template <typename VectorType,
3859 std::enable_if_t<!has_exchange_on_subset<VectorType> &&
3860 !is_not_parallel_vector<VectorType>,
3861 VectorType> * = nullptr>
3862 void
3863 reset_ghost_values(const VectorType &vec) const
3864 {
3865 if (ghosts_were_set == true)
3866 return;
3867
3868 vec.zero_out_ghost_values();
3869 }
3870
3871
3872
3878 template <typename VectorType,
3879 std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3880 * = nullptr>
3881 void
3882 reset_ghost_values(const VectorType &vec) const
3883 {
3884 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3885 "Type mismatch between VectorType and VectorDataExchange");
3886 if (ghosts_were_set == true)
3887 return;
3888
3889 if (vec.size() != 0)
3890 {
3891# ifdef DEAL_II_WITH_MPI
3892 AssertDimension(requests.size(), tmp_data.size());
3893
3894 const unsigned int mf_component = find_vector_in_mf(vec);
3895
3896 const auto &part = get_partitioner(mf_component);
3897
3898 if (part.n_ghost_indices() > 0)
3899 {
3900 part.reset_ghost_values(
3901 ArrayView<Number>(const_cast<VectorType &>(vec).begin() +
3902 part.locally_owned_size(),
3903 matrix_free.get_dof_info(mf_component)
3904 .vector_partitioner->n_ghost_indices()));
3905 }
3906
3907# endif
3908 }
3909 // let vector know that it's not ghosted anymore
3910 vec.set_ghost_state(false);
3911 }
3912
3913
3914
3920 template <typename VectorType,
3921 std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3922 * = nullptr>
3923 void
3924 zero_vector_region(const unsigned int range_index, VectorType &vec) const
3925 {
3926 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3927 "Type mismatch between VectorType and VectorDataExchange");
3928 if (range_index == numbers::invalid_unsigned_int)
3929 vec = Number();
3930 else
3931 {
3932 const unsigned int mf_component = find_vector_in_mf(vec, false);
3934 matrix_free.get_dof_info(mf_component);
3935 Assert(dof_info.vector_zero_range_list_index.empty() == false,
3937
3938 Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3940 AssertIndexRange(range_index,
3941 dof_info.vector_zero_range_list_index.size() - 1);
3942 for (unsigned int id =
3943 dof_info.vector_zero_range_list_index[range_index];
3944 id != dof_info.vector_zero_range_list_index[range_index + 1];
3945 ++id)
3946 std::memset(vec.begin() + dof_info.vector_zero_range_list[id].first,
3947 0,
3948 (dof_info.vector_zero_range_list[id].second -
3949 dof_info.vector_zero_range_list[id].first) *
3950 sizeof(Number));
3951 }
3952 }
3953
3954
3955
3961 template <typename VectorType,
3962 std::enable_if_t<!has_exchange_on_subset<VectorType>, VectorType>
3963 * = nullptr,
3964 std::enable_if_t<has_assignment_operator<VectorType>, VectorType>
3965 * = nullptr>
3966 void
3967 zero_vector_region(const unsigned int range_index, VectorType &vec) const
3968 {
3969 if (range_index == numbers::invalid_unsigned_int || range_index == 0)
3970 {
3971 if constexpr (std::is_same_v<
3973 VectorType>)
3974 {
3975 for (unsigned int i = 0; i < vec.size(); ++i)
3976 vec[i] = typename VectorType::value_type();
3977 }
3978 else
3979 vec = typename VectorType::value_type();
3980 }
3981 }
3982
3983
3984
3989 void
3990 zero_vector_region(const unsigned int, ...) const
3991 {
3992 Assert(false,
3994 "Zeroing is only implemented for vector types "
3995 "which provide operator=(const VectorType::value_type)"));
3996 }
3997
3998
3999
4000 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
4001 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
4002 DataAccessOnFaces vector_face_access;
4003 bool ghosts_were_set;
4004# ifdef DEAL_II_WITH_MPI
4005 std::vector<AlignedVector<Number> *> tmp_data;
4006 std::vector<std::vector<MPI_Request>> requests;
4007# endif
4008 }; // VectorDataExchange
4009
4010 template <typename VectorStruct>
4011 unsigned int
4012 n_components(const VectorStruct &vec);
4013
4014 template <typename VectorStruct>
4015 unsigned int
4016 n_components_block(const VectorStruct &vec, const std::bool_constant<true>)
4017 {
4018 unsigned int components = 0;
4019 for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
4020 components += n_components(vec.block(bl));
4021 return components;
4022 }
4023
4024 template <typename VectorStruct>
4025 unsigned int
4026 n_components_block(const VectorStruct &, const std::bool_constant<false>)
4027 {
4028 return 1;
4029 }
4030
4031 template <typename VectorStruct>
4032 unsigned int
4033 n_components(const VectorStruct &vec)
4034 {
4035 return n_components_block(
4036 vec, std::bool_constant<IsBlockVector<VectorStruct>::value>());
4037 }
4038
4039 template <typename VectorStruct>
4040 inline unsigned int
4041 n_components(const std::vector<VectorStruct> &vec)
4042 {
4043 unsigned int components = 0;
4044 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4045 components += n_components_block(
4046 vec[comp], std::bool_constant<IsBlockVector<VectorStruct>::value>());
4047 return components;
4048 }
4049
4050 template <typename VectorStruct>
4051 inline unsigned int
4052 n_components(const std::vector<VectorStruct *> &vec)
4053 {
4054 unsigned int components = 0;
4055 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4056 components += n_components_block(
4057 *vec[comp], std::bool_constant<IsBlockVector<VectorStruct>::value>());
4058 return components;
4059 }
4060
4061
4062
4063 // A helper function to identify block vectors with many components where we
4064 // should not try to overlap computations and communication because there
4065 // would be too many outstanding communication requests.
4066
4067 // default value for vectors that do not have communication_block_size
4068 template <typename VectorStruct,
4069 std::enable_if_t<!has_communication_block_size<VectorStruct>,
4070 VectorStruct> * = nullptr>
4071 constexpr unsigned int
4072 get_communication_block_size(const VectorStruct &)
4073 {
4075 }
4076
4077
4078
4079 template <typename VectorStruct,
4080 std::enable_if_t<has_communication_block_size<VectorStruct>,
4081 VectorStruct> * = nullptr>
4082 constexpr unsigned int
4083 get_communication_block_size(const VectorStruct &)
4084 {
4085 return VectorStruct::communication_block_size;
4086 }
4087
4088
4089
4090 template <typename VectorType,
4091 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType> * =
4092 nullptr>
4093 bool
4094 has_ghost_elements(const VectorType &vec)
4095 {
4096 (void)vec;
4097 return false;
4098 }
4099
4100
4101
4102 template <typename VectorType,
4103 std::enable_if_t<!is_not_parallel_vector<VectorType>, VectorType>
4104 * = nullptr>
4105 bool
4106 has_ghost_elements(const VectorType &vec)
4107 {
4108 return vec.has_ghost_elements();
4109 }
4110
4111
4112
4113 // --------------------------------------------------------------------------
4114 // below we have wrappers to distinguish between block and non-block vectors.
4115 // --------------------------------------------------------------------------
4116
4117 //
4118 // update_ghost_values_start
4119 //
4120
4121 // update_ghost_values for block vectors
4122 template <int dim,
4123 typename VectorStruct,
4124 typename Number,
4125 typename VectorizedArrayType,
4126 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4127 * = nullptr>
4128 void
4129 update_ghost_values_start(
4130 const VectorStruct &vec,
4131 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4132 const unsigned int channel = 0)
4133 {
4134 if (get_communication_block_size(vec) < vec.n_blocks())
4135 {
4136 const bool ghosts_set = vec.has_ghost_elements();
4137
4138 Assert(exchanger.matrix_free.get_task_info()
4139 .allow_ghosted_vectors_in_loops ||
4140 ghosts_set == false,
4142
4143 if (ghosts_set)
4144 {
4145 exchanger.ghosts_were_set = true;
4146 return;
4147 }
4148
4149 vec.update_ghost_values();
4150 }
4151 else
4152 {
4153 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4154 update_ghost_values_start(vec.block(i), exchanger, channel + i);
4155 }
4156 }
4157
4158
4159
4160 // update_ghost_values for non-block vectors
4161 template <int dim,
4162 typename VectorStruct,
4163 typename Number,
4164 typename VectorizedArrayType,
4165 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4166 * = nullptr>
4167 void
4168 update_ghost_values_start(
4169 const VectorStruct &vec,
4170 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4171 const unsigned int channel = 0)
4172 {
4173 exchanger.update_ghost_values_start(channel, vec);
4174 }
4175
4176
4177
4178 // update_ghost_values_start() for vector of vectors
4179 template <int dim,
4180 typename VectorStruct,
4181 typename Number,
4182 typename VectorizedArrayType>
4183 inline void
4184 update_ghost_values_start(
4185 const std::vector<VectorStruct> &vec,
4186 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4187 {
4188 unsigned int component_index = 0;
4189 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4190 {
4191 update_ghost_values_start(vec[comp], exchanger, component_index);
4192 component_index += n_components(vec[comp]);
4193 }
4194 }
4195
4196
4197
4198 // update_ghost_values_start() for vector of pointers to vectors
4199 template <int dim,
4200 typename VectorStruct,
4201 typename Number,
4202 typename VectorizedArrayType>
4203 inline void
4204 update_ghost_values_start(
4205 const std::vector<VectorStruct *> &vec,
4206 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4207 {
4208 unsigned int component_index = 0;
4209 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4210 {
4211 update_ghost_values_start(*vec[comp], exchanger, component_index);
4212 component_index += n_components(*vec[comp]);
4213 }
4214 }
4215
4216
4217
4218 //
4219 // update_ghost_values_finish
4220 //
4221
4222 // for block vectors
4223 template <int dim,
4224 typename VectorStruct,
4225 typename Number,
4226 typename VectorizedArrayType,
4227 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4228 * = nullptr>
4229 void
4230 update_ghost_values_finish(
4231 const VectorStruct &vec,
4232 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4233 const unsigned int channel = 0)
4234 {
4235 if (get_communication_block_size(vec) < vec.n_blocks())
4236 {
4237 // do nothing, everything has already been completed in the _start()
4238 // call
4239 }
4240 else
4241 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4242 update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4243 }
4244
4245
4246
4247 // for non-block vectors
4248 template <int dim,
4249 typename VectorStruct,
4250 typename Number,
4251 typename VectorizedArrayType,
4252 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4253 * = nullptr>
4254 void
4255 update_ghost_values_finish(
4256 const VectorStruct &vec,
4257 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4258 const unsigned int channel = 0)
4259 {
4260 exchanger.update_ghost_values_finish(channel, vec);
4261 }
4262
4263
4264
4265 // for vector of vectors
4266 template <int dim,
4267 typename VectorStruct,
4268 typename Number,
4269 typename VectorizedArrayType>
4270 inline void
4271 update_ghost_values_finish(
4272 const std::vector<VectorStruct> &vec,
4273 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4274 {
4275 unsigned int component_index = 0;
4276 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4277 {
4278 update_ghost_values_finish(vec[comp], exchanger, component_index);
4279 component_index += n_components(vec[comp]);
4280 }
4281 }
4282
4283
4284
4285 // for vector of pointers to vectors
4286 template <int dim,
4287 typename VectorStruct,
4288 typename Number,
4289 typename VectorizedArrayType>
4290 inline void
4291 update_ghost_values_finish(
4292 const std::vector<VectorStruct *> &vec,
4293 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4294 {
4295 unsigned int component_index = 0;
4296 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4297 {
4298 update_ghost_values_finish(*vec[comp], exchanger, component_index);
4299 component_index += n_components(*vec[comp]);
4300 }
4301 }
4302
4303
4304
4305 //
4306 // compress_start
4307 //
4308
4309 // for block vectors
4310 template <int dim,
4311 typename VectorStruct,
4312 typename Number,
4313 typename VectorizedArrayType,
4314 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4315 * = nullptr>
4316 inline void
4317 compress_start(
4318 VectorStruct &vec,
4319 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4320 const unsigned int channel = 0)
4321 {
4322 if (get_communication_block_size(vec) < vec.n_blocks())
4323 vec.compress(VectorOperation::add);
4324 else
4325 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4326 compress_start(vec.block(i), exchanger, channel + i);
4327 }
4328
4329
4330
4331 // for non-block vectors
4332 template <int dim,
4333 typename VectorStruct,
4334 typename Number,
4335 typename VectorizedArrayType,
4336 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4337 * = nullptr>
4338 inline void
4339 compress_start(
4340 VectorStruct &vec,
4341 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4342 const unsigned int channel = 0)
4343 {
4344 exchanger.compress_start(channel, vec);
4345 }
4346
4347
4348
4349 // for std::vector of vectors
4350 template <int dim,
4351 typename VectorStruct,
4352 typename Number,
4353 typename VectorizedArrayType>
4354 inline void
4355 compress_start(
4356 std::vector<VectorStruct> &vec,
4357 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4358 {
4359 unsigned int component_index = 0;
4360 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4361 {
4362 compress_start(vec[comp], exchanger, component_index);
4363 component_index += n_components(vec[comp]);
4364 }
4365 }
4366
4367
4368
4369 // for std::vector of pointer to vectors
4370 template <int dim,
4371 typename VectorStruct,
4372 typename Number,
4373 typename VectorizedArrayType>
4374 inline void
4375 compress_start(
4376 std::vector<VectorStruct *> &vec,
4377 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4378 {
4379 unsigned int component_index = 0;
4380 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4381 {
4382 compress_start(*vec[comp], exchanger, component_index);
4383 component_index += n_components(*vec[comp]);
4384 }
4385 }
4386
4387
4388
4389 //
4390 // compress_finish
4391 //
4392
4393 // for block vectors
4394 template <int dim,
4395 typename VectorStruct,
4396 typename Number,
4397 typename VectorizedArrayType,
4398 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4399 * = nullptr>
4400 inline void
4401 compress_finish(
4402 VectorStruct &vec,
4403 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4404 const unsigned int channel = 0)
4405 {
4406 if (get_communication_block_size(vec) < vec.n_blocks())
4407 {
4408 // do nothing, everything has already been completed in the _start()
4409 // call
4410 }
4411 else
4412 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4413 compress_finish(vec.block(i), exchanger, channel + i);
4414 }
4415
4416
4417
4418 // for non-block vectors
4419 template <int dim,
4420 typename VectorStruct,
4421 typename Number,
4422 typename VectorizedArrayType,
4423 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4424 * = nullptr>
4425 inline void
4426 compress_finish(
4427 VectorStruct &vec,
4428 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4429 const unsigned int channel = 0)
4430 {
4431 exchanger.compress_finish(channel, vec);
4432 }
4433
4434
4435
4436 // for std::vector of vectors
4437 template <int dim,
4438 typename VectorStruct,
4439 typename Number,
4440 typename VectorizedArrayType>
4441 inline void
4442 compress_finish(
4443 std::vector<VectorStruct> &vec,
4444 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4445 {
4446 unsigned int component_index = 0;
4447 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4448 {
4449 compress_finish(vec[comp], exchanger, component_index);
4450 component_index += n_components(vec[comp]);
4451 }
4452 }
4453
4454
4455
4456 // for std::vector of pointer to vectors
4457 template <int dim,
4458 typename VectorStruct,
4459 typename Number,
4460 typename VectorizedArrayType>
4461 inline void
4462 compress_finish(
4463 std::vector<VectorStruct *> &vec,
4464 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4465 {
4466 unsigned int component_index = 0;
4467 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4468 {
4469 compress_finish(*vec[comp], exchanger, component_index);
4470 component_index += n_components(*vec[comp]);
4471 }
4472 }
4473
4474
4475
4476 //
4477 // reset_ghost_values:
4478 //
4479 // if the input vector did not have ghosts imported, clear them here again
4480 // in order to avoid subsequent operations e.g. in linear solvers to work
4481 // with ghosts all the time
4482 //
4483
4484 // for block vectors
4485 template <int dim,
4486 typename VectorStruct,
4487 typename Number,
4488 typename VectorizedArrayType,
4489 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4490 * = nullptr>
4491 inline void
4492 reset_ghost_values(
4493 const VectorStruct &vec,
4494 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4495 {
4496 // return immediately if there is nothing to do.
4497 if (exchanger.ghosts_were_set == true)
4498 return;
4499
4500 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4501 reset_ghost_values(vec.block(i), exchanger);
4502 }
4503
4504
4505
4506 // for non-block vectors
4507 template <int dim,
4508 typename VectorStruct,
4509 typename Number,
4510 typename VectorizedArrayType,
4511 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4512 * = nullptr>
4513 inline void
4514 reset_ghost_values(
4515 const VectorStruct &vec,
4516 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4517 {
4518 // return immediately if there is nothing to do.
4519 if (exchanger.ghosts_were_set == true)
4520 return;
4521
4522 exchanger.reset_ghost_values(vec);
4523 }
4524
4525
4526
4527 // for std::vector of vectors
4528 template <int dim,
4529 typename VectorStruct,
4530 typename Number,
4531 typename VectorizedArrayType>
4532 inline void
4533 reset_ghost_values(
4534 const std::vector<VectorStruct> &vec,
4535 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4536 {
4537 // return immediately if there is nothing to do.
4538 if (exchanger.ghosts_were_set == true)
4539 return;
4540
4541 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4542 reset_ghost_values(vec[comp], exchanger);
4543 }
4544
4545
4546
4547 // for std::vector of pointer to vectors
4548 template <int dim,
4549 typename VectorStruct,
4550 typename Number,
4551 typename VectorizedArrayType>
4552 inline void
4553 reset_ghost_values(
4554 const std::vector<VectorStruct *> &vec,
4555 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4556 {
4557 // return immediately if there is nothing to do.
4558 if (exchanger.ghosts_were_set == true)
4559 return;
4560
4561 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4562 reset_ghost_values(*vec[comp], exchanger);
4563 }
4564
4565
4566
4567 //
4568 // zero_vector_region
4569 //
4570
4571 // for block vectors
4572 template <int dim,
4573 typename VectorStruct,
4574 typename Number,
4575 typename VectorizedArrayType,
4576 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4577 * = nullptr>
4578 inline void
4579 zero_vector_region(
4580 const unsigned int range_index,
4581 VectorStruct &vec,
4582 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4583 {
4584 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4585 exchanger.zero_vector_region(range_index, vec.block(i));
4586 }
4587
4588
4589
4590 // for non-block vectors
4591 template <int dim,
4592 typename VectorStruct,
4593 typename Number,
4594 typename VectorizedArrayType,
4595 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4596 * = nullptr>
4597 inline void
4598 zero_vector_region(
4599 const unsigned int range_index,
4600 VectorStruct &vec,
4601 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4602 {
4603 exchanger.zero_vector_region(range_index, vec);
4604 }
4605
4606
4607
4608 // for std::vector of vectors
4609 template <int dim,
4610 typename VectorStruct,
4611 typename Number,
4612 typename VectorizedArrayType>
4613 inline void
4614 zero_vector_region(
4615 const unsigned int range_index,
4616 std::vector<VectorStruct> &vec,
4617 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4618 {
4619 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4620 zero_vector_region(range_index, vec[comp], exchanger);
4621 }
4622
4623
4624
4625 // for std::vector of pointers to vectors
4626 template <int dim,
4627 typename VectorStruct,
4628 typename Number,
4629 typename VectorizedArrayType>
4630 inline void
4631 zero_vector_region(
4632 const unsigned int range_index,
4633 std::vector<VectorStruct *> &vec,
4634 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4635 {
4636 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4637 zero_vector_region(range_index, *vec[comp], exchanger);
4638 }
4639
4640
4641
4642 // Apply a unit matrix operation to constrained DoFs: Default cases where we
4643 // cannot detect a LinearAlgebra::distributed::Vector, we do not do
4644 // anything, else we apply the constraints as a unit operation
4645 template <typename VectorStruct1, typename VectorStruct2>
4646 inline void
4647 apply_operation_to_constrained_dofs(const std::vector<unsigned int> &,
4648 const VectorStruct1 &,
4649 VectorStruct2 &)
4650 {}
4651
4652 template <typename Number>
4653 inline void
4654 apply_operation_to_constrained_dofs(
4655 const std::vector<unsigned int> &constrained_dofs,
4658 {
4659 for (const unsigned int i : constrained_dofs)
4660 dst.local_element(i) = src.local_element(i);
4661 }
4662
4663
4664 namespace MatrixFreeFunctions
4665 {
4666 // struct to select between a const interface and a non-const interface
4667 // for MFWorker
4668 template <typename, typename, typename, typename, bool>
4669 struct InterfaceSelector
4670 {};
4671
4672 // Version for constant functions
4673 template <typename MF,
4674 typename InVector,
4675 typename OutVector,
4676 typename Container>
4677 struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4678 {
4679 using function_type = void (Container::*)(
4680 const MF &,
4681 OutVector &,
4682 const InVector &,
4683 const std::pair<unsigned int, unsigned int> &) const;
4684 };
4685
4686 // Version for non-constant functions
4687 template <typename MF,
4688 typename InVector,
4689 typename OutVector,
4690 typename Container>
4691 struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4692 {
4693 using function_type =
4694 void (Container::*)(const MF &,
4695 OutVector &,
4696 const InVector &,
4697 const std::pair<unsigned int, unsigned int> &);
4698 };
4699 } // namespace MatrixFreeFunctions
4700
4701
4702
4703 // A implementation class for the worker object that runs the various
4704 // operations we want to perform during the matrix-free loop
4705 template <typename MF,
4706 typename InVector,
4707 typename OutVector,
4708 typename Container,
4709 bool is_constant>
4710 class MFWorker : public MFWorkerInterface
4711 {
4712 public:
4713 // An alias to make the arguments further down more readable
4714 using function_type = typename MatrixFreeFunctions::
4715 InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4716 function_type;
4717
4718 // constructor, binds all the arguments to this class
4719 MFWorker(const MF &matrix_free,
4720 const InVector &src,
4721 OutVector &dst,
4722 const bool zero_dst_vector_setting,
4723 const Container &container,
4724 function_type cell_function,
4725 function_type face_function,
4726 function_type boundary_function,
4727 const typename MF::DataAccessOnFaces src_vector_face_access =
4728 MF::DataAccessOnFaces::none,
4729 const typename MF::DataAccessOnFaces dst_vector_face_access =
4730 MF::DataAccessOnFaces::none,
4731 const std::function<void(const unsigned int, const unsigned int)>
4732 &operation_before_loop = {},
4733 const std::function<void(const unsigned int, const unsigned int)>
4734 &operation_after_loop = {},
4735 const unsigned int dof_handler_index_pre_post = 0)
4736 : matrix_free(matrix_free)
4737 , container(const_cast<Container &>(container))
4738 , cell_function(cell_function)
4739 , face_function(face_function)
4740 , boundary_function(boundary_function)
4741 , src(src)
4742 , dst(dst)
4743 , src_data_exchanger(matrix_free,
4744 src_vector_face_access,
4745 n_components(src))
4746 , dst_data_exchanger(matrix_free,
4747 dst_vector_face_access,
4748 n_components(dst))
4749 , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4750 , zero_dst_vector_setting(zero_dst_vector_setting &&
4751 !src_and_dst_are_same)
4752 , operation_before_loop(operation_before_loop)
4753 , operation_after_loop(operation_after_loop)
4754 , dof_handler_index_pre_post(dof_handler_index_pre_post)
4755 {
4756 Assert(!has_ghost_elements(dst),
4757 ExcMessage("The destination vector passed to the matrix-free "
4758 "loop is ghosted. This is not allowed."));
4759 }
4760
4761 // Runs the cell work. If no function is given, nothing is done
4762 virtual void
4763 cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4764 {
4765 if (cell_function != nullptr && cell_range.second > cell_range.first)
4766 for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4767 {
4768 const auto cell_subrange =
4769 matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4770
4771 if (cell_subrange.second <= cell_subrange.first)
4772 continue;
4773
4774 (container.*
4775 cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4776 }
4777 }
4778
4779 virtual void
4780 cell(const unsigned int range_index) override
4781 {
4782 process_range(cell_function,
4783 matrix_free.get_task_info().cell_partition_data_hp_ptr,
4784 matrix_free.get_task_info().cell_partition_data_hp,
4785 range_index);
4786 }
4787
4788 virtual void
4789 face(const unsigned int range_index) override
4790 {
4791 process_range(face_function,
4792 matrix_free.get_task_info().face_partition_data_hp_ptr,
4793 matrix_free.get_task_info().face_partition_data_hp,
4794 range_index);
4795 }
4796
4797 virtual void
4798 boundary(const unsigned int range_index) override
4799 {
4800 process_range(boundary_function,
4801 matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4802 matrix_free.get_task_info().boundary_partition_data_hp,
4803 range_index);
4804 }
4805
4806 private:
4807 void
4808 process_range(const function_type &fu,
4809 const std::vector<unsigned int> &ptr,
4810 const std::vector<unsigned int> &data,
4811 const unsigned int range_index)
4812 {
4813 if (fu == nullptr)
4814 return;
4815
4816 AssertIndexRange(range_index + 1, ptr.size());
4817 for (unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4818 {
4819 AssertIndexRange(2 * i + 1, data.size());
4820 (container.*fu)(matrix_free,
4821 this->dst,
4822 this->src,
4823 std::make_pair(data[2 * i], data[2 * i + 1]));
4824 }
4825 }
4826
4827 public:
4828 // Starts the communication for the update ghost values operation. We
4829 // cannot call this update if ghost and destination are the same because
4830 // that would introduce spurious entries in the destination (there is also
4831 // the problem that reading from a vector that we also write to is usually
4832 // not intended in case there is overlap, but this is up to the
4833 // application code to decide and we cannot catch this case here).
4834 virtual void
4835 vector_update_ghosts_start() override
4836 {
4837 if (!src_and_dst_are_same)
4838 internal::update_ghost_values_start(src, src_data_exchanger);
4839 }
4840
4841 // Finishes the communication for the update ghost values operation
4842 virtual void
4843 vector_update_ghosts_finish() override
4844 {
4845 if (!src_and_dst_are_same)
4846 internal::update_ghost_values_finish(src, src_data_exchanger);
4847 }
4848
4849 // Starts the communication for the vector compress operation
4850 virtual void
4851 vector_compress_start() override
4852 {
4853 internal::compress_start(dst, dst_data_exchanger);
4854 }
4855
4856 // Finishes the communication for the vector compress operation
4857 virtual void
4858 vector_compress_finish() override
4859 {
4860 internal::compress_finish(dst, dst_data_exchanger);
4861 if (!src_and_dst_are_same)
4862 internal::reset_ghost_values(src, src_data_exchanger);
4863 }
4864
4865 // Zeros the given input vector
4866 virtual void
4867 zero_dst_vector_range(const unsigned int range_index) override
4868 {
4869 if (zero_dst_vector_setting)
4870 internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4871 }
4872
4873 virtual void
4874 cell_loop_pre_range(const unsigned int range_index) override
4875 {
4876 if (operation_before_loop)
4877 {
4879 matrix_free.get_dof_info(dof_handler_index_pre_post);
4880 if (range_index == numbers::invalid_unsigned_int)
4881 {
4882 // Case with threaded loop -> currently no overlap implemented
4884 0U,
4885 dof_info.vector_partitioner->locally_owned_size(),
4886 operation_before_loop,
4888 }
4889 else
4890 {
4891 AssertIndexRange(range_index,
4892 dof_info.cell_loop_pre_list_index.size() - 1);
4893 for (unsigned int id =
4894 dof_info.cell_loop_pre_list_index[range_index];
4895 id != dof_info.cell_loop_pre_list_index[range_index + 1];
4896 ++id)
4897 operation_before_loop(dof_info.cell_loop_pre_list[id].first,
4898 dof_info.cell_loop_pre_list[id].second);
4899 }
4900 }
4901 }
4902
4903 virtual void
4904 cell_loop_post_range(const unsigned int range_index) override
4905 {
4906 if (operation_after_loop)
4907 {
4908 // Run unit matrix operation on constrained dofs if we are at the
4909 // last range
4910 const std::vector<unsigned int> &partition_row_index =
4911 matrix_free.get_task_info().partition_row_index;
4912 if (range_index ==
4913 partition_row_index[partition_row_index.size() - 2] - 1)
4914 apply_operation_to_constrained_dofs(
4915 matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
4916 src,
4917 dst);
4918
4920 matrix_free.get_dof_info(dof_handler_index_pre_post);
4921 if (range_index == numbers::invalid_unsigned_int)
4922 {
4923 // Case with threaded loop -> currently no overlap implemented
4925 0U,
4926 dof_info.vector_partitioner->locally_owned_size(),
4927 operation_after_loop,
4929 }
4930 else
4931 {
4932 AssertIndexRange(range_index,
4933 dof_info.cell_loop_post_list_index.size() - 1);
4934 for (unsigned int id =
4935 dof_info.cell_loop_post_list_index[range_index];
4936 id != dof_info.cell_loop_post_list_index[range_index + 1];
4937 ++id)
4938 operation_after_loop(dof_info.cell_loop_post_list[id].first,
4939 dof_info.cell_loop_post_list[id].second);
4940 }
4941 }
4942 }
4943
4944 private:
4945 const MF &matrix_free;
4946 Container &container;
4947 function_type cell_function;
4948 function_type face_function;
4949 function_type boundary_function;
4950
4951 const InVector &src;
4952 OutVector &dst;
4953 VectorDataExchange<MF::dimension,
4954 typename MF::value_type,
4955 typename MF::vectorized_value_type>
4956 src_data_exchanger;
4957 VectorDataExchange<MF::dimension,
4958 typename MF::value_type,
4959 typename MF::vectorized_value_type>
4960 dst_data_exchanger;
4961 const bool src_and_dst_are_same;
4962 const bool zero_dst_vector_setting;
4963 const std::function<void(const unsigned int, const unsigned int)>
4964 operation_before_loop;
4965 const std::function<void(const unsigned int, const unsigned int)>
4966 operation_after_loop;
4967 const unsigned int dof_handler_index_pre_post;
4968 };
4969
4970
4971
4976 template <class MF, typename InVector, typename OutVector>
4977 struct MFClassWrapper
4978 {
4979 using function_type =
4980 std::function<void(const MF &,
4981 OutVector &,
4982 const InVector &,
4983 const std::pair<unsigned int, unsigned int> &)>;
4984
4985 MFClassWrapper(const function_type cell,
4986 const function_type face,
4987 const function_type boundary)
4988 : cell(cell)
4989 , face(face)
4990 , boundary(boundary)
4991 {}
4992
4993 void
4994 cell_integrator(const MF &mf,
4995 OutVector &dst,
4996 const InVector &src,
4997 const std::pair<unsigned int, unsigned int> &range) const
4998 {
4999 if (cell)
5000 cell(mf, dst, src, range);
5001 }
5002
5003 void
5004 face_integrator(const MF &mf,
5005 OutVector &dst,
5006 const InVector &src,
5007 const std::pair<unsigned int, unsigned int> &range) const
5008 {
5009 if (face)
5010 face(mf, dst, src, range);
5011 }
5012
5013 void
5014 boundary_integrator(
5015 const MF &mf,
5016 OutVector &dst,
5017 const InVector &src,
5018 const std::pair<unsigned int, unsigned int> &range) const
5019 {
5020 if (boundary)
5021 boundary(mf, dst, src, range);
5022 }
5023
5024 const function_type cell;
5025 const function_type face;
5026 const function_type boundary;
5027 };
5028
5029} // end of namespace internal
5030
5031
5032
5033template <int dim, typename Number, typename VectorizedArrayType>
5034template <typename OutVector, typename InVector>
5035inline void
5037 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5038 OutVector &,
5039 const InVector &,
5040 const std::pair<unsigned int, unsigned int> &)>
5041 &cell_operation,
5042 OutVector &dst,
5043 const InVector &src,
5044 const bool zero_dst_vector) const
5045{
5046 using Wrapper =
5047 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5048 InVector,
5049 OutVector>;
5050 Wrapper wrap(cell_operation, nullptr, nullptr);
5051 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5052 InVector,
5053 OutVector,
5054 Wrapper,
5055 true>
5056 worker(*this,
5057 src,
5058 dst,
5059 zero_dst_vector,
5060 wrap,
5061 &Wrapper::cell_integrator,
5062 &Wrapper::face_integrator,
5063 &Wrapper::boundary_integrator);
5064
5065 task_info.loop(worker);
5066}
5067
5068
5069
5070template <int dim, typename Number, typename VectorizedArrayType>
5071template <typename OutVector, typename InVector>
5072inline void
5074 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5075 OutVector &,
5076 const InVector &,
5077 const std::pair<unsigned int, unsigned int> &)>
5078 &cell_operation,
5079 OutVector &dst,
5080 const InVector &src,
5081 const std::function<void(const unsigned int, const unsigned int)>
5082 &operation_before_loop,
5083 const std::function<void(const unsigned int, const unsigned int)>
5084 &operation_after_loop,
5085 const unsigned int dof_handler_index_pre_post) const
5086{
5087 using Wrapper =
5088 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5089 InVector,
5090 OutVector>;
5091 Wrapper wrap(cell_operation, nullptr, nullptr);
5092 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5093 InVector,
5094 OutVector,
5095 Wrapper,
5096 true>
5097 worker(*this,
5098 src,
5099 dst,
5100 false,
5101 wrap,
5102 &Wrapper::cell_integrator,
5103 &Wrapper::face_integrator,
5104 &Wrapper::boundary_integrator,
5105 DataAccessOnFaces::none,
5106 DataAccessOnFaces::none,
5107 operation_before_loop,
5108 operation_after_loop,
5109 dof_handler_index_pre_post);
5110
5111 task_info.loop(worker);
5112}
5113
5114
5115
5116template <int dim, typename Number, typename VectorizedArrayType>
5117template <typename OutVector, typename InVector>
5118inline void
5120 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5121 OutVector &,
5122 const InVector &,
5123 const std::pair<unsigned int, unsigned int> &)>
5124 &cell_operation,
5125 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5126 OutVector &,
5127 const InVector &,
5128 const std::pair<unsigned int, unsigned int> &)>
5129 &inner_face_operation,
5130 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5131 OutVector &,
5132 const InVector &,
5133 const std::pair<unsigned int, unsigned int> &)>
5134 &boundary_face_operation,
5135 OutVector &dst,
5136 const InVector &src,
5137 const bool zero_dst_vector,
5138 const DataAccessOnFaces dst_vector_face_access,
5139 const DataAccessOnFaces src_vector_face_access) const
5140{
5141 using Wrapper =
5142 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5143 InVector,
5144 OutVector>;
5145 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5146 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5147 InVector,
5148 OutVector,
5149 Wrapper,
5150 true>
5151 worker(*this,
5152 src,
5153 dst,
5154 zero_dst_vector,
5155 wrap,
5156 &Wrapper::cell_integrator,
5157 &Wrapper::face_integrator,
5158 &Wrapper::boundary_integrator,
5159 src_vector_face_access,
5160 dst_vector_face_access);
5161
5162 task_info.loop(worker);
5163}
5164
5165
5166
5167template <int dim, typename Number, typename VectorizedArrayType>
5168template <typename CLASS, typename OutVector, typename InVector>
5169inline void
5171 void (CLASS::*function_pointer)(
5172 const MatrixFree<dim, Number, VectorizedArrayType> &,
5173 OutVector &,
5174 const InVector &,
5175 const std::pair<unsigned int, unsigned int> &) const,
5176 const CLASS *owning_class,
5177 OutVector &dst,
5178 const InVector &src,
5179 const bool zero_dst_vector) const
5180{
5181 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5182 InVector,
5183 OutVector,
5184 CLASS,
5185 true>
5186 worker(*this,
5187 src,
5188 dst,
5189 zero_dst_vector,
5190 *owning_class,
5191 function_pointer,
5192 nullptr,
5193 nullptr);
5194 task_info.loop(worker);
5195}
5196
5197
5198
5199template <int dim, typename Number, typename VectorizedArrayType>
5200template <typename CLASS, typename OutVector, typename InVector>
5201inline void
5203 void (CLASS::*function_pointer)(
5204 const MatrixFree<dim, Number, VectorizedArrayType> &,
5205 OutVector &,
5206 const InVector &,
5207 const std::pair<unsigned int, unsigned int> &) const,
5208 const CLASS *owning_class,
5209 OutVector &dst,
5210 const InVector &src,
5211 const std::function<void(const unsigned int, const unsigned int)>
5212 &operation_before_loop,
5213 const std::function<void(const unsigned int, const unsigned int)>
5214 &operation_after_loop,
5215 const unsigned int dof_handler_index_pre_post) const
5216{
5217 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5218 InVector,
5219 OutVector,
5220 CLASS,
5221 true>
5222 worker(*this,
5223 src,
5224 dst,
5225 false,
5226 *owning_class,
5227 function_pointer,
5228 nullptr,
5229 nullptr,
5232 operation_before_loop,
5233 operation_after_loop,
5234 dof_handler_index_pre_post);
5235 task_info.loop(worker);
5236}
5237
5238
5239
5240template <int dim, typename Number, typename VectorizedArrayType>
5241template <typename CLASS, typename OutVector, typename InVector>
5242inline void
5244 void (CLASS::*cell_operation)(
5245 const MatrixFree<dim, Number, VectorizedArrayType> &,
5246 OutVector &,
5247 const InVector &,
5248 const std::pair<unsigned int, unsigned int> &) const,
5249 void (CLASS::*inner_face_operation)(
5250 const MatrixFree<dim, Number, VectorizedArrayType> &,
5251 OutVector &,
5252 const InVector &,
5253 const std::pair<unsigned int, unsigned int> &) const,
5254 void (CLASS::*boundary_face_operation)(
5255 const MatrixFree<dim, Number, VectorizedArrayType> &,
5256 OutVector &,
5257 const InVector &,
5258 const std::pair<unsigned int, unsigned int> &) const,
5259 const CLASS *owning_class,
5260 OutVector &dst,
5261 const InVector &src,
5262 const bool zero_dst_vector,
5263 const DataAccessOnFaces dst_vector_face_access,
5264 const DataAccessOnFaces src_vector_face_access) const
5265{
5266 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5267 InVector,
5268 OutVector,
5269 CLASS,
5270 true>
5271 worker(*this,
5272 src,
5273 dst,
5274 zero_dst_vector,
5275 *owning_class,
5276 cell_operation,
5277 inner_face_operation,
5278 boundary_face_operation,
5279 src_vector_face_access,
5280 dst_vector_face_access);
5281 task_info.loop(worker);
5282}
5283
5284
5285
5286template <int dim, typename Number, typename VectorizedArrayType>
5287template <typename CLASS, typename OutVector, typename InVector>
5288inline void
5290 void (CLASS::*function_pointer)(
5291 const MatrixFree<dim, Number, VectorizedArrayType> &,
5292 OutVector &,
5293 const InVector &,
5294 const std::pair<unsigned int, unsigned int> &),
5295 CLASS *owning_class,
5296 OutVector &dst,
5297 const InVector &src,
5298 const bool zero_dst_vector) const
5299{
5300 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5301 InVector,
5302 OutVector,
5303 CLASS,
5304 false>
5305 worker(*this,
5306 src,
5307 dst,
5308 zero_dst_vector,
5309 *owning_class,
5310 function_pointer,
5311 nullptr,
5312 nullptr);
5313 task_info.loop(worker);
5314}
5315
5316
5317
5318template <int dim, typename Number, typename VectorizedArrayType>
5319template <typename CLASS, typename OutVector, typename InVector>
5320inline void
5322 void (CLASS::*function_pointer)(
5323 const MatrixFree<dim, Number, VectorizedArrayType> &,
5324 OutVector &,
5325 const InVector &,
5326 const std::pair<unsigned int, unsigned int> &),
5327 CLASS *owning_class,
5328 OutVector &dst,
5329 const InVector &src,
5330 const std::function<void(const unsigned int, const unsigned int)>
5331 &operation_before_loop,
5332 const std::function<void(const unsigned int, const unsigned int)>
5333 &operation_after_loop,
5334 const unsigned int dof_handler_index_pre_post) const
5335{
5336 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5337 InVector,
5338 OutVector,
5339 CLASS,
5340 false>
5341 worker(*this,
5342 src,
5343 dst,
5344 false,
5345 *owning_class,
5346 function_pointer,
5347 nullptr,
5348 nullptr,
5351 operation_before_loop,
5352 operation_after_loop,
5353 dof_handler_index_pre_post);
5354 task_info.loop(worker);
5355}
5356
5357
5358
5359template <int dim, typename Number, typename VectorizedArrayType>
5360template <typename CLASS, typename OutVector, typename InVector>
5361inline void
5363 void (CLASS::*cell_operation)(
5364 const MatrixFree<dim, Number, VectorizedArrayType> &,
5365 OutVector &,
5366 const InVector &,
5367 const std::pair<unsigned int, unsigned int> &),
5368 void (CLASS::*inner_face_operation)(
5369 const MatrixFree<dim, Number, VectorizedArrayType> &,
5370 OutVector &,
5371 const InVector &,
5372 const std::pair<unsigned int, unsigned int> &),
5373 void (CLASS::*boundary_face_operation)(
5374 const MatrixFree<dim, Number, VectorizedArrayType> &,
5375 OutVector &,
5376 const InVector &,
5377 const std::pair<unsigned int, unsigned int> &),
5378 CLASS *owning_class,
5379 OutVector &dst,
5380 const InVector &src,
5381 const bool zero_dst_vector,
5382 const DataAccessOnFaces dst_vector_face_access,
5383 const DataAccessOnFaces src_vector_face_access) const
5384{
5385 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5386 InVector,
5387 OutVector,
5388 CLASS,
5389 false>
5390 worker(*this,
5391 src,
5392 dst,
5393 zero_dst_vector,
5394 *owning_class,
5395 cell_operation,
5396 inner_face_operation,
5397 boundary_face_operation,
5398 src_vector_face_access,
5399 dst_vector_face_access);
5400 task_info.loop(worker);
5401}
5402
5403
5404
5405template <int dim, typename Number, typename VectorizedArrayType>
5406template <typename OutVector, typename InVector>
5407inline void
5409 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5410 OutVector &,
5411 const InVector &,
5412 const std::pair<unsigned int, unsigned int> &)>
5413 &cell_operation,
5414 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5415 OutVector &,
5416 const InVector &,
5417 const std::pair<unsigned int, unsigned int> &)>
5418 &inner_face_operation,
5419 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5420 OutVector &,
5421 const InVector &,
5422 const std::pair<unsigned int, unsigned int> &)>
5423 &boundary_face_operation,
5424 OutVector &dst,
5425 const InVector &src,
5426 const std::function<void(const unsigned int, const unsigned int)>
5427 &operation_before_loop,
5428 const std::function<void(const unsigned int, const unsigned int)>
5429 &operation_after_loop,
5430 const unsigned int dof_handler_index_pre_post,
5431 const DataAccessOnFaces dst_vector_face_access,
5432 const DataAccessOnFaces src_vector_face_access) const
5433{
5434 using Wrapper =
5435 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5436 InVector,
5437 OutVector>;
5438 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5439 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5440 InVector,
5441 OutVector,
5442 Wrapper,
5443 true>
5444 worker(*this,
5445 src,
5446 dst,
5447 false,
5448 wrap,
5449 &Wrapper::cell_integrator,
5450 &Wrapper::face_integrator,
5451 &Wrapper::boundary_integrator,
5452 src_vector_face_access,
5453 dst_vector_face_access,
5454 operation_before_loop,
5455 operation_after_loop,
5456 dof_handler_index_pre_post);
5457
5458 task_info.loop(worker);
5459}
5460
5461
5462
5463template <int dim, typename Number, typename VectorizedArrayType>
5464template <typename CLASS, typename OutVector, typename InVector>
5465inline void
5467 void (CLASS::*cell_operation)(const MatrixFree &,
5468 OutVector &,
5469 const InVector &,
5470 const std::pair<unsigned int, unsigned int> &)
5471 const,
5472 void (CLASS::*inner_face_operation)(
5473 const MatrixFree &,
5474 OutVector &,
5475 const InVector &,
5476 const std::pair<unsigned int, unsigned int> &) const,
5477 void (CLASS::*boundary_face_operation)(
5478 const MatrixFree &,
5479 OutVector &,
5480 const InVector &,
5481 const std::pair<unsigned int, unsigned int> &) const,
5482 const CLASS *owning_class,
5483 OutVector &dst,
5484 const InVector &src,
5485 const std::function<void(const unsigned int, const unsigned int)>
5486 &operation_before_loop,
5487 const std::function<void(const unsigned int, const unsigned int)>
5488 &operation_after_loop,
5489 const unsigned int dof_handler_index_pre_post,
5490 const DataAccessOnFaces dst_vector_face_access,
5491 const DataAccessOnFaces src_vector_face_access) const
5492{
5493 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5494 InVector,
5495 OutVector,
5496 CLASS,
5497 true>
5498 worker(*this,
5499 src,
5500 dst,
5501 false,
5502 *owning_class,
5503 cell_operation,
5504 inner_face_operation,
5505 boundary_face_operation,
5506 src_vector_face_access,
5507 dst_vector_face_access,
5508 operation_before_loop,
5509 operation_after_loop,
5510 dof_handler_index_pre_post);
5511 task_info.loop(worker);
5512}
5513
5514
5515
5516template <int dim, typename Number, typename VectorizedArrayType>
5517template <typename CLASS, typename OutVector, typename InVector>
5518inline void
5520 void (CLASS::*cell_operation)(const MatrixFree &,
5521 OutVector &,
5522 const InVector &,
5523 const std::pair<unsigned int, unsigned int> &),
5524 void (CLASS::*inner_face_operation)(
5525 const MatrixFree &,
5526 OutVector &,
5527 const InVector &,
5528 const std::pair<unsigned int, unsigned int> &),
5529 void (CLASS::*boundary_face_operation)(
5530 const MatrixFree &,
5531 OutVector &,
5532 const InVector &,
5533 const std::pair<unsigned int, unsigned int> &),
5534 const CLASS *owning_class,
5535 OutVector &dst,
5536 const InVector &src,
5537 const std::function<void(const unsigned int, const unsigned int)>
5538 &operation_before_loop,
5539 const std::function<void(const unsigned int, const unsigned int)>
5540 &operation_after_loop,
5541 const unsigned int dof_handler_index_pre_post,
5542 const DataAccessOnFaces dst_vector_face_access,
5543 const DataAccessOnFaces src_vector_face_access) const
5544{
5545 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5546 InVector,
5547 OutVector,
5548 CLASS,
5549 false>
5550 worker(*this,
5551 src,
5552 dst,
5553 false,
5554 *owning_class,
5555 cell_operation,
5556 inner_face_operation,
5557 boundary_face_operation,
5558 src_vector_face_access,
5559 dst_vector_face_access,
5560 operation_before_loop,
5561 operation_after_loop,
5562 dof_handler_index_pre_post);
5563 task_info.loop(worker);
5564}
5565
5566
5567
5568template <int dim, typename Number, typename VectorizedArrayType>
5569template <typename CLASS, typename OutVector, typename InVector>
5570inline void
5572 void (CLASS::*function_pointer)(
5573 const MatrixFree<dim, Number, VectorizedArrayType> &,
5574 OutVector &,
5575 const InVector &,
5576 const std::pair<unsigned int, unsigned int> &) const,
5577 const CLASS *owning_class,
5578 OutVector &dst,
5579 const InVector &src,
5580 const bool zero_dst_vector,
5581 const DataAccessOnFaces src_vector_face_access) const
5582{
5583 auto src_vector_face_access_temp = src_vector_face_access;
5584 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5585 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5586 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5587 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5588
5589 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5590 InVector,
5591 OutVector,
5592 CLASS,
5593 true>
5594 worker(*this,
5595 src,
5596 dst,
5597 zero_dst_vector,
5598 *owning_class,
5599 function_pointer,
5600 nullptr,
5601 nullptr,
5602 src_vector_face_access_temp,
5604 task_info.loop(worker);
5605}
5606
5607
5608
5609template <int dim, typename Number, typename VectorizedArrayType>
5610template <typename CLASS, typename OutVector, typename InVector>
5611inline void
5613 void (CLASS::*function_pointer)(
5614 const MatrixFree<dim, Number, VectorizedArrayType> &,
5615 OutVector &,
5616 const InVector &,
5617 const std::pair<unsigned int, unsigned int> &),
5618 CLASS *owning_class,
5619 OutVector &dst,
5620 const InVector &src,
5621 const bool zero_dst_vector,
5622 const DataAccessOnFaces src_vector_face_access) const
5623{
5624 auto src_vector_face_access_temp = src_vector_face_access;
5625 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5626 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5627 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5628 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5629
5630 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5631 InVector,
5632 OutVector,
5633 CLASS,
5634 false>
5635 worker(*this,
5636 src,
5637 dst,
5638 zero_dst_vector,
5639 *owning_class,
5640 function_pointer,
5641 nullptr,
5642 nullptr,
5643 src_vector_face_access_temp,
5645 task_info.loop(worker);
5646}
5647
5648
5649
5650template <int dim, typename Number, typename VectorizedArrayType>
5651template <typename OutVector, typename InVector>
5652inline void
5654 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5655 OutVector &,
5656 const InVector &,
5657 const std::pair<unsigned int, unsigned int> &)>
5658 &cell_operation,
5659 OutVector &dst,
5660 const InVector &src,
5661 const bool zero_dst_vector,
5662 const DataAccessOnFaces src_vector_face_access) const
5663{
5664 auto src_vector_face_access_temp = src_vector_face_access;
5665 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5666 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5667 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5668 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5669
5670 using Wrapper =
5671 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5672 InVector,
5673 OutVector>;
5674 Wrapper wrap(cell_operation, nullptr, nullptr);
5675
5676 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5677 InVector,
5678 OutVector,
5679 Wrapper,
5680 true>
5681 worker(*this,
5682 src,
5683 dst,
5684 zero_dst_vector,
5685 wrap,
5686 &Wrapper::cell_integrator,
5687 &Wrapper::face_integrator,
5688 &Wrapper::boundary_integrator,
5689 src_vector_face_access_temp,
5690 DataAccessOnFaces::none);
5691 task_info.loop(worker);
5692}
5693
5694
5695#endif // ifndef DOXYGEN
5696
5697
5698
5700
5701#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:320
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
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
unsigned int get_cell_category(const unsigned int cell_batch_index, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) 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
std::pair< unsigned int, unsigned int > get_face_range_category(const std::pair< unsigned int, unsigned int > face_batch_range, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) 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
std::vector< ObserverPointer< const DoFHandler< dim > > > dof_handlers
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
unsigned int first_hp_dof_handler_index
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
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
std::vector< ObserverPointer< const AffineConstraints< Number > > > affine_constraints
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_face_active_fe_index(const std::pair< unsigned int, unsigned int > range, const bool is_interior_face=true, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) const
unsigned int get_cell_active_fe_index(const std::pair< unsigned int, unsigned int > range, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) 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 dof_handler_index=numbers::invalid_unsigned_int) 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
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)
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
unsigned int get_cell_range_category(const std::pair< unsigned int, unsigned int > cell_batch_range, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) 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 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())
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:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
Point< 2 > second
Definition grid_out.cc:4633
Point< 2 > first
Definition grid_out.cc:4632
unsigned int level
Definition grid_out.cc:4635
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.
std::vector< index_type > data
Definition mpi.cc:746
std::size_t size
Definition mpi.cc:745
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
VectorType::value_type * begin(VectorType &V)
bool job_supports_mpi()
Definition mpi.cc:692
unsigned int minimum_parallel_grain_size
Definition parallel.cc:50
constexpr unsigned int invalid_unsigned_int
Definition types.h:238
constexpr types::boundary_id invalid_boundary_id
Definition types.h:309
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:540
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
unsigned short int fe_index
Definition types.h:72
unsigned int boundary_id
Definition types.h:161
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
UpdateFlags mapping_update_flags_boundary_faces
UpdateFlags mapping_update_flags_faces_by_cells
AdditionalData & operator=(const AdditionalData &other)=default
AdditionalData(const AdditionalData &other)
std::vector< unsigned int > cell_loop_pre_list_index
Definition dof_info.h:749
std::vector< unsigned int > cell_loop_post_list_index
Definition dof_info.h:762
std::vector< std::pair< unsigned int, unsigned int > > vector_zero_range_list
Definition dof_info.h:742
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition dof_info.h:593
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_pre_list
Definition dof_info.h:755
std::vector< unsigned int > vector_zero_range_list_index
Definition dof_info.h:737
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_post_list
Definition dof_info.h:768
void loop(MFWorkerInterface &worker) const
Definition task_info.cc:350