Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-3063-g10967a192b 2025-04-12 13:30:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
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
27
29
30#include <deal.II/fe/fe.h>
31#include <deal.II/fe/mapping.h>
32
35
40
47
48#include <cstdlib>
49#include <limits>
50#include <list>
51#include <memory>
52
53
55
56
57
109template <int dim,
110 typename Number = double,
111 typename VectorizedArrayType = VectorizedArray<Number>>
113{
114 static_assert(
115 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
116 "Type of Number and of VectorizedArrayType do not match.");
117
118public:
123 using value_type = Number;
124 using vectorized_value_type = VectorizedArrayType;
125
129 static constexpr unsigned int dimension = dim;
130
184 {
189
217
223 const unsigned int tasks_block_size = 0,
229 const unsigned int mg_level = numbers::invalid_unsigned_int,
230 const bool store_plain_indices = true,
231 const bool initialize_indices = true,
232 const bool initialize_mapping = true,
233 const bool overlap_communication_computation = true,
234 const bool hold_all_faces_to_owned_cells = false,
235 const bool cell_vectorization_categories_strict = false,
236 const bool allow_ghosted_vectors_in_loops = true)
252 , store_ghost_cells(false)
253 , communicator_sm(MPI_COMM_SELF)
254 {}
255
282
287 operator=(const AdditionalData &other) = default;
288
326
336 unsigned int tasks_block_size;
337
352
372
392
420
429 unsigned int mg_level;
430
438
449
458
472
481
503 std::vector<unsigned int> cell_vectorization_category;
504
515
527
534
539 };
540
549
554
558 ~MatrixFree() override = default;
559
572 template <typename QuadratureType, typename number2, typename MappingType>
573 void
574 reinit(const MappingType &mapping,
575 const DoFHandler<dim> &dof_handler,
576 const AffineConstraints<number2> &constraint,
577 const QuadratureType &quad,
578 const AdditionalData &additional_data = AdditionalData());
579
601 template <typename QuadratureType, typename number2, typename MappingType>
602 void
603 reinit(const MappingType &mapping,
604 const std::vector<const DoFHandler<dim> *> &dof_handler,
605 const std::vector<const AffineConstraints<number2> *> &constraint,
606 const std::vector<QuadratureType> &quad,
607 const AdditionalData &additional_data = AdditionalData());
608
616 template <typename QuadratureType, typename number2, typename MappingType>
617 void
618 reinit(const MappingType &mapping,
619 const std::vector<const DoFHandler<dim> *> &dof_handler,
620 const std::vector<const AffineConstraints<number2> *> &constraint,
621 const QuadratureType &quad,
622 const AdditionalData &additional_data = AdditionalData());
623
629 void
631 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free_base);
632
642 void
644
648 void
649 update_mapping(const std::shared_ptr<hp::MappingCollection<dim>> &mapping);
650
655 void
657
672 {
679 none,
680
691 values,
692
702
713 gradients,
714
724
732 };
733
778 template <typename OutVector, typename InVector>
779 void
780 cell_loop(const std::function<void(
782 OutVector &,
783 const InVector &,
784 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
785 OutVector &dst,
786 const InVector &src,
787 const bool zero_dst_vector = false) const;
788
835 template <typename CLASS, typename OutVector, typename InVector>
836 void
837 cell_loop(void (CLASS::*cell_operation)(
838 const MatrixFree &,
839 OutVector &,
840 const InVector &,
841 const std::pair<unsigned int, unsigned int> &) const,
842 const CLASS *owning_class,
843 OutVector &dst,
844 const InVector &src,
845 const bool zero_dst_vector = false) const;
846
850 template <typename CLASS, typename OutVector, typename InVector>
851 void
852 cell_loop(void (CLASS::*cell_operation)(
853 const MatrixFree &,
854 OutVector &,
855 const InVector &,
856 const std::pair<unsigned int, unsigned int> &),
857 CLASS *owning_class,
858 OutVector &dst,
859 const InVector &src,
860 const bool zero_dst_vector = false) const;
861
946 template <typename CLASS, typename OutVector, typename InVector>
947 void
948 cell_loop(void (CLASS::*cell_operation)(
949 const MatrixFree &,
950 OutVector &,
951 const InVector &,
952 const std::pair<unsigned int, unsigned int> &) const,
953 const CLASS *owning_class,
954 OutVector &dst,
955 const InVector &src,
956 const std::function<void(const unsigned int, const unsigned int)>
957 &operation_before_loop,
958 const std::function<void(const unsigned int, const unsigned int)>
959 &operation_after_loop,
960 const unsigned int dof_handler_index_pre_post = 0) const;
961
965 template <typename CLASS, typename OutVector, typename InVector>
966 void
967 cell_loop(void (CLASS::*cell_operation)(
968 const MatrixFree &,
969 OutVector &,
970 const InVector &,
971 const std::pair<unsigned int, unsigned int> &),
972 CLASS *owning_class,
973 OutVector &dst,
974 const InVector &src,
975 const std::function<void(const unsigned int, const unsigned int)>
976 &operation_before_loop,
977 const std::function<void(const unsigned int, const unsigned int)>
978 &operation_after_loop,
979 const unsigned int dof_handler_index_pre_post = 0) const;
980
985 template <typename OutVector, typename InVector>
986 void
987 cell_loop(const std::function<void(
989 OutVector &,
990 const InVector &,
991 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
992 OutVector &dst,
993 const InVector &src,
994 const std::function<void(const unsigned int, const unsigned int)>
995 &operation_before_loop,
996 const std::function<void(const unsigned int, const unsigned int)>
997 &operation_after_loop,
998 const unsigned int dof_handler_index_pre_post = 0) const;
999
1075 template <typename OutVector, typename InVector>
1076 void
1078 const std::function<
1080 OutVector &,
1081 const InVector &,
1082 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1083 const std::function<void(
1085 OutVector &,
1086 const InVector &,
1087 const std::pair<unsigned int, unsigned int> &)> &inner_face_operation,
1088 const std::function<void(
1090 OutVector &,
1091 const InVector &,
1092 const std::pair<unsigned int, unsigned int> &)> &boundary_face_operation,
1093 OutVector &dst,
1094 const InVector &src,
1095 const bool zero_dst_vector = false,
1096 const DataAccessOnFaces dst_vector_face_access =
1098 const DataAccessOnFaces src_vector_face_access =
1100
1189 template <typename CLASS, typename OutVector, typename InVector>
1190 void
1191 loop(void (CLASS::*cell_operation)(
1192 const MatrixFree &,
1193 OutVector &,
1194 const InVector &,
1195 const std::pair<unsigned int, unsigned int> &) const,
1196 void (CLASS::*inner_face_operation)(
1197 const MatrixFree &,
1198 OutVector &,
1199 const InVector &,
1200 const std::pair<unsigned int, unsigned int> &) const,
1201 void (CLASS::*boundary_face_operation)(
1202 const MatrixFree &,
1203 OutVector &,
1204 const InVector &,
1205 const std::pair<unsigned int, unsigned int> &) const,
1206 const CLASS *owning_class,
1207 OutVector &dst,
1208 const InVector &src,
1209 const bool zero_dst_vector = false,
1210 const DataAccessOnFaces dst_vector_face_access =
1212 const DataAccessOnFaces src_vector_face_access =
1214
1218 template <typename CLASS, typename OutVector, typename InVector>
1219 void
1220 loop(void (CLASS::*cell_operation)(
1221 const MatrixFree &,
1222 OutVector &,
1223 const InVector &,
1224 const std::pair<unsigned int, unsigned int> &),
1225 void (CLASS::*inner_face_operation)(
1226 const MatrixFree &,
1227 OutVector &,
1228 const InVector &,
1229 const std::pair<unsigned int, unsigned int> &),
1230 void (CLASS::*boundary_face_operation)(
1231 const MatrixFree &,
1232 OutVector &,
1233 const InVector &,
1234 const std::pair<unsigned int, unsigned int> &),
1235 CLASS *owning_class,
1236 OutVector &dst,
1237 const InVector &src,
1238 const bool zero_dst_vector = false,
1239 const DataAccessOnFaces dst_vector_face_access =
1241 const DataAccessOnFaces src_vector_face_access =
1243
1365 template <typename CLASS, typename OutVector, typename InVector>
1366 void
1367 loop(void (CLASS::*cell_operation)(
1368 const MatrixFree &,
1369 OutVector &,
1370 const InVector &,
1371 const std::pair<unsigned int, unsigned int> &) const,
1372 void (CLASS::*inner_face_operation)(
1373 const MatrixFree &,
1374 OutVector &,
1375 const InVector &,
1376 const std::pair<unsigned int, unsigned int> &) const,
1377 void (CLASS::*boundary_face_operation)(
1378 const MatrixFree &,
1379 OutVector &,
1380 const InVector &,
1381 const std::pair<unsigned int, unsigned int> &) const,
1382 const CLASS *owning_class,
1383 OutVector &dst,
1384 const InVector &src,
1385 const std::function<void(const unsigned int, const unsigned int)>
1386 &operation_before_loop,
1387 const std::function<void(const unsigned int, const unsigned int)>
1388 &operation_after_loop,
1389 const unsigned int dof_handler_index_pre_post = 0,
1390 const DataAccessOnFaces dst_vector_face_access =
1392 const DataAccessOnFaces src_vector_face_access =
1394
1398 template <typename CLASS, typename OutVector, typename InVector>
1399 void
1400 loop(void (CLASS::*cell_operation)(
1401 const MatrixFree &,
1402 OutVector &,
1403 const InVector &,
1404 const std::pair<unsigned int, unsigned int> &),
1405 void (CLASS::*inner_face_operation)(
1406 const MatrixFree &,
1407 OutVector &,
1408 const InVector &,
1409 const std::pair<unsigned int, unsigned int> &),
1410 void (CLASS::*boundary_face_operation)(
1411 const MatrixFree &,
1412 OutVector &,
1413 const InVector &,
1414 const std::pair<unsigned int, unsigned int> &),
1415 const CLASS *owning_class,
1416 OutVector &dst,
1417 const InVector &src,
1418 const std::function<void(const unsigned int, const unsigned int)>
1419 &operation_before_loop,
1420 const std::function<void(const unsigned int, const unsigned int)>
1421 &operation_after_loop,
1422 const unsigned int dof_handler_index_pre_post = 0,
1423 const DataAccessOnFaces dst_vector_face_access =
1425 const DataAccessOnFaces src_vector_face_access =
1427
1433 template <typename OutVector, typename InVector>
1434 void
1436 const std::function<
1438 OutVector &,
1439 const InVector &,
1440 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1441 const std::function<void(
1443 OutVector &,
1444 const InVector &,
1445 const std::pair<unsigned int, unsigned int> &)> &inner_face_operation,
1446 const std::function<void(
1448 OutVector &,
1449 const InVector &,
1450 const std::pair<unsigned int, unsigned int> &)> &boundary_face_operation,
1451 OutVector &dst,
1452 const InVector &src,
1453 const std::function<void(const unsigned int, const unsigned int)>
1454 &operation_before_loop,
1455 const std::function<void(const unsigned int, const unsigned int)>
1456 &operation_after_loop,
1457 const unsigned int dof_handler_index_pre_post = 0,
1458 const DataAccessOnFaces dst_vector_face_access =
1460 const DataAccessOnFaces src_vector_face_access =
1462
1527 template <typename CLASS, typename OutVector, typename InVector>
1528 void
1529 loop_cell_centric(void (CLASS::*cell_operation)(
1530 const MatrixFree &,
1531 OutVector &,
1532 const InVector &,
1533 const std::pair<unsigned int, unsigned int> &) const,
1534 const CLASS *owning_class,
1535 OutVector &dst,
1536 const InVector &src,
1537 const bool zero_dst_vector = false,
1538 const DataAccessOnFaces src_vector_face_access =
1540
1544 template <typename CLASS, typename OutVector, typename InVector>
1545 void
1546 loop_cell_centric(void (CLASS::*cell_operation)(
1547 const MatrixFree &,
1548 OutVector &,
1549 const InVector &,
1550 const std::pair<unsigned int, unsigned int> &),
1551 CLASS *owning_class,
1552 OutVector &dst,
1553 const InVector &src,
1554 const bool zero_dst_vector = false,
1555 const DataAccessOnFaces src_vector_face_access =
1557
1561 template <typename OutVector, typename InVector>
1562 void
1564 const std::function<void(const MatrixFree &,
1565 OutVector &,
1566 const InVector &,
1567 const std::pair<unsigned int, unsigned int> &)>
1568 &cell_operation,
1569 OutVector &dst,
1570 const InVector &src,
1571 const bool zero_dst_vector = false,
1572 const DataAccessOnFaces src_vector_face_access =
1574
1582 std::pair<unsigned int, unsigned int>
1583 create_cell_subrange_hp(const std::pair<unsigned int, unsigned int> &range,
1584 const unsigned int fe_degree,
1585 const unsigned int dof_handler_index = 0) const;
1586
1593 std::pair<unsigned int, unsigned int>
1595 const std::pair<unsigned int, unsigned int> &range,
1596 const unsigned int fe_index,
1597 const unsigned int dof_handler_index = 0) const;
1598
1602 unsigned int
1604
1608 unsigned int
1610 const std::pair<unsigned int, unsigned int> range) const;
1611
1615 unsigned int
1616 get_face_active_fe_index(const std::pair<unsigned int, unsigned int> range,
1617 const bool is_interior_face = true) const;
1618
1630 template <typename T>
1631 void
1633
1639 template <typename T>
1640 void
1642
1656 template <typename VectorType>
1657 void
1658 initialize_dof_vector(VectorType &vec,
1659 const unsigned int dof_handler_index = 0) const;
1660
1678 template <typename Number2>
1679 void
1681 const unsigned int dof_handler_index = 0) const;
1682
1693 const std::shared_ptr<const Utilities::MPI::Partitioner> &
1694 get_vector_partitioner(const unsigned int dof_handler_index = 0) const;
1695
1699 const IndexSet &
1700 get_locally_owned_set(const unsigned int dof_handler_index = 0) const;
1701
1705 const IndexSet &
1706 get_ghost_set(const unsigned int dof_handler_index = 0) const;
1707
1717 const std::vector<unsigned int> &
1718 get_constrained_dofs(const unsigned int dof_handler_index = 0) const;
1719
1730 void
1731 renumber_dofs(std::vector<types::global_dof_index> &renumbering,
1732 const unsigned int dof_handler_index = 0);
1733
1743 template <int spacedim>
1744 static bool
1746
1750 unsigned int
1752
1757 unsigned int
1758 n_base_elements(const unsigned int dof_handler_index) const;
1759
1767 unsigned int
1769
1779 unsigned int
1781
1788 unsigned int
1790
1800 unsigned int
1802
1813 unsigned int
1815
1820 unsigned int
1822
1834 get_boundary_id(const unsigned int face_batch_index) const;
1835
1840 std::array<types::boundary_id, VectorizedArrayType::size()>
1841 get_faces_by_cells_boundary_id(const unsigned int cell_batch_index,
1842 const unsigned int face_number) const;
1843
1848 const DoFHandler<dim> &
1849 get_dof_handler(const unsigned int dof_handler_index = 0) const;
1850
1858 get_affine_constraints(const unsigned int dof_handler_index = 0) const;
1859
1873 get_cell_iterator(const unsigned int cell_batch_index,
1874 const unsigned int lane_index,
1875 const unsigned int dof_handler_index = 0) const;
1876
1882 std::pair<int, int>
1883 get_cell_level_and_index(const unsigned int cell_batch_index,
1884 const unsigned int lane_index) const;
1885
1892 unsigned int
1894 const typename Triangulation<dim>::cell_iterator &cell) const;
1895
1908 std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>
1909 get_face_iterator(const unsigned int face_batch_index,
1910 const unsigned int lane_index,
1911 const bool interior = true,
1912 const unsigned int fe_component = 0) const;
1913
1926 bool
1927 at_irregular_cell(const unsigned int cell_batch_index) const;
1928
1938 unsigned int
1939 n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const;
1940
1950 unsigned int
1951 n_active_entries_per_face_batch(const unsigned int face_batch_index) const;
1952
1956 unsigned int
1957 get_dofs_per_cell(const unsigned int dof_handler_index = 0,
1958 const unsigned int hp_active_fe_index = 0) const;
1959
1963 unsigned int
1964 get_n_q_points(const unsigned int quad_index = 0,
1965 const unsigned int hp_active_fe_index = 0) const;
1966
1971 unsigned int
1972 get_dofs_per_face(const unsigned int dof_handler_index = 0,
1973 const unsigned int hp_active_fe_index = 0) const;
1974
1979 unsigned int
1980 get_n_q_points_face(const unsigned int quad_index = 0,
1981 const unsigned int hp_active_fe_index = 0) const;
1982
1986 const Quadrature<dim> &
1987 get_quadrature(const unsigned int quad_index = 0,
1988 const unsigned int hp_active_fe_index = 0) const;
1989
1993 const Quadrature<dim - 1> &
1994 get_face_quadrature(const unsigned int quad_index = 0,
1995 const unsigned int hp_active_fe_index = 0) const;
1996
2009 unsigned int
2011 const std::pair<unsigned int, unsigned int> cell_batch_range) const;
2012
2017 std::pair<unsigned int, unsigned int>
2019 const std::pair<unsigned int, unsigned int> face_batch_range) const;
2020
2032 unsigned int
2033 get_cell_category(const unsigned int cell_batch_index) const;
2034
2039 std::pair<unsigned int, unsigned int>
2040 get_face_category(const unsigned int face_batch_index) const;
2041
2045 bool
2047
2052 bool
2054
2059 unsigned int
2061
2066 std::size_t
2068
2073 template <typename StreamType>
2074 void
2075 print_memory_consumption(StreamType &out) const;
2076
2081 void
2082 print(std::ostream &out) const;
2083
2097
2098 /*
2099 * Return geometry-dependent information on the cells.
2100 */
2101 const internal::MatrixFreeFunctions::
2102 MappingInfo<dim, Number, VectorizedArrayType> &
2104
2109 get_dof_info(const unsigned int dof_handler_index_component = 0) const;
2110
2114 unsigned int
2116
2121 const Number *
2122 constraint_pool_begin(const unsigned int pool_index) const;
2123
2129 const Number *
2130 constraint_pool_end(const unsigned int pool_index) const;
2131
2136 get_shape_info(const unsigned int dof_handler_index_component = 0,
2137 const unsigned int quad_index = 0,
2138 const unsigned int fe_base_element = 0,
2139 const unsigned int hp_active_fe_index = 0,
2140 const unsigned int hp_active_quad_index = 0) const;
2141
2146 VectorizedArrayType::size()> &
2147 get_face_info(const unsigned int face_batch_index) const;
2148
2149
2157
2173
2177 void
2179
2191
2195 void
2197 const AlignedVector<Number> *memory) const;
2198
2201private:
2206 template <typename number2, int q_dim>
2207 void
2209 const std::shared_ptr<hp::MappingCollection<dim>> &mapping,
2210 const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
2211 const std::vector<const AffineConstraints<number2> *> &constraint,
2212 const std::vector<IndexSet> &locally_owned_set,
2213 const std::vector<hp::QCollection<q_dim>> &quad,
2214 const AdditionalData &additional_data);
2215
2222 template <typename number2>
2223 void
2225 const std::vector<const AffineConstraints<number2> *> &constraint,
2226 const std::vector<IndexSet> &locally_owned_set,
2227 const AdditionalData &additional_data);
2228
2232 void
2234 const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
2235 const AdditionalData &additional_data);
2236
2240 std::vector<ObserverPointer<const DoFHandler<dim>>> dof_handlers;
2241
2248 std::vector<ObserverPointer<const AffineConstraints<Number>>>
2250
2255 std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
2256
2263 std::vector<Number> constraint_pool_data;
2264
2269 std::vector<unsigned int> constraint_pool_row_index;
2270
2277
2282
2289 std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
2290
2295 std::vector<unsigned int> mf_cell_indices;
2296
2304
2311
2316 internal::MatrixFreeFunctions::FaceInfo<VectorizedArrayType::size()>
2318
2323
2328
2337 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2339
2344 mutable std::list<std::pair<bool, AlignedVector<Number>>>
2346
2350 unsigned int mg_level;
2351
2357};
2358
2359
2360
2361/*----------------------- Inline functions ----------------------------------*/
2362
2363#ifndef DOXYGEN
2364
2365
2366
2367template <int dim, typename Number, typename VectorizedArrayType>
2368template <typename T>
2369inline void
2371 AlignedVector<T> &vec) const
2372{
2373 vec.resize(this->n_cell_batches() + this->n_ghost_cell_batches());
2374}
2375
2376
2377
2378template <int dim, typename Number, typename VectorizedArrayType>
2379template <typename T>
2380inline void
2382 AlignedVector<T> &vec) const
2383{
2384 vec.resize(this->n_inner_face_batches() + this->n_boundary_face_batches() +
2385 this->n_ghost_inner_face_batches());
2386}
2387
2388
2389
2390template <int dim, typename Number, typename VectorizedArrayType>
2391template <typename VectorType>
2392inline void
2394 VectorType &vec,
2395 const unsigned int comp) const
2396{
2397 static_assert(IsBlockVector<VectorType>::value == false,
2398 "This function is not supported for block vectors.");
2399
2400 Assert(task_info.n_procs == 1,
2401 ExcMessage("This function can only be used in serial."));
2402
2403 AssertIndexRange(comp, n_components());
2404 vec.reinit(dof_info[comp].vector_partitioner->size());
2405}
2406
2407
2408
2409template <int dim, typename Number, typename VectorizedArrayType>
2410template <typename Number2>
2411inline void
2414 const unsigned int comp) const
2415{
2416 AssertIndexRange(comp, n_components());
2417 vec.reinit(dof_info[comp].vector_partitioner, task_info.communicator_sm);
2418}
2419
2420
2421
2422template <int dim, typename Number, typename VectorizedArrayType>
2423inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2425 const unsigned int comp) const
2426{
2427 AssertIndexRange(comp, n_components());
2428 return dof_info[comp].vector_partitioner;
2429}
2430
2431
2432
2433template <int dim, typename Number, typename VectorizedArrayType>
2434inline const std::vector<unsigned int> &
2436 const unsigned int comp) const
2437{
2438 AssertIndexRange(comp, n_components());
2439 return dof_info[comp].constrained_dofs;
2440}
2441
2442
2443
2444template <int dim, typename Number, typename VectorizedArrayType>
2445inline unsigned int
2447{
2448 AssertDimension(dof_handlers.size(), dof_info.size());
2449 return dof_handlers.size();
2450}
2451
2452
2453
2454template <int dim, typename Number, typename VectorizedArrayType>
2455inline unsigned int
2457 const unsigned int dof_no) const
2458{
2459 AssertDimension(dof_handlers.size(), dof_info.size());
2460 AssertIndexRange(dof_no, dof_handlers.size());
2461 return dof_handlers[dof_no]->get_fe().n_base_elements();
2462}
2463
2464
2465
2466template <int dim, typename Number, typename VectorizedArrayType>
2469{
2470 return task_info;
2471}
2472
2473
2474
2475template <int dim, typename Number, typename VectorizedArrayType>
2476inline unsigned int
2478{
2479 return task_info.n_active_cells;
2480}
2481
2482
2483
2484template <int dim, typename Number, typename VectorizedArrayType>
2485inline unsigned int
2487{
2488 return *(task_info.cell_partition_data.end() - 2);
2489}
2490
2491
2492
2493template <int dim, typename Number, typename VectorizedArrayType>
2494inline unsigned int
2496{
2497 return *(task_info.cell_partition_data.end() - 1) -
2498 *(task_info.cell_partition_data.end() - 2);
2499}
2500
2501
2502
2503template <int dim, typename Number, typename VectorizedArrayType>
2504inline unsigned int
2506{
2507 if (task_info.face_partition_data.empty())
2508 return 0;
2509 return task_info.face_partition_data.back();
2510}
2511
2512
2513
2514template <int dim, typename Number, typename VectorizedArrayType>
2515inline unsigned int
2517{
2518 if (task_info.face_partition_data.empty())
2519 return 0;
2520 return task_info.boundary_partition_data.back() -
2521 task_info.face_partition_data.back();
2522}
2523
2524
2525
2526template <int dim, typename Number, typename VectorizedArrayType>
2527inline unsigned int
2529{
2530 if (task_info.face_partition_data.empty())
2531 return 0;
2532 return face_info.faces.size() - task_info.boundary_partition_data.back();
2533}
2534
2535
2536
2537template <int dim, typename Number, typename VectorizedArrayType>
2538inline types::boundary_id
2540 const unsigned int face_batch_index) const
2541{
2542 Assert(face_batch_index >= task_info.boundary_partition_data[0] &&
2543 face_batch_index < task_info.boundary_partition_data.back(),
2544 ExcIndexRange(face_batch_index,
2545 task_info.boundary_partition_data[0],
2546 task_info.boundary_partition_data.back()));
2547 return types::boundary_id(face_info.faces[face_batch_index].exterior_face_no);
2548}
2549
2550
2551
2552template <int dim, typename Number, typename VectorizedArrayType>
2553inline std::array<types::boundary_id, VectorizedArrayType::size()>
2555 const unsigned int cell_batch_index,
2556 const unsigned int face_number) const
2557{
2558 AssertIndexRange(cell_batch_index, n_cell_batches());
2559 AssertIndexRange(face_number, ReferenceCells::max_n_faces<dim>());
2560 Assert(face_info.cell_and_face_boundary_id.size(0) >= n_cell_batches(),
2562 std::array<types::boundary_id, VectorizedArrayType::size()> result;
2563 result.fill(numbers::invalid_boundary_id);
2564 for (unsigned int v = 0;
2565 v < n_active_entries_per_cell_batch(cell_batch_index);
2566 ++v)
2567 result[v] =
2568 face_info.cell_and_face_boundary_id(cell_batch_index, face_number, v);
2569 return result;
2570}
2571
2572
2573
2574template <int dim, typename Number, typename VectorizedArrayType>
2575inline const internal::MatrixFreeFunctions::
2576 MappingInfo<dim, Number, VectorizedArrayType> &
2578{
2579 return mapping_info;
2580}
2581
2582
2583
2584template <int dim, typename Number, typename VectorizedArrayType>
2587 const unsigned int dof_index) const
2588{
2589 AssertIndexRange(dof_index, n_components());
2590 return dof_info[dof_index];
2591}
2592
2593
2594
2595template <int dim, typename Number, typename VectorizedArrayType>
2596inline unsigned int
2598{
2599 return constraint_pool_row_index.size() - 1;
2600}
2601
2602
2603
2604template <int dim, typename Number, typename VectorizedArrayType>
2605inline const Number *
2607 const unsigned int row) const
2608{
2609 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2610 return constraint_pool_data.empty() ?
2611 nullptr :
2612 constraint_pool_data.data() + constraint_pool_row_index[row];
2613}
2614
2615
2616
2617template <int dim, typename Number, typename VectorizedArrayType>
2618inline const Number *
2620 const unsigned int row) const
2621{
2622 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2623 return constraint_pool_data.empty() ?
2624 nullptr :
2625 constraint_pool_data.data() + constraint_pool_row_index[row + 1];
2626}
2627
2628
2629
2630template <int dim, typename Number, typename VectorizedArrayType>
2631inline std::pair<unsigned int, unsigned int>
2633 const std::pair<unsigned int, unsigned int> &range,
2634 const unsigned int degree,
2635 const unsigned int dof_handler_component) const
2636{
2637 if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2638 {
2640 dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2642 dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2643 if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2644 return range;
2645 else
2646 return {range.second, range.second};
2647 }
2648
2649 const unsigned int fe_index =
2650 dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2651 if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2652 return {range.second, range.second};
2653 else
2654 return create_cell_subrange_hp_by_index(range,
2655 fe_index,
2656 dof_handler_component);
2657}
2658
2659
2660
2661template <int dim, typename Number, typename VectorizedArrayType>
2662inline bool
2664 const unsigned int cell_batch_index) const
2665{
2666 AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2667 return VectorizedArrayType::size() > 1 &&
2668 cell_level_index[(cell_batch_index + 1) * VectorizedArrayType::size() -
2669 1] == cell_level_index[(cell_batch_index + 1) *
2670 VectorizedArrayType::size() -
2671 2];
2672}
2673
2674
2675
2676template <int dim, typename Number, typename VectorizedArrayType>
2677unsigned int
2679{
2680 return shape_info.size(2);
2681}
2682
2683
2684template <int dim, typename Number, typename VectorizedArrayType>
2685unsigned int
2687 const std::pair<unsigned int, unsigned int> range) const
2688{
2689 const auto &fe_indices =
2690 dof_info[first_hp_dof_handler_index].cell_active_fe_index;
2691
2692 if (fe_indices.empty() == true ||
2693 dof_handlers[first_hp_dof_handler_index]->get_fe_collection().size() == 1)
2694 return 0;
2695
2696 const auto index = fe_indices[range.first];
2697
2698 for (unsigned int i = range.first; i < range.second; ++i)
2699 AssertDimension(index, fe_indices[i]);
2700
2701 return index;
2702}
2703
2704
2705
2706template <int dim, typename Number, typename VectorizedArrayType>
2707unsigned int
2709 const std::pair<unsigned int, unsigned int> range,
2710 const bool is_interior_face) const
2711{
2712 const auto &fe_indices =
2713 dof_info[first_hp_dof_handler_index].cell_active_fe_index;
2714
2715 if (fe_indices.empty() == true)
2716 return 0;
2717
2718 if (is_interior_face)
2719 {
2720 const unsigned int index =
2721 fe_indices[face_info.faces[range.first].cells_interior[0] /
2722 VectorizedArrayType::size()];
2723
2724 for (unsigned int i = range.first; i < range.second; ++i)
2725 AssertDimension(index,
2726 fe_indices[face_info.faces[i].cells_interior[0] /
2727 VectorizedArrayType::size()]);
2728
2729 return index;
2730 }
2731 else
2732 {
2733 const unsigned int index =
2734 fe_indices[face_info.faces[range.first].cells_exterior[0] /
2735 VectorizedArrayType::size()];
2736
2737 for (unsigned int i = range.first; i < range.second; ++i)
2738 AssertDimension(index,
2739 fe_indices[face_info.faces[i].cells_exterior[0] /
2740 VectorizedArrayType::size()]);
2741
2742 return index;
2743 }
2744}
2745
2746
2747
2748template <int dim, typename Number, typename VectorizedArrayType>
2749inline unsigned int
2751 const unsigned int cell_batch_index) const
2752{
2753 Assert(!dof_info.empty(), ExcNotInitialized());
2754 AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2755 const std::vector<unsigned char> &n_lanes_filled =
2756 dof_info[0].n_vectorization_lanes_filled
2758 AssertIndexRange(cell_batch_index, n_lanes_filled.size());
2759
2760 return n_lanes_filled[cell_batch_index];
2761}
2762
2763
2764
2765template <int dim, typename Number, typename VectorizedArrayType>
2766inline unsigned int
2768 const unsigned int face_batch_index) const
2769{
2770 AssertIndexRange(face_batch_index, face_info.faces.size());
2771 Assert(!dof_info.empty(), ExcNotInitialized());
2772 const std::vector<unsigned char> &n_lanes_filled =
2773 dof_info[0].n_vectorization_lanes_filled
2775 AssertIndexRange(face_batch_index, n_lanes_filled.size());
2776 return n_lanes_filled[face_batch_index];
2777}
2778
2779
2780
2781template <int dim, typename Number, typename VectorizedArrayType>
2782inline unsigned int
2784 const unsigned int dof_handler_index,
2785 const unsigned int active_fe_index) const
2786{
2787 return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2788}
2789
2790
2791
2792template <int dim, typename Number, typename VectorizedArrayType>
2793inline unsigned int
2795 const unsigned int quad_index,
2796 const unsigned int active_fe_index) const
2797{
2798 AssertIndexRange(quad_index, mapping_info.cell_data.size());
2799 return mapping_info.cell_data[quad_index]
2800 .descriptor[active_fe_index]
2801 .n_q_points;
2802}
2803
2804
2805
2806template <int dim, typename Number, typename VectorizedArrayType>
2807inline unsigned int
2809 const unsigned int dof_handler_index,
2810 const unsigned int active_fe_index) const
2811{
2812 return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2813}
2814
2815
2816
2817template <int dim, typename Number, typename VectorizedArrayType>
2818inline unsigned int
2820 const unsigned int quad_index,
2821 const unsigned int active_fe_index) const
2822{
2823 AssertIndexRange(quad_index, mapping_info.face_data.size());
2824 return mapping_info.face_data[quad_index]
2825 .descriptor[active_fe_index]
2826 .n_q_points;
2827}
2828
2829
2830
2831template <int dim, typename Number, typename VectorizedArrayType>
2832inline const IndexSet &
2834 const unsigned int dof_handler_index) const
2835{
2836 return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2837}
2838
2839
2840
2841template <int dim, typename Number, typename VectorizedArrayType>
2842inline const IndexSet &
2844 const unsigned int dof_handler_index) const
2845{
2846 return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2847}
2848
2849
2850
2851template <int dim, typename Number, typename VectorizedArrayType>
2854 const unsigned int dof_handler_index,
2855 const unsigned int index_quad,
2856 const unsigned int index_fe,
2857 const unsigned int active_fe_index,
2858 const unsigned int active_quad_index) const
2859{
2860 AssertIndexRange(dof_handler_index, dof_info.size());
2861 const unsigned int ind =
2862 dof_info[dof_handler_index].global_base_element_offset + index_fe;
2863 AssertIndexRange(ind, shape_info.size(0));
2864 AssertIndexRange(index_quad, shape_info.size(1));
2865 AssertIndexRange(active_fe_index, shape_info.size(2));
2866 AssertIndexRange(active_quad_index, shape_info.size(3));
2867 return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2868}
2869
2870
2871
2872template <int dim, typename Number, typename VectorizedArrayType>
2874 VectorizedArrayType::size()> &
2876 const unsigned int face_batch_index) const
2877{
2878 AssertIndexRange(face_batch_index, face_info.faces.size());
2879 return face_info.faces[face_batch_index];
2880}
2881
2882
2883
2884template <int dim, typename Number, typename VectorizedArrayType>
2885inline const Table<3, unsigned int> &
2887 const
2888{
2889 return face_info.cell_and_face_to_plain_faces;
2890}
2891
2892
2893
2894template <int dim, typename Number, typename VectorizedArrayType>
2895inline const Quadrature<dim> &
2897 const unsigned int quad_index,
2898 const unsigned int active_fe_index) const
2899{
2900 AssertIndexRange(quad_index, mapping_info.cell_data.size());
2901 return mapping_info.cell_data[quad_index]
2902 .descriptor[active_fe_index]
2903 .quadrature;
2904}
2905
2906
2907
2908template <int dim, typename Number, typename VectorizedArrayType>
2909inline const Quadrature<dim - 1> &
2911 const unsigned int quad_index,
2912 const unsigned int active_fe_index) const
2913{
2914 AssertIndexRange(quad_index, mapping_info.face_data.size());
2915 return mapping_info.face_data[quad_index]
2916 .descriptor[active_fe_index]
2917 .quadrature;
2918}
2919
2920
2921
2922template <int dim, typename Number, typename VectorizedArrayType>
2923inline unsigned int
2925 const std::pair<unsigned int, unsigned int> range) const
2926{
2927 auto result = get_cell_category(range.first);
2928
2929 for (unsigned int i = range.first; i < range.second; ++i)
2930 result = std::max(result, get_cell_category(i));
2931
2932 return result;
2933}
2934
2935
2936
2937template <int dim, typename Number, typename VectorizedArrayType>
2938inline std::pair<unsigned int, unsigned int>
2940 const std::pair<unsigned int, unsigned int> range) const
2941{
2942 auto result = get_face_category(range.first);
2943
2944 for (unsigned int i = range.first; i < range.second; ++i)
2945 {
2946 result.first = std::max(result.first, get_face_category(i).first);
2947 result.second = std::max(result.second, get_face_category(i).second);
2948 }
2949
2950 return result;
2951}
2952
2953
2954
2955template <int dim, typename Number, typename VectorizedArrayType>
2956inline unsigned int
2958 const unsigned int cell_batch_index) const
2959{
2960 AssertIndexRange(0, dof_info.size());
2962 cell_batch_index,
2963 dof_info[first_hp_dof_handler_index].cell_active_fe_index.size());
2964 if (dof_info[first_hp_dof_handler_index].cell_active_fe_index.empty())
2965 return 0;
2966 else
2967 return dof_info[first_hp_dof_handler_index]
2968 .cell_active_fe_index[cell_batch_index];
2969}
2970
2971
2972
2973template <int dim, typename Number, typename VectorizedArrayType>
2974inline std::pair<unsigned int, unsigned int>
2976 const unsigned int face_batch_index) const
2977{
2978 AssertIndexRange(face_batch_index, face_info.faces.size());
2979 if (dof_info[first_hp_dof_handler_index].cell_active_fe_index.empty())
2980 return std::make_pair(0U, 0U);
2981
2982 std::pair<unsigned int, unsigned int> result = std::make_pair(0U, 0U);
2983 for (unsigned int v = 0;
2984 v < VectorizedArrayType::size() &&
2985 face_info.faces[face_batch_index].cells_interior[v] !=
2987 ++v)
2988 result.first =
2989 std::max(result.first,
2990 dof_info[first_hp_dof_handler_index].cell_active_fe_index
2991 [face_info.faces[face_batch_index].cells_interior[v] /
2992 VectorizedArrayType::size()]);
2993 if (face_info.faces[face_batch_index].cells_exterior[0] !=
2995 for (unsigned int v = 0;
2996 v < VectorizedArrayType::size() &&
2997 face_info.faces[face_batch_index].cells_exterior[v] !=
2999 ++v)
3000 result.second =
3001 std::max(result.second,
3002 dof_info[first_hp_dof_handler_index].cell_active_fe_index
3003 [face_info.faces[face_batch_index].cells_exterior[v] /
3004 VectorizedArrayType::size()]);
3005 else
3006 result.second = numbers::invalid_unsigned_int;
3007 return result;
3008}
3009
3010
3011
3012template <int dim, typename Number, typename VectorizedArrayType>
3013inline bool
3015{
3016 return indices_are_initialized;
3017}
3018
3019
3020
3021template <int dim, typename Number, typename VectorizedArrayType>
3022inline bool
3024{
3025 return mapping_is_initialized;
3026}
3027
3028
3029template <int dim, typename Number, typename VectorizedArrayType>
3030inline unsigned int
3032{
3033 return mg_level;
3034}
3035
3036
3037
3038template <int dim, typename Number, typename VectorizedArrayType>
3041{
3042 using list_type =
3043 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3044 list_type &data = scratch_pad.get();
3045 for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
3046 if (it->first == false)
3047 {
3048 it->first = true;
3049 return &it->second;
3050 }
3051 data.emplace_front(true, AlignedVector<VectorizedArrayType>());
3052 return &data.front().second;
3053}
3054
3055
3056
3057template <int dim, typename Number, typename VectorizedArrayType>
3058void
3060 const AlignedVector<VectorizedArrayType> *scratch) const
3061{
3062 using list_type =
3063 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3064 list_type &data = scratch_pad.get();
3065 for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
3066 if (&it->second == scratch)
3067 {
3068 Assert(it->first == true, ExcInternalError());
3069 it->first = false;
3070 return;
3071 }
3072 AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
3073}
3074
3075
3076
3077template <int dim, typename Number, typename VectorizedArrayType>
3081{
3082 for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
3083 scratch_pad_non_threadsafe.begin();
3084 it != scratch_pad_non_threadsafe.end();
3085 ++it)
3086 if (it->first == false)
3087 {
3088 it->first = true;
3089 return &it->second;
3090 }
3091 scratch_pad_non_threadsafe.push_front(
3092 std::make_pair(true, AlignedVector<Number>()));
3093 return &scratch_pad_non_threadsafe.front().second;
3094}
3095
3096
3097
3098template <int dim, typename Number, typename VectorizedArrayType>
3099void
3102 const AlignedVector<Number> *scratch) const
3103{
3104 for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
3105 scratch_pad_non_threadsafe.begin();
3106 it != scratch_pad_non_threadsafe.end();
3107 ++it)
3108 if (&it->second == scratch)
3109 {
3110 Assert(it->first == true, ExcInternalError());
3111 it->first = false;
3112 return;
3113 }
3114 AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
3115}
3116
3117
3118
3119// ------------------------------ reinit functions ---------------------------
3120
3121namespace internal
3122{
3123 namespace MatrixFreeImplementation
3124 {
3125 template <int dim, int spacedim>
3126 inline std::vector<IndexSet>
3127 extract_locally_owned_index_sets(
3128 const std::vector<const DoFHandler<dim, spacedim> *> &dofh,
3129 const unsigned int level)
3130 {
3131 std::vector<IndexSet> locally_owned_set;
3132 locally_owned_set.reserve(dofh.size());
3133 for (unsigned int j = 0; j < dofh.size(); ++j)
3135 locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
3136 else
3137 locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
3138 return locally_owned_set;
3139 }
3140 } // namespace MatrixFreeImplementation
3141} // namespace internal
3142
3143
3144
3145template <int dim, typename Number, typename VectorizedArrayType>
3146template <typename QuadratureType, typename number2, typename MappingType>
3147void
3149 const MappingType &mapping,
3150 const DoFHandler<dim> &dof_handler,
3151 const AffineConstraints<number2> &constraints_in,
3152 const QuadratureType &quad,
3154 &additional_data)
3155{
3156 std::vector<const DoFHandler<dim, dim> *> dof_handlers;
3157 std::vector<const AffineConstraints<number2> *> constraints;
3158
3159 dof_handlers.push_back(&dof_handler);
3160 constraints.push_back(&constraints_in);
3161
3162 std::vector<IndexSet> locally_owned_sets =
3163 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3164 dof_handlers, additional_data.mg_level);
3165
3166 std::vector<hp::QCollection<dim>> quad_hp;
3167 quad_hp.emplace_back(quad);
3168
3169 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3170 dof_handlers,
3171 constraints,
3172 locally_owned_sets,
3173 quad_hp,
3174 additional_data);
3175}
3176
3177
3178
3179template <int dim, typename Number, typename VectorizedArrayType>
3180template <typename QuadratureType, typename number2, typename MappingType>
3181void
3183 const MappingType &mapping,
3184 const std::vector<const DoFHandler<dim> *> &dof_handler,
3185 const std::vector<const AffineConstraints<number2> *> &constraint,
3186 const QuadratureType &quad,
3188 &additional_data)
3189{
3190 std::vector<IndexSet> locally_owned_set =
3191 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3192 dof_handler, additional_data.mg_level);
3193 std::vector<hp::QCollection<dim>> quad_hp;
3194 quad_hp.emplace_back(quad);
3195
3196 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3197 dof_handler,
3198 constraint,
3199 locally_owned_set,
3200 quad_hp,
3201 additional_data);
3202}
3203
3204
3205
3206template <int dim, typename Number, typename VectorizedArrayType>
3207template <typename QuadratureType, typename number2, typename MappingType>
3208void
3210 const MappingType &mapping,
3211 const std::vector<const DoFHandler<dim> *> &dof_handler,
3212 const std::vector<const AffineConstraints<number2> *> &constraint,
3213 const std::vector<QuadratureType> &quad,
3215 &additional_data)
3216{
3217 std::vector<IndexSet> locally_owned_set =
3218 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3219 dof_handler, additional_data.mg_level);
3220 std::vector<hp::QCollection<dim>> quad_hp;
3221 for (unsigned int q = 0; q < quad.size(); ++q)
3222 quad_hp.emplace_back(quad[q]);
3223
3224 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3225 dof_handler,
3226 constraint,
3227 locally_owned_set,
3228 quad_hp,
3229 additional_data);
3230}
3231
3232
3233
3234// ------------------------------ implementation of loops --------------------
3235
3236// internal helper functions that define how to call MPI data exchange
3237// functions: for generic vectors, do nothing at all. For distributed vectors,
3238// call update_ghost_values_start function and so on. If we have collections
3239// of vectors, just do the individual functions of the components. In order to
3240// keep ghost values consistent (whether we are in read or write mode), we
3241// also reset the values at the end. the whole situation is a bit complicated
3242// by the fact that we need to treat block vectors differently, which use some
3243// additional helper functions to select the blocks and template magic.
3244namespace internal
3245{
3249 template <int dim, typename Number, typename VectorizedArrayType>
3250 struct VectorDataExchange
3251 {
3252 // A shift for the MPI messages to reduce the risk for accidental
3253 // interaction with other open communications that a user program might
3254 // set up (parallel vectors support unfinished communication). We let
3255 // the other vectors use the first 20 assigned numbers and start the
3256 // matrix-free communication.
3257 static constexpr unsigned int channel_shift = 20;
3258
3259
3260
3265 VectorDataExchange(
3266 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3267 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3268 DataAccessOnFaces vector_face_access,
3269 const unsigned int n_components)
3270 : matrix_free(matrix_free)
3271 , vector_face_access(
3272 matrix_free.get_task_info().face_partition_data.empty() ?
3273 ::MatrixFree<dim, Number, VectorizedArrayType>::
3274 DataAccessOnFaces::unspecified :
3275 vector_face_access)
3276 , ghosts_were_set(false)
3277# ifdef DEAL_II_WITH_MPI
3278 , tmp_data(n_components)
3279 , requests(n_components)
3280# endif
3281 {
3282 (void)n_components;
3283 if (this->vector_face_access !=
3285 DataAccessOnFaces::unspecified)
3286 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3288 matrix_free.get_dof_info(c).vector_exchanger_face_variants.size(),
3289 5);
3290 }
3291
3292
3293
3297 ~VectorDataExchange() // NOLINT
3298 {
3299# ifdef DEAL_II_WITH_MPI
3300 for (unsigned int i = 0; i < tmp_data.size(); ++i)
3301 if (tmp_data[i] != nullptr)
3302 matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3303# endif
3304 }
3305
3306
3307
3312 template <typename VectorType>
3313 unsigned int
3314 find_vector_in_mf(const VectorType &vec,
3315 const bool check_global_compatibility = true) const
3316 {
3317 // case 1: vector was set up with MatrixFree::initialize_dof_vector()
3318 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3319 if (vec.get_partitioner().get() ==
3320 matrix_free.get_dof_info(c).vector_partitioner.get())
3321 return c;
3322
3323 // case 2: user provided own partitioner (compatibility mode)
3324 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3325 if (check_global_compatibility ?
3326 vec.get_partitioner()->is_globally_compatible(
3327 *matrix_free.get_dof_info(c).vector_partitioner) :
3328 vec.get_partitioner()->is_compatible(
3329 *matrix_free.get_dof_info(c).vector_partitioner))
3330 return c;
3331
3332 Assert(false,
3333 ExcNotImplemented("Could not find partitioner that fits vector"));
3334
3336 }
3337
3338
3339
3345 get_partitioner(const unsigned int mf_component) const
3346 {
3347 AssertDimension(matrix_free.get_dof_info(mf_component)
3348 .vector_exchanger_face_variants.size(),
3349 5);
3350 if (vector_face_access ==
3352 DataAccessOnFaces::none)
3353 return *matrix_free.get_dof_info(mf_component)
3354 .vector_exchanger_face_variants[0];
3355 else if (vector_face_access ==
3357 DataAccessOnFaces::values)
3358 return *matrix_free.get_dof_info(mf_component)
3359 .vector_exchanger_face_variants[1];
3360 else if (vector_face_access ==
3362 DataAccessOnFaces::gradients)
3363 return *matrix_free.get_dof_info(mf_component)
3364 .vector_exchanger_face_variants[2];
3365 else if (vector_face_access ==
3367 DataAccessOnFaces::values_all_faces)
3368 return *matrix_free.get_dof_info(mf_component)
3369 .vector_exchanger_face_variants[3];
3370 else if (vector_face_access ==
3372 DataAccessOnFaces::gradients_all_faces)
3373 return *matrix_free.get_dof_info(mf_component)
3374 .vector_exchanger_face_variants[4];
3375 else
3376 return *matrix_free.get_dof_info(mf_component).vector_exchanger.get();
3377 }
3378
3379
3380
3384 template <typename VectorType,
3385 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3386 * = nullptr>
3387 void
3388 update_ghost_values_start(const unsigned int /*component_in_block_vector*/,
3389 const VectorType & /*vec*/)
3390 {}
3391
3392
3397 template <typename VectorType,
3398 std::enable_if_t<!has_update_ghost_values_start<VectorType> &&
3399 !is_not_parallel_vector<VectorType>,
3400 VectorType> * = nullptr>
3401 void
3402 update_ghost_values_start(const unsigned int component_in_block_vector,
3403 const VectorType &vec)
3404 {
3405 (void)component_in_block_vector;
3406 const bool ghosts_set = vec.has_ghost_elements();
3407
3408 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3409 ghosts_set == false,
3411
3412 if (ghosts_set)
3413 {
3414 ghosts_were_set = true;
3415 return;
3416 }
3417
3418 vec.update_ghost_values();
3419 }
3420
3421
3422
3428 template <typename VectorType,
3429 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3430 !has_exchange_on_subset<VectorType>,
3431 VectorType> * = nullptr>
3432 void
3433 update_ghost_values_start(const unsigned int component_in_block_vector,
3434 const VectorType &vec)
3435 {
3436 (void)component_in_block_vector;
3437 const bool ghosts_set = vec.has_ghost_elements();
3438
3439 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3440 ghosts_set == false,
3442
3443 if (ghosts_set)
3444 {
3445 ghosts_were_set = true;
3446 return;
3447 }
3448
3449 vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3450 }
3451
3452
3453
3460 template <typename VectorType,
3461 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3462 has_exchange_on_subset<VectorType>,
3463 VectorType> * = nullptr>
3464 void
3465 update_ghost_values_start(const unsigned int component_in_block_vector,
3466 const VectorType &vec)
3467 {
3468 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3469 "Type mismatch between VectorType and VectorDataExchange");
3470 (void)component_in_block_vector;
3471 const bool ghosts_set = vec.has_ghost_elements();
3472
3473 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3474 ghosts_set == false,
3476
3477 if (ghosts_set)
3478 {
3479 ghosts_were_set = true;
3480 return;
3481 }
3482
3483 if (vec.size() != 0)
3484 {
3485# ifdef DEAL_II_WITH_MPI
3486 const unsigned int mf_component = find_vector_in_mf(vec);
3487
3488 const auto &part = get_partitioner(mf_component);
3489
3490 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3491 part.n_import_sm_procs() == 0)
3492 return;
3493
3494 tmp_data[component_in_block_vector] =
3495 matrix_free.acquire_scratch_data_non_threadsafe();
3496 tmp_data[component_in_block_vector]->resize_fast(
3497 part.n_import_indices());
3498 AssertDimension(requests.size(), tmp_data.size());
3499
3500 part.export_to_ghosted_array_start(
3501 component_in_block_vector * 2 + channel_shift,
3502 ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3503 vec.shared_vector_data(),
3504 ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3505 part.locally_owned_size(),
3506 matrix_free.get_dof_info(mf_component)
3507 .vector_partitioner->n_ghost_indices()),
3508 ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3509 part.n_import_indices()),
3510 this->requests[component_in_block_vector]);
3511# endif
3512 }
3513 }
3514
3515
3516
3521 template <typename VectorType,
3522 std::enable_if_t<!has_update_ghost_values_start<VectorType>,
3523 VectorType> * = nullptr>
3524 void
3525 update_ghost_values_finish(const unsigned int /*component_in_block_vector*/,
3526 const VectorType & /*vec*/)
3527 {}
3528
3529
3530
3536 template <typename VectorType,
3537 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3538 !has_exchange_on_subset<VectorType>,
3539 VectorType> * = nullptr>
3540 void
3541 update_ghost_values_finish(const unsigned int component_in_block_vector,
3542 const VectorType &vec)
3543 {
3544 (void)component_in_block_vector;
3545
3546 if (ghosts_were_set)
3547 return;
3548
3549 vec.update_ghost_values_finish();
3550 }
3551
3552
3553
3560 template <typename VectorType,
3561 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3562 has_exchange_on_subset<VectorType>,
3563 VectorType> * = nullptr>
3564 void
3565 update_ghost_values_finish(const unsigned int component_in_block_vector,
3566 const VectorType &vec)
3567 {
3568 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3569 "Type mismatch between VectorType and VectorDataExchange");
3570 (void)component_in_block_vector;
3571
3572 if (ghosts_were_set)
3573 return;
3574
3575 if (vec.size() != 0)
3576 {
3577# ifdef DEAL_II_WITH_MPI
3578 AssertIndexRange(component_in_block_vector, tmp_data.size());
3579 AssertDimension(requests.size(), tmp_data.size());
3580
3581 const unsigned int mf_component = find_vector_in_mf(vec);
3582
3583 const auto &part = get_partitioner(mf_component);
3584
3585 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3586 part.n_import_sm_procs() != 0)
3587 {
3588 part.export_to_ghosted_array_finish(
3589 ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3590 vec.shared_vector_data(),
3591 ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3592 part.locally_owned_size(),
3593 matrix_free.get_dof_info(mf_component)
3594 .vector_partitioner->n_ghost_indices()),
3595 this->requests[component_in_block_vector]);
3596
3597 matrix_free.release_scratch_data_non_threadsafe(
3598 tmp_data[component_in_block_vector]);
3599 tmp_data[component_in_block_vector] = nullptr;
3600 }
3601# endif
3602 }
3603 // let vector know that ghosts are being updated and we can read from
3604 // them
3605 vec.set_ghost_state(true);
3606 }
3607
3608
3609
3613 template <typename VectorType,
3614 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3615 * = nullptr>
3616 void
3617 compress_start(const unsigned int /*component_in_block_vector*/,
3618 VectorType & /*vec*/)
3619 {}
3620
3621
3622
3627 template <typename VectorType,
3628 std::enable_if_t<!has_compress_start<VectorType> &&
3629 !is_not_parallel_vector<VectorType>,
3630 VectorType> * = nullptr>
3631 void
3632 compress_start(const unsigned int component_in_block_vector,
3633 VectorType &vec)
3634 {
3635 (void)component_in_block_vector;
3636 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3637 vec.compress(VectorOperation::add);
3638 }
3639
3640
3641
3647 template <typename VectorType,
3648 std::enable_if_t<has_compress_start<VectorType> &&
3649 !has_exchange_on_subset<VectorType>,
3650 VectorType> * = nullptr>
3651 void
3652 compress_start(const unsigned int component_in_block_vector,
3653 VectorType &vec)
3654 {
3655 (void)component_in_block_vector;
3656 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3657 vec.compress_start(component_in_block_vector + channel_shift);
3658 }
3659
3660
3661
3668 template <typename VectorType,
3669 std::enable_if_t<has_compress_start<VectorType> &&
3670 has_exchange_on_subset<VectorType>,
3671 VectorType> * = nullptr>
3672 void
3673 compress_start(const unsigned int component_in_block_vector,
3674 VectorType &vec)
3675 {
3676 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3677 "Type mismatch between VectorType and VectorDataExchange");
3678 (void)component_in_block_vector;
3679 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3680
3681 if (vec.size() != 0)
3682 {
3683# ifdef DEAL_II_WITH_MPI
3684 const unsigned int mf_component = find_vector_in_mf(vec);
3685
3686 const auto &part = get_partitioner(mf_component);
3687
3688 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3689 part.n_import_sm_procs() == 0)
3690 return;
3691
3692 tmp_data[component_in_block_vector] =
3693 matrix_free.acquire_scratch_data_non_threadsafe();
3694 tmp_data[component_in_block_vector]->resize_fast(
3695 part.n_import_indices());
3696 AssertDimension(requests.size(), tmp_data.size());
3697
3698 part.import_from_ghosted_array_start(
3700 component_in_block_vector * 2 + channel_shift,
3701 ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3702 vec.shared_vector_data(),
3703 ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3704 matrix_free.get_dof_info(mf_component)
3705 .vector_partitioner->n_ghost_indices()),
3706 ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3707 part.n_import_indices()),
3708 this->requests[component_in_block_vector]);
3709# endif
3710 }
3711 }
3712
3713
3714
3719 template <
3720 typename VectorType,
3721 std::enable_if_t<!has_compress_start<VectorType>, VectorType> * = nullptr>
3722 void
3723 compress_finish(const unsigned int /*component_in_block_vector*/,
3724 VectorType & /*vec*/)
3725 {}
3726
3727
3728
3734 template <typename VectorType,
3735 std::enable_if_t<has_compress_start<VectorType> &&
3736 !has_exchange_on_subset<VectorType>,
3737 VectorType> * = nullptr>
3738 void
3739 compress_finish(const unsigned int component_in_block_vector,
3740 VectorType &vec)
3741 {
3742 (void)component_in_block_vector;
3743 vec.compress_finish(VectorOperation::add);
3744 }
3745
3746
3747
3754 template <typename VectorType,
3755 std::enable_if_t<has_compress_start<VectorType> &&
3756 has_exchange_on_subset<VectorType>,
3757 VectorType> * = nullptr>
3758 void
3759 compress_finish(const unsigned int component_in_block_vector,
3760 VectorType &vec)
3761 {
3762 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3763 "Type mismatch between VectorType and VectorDataExchange");
3764 (void)component_in_block_vector;
3765 if (vec.size() != 0)
3766 {
3767# ifdef DEAL_II_WITH_MPI
3768 AssertIndexRange(component_in_block_vector, tmp_data.size());
3769 AssertDimension(requests.size(), tmp_data.size());
3770
3771 const unsigned int mf_component = find_vector_in_mf(vec);
3772
3773 const auto &part = get_partitioner(mf_component);
3774
3775 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3776 part.n_import_sm_procs() != 0)
3777 {
3778 part.import_from_ghosted_array_finish(
3780 ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3781 vec.shared_vector_data(),
3782 ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3783 matrix_free.get_dof_info(mf_component)
3784 .vector_partitioner->n_ghost_indices()),
3786 tmp_data[component_in_block_vector]->begin(),
3787 part.n_import_indices()),
3788 this->requests[component_in_block_vector]);
3789
3790 matrix_free.release_scratch_data_non_threadsafe(
3791 tmp_data[component_in_block_vector]);
3792 tmp_data[component_in_block_vector] = nullptr;
3793 }
3794
3796 {
3797 const int ierr =
3798 MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3799 AssertThrowMPI(ierr);
3800 }
3801# endif
3802 }
3803 }
3804
3805
3806
3810 template <typename VectorType,
3811 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3812 * = nullptr>
3813 void
3814 reset_ghost_values(const VectorType & /*vec*/) const
3815 {}
3816
3817
3818
3823 template <typename VectorType,
3824 std::enable_if_t<!has_exchange_on_subset<VectorType> &&
3825 !is_not_parallel_vector<VectorType>,
3826 VectorType> * = nullptr>
3827 void
3828 reset_ghost_values(const VectorType &vec) const
3829 {
3830 if (ghosts_were_set == true)
3831 return;
3832
3833 vec.zero_out_ghost_values();
3834 }
3835
3836
3837
3843 template <typename VectorType,
3844 std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3845 * = nullptr>
3846 void
3847 reset_ghost_values(const VectorType &vec) const
3848 {
3849 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3850 "Type mismatch between VectorType and VectorDataExchange");
3851 if (ghosts_were_set == true)
3852 return;
3853
3854 if (vec.size() != 0)
3855 {
3856# ifdef DEAL_II_WITH_MPI
3857 AssertDimension(requests.size(), tmp_data.size());
3858
3859 const unsigned int mf_component = find_vector_in_mf(vec);
3860
3861 const auto &part = get_partitioner(mf_component);
3862
3863 if (part.n_ghost_indices() > 0)
3864 {
3865 part.reset_ghost_values(
3866 ArrayView<Number>(const_cast<VectorType &>(vec).begin() +
3867 part.locally_owned_size(),
3868 matrix_free.get_dof_info(mf_component)
3869 .vector_partitioner->n_ghost_indices()));
3870 }
3871
3872# endif
3873 }
3874 // let vector know that it's not ghosted anymore
3875 vec.set_ghost_state(false);
3876 }
3877
3878
3879
3885 template <typename VectorType,
3886 std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3887 * = nullptr>
3888 void
3889 zero_vector_region(const unsigned int range_index, VectorType &vec) const
3890 {
3891 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3892 "Type mismatch between VectorType and VectorDataExchange");
3893 if (range_index == numbers::invalid_unsigned_int)
3894 vec = Number();
3895 else
3896 {
3897 const unsigned int mf_component = find_vector_in_mf(vec, false);
3899 matrix_free.get_dof_info(mf_component);
3900 Assert(dof_info.vector_zero_range_list_index.empty() == false,
3902
3903 Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3905 AssertIndexRange(range_index,
3906 dof_info.vector_zero_range_list_index.size() - 1);
3907 for (unsigned int id =
3908 dof_info.vector_zero_range_list_index[range_index];
3909 id != dof_info.vector_zero_range_list_index[range_index + 1];
3910 ++id)
3911 std::memset(vec.begin() + dof_info.vector_zero_range_list[id].first,
3912 0,
3913 (dof_info.vector_zero_range_list[id].second -
3914 dof_info.vector_zero_range_list[id].first) *
3915 sizeof(Number));
3916 }
3917 }
3918
3919
3920
3926 template <typename VectorType,
3927 std::enable_if_t<!has_exchange_on_subset<VectorType>, VectorType>
3928 * = nullptr,
3929 std::enable_if_t<has_assignment_operator<VectorType>, VectorType>
3930 * = nullptr>
3931 void
3932 zero_vector_region(const unsigned int range_index, VectorType &vec) const
3933 {
3934 if (range_index == numbers::invalid_unsigned_int || range_index == 0)
3935 {
3936 if constexpr (std::is_same_v<
3938 VectorType>)
3939 {
3940 for (unsigned int i = 0; i < vec.size(); ++i)
3941 vec[i] = typename VectorType::value_type();
3942 }
3943 else
3944 vec = typename VectorType::value_type();
3945 }
3946 }
3947
3948
3949
3954 void
3955 zero_vector_region(const unsigned int, ...) const
3956 {
3957 Assert(false,
3959 "Zeroing is only implemented for vector types "
3960 "which provide operator=(const VectorType::value_type)"));
3961 }
3962
3963
3964
3965 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
3966 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3967 DataAccessOnFaces vector_face_access;
3968 bool ghosts_were_set;
3969# ifdef DEAL_II_WITH_MPI
3970 std::vector<AlignedVector<Number> *> tmp_data;
3971 std::vector<std::vector<MPI_Request>> requests;
3972# endif
3973 }; // VectorDataExchange
3974
3975 template <typename VectorStruct>
3976 unsigned int
3977 n_components(const VectorStruct &vec);
3978
3979 template <typename VectorStruct>
3980 unsigned int
3981 n_components_block(const VectorStruct &vec, const std::bool_constant<true>)
3982 {
3983 unsigned int components = 0;
3984 for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
3985 components += n_components(vec.block(bl));
3986 return components;
3987 }
3988
3989 template <typename VectorStruct>
3990 unsigned int
3991 n_components_block(const VectorStruct &, const std::bool_constant<false>)
3992 {
3993 return 1;
3994 }
3995
3996 template <typename VectorStruct>
3997 unsigned int
3998 n_components(const VectorStruct &vec)
3999 {
4000 return n_components_block(
4001 vec, std::bool_constant<IsBlockVector<VectorStruct>::value>());
4002 }
4003
4004 template <typename VectorStruct>
4005 inline unsigned int
4006 n_components(const std::vector<VectorStruct> &vec)
4007 {
4008 unsigned int components = 0;
4009 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4010 components += n_components_block(
4011 vec[comp], std::bool_constant<IsBlockVector<VectorStruct>::value>());
4012 return components;
4013 }
4014
4015 template <typename VectorStruct>
4016 inline unsigned int
4017 n_components(const std::vector<VectorStruct *> &vec)
4018 {
4019 unsigned int components = 0;
4020 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4021 components += n_components_block(
4022 *vec[comp], std::bool_constant<IsBlockVector<VectorStruct>::value>());
4023 return components;
4024 }
4025
4026
4027
4028 // A helper function to identify block vectors with many components where we
4029 // should not try to overlap computations and communication because there
4030 // would be too many outstanding communication requests.
4031
4032 // default value for vectors that do not have communication_block_size
4033 template <typename VectorStruct,
4034 std::enable_if_t<!has_communication_block_size<VectorStruct>,
4035 VectorStruct> * = nullptr>
4036 constexpr unsigned int
4037 get_communication_block_size(const VectorStruct &)
4038 {
4040 }
4041
4042
4043
4044 template <typename VectorStruct,
4045 std::enable_if_t<has_communication_block_size<VectorStruct>,
4046 VectorStruct> * = nullptr>
4047 constexpr unsigned int
4048 get_communication_block_size(const VectorStruct &)
4049 {
4050 return VectorStruct::communication_block_size;
4051 }
4052
4053
4054
4055 template <typename VectorType,
4056 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType> * =
4057 nullptr>
4058 bool
4059 has_ghost_elements(const VectorType &vec)
4060 {
4061 (void)vec;
4062 return false;
4063 }
4064
4065
4066
4067 template <typename VectorType,
4068 std::enable_if_t<!is_not_parallel_vector<VectorType>, VectorType>
4069 * = nullptr>
4070 bool
4071 has_ghost_elements(const VectorType &vec)
4072 {
4073 return vec.has_ghost_elements();
4074 }
4075
4076
4077
4078 // --------------------------------------------------------------------------
4079 // below we have wrappers to distinguish between block and non-block vectors.
4080 // --------------------------------------------------------------------------
4081
4082 //
4083 // update_ghost_values_start
4084 //
4085
4086 // update_ghost_values for block vectors
4087 template <int dim,
4088 typename VectorStruct,
4089 typename Number,
4090 typename VectorizedArrayType,
4091 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4092 * = nullptr>
4093 void
4094 update_ghost_values_start(
4095 const VectorStruct &vec,
4096 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4097 const unsigned int channel = 0)
4098 {
4099 if (get_communication_block_size(vec) < vec.n_blocks())
4100 {
4101 const bool ghosts_set = vec.has_ghost_elements();
4102
4103 Assert(exchanger.matrix_free.get_task_info()
4104 .allow_ghosted_vectors_in_loops ||
4105 ghosts_set == false,
4107
4108 if (ghosts_set)
4109 {
4110 exchanger.ghosts_were_set = true;
4111 return;
4112 }
4113
4114 vec.update_ghost_values();
4115 }
4116 else
4117 {
4118 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4119 update_ghost_values_start(vec.block(i), exchanger, channel + i);
4120 }
4121 }
4122
4123
4124
4125 // update_ghost_values for non-block vectors
4126 template <int dim,
4127 typename VectorStruct,
4128 typename Number,
4129 typename VectorizedArrayType,
4130 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4131 * = nullptr>
4132 void
4133 update_ghost_values_start(
4134 const VectorStruct &vec,
4135 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4136 const unsigned int channel = 0)
4137 {
4138 exchanger.update_ghost_values_start(channel, vec);
4139 }
4140
4141
4142
4143 // update_ghost_values_start() for vector of vectors
4144 template <int dim,
4145 typename VectorStruct,
4146 typename Number,
4147 typename VectorizedArrayType>
4148 inline void
4149 update_ghost_values_start(
4150 const std::vector<VectorStruct> &vec,
4151 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4152 {
4153 unsigned int component_index = 0;
4154 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4155 {
4156 update_ghost_values_start(vec[comp], exchanger, component_index);
4157 component_index += n_components(vec[comp]);
4158 }
4159 }
4160
4161
4162
4163 // update_ghost_values_start() for vector of pointers to vectors
4164 template <int dim,
4165 typename VectorStruct,
4166 typename Number,
4167 typename VectorizedArrayType>
4168 inline void
4169 update_ghost_values_start(
4170 const std::vector<VectorStruct *> &vec,
4171 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4172 {
4173 unsigned int component_index = 0;
4174 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4175 {
4176 update_ghost_values_start(*vec[comp], exchanger, component_index);
4177 component_index += n_components(*vec[comp]);
4178 }
4179 }
4180
4181
4182
4183 //
4184 // update_ghost_values_finish
4185 //
4186
4187 // for block vectors
4188 template <int dim,
4189 typename VectorStruct,
4190 typename Number,
4191 typename VectorizedArrayType,
4192 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4193 * = nullptr>
4194 void
4195 update_ghost_values_finish(
4196 const VectorStruct &vec,
4197 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4198 const unsigned int channel = 0)
4199 {
4200 if (get_communication_block_size(vec) < vec.n_blocks())
4201 {
4202 // do nothing, everything has already been completed in the _start()
4203 // call
4204 }
4205 else
4206 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4207 update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4208 }
4209
4210
4211
4212 // for non-block vectors
4213 template <int dim,
4214 typename VectorStruct,
4215 typename Number,
4216 typename VectorizedArrayType,
4217 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4218 * = nullptr>
4219 void
4220 update_ghost_values_finish(
4221 const VectorStruct &vec,
4222 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4223 const unsigned int channel = 0)
4224 {
4225 exchanger.update_ghost_values_finish(channel, vec);
4226 }
4227
4228
4229
4230 // for vector of vectors
4231 template <int dim,
4232 typename VectorStruct,
4233 typename Number,
4234 typename VectorizedArrayType>
4235 inline void
4236 update_ghost_values_finish(
4237 const std::vector<VectorStruct> &vec,
4238 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4239 {
4240 unsigned int component_index = 0;
4241 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4242 {
4243 update_ghost_values_finish(vec[comp], exchanger, component_index);
4244 component_index += n_components(vec[comp]);
4245 }
4246 }
4247
4248
4249
4250 // for vector of pointers to vectors
4251 template <int dim,
4252 typename VectorStruct,
4253 typename Number,
4254 typename VectorizedArrayType>
4255 inline void
4256 update_ghost_values_finish(
4257 const std::vector<VectorStruct *> &vec,
4258 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4259 {
4260 unsigned int component_index = 0;
4261 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4262 {
4263 update_ghost_values_finish(*vec[comp], exchanger, component_index);
4264 component_index += n_components(*vec[comp]);
4265 }
4266 }
4267
4268
4269
4270 //
4271 // compress_start
4272 //
4273
4274 // for block vectors
4275 template <int dim,
4276 typename VectorStruct,
4277 typename Number,
4278 typename VectorizedArrayType,
4279 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4280 * = nullptr>
4281 inline void
4282 compress_start(
4283 VectorStruct &vec,
4284 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4285 const unsigned int channel = 0)
4286 {
4287 if (get_communication_block_size(vec) < vec.n_blocks())
4288 vec.compress(VectorOperation::add);
4289 else
4290 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4291 compress_start(vec.block(i), exchanger, channel + i);
4292 }
4293
4294
4295
4296 // for non-block vectors
4297 template <int dim,
4298 typename VectorStruct,
4299 typename Number,
4300 typename VectorizedArrayType,
4301 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4302 * = nullptr>
4303 inline void
4304 compress_start(
4305 VectorStruct &vec,
4306 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4307 const unsigned int channel = 0)
4308 {
4309 exchanger.compress_start(channel, vec);
4310 }
4311
4312
4313
4314 // for std::vector of vectors
4315 template <int dim,
4316 typename VectorStruct,
4317 typename Number,
4318 typename VectorizedArrayType>
4319 inline void
4320 compress_start(
4321 std::vector<VectorStruct> &vec,
4322 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4323 {
4324 unsigned int component_index = 0;
4325 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4326 {
4327 compress_start(vec[comp], exchanger, component_index);
4328 component_index += n_components(vec[comp]);
4329 }
4330 }
4331
4332
4333
4334 // for std::vector of pointer to vectors
4335 template <int dim,
4336 typename VectorStruct,
4337 typename Number,
4338 typename VectorizedArrayType>
4339 inline void
4340 compress_start(
4341 std::vector<VectorStruct *> &vec,
4342 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4343 {
4344 unsigned int component_index = 0;
4345 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4346 {
4347 compress_start(*vec[comp], exchanger, component_index);
4348 component_index += n_components(*vec[comp]);
4349 }
4350 }
4351
4352
4353
4354 //
4355 // compress_finish
4356 //
4357
4358 // for block vectors
4359 template <int dim,
4360 typename VectorStruct,
4361 typename Number,
4362 typename VectorizedArrayType,
4363 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4364 * = nullptr>
4365 inline void
4366 compress_finish(
4367 VectorStruct &vec,
4368 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4369 const unsigned int channel = 0)
4370 {
4371 if (get_communication_block_size(vec) < vec.n_blocks())
4372 {
4373 // do nothing, everything has already been completed in the _start()
4374 // call
4375 }
4376 else
4377 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4378 compress_finish(vec.block(i), exchanger, channel + i);
4379 }
4380
4381
4382
4383 // for non-block vectors
4384 template <int dim,
4385 typename VectorStruct,
4386 typename Number,
4387 typename VectorizedArrayType,
4388 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4389 * = nullptr>
4390 inline void
4391 compress_finish(
4392 VectorStruct &vec,
4393 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4394 const unsigned int channel = 0)
4395 {
4396 exchanger.compress_finish(channel, vec);
4397 }
4398
4399
4400
4401 // for std::vector of vectors
4402 template <int dim,
4403 typename VectorStruct,
4404 typename Number,
4405 typename VectorizedArrayType>
4406 inline void
4407 compress_finish(
4408 std::vector<VectorStruct> &vec,
4409 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4410 {
4411 unsigned int component_index = 0;
4412 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4413 {
4414 compress_finish(vec[comp], exchanger, component_index);
4415 component_index += n_components(vec[comp]);
4416 }
4417 }
4418
4419
4420
4421 // for std::vector of pointer to vectors
4422 template <int dim,
4423 typename VectorStruct,
4424 typename Number,
4425 typename VectorizedArrayType>
4426 inline void
4427 compress_finish(
4428 std::vector<VectorStruct *> &vec,
4429 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4430 {
4431 unsigned int component_index = 0;
4432 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4433 {
4434 compress_finish(*vec[comp], exchanger, component_index);
4435 component_index += n_components(*vec[comp]);
4436 }
4437 }
4438
4439
4440
4441 //
4442 // reset_ghost_values:
4443 //
4444 // if the input vector did not have ghosts imported, clear them here again
4445 // in order to avoid subsequent operations e.g. in linear solvers to work
4446 // with ghosts all the time
4447 //
4448
4449 // for block vectors
4450 template <int dim,
4451 typename VectorStruct,
4452 typename Number,
4453 typename VectorizedArrayType,
4454 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4455 * = nullptr>
4456 inline void
4457 reset_ghost_values(
4458 const VectorStruct &vec,
4459 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4460 {
4461 // return immediately if there is nothing to do.
4462 if (exchanger.ghosts_were_set == true)
4463 return;
4464
4465 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4466 reset_ghost_values(vec.block(i), exchanger);
4467 }
4468
4469
4470
4471 // for non-block vectors
4472 template <int dim,
4473 typename VectorStruct,
4474 typename Number,
4475 typename VectorizedArrayType,
4476 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4477 * = nullptr>
4478 inline void
4479 reset_ghost_values(
4480 const VectorStruct &vec,
4481 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4482 {
4483 // return immediately if there is nothing to do.
4484 if (exchanger.ghosts_were_set == true)
4485 return;
4486
4487 exchanger.reset_ghost_values(vec);
4488 }
4489
4490
4491
4492 // for std::vector of vectors
4493 template <int dim,
4494 typename VectorStruct,
4495 typename Number,
4496 typename VectorizedArrayType>
4497 inline void
4498 reset_ghost_values(
4499 const std::vector<VectorStruct> &vec,
4500 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4501 {
4502 // return immediately if there is nothing to do.
4503 if (exchanger.ghosts_were_set == true)
4504 return;
4505
4506 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4507 reset_ghost_values(vec[comp], exchanger);
4508 }
4509
4510
4511
4512 // for std::vector of pointer to vectors
4513 template <int dim,
4514 typename VectorStruct,
4515 typename Number,
4516 typename VectorizedArrayType>
4517 inline void
4518 reset_ghost_values(
4519 const std::vector<VectorStruct *> &vec,
4520 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4521 {
4522 // return immediately if there is nothing to do.
4523 if (exchanger.ghosts_were_set == true)
4524 return;
4525
4526 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4527 reset_ghost_values(*vec[comp], exchanger);
4528 }
4529
4530
4531
4532 //
4533 // zero_vector_region
4534 //
4535
4536 // for block vectors
4537 template <int dim,
4538 typename VectorStruct,
4539 typename Number,
4540 typename VectorizedArrayType,
4541 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4542 * = nullptr>
4543 inline void
4544 zero_vector_region(
4545 const unsigned int range_index,
4546 VectorStruct &vec,
4547 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4548 {
4549 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4550 exchanger.zero_vector_region(range_index, vec.block(i));
4551 }
4552
4553
4554
4555 // for non-block vectors
4556 template <int dim,
4557 typename VectorStruct,
4558 typename Number,
4559 typename VectorizedArrayType,
4560 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4561 * = nullptr>
4562 inline void
4563 zero_vector_region(
4564 const unsigned int range_index,
4565 VectorStruct &vec,
4566 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4567 {
4568 exchanger.zero_vector_region(range_index, vec);
4569 }
4570
4571
4572
4573 // for std::vector of vectors
4574 template <int dim,
4575 typename VectorStruct,
4576 typename Number,
4577 typename VectorizedArrayType>
4578 inline void
4579 zero_vector_region(
4580 const unsigned int range_index,
4581 std::vector<VectorStruct> &vec,
4582 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4583 {
4584 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4585 zero_vector_region(range_index, vec[comp], exchanger);
4586 }
4587
4588
4589
4590 // for std::vector of pointers to vectors
4591 template <int dim,
4592 typename VectorStruct,
4593 typename Number,
4594 typename VectorizedArrayType>
4595 inline void
4596 zero_vector_region(
4597 const unsigned int range_index,
4598 std::vector<VectorStruct *> &vec,
4599 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4600 {
4601 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4602 zero_vector_region(range_index, *vec[comp], exchanger);
4603 }
4604
4605
4606
4607 // Apply a unit matrix operation to constrained DoFs: Default cases where we
4608 // cannot detect a LinearAlgebra::distributed::Vector, we do not do
4609 // anything, else we apply the constraints as a unit operation
4610 template <typename VectorStruct1, typename VectorStruct2>
4611 inline void
4612 apply_operation_to_constrained_dofs(const std::vector<unsigned int> &,
4613 const VectorStruct1 &,
4614 VectorStruct2 &)
4615 {}
4616
4617 template <typename Number>
4618 inline void
4619 apply_operation_to_constrained_dofs(
4620 const std::vector<unsigned int> &constrained_dofs,
4623 {
4624 for (const unsigned int i : constrained_dofs)
4625 dst.local_element(i) = src.local_element(i);
4626 }
4627
4628
4629 namespace MatrixFreeFunctions
4630 {
4631 // struct to select between a const interface and a non-const interface
4632 // for MFWorker
4633 template <typename, typename, typename, typename, bool>
4634 struct InterfaceSelector
4635 {};
4636
4637 // Version for constant functions
4638 template <typename MF,
4639 typename InVector,
4640 typename OutVector,
4641 typename Container>
4642 struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4643 {
4644 using function_type = void (Container::*)(
4645 const MF &,
4646 OutVector &,
4647 const InVector &,
4648 const std::pair<unsigned int, unsigned int> &) const;
4649 };
4650
4651 // Version for non-constant functions
4652 template <typename MF,
4653 typename InVector,
4654 typename OutVector,
4655 typename Container>
4656 struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4657 {
4658 using function_type =
4659 void (Container::*)(const MF &,
4660 OutVector &,
4661 const InVector &,
4662 const std::pair<unsigned int, unsigned int> &);
4663 };
4664 } // namespace MatrixFreeFunctions
4665
4666
4667
4668 // A implementation class for the worker object that runs the various
4669 // operations we want to perform during the matrix-free loop
4670 template <typename MF,
4671 typename InVector,
4672 typename OutVector,
4673 typename Container,
4674 bool is_constant>
4675 class MFWorker : public MFWorkerInterface
4676 {
4677 public:
4678 // An alias to make the arguments further down more readable
4679 using function_type = typename MatrixFreeFunctions::
4680 InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4681 function_type;
4682
4683 // constructor, binds all the arguments to this class
4684 MFWorker(const MF &matrix_free,
4685 const InVector &src,
4686 OutVector &dst,
4687 const bool zero_dst_vector_setting,
4688 const Container &container,
4689 function_type cell_function,
4690 function_type face_function,
4691 function_type boundary_function,
4692 const typename MF::DataAccessOnFaces src_vector_face_access =
4693 MF::DataAccessOnFaces::none,
4694 const typename MF::DataAccessOnFaces dst_vector_face_access =
4695 MF::DataAccessOnFaces::none,
4696 const std::function<void(const unsigned int, const unsigned int)>
4697 &operation_before_loop = {},
4698 const std::function<void(const unsigned int, const unsigned int)>
4699 &operation_after_loop = {},
4700 const unsigned int dof_handler_index_pre_post = 0)
4701 : matrix_free(matrix_free)
4702 , container(const_cast<Container &>(container))
4703 , cell_function(cell_function)
4704 , face_function(face_function)
4705 , boundary_function(boundary_function)
4706 , src(src)
4707 , dst(dst)
4708 , src_data_exchanger(matrix_free,
4709 src_vector_face_access,
4710 n_components(src))
4711 , dst_data_exchanger(matrix_free,
4712 dst_vector_face_access,
4713 n_components(dst))
4714 , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4715 , zero_dst_vector_setting(zero_dst_vector_setting &&
4716 !src_and_dst_are_same)
4717 , operation_before_loop(operation_before_loop)
4718 , operation_after_loop(operation_after_loop)
4719 , dof_handler_index_pre_post(dof_handler_index_pre_post)
4720 {
4721 Assert(!has_ghost_elements(dst),
4722 ExcMessage("The destination vector passed to the matrix-free "
4723 "loop is ghosted. This is not allowed."));
4724 }
4725
4726 // Runs the cell work. If no function is given, nothing is done
4727 virtual void
4728 cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4729 {
4730 if (cell_function != nullptr && cell_range.second > cell_range.first)
4731 for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4732 {
4733 const auto cell_subrange =
4734 matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4735
4736 if (cell_subrange.second <= cell_subrange.first)
4737 continue;
4738
4739 (container.*
4740 cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4741 }
4742 }
4743
4744 virtual void
4745 cell(const unsigned int range_index) override
4746 {
4747 process_range(cell_function,
4748 matrix_free.get_task_info().cell_partition_data_hp_ptr,
4749 matrix_free.get_task_info().cell_partition_data_hp,
4750 range_index);
4751 }
4752
4753 virtual void
4754 face(const unsigned int range_index) override
4755 {
4756 process_range(face_function,
4757 matrix_free.get_task_info().face_partition_data_hp_ptr,
4758 matrix_free.get_task_info().face_partition_data_hp,
4759 range_index);
4760 }
4761
4762 virtual void
4763 boundary(const unsigned int range_index) override
4764 {
4765 process_range(boundary_function,
4766 matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4767 matrix_free.get_task_info().boundary_partition_data_hp,
4768 range_index);
4769 }
4770
4771 private:
4772 void
4773 process_range(const function_type &fu,
4774 const std::vector<unsigned int> &ptr,
4775 const std::vector<unsigned int> &data,
4776 const unsigned int range_index)
4777 {
4778 if (fu == nullptr)
4779 return;
4780
4781 AssertIndexRange(range_index + 1, ptr.size());
4782 for (unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4783 {
4784 AssertIndexRange(2 * i + 1, data.size());
4785 (container.*fu)(matrix_free,
4786 this->dst,
4787 this->src,
4788 std::make_pair(data[2 * i], data[2 * i + 1]));
4789 }
4790 }
4791
4792 public:
4793 // Starts the communication for the update ghost values operation. We
4794 // cannot call this update if ghost and destination are the same because
4795 // that would introduce spurious entries in the destination (there is also
4796 // the problem that reading from a vector that we also write to is usually
4797 // not intended in case there is overlap, but this is up to the
4798 // application code to decide and we cannot catch this case here).
4799 virtual void
4800 vector_update_ghosts_start() override
4801 {
4802 if (!src_and_dst_are_same)
4803 internal::update_ghost_values_start(src, src_data_exchanger);
4804 }
4805
4806 // Finishes the communication for the update ghost values operation
4807 virtual void
4808 vector_update_ghosts_finish() override
4809 {
4810 if (!src_and_dst_are_same)
4811 internal::update_ghost_values_finish(src, src_data_exchanger);
4812 }
4813
4814 // Starts the communication for the vector compress operation
4815 virtual void
4816 vector_compress_start() override
4817 {
4818 internal::compress_start(dst, dst_data_exchanger);
4819 }
4820
4821 // Finishes the communication for the vector compress operation
4822 virtual void
4823 vector_compress_finish() override
4824 {
4825 internal::compress_finish(dst, dst_data_exchanger);
4826 if (!src_and_dst_are_same)
4827 internal::reset_ghost_values(src, src_data_exchanger);
4828 }
4829
4830 // Zeros the given input vector
4831 virtual void
4832 zero_dst_vector_range(const unsigned int range_index) override
4833 {
4834 if (zero_dst_vector_setting)
4835 internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4836 }
4837
4838 virtual void
4839 cell_loop_pre_range(const unsigned int range_index) override
4840 {
4841 if (operation_before_loop)
4842 {
4844 matrix_free.get_dof_info(dof_handler_index_pre_post);
4845 if (range_index == numbers::invalid_unsigned_int)
4846 {
4847 // Case with threaded loop -> currently no overlap implemented
4849 0U,
4850 dof_info.vector_partitioner->locally_owned_size(),
4851 operation_before_loop,
4853 }
4854 else
4855 {
4856 AssertIndexRange(range_index,
4857 dof_info.cell_loop_pre_list_index.size() - 1);
4858 for (unsigned int id =
4859 dof_info.cell_loop_pre_list_index[range_index];
4860 id != dof_info.cell_loop_pre_list_index[range_index + 1];
4861 ++id)
4862 operation_before_loop(dof_info.cell_loop_pre_list[id].first,
4863 dof_info.cell_loop_pre_list[id].second);
4864 }
4865 }
4866 }
4867
4868 virtual void
4869 cell_loop_post_range(const unsigned int range_index) override
4870 {
4871 if (operation_after_loop)
4872 {
4873 // Run unit matrix operation on constrained dofs if we are at the
4874 // last range
4875 const std::vector<unsigned int> &partition_row_index =
4876 matrix_free.get_task_info().partition_row_index;
4877 if (range_index ==
4878 partition_row_index[partition_row_index.size() - 2] - 1)
4879 apply_operation_to_constrained_dofs(
4880 matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
4881 src,
4882 dst);
4883
4885 matrix_free.get_dof_info(dof_handler_index_pre_post);
4886 if (range_index == numbers::invalid_unsigned_int)
4887 {
4888 // Case with threaded loop -> currently no overlap implemented
4890 0U,
4891 dof_info.vector_partitioner->locally_owned_size(),
4892 operation_after_loop,
4894 }
4895 else
4896 {
4897 AssertIndexRange(range_index,
4898 dof_info.cell_loop_post_list_index.size() - 1);
4899 for (unsigned int id =
4900 dof_info.cell_loop_post_list_index[range_index];
4901 id != dof_info.cell_loop_post_list_index[range_index + 1];
4902 ++id)
4903 operation_after_loop(dof_info.cell_loop_post_list[id].first,
4904 dof_info.cell_loop_post_list[id].second);
4905 }
4906 }
4907 }
4908
4909 private:
4910 const MF &matrix_free;
4911 Container &container;
4912 function_type cell_function;
4913 function_type face_function;
4914 function_type boundary_function;
4915
4916 const InVector &src;
4917 OutVector &dst;
4918 VectorDataExchange<MF::dimension,
4919 typename MF::value_type,
4920 typename MF::vectorized_value_type>
4921 src_data_exchanger;
4922 VectorDataExchange<MF::dimension,
4923 typename MF::value_type,
4924 typename MF::vectorized_value_type>
4925 dst_data_exchanger;
4926 const bool src_and_dst_are_same;
4927 const bool zero_dst_vector_setting;
4928 const std::function<void(const unsigned int, const unsigned int)>
4929 operation_before_loop;
4930 const std::function<void(const unsigned int, const unsigned int)>
4931 operation_after_loop;
4932 const unsigned int dof_handler_index_pre_post;
4933 };
4934
4935
4936
4941 template <class MF, typename InVector, typename OutVector>
4942 struct MFClassWrapper
4943 {
4944 using function_type =
4945 std::function<void(const MF &,
4946 OutVector &,
4947 const InVector &,
4948 const std::pair<unsigned int, unsigned int> &)>;
4949
4950 MFClassWrapper(const function_type cell,
4951 const function_type face,
4952 const function_type boundary)
4953 : cell(cell)
4954 , face(face)
4955 , boundary(boundary)
4956 {}
4957
4958 void
4959 cell_integrator(const MF &mf,
4960 OutVector &dst,
4961 const InVector &src,
4962 const std::pair<unsigned int, unsigned int> &range) const
4963 {
4964 if (cell)
4965 cell(mf, dst, src, range);
4966 }
4967
4968 void
4969 face_integrator(const MF &mf,
4970 OutVector &dst,
4971 const InVector &src,
4972 const std::pair<unsigned int, unsigned int> &range) const
4973 {
4974 if (face)
4975 face(mf, dst, src, range);
4976 }
4977
4978 void
4979 boundary_integrator(
4980 const MF &mf,
4981 OutVector &dst,
4982 const InVector &src,
4983 const std::pair<unsigned int, unsigned int> &range) const
4984 {
4985 if (boundary)
4986 boundary(mf, dst, src, range);
4987 }
4988
4989 const function_type cell;
4990 const function_type face;
4991 const function_type boundary;
4992 };
4993
4994} // end of namespace internal
4995
4996
4997
4998template <int dim, typename Number, typename VectorizedArrayType>
4999template <typename OutVector, typename InVector>
5000inline void
5002 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5003 OutVector &,
5004 const InVector &,
5005 const std::pair<unsigned int, unsigned int> &)>
5006 &cell_operation,
5007 OutVector &dst,
5008 const InVector &src,
5009 const bool zero_dst_vector) const
5010{
5011 using Wrapper =
5012 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5013 InVector,
5014 OutVector>;
5015 Wrapper wrap(cell_operation, nullptr, nullptr);
5016 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5017 InVector,
5018 OutVector,
5019 Wrapper,
5020 true>
5021 worker(*this,
5022 src,
5023 dst,
5024 zero_dst_vector,
5025 wrap,
5026 &Wrapper::cell_integrator,
5027 &Wrapper::face_integrator,
5028 &Wrapper::boundary_integrator);
5029
5030 task_info.loop(worker);
5031}
5032
5033
5034
5035template <int dim, typename Number, typename VectorizedArrayType>
5036template <typename OutVector, typename InVector>
5037inline void
5039 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5040 OutVector &,
5041 const InVector &,
5042 const std::pair<unsigned int, unsigned int> &)>
5043 &cell_operation,
5044 OutVector &dst,
5045 const InVector &src,
5046 const std::function<void(const unsigned int, const unsigned int)>
5047 &operation_before_loop,
5048 const std::function<void(const unsigned int, const unsigned int)>
5049 &operation_after_loop,
5050 const unsigned int dof_handler_index_pre_post) const
5051{
5052 using Wrapper =
5053 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5054 InVector,
5055 OutVector>;
5056 Wrapper wrap(cell_operation, nullptr, nullptr);
5057 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5058 InVector,
5059 OutVector,
5060 Wrapper,
5061 true>
5062 worker(*this,
5063 src,
5064 dst,
5065 false,
5066 wrap,
5067 &Wrapper::cell_integrator,
5068 &Wrapper::face_integrator,
5069 &Wrapper::boundary_integrator,
5070 DataAccessOnFaces::none,
5071 DataAccessOnFaces::none,
5072 operation_before_loop,
5073 operation_after_loop,
5074 dof_handler_index_pre_post);
5075
5076 task_info.loop(worker);
5077}
5078
5079
5080
5081template <int dim, typename Number, typename VectorizedArrayType>
5082template <typename OutVector, typename InVector>
5083inline void
5085 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5086 OutVector &,
5087 const InVector &,
5088 const std::pair<unsigned int, unsigned int> &)>
5089 &cell_operation,
5090 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5091 OutVector &,
5092 const InVector &,
5093 const std::pair<unsigned int, unsigned int> &)>
5094 &inner_face_operation,
5095 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5096 OutVector &,
5097 const InVector &,
5098 const std::pair<unsigned int, unsigned int> &)>
5099 &boundary_face_operation,
5100 OutVector &dst,
5101 const InVector &src,
5102 const bool zero_dst_vector,
5103 const DataAccessOnFaces dst_vector_face_access,
5104 const DataAccessOnFaces src_vector_face_access) const
5105{
5106 using Wrapper =
5107 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5108 InVector,
5109 OutVector>;
5110 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5111 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5112 InVector,
5113 OutVector,
5114 Wrapper,
5115 true>
5116 worker(*this,
5117 src,
5118 dst,
5119 zero_dst_vector,
5120 wrap,
5121 &Wrapper::cell_integrator,
5122 &Wrapper::face_integrator,
5123 &Wrapper::boundary_integrator,
5124 src_vector_face_access,
5125 dst_vector_face_access);
5126
5127 task_info.loop(worker);
5128}
5129
5130
5131
5132template <int dim, typename Number, typename VectorizedArrayType>
5133template <typename CLASS, typename OutVector, typename InVector>
5134inline void
5136 void (CLASS::*function_pointer)(
5137 const MatrixFree<dim, Number, VectorizedArrayType> &,
5138 OutVector &,
5139 const InVector &,
5140 const std::pair<unsigned int, unsigned int> &) const,
5141 const CLASS *owning_class,
5142 OutVector &dst,
5143 const InVector &src,
5144 const bool zero_dst_vector) const
5145{
5146 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5147 InVector,
5148 OutVector,
5149 CLASS,
5150 true>
5151 worker(*this,
5152 src,
5153 dst,
5154 zero_dst_vector,
5155 *owning_class,
5156 function_pointer,
5157 nullptr,
5158 nullptr);
5159 task_info.loop(worker);
5160}
5161
5162
5163
5164template <int dim, typename Number, typename VectorizedArrayType>
5165template <typename CLASS, typename OutVector, typename InVector>
5166inline void
5168 void (CLASS::*function_pointer)(
5169 const MatrixFree<dim, Number, VectorizedArrayType> &,
5170 OutVector &,
5171 const InVector &,
5172 const std::pair<unsigned int, unsigned int> &) const,
5173 const CLASS *owning_class,
5174 OutVector &dst,
5175 const InVector &src,
5176 const std::function<void(const unsigned int, const unsigned int)>
5177 &operation_before_loop,
5178 const std::function<void(const unsigned int, const unsigned int)>
5179 &operation_after_loop,
5180 const unsigned int dof_handler_index_pre_post) const
5181{
5182 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5183 InVector,
5184 OutVector,
5185 CLASS,
5186 true>
5187 worker(*this,
5188 src,
5189 dst,
5190 false,
5191 *owning_class,
5192 function_pointer,
5193 nullptr,
5194 nullptr,
5197 operation_before_loop,
5198 operation_after_loop,
5199 dof_handler_index_pre_post);
5200 task_info.loop(worker);
5201}
5202
5203
5204
5205template <int dim, typename Number, typename VectorizedArrayType>
5206template <typename CLASS, typename OutVector, typename InVector>
5207inline void
5209 void (CLASS::*cell_operation)(
5210 const MatrixFree<dim, Number, VectorizedArrayType> &,
5211 OutVector &,
5212 const InVector &,
5213 const std::pair<unsigned int, unsigned int> &) const,
5214 void (CLASS::*inner_face_operation)(
5215 const MatrixFree<dim, Number, VectorizedArrayType> &,
5216 OutVector &,
5217 const InVector &,
5218 const std::pair<unsigned int, unsigned int> &) const,
5219 void (CLASS::*boundary_face_operation)(
5220 const MatrixFree<dim, Number, VectorizedArrayType> &,
5221 OutVector &,
5222 const InVector &,
5223 const std::pair<unsigned int, unsigned int> &) const,
5224 const CLASS *owning_class,
5225 OutVector &dst,
5226 const InVector &src,
5227 const bool zero_dst_vector,
5228 const DataAccessOnFaces dst_vector_face_access,
5229 const DataAccessOnFaces src_vector_face_access) const
5230{
5231 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5232 InVector,
5233 OutVector,
5234 CLASS,
5235 true>
5236 worker(*this,
5237 src,
5238 dst,
5239 zero_dst_vector,
5240 *owning_class,
5241 cell_operation,
5242 inner_face_operation,
5243 boundary_face_operation,
5244 src_vector_face_access,
5245 dst_vector_face_access);
5246 task_info.loop(worker);
5247}
5248
5249
5250
5251template <int dim, typename Number, typename VectorizedArrayType>
5252template <typename CLASS, typename OutVector, typename InVector>
5253inline void
5255 void (CLASS::*function_pointer)(
5256 const MatrixFree<dim, Number, VectorizedArrayType> &,
5257 OutVector &,
5258 const InVector &,
5259 const std::pair<unsigned int, unsigned int> &),
5260 CLASS *owning_class,
5261 OutVector &dst,
5262 const InVector &src,
5263 const bool zero_dst_vector) const
5264{
5265 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5266 InVector,
5267 OutVector,
5268 CLASS,
5269 false>
5270 worker(*this,
5271 src,
5272 dst,
5273 zero_dst_vector,
5274 *owning_class,
5275 function_pointer,
5276 nullptr,
5277 nullptr);
5278 task_info.loop(worker);
5279}
5280
5281
5282
5283template <int dim, typename Number, typename VectorizedArrayType>
5284template <typename CLASS, typename OutVector, typename InVector>
5285inline void
5287 void (CLASS::*function_pointer)(
5288 const MatrixFree<dim, Number, VectorizedArrayType> &,
5289 OutVector &,
5290 const InVector &,
5291 const std::pair<unsigned int, unsigned int> &),
5292 CLASS *owning_class,
5293 OutVector &dst,
5294 const InVector &src,
5295 const std::function<void(const unsigned int, const unsigned int)>
5296 &operation_before_loop,
5297 const std::function<void(const unsigned int, const unsigned int)>
5298 &operation_after_loop,
5299 const unsigned int dof_handler_index_pre_post) const
5300{
5301 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5302 InVector,
5303 OutVector,
5304 CLASS,
5305 false>
5306 worker(*this,
5307 src,
5308 dst,
5309 false,
5310 *owning_class,
5311 function_pointer,
5312 nullptr,
5313 nullptr,
5316 operation_before_loop,
5317 operation_after_loop,
5318 dof_handler_index_pre_post);
5319 task_info.loop(worker);
5320}
5321
5322
5323
5324template <int dim, typename Number, typename VectorizedArrayType>
5325template <typename CLASS, typename OutVector, typename InVector>
5326inline void
5328 void (CLASS::*cell_operation)(
5329 const MatrixFree<dim, Number, VectorizedArrayType> &,
5330 OutVector &,
5331 const InVector &,
5332 const std::pair<unsigned int, unsigned int> &),
5333 void (CLASS::*inner_face_operation)(
5334 const MatrixFree<dim, Number, VectorizedArrayType> &,
5335 OutVector &,
5336 const InVector &,
5337 const std::pair<unsigned int, unsigned int> &),
5338 void (CLASS::*boundary_face_operation)(
5339 const MatrixFree<dim, Number, VectorizedArrayType> &,
5340 OutVector &,
5341 const InVector &,
5342 const std::pair<unsigned int, unsigned int> &),
5343 CLASS *owning_class,
5344 OutVector &dst,
5345 const InVector &src,
5346 const bool zero_dst_vector,
5347 const DataAccessOnFaces dst_vector_face_access,
5348 const DataAccessOnFaces src_vector_face_access) const
5349{
5350 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5351 InVector,
5352 OutVector,
5353 CLASS,
5354 false>
5355 worker(*this,
5356 src,
5357 dst,
5358 zero_dst_vector,
5359 *owning_class,
5360 cell_operation,
5361 inner_face_operation,
5362 boundary_face_operation,
5363 src_vector_face_access,
5364 dst_vector_face_access);
5365 task_info.loop(worker);
5366}
5367
5368
5369
5370template <int dim, typename Number, typename VectorizedArrayType>
5371template <typename OutVector, typename InVector>
5372inline void
5374 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5375 OutVector &,
5376 const InVector &,
5377 const std::pair<unsigned int, unsigned int> &)>
5378 &cell_operation,
5379 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5380 OutVector &,
5381 const InVector &,
5382 const std::pair<unsigned int, unsigned int> &)>
5383 &inner_face_operation,
5384 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5385 OutVector &,
5386 const InVector &,
5387 const std::pair<unsigned int, unsigned int> &)>
5388 &boundary_face_operation,
5389 OutVector &dst,
5390 const InVector &src,
5391 const std::function<void(const unsigned int, const unsigned int)>
5392 &operation_before_loop,
5393 const std::function<void(const unsigned int, const unsigned int)>
5394 &operation_after_loop,
5395 const unsigned int dof_handler_index_pre_post,
5396 const DataAccessOnFaces dst_vector_face_access,
5397 const DataAccessOnFaces src_vector_face_access) const
5398{
5399 using Wrapper =
5400 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5401 InVector,
5402 OutVector>;
5403 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5404 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5405 InVector,
5406 OutVector,
5407 Wrapper,
5408 true>
5409 worker(*this,
5410 src,
5411 dst,
5412 false,
5413 wrap,
5414 &Wrapper::cell_integrator,
5415 &Wrapper::face_integrator,
5416 &Wrapper::boundary_integrator,
5417 src_vector_face_access,
5418 dst_vector_face_access,
5419 operation_before_loop,
5420 operation_after_loop,
5421 dof_handler_index_pre_post);
5422
5423 task_info.loop(worker);
5424}
5425
5426
5427
5428template <int dim, typename Number, typename VectorizedArrayType>
5429template <typename CLASS, typename OutVector, typename InVector>
5430inline void
5432 void (CLASS::*cell_operation)(const MatrixFree &,
5433 OutVector &,
5434 const InVector &,
5435 const std::pair<unsigned int, unsigned int> &)
5436 const,
5437 void (CLASS::*inner_face_operation)(
5438 const MatrixFree &,
5439 OutVector &,
5440 const InVector &,
5441 const std::pair<unsigned int, unsigned int> &) const,
5442 void (CLASS::*boundary_face_operation)(
5443 const MatrixFree &,
5444 OutVector &,
5445 const InVector &,
5446 const std::pair<unsigned int, unsigned int> &) const,
5447 const CLASS *owning_class,
5448 OutVector &dst,
5449 const InVector &src,
5450 const std::function<void(const unsigned int, const unsigned int)>
5451 &operation_before_loop,
5452 const std::function<void(const unsigned int, const unsigned int)>
5453 &operation_after_loop,
5454 const unsigned int dof_handler_index_pre_post,
5455 const DataAccessOnFaces dst_vector_face_access,
5456 const DataAccessOnFaces src_vector_face_access) const
5457{
5458 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5459 InVector,
5460 OutVector,
5461 CLASS,
5462 true>
5463 worker(*this,
5464 src,
5465 dst,
5466 false,
5467 *owning_class,
5468 cell_operation,
5469 inner_face_operation,
5470 boundary_face_operation,
5471 src_vector_face_access,
5472 dst_vector_face_access,
5473 operation_before_loop,
5474 operation_after_loop,
5475 dof_handler_index_pre_post);
5476 task_info.loop(worker);
5477}
5478
5479
5480
5481template <int dim, typename Number, typename VectorizedArrayType>
5482template <typename CLASS, typename OutVector, typename InVector>
5483inline void
5485 void (CLASS::*cell_operation)(const MatrixFree &,
5486 OutVector &,
5487 const InVector &,
5488 const std::pair<unsigned int, unsigned int> &),
5489 void (CLASS::*inner_face_operation)(
5490 const MatrixFree &,
5491 OutVector &,
5492 const InVector &,
5493 const std::pair<unsigned int, unsigned int> &),
5494 void (CLASS::*boundary_face_operation)(
5495 const MatrixFree &,
5496 OutVector &,
5497 const InVector &,
5498 const std::pair<unsigned int, unsigned int> &),
5499 const CLASS *owning_class,
5500 OutVector &dst,
5501 const InVector &src,
5502 const std::function<void(const unsigned int, const unsigned int)>
5503 &operation_before_loop,
5504 const std::function<void(const unsigned int, const unsigned int)>
5505 &operation_after_loop,
5506 const unsigned int dof_handler_index_pre_post,
5507 const DataAccessOnFaces dst_vector_face_access,
5508 const DataAccessOnFaces src_vector_face_access) const
5509{
5510 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5511 InVector,
5512 OutVector,
5513 CLASS,
5514 false>
5515 worker(*this,
5516 src,
5517 dst,
5518 false,
5519 *owning_class,
5520 cell_operation,
5521 inner_face_operation,
5522 boundary_face_operation,
5523 src_vector_face_access,
5524 dst_vector_face_access,
5525 operation_before_loop,
5526 operation_after_loop,
5527 dof_handler_index_pre_post);
5528 task_info.loop(worker);
5529}
5530
5531
5532
5533template <int dim, typename Number, typename VectorizedArrayType>
5534template <typename CLASS, typename OutVector, typename InVector>
5535inline void
5537 void (CLASS::*function_pointer)(
5538 const MatrixFree<dim, Number, VectorizedArrayType> &,
5539 OutVector &,
5540 const InVector &,
5541 const std::pair<unsigned int, unsigned int> &) const,
5542 const CLASS *owning_class,
5543 OutVector &dst,
5544 const InVector &src,
5545 const bool zero_dst_vector,
5546 const DataAccessOnFaces src_vector_face_access) const
5547{
5548 auto src_vector_face_access_temp = src_vector_face_access;
5549 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5550 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5551 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5552 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5553
5554 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5555 InVector,
5556 OutVector,
5557 CLASS,
5558 true>
5559 worker(*this,
5560 src,
5561 dst,
5562 zero_dst_vector,
5563 *owning_class,
5564 function_pointer,
5565 nullptr,
5566 nullptr,
5567 src_vector_face_access_temp,
5569 task_info.loop(worker);
5570}
5571
5572
5573
5574template <int dim, typename Number, typename VectorizedArrayType>
5575template <typename CLASS, typename OutVector, typename InVector>
5576inline void
5578 void (CLASS::*function_pointer)(
5579 const MatrixFree<dim, Number, VectorizedArrayType> &,
5580 OutVector &,
5581 const InVector &,
5582 const std::pair<unsigned int, unsigned int> &),
5583 CLASS *owning_class,
5584 OutVector &dst,
5585 const InVector &src,
5586 const bool zero_dst_vector,
5587 const DataAccessOnFaces src_vector_face_access) const
5588{
5589 auto src_vector_face_access_temp = src_vector_face_access;
5590 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5591 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5592 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5593 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5594
5595 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5596 InVector,
5597 OutVector,
5598 CLASS,
5599 false>
5600 worker(*this,
5601 src,
5602 dst,
5603 zero_dst_vector,
5604 *owning_class,
5605 function_pointer,
5606 nullptr,
5607 nullptr,
5608 src_vector_face_access_temp,
5610 task_info.loop(worker);
5611}
5612
5613
5614
5615template <int dim, typename Number, typename VectorizedArrayType>
5616template <typename OutVector, typename InVector>
5617inline void
5619 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5620 OutVector &,
5621 const InVector &,
5622 const std::pair<unsigned int, unsigned int> &)>
5623 &cell_operation,
5624 OutVector &dst,
5625 const InVector &src,
5626 const bool zero_dst_vector,
5627 const DataAccessOnFaces src_vector_face_access) const
5628{
5629 auto src_vector_face_access_temp = src_vector_face_access;
5630 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5631 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5632 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5633 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5634
5635 using Wrapper =
5636 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5637 InVector,
5638 OutVector>;
5639 Wrapper wrap(cell_operation, nullptr, nullptr);
5640
5641 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5642 InVector,
5643 OutVector,
5644 Wrapper,
5645 true>
5646 worker(*this,
5647 src,
5648 dst,
5649 zero_dst_vector,
5650 wrap,
5651 &Wrapper::cell_integrator,
5652 &Wrapper::face_integrator,
5653 &Wrapper::boundary_integrator,
5654 src_vector_face_access_temp,
5655 DataAccessOnFaces::none);
5656 task_info.loop(worker);
5657}
5658
5659
5660#endif // ifndef DOXYGEN
5661
5662
5663
5665
5666#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
unsigned int get_cell_active_fe_index(const std::pair< unsigned int, unsigned int > range) const
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
unsigned int n_ghost_cell_batches() const
void initialize_dof_handlers(const std::vector< const DoFHandler< dim, dim > * > &dof_handlers, const AdditionalData &additional_data)
types::boundary_id get_boundary_id(const unsigned int face_batch_index) const
void print(std::ostream &out) const
void update_mapping(const std::shared_ptr< hp::MappingCollection< dim > > &mapping)
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
unsigned int n_inner_face_batches() const
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
void update_mapping(const Mapping< dim > &mapping)
const Quadrature< dim > & get_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_index) const
unsigned int get_mg_level() const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*boundary_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void print_memory_consumption(StreamType &out) const
bool mapping_initialized() const
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
void clear()
void internal_reinit(const std::shared_ptr< hp::MappingCollection< dim > > &mapping, const std::vector< const DoFHandler< dim, dim > * > &dof_handlers, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const std::vector< hp::QCollection< q_dim > > &quad, const AdditionalData &additional_data)
unsigned int get_n_q_points_face(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &inner_face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_face_operation, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
bool at_irregular_cell(const unsigned int cell_batch_index) const
~MatrixFree() override=default
const AffineConstraints< Number > & get_affine_constraints(const unsigned int dof_handler_index=0) const
void copy_from(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free_base)
void initialize_face_data_vector(AlignedVector< T > &vec) const
unsigned int get_dofs_per_cell(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int n_constraint_pool_entries() const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
const Number * constraint_pool_begin(const unsigned int pool_index) const
AlignedVector< Number > * acquire_scratch_data_non_threadsafe() const
const IndexSet & get_locally_owned_set(const unsigned int dof_handler_index=0) const
unsigned int get_dofs_per_face(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
internal::MatrixFreeFunctions::TaskInfo task_info
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
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
std::pair< unsigned int, unsigned int > get_face_range_category(const std::pair< unsigned int, unsigned int > face_batch_range) const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
std::vector< std::pair< unsigned int, unsigned int > > cell_level_index
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner(const unsigned int dof_handler_index=0) const
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_cell_category(const unsigned int cell_batch_index) const
Threads::ThreadLocalStorage< std::list< std::pair< bool, AlignedVector< VectorizedArrayType > > > > scratch_pad
unsigned int get_matrix_free_cell_index(const typename Triangulation< dim >::cell_iterator &cell) const
Number value_type
std::size_t memory_consumption() const
std::pair< unsigned int, unsigned int > get_face_category(const unsigned int face_batch_index) const
unsigned int n_boundary_face_batches() const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*boundary_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
std::list< std::pair< bool, AlignedVector< Number > > > scratch_pad_non_threadsafe
const Number * constraint_pool_end(const unsigned int pool_index) const
void loop_cell_centric(const std::function< void(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
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
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*boundary_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
std::vector< internal::MatrixFreeFunctions::DoFInfo > dof_info
void initialize_cell_data_vector(AlignedVector< T > &vec) const
unsigned int n_ghost_inner_face_batches() const
void release_scratch_data_non_threadsafe(const AlignedVector< Number > *memory) const
void renumber_dofs(std::vector< types::global_dof_index > &renumbering, const unsigned int dof_handler_index=0)
unsigned int get_face_active_fe_index(const std::pair< unsigned int, unsigned int > range, const bool is_interior_face=true) const
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
bool indices_are_initialized
internal::MatrixFreeFunctions::FaceInfo< VectorizedArrayType::size()> face_info
void reinit(const MappingType &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< number2 > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
unsigned int get_cell_range_category(const std::pair< unsigned int, unsigned int > cell_batch_range) const
void reinit(const MappingType &mapping, const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< QuadratureType > &quad, const AdditionalData &additional_data=AdditionalData())
std::array< types::boundary_id, VectorizedArrayType::size()> get_faces_by_cells_boundary_id(const unsigned int cell_batch_index, const unsigned int face_number) const
std::vector< unsigned int > constraint_pool_row_index
const internal::MatrixFreeFunctions::ShapeInfo< Number > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
const std::vector< unsigned int > & get_constrained_dofs(const unsigned int dof_handler_index=0) const
unsigned int cell_level_index_end_local
unsigned int n_base_elements(const unsigned int dof_handler_index) const
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< Number > > shape_info
static bool is_supported(const FiniteElement< dim, spacedim > &fe)
static constexpr unsigned int dimension
A class that provides a separate storage location on each thread that accesses the object.
#define DEAL_II_NAMESPACE_OPEN
Definition config.h: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