Reference documentation for deal.II version GIT relicensing-426-g7976cfd195 2024-04-18 21:10:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
matrix_free.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2012 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
16#ifndef dealii_matrix_free_h
17#define dealii_matrix_free_h
18
19#include <deal.II/base/config.h>
20
27
29
30#include <deal.II/fe/fe.h>
31#include <deal.II/fe/mapping.h>
32
35
40
47
48#include <cstdlib>
49#include <limits>
50#include <list>
51#include <memory>
52
53
55
56
57
109template <int dim,
110 typename Number = double,
111 typename VectorizedArrayType = VectorizedArray<Number>>
113{
114 static_assert(
115 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
116 "Type of Number and of VectorizedArrayType do not match.");
117
118public:
123 using value_type = Number;
124 using vectorized_value_type = VectorizedArrayType;
125
129 static constexpr unsigned int dimension = dim;
130
184 {
189
217
223 const unsigned int tasks_block_size = 0,
229 const unsigned int mg_level = numbers::invalid_unsigned_int,
230 const bool store_plain_indices = true,
231 const bool initialize_indices = true,
232 const bool initialize_mapping = true,
233 const bool overlap_communication_computation = true,
234 const bool hold_all_faces_to_owned_cells = false,
235 const bool cell_vectorization_categories_strict = false,
236 const bool allow_ghosted_vectors_in_loops = true)
252 , communicator_sm(MPI_COMM_SELF)
253 {}
254
280
310
348
358 unsigned int tasks_block_size;
359
374
394
414
442
451 unsigned int mg_level;
452
460
471
480
494
503
525 std::vector<unsigned int> cell_vectorization_category;
526
537
549
554 };
555
564
569
573 ~MatrixFree() override = default;
574
587 template <typename QuadratureType, typename number2, typename MappingType>
588 void
589 reinit(const MappingType &mapping,
590 const DoFHandler<dim> &dof_handler,
591 const AffineConstraints<number2> &constraint,
592 const QuadratureType &quad,
593 const AdditionalData &additional_data = AdditionalData());
594
616 template <typename QuadratureType, typename number2, typename MappingType>
617 void
618 reinit(const MappingType &mapping,
619 const std::vector<const DoFHandler<dim> *> &dof_handler,
620 const std::vector<const AffineConstraints<number2> *> &constraint,
621 const std::vector<QuadratureType> &quad,
622 const AdditionalData &additional_data = AdditionalData());
623
631 template <typename QuadratureType, typename number2, typename MappingType>
632 void
633 reinit(const MappingType &mapping,
634 const std::vector<const DoFHandler<dim> *> &dof_handler,
635 const std::vector<const AffineConstraints<number2> *> &constraint,
636 const QuadratureType &quad,
637 const AdditionalData &additional_data = AdditionalData());
638
644 void
646 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free_base);
647
657 void
659
663 void
664 update_mapping(const std::shared_ptr<hp::MappingCollection<dim>> &mapping);
665
670 void
672
687 {
694 none,
695
706 values,
707
717
728 gradients,
729
739
747 };
748
793 template <typename OutVector, typename InVector>
794 void
795 cell_loop(const std::function<void(
797 OutVector &,
798 const InVector &,
799 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
800 OutVector &dst,
801 const InVector &src,
802 const bool zero_dst_vector = false) const;
803
850 template <typename CLASS, typename OutVector, typename InVector>
851 void
852 cell_loop(void (CLASS::*cell_operation)(
853 const MatrixFree &,
854 OutVector &,
855 const InVector &,
856 const std::pair<unsigned int, unsigned int> &) const,
857 const CLASS *owning_class,
858 OutVector &dst,
859 const InVector &src,
860 const bool zero_dst_vector = false) const;
861
865 template <typename CLASS, typename OutVector, typename InVector>
866 void
867 cell_loop(void (CLASS::*cell_operation)(
868 const MatrixFree &,
869 OutVector &,
870 const InVector &,
871 const std::pair<unsigned int, unsigned int> &),
872 CLASS *owning_class,
873 OutVector &dst,
874 const InVector &src,
875 const bool zero_dst_vector = false) const;
876
961 template <typename CLASS, typename OutVector, typename InVector>
962 void
963 cell_loop(void (CLASS::*cell_operation)(
964 const MatrixFree &,
965 OutVector &,
966 const InVector &,
967 const std::pair<unsigned int, unsigned int> &) const,
968 const CLASS *owning_class,
969 OutVector &dst,
970 const InVector &src,
971 const std::function<void(const unsigned int, const unsigned int)>
972 &operation_before_loop,
973 const std::function<void(const unsigned int, const unsigned int)>
974 &operation_after_loop,
975 const unsigned int dof_handler_index_pre_post = 0) const;
976
980 template <typename CLASS, typename OutVector, typename InVector>
981 void
982 cell_loop(void (CLASS::*cell_operation)(
983 const MatrixFree &,
984 OutVector &,
985 const InVector &,
986 const std::pair<unsigned int, unsigned int> &),
987 CLASS *owning_class,
988 OutVector &dst,
989 const InVector &src,
990 const std::function<void(const unsigned int, const unsigned int)>
991 &operation_before_loop,
992 const std::function<void(const unsigned int, const unsigned int)>
993 &operation_after_loop,
994 const unsigned int dof_handler_index_pre_post = 0) const;
995
1000 template <typename OutVector, typename InVector>
1001 void
1002 cell_loop(const std::function<void(
1004 OutVector &,
1005 const InVector &,
1006 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1007 OutVector &dst,
1008 const InVector &src,
1009 const std::function<void(const unsigned int, const unsigned int)>
1010 &operation_before_loop,
1011 const std::function<void(const unsigned int, const unsigned int)>
1012 &operation_after_loop,
1013 const unsigned int dof_handler_index_pre_post = 0) const;
1014
1090 template <typename OutVector, typename InVector>
1091 void
1092 loop(const std::function<
1094 OutVector &,
1095 const InVector &,
1096 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1097 const std::function<
1099 OutVector &,
1100 const InVector &,
1101 const std::pair<unsigned int, unsigned int> &)> &face_operation,
1102 const std::function<void(
1104 OutVector &,
1105 const InVector &,
1106 const std::pair<unsigned int, unsigned int> &)> &boundary_operation,
1107 OutVector &dst,
1108 const InVector &src,
1109 const bool zero_dst_vector = false,
1110 const DataAccessOnFaces dst_vector_face_access =
1112 const DataAccessOnFaces src_vector_face_access =
1114
1200 template <typename CLASS, typename OutVector, typename InVector>
1201 void
1203 void (CLASS::*cell_operation)(const MatrixFree &,
1204 OutVector &,
1205 const InVector &,
1206 const std::pair<unsigned int, unsigned int> &)
1207 const,
1208 void (CLASS::*face_operation)(const MatrixFree &,
1209 OutVector &,
1210 const InVector &,
1211 const std::pair<unsigned int, unsigned int> &)
1212 const,
1213 void (CLASS::*boundary_operation)(
1214 const MatrixFree &,
1215 OutVector &,
1216 const InVector &,
1217 const std::pair<unsigned int, unsigned int> &) const,
1218 const CLASS *owning_class,
1219 OutVector &dst,
1220 const InVector &src,
1221 const bool zero_dst_vector = false,
1222 const DataAccessOnFaces dst_vector_face_access =
1224 const DataAccessOnFaces src_vector_face_access =
1226
1230 template <typename CLASS, typename OutVector, typename InVector>
1231 void
1232 loop(void (CLASS::*cell_operation)(
1233 const MatrixFree &,
1234 OutVector &,
1235 const InVector &,
1236 const std::pair<unsigned int, unsigned int> &),
1237 void (CLASS::*face_operation)(
1238 const MatrixFree &,
1239 OutVector &,
1240 const InVector &,
1241 const std::pair<unsigned int, unsigned int> &),
1242 void (CLASS::*boundary_operation)(
1243 const MatrixFree &,
1244 OutVector &,
1245 const InVector &,
1246 const std::pair<unsigned int, unsigned int> &),
1247 CLASS *owning_class,
1248 OutVector &dst,
1249 const InVector &src,
1250 const bool zero_dst_vector = false,
1251 const DataAccessOnFaces dst_vector_face_access =
1253 const DataAccessOnFaces src_vector_face_access =
1255
1375 template <typename CLASS, typename OutVector, typename InVector>
1376 void
1378 void (CLASS::*cell_operation)(const MatrixFree &,
1379 OutVector &,
1380 const InVector &,
1381 const std::pair<unsigned int, unsigned int> &)
1382 const,
1383 void (CLASS::*face_operation)(const MatrixFree &,
1384 OutVector &,
1385 const InVector &,
1386 const std::pair<unsigned int, unsigned int> &)
1387 const,
1388 void (CLASS::*boundary_operation)(
1389 const MatrixFree &,
1390 OutVector &,
1391 const InVector &,
1392 const std::pair<unsigned int, unsigned int> &) const,
1393 const CLASS *owning_class,
1394 OutVector &dst,
1395 const InVector &src,
1396 const std::function<void(const unsigned int, const unsigned int)>
1397 &operation_before_loop,
1398 const std::function<void(const unsigned int, const unsigned int)>
1399 &operation_after_loop,
1400 const unsigned int dof_handler_index_pre_post = 0,
1401 const DataAccessOnFaces dst_vector_face_access =
1403 const DataAccessOnFaces src_vector_face_access =
1405
1409 template <typename CLASS, typename OutVector, typename InVector>
1410 void
1411 loop(void (CLASS::*cell_operation)(
1412 const MatrixFree &,
1413 OutVector &,
1414 const InVector &,
1415 const std::pair<unsigned int, unsigned int> &),
1416 void (CLASS::*face_operation)(
1417 const MatrixFree &,
1418 OutVector &,
1419 const InVector &,
1420 const std::pair<unsigned int, unsigned int> &),
1421 void (CLASS::*boundary_operation)(
1422 const MatrixFree &,
1423 OutVector &,
1424 const InVector &,
1425 const std::pair<unsigned int, unsigned int> &),
1426 const CLASS *owning_class,
1427 OutVector &dst,
1428 const InVector &src,
1429 const std::function<void(const unsigned int, const unsigned int)>
1430 &operation_before_loop,
1431 const std::function<void(const unsigned int, const unsigned int)>
1432 &operation_after_loop,
1433 const unsigned int dof_handler_index_pre_post = 0,
1434 const DataAccessOnFaces dst_vector_face_access =
1436 const DataAccessOnFaces src_vector_face_access =
1438
1444 template <typename OutVector, typename InVector>
1445 void
1446 loop(const std::function<
1448 OutVector &,
1449 const InVector &,
1450 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1451 const std::function<
1453 OutVector &,
1454 const InVector &,
1455 const std::pair<unsigned int, unsigned int> &)> &face_operation,
1456 const std::function<void(
1458 OutVector &,
1459 const InVector &,
1460 const std::pair<unsigned int, unsigned int> &)> &boundary_operation,
1461 OutVector &dst,
1462 const InVector &src,
1463 const std::function<void(const unsigned int, const unsigned int)>
1464 &operation_before_loop,
1465 const std::function<void(const unsigned int, const unsigned int)>
1466 &operation_after_loop,
1467 const unsigned int dof_handler_index_pre_post = 0,
1468 const DataAccessOnFaces dst_vector_face_access =
1470 const DataAccessOnFaces src_vector_face_access =
1472
1537 template <typename CLASS, typename OutVector, typename InVector>
1538 void
1539 loop_cell_centric(void (CLASS::*cell_operation)(
1540 const MatrixFree &,
1541 OutVector &,
1542 const InVector &,
1543 const std::pair<unsigned int, unsigned int> &) const,
1544 const CLASS *owning_class,
1545 OutVector &dst,
1546 const InVector &src,
1547 const bool zero_dst_vector = false,
1548 const DataAccessOnFaces src_vector_face_access =
1550
1554 template <typename CLASS, typename OutVector, typename InVector>
1555 void
1556 loop_cell_centric(void (CLASS::*cell_operation)(
1557 const MatrixFree &,
1558 OutVector &,
1559 const InVector &,
1560 const std::pair<unsigned int, unsigned int> &),
1561 CLASS *owning_class,
1562 OutVector &dst,
1563 const InVector &src,
1564 const bool zero_dst_vector = false,
1565 const DataAccessOnFaces src_vector_face_access =
1567
1571 template <typename OutVector, typename InVector>
1572 void
1574 const std::function<void(const MatrixFree &,
1575 OutVector &,
1576 const InVector &,
1577 const std::pair<unsigned int, unsigned int> &)>
1578 &cell_operation,
1579 OutVector &dst,
1580 const InVector &src,
1581 const bool zero_dst_vector = false,
1582 const DataAccessOnFaces src_vector_face_access =
1584
1592 std::pair<unsigned int, unsigned int>
1593 create_cell_subrange_hp(const std::pair<unsigned int, unsigned int> &range,
1594 const unsigned int fe_degree,
1595 const unsigned int dof_handler_index = 0) const;
1596
1603 std::pair<unsigned int, unsigned int>
1605 const std::pair<unsigned int, unsigned int> &range,
1606 const unsigned int fe_index,
1607 const unsigned int dof_handler_index = 0) const;
1608
1612 unsigned int
1614
1618 unsigned int
1620 const std::pair<unsigned int, unsigned int> range) const;
1621
1625 unsigned int
1626 get_face_active_fe_index(const std::pair<unsigned int, unsigned int> range,
1627 const bool is_interior_face = true) const;
1628
1640 template <typename T>
1641 void
1643
1649 template <typename T>
1650 void
1652
1666 template <typename VectorType>
1667 void
1668 initialize_dof_vector(VectorType &vec,
1669 const unsigned int dof_handler_index = 0) const;
1670
1688 template <typename Number2>
1689 void
1691 const unsigned int dof_handler_index = 0) const;
1692
1703 const std::shared_ptr<const Utilities::MPI::Partitioner> &
1704 get_vector_partitioner(const unsigned int dof_handler_index = 0) const;
1705
1709 const IndexSet &
1710 get_locally_owned_set(const unsigned int dof_handler_index = 0) const;
1711
1715 const IndexSet &
1716 get_ghost_set(const unsigned int dof_handler_index = 0) const;
1717
1727 const std::vector<unsigned int> &
1728 get_constrained_dofs(const unsigned int dof_handler_index = 0) const;
1729
1740 void
1741 renumber_dofs(std::vector<types::global_dof_index> &renumbering,
1742 const unsigned int dof_handler_index = 0);
1743
1753 template <int spacedim>
1754 static bool
1756
1760 unsigned int
1762
1767 unsigned int
1768 n_base_elements(const unsigned int dof_handler_index) const;
1769
1777 unsigned int
1779
1789 unsigned int
1791
1798 unsigned int
1800
1810 unsigned int
1812
1823 unsigned int
1825
1830 unsigned int
1832
1844 get_boundary_id(const unsigned int face_batch_index) const;
1845
1850 std::array<types::boundary_id, VectorizedArrayType::size()>
1851 get_faces_by_cells_boundary_id(const unsigned int cell_batch_index,
1852 const unsigned int face_number) const;
1853
1858 const DoFHandler<dim> &
1859 get_dof_handler(const unsigned int dof_handler_index = 0) const;
1860
1868 get_affine_constraints(const unsigned int dof_handler_index = 0) const;
1869
1883 get_cell_iterator(const unsigned int cell_batch_index,
1884 const unsigned int lane_index,
1885 const unsigned int dof_handler_index = 0) const;
1886
1892 std::pair<int, int>
1893 get_cell_level_and_index(const unsigned int cell_batch_index,
1894 const unsigned int lane_index) const;
1895
1902 unsigned int
1904 const typename Triangulation<dim>::cell_iterator &cell) const;
1905
1918 std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>
1919 get_face_iterator(const unsigned int face_batch_index,
1920 const unsigned int lane_index,
1921 const bool interior = true,
1922 const unsigned int fe_component = 0) const;
1923
1936 bool
1937 at_irregular_cell(const unsigned int cell_batch_index) const;
1938
1948 unsigned int
1949 n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const;
1950
1960 unsigned int
1961 n_active_entries_per_face_batch(const unsigned int face_batch_index) const;
1962
1966 unsigned int
1967 get_dofs_per_cell(const unsigned int dof_handler_index = 0,
1968 const unsigned int hp_active_fe_index = 0) const;
1969
1973 unsigned int
1974 get_n_q_points(const unsigned int quad_index = 0,
1975 const unsigned int hp_active_fe_index = 0) const;
1976
1981 unsigned int
1982 get_dofs_per_face(const unsigned int dof_handler_index = 0,
1983 const unsigned int hp_active_fe_index = 0) const;
1984
1989 unsigned int
1990 get_n_q_points_face(const unsigned int quad_index = 0,
1991 const unsigned int hp_active_fe_index = 0) const;
1992
1996 const Quadrature<dim> &
1997 get_quadrature(const unsigned int quad_index = 0,
1998 const unsigned int hp_active_fe_index = 0) const;
1999
2003 const Quadrature<dim - 1> &
2004 get_face_quadrature(const unsigned int quad_index = 0,
2005 const unsigned int hp_active_fe_index = 0) const;
2006
2019 unsigned int
2021 const std::pair<unsigned int, unsigned int> cell_batch_range) const;
2022
2027 std::pair<unsigned int, unsigned int>
2029 const std::pair<unsigned int, unsigned int> face_batch_range) const;
2030
2042 unsigned int
2043 get_cell_category(const unsigned int cell_batch_index) const;
2044
2049 std::pair<unsigned int, unsigned int>
2050 get_face_category(const unsigned int face_batch_index) const;
2051
2055 bool
2057
2062 bool
2064
2069 unsigned int
2071
2076 std::size_t
2078
2083 template <typename StreamType>
2084 void
2085 print_memory_consumption(StreamType &out) const;
2086
2091 void
2092 print(std::ostream &out) const;
2093
2107
2108 /*
2109 * Return geometry-dependent information on the cells.
2110 */
2111 const internal::MatrixFreeFunctions::
2112 MappingInfo<dim, Number, VectorizedArrayType> &
2114
2119 get_dof_info(const unsigned int dof_handler_index_component = 0) const;
2120
2124 unsigned int
2126
2131 const Number *
2132 constraint_pool_begin(const unsigned int pool_index) const;
2133
2139 const Number *
2140 constraint_pool_end(const unsigned int pool_index) const;
2141
2146 get_shape_info(const unsigned int dof_handler_index_component = 0,
2147 const unsigned int quad_index = 0,
2148 const unsigned int fe_base_element = 0,
2149 const unsigned int hp_active_fe_index = 0,
2150 const unsigned int hp_active_quad_index = 0) const;
2151
2156 VectorizedArrayType::size()> &
2157 get_face_info(const unsigned int face_batch_index) const;
2158
2159
2167
2183
2187 void
2189
2201
2205 void
2207 const AlignedVector<Number> *memory) const;
2208
2211private:
2216 template <typename number2, int q_dim>
2217 void
2219 const std::shared_ptr<hp::MappingCollection<dim>> &mapping,
2220 const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
2221 const std::vector<const AffineConstraints<number2> *> &constraint,
2222 const std::vector<IndexSet> &locally_owned_set,
2223 const std::vector<hp::QCollection<q_dim>> &quad,
2224 const AdditionalData &additional_data);
2225
2232 template <typename number2>
2233 void
2235 const std::vector<const AffineConstraints<number2> *> &constraint,
2236 const std::vector<IndexSet> &locally_owned_set,
2237 const AdditionalData &additional_data);
2238
2242 void
2244 const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
2245 const AdditionalData &additional_data);
2246
2250 std::vector<SmartPointer<const DoFHandler<dim>>> dof_handlers;
2251
2258 std::vector<SmartPointer<const AffineConstraints<Number>>> affine_constraints;
2259
2264 std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
2265
2272 std::vector<Number> constraint_pool_data;
2273
2278 std::vector<unsigned int> constraint_pool_row_index;
2279
2286
2292
2299 std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
2300
2305 std::vector<unsigned int> mf_cell_indices;
2306
2314
2321
2326 internal::MatrixFreeFunctions::FaceInfo<VectorizedArrayType::size()>
2328
2333
2338
2347 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2349
2354 mutable std::list<std::pair<bool, AlignedVector<Number>>>
2356
2360 unsigned int mg_level;
2361};
2362
2363
2364
2365/*----------------------- Inline functions ----------------------------------*/
2366
2367#ifndef DOXYGEN
2368
2369
2370
2371template <int dim, typename Number, typename VectorizedArrayType>
2372template <typename T>
2373inline void
2375 AlignedVector<T> &vec) const
2376{
2377 vec.resize(this->n_cell_batches() + this->n_ghost_cell_batches());
2378}
2379
2380
2381
2382template <int dim, typename Number, typename VectorizedArrayType>
2383template <typename T>
2384inline void
2386 AlignedVector<T> &vec) const
2387{
2388 vec.resize(this->n_inner_face_batches() + this->n_boundary_face_batches() +
2389 this->n_ghost_inner_face_batches());
2390}
2391
2392
2393
2394template <int dim, typename Number, typename VectorizedArrayType>
2395template <typename VectorType>
2396inline void
2398 VectorType &vec,
2399 const unsigned int comp) const
2400{
2401 static_assert(IsBlockVector<VectorType>::value == false,
2402 "This function is not supported for block vectors.");
2403
2404 Assert(task_info.n_procs == 1,
2405 ExcMessage("This function can only be used in serial."));
2406
2407 AssertIndexRange(comp, n_components());
2408 vec.reinit(dof_info[comp].vector_partitioner->size());
2409}
2410
2411
2412
2413template <int dim, typename Number, typename VectorizedArrayType>
2414template <typename Number2>
2415inline void
2418 const unsigned int comp) const
2419{
2420 AssertIndexRange(comp, n_components());
2421 vec.reinit(dof_info[comp].vector_partitioner, task_info.communicator_sm);
2422}
2423
2424
2425
2426template <int dim, typename Number, typename VectorizedArrayType>
2427inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2429 const unsigned int comp) const
2430{
2431 AssertIndexRange(comp, n_components());
2432 return dof_info[comp].vector_partitioner;
2433}
2434
2435
2436
2437template <int dim, typename Number, typename VectorizedArrayType>
2438inline const std::vector<unsigned int> &
2440 const unsigned int comp) const
2441{
2442 AssertIndexRange(comp, n_components());
2443 return dof_info[comp].constrained_dofs;
2444}
2445
2446
2447
2448template <int dim, typename Number, typename VectorizedArrayType>
2449inline unsigned int
2451{
2452 AssertDimension(dof_handlers.size(), dof_info.size());
2453 return dof_handlers.size();
2454}
2455
2456
2457
2458template <int dim, typename Number, typename VectorizedArrayType>
2459inline unsigned int
2461 const unsigned int dof_no) const
2462{
2463 AssertDimension(dof_handlers.size(), dof_info.size());
2464 AssertIndexRange(dof_no, dof_handlers.size());
2465 return dof_handlers[dof_no]->get_fe().n_base_elements();
2466}
2467
2468
2469
2470template <int dim, typename Number, typename VectorizedArrayType>
2473{
2474 return task_info;
2475}
2476
2477
2478
2479template <int dim, typename Number, typename VectorizedArrayType>
2480inline unsigned int
2482{
2483 return task_info.n_active_cells;
2484}
2485
2486
2487
2488template <int dim, typename Number, typename VectorizedArrayType>
2489inline unsigned int
2491{
2492 return *(task_info.cell_partition_data.end() - 2);
2493}
2494
2495
2496
2497template <int dim, typename Number, typename VectorizedArrayType>
2498inline unsigned int
2500{
2501 return *(task_info.cell_partition_data.end() - 1) -
2502 *(task_info.cell_partition_data.end() - 2);
2503}
2504
2505
2506
2507template <int dim, typename Number, typename VectorizedArrayType>
2508inline unsigned int
2510{
2511 if (task_info.face_partition_data.empty())
2512 return 0;
2513 return task_info.face_partition_data.back();
2514}
2515
2516
2517
2518template <int dim, typename Number, typename VectorizedArrayType>
2519inline unsigned int
2521{
2522 if (task_info.face_partition_data.empty())
2523 return 0;
2524 return task_info.boundary_partition_data.back() -
2525 task_info.face_partition_data.back();
2526}
2527
2528
2529
2530template <int dim, typename Number, typename VectorizedArrayType>
2531inline unsigned int
2533{
2534 if (task_info.face_partition_data.empty())
2535 return 0;
2536 return face_info.faces.size() - task_info.boundary_partition_data.back();
2537}
2538
2539
2540
2541template <int dim, typename Number, typename VectorizedArrayType>
2542inline types::boundary_id
2544 const unsigned int face_batch_index) const
2545{
2546 Assert(face_batch_index >= task_info.boundary_partition_data[0] &&
2547 face_batch_index < task_info.boundary_partition_data.back(),
2548 ExcIndexRange(face_batch_index,
2549 task_info.boundary_partition_data[0],
2550 task_info.boundary_partition_data.back()));
2551 return types::boundary_id(face_info.faces[face_batch_index].exterior_face_no);
2552}
2553
2554
2555
2556template <int dim, typename Number, typename VectorizedArrayType>
2557inline std::array<types::boundary_id, VectorizedArrayType::size()>
2559 const unsigned int cell_batch_index,
2560 const unsigned int face_number) const
2561{
2562 AssertIndexRange(cell_batch_index, n_cell_batches());
2564 Assert(face_info.cell_and_face_boundary_id.size(0) >= n_cell_batches(),
2566 std::array<types::boundary_id, VectorizedArrayType::size()> result;
2567 result.fill(numbers::invalid_boundary_id);
2568 for (unsigned int v = 0;
2569 v < n_active_entries_per_cell_batch(cell_batch_index);
2570 ++v)
2571 result[v] =
2572 face_info.cell_and_face_boundary_id(cell_batch_index, face_number, v);
2573 return result;
2574}
2575
2576
2577
2578template <int dim, typename Number, typename VectorizedArrayType>
2579inline const internal::MatrixFreeFunctions::
2580 MappingInfo<dim, Number, VectorizedArrayType> &
2582{
2583 return mapping_info;
2584}
2585
2586
2587
2588template <int dim, typename Number, typename VectorizedArrayType>
2591 const unsigned int dof_index) const
2592{
2593 AssertIndexRange(dof_index, n_components());
2594 return dof_info[dof_index];
2595}
2596
2597
2598
2599template <int dim, typename Number, typename VectorizedArrayType>
2600inline unsigned int
2602{
2603 return constraint_pool_row_index.size() - 1;
2604}
2605
2606
2607
2608template <int dim, typename Number, typename VectorizedArrayType>
2609inline const Number *
2611 const unsigned int row) const
2612{
2613 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2614 return constraint_pool_data.empty() ?
2615 nullptr :
2616 constraint_pool_data.data() + constraint_pool_row_index[row];
2617}
2618
2619
2620
2621template <int dim, typename Number, typename VectorizedArrayType>
2622inline const Number *
2624 const unsigned int row) const
2625{
2626 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2627 return constraint_pool_data.empty() ?
2628 nullptr :
2629 constraint_pool_data.data() + constraint_pool_row_index[row + 1];
2630}
2631
2632
2633
2634template <int dim, typename Number, typename VectorizedArrayType>
2635inline std::pair<unsigned int, unsigned int>
2637 const std::pair<unsigned int, unsigned int> &range,
2638 const unsigned int degree,
2639 const unsigned int dof_handler_component) const
2640{
2641 if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2642 {
2644 dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2646 dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2647 if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2648 return range;
2649 else
2650 return {range.second, range.second};
2651 }
2652
2653 const unsigned int fe_index =
2654 dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2655 if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2656 return {range.second, range.second};
2657 else
2658 return create_cell_subrange_hp_by_index(range,
2659 fe_index,
2660 dof_handler_component);
2661}
2662
2663
2664
2665template <int dim, typename Number, typename VectorizedArrayType>
2666inline bool
2668 const unsigned int cell_batch_index) const
2669{
2670 AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2671 return VectorizedArrayType::size() > 1 &&
2672 cell_level_index[(cell_batch_index + 1) * VectorizedArrayType::size() -
2673 1] == cell_level_index[(cell_batch_index + 1) *
2674 VectorizedArrayType::size() -
2675 2];
2676}
2677
2678
2679
2680template <int dim, typename Number, typename VectorizedArrayType>
2681unsigned int
2683{
2684 return shape_info.size(2);
2685}
2686
2687
2688template <int dim, typename Number, typename VectorizedArrayType>
2689unsigned int
2691 const std::pair<unsigned int, unsigned int> range) const
2692{
2693 const auto &fe_indices = dof_info[0].cell_active_fe_index;
2694
2695 if (fe_indices.empty() == true ||
2696 dof_handlers[0]->get_fe_collection().size() == 1)
2697 return 0;
2698
2699 const auto index = fe_indices[range.first];
2700
2701 for (unsigned int i = range.first; i < range.second; ++i)
2702 AssertDimension(index, fe_indices[i]);
2703
2704 return index;
2705}
2706
2707
2708
2709template <int dim, typename Number, typename VectorizedArrayType>
2710unsigned int
2712 const std::pair<unsigned int, unsigned int> range,
2713 const bool is_interior_face) const
2714{
2715 const auto &fe_indices = dof_info[0].cell_active_fe_index;
2716
2717 if (fe_indices.empty() == true)
2718 return 0;
2719
2720 if (is_interior_face)
2721 {
2722 const unsigned int index =
2723 fe_indices[face_info.faces[range.first].cells_interior[0] /
2724 VectorizedArrayType::size()];
2725
2726 for (unsigned int i = range.first; i < range.second; ++i)
2727 AssertDimension(index,
2728 fe_indices[face_info.faces[i].cells_interior[0] /
2729 VectorizedArrayType::size()]);
2730
2731 return index;
2732 }
2733 else
2734 {
2735 const unsigned int index =
2736 fe_indices[face_info.faces[range.first].cells_exterior[0] /
2737 VectorizedArrayType::size()];
2738
2739 for (unsigned int i = range.first; i < range.second; ++i)
2740 AssertDimension(index,
2741 fe_indices[face_info.faces[i].cells_exterior[0] /
2742 VectorizedArrayType::size()]);
2743
2744 return index;
2745 }
2746}
2747
2748
2749
2750template <int dim, typename Number, typename VectorizedArrayType>
2751inline unsigned int
2753 const unsigned int cell_batch_index) const
2754{
2755 Assert(!dof_info.empty(), ExcNotInitialized());
2756 AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2757 const std::vector<unsigned char> &n_lanes_filled =
2758 dof_info[0].n_vectorization_lanes_filled
2760 AssertIndexRange(cell_batch_index, n_lanes_filled.size());
2761
2762 return n_lanes_filled[cell_batch_index];
2763}
2764
2765
2766
2767template <int dim, typename Number, typename VectorizedArrayType>
2768inline unsigned int
2770 const unsigned int face_batch_index) const
2771{
2772 AssertIndexRange(face_batch_index, face_info.faces.size());
2773 Assert(!dof_info.empty(), ExcNotInitialized());
2774 const std::vector<unsigned char> &n_lanes_filled =
2775 dof_info[0].n_vectorization_lanes_filled
2777 AssertIndexRange(face_batch_index, n_lanes_filled.size());
2778 return n_lanes_filled[face_batch_index];
2779}
2780
2781
2782
2783template <int dim, typename Number, typename VectorizedArrayType>
2784inline unsigned int
2786 const unsigned int dof_handler_index,
2787 const unsigned int active_fe_index) const
2788{
2789 return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2790}
2791
2792
2793
2794template <int dim, typename Number, typename VectorizedArrayType>
2795inline unsigned int
2797 const unsigned int quad_index,
2798 const unsigned int active_fe_index) const
2799{
2800 AssertIndexRange(quad_index, mapping_info.cell_data.size());
2801 return mapping_info.cell_data[quad_index]
2802 .descriptor[active_fe_index]
2803 .n_q_points;
2804}
2805
2806
2807
2808template <int dim, typename Number, typename VectorizedArrayType>
2809inline unsigned int
2811 const unsigned int dof_handler_index,
2812 const unsigned int active_fe_index) const
2813{
2814 return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2815}
2816
2817
2818
2819template <int dim, typename Number, typename VectorizedArrayType>
2820inline unsigned int
2822 const unsigned int quad_index,
2823 const unsigned int active_fe_index) const
2824{
2825 AssertIndexRange(quad_index, mapping_info.face_data.size());
2826 return mapping_info.face_data[quad_index]
2827 .descriptor[active_fe_index]
2828 .n_q_points;
2829}
2830
2831
2832
2833template <int dim, typename Number, typename VectorizedArrayType>
2834inline const IndexSet &
2836 const unsigned int dof_handler_index) const
2837{
2838 return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2839}
2840
2841
2842
2843template <int dim, typename Number, typename VectorizedArrayType>
2844inline const IndexSet &
2846 const unsigned int dof_handler_index) const
2847{
2848 return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2849}
2850
2851
2852
2853template <int dim, typename Number, typename VectorizedArrayType>
2856 const unsigned int dof_handler_index,
2857 const unsigned int index_quad,
2858 const unsigned int index_fe,
2859 const unsigned int active_fe_index,
2860 const unsigned int active_quad_index) const
2861{
2862 AssertIndexRange(dof_handler_index, dof_info.size());
2863 const unsigned int ind =
2864 dof_info[dof_handler_index].global_base_element_offset + index_fe;
2865 AssertIndexRange(ind, shape_info.size(0));
2866 AssertIndexRange(index_quad, shape_info.size(1));
2867 AssertIndexRange(active_fe_index, shape_info.size(2));
2868 AssertIndexRange(active_quad_index, shape_info.size(3));
2869 return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2870}
2871
2872
2873
2874template <int dim, typename Number, typename VectorizedArrayType>
2876 VectorizedArrayType::size()> &
2878 const unsigned int face_batch_index) const
2879{
2880 AssertIndexRange(face_batch_index, face_info.faces.size());
2881 return face_info.faces[face_batch_index];
2882}
2883
2884
2885
2886template <int dim, typename Number, typename VectorizedArrayType>
2887inline const Table<3, unsigned int> &
2889 const
2890{
2891 return face_info.cell_and_face_to_plain_faces;
2892}
2893
2894
2895
2896template <int dim, typename Number, typename VectorizedArrayType>
2897inline const Quadrature<dim> &
2899 const unsigned int quad_index,
2900 const unsigned int active_fe_index) const
2901{
2902 AssertIndexRange(quad_index, mapping_info.cell_data.size());
2903 return mapping_info.cell_data[quad_index]
2904 .descriptor[active_fe_index]
2905 .quadrature;
2906}
2907
2908
2909
2910template <int dim, typename Number, typename VectorizedArrayType>
2911inline const Quadrature<dim - 1> &
2913 const unsigned int quad_index,
2914 const unsigned int active_fe_index) const
2915{
2916 AssertIndexRange(quad_index, mapping_info.face_data.size());
2917 return mapping_info.face_data[quad_index]
2918 .descriptor[active_fe_index]
2919 .quadrature;
2920}
2921
2922
2923
2924template <int dim, typename Number, typename VectorizedArrayType>
2925inline unsigned int
2927 const std::pair<unsigned int, unsigned int> range) const
2928{
2929 auto result = get_cell_category(range.first);
2930
2931 for (unsigned int i = range.first; i < range.second; ++i)
2932 result = std::max(result, get_cell_category(i));
2933
2934 return result;
2935}
2936
2937
2938
2939template <int dim, typename Number, typename VectorizedArrayType>
2940inline std::pair<unsigned int, unsigned int>
2942 const std::pair<unsigned int, unsigned int> range) const
2943{
2944 auto result = get_face_category(range.first);
2945
2946 for (unsigned int i = range.first; i < range.second; ++i)
2947 {
2948 result.first = std::max(result.first, get_face_category(i).first);
2949 result.second = std::max(result.second, get_face_category(i).second);
2950 }
2951
2952 return result;
2953}
2954
2955
2956
2957template <int dim, typename Number, typename VectorizedArrayType>
2958inline unsigned int
2960 const unsigned int cell_batch_index) const
2961{
2962 AssertIndexRange(0, dof_info.size());
2963 AssertIndexRange(cell_batch_index, dof_info[0].cell_active_fe_index.size());
2964 if (dof_info[0].cell_active_fe_index.empty())
2965 return 0;
2966 else
2967 return dof_info[0].cell_active_fe_index[cell_batch_index];
2968}
2969
2970
2971
2972template <int dim, typename Number, typename VectorizedArrayType>
2973inline std::pair<unsigned int, unsigned int>
2975 const unsigned int face_batch_index) const
2976{
2977 AssertIndexRange(face_batch_index, face_info.faces.size());
2978 if (dof_info[0].cell_active_fe_index.empty())
2979 return std::make_pair(0U, 0U);
2980
2981 std::pair<unsigned int, unsigned int> result = std::make_pair(0U, 0U);
2982 for (unsigned int v = 0;
2983 v < VectorizedArrayType::size() &&
2984 face_info.faces[face_batch_index].cells_interior[v] !=
2986 ++v)
2987 result.first = std::max(
2988 result.first,
2989 dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
2990 .cells_interior[v] /
2991 VectorizedArrayType::size()]);
2992 if (face_info.faces[face_batch_index].cells_exterior[0] !=
2994 for (unsigned int v = 0;
2995 v < VectorizedArrayType::size() &&
2996 face_info.faces[face_batch_index].cells_exterior[v] !=
2998 ++v)
2999 result.second = std::max(
3000 result.second,
3001 dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
3002 .cells_exterior[v] /
3003 VectorizedArrayType::size()]);
3004 else
3005 result.second = numbers::invalid_unsigned_int;
3006 return result;
3007}
3008
3009
3010
3011template <int dim, typename Number, typename VectorizedArrayType>
3012inline bool
3014{
3015 return indices_are_initialized;
3016}
3017
3018
3019
3020template <int dim, typename Number, typename VectorizedArrayType>
3021inline bool
3023{
3024 return mapping_is_initialized;
3025}
3026
3027
3028template <int dim, typename Number, typename VectorizedArrayType>
3029inline unsigned int
3031{
3032 return mg_level;
3033}
3034
3035
3036
3037template <int dim, typename Number, typename VectorizedArrayType>
3040{
3041 using list_type =
3042 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3043 list_type &data = scratch_pad.get();
3044 for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
3045 if (it->first == false)
3046 {
3047 it->first = true;
3048 return &it->second;
3049 }
3050 data.emplace_front(true, AlignedVector<VectorizedArrayType>());
3051 return &data.front().second;
3052}
3053
3054
3055
3056template <int dim, typename Number, typename VectorizedArrayType>
3057void
3059 const AlignedVector<VectorizedArrayType> *scratch) const
3060{
3061 using list_type =
3062 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3063 list_type &data = scratch_pad.get();
3064 for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
3065 if (&it->second == scratch)
3066 {
3067 Assert(it->first == true, ExcInternalError());
3068 it->first = false;
3069 return;
3070 }
3071 AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
3072}
3073
3074
3075
3076template <int dim, typename Number, typename VectorizedArrayType>
3080{
3081 for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
3082 scratch_pad_non_threadsafe.begin();
3083 it != scratch_pad_non_threadsafe.end();
3084 ++it)
3085 if (it->first == false)
3086 {
3087 it->first = true;
3088 return &it->second;
3089 }
3090 scratch_pad_non_threadsafe.push_front(
3091 std::make_pair(true, AlignedVector<Number>()));
3092 return &scratch_pad_non_threadsafe.front().second;
3093}
3094
3095
3096
3097template <int dim, typename Number, typename VectorizedArrayType>
3098void
3101 const AlignedVector<Number> *scratch) const
3102{
3103 for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
3104 scratch_pad_non_threadsafe.begin();
3105 it != scratch_pad_non_threadsafe.end();
3106 ++it)
3107 if (&it->second == scratch)
3108 {
3109 Assert(it->first == true, ExcInternalError());
3110 it->first = false;
3111 return;
3112 }
3113 AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
3114}
3115
3116
3117
3118// ------------------------------ reinit functions ---------------------------
3119
3120namespace internal
3121{
3122 namespace MatrixFreeImplementation
3123 {
3124 template <int dim, int spacedim>
3125 inline std::vector<IndexSet>
3126 extract_locally_owned_index_sets(
3127 const std::vector<const DoFHandler<dim, spacedim> *> &dofh,
3128 const unsigned int level)
3129 {
3130 std::vector<IndexSet> locally_owned_set;
3131 locally_owned_set.reserve(dofh.size());
3132 for (unsigned int j = 0; j < dofh.size(); ++j)
3134 locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
3135 else
3136 locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
3137 return locally_owned_set;
3138 }
3139 } // namespace MatrixFreeImplementation
3140} // namespace internal
3141
3142
3143
3144template <int dim, typename Number, typename VectorizedArrayType>
3145template <typename QuadratureType, typename number2, typename MappingType>
3146void
3148 const MappingType &mapping,
3149 const DoFHandler<dim> &dof_handler,
3150 const AffineConstraints<number2> &constraints_in,
3151 const QuadratureType &quad,
3153 &additional_data)
3154{
3155 std::vector<const DoFHandler<dim, dim> *> dof_handlers;
3156 std::vector<const AffineConstraints<number2> *> constraints;
3157
3158 dof_handlers.push_back(&dof_handler);
3159 constraints.push_back(&constraints_in);
3160
3161 std::vector<IndexSet> locally_owned_sets =
3162 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3163 dof_handlers, additional_data.mg_level);
3164
3165 std::vector<hp::QCollection<dim>> quad_hp;
3166 quad_hp.emplace_back(quad);
3167
3168 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3169 dof_handlers,
3170 constraints,
3171 locally_owned_sets,
3172 quad_hp,
3173 additional_data);
3174}
3175
3176
3177
3178template <int dim, typename Number, typename VectorizedArrayType>
3179template <typename QuadratureType, typename number2, typename MappingType>
3180void
3182 const MappingType &mapping,
3183 const std::vector<const DoFHandler<dim> *> &dof_handler,
3184 const std::vector<const AffineConstraints<number2> *> &constraint,
3185 const QuadratureType &quad,
3187 &additional_data)
3188{
3189 std::vector<IndexSet> locally_owned_set =
3190 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3191 dof_handler, additional_data.mg_level);
3192 std::vector<hp::QCollection<dim>> quad_hp;
3193 quad_hp.emplace_back(quad);
3194
3195 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3196 dof_handler,
3197 constraint,
3198 locally_owned_set,
3199 quad_hp,
3200 additional_data);
3201}
3202
3203
3204
3205template <int dim, typename Number, typename VectorizedArrayType>
3206template <typename QuadratureType, typename number2, typename MappingType>
3207void
3209 const MappingType &mapping,
3210 const std::vector<const DoFHandler<dim> *> &dof_handler,
3211 const std::vector<const AffineConstraints<number2> *> &constraint,
3212 const std::vector<QuadratureType> &quad,
3214 &additional_data)
3215{
3216 std::vector<IndexSet> locally_owned_set =
3217 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3218 dof_handler, additional_data.mg_level);
3219 std::vector<hp::QCollection<dim>> quad_hp;
3220 for (unsigned int q = 0; q < quad.size(); ++q)
3221 quad_hp.emplace_back(quad[q]);
3222
3223 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3224 dof_handler,
3225 constraint,
3226 locally_owned_set,
3227 quad_hp,
3228 additional_data);
3229}
3230
3231
3232
3233// ------------------------------ implementation of loops --------------------
3234
3235// internal helper functions that define how to call MPI data exchange
3236// functions: for generic vectors, do nothing at all. For distributed vectors,
3237// call update_ghost_values_start function and so on. If we have collections
3238// of vectors, just do the individual functions of the components. In order to
3239// keep ghost values consistent (whether we are in read or write mode), we
3240// also reset the values at the end. the whole situation is a bit complicated
3241// by the fact that we need to treat block vectors differently, which use some
3242// additional helper functions to select the blocks and template magic.
3243namespace internal
3244{
3248 template <int dim, typename Number, typename VectorizedArrayType>
3249 struct VectorDataExchange
3250 {
3251 // A shift for the MPI messages to reduce the risk for accidental
3252 // interaction with other open communications that a user program might
3253 // set up (parallel vectors support unfinished communication). We let
3254 // the other vectors use the first 20 assigned numbers and start the
3255 // matrix-free communication.
3256 static constexpr unsigned int channel_shift = 20;
3257
3258
3259
3264 VectorDataExchange(
3265 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3266 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3267 DataAccessOnFaces vector_face_access,
3268 const unsigned int n_components)
3269 : matrix_free(matrix_free)
3270 , vector_face_access(
3271 matrix_free.get_task_info().face_partition_data.empty() ?
3272 ::MatrixFree<dim, Number, VectorizedArrayType>::
3273 DataAccessOnFaces::unspecified :
3274 vector_face_access)
3275 , ghosts_were_set(false)
3276# ifdef DEAL_II_WITH_MPI
3277 , tmp_data(n_components)
3278 , requests(n_components)
3279# endif
3280 {
3281 (void)n_components;
3282 if (this->vector_face_access !=
3284 DataAccessOnFaces::unspecified)
3285 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3287 matrix_free.get_dof_info(c).vector_exchanger_face_variants.size(),
3288 5);
3289 }
3290
3291
3292
3296 ~VectorDataExchange() // NOLINT
3297 {
3298# ifdef DEAL_II_WITH_MPI
3299 for (unsigned int i = 0; i < tmp_data.size(); ++i)
3300 if (tmp_data[i] != nullptr)
3301 matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3302# endif
3303 }
3304
3305
3306
3311 template <typename VectorType>
3312 unsigned int
3313 find_vector_in_mf(const VectorType &vec,
3314 const bool check_global_compatibility = true) const
3315 {
3316 // case 1: vector was set up with MatrixFree::initialize_dof_vector()
3317 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3318 if (vec.get_partitioner().get() ==
3319 matrix_free.get_dof_info(c).vector_partitioner.get())
3320 return c;
3321
3322 // case 2: user provided own partitioner (compatibility mode)
3323 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3324 if (check_global_compatibility ?
3325 vec.get_partitioner()->is_globally_compatible(
3326 *matrix_free.get_dof_info(c).vector_partitioner) :
3327 vec.get_partitioner()->is_compatible(
3328 *matrix_free.get_dof_info(c).vector_partitioner))
3329 return c;
3330
3331 Assert(false,
3332 ExcNotImplemented("Could not find partitioner that fits vector"));
3333
3335 }
3336
3337
3338
3344 get_partitioner(const unsigned int mf_component) const
3345 {
3346 AssertDimension(matrix_free.get_dof_info(mf_component)
3347 .vector_exchanger_face_variants.size(),
3348 5);
3349 if (vector_face_access ==
3351 DataAccessOnFaces::none)
3352 return *matrix_free.get_dof_info(mf_component)
3353 .vector_exchanger_face_variants[0];
3354 else if (vector_face_access ==
3356 DataAccessOnFaces::values)
3357 return *matrix_free.get_dof_info(mf_component)
3358 .vector_exchanger_face_variants[1];
3359 else if (vector_face_access ==
3361 DataAccessOnFaces::gradients)
3362 return *matrix_free.get_dof_info(mf_component)
3363 .vector_exchanger_face_variants[2];
3364 else if (vector_face_access ==
3366 DataAccessOnFaces::values_all_faces)
3367 return *matrix_free.get_dof_info(mf_component)
3368 .vector_exchanger_face_variants[3];
3369 else if (vector_face_access ==
3371 DataAccessOnFaces::gradients_all_faces)
3372 return *matrix_free.get_dof_info(mf_component)
3373 .vector_exchanger_face_variants[4];
3374 else
3375 return *matrix_free.get_dof_info(mf_component).vector_exchanger.get();
3376 }
3377
3378
3379
3383 template <typename VectorType,
3384 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3385 * = nullptr>
3386 void
3387 update_ghost_values_start(const unsigned int /*component_in_block_vector*/,
3388 const VectorType & /*vec*/)
3389 {}
3390
3391
3396 template <typename VectorType,
3397 std::enable_if_t<!has_update_ghost_values_start<VectorType> &&
3398 !is_not_parallel_vector<VectorType>,
3399 VectorType> * = nullptr>
3400 void
3401 update_ghost_values_start(const unsigned int component_in_block_vector,
3402 const VectorType &vec)
3403 {
3404 (void)component_in_block_vector;
3405 const bool ghosts_set = vec.has_ghost_elements();
3406
3407 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3408 ghosts_set == false,
3410
3411 if (ghosts_set)
3412 ghosts_were_set = true;
3413
3414 vec.update_ghost_values();
3415 }
3416
3417
3418
3424 template <typename VectorType,
3425 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3426 !has_exchange_on_subset<VectorType>,
3427 VectorType> * = nullptr>
3428 void
3429 update_ghost_values_start(const unsigned int component_in_block_vector,
3430 const VectorType &vec)
3431 {
3432 (void)component_in_block_vector;
3433 const bool ghosts_set = vec.has_ghost_elements();
3434
3435 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3436 ghosts_set == false,
3438
3439 if (ghosts_set)
3440 ghosts_were_set = true;
3441
3442 vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3443 }
3444
3445
3446
3453 template <typename VectorType,
3454 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3455 has_exchange_on_subset<VectorType>,
3456 VectorType> * = nullptr>
3457 void
3458 update_ghost_values_start(const unsigned int component_in_block_vector,
3459 const VectorType &vec)
3460 {
3461 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3462 "Type mismatch between VectorType and VectorDataExchange");
3463 (void)component_in_block_vector;
3464 const bool ghosts_set = vec.has_ghost_elements();
3465
3466 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3467 ghosts_set == false,
3469
3470 if (ghosts_set)
3471 ghosts_were_set = true;
3472
3473 if (vec.size() != 0)
3474 {
3475# ifdef DEAL_II_WITH_MPI
3476 const unsigned int mf_component = find_vector_in_mf(vec);
3477
3478 const auto &part = get_partitioner(mf_component);
3479
3480 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3481 part.n_import_sm_procs() == 0)
3482 return;
3483
3484 tmp_data[component_in_block_vector] =
3485 matrix_free.acquire_scratch_data_non_threadsafe();
3486 tmp_data[component_in_block_vector]->resize_fast(
3487 part.n_import_indices());
3488 AssertDimension(requests.size(), tmp_data.size());
3489
3490 part.export_to_ghosted_array_start(
3491 component_in_block_vector * 2 + channel_shift,
3492 ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3493 vec.shared_vector_data(),
3494 ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3495 part.locally_owned_size(),
3496 matrix_free.get_dof_info(mf_component)
3497 .vector_partitioner->n_ghost_indices()),
3498 ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3499 part.n_import_indices()),
3500 this->requests[component_in_block_vector]);
3501# endif
3502 }
3503 }
3504
3505
3506
3511 template <typename VectorType,
3512 std::enable_if_t<!has_update_ghost_values_start<VectorType>,
3513 VectorType> * = nullptr>
3514 void
3515 update_ghost_values_finish(const unsigned int /*component_in_block_vector*/,
3516 const VectorType & /*vec*/)
3517 {}
3518
3519
3520
3526 template <typename VectorType,
3527 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3528 !has_exchange_on_subset<VectorType>,
3529 VectorType> * = nullptr>
3530 void
3531 update_ghost_values_finish(const unsigned int component_in_block_vector,
3532 const VectorType &vec)
3533 {
3534 (void)component_in_block_vector;
3535 vec.update_ghost_values_finish();
3536 }
3537
3538
3539
3546 template <typename VectorType,
3547 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3548 has_exchange_on_subset<VectorType>,
3549 VectorType> * = nullptr>
3550 void
3551 update_ghost_values_finish(const unsigned int component_in_block_vector,
3552 const VectorType &vec)
3553 {
3554 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3555 "Type mismatch between VectorType and VectorDataExchange");
3556 (void)component_in_block_vector;
3557
3558 if (vec.size() != 0)
3559 {
3560# ifdef DEAL_II_WITH_MPI
3561 AssertIndexRange(component_in_block_vector, tmp_data.size());
3562 AssertDimension(requests.size(), tmp_data.size());
3563
3564 const unsigned int mf_component = find_vector_in_mf(vec);
3565
3566 const auto &part = get_partitioner(mf_component);
3567
3568 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3569 part.n_import_sm_procs() != 0)
3570 {
3571 part.export_to_ghosted_array_finish(
3572 ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3573 vec.shared_vector_data(),
3574 ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3575 part.locally_owned_size(),
3576 matrix_free.get_dof_info(mf_component)
3577 .vector_partitioner->n_ghost_indices()),
3578 this->requests[component_in_block_vector]);
3579
3580 matrix_free.release_scratch_data_non_threadsafe(
3581 tmp_data[component_in_block_vector]);
3582 tmp_data[component_in_block_vector] = nullptr;
3583 }
3584# endif
3585 }
3586 // let vector know that ghosts are being updated and we can read from
3587 // them
3588 vec.set_ghost_state(true);
3589 }
3590
3591
3592
3596 template <typename VectorType,
3597 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3598 * = nullptr>
3599 void
3600 compress_start(const unsigned int /*component_in_block_vector*/,
3601 VectorType & /*vec*/)
3602 {}
3603
3604
3605
3610 template <typename VectorType,
3611 std::enable_if_t<!has_compress_start<VectorType> &&
3612 !is_not_parallel_vector<VectorType>,
3613 VectorType> * = nullptr>
3614 void
3615 compress_start(const unsigned int component_in_block_vector,
3616 VectorType &vec)
3617 {
3618 (void)component_in_block_vector;
3619 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3620 vec.compress(VectorOperation::add);
3621 }
3622
3623
3624
3630 template <typename VectorType,
3631 std::enable_if_t<has_compress_start<VectorType> &&
3632 !has_exchange_on_subset<VectorType>,
3633 VectorType> * = nullptr>
3634 void
3635 compress_start(const unsigned int component_in_block_vector,
3636 VectorType &vec)
3637 {
3638 (void)component_in_block_vector;
3639 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3640 vec.compress_start(component_in_block_vector + channel_shift);
3641 }
3642
3643
3644
3651 template <typename VectorType,
3652 std::enable_if_t<has_compress_start<VectorType> &&
3653 has_exchange_on_subset<VectorType>,
3654 VectorType> * = nullptr>
3655 void
3656 compress_start(const unsigned int component_in_block_vector,
3657 VectorType &vec)
3658 {
3659 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3660 "Type mismatch between VectorType and VectorDataExchange");
3661 (void)component_in_block_vector;
3662 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3663
3664 if (vec.size() != 0)
3665 {
3666# ifdef DEAL_II_WITH_MPI
3667 const unsigned int mf_component = find_vector_in_mf(vec);
3668
3669 const auto &part = get_partitioner(mf_component);
3670
3671 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3672 part.n_import_sm_procs() == 0)
3673 return;
3674
3675 tmp_data[component_in_block_vector] =
3676 matrix_free.acquire_scratch_data_non_threadsafe();
3677 tmp_data[component_in_block_vector]->resize_fast(
3678 part.n_import_indices());
3679 AssertDimension(requests.size(), tmp_data.size());
3680
3681 part.import_from_ghosted_array_start(
3683 component_in_block_vector * 2 + channel_shift,
3684 ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3685 vec.shared_vector_data(),
3686 ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3687 matrix_free.get_dof_info(mf_component)
3688 .vector_partitioner->n_ghost_indices()),
3689 ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3690 part.n_import_indices()),
3691 this->requests[component_in_block_vector]);
3692# endif
3693 }
3694 }
3695
3696
3697
3702 template <
3703 typename VectorType,
3704 std::enable_if_t<!has_compress_start<VectorType>, VectorType> * = nullptr>
3705 void
3706 compress_finish(const unsigned int /*component_in_block_vector*/,
3707 VectorType & /*vec*/)
3708 {}
3709
3710
3711
3717 template <typename VectorType,
3718 std::enable_if_t<has_compress_start<VectorType> &&
3719 !has_exchange_on_subset<VectorType>,
3720 VectorType> * = nullptr>
3721 void
3722 compress_finish(const unsigned int component_in_block_vector,
3723 VectorType &vec)
3724 {
3725 (void)component_in_block_vector;
3726 vec.compress_finish(VectorOperation::add);
3727 }
3728
3729
3730
3737 template <typename VectorType,
3738 std::enable_if_t<has_compress_start<VectorType> &&
3739 has_exchange_on_subset<VectorType>,
3740 VectorType> * = nullptr>
3741 void
3742 compress_finish(const unsigned int component_in_block_vector,
3743 VectorType &vec)
3744 {
3745 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3746 "Type mismatch between VectorType and VectorDataExchange");
3747 (void)component_in_block_vector;
3748 if (vec.size() != 0)
3749 {
3750# ifdef DEAL_II_WITH_MPI
3751 AssertIndexRange(component_in_block_vector, tmp_data.size());
3752 AssertDimension(requests.size(), tmp_data.size());
3753
3754 const unsigned int mf_component = find_vector_in_mf(vec);
3755
3756 const auto &part = get_partitioner(mf_component);
3757
3758 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3759 part.n_import_sm_procs() != 0)
3760 {
3761 part.import_from_ghosted_array_finish(
3763 ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3764 vec.shared_vector_data(),
3765 ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3766 matrix_free.get_dof_info(mf_component)
3767 .vector_partitioner->n_ghost_indices()),
3769 tmp_data[component_in_block_vector]->begin(),
3770 part.n_import_indices()),
3771 this->requests[component_in_block_vector]);
3772
3773 matrix_free.release_scratch_data_non_threadsafe(
3774 tmp_data[component_in_block_vector]);
3775 tmp_data[component_in_block_vector] = nullptr;
3776 }
3777
3779 {
3780 const int ierr =
3781 MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3782 AssertThrowMPI(ierr);
3783 }
3784# endif
3785 }
3786 }
3787
3788
3789
3793 template <typename VectorType,
3794 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3795 * = nullptr>
3796 void
3797 reset_ghost_values(const VectorType & /*vec*/) const
3798 {}
3799
3800
3801
3806 template <typename VectorType,
3807 std::enable_if_t<!has_exchange_on_subset<VectorType> &&
3808 !is_not_parallel_vector<VectorType>,
3809 VectorType> * = nullptr>
3810 void
3811 reset_ghost_values(const VectorType &vec) const
3812 {
3813 if (ghosts_were_set == true)
3814 return;
3815
3816 vec.zero_out_ghost_values();
3817 }
3818
3819
3820
3826 template <typename VectorType,
3827 std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3828 * = nullptr>
3829 void
3830 reset_ghost_values(const VectorType &vec) const
3831 {
3832 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3833 "Type mismatch between VectorType and VectorDataExchange");
3834 if (ghosts_were_set == true)
3835 return;
3836
3837 if (vec.size() != 0)
3838 {
3839# ifdef DEAL_II_WITH_MPI
3840 AssertDimension(requests.size(), tmp_data.size());
3841
3842 const unsigned int mf_component = find_vector_in_mf(vec);
3843
3844 const auto &part = get_partitioner(mf_component);
3845
3846 if (part.n_ghost_indices() > 0)
3847 {
3848 part.reset_ghost_values(ArrayView<Number>(
3850 .begin() +
3851 part.locally_owned_size(),
3852 matrix_free.get_dof_info(mf_component)
3853 .vector_partitioner->n_ghost_indices()));
3854 }
3855
3856# endif
3857 }
3858 // let vector know that it's not ghosted anymore
3859 vec.set_ghost_state(false);
3860 }
3861
3862
3863
3869 template <typename VectorType,
3870 std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3871 * = nullptr>
3872 void
3873 zero_vector_region(const unsigned int range_index, VectorType &vec) const
3874 {
3875 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3876 "Type mismatch between VectorType and VectorDataExchange");
3877 if (range_index == numbers::invalid_unsigned_int)
3878 vec = Number();
3879 else
3880 {
3881 const unsigned int mf_component = find_vector_in_mf(vec, false);
3883 matrix_free.get_dof_info(mf_component);
3884 Assert(dof_info.vector_zero_range_list_index.empty() == false,
3886
3887 Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3889 AssertIndexRange(range_index,
3890 dof_info.vector_zero_range_list_index.size() - 1);
3891 for (unsigned int id =
3892 dof_info.vector_zero_range_list_index[range_index];
3893 id != dof_info.vector_zero_range_list_index[range_index + 1];
3894 ++id)
3895 std::memset(vec.begin() + dof_info.vector_zero_range_list[id].first,
3896 0,
3897 (dof_info.vector_zero_range_list[id].second -
3898 dof_info.vector_zero_range_list[id].first) *
3899 sizeof(Number));
3900 }
3901 }
3902
3903
3904
3910 template <typename VectorType,
3911 std::enable_if_t<!has_exchange_on_subset<VectorType>, VectorType>
3912 * = nullptr,
3913 typename VectorType::value_type * = nullptr>
3914 void
3915 zero_vector_region(const unsigned int range_index, VectorType &vec) const
3916 {
3917 if (range_index == numbers::invalid_unsigned_int || range_index == 0)
3918 vec = typename VectorType::value_type();
3919 }
3920
3921
3922
3927 void
3928 zero_vector_region(const unsigned int, ...) const
3929 {
3930 Assert(false,
3931 ExcNotImplemented("Zeroing is only implemented for vector types "
3932 "which provide VectorType::value_type"));
3933 }
3934
3935
3936
3937 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
3938 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3939 DataAccessOnFaces vector_face_access;
3940 bool ghosts_were_set;
3941# ifdef DEAL_II_WITH_MPI
3942 std::vector<AlignedVector<Number> *> tmp_data;
3943 std::vector<std::vector<MPI_Request>> requests;
3944# endif
3945 }; // VectorDataExchange
3946
3947 template <typename VectorStruct>
3948 unsigned int
3949 n_components(const VectorStruct &vec);
3950
3951 template <typename VectorStruct>
3952 unsigned int
3953 n_components_block(const VectorStruct &vec, const std::bool_constant<true>)
3954 {
3955 unsigned int components = 0;
3956 for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
3957 components += n_components(vec.block(bl));
3958 return components;
3959 }
3960
3961 template <typename VectorStruct>
3962 unsigned int
3963 n_components_block(const VectorStruct &, const std::bool_constant<false>)
3964 {
3965 return 1;
3966 }
3967
3968 template <typename VectorStruct>
3969 unsigned int
3970 n_components(const VectorStruct &vec)
3971 {
3972 return n_components_block(
3973 vec, std::bool_constant<IsBlockVector<VectorStruct>::value>());
3974 }
3975
3976 template <typename VectorStruct>
3977 inline unsigned int
3978 n_components(const std::vector<VectorStruct> &vec)
3979 {
3980 unsigned int components = 0;
3981 for (unsigned int comp = 0; comp < vec.size(); ++comp)
3982 components += n_components_block(
3983 vec[comp], std::bool_constant<IsBlockVector<VectorStruct>::value>());
3984 return components;
3985 }
3986
3987 template <typename VectorStruct>
3988 inline unsigned int
3989 n_components(const std::vector<VectorStruct *> &vec)
3990 {
3991 unsigned int components = 0;
3992 for (unsigned int comp = 0; comp < vec.size(); ++comp)
3993 components += n_components_block(
3994 *vec[comp], std::bool_constant<IsBlockVector<VectorStruct>::value>());
3995 return components;
3996 }
3997
3998
3999
4000 // A helper function to identify block vectors with many components where we
4001 // should not try to overlap computations and communication because there
4002 // would be too many outstanding communication requests.
4003
4004 // default value for vectors that do not have communication_block_size
4005 template <typename VectorStruct,
4006 std::enable_if_t<!has_communication_block_size<VectorStruct>,
4007 VectorStruct> * = nullptr>
4008 constexpr unsigned int
4009 get_communication_block_size(const VectorStruct &)
4010 {
4012 }
4013
4014
4015
4016 template <typename VectorStruct,
4017 std::enable_if_t<has_communication_block_size<VectorStruct>,
4018 VectorStruct> * = nullptr>
4019 constexpr unsigned int
4020 get_communication_block_size(const VectorStruct &)
4021 {
4022 return VectorStruct::communication_block_size;
4023 }
4024
4025
4026
4027 template <typename VectorType,
4028 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType> * =
4029 nullptr>
4030 bool
4031 has_ghost_elements(const VectorType &vec)
4032 {
4033 (void)vec;
4034 return false;
4035 }
4036
4037
4038
4039 template <typename VectorType,
4040 std::enable_if_t<!is_not_parallel_vector<VectorType>, VectorType>
4041 * = nullptr>
4042 bool
4043 has_ghost_elements(const VectorType &vec)
4044 {
4045 return vec.has_ghost_elements();
4046 }
4047
4048
4049
4050 // --------------------------------------------------------------------------
4051 // below we have wrappers to distinguish between block and non-block vectors.
4052 // --------------------------------------------------------------------------
4053
4054 //
4055 // update_ghost_values_start
4056 //
4057
4058 // update_ghost_values for block vectors
4059 template <int dim,
4060 typename VectorStruct,
4061 typename Number,
4062 typename VectorizedArrayType,
4063 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4064 * = nullptr>
4065 void
4066 update_ghost_values_start(
4067 const VectorStruct &vec,
4068 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4069 const unsigned int channel = 0)
4070 {
4071 if (get_communication_block_size(vec) < vec.n_blocks())
4072 {
4073 const bool ghosts_set = vec.has_ghost_elements();
4074
4075 Assert(exchanger.matrix_free.get_task_info()
4076 .allow_ghosted_vectors_in_loops ||
4077 ghosts_set == false,
4079
4080 if (ghosts_set)
4081 exchanger.ghosts_were_set = true;
4082
4083 vec.update_ghost_values();
4084 }
4085 else
4086 {
4087 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4088 update_ghost_values_start(vec.block(i), exchanger, channel + i);
4089 }
4090 }
4091
4092
4093
4094 // update_ghost_values for non-block vectors
4095 template <int dim,
4096 typename VectorStruct,
4097 typename Number,
4098 typename VectorizedArrayType,
4099 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4100 * = nullptr>
4101 void
4102 update_ghost_values_start(
4103 const VectorStruct &vec,
4104 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4105 const unsigned int channel = 0)
4106 {
4107 exchanger.update_ghost_values_start(channel, vec);
4108 }
4109
4110
4111
4112 // update_ghost_values_start() for vector of vectors
4113 template <int dim,
4114 typename VectorStruct,
4115 typename Number,
4116 typename VectorizedArrayType>
4117 inline void
4118 update_ghost_values_start(
4119 const std::vector<VectorStruct> &vec,
4120 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4121 {
4122 unsigned int component_index = 0;
4123 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4124 {
4125 update_ghost_values_start(vec[comp], exchanger, component_index);
4126 component_index += n_components(vec[comp]);
4127 }
4128 }
4129
4130
4131
4132 // update_ghost_values_start() for vector of pointers to vectors
4133 template <int dim,
4134 typename VectorStruct,
4135 typename Number,
4136 typename VectorizedArrayType>
4137 inline void
4138 update_ghost_values_start(
4139 const std::vector<VectorStruct *> &vec,
4140 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4141 {
4142 unsigned int component_index = 0;
4143 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4144 {
4145 update_ghost_values_start(*vec[comp], exchanger, component_index);
4146 component_index += n_components(*vec[comp]);
4147 }
4148 }
4149
4150
4151
4152 //
4153 // update_ghost_values_finish
4154 //
4155
4156 // for block vectors
4157 template <int dim,
4158 typename VectorStruct,
4159 typename Number,
4160 typename VectorizedArrayType,
4161 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4162 * = nullptr>
4163 void
4164 update_ghost_values_finish(
4165 const VectorStruct &vec,
4166 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4167 const unsigned int channel = 0)
4168 {
4169 if (get_communication_block_size(vec) < vec.n_blocks())
4170 {
4171 // do nothing, everything has already been completed in the _start()
4172 // call
4173 }
4174 else
4175 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4176 update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4177 }
4178
4179
4180
4181 // for non-block vectors
4182 template <int dim,
4183 typename VectorStruct,
4184 typename Number,
4185 typename VectorizedArrayType,
4186 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4187 * = nullptr>
4188 void
4189 update_ghost_values_finish(
4190 const VectorStruct &vec,
4191 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4192 const unsigned int channel = 0)
4193 {
4194 exchanger.update_ghost_values_finish(channel, vec);
4195 }
4196
4197
4198
4199 // for vector of vectors
4200 template <int dim,
4201 typename VectorStruct,
4202 typename Number,
4203 typename VectorizedArrayType>
4204 inline void
4205 update_ghost_values_finish(
4206 const std::vector<VectorStruct> &vec,
4207 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4208 {
4209 unsigned int component_index = 0;
4210 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4211 {
4212 update_ghost_values_finish(vec[comp], exchanger, component_index);
4213 component_index += n_components(vec[comp]);
4214 }
4215 }
4216
4217
4218
4219 // for vector of pointers to vectors
4220 template <int dim,
4221 typename VectorStruct,
4222 typename Number,
4223 typename VectorizedArrayType>
4224 inline void
4225 update_ghost_values_finish(
4226 const std::vector<VectorStruct *> &vec,
4227 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4228 {
4229 unsigned int component_index = 0;
4230 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4231 {
4232 update_ghost_values_finish(*vec[comp], exchanger, component_index);
4233 component_index += n_components(*vec[comp]);
4234 }
4235 }
4236
4237
4238
4239 //
4240 // compress_start
4241 //
4242
4243 // for block vectors
4244 template <int dim,
4245 typename VectorStruct,
4246 typename Number,
4247 typename VectorizedArrayType,
4248 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4249 * = nullptr>
4250 inline void
4251 compress_start(
4252 VectorStruct &vec,
4253 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4254 const unsigned int channel = 0)
4255 {
4256 if (get_communication_block_size(vec) < vec.n_blocks())
4257 vec.compress(VectorOperation::add);
4258 else
4259 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4260 compress_start(vec.block(i), exchanger, channel + i);
4261 }
4262
4263
4264
4265 // for non-block vectors
4266 template <int dim,
4267 typename VectorStruct,
4268 typename Number,
4269 typename VectorizedArrayType,
4270 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4271 * = nullptr>
4272 inline void
4273 compress_start(
4274 VectorStruct &vec,
4275 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4276 const unsigned int channel = 0)
4277 {
4278 exchanger.compress_start(channel, vec);
4279 }
4280
4281
4282
4283 // for std::vector of vectors
4284 template <int dim,
4285 typename VectorStruct,
4286 typename Number,
4287 typename VectorizedArrayType>
4288 inline void
4289 compress_start(
4290 std::vector<VectorStruct> &vec,
4291 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4292 {
4293 unsigned int component_index = 0;
4294 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4295 {
4296 compress_start(vec[comp], exchanger, component_index);
4297 component_index += n_components(vec[comp]);
4298 }
4299 }
4300
4301
4302
4303 // for std::vector of pointer to vectors
4304 template <int dim,
4305 typename VectorStruct,
4306 typename Number,
4307 typename VectorizedArrayType>
4308 inline void
4309 compress_start(
4310 std::vector<VectorStruct *> &vec,
4311 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4312 {
4313 unsigned int component_index = 0;
4314 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4315 {
4316 compress_start(*vec[comp], exchanger, component_index);
4317 component_index += n_components(*vec[comp]);
4318 }
4319 }
4320
4321
4322
4323 //
4324 // compress_finish
4325 //
4326
4327 // for block vectors
4328 template <int dim,
4329 typename VectorStruct,
4330 typename Number,
4331 typename VectorizedArrayType,
4332 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4333 * = nullptr>
4334 inline void
4335 compress_finish(
4336 VectorStruct &vec,
4337 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4338 const unsigned int channel = 0)
4339 {
4340 if (get_communication_block_size(vec) < vec.n_blocks())
4341 {
4342 // do nothing, everything has already been completed in the _start()
4343 // call
4344 }
4345 else
4346 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4347 compress_finish(vec.block(i), exchanger, channel + i);
4348 }
4349
4350
4351
4352 // for non-block vectors
4353 template <int dim,
4354 typename VectorStruct,
4355 typename Number,
4356 typename VectorizedArrayType,
4357 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4358 * = nullptr>
4359 inline void
4360 compress_finish(
4361 VectorStruct &vec,
4362 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4363 const unsigned int channel = 0)
4364 {
4365 exchanger.compress_finish(channel, vec);
4366 }
4367
4368
4369
4370 // for std::vector of vectors
4371 template <int dim,
4372 typename VectorStruct,
4373 typename Number,
4374 typename VectorizedArrayType>
4375 inline void
4376 compress_finish(
4377 std::vector<VectorStruct> &vec,
4378 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4379 {
4380 unsigned int component_index = 0;
4381 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4382 {
4383 compress_finish(vec[comp], exchanger, component_index);
4384 component_index += n_components(vec[comp]);
4385 }
4386 }
4387
4388
4389
4390 // for std::vector of pointer to vectors
4391 template <int dim,
4392 typename VectorStruct,
4393 typename Number,
4394 typename VectorizedArrayType>
4395 inline void
4396 compress_finish(
4397 std::vector<VectorStruct *> &vec,
4398 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4399 {
4400 unsigned int component_index = 0;
4401 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4402 {
4403 compress_finish(*vec[comp], exchanger, component_index);
4404 component_index += n_components(*vec[comp]);
4405 }
4406 }
4407
4408
4409
4410 //
4411 // reset_ghost_values:
4412 //
4413 // if the input vector did not have ghosts imported, clear them here again
4414 // in order to avoid subsequent operations e.g. in linear solvers to work
4415 // with ghosts all the time
4416 //
4417
4418 // for block vectors
4419 template <int dim,
4420 typename VectorStruct,
4421 typename Number,
4422 typename VectorizedArrayType,
4423 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4424 * = nullptr>
4425 inline void
4426 reset_ghost_values(
4427 const VectorStruct &vec,
4428 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4429 {
4430 // return immediately if there is nothing to do.
4431 if (exchanger.ghosts_were_set == true)
4432 return;
4433
4434 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4435 reset_ghost_values(vec.block(i), exchanger);
4436 }
4437
4438
4439
4440 // for non-block vectors
4441 template <int dim,
4442 typename VectorStruct,
4443 typename Number,
4444 typename VectorizedArrayType,
4445 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4446 * = nullptr>
4447 inline void
4448 reset_ghost_values(
4449 const VectorStruct &vec,
4450 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4451 {
4452 exchanger.reset_ghost_values(vec);
4453 }
4454
4455
4456
4457 // for std::vector of vectors
4458 template <int dim,
4459 typename VectorStruct,
4460 typename Number,
4461 typename VectorizedArrayType>
4462 inline void
4463 reset_ghost_values(
4464 const std::vector<VectorStruct> &vec,
4465 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4466 {
4467 // return immediately if there is nothing to do.
4468 if (exchanger.ghosts_were_set == true)
4469 return;
4470
4471 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4472 reset_ghost_values(vec[comp], exchanger);
4473 }
4474
4475
4476
4477 // for std::vector of pointer to vectors
4478 template <int dim,
4479 typename VectorStruct,
4480 typename Number,
4481 typename VectorizedArrayType>
4482 inline void
4483 reset_ghost_values(
4484 const std::vector<VectorStruct *> &vec,
4485 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4486 {
4487 // return immediately if there is nothing to do.
4488 if (exchanger.ghosts_were_set == true)
4489 return;
4490
4491 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4492 reset_ghost_values(*vec[comp], exchanger);
4493 }
4494
4495
4496
4497 //
4498 // zero_vector_region
4499 //
4500
4501 // for block vectors
4502 template <int dim,
4503 typename VectorStruct,
4504 typename Number,
4505 typename VectorizedArrayType,
4506 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4507 * = nullptr>
4508 inline void
4509 zero_vector_region(
4510 const unsigned int range_index,
4511 VectorStruct &vec,
4512 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4513 {
4514 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4515 exchanger.zero_vector_region(range_index, vec.block(i));
4516 }
4517
4518
4519
4520 // for non-block vectors
4521 template <int dim,
4522 typename VectorStruct,
4523 typename Number,
4524 typename VectorizedArrayType,
4525 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4526 * = nullptr>
4527 inline void
4528 zero_vector_region(
4529 const unsigned int range_index,
4530 VectorStruct &vec,
4531 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4532 {
4533 exchanger.zero_vector_region(range_index, vec);
4534 }
4535
4536
4537
4538 // for std::vector of vectors
4539 template <int dim,
4540 typename VectorStruct,
4541 typename Number,
4542 typename VectorizedArrayType>
4543 inline void
4544 zero_vector_region(
4545 const unsigned int range_index,
4546 std::vector<VectorStruct> &vec,
4547 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4548 {
4549 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4550 zero_vector_region(range_index, vec[comp], exchanger);
4551 }
4552
4553
4554
4555 // for std::vector of pointers to vectors
4556 template <int dim,
4557 typename VectorStruct,
4558 typename Number,
4559 typename VectorizedArrayType>
4560 inline void
4561 zero_vector_region(
4562 const unsigned int range_index,
4563 std::vector<VectorStruct *> &vec,
4564 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4565 {
4566 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4567 zero_vector_region(range_index, *vec[comp], exchanger);
4568 }
4569
4570
4571
4572 // Apply a unit matrix operation to constrained DoFs: Default cases where we
4573 // cannot detect a LinearAlgebra::distributed::Vector, we do not do
4574 // anything, else we apply the constraints as a unit operation
4575 template <typename VectorStruct1, typename VectorStruct2>
4576 inline void
4577 apply_operation_to_constrained_dofs(const std::vector<unsigned int> &,
4578 const VectorStruct1 &,
4579 VectorStruct2 &)
4580 {}
4581
4582 template <typename Number>
4583 inline void
4584 apply_operation_to_constrained_dofs(
4585 const std::vector<unsigned int> &constrained_dofs,
4588 {
4589 for (const unsigned int i : constrained_dofs)
4590 dst.local_element(i) = src.local_element(i);
4591 }
4592
4593
4594 namespace MatrixFreeFunctions
4595 {
4596 // struct to select between a const interface and a non-const interface
4597 // for MFWorker
4598 template <typename, typename, typename, typename, bool>
4599 struct InterfaceSelector
4600 {};
4601
4602 // Version for constant functions
4603 template <typename MF,
4604 typename InVector,
4605 typename OutVector,
4606 typename Container>
4607 struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4608 {
4609 using function_type = void (Container::*)(
4610 const MF &,
4611 OutVector &,
4612 const InVector &,
4613 const std::pair<unsigned int, unsigned int> &) const;
4614 };
4615
4616 // Version for non-constant functions
4617 template <typename MF,
4618 typename InVector,
4619 typename OutVector,
4620 typename Container>
4621 struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4622 {
4623 using function_type =
4624 void (Container::*)(const MF &,
4625 OutVector &,
4626 const InVector &,
4627 const std::pair<unsigned int, unsigned int> &);
4628 };
4629 } // namespace MatrixFreeFunctions
4630
4631
4632
4633 // A implementation class for the worker object that runs the various
4634 // operations we want to perform during the matrix-free loop
4635 template <typename MF,
4636 typename InVector,
4637 typename OutVector,
4638 typename Container,
4639 bool is_constant>
4640 class MFWorker : public MFWorkerInterface
4641 {
4642 public:
4643 // An alias to make the arguments further down more readable
4644 using function_type = typename MatrixFreeFunctions::
4645 InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4646 function_type;
4647
4648 // constructor, binds all the arguments to this class
4649 MFWorker(const MF &matrix_free,
4650 const InVector &src,
4651 OutVector &dst,
4652 const bool zero_dst_vector_setting,
4653 const Container &container,
4654 function_type cell_function,
4655 function_type face_function,
4656 function_type boundary_function,
4657 const typename MF::DataAccessOnFaces src_vector_face_access =
4658 MF::DataAccessOnFaces::none,
4659 const typename MF::DataAccessOnFaces dst_vector_face_access =
4660 MF::DataAccessOnFaces::none,
4661 const std::function<void(const unsigned int, const unsigned int)>
4662 &operation_before_loop = {},
4663 const std::function<void(const unsigned int, const unsigned int)>
4664 &operation_after_loop = {},
4665 const unsigned int dof_handler_index_pre_post = 0)
4666 : matrix_free(matrix_free)
4667 , container(const_cast<Container &>(container))
4668 , cell_function(cell_function)
4669 , face_function(face_function)
4670 , boundary_function(boundary_function)
4671 , src(src)
4672 , dst(dst)
4673 , src_data_exchanger(matrix_free,
4674 src_vector_face_access,
4675 n_components(src))
4676 , dst_data_exchanger(matrix_free,
4677 dst_vector_face_access,
4678 n_components(dst))
4679 , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4680 , zero_dst_vector_setting(zero_dst_vector_setting &&
4681 !src_and_dst_are_same)
4682 , operation_before_loop(operation_before_loop)
4683 , operation_after_loop(operation_after_loop)
4684 , dof_handler_index_pre_post(dof_handler_index_pre_post)
4685 {
4686 Assert(!has_ghost_elements(dst),
4687 ExcMessage("The destination vector passed to the matrix-free "
4688 "loop is ghosted. This is not allowed."));
4689 }
4690
4691 // Runs the cell work. If no function is given, nothing is done
4692 virtual void
4693 cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4694 {
4695 if (cell_function != nullptr && cell_range.second > cell_range.first)
4696 for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4697 {
4698 const auto cell_subrange =
4699 matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4700
4701 if (cell_subrange.second <= cell_subrange.first)
4702 continue;
4703
4704 (container.*
4705 cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4706 }
4707 }
4708
4709 virtual void
4710 cell(const unsigned int range_index) override
4711 {
4712 process_range(cell_function,
4713 matrix_free.get_task_info().cell_partition_data_hp_ptr,
4714 matrix_free.get_task_info().cell_partition_data_hp,
4715 range_index);
4716 }
4717
4718 virtual void
4719 face(const unsigned int range_index) override
4720 {
4721 process_range(face_function,
4722 matrix_free.get_task_info().face_partition_data_hp_ptr,
4723 matrix_free.get_task_info().face_partition_data_hp,
4724 range_index);
4725 }
4726
4727 virtual void
4728 boundary(const unsigned int range_index) override
4729 {
4730 process_range(boundary_function,
4731 matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4732 matrix_free.get_task_info().boundary_partition_data_hp,
4733 range_index);
4734 }
4735
4736 private:
4737 void
4738 process_range(const function_type &fu,
4739 const std::vector<unsigned int> &ptr,
4740 const std::vector<unsigned int> &data,
4741 const unsigned int range_index)
4742 {
4743 if (fu == nullptr)
4744 return;
4745
4746 AssertIndexRange(range_index + 1, ptr.size());
4747 for (unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4748 {
4749 AssertIndexRange(2 * i + 1, data.size());
4750 (container.*fu)(matrix_free,
4751 this->dst,
4752 this->src,
4753 std::make_pair(data[2 * i], data[2 * i + 1]));
4754 }
4755 }
4756
4757 public:
4758 // Starts the communication for the update ghost values operation. We
4759 // cannot call this update if ghost and destination are the same because
4760 // that would introduce spurious entries in the destination (there is also
4761 // the problem that reading from a vector that we also write to is usually
4762 // not intended in case there is overlap, but this is up to the
4763 // application code to decide and we cannot catch this case here).
4764 virtual void
4765 vector_update_ghosts_start() override
4766 {
4767 if (!src_and_dst_are_same)
4768 internal::update_ghost_values_start(src, src_data_exchanger);
4769 }
4770
4771 // Finishes the communication for the update ghost values operation
4772 virtual void
4773 vector_update_ghosts_finish() override
4774 {
4775 if (!src_and_dst_are_same)
4776 internal::update_ghost_values_finish(src, src_data_exchanger);
4777 }
4778
4779 // Starts the communication for the vector compress operation
4780 virtual void
4781 vector_compress_start() override
4782 {
4783 internal::compress_start(dst, dst_data_exchanger);
4784 }
4785
4786 // Finishes the communication for the vector compress operation
4787 virtual void
4788 vector_compress_finish() override
4789 {
4790 internal::compress_finish(dst, dst_data_exchanger);
4791 if (!src_and_dst_are_same)
4792 internal::reset_ghost_values(src, src_data_exchanger);
4793 }
4794
4795 // Zeros the given input vector
4796 virtual void
4797 zero_dst_vector_range(const unsigned int range_index) override
4798 {
4799 if (zero_dst_vector_setting)
4800 internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4801 }
4802
4803 virtual void
4804 cell_loop_pre_range(const unsigned int range_index) override
4805 {
4806 if (operation_before_loop)
4807 {
4809 matrix_free.get_dof_info(dof_handler_index_pre_post);
4810 if (range_index == numbers::invalid_unsigned_int)
4811 {
4812 // Case with threaded loop -> currently no overlap implemented
4814 0U,
4815 dof_info.vector_partitioner->locally_owned_size(),
4816 operation_before_loop,
4818 }
4819 else
4820 {
4821 AssertIndexRange(range_index,
4822 dof_info.cell_loop_pre_list_index.size() - 1);
4823 for (unsigned int id =
4824 dof_info.cell_loop_pre_list_index[range_index];
4825 id != dof_info.cell_loop_pre_list_index[range_index + 1];
4826 ++id)
4827 operation_before_loop(dof_info.cell_loop_pre_list[id].first,
4828 dof_info.cell_loop_pre_list[id].second);
4829 }
4830 }
4831 }
4832
4833 virtual void
4834 cell_loop_post_range(const unsigned int range_index) override
4835 {
4836 if (operation_after_loop)
4837 {
4838 // Run unit matrix operation on constrained dofs if we are at the
4839 // last range
4840 const std::vector<unsigned int> &partition_row_index =
4841 matrix_free.get_task_info().partition_row_index;
4842 if (range_index ==
4843 partition_row_index[partition_row_index.size() - 2] - 1)
4844 apply_operation_to_constrained_dofs(
4845 matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
4846 src,
4847 dst);
4848
4850 matrix_free.get_dof_info(dof_handler_index_pre_post);
4851 if (range_index == numbers::invalid_unsigned_int)
4852 {
4853 // Case with threaded loop -> currently no overlap implemented
4855 0U,
4856 dof_info.vector_partitioner->locally_owned_size(),
4857 operation_after_loop,
4859 }
4860 else
4861 {
4862 AssertIndexRange(range_index,
4863 dof_info.cell_loop_post_list_index.size() - 1);
4864 for (unsigned int id =
4865 dof_info.cell_loop_post_list_index[range_index];
4866 id != dof_info.cell_loop_post_list_index[range_index + 1];
4867 ++id)
4868 operation_after_loop(dof_info.cell_loop_post_list[id].first,
4869 dof_info.cell_loop_post_list[id].second);
4870 }
4871 }
4872 }
4873
4874 private:
4875 const MF &matrix_free;
4876 Container &container;
4877 function_type cell_function;
4878 function_type face_function;
4879 function_type boundary_function;
4880
4881 const InVector &src;
4882 OutVector &dst;
4883 VectorDataExchange<MF::dimension,
4884 typename MF::value_type,
4885 typename MF::vectorized_value_type>
4886 src_data_exchanger;
4887 VectorDataExchange<MF::dimension,
4888 typename MF::value_type,
4889 typename MF::vectorized_value_type>
4890 dst_data_exchanger;
4891 const bool src_and_dst_are_same;
4892 const bool zero_dst_vector_setting;
4893 const std::function<void(const unsigned int, const unsigned int)>
4894 operation_before_loop;
4895 const std::function<void(const unsigned int, const unsigned int)>
4896 operation_after_loop;
4897 const unsigned int dof_handler_index_pre_post;
4898 };
4899
4900
4901
4906 template <class MF, typename InVector, typename OutVector>
4907 struct MFClassWrapper
4908 {
4909 using function_type =
4910 std::function<void(const MF &,
4911 OutVector &,
4912 const InVector &,
4913 const std::pair<unsigned int, unsigned int> &)>;
4914
4915 MFClassWrapper(const function_type cell,
4916 const function_type face,
4917 const function_type boundary)
4918 : cell(cell)
4919 , face(face)
4920 , boundary(boundary)
4921 {}
4922
4923 void
4924 cell_integrator(const MF &mf,
4925 OutVector &dst,
4926 const InVector &src,
4927 const std::pair<unsigned int, unsigned int> &range) const
4928 {
4929 if (cell)
4930 cell(mf, dst, src, range);
4931 }
4932
4933 void
4934 face_integrator(const MF &mf,
4935 OutVector &dst,
4936 const InVector &src,
4937 const std::pair<unsigned int, unsigned int> &range) const
4938 {
4939 if (face)
4940 face(mf, dst, src, range);
4941 }
4942
4943 void
4944 boundary_integrator(
4945 const MF &mf,
4946 OutVector &dst,
4947 const InVector &src,
4948 const std::pair<unsigned int, unsigned int> &range) const
4949 {
4950 if (boundary)
4951 boundary(mf, dst, src, range);
4952 }
4953
4954 const function_type cell;
4955 const function_type face;
4956 const function_type boundary;
4957 };
4958
4959} // end of namespace internal
4960
4961
4962
4963template <int dim, typename Number, typename VectorizedArrayType>
4964template <typename OutVector, typename InVector>
4965inline void
4967 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4968 OutVector &,
4969 const InVector &,
4970 const std::pair<unsigned int, unsigned int> &)>
4971 &cell_operation,
4972 OutVector &dst,
4973 const InVector &src,
4974 const bool zero_dst_vector) const
4975{
4976 using Wrapper =
4977 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4978 InVector,
4979 OutVector>;
4980 Wrapper wrap(cell_operation, nullptr, nullptr);
4981 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4982 InVector,
4983 OutVector,
4984 Wrapper,
4985 true>
4986 worker(*this,
4987 src,
4988 dst,
4989 zero_dst_vector,
4990 wrap,
4991 &Wrapper::cell_integrator,
4992 &Wrapper::face_integrator,
4993 &Wrapper::boundary_integrator);
4994
4995 task_info.loop(worker);
4996}
4997
4998
4999
5000template <int dim, typename Number, typename VectorizedArrayType>
5001template <typename OutVector, typename InVector>
5002inline void
5004 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5005 OutVector &,
5006 const InVector &,
5007 const std::pair<unsigned int, unsigned int> &)>
5008 &cell_operation,
5009 OutVector &dst,
5010 const InVector &src,
5011 const std::function<void(const unsigned int, const unsigned int)>
5012 &operation_before_loop,
5013 const std::function<void(const unsigned int, const unsigned int)>
5014 &operation_after_loop,
5015 const unsigned int dof_handler_index_pre_post) const
5016{
5017 using Wrapper =
5018 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5019 InVector,
5020 OutVector>;
5021 Wrapper wrap(cell_operation, nullptr, nullptr);
5022 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5023 InVector,
5024 OutVector,
5025 Wrapper,
5026 true>
5027 worker(*this,
5028 src,
5029 dst,
5030 false,
5031 wrap,
5032 &Wrapper::cell_integrator,
5033 &Wrapper::face_integrator,
5034 &Wrapper::boundary_integrator,
5035 DataAccessOnFaces::none,
5036 DataAccessOnFaces::none,
5037 operation_before_loop,
5038 operation_after_loop,
5039 dof_handler_index_pre_post);
5040
5041 task_info.loop(worker);
5042}
5043
5044
5045
5046template <int dim, typename Number, typename VectorizedArrayType>
5047template <typename OutVector, typename InVector>
5048inline void
5050 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5051 OutVector &,
5052 const InVector &,
5053 const std::pair<unsigned int, unsigned int> &)>
5054 &cell_operation,
5055 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5056 OutVector &,
5057 const InVector &,
5058 const std::pair<unsigned int, unsigned int> &)>
5059 &face_operation,
5060 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5061 OutVector &,
5062 const InVector &,
5063 const std::pair<unsigned int, unsigned int> &)>
5064 &boundary_operation,
5065 OutVector &dst,
5066 const InVector &src,
5067 const bool zero_dst_vector,
5068 const DataAccessOnFaces dst_vector_face_access,
5069 const DataAccessOnFaces src_vector_face_access) const
5070{
5071 using Wrapper =
5072 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5073 InVector,
5074 OutVector>;
5075 Wrapper wrap(cell_operation, face_operation, boundary_operation);
5076 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5077 InVector,
5078 OutVector,
5079 Wrapper,
5080 true>
5081 worker(*this,
5082 src,
5083 dst,
5084 zero_dst_vector,
5085 wrap,
5086 &Wrapper::cell_integrator,
5087 &Wrapper::face_integrator,
5088 &Wrapper::boundary_integrator,
5089 src_vector_face_access,
5090 dst_vector_face_access);
5091
5092 task_info.loop(worker);
5093}
5094
5095
5096
5097template <int dim, typename Number, typename VectorizedArrayType>
5098template <typename CLASS, typename OutVector, typename InVector>
5099inline void
5101 void (CLASS::*function_pointer)(
5102 const MatrixFree<dim, Number, VectorizedArrayType> &,
5103 OutVector &,
5104 const InVector &,
5105 const std::pair<unsigned int, unsigned int> &) const,
5106 const CLASS *owning_class,
5107 OutVector &dst,
5108 const InVector &src,
5109 const bool zero_dst_vector) const
5110{
5111 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5112 InVector,
5113 OutVector,
5114 CLASS,
5115 true>
5116 worker(*this,
5117 src,
5118 dst,
5119 zero_dst_vector,
5120 *owning_class,
5121 function_pointer,
5122 nullptr,
5123 nullptr);
5124 task_info.loop(worker);
5125}
5126
5127
5128
5129template <int dim, typename Number, typename VectorizedArrayType>
5130template <typename CLASS, typename OutVector, typename InVector>
5131inline void
5133 void (CLASS::*function_pointer)(
5134 const MatrixFree<dim, Number, VectorizedArrayType> &,
5135 OutVector &,
5136 const InVector &,
5137 const std::pair<unsigned int, unsigned int> &) const,
5138 const CLASS *owning_class,
5139 OutVector &dst,
5140 const InVector &src,
5141 const std::function<void(const unsigned int, const unsigned int)>
5142 &operation_before_loop,
5143 const std::function<void(const unsigned int, const unsigned int)>
5144 &operation_after_loop,
5145 const unsigned int dof_handler_index_pre_post) const
5146{
5147 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5148 InVector,
5149 OutVector,
5150 CLASS,
5151 true>
5152 worker(*this,
5153 src,
5154 dst,
5155 false,
5156 *owning_class,
5157 function_pointer,
5158 nullptr,
5159 nullptr,
5162 operation_before_loop,
5163 operation_after_loop,
5164 dof_handler_index_pre_post);
5165 task_info.loop(worker);
5166}
5167
5168
5169
5170template <int dim, typename Number, typename VectorizedArrayType>
5171template <typename CLASS, typename OutVector, typename InVector>
5172inline void
5174 void (CLASS::*cell_operation)(
5175 const MatrixFree<dim, Number, VectorizedArrayType> &,
5176 OutVector &,
5177 const InVector &,
5178 const std::pair<unsigned int, unsigned int> &) const,
5179 void (CLASS::*face_operation)(
5180 const MatrixFree<dim, Number, VectorizedArrayType> &,
5181 OutVector &,
5182 const InVector &,
5183 const std::pair<unsigned int, unsigned int> &) const,
5184 void (CLASS::*boundary_operation)(
5185 const MatrixFree<dim, Number, VectorizedArrayType> &,
5186 OutVector &,
5187 const InVector &,
5188 const std::pair<unsigned int, unsigned int> &) const,
5189 const CLASS *owning_class,
5190 OutVector &dst,
5191 const InVector &src,
5192 const bool zero_dst_vector,
5193 const DataAccessOnFaces dst_vector_face_access,
5194 const DataAccessOnFaces src_vector_face_access) const
5195{
5196 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5197 InVector,
5198 OutVector,
5199 CLASS,
5200 true>
5201 worker(*this,
5202 src,
5203 dst,
5204 zero_dst_vector,
5205 *owning_class,
5206 cell_operation,
5207 face_operation,
5208 boundary_operation,
5209 src_vector_face_access,
5210 dst_vector_face_access);
5211 task_info.loop(worker);
5212}
5213
5214
5215
5216template <int dim, typename Number, typename VectorizedArrayType>
5217template <typename CLASS, typename OutVector, typename InVector>
5218inline void
5220 void (CLASS::*function_pointer)(
5221 const MatrixFree<dim, Number, VectorizedArrayType> &,
5222 OutVector &,
5223 const InVector &,
5224 const std::pair<unsigned int, unsigned int> &),
5225 CLASS *owning_class,
5226 OutVector &dst,
5227 const InVector &src,
5228 const bool zero_dst_vector) const
5229{
5230 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5231 InVector,
5232 OutVector,
5233 CLASS,
5234 false>
5235 worker(*this,
5236 src,
5237 dst,
5238 zero_dst_vector,
5239 *owning_class,
5240 function_pointer,
5241 nullptr,
5242 nullptr);
5243 task_info.loop(worker);
5244}
5245
5246
5247
5248template <int dim, typename Number, typename VectorizedArrayType>
5249template <typename CLASS, typename OutVector, typename InVector>
5250inline void
5252 void (CLASS::*function_pointer)(
5253 const MatrixFree<dim, Number, VectorizedArrayType> &,
5254 OutVector &,
5255 const InVector &,
5256 const std::pair<unsigned int, unsigned int> &),
5257 CLASS *owning_class,
5258 OutVector &dst,
5259 const InVector &src,
5260 const std::function<void(const unsigned int, const unsigned int)>
5261 &operation_before_loop,
5262 const std::function<void(const unsigned int, const unsigned int)>
5263 &operation_after_loop,
5264 const unsigned int dof_handler_index_pre_post) const
5265{
5266 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5267 InVector,
5268 OutVector,
5269 CLASS,
5270 false>
5271 worker(*this,
5272 src,
5273 dst,
5274 false,
5275 *owning_class,
5276 function_pointer,
5277 nullptr,
5278 nullptr,
5281 operation_before_loop,
5282 operation_after_loop,
5283 dof_handler_index_pre_post);
5284 task_info.loop(worker);
5285}
5286
5287
5288
5289template <int dim, typename Number, typename VectorizedArrayType>
5290template <typename CLASS, typename OutVector, typename InVector>
5291inline void
5293 void (CLASS::*cell_operation)(
5294 const MatrixFree<dim, Number, VectorizedArrayType> &,
5295 OutVector &,
5296 const InVector &,
5297 const std::pair<unsigned int, unsigned int> &),
5298 void (CLASS::*face_operation)(
5299 const MatrixFree<dim, Number, VectorizedArrayType> &,
5300 OutVector &,
5301 const InVector &,
5302 const std::pair<unsigned int, unsigned int> &),
5303 void (CLASS::*boundary_operation)(
5304 const MatrixFree<dim, Number, VectorizedArrayType> &,
5305 OutVector &,
5306 const InVector &,
5307 const std::pair<unsigned int, unsigned int> &),
5308 CLASS *owning_class,
5309 OutVector &dst,
5310 const InVector &src,
5311 const bool zero_dst_vector,
5312 const DataAccessOnFaces dst_vector_face_access,
5313 const DataAccessOnFaces src_vector_face_access) const
5314{
5315 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5316 InVector,
5317 OutVector,
5318 CLASS,
5319 false>
5320 worker(*this,
5321 src,
5322 dst,
5323 zero_dst_vector,
5324 *owning_class,
5325 cell_operation,
5326 face_operation,
5327 boundary_operation,
5328 src_vector_face_access,
5329 dst_vector_face_access);
5330 task_info.loop(worker);
5331}
5332
5333
5334
5335template <int dim, typename Number, typename VectorizedArrayType>
5336template <typename OutVector, typename InVector>
5337inline void
5339 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5340 OutVector &,
5341 const InVector &,
5342 const std::pair<unsigned int, unsigned int> &)>
5343 &cell_operation,
5344 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5345 OutVector &,
5346 const InVector &,
5347 const std::pair<unsigned int, unsigned int> &)>
5348 &face_operation,
5349 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5350 OutVector &,
5351 const InVector &,
5352 const std::pair<unsigned int, unsigned int> &)>
5353 &boundary_operation,
5354 OutVector &dst,
5355 const InVector &src,
5356 const std::function<void(const unsigned int, const unsigned int)>
5357 &operation_before_loop,
5358 const std::function<void(const unsigned int, const unsigned int)>
5359 &operation_after_loop,
5360 const unsigned int dof_handler_index_pre_post,
5361 const DataAccessOnFaces dst_vector_face_access,
5362 const DataAccessOnFaces src_vector_face_access) const
5363{
5364 using Wrapper =
5365 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5366 InVector,
5367 OutVector>;
5368 Wrapper wrap(cell_operation, face_operation, boundary_operation);
5369 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5370 InVector,
5371 OutVector,
5372 Wrapper,
5373 true>
5374 worker(*this,
5375 src,
5376 dst,
5377 false,
5378 wrap,
5379 &Wrapper::cell_integrator,
5380 &Wrapper::face_integrator,
5381 &Wrapper::boundary_integrator,
5382 src_vector_face_access,
5383 dst_vector_face_access,
5384 operation_before_loop,
5385 operation_after_loop,
5386 dof_handler_index_pre_post);
5387
5388 task_info.loop(worker);
5389}
5390
5391
5392
5393template <int dim, typename Number, typename VectorizedArrayType>
5394template <typename CLASS, typename OutVector, typename InVector>
5395inline void
5397 void (CLASS::*cell_operation)(const MatrixFree &,
5398 OutVector &,
5399 const InVector &,
5400 const std::pair<unsigned int, unsigned int> &)
5401 const,
5402 void (CLASS::*face_operation)(const MatrixFree &,
5403 OutVector &,
5404 const InVector &,
5405 const std::pair<unsigned int, unsigned int> &)
5406 const,
5407 void (CLASS::*boundary_operation)(
5408 const MatrixFree &,
5409 OutVector &,
5410 const InVector &,
5411 const std::pair<unsigned int, unsigned int> &) const,
5412 const CLASS *owning_class,
5413 OutVector &dst,
5414 const InVector &src,
5415 const std::function<void(const unsigned int, const unsigned int)>
5416 &operation_before_loop,
5417 const std::function<void(const unsigned int, const unsigned int)>
5418 &operation_after_loop,
5419 const unsigned int dof_handler_index_pre_post,
5420 const DataAccessOnFaces dst_vector_face_access,
5421 const DataAccessOnFaces src_vector_face_access) const
5422{
5423 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5424 InVector,
5425 OutVector,
5426 CLASS,
5427 true>
5428 worker(*this,
5429 src,
5430 dst,
5431 false,
5432 *owning_class,
5433 cell_operation,
5434 face_operation,
5435 boundary_operation,
5436 src_vector_face_access,
5437 dst_vector_face_access,
5438 operation_before_loop,
5439 operation_after_loop,
5440 dof_handler_index_pre_post);
5441 task_info.loop(worker);
5442}
5443
5444
5445
5446template <int dim, typename Number, typename VectorizedArrayType>
5447template <typename CLASS, typename OutVector, typename InVector>
5448inline void
5450 void (CLASS::*cell_operation)(const MatrixFree &,
5451 OutVector &,
5452 const InVector &,
5453 const std::pair<unsigned int, unsigned int> &),
5454 void (CLASS::*face_operation)(const MatrixFree &,
5455 OutVector &,
5456 const InVector &,
5457 const std::pair<unsigned int, unsigned int> &),
5458 void (CLASS::*boundary_operation)(
5459 const MatrixFree &,
5460 OutVector &,
5461 const InVector &,
5462 const std::pair<unsigned int, unsigned int> &),
5463 const CLASS *owning_class,
5464 OutVector &dst,
5465 const InVector &src,
5466 const std::function<void(const unsigned int, const unsigned int)>
5467 &operation_before_loop,
5468 const std::function<void(const unsigned int, const unsigned int)>
5469 &operation_after_loop,
5470 const unsigned int dof_handler_index_pre_post,
5471 const DataAccessOnFaces dst_vector_face_access,
5472 const DataAccessOnFaces src_vector_face_access) const
5473{
5474 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5475 InVector,
5476 OutVector,
5477 CLASS,
5478 false>
5479 worker(*this,
5480 src,
5481 dst,
5482 false,
5483 *owning_class,
5484 cell_operation,
5485 face_operation,
5486 boundary_operation,
5487 src_vector_face_access,
5488 dst_vector_face_access,
5489 operation_before_loop,
5490 operation_after_loop,
5491 dof_handler_index_pre_post);
5492 task_info.loop(worker);
5493}
5494
5495
5496
5497template <int dim, typename Number, typename VectorizedArrayType>
5498template <typename CLASS, typename OutVector, typename InVector>
5499inline void
5501 void (CLASS::*function_pointer)(
5502 const MatrixFree<dim, Number, VectorizedArrayType> &,
5503 OutVector &,
5504 const InVector &,
5505 const std::pair<unsigned int, unsigned int> &) const,
5506 const CLASS *owning_class,
5507 OutVector &dst,
5508 const InVector &src,
5509 const bool zero_dst_vector,
5510 const DataAccessOnFaces src_vector_face_access) const
5511{
5512 auto src_vector_face_access_temp = src_vector_face_access;
5513 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5514 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5515 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5516 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5517
5518 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5519 InVector,
5520 OutVector,
5521 CLASS,
5522 true>
5523 worker(*this,
5524 src,
5525 dst,
5526 zero_dst_vector,
5527 *owning_class,
5528 function_pointer,
5529 nullptr,
5530 nullptr,
5531 src_vector_face_access_temp,
5533 task_info.loop(worker);
5534}
5535
5536
5537
5538template <int dim, typename Number, typename VectorizedArrayType>
5539template <typename CLASS, typename OutVector, typename InVector>
5540inline void
5542 void (CLASS::*function_pointer)(
5543 const MatrixFree<dim, Number, VectorizedArrayType> &,
5544 OutVector &,
5545 const InVector &,
5546 const std::pair<unsigned int, unsigned int> &),
5547 CLASS *owning_class,
5548 OutVector &dst,
5549 const InVector &src,
5550 const bool zero_dst_vector,
5551 const DataAccessOnFaces src_vector_face_access) const
5552{
5553 auto src_vector_face_access_temp = src_vector_face_access;
5554 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5555 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5556 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5557 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5558
5559 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5560 InVector,
5561 OutVector,
5562 CLASS,
5563 false>
5564 worker(*this,
5565 src,
5566 dst,
5567 zero_dst_vector,
5568 *owning_class,
5569 function_pointer,
5570 nullptr,
5571 nullptr,
5572 src_vector_face_access_temp,
5574 task_info.loop(worker);
5575}
5576
5577
5578
5579template <int dim, typename Number, typename VectorizedArrayType>
5580template <typename OutVector, typename InVector>
5581inline void
5583 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5584 OutVector &,
5585 const InVector &,
5586 const std::pair<unsigned int, unsigned int> &)>
5587 &cell_operation,
5588 OutVector &dst,
5589 const InVector &src,
5590 const bool zero_dst_vector,
5591 const DataAccessOnFaces src_vector_face_access) const
5592{
5593 auto src_vector_face_access_temp = src_vector_face_access;
5594 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5595 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5596 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5597 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5598
5599 using Wrapper =
5600 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5601 InVector,
5602 OutVector>;
5603 Wrapper wrap(cell_operation, nullptr, nullptr);
5604
5605 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5606 InVector,
5607 OutVector,
5608 Wrapper,
5609 true>
5610 worker(*this,
5611 src,
5612 dst,
5613 zero_dst_vector,
5614 wrap,
5615 &Wrapper::cell_integrator,
5616 &Wrapper::face_integrator,
5617 &Wrapper::boundary_integrator,
5618 src_vector_face_access_temp,
5619 DataAccessOnFaces::none);
5620 task_info.loop(worker);
5621}
5622
5623
5624#endif // ifndef DOXYGEN
5625
5626
5627
5629
5630#endif
iterator begin()
void resize(const size_type new_size)
void reinit(const size_type size, const bool omit_zeroing_entries=false)
Abstract base class for mapping classes.
Definition mapping.h:318
unsigned int n_active_fe_indices() const
std::pair< typename DoFHandler< dim >::cell_iterator, unsigned int > get_face_iterator(const unsigned int face_batch_index, const unsigned int lane_index, const bool interior=true, const unsigned int fe_component=0) const
void loop_cell_centric(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > mapping_info
unsigned int get_cell_active_fe_index(const std::pair< unsigned int, unsigned int > range) const
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
unsigned int n_ghost_cell_batches() const
void 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 > &)> &face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_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
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & 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
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
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > > shape_info
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 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
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
bool at_irregular_cell(const unsigned int cell_batch_index) const
~MatrixFree() override=default
const AffineConstraints< Number > & get_affine_constraints(const unsigned int dof_handler_index=0) const
void copy_from(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free_base)
void initialize_face_data_vector(AlignedVector< T > &vec) const
unsigned int get_dofs_per_cell(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int n_constraint_pool_entries() const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
const Number * constraint_pool_begin(const unsigned int pool_index) const
AlignedVector< Number > * acquire_scratch_data_non_threadsafe() const
const IndexSet & get_locally_owned_set(const unsigned int dof_handler_index=0) const
unsigned int get_dofs_per_face(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
internal::MatrixFreeFunctions::TaskInfo task_info
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*boundary_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
bool mapping_is_initialized
void loop_cell_centric(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void initialize_dof_vector(LinearAlgebra::distributed::Vector< Number2 > &vec, const unsigned int dof_handler_index=0) const
DoFHandler< dim >::cell_iterator get_cell_iterator(const unsigned int cell_batch_index, const unsigned int lane_index, const unsigned int dof_handler_index=0) const
MatrixFree(const MatrixFree< dim, Number, VectorizedArrayType > &other)
std::pair< unsigned int, unsigned int > create_cell_subrange_hp(const std::pair< unsigned int, unsigned int > &range, const unsigned int fe_degree, const unsigned int dof_handler_index=0) const
std::pair< unsigned int, unsigned int > create_cell_subrange_hp_by_index(const std::pair< unsigned int, unsigned int > &range, const unsigned int fe_index, const unsigned int dof_handler_index=0) const
unsigned int mg_level
std::pair< unsigned int, unsigned int > get_face_range_category(const std::pair< unsigned int, unsigned int > face_batch_range) const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
std::vector< std::pair< unsigned int, unsigned int > > cell_level_index
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner(const unsigned int dof_handler_index=0) const
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
std::vector< Number > constraint_pool_data
VectorizedArrayType vectorized_value_type
unsigned int n_cell_batches() const
const IndexSet & get_ghost_set(const unsigned int dof_handler_index=0) const
std::pair< int, int > get_cell_level_and_index(const unsigned int cell_batch_index, const unsigned int lane_index) const
bool indices_initialized() const
const Quadrature< dim - 1 > & get_face_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int get_cell_category(const unsigned int cell_batch_index) const
Threads::ThreadLocalStorage< std::list< std::pair< bool, AlignedVector< VectorizedArrayType > > > > scratch_pad
unsigned int get_matrix_free_cell_index(const typename Triangulation< dim >::cell_iterator &cell) const
Number value_type
std::size_t memory_consumption() const
std::pair< unsigned int, unsigned int > get_face_category(const unsigned int face_batch_index) const
unsigned int n_boundary_face_batches() const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
std::list< std::pair< bool, AlignedVector< Number > > > scratch_pad_non_threadsafe
const Number * constraint_pool_end(const unsigned int pool_index) const
void loop_cell_centric(const std::function< void(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
std::vector< SmartPointer< const DoFHandler< dim > > > dof_handlers
void initialize_indices(const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const AdditionalData &additional_data)
std::vector< SmartPointer< const AffineConstraints< Number > > > affine_constraints
unsigned int n_components() const
unsigned int n_physical_cells() const
void initialize_dof_vector(VectorType &vec, const unsigned int dof_handler_index=0) const
unsigned int get_n_q_points(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
void reinit(const MappingType &mapping, const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
std::vector< unsigned int > mf_cell_indices
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*boundary_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
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 loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*boundary_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 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
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 > &)> &face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_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 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
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*boundary_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
static bool is_supported(const FiniteElement< dim, spacedim > &fe)
static constexpr unsigned int dimension
A class that provides a separate storage location on each thread that accesses the object.
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
unsigned int level
Definition grid_out.cc:4626
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ActiveSelector::cell_iterator cell_iterator
UpdateFlags
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_default
No update.
for(unsigned int j=best_vertex+1;j< vertices.size();++j) if(vertices_to_use[j])
bool job_supports_mpi()
Definition mpi.cc:697
unsigned int minimum_parallel_grain_size
Definition parallel.cc:33
const types::boundary_id invalid_boundary_id
Definition types.h:292
static const unsigned int invalid_unsigned_int
Definition types.h:220
void apply_to_subranges(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, const Function &f, const unsigned int grainsize)
Definition parallel.h:452
const InputIterator OutputIterator const Function & function
Definition parallel.h:168
const Iterator & begin
Definition parallel.h:609
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
unsigned short int fe_index
Definition types.h:59
unsigned int boundary_id
Definition types.h:144
TasksParallelScheme tasks_parallel_scheme
UpdateFlags mapping_update_flags_inner_faces
AdditionalData(const TasksParallelScheme tasks_parallel_scheme=partition_partition, const unsigned int tasks_block_size=0, const UpdateFlags mapping_update_flags=update_gradients|update_JxW_values, const UpdateFlags mapping_update_flags_boundary_faces=update_default, const UpdateFlags mapping_update_flags_inner_faces=update_default, const UpdateFlags mapping_update_flags_faces_by_cells=update_default, const unsigned int mg_level=numbers::invalid_unsigned_int, const bool store_plain_indices=true, const bool initialize_indices=true, const bool initialize_mapping=true, const bool overlap_communication_computation=true, const bool hold_all_faces_to_owned_cells=false, const bool cell_vectorization_categories_strict=false, const bool allow_ghosted_vectors_in_loops=true)
std::vector< unsigned int > cell_vectorization_category
AdditionalData & operator=(const AdditionalData &other)
UpdateFlags mapping_update_flags_boundary_faces
UpdateFlags mapping_update_flags_faces_by_cells
AdditionalData(const AdditionalData &other)
std::vector< unsigned int > cell_loop_pre_list_index
Definition dof_info.h:748
std::vector< unsigned int > cell_loop_post_list_index
Definition dof_info.h:761
std::vector< std::pair< unsigned int, unsigned int > > vector_zero_range_list
Definition dof_info.h:741
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition dof_info.h:592
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_pre_list
Definition dof_info.h:754
std::vector< unsigned int > vector_zero_range_list_index
Definition dof_info.h:736
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_post_list
Definition dof_info.h:767
void loop(MFWorkerInterface &worker) const
Definition task_info.cc:350