deal.II version GIT relicensing-6917-gb466a1809f 2026-10-16 06:30: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: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2012 - 2026 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13
14#ifndef dealii_matrix_free_h
15#define dealii_matrix_free_h
16
17#include <deal.II/base/config.h>
18
26
28
29#include <deal.II/fe/fe.h>
30#include <deal.II/fe/mapping.h>
31
34
39
46
47#include <boost/container/small_vector.hpp>
48
49#include <cstdlib>
50#include <limits>
51#include <list>
52#include <memory>
53
54
56
57
58
110template <int dim,
111 typename Number = double,
112 typename VectorizedArrayType = VectorizedArray<Number>>
114{
115 static_assert(
116 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
117 "Type of Number and of VectorizedArrayType do not match.");
118
119public:
124 using value_type = Number;
125 using vectorized_value_type = VectorizedArrayType;
126
130 static constexpr unsigned int dimension = dim;
131
198 {
203
231
237 const unsigned int tasks_block_size = 0,
243 const unsigned int mg_level = numbers::invalid_unsigned_int,
244 const bool store_plain_indices = true,
245 const bool initialize_indices = true,
246 const bool initialize_mapping = true,
247 const bool overlap_communication_computation = true,
248 const bool hold_all_faces_to_owned_cells = false,
249 const bool cell_vectorization_categories_strict = false,
250 const bool allow_ghosted_vectors_in_loops = true)
266 , store_ghost_cells(false)
267 , communicator_sm(MPI_COMM_SELF)
268 {}
269
296
301 operator=(const AdditionalData &other) = default;
302
340
350 unsigned int tasks_block_size;
351
366
386
406
434
443 unsigned int mg_level;
444
452
463
472
486
495
531 std::vector<unsigned int> cell_vectorization_category;
532
543
555
562
567 };
568
577
582
586 ~MatrixFree() override = default;
587
600 template <typename QuadratureType, typename number2, typename MappingType>
601 void
602 reinit(const MappingType &mapping,
603 const DoFHandler<dim> &dof_handler,
604 const AffineConstraints<number2> &constraint,
605 const QuadratureType &quad,
606 const AdditionalData &additional_data = AdditionalData());
607
629 template <typename QuadratureType, typename number2, typename MappingType>
630 void
631 reinit(const MappingType &mapping,
632 const std::vector<const DoFHandler<dim> *> &dof_handler,
633 const std::vector<const AffineConstraints<number2> *> &constraint,
634 const std::vector<QuadratureType> &quad,
635 const AdditionalData &additional_data = AdditionalData());
636
644 template <typename QuadratureType, typename number2, typename MappingType>
645 void
646 reinit(const MappingType &mapping,
647 const std::vector<const DoFHandler<dim> *> &dof_handler,
648 const std::vector<const AffineConstraints<number2> *> &constraint,
649 const QuadratureType &quad,
650 const AdditionalData &additional_data = AdditionalData());
651
657 void
659 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free_base);
660
670 void
672
676 void
677 update_mapping(const std::shared_ptr<hp::MappingCollection<dim>> &mapping);
678
683 void
685
700 {
707 none,
708
719 values,
720
730
741 gradients,
742
752
760 };
761
806 template <typename OutVector, typename InVector>
807 void
808 cell_loop(const std::function<void(
810 OutVector &,
811 const InVector &,
812 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
813 OutVector &dst,
814 const InVector &src,
815 const bool zero_dst_vector = false) const;
816
863 template <typename CLASS, typename OutVector, typename InVector>
864 void
865 cell_loop(void (CLASS::*cell_operation)(
866 const MatrixFree &,
867 OutVector &,
868 const InVector &,
869 const std::pair<unsigned int, unsigned int> &) const,
870 const CLASS *owning_class,
871 OutVector &dst,
872 const InVector &src,
873 const bool zero_dst_vector = false) const;
874
878 template <typename CLASS, typename OutVector, typename InVector>
879 void
880 cell_loop(void (CLASS::*cell_operation)(
881 const MatrixFree &,
882 OutVector &,
883 const InVector &,
884 const std::pair<unsigned int, unsigned int> &),
885 CLASS *owning_class,
886 OutVector &dst,
887 const InVector &src,
888 const bool zero_dst_vector = false) const;
889
974 template <typename CLASS, typename OutVector, typename InVector>
975 void
976 cell_loop(void (CLASS::*cell_operation)(
977 const MatrixFree &,
978 OutVector &,
979 const InVector &,
980 const std::pair<unsigned int, unsigned int> &) const,
981 const CLASS *owning_class,
982 OutVector &dst,
983 const InVector &src,
984 const std::function<void(const unsigned int, const unsigned int)>
985 &operation_before_loop,
986 const std::function<void(const unsigned int, const unsigned int)>
987 &operation_after_loop,
988 const unsigned int dof_handler_index_pre_post = 0) const;
989
993 template <typename CLASS, typename OutVector, typename InVector>
994 void
995 cell_loop(void (CLASS::*cell_operation)(
996 const MatrixFree &,
997 OutVector &,
998 const InVector &,
999 const std::pair<unsigned int, unsigned int> &),
1000 CLASS *owning_class,
1001 OutVector &dst,
1002 const InVector &src,
1003 const std::function<void(const unsigned int, const unsigned int)>
1004 &operation_before_loop,
1005 const std::function<void(const unsigned int, const unsigned int)>
1006 &operation_after_loop,
1007 const unsigned int dof_handler_index_pre_post = 0) const;
1008
1013 template <typename OutVector, typename InVector>
1014 void
1015 cell_loop(const std::function<void(
1017 OutVector &,
1018 const InVector &,
1019 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1020 OutVector &dst,
1021 const InVector &src,
1022 const std::function<void(const unsigned int, const unsigned int)>
1023 &operation_before_loop,
1024 const std::function<void(const unsigned int, const unsigned int)>
1025 &operation_after_loop,
1026 const unsigned int dof_handler_index_pre_post = 0) const;
1027
1103 template <typename OutVector, typename InVector>
1104 void
1106 const std::function<
1108 OutVector &,
1109 const InVector &,
1110 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1111 const std::function<void(
1113 OutVector &,
1114 const InVector &,
1115 const std::pair<unsigned int, unsigned int> &)> &inner_face_operation,
1116 const std::function<void(
1118 OutVector &,
1119 const InVector &,
1120 const std::pair<unsigned int, unsigned int> &)> &boundary_face_operation,
1121 OutVector &dst,
1122 const InVector &src,
1123 const bool zero_dst_vector = false,
1124 const DataAccessOnFaces dst_vector_face_access =
1126 const DataAccessOnFaces src_vector_face_access =
1128
1217 template <typename CLASS, typename OutVector, typename InVector>
1218 void
1219 loop(void (CLASS::*cell_operation)(
1220 const MatrixFree &,
1221 OutVector &,
1222 const InVector &,
1223 const std::pair<unsigned int, unsigned int> &) const,
1224 void (CLASS::*inner_face_operation)(
1225 const MatrixFree &,
1226 OutVector &,
1227 const InVector &,
1228 const std::pair<unsigned int, unsigned int> &) const,
1229 void (CLASS::*boundary_face_operation)(
1230 const MatrixFree &,
1231 OutVector &,
1232 const InVector &,
1233 const std::pair<unsigned int, unsigned int> &) const,
1234 const CLASS *owning_class,
1235 OutVector &dst,
1236 const InVector &src,
1237 const bool zero_dst_vector = false,
1238 const DataAccessOnFaces dst_vector_face_access =
1240 const DataAccessOnFaces src_vector_face_access =
1242
1246 template <typename CLASS, typename OutVector, typename InVector>
1247 void
1248 loop(void (CLASS::*cell_operation)(
1249 const MatrixFree &,
1250 OutVector &,
1251 const InVector &,
1252 const std::pair<unsigned int, unsigned int> &),
1253 void (CLASS::*inner_face_operation)(
1254 const MatrixFree &,
1255 OutVector &,
1256 const InVector &,
1257 const std::pair<unsigned int, unsigned int> &),
1258 void (CLASS::*boundary_face_operation)(
1259 const MatrixFree &,
1260 OutVector &,
1261 const InVector &,
1262 const std::pair<unsigned int, unsigned int> &),
1263 CLASS *owning_class,
1264 OutVector &dst,
1265 const InVector &src,
1266 const bool zero_dst_vector = false,
1267 const DataAccessOnFaces dst_vector_face_access =
1269 const DataAccessOnFaces src_vector_face_access =
1271
1393 template <typename CLASS, typename OutVector, typename InVector>
1394 void
1395 loop(void (CLASS::*cell_operation)(
1396 const MatrixFree &,
1397 OutVector &,
1398 const InVector &,
1399 const std::pair<unsigned int, unsigned int> &) const,
1400 void (CLASS::*inner_face_operation)(
1401 const MatrixFree &,
1402 OutVector &,
1403 const InVector &,
1404 const std::pair<unsigned int, unsigned int> &) const,
1405 void (CLASS::*boundary_face_operation)(
1406 const MatrixFree &,
1407 OutVector &,
1408 const InVector &,
1409 const std::pair<unsigned int, unsigned int> &) const,
1410 const CLASS *owning_class,
1411 OutVector &dst,
1412 const InVector &src,
1413 const std::function<void(const unsigned int, const unsigned int)>
1414 &operation_before_loop,
1415 const std::function<void(const unsigned int, const unsigned int)>
1416 &operation_after_loop,
1417 const unsigned int dof_handler_index_pre_post = 0,
1418 const DataAccessOnFaces dst_vector_face_access =
1420 const DataAccessOnFaces src_vector_face_access =
1422
1426 template <typename CLASS, typename OutVector, typename InVector>
1427 void
1428 loop(void (CLASS::*cell_operation)(
1429 const MatrixFree &,
1430 OutVector &,
1431 const InVector &,
1432 const std::pair<unsigned int, unsigned int> &),
1433 void (CLASS::*inner_face_operation)(
1434 const MatrixFree &,
1435 OutVector &,
1436 const InVector &,
1437 const std::pair<unsigned int, unsigned int> &),
1438 void (CLASS::*boundary_face_operation)(
1439 const MatrixFree &,
1440 OutVector &,
1441 const InVector &,
1442 const std::pair<unsigned int, unsigned int> &),
1443 const CLASS *owning_class,
1444 OutVector &dst,
1445 const InVector &src,
1446 const std::function<void(const unsigned int, const unsigned int)>
1447 &operation_before_loop,
1448 const std::function<void(const unsigned int, const unsigned int)>
1449 &operation_after_loop,
1450 const unsigned int dof_handler_index_pre_post = 0,
1451 const DataAccessOnFaces dst_vector_face_access =
1453 const DataAccessOnFaces src_vector_face_access =
1455
1461 template <typename OutVector, typename InVector>
1462 void
1464 const std::function<
1466 OutVector &,
1467 const InVector &,
1468 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1469 const std::function<void(
1471 OutVector &,
1472 const InVector &,
1473 const std::pair<unsigned int, unsigned int> &)> &inner_face_operation,
1474 const std::function<void(
1476 OutVector &,
1477 const InVector &,
1478 const std::pair<unsigned int, unsigned int> &)> &boundary_face_operation,
1479 OutVector &dst,
1480 const InVector &src,
1481 const std::function<void(const unsigned int, const unsigned int)>
1482 &operation_before_loop,
1483 const std::function<void(const unsigned int, const unsigned int)>
1484 &operation_after_loop,
1485 const unsigned int dof_handler_index_pre_post = 0,
1486 const DataAccessOnFaces dst_vector_face_access =
1488 const DataAccessOnFaces src_vector_face_access =
1490
1555 template <typename CLASS, typename OutVector, typename InVector>
1556 void
1557 loop_cell_centric(void (CLASS::*cell_operation)(
1558 const MatrixFree &,
1559 OutVector &,
1560 const InVector &,
1561 const std::pair<unsigned int, unsigned int> &) const,
1562 const CLASS *owning_class,
1563 OutVector &dst,
1564 const InVector &src,
1565 const bool zero_dst_vector = false,
1566 const DataAccessOnFaces src_vector_face_access =
1568
1572 template <typename CLASS, typename OutVector, typename InVector>
1573 void
1574 loop_cell_centric(void (CLASS::*cell_operation)(
1575 const MatrixFree &,
1576 OutVector &,
1577 const InVector &,
1578 const std::pair<unsigned int, unsigned int> &),
1579 CLASS *owning_class,
1580 OutVector &dst,
1581 const InVector &src,
1582 const bool zero_dst_vector = false,
1583 const DataAccessOnFaces src_vector_face_access =
1585
1589 template <typename OutVector, typename InVector>
1590 void
1592 const std::function<void(const MatrixFree &,
1593 OutVector &,
1594 const InVector &,
1595 const std::pair<unsigned int, unsigned int> &)>
1596 &cell_operation,
1597 OutVector &dst,
1598 const InVector &src,
1599 const bool zero_dst_vector = false,
1600 const DataAccessOnFaces src_vector_face_access =
1602
1610 std::pair<unsigned int, unsigned int>
1611 create_cell_subrange_hp(const std::pair<unsigned int, unsigned int> &range,
1612 const unsigned int fe_degree,
1613 const unsigned int dof_handler_index = 0) const;
1614
1621 std::pair<unsigned int, unsigned int>
1623 const std::pair<unsigned int, unsigned int> &range,
1624 const unsigned int fe_index,
1625 const unsigned int dof_handler_index = 0) const;
1626
1630 unsigned int
1632
1636 unsigned int
1638 const std::pair<unsigned int, unsigned int> range,
1639 const unsigned int dof_handler_index = numbers::invalid_unsigned_int) const;
1640
1644 unsigned int
1646 const std::pair<unsigned int, unsigned int> range,
1647 const bool is_interior_face = true,
1648 const unsigned int dof_handler_index = numbers::invalid_unsigned_int) const;
1649
1661 template <typename T>
1662 void
1664
1670 template <typename T>
1671 void
1673
1687 template <typename VectorType>
1688 void
1689 initialize_dof_vector(VectorType &vec,
1690 const unsigned int dof_handler_index = 0) const;
1691
1709 template <typename Number2, typename MemorySpace>
1710 void
1713 const unsigned int dof_handler_index = 0) const;
1714
1723 template <typename Number2, typename MemorySpace>
1724 void
1727
1738 const std::shared_ptr<const Utilities::MPI::Partitioner> &
1739 get_vector_partitioner(const unsigned int dof_handler_index = 0) const;
1740
1744 const IndexSet &
1745 get_locally_owned_set(const unsigned int dof_handler_index = 0) const;
1746
1750 const IndexSet &
1751 get_ghost_set(const unsigned int dof_handler_index = 0) const;
1752
1762 const std::vector<unsigned int> &
1763 get_constrained_dofs(const unsigned int dof_handler_index = 0) const;
1764
1775 void
1776 renumber_dofs(std::vector<types::global_dof_index> &renumbering,
1777 const unsigned int dof_handler_index = 0);
1778
1788 template <int spacedim>
1789 static bool
1791
1795 unsigned int
1797
1802 unsigned int
1803 n_base_elements(const unsigned int dof_handler_index) const;
1804
1812 unsigned int
1814
1824 unsigned int
1826
1833 unsigned int
1835
1845 unsigned int
1847
1858 unsigned int
1860
1865 unsigned int
1867
1879 get_boundary_id(const unsigned int face_batch_index) const;
1880
1885 std::array<types::boundary_id, VectorizedArrayType::size()>
1886 get_faces_by_cells_boundary_id(const unsigned int cell_batch_index,
1887 const unsigned int face_number) const;
1888
1893 const DoFHandler<dim> &
1894 get_dof_handler(const unsigned int dof_handler_index = 0) const;
1895
1903 get_affine_constraints(const unsigned int dof_handler_index = 0) const;
1904
1918 get_cell_iterator(const unsigned int cell_batch_index,
1919 const unsigned int lane_index,
1920 const unsigned int dof_handler_index = 0) const;
1921
1927 std::pair<int, int>
1928 get_cell_level_and_index(const unsigned int cell_batch_index,
1929 const unsigned int lane_index) const;
1930
1937 unsigned int
1939 const typename Triangulation<dim>::cell_iterator &cell) const;
1940
1953 std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>
1954 get_face_iterator(const unsigned int face_batch_index,
1955 const unsigned int lane_index,
1956 const bool interior = true,
1957 const unsigned int fe_component = 0) const;
1958
1971 bool
1972 at_irregular_cell(const unsigned int cell_batch_index) const;
1973
1983 unsigned int
1984 n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const;
1985
1995 unsigned int
1996 n_active_entries_per_face_batch(const unsigned int face_batch_index) const;
1997
2001 unsigned int
2002 get_dofs_per_cell(const unsigned int dof_handler_index = 0,
2003 const unsigned int hp_active_fe_index = 0) const;
2004
2008 unsigned int
2009 get_n_q_points(const unsigned int quad_index = 0,
2010 const unsigned int hp_active_fe_index = 0) const;
2011
2016 unsigned int
2017 get_dofs_per_face(const unsigned int dof_handler_index = 0,
2018 const unsigned int hp_active_fe_index = 0) const;
2019
2024 unsigned int
2025 get_n_q_points_face(const unsigned int quad_index = 0,
2026 const unsigned int hp_active_fe_index = 0) const;
2027
2031 const Quadrature<dim> &
2032 get_quadrature(const unsigned int quad_index = 0,
2033 const unsigned int hp_active_fe_index = 0) const;
2034
2038 const Quadrature<dim - 1> &
2039 get_face_quadrature(const unsigned int quad_index = 0,
2040 const unsigned int hp_active_fe_index = 0) const;
2041
2054 unsigned int
2056 const std::pair<unsigned int, unsigned int> cell_batch_range,
2057 const unsigned int dof_handler_index = numbers::invalid_unsigned_int) const;
2058
2063 std::pair<unsigned int, unsigned int>
2065 const std::pair<unsigned int, unsigned int> face_batch_range,
2066 const unsigned int dof_handler_index = numbers::invalid_unsigned_int) const;
2067
2079 unsigned int
2081 const unsigned int cell_batch_index,
2082 const unsigned int dof_handler_index = numbers::invalid_unsigned_int) const;
2083
2088 std::pair<unsigned int, unsigned int>
2090 const unsigned int face_batch_index,
2091 const unsigned int dof_handler_index = numbers::invalid_unsigned_int) const;
2092
2096 bool
2098
2103 bool
2105
2110 unsigned int
2112
2117 std::size_t
2119
2124 template <typename StreamType>
2125 void
2126 print_memory_consumption(StreamType &out) const;
2127
2132 void
2133 print(std::ostream &out) const;
2134
2148
2149 /*
2150 * Return geometry-dependent information on the cells.
2151 */
2152 const internal::MatrixFreeFunctions::
2153 MappingInfo<dim, Number, VectorizedArrayType> &
2155
2160 get_dof_info(const unsigned int dof_handler_index_component = 0) const;
2161
2165 unsigned int
2167
2172 const Number *
2173 constraint_pool_begin(const unsigned int pool_index) const;
2174
2180 const Number *
2181 constraint_pool_end(const unsigned int pool_index) const;
2182
2187 get_shape_info(const unsigned int dof_handler_index_component = 0,
2188 const unsigned int quad_index = 0,
2189 const unsigned int fe_base_element = 0,
2190 const unsigned int hp_active_fe_index = 0,
2191 const unsigned int hp_active_quad_index = 0) const;
2192
2197 VectorizedArrayType::size()> &
2198 get_face_info(const unsigned int face_batch_index) const;
2199
2200
2208
2224
2228 void
2230
2242
2246 void
2248 const AlignedVector<Number> *memory) const;
2249
2252private:
2257 template <typename number2, int q_dim>
2258 void
2260 const std::shared_ptr<hp::MappingCollection<dim>> &mapping,
2261 const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
2262 const std::vector<const AffineConstraints<number2> *> &constraint,
2263 const std::vector<IndexSet> &locally_owned_set,
2264 const std::vector<hp::QCollection<q_dim>> &quad,
2265 const AdditionalData &additional_data);
2266
2273 template <typename number2>
2274 void
2276 const std::vector<const AffineConstraints<number2> *> &constraint,
2277 const std::vector<IndexSet> &locally_owned_set,
2278 const AdditionalData &additional_data);
2279
2283 void
2285 const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
2286 const AdditionalData &additional_data);
2287
2291 std::vector<ObserverPointer<const DoFHandler<dim>>> dof_handlers;
2292
2299 std::vector<ObserverPointer<const AffineConstraints<Number>>>
2301
2306 std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
2307
2314 std::vector<Number> constraint_pool_data;
2315
2320 std::vector<unsigned int> constraint_pool_row_index;
2321
2328
2333
2340 std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
2341
2346 std::vector<unsigned int> mf_cell_indices;
2347
2355
2362
2367 internal::MatrixFreeFunctions::FaceInfo<VectorizedArrayType::size()>
2369
2374
2379
2388 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2390
2395 mutable std::list<std::pair<bool, AlignedVector<Number>>>
2397
2401 unsigned int mg_level;
2402
2408};
2409
2410
2411
2412/*----------------------- Inline functions ----------------------------------*/
2413
2414#ifndef DOXYGEN
2415
2416
2417
2418template <int dim, typename Number, typename VectorizedArrayType>
2419template <typename T>
2420inline void
2422 AlignedVector<T> &vec) const
2423{
2424 vec.resize(this->n_cell_batches() + this->n_ghost_cell_batches());
2425}
2426
2427
2428
2429template <int dim, typename Number, typename VectorizedArrayType>
2430template <typename T>
2431inline void
2433 AlignedVector<T> &vec) const
2434{
2435 vec.resize(this->n_inner_face_batches() + this->n_boundary_face_batches() +
2436 this->n_ghost_inner_face_batches());
2437}
2438
2439
2440
2441template <int dim, typename Number, typename VectorizedArrayType>
2442template <typename VectorType>
2443inline void
2445 VectorType &vec,
2446 const unsigned int comp) const
2447{
2448 static_assert(IsBlockVector<VectorType>::value == false,
2449 "This function is not supported for block vectors.");
2450
2451 Assert(task_info.n_procs == 1,
2452 ExcMessage("This function can only be used in serial."));
2453
2454 AssertIndexRange(comp, n_components());
2455 vec.reinit(dof_info[comp].vector_partitioner->size());
2456}
2457
2458
2459
2460template <int dim, typename Number, typename VectorizedArrayType>
2461template <typename Number2, typename MemorySpace>
2462inline void
2465 const unsigned int comp) const
2466{
2467 AssertIndexRange(comp, n_components());
2468 vec.reinit(dof_info[comp].vector_partitioner, task_info.communicator_sm);
2469}
2470
2471
2472
2473template <int dim, typename Number, typename VectorizedArrayType>
2474template <typename Number2, typename MemorySpace>
2475inline void
2478{
2479 vec.reinit(n_components());
2480 for (unsigned int c = 0; c < n_components(); ++c)
2481 this->initialize_dof_vector(vec.block(c), c);
2482
2483 vec.collect_sizes();
2484}
2485
2486
2487
2488template <int dim, typename Number, typename VectorizedArrayType>
2489inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2491 const unsigned int comp) const
2492{
2493 AssertIndexRange(comp, n_components());
2494 return dof_info[comp].vector_partitioner;
2495}
2496
2497
2498
2499template <int dim, typename Number, typename VectorizedArrayType>
2500inline const std::vector<unsigned int> &
2502 const unsigned int comp) const
2503{
2504 AssertIndexRange(comp, n_components());
2505 return dof_info[comp].constrained_dofs;
2506}
2507
2508
2509
2510template <int dim, typename Number, typename VectorizedArrayType>
2511inline unsigned int
2513{
2514 AssertDimension(dof_handlers.size(), dof_info.size());
2515 return dof_handlers.size();
2516}
2517
2518
2519
2520template <int dim, typename Number, typename VectorizedArrayType>
2521inline unsigned int
2523 const unsigned int dof_handler_index) const
2524{
2525 AssertDimension(dof_handlers.size(), dof_info.size());
2526 AssertIndexRange(dof_handler_index, dof_handlers.size());
2527 return dof_handlers[dof_handler_index]->get_fe().n_base_elements();
2528}
2529
2530
2531
2532template <int dim, typename Number, typename VectorizedArrayType>
2535{
2536 return task_info;
2537}
2538
2539
2540
2541template <int dim, typename Number, typename VectorizedArrayType>
2542inline unsigned int
2544{
2545 return task_info.n_active_cells;
2546}
2547
2548
2549
2550template <int dim, typename Number, typename VectorizedArrayType>
2551inline unsigned int
2553{
2554 return *(task_info.cell_partition_data.end() - 2);
2555}
2556
2557
2558
2559template <int dim, typename Number, typename VectorizedArrayType>
2560inline unsigned int
2562{
2563 return *(task_info.cell_partition_data.end() - 1) -
2564 *(task_info.cell_partition_data.end() - 2);
2565}
2566
2567
2568
2569template <int dim, typename Number, typename VectorizedArrayType>
2570inline unsigned int
2572{
2573 if (task_info.face_partition_data.empty())
2574 return 0;
2575 return task_info.face_partition_data.back();
2576}
2577
2578
2579
2580template <int dim, typename Number, typename VectorizedArrayType>
2581inline unsigned int
2583{
2584 if (task_info.face_partition_data.empty())
2585 return 0;
2586 return task_info.boundary_partition_data.back() -
2587 task_info.face_partition_data.back();
2588}
2589
2590
2591
2592template <int dim, typename Number, typename VectorizedArrayType>
2593inline unsigned int
2595{
2596 if (task_info.face_partition_data.empty())
2597 return 0;
2598 return face_info.faces.size() - task_info.boundary_partition_data.back();
2599}
2600
2601
2602
2603template <int dim, typename Number, typename VectorizedArrayType>
2604inline types::boundary_id
2606 const unsigned int face_batch_index) const
2607{
2608 Assert(face_batch_index >= task_info.boundary_partition_data[0] &&
2609 face_batch_index < task_info.boundary_partition_data.back(),
2610 ExcIndexRange(face_batch_index,
2611 task_info.boundary_partition_data[0],
2612 task_info.boundary_partition_data.back()));
2613 return types::boundary_id(face_info.faces[face_batch_index].exterior_face_no);
2614}
2615
2616
2617
2618template <int dim, typename Number, typename VectorizedArrayType>
2619inline std::array<types::boundary_id, VectorizedArrayType::size()>
2621 const unsigned int cell_batch_index,
2622 const unsigned int face_number) const
2623{
2624 AssertIndexRange(cell_batch_index, n_cell_batches());
2625 AssertIndexRange(face_number, ReferenceCells::max_n_faces<dim>());
2626 Assert(face_info.cell_and_face_boundary_id.size(0) >= n_cell_batches(),
2628 std::array<types::boundary_id, VectorizedArrayType::size()> result;
2629 result.fill(numbers::invalid_boundary_id);
2630 for (unsigned int v = 0;
2631 v < n_active_entries_per_cell_batch(cell_batch_index);
2632 ++v)
2633 result[v] =
2634 face_info.cell_and_face_boundary_id(cell_batch_index, face_number, v);
2635 return result;
2636}
2637
2638
2639
2640template <int dim, typename Number, typename VectorizedArrayType>
2641inline const internal::MatrixFreeFunctions::
2642 MappingInfo<dim, Number, VectorizedArrayType> &
2644{
2645 return mapping_info;
2646}
2647
2648
2649
2650template <int dim, typename Number, typename VectorizedArrayType>
2653 const unsigned int dof_index) const
2654{
2655 AssertIndexRange(dof_index, n_components());
2656 return dof_info[dof_index];
2657}
2658
2659
2660
2661template <int dim, typename Number, typename VectorizedArrayType>
2662inline unsigned int
2664{
2665 return constraint_pool_row_index.size() - 1;
2666}
2667
2668
2669
2670template <int dim, typename Number, typename VectorizedArrayType>
2671inline const Number *
2673 const unsigned int row) const
2674{
2675 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2676 return constraint_pool_data.empty() ?
2677 nullptr :
2678 constraint_pool_data.data() + constraint_pool_row_index[row];
2679}
2680
2681
2682
2683template <int dim, typename Number, typename VectorizedArrayType>
2684inline const Number *
2686 const unsigned int row) const
2687{
2688 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2689 return constraint_pool_data.empty() ?
2690 nullptr :
2691 constraint_pool_data.data() + constraint_pool_row_index[row + 1];
2692}
2693
2694
2695
2696template <int dim, typename Number, typename VectorizedArrayType>
2697inline std::pair<unsigned int, unsigned int>
2699 const std::pair<unsigned int, unsigned int> &range,
2700 const unsigned int degree,
2701 const unsigned int dof_handler_component) const
2702{
2703 if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2704 {
2706 dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2708 dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2709 if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2710 return range;
2711 else
2712 return {range.second, range.second};
2713 }
2714
2715 const unsigned int fe_index =
2716 dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2717 if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2718 return {range.second, range.second};
2719 else
2720 return create_cell_subrange_hp_by_index(range,
2721 fe_index,
2722 dof_handler_component);
2723}
2724
2725
2726
2727template <int dim, typename Number, typename VectorizedArrayType>
2728inline bool
2730 const unsigned int cell_batch_index) const
2731{
2732 AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2733 return VectorizedArrayType::size() > 1 &&
2734 cell_level_index[(cell_batch_index + 1) * VectorizedArrayType::size() -
2735 1] == cell_level_index[(cell_batch_index + 1) *
2736 VectorizedArrayType::size() -
2737 2];
2738}
2739
2740
2741
2742template <int dim, typename Number, typename VectorizedArrayType>
2743unsigned int
2745{
2746 return shape_info.size(2);
2747}
2748
2749
2750template <int dim, typename Number, typename VectorizedArrayType>
2751unsigned int
2753 const std::pair<unsigned int, unsigned int> range,
2754 const unsigned int dof_handler_index) const
2755{
2756 const unsigned int dof_handler_index_local =
2757 dof_handler_index == numbers::invalid_unsigned_int ?
2758 first_hp_dof_handler_index :
2759 dof_handler_index;
2760
2761 const auto &fe_indices =
2762 dof_info[dof_handler_index_local].cell_active_fe_index;
2763
2764 if (fe_indices.empty() == true ||
2765 dof_handlers[dof_handler_index_local]->get_fe_collection().size() == 1)
2766 return 0;
2767
2768 const auto index = fe_indices[range.first];
2769
2770 for (unsigned int i = range.first; i < range.second; ++i)
2771 AssertDimension(index, fe_indices[i]);
2772
2773 return index;
2774}
2775
2776
2777
2778template <int dim, typename Number, typename VectorizedArrayType>
2779unsigned int
2781 const std::pair<unsigned int, unsigned int> range,
2782 const bool is_interior_face,
2783 const unsigned int dof_handler_index) const
2784{
2785 const unsigned int dof_handler_index_local =
2786 dof_handler_index == numbers::invalid_unsigned_int ?
2787 first_hp_dof_handler_index :
2788 dof_handler_index;
2789
2790 const auto &fe_indices =
2791 dof_info[dof_handler_index_local].cell_active_fe_index;
2792
2793 if (fe_indices.empty() == true)
2794 return 0;
2795
2796 if (is_interior_face)
2797 {
2798 const unsigned int index =
2799 fe_indices[face_info.faces[range.first].cells_interior[0] /
2800 VectorizedArrayType::size()];
2801
2802 for (unsigned int i = range.first; i < range.second; ++i)
2803 AssertDimension(index,
2804 fe_indices[face_info.faces[i].cells_interior[0] /
2805 VectorizedArrayType::size()]);
2806
2807 return index;
2808 }
2809 else
2810 {
2811 const unsigned int index =
2812 fe_indices[face_info.faces[range.first].cells_exterior[0] /
2813 VectorizedArrayType::size()];
2814
2815 for (unsigned int i = range.first; i < range.second; ++i)
2816 AssertDimension(index,
2817 fe_indices[face_info.faces[i].cells_exterior[0] /
2818 VectorizedArrayType::size()]);
2819
2820 return index;
2821 }
2822}
2823
2824
2825
2826template <int dim, typename Number, typename VectorizedArrayType>
2827inline unsigned int
2829 const unsigned int cell_batch_index) const
2830{
2831 Assert(!dof_info.empty(), ExcNotInitialized());
2832 AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2833 const std::vector<unsigned char> &n_lanes_filled =
2834 dof_info[0].n_vectorization_lanes_filled
2836 AssertIndexRange(cell_batch_index, n_lanes_filled.size());
2837
2838 return n_lanes_filled[cell_batch_index];
2839}
2840
2841
2842
2843template <int dim, typename Number, typename VectorizedArrayType>
2844inline unsigned int
2846 const unsigned int face_batch_index) const
2847{
2848 AssertIndexRange(face_batch_index, face_info.faces.size());
2849 Assert(!dof_info.empty(), ExcNotInitialized());
2850 const std::vector<unsigned char> &n_lanes_filled =
2851 dof_info[0].n_vectorization_lanes_filled
2853 AssertIndexRange(face_batch_index, n_lanes_filled.size());
2854 return n_lanes_filled[face_batch_index];
2855}
2856
2857
2858
2859template <int dim, typename Number, typename VectorizedArrayType>
2860inline unsigned int
2862 const unsigned int dof_handler_index,
2863 const unsigned int active_fe_index) const
2864{
2865 return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2866}
2867
2868
2869
2870template <int dim, typename Number, typename VectorizedArrayType>
2871inline unsigned int
2873 const unsigned int quad_index,
2874 const unsigned int active_fe_index) const
2875{
2876 AssertIndexRange(quad_index, mapping_info.cell_data.size());
2877 return mapping_info.cell_data[quad_index]
2878 .descriptor[active_fe_index]
2879 .n_q_points;
2880}
2881
2882
2883
2884template <int dim, typename Number, typename VectorizedArrayType>
2885inline unsigned int
2887 const unsigned int dof_handler_index,
2888 const unsigned int active_fe_index) const
2889{
2890 return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2891}
2892
2893
2894
2895template <int dim, typename Number, typename VectorizedArrayType>
2896inline unsigned int
2898 const unsigned int quad_index,
2899 const unsigned int active_fe_index) const
2900{
2901 AssertIndexRange(quad_index, mapping_info.face_data.size());
2902 return mapping_info.face_data[quad_index]
2903 .descriptor[active_fe_index]
2904 .n_q_points;
2905}
2906
2907
2908
2909template <int dim, typename Number, typename VectorizedArrayType>
2910inline const IndexSet &
2912 const unsigned int dof_handler_index) const
2913{
2914 return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2915}
2916
2917
2918
2919template <int dim, typename Number, typename VectorizedArrayType>
2920inline const IndexSet &
2922 const unsigned int dof_handler_index) const
2923{
2924 return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2925}
2926
2927
2928
2929template <int dim, typename Number, typename VectorizedArrayType>
2932 const unsigned int dof_handler_index,
2933 const unsigned int index_quad,
2934 const unsigned int index_fe,
2935 const unsigned int active_fe_index,
2936 const unsigned int active_quad_index) const
2937{
2938 AssertIndexRange(dof_handler_index, dof_info.size());
2939 const unsigned int ind =
2940 dof_info[dof_handler_index].global_base_element_offset + index_fe;
2941 AssertIndexRange(ind, shape_info.size(0));
2942 AssertIndexRange(index_quad, shape_info.size(1));
2943 AssertIndexRange(active_fe_index, shape_info.size(2));
2944 AssertIndexRange(active_quad_index, shape_info.size(3));
2945 return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2946}
2947
2948
2949
2950template <int dim, typename Number, typename VectorizedArrayType>
2952 VectorizedArrayType::size()> &
2954 const unsigned int face_batch_index) const
2955{
2956 AssertIndexRange(face_batch_index, face_info.faces.size());
2957 return face_info.faces[face_batch_index];
2958}
2959
2960
2961
2962template <int dim, typename Number, typename VectorizedArrayType>
2963inline const Table<3, unsigned int> &
2965 const
2966{
2967 return face_info.cell_and_face_to_plain_faces;
2968}
2969
2970
2971
2972template <int dim, typename Number, typename VectorizedArrayType>
2973inline const Quadrature<dim> &
2975 const unsigned int quad_index,
2976 const unsigned int active_fe_index) const
2977{
2978 AssertIndexRange(quad_index, mapping_info.cell_data.size());
2979 return mapping_info.cell_data[quad_index]
2980 .descriptor[active_fe_index]
2981 .quadrature;
2982}
2983
2984
2985
2986template <int dim, typename Number, typename VectorizedArrayType>
2987inline const Quadrature<dim - 1> &
2989 const unsigned int quad_index,
2990 const unsigned int active_fe_index) const
2991{
2992 AssertIndexRange(quad_index, mapping_info.face_data.size());
2993 return mapping_info.face_data[quad_index]
2994 .descriptor[active_fe_index]
2995 .quadrature;
2996}
2997
2998
2999
3000template <int dim, typename Number, typename VectorizedArrayType>
3001inline unsigned int
3003 const std::pair<unsigned int, unsigned int> range,
3004 const unsigned int dof_handler_index) const
3005{
3006 auto result = get_cell_category(range.first, dof_handler_index);
3007
3008 for (unsigned int i = range.first; i < range.second; ++i)
3009 result = std::max(result, get_cell_category(i, dof_handler_index));
3010
3011 return result;
3012}
3013
3014
3015
3016template <int dim, typename Number, typename VectorizedArrayType>
3017inline std::pair<unsigned int, unsigned int>
3019 const std::pair<unsigned int, unsigned int> range,
3020 const unsigned int dof_handler_index) const
3021{
3022 auto result = get_face_category(range.first, dof_handler_index);
3023
3024 for (unsigned int i = range.first; i < range.second; ++i)
3025 {
3026 result.first =
3027 std::max(result.first, get_face_category(i, dof_handler_index).first);
3028 result.second =
3029 std::max(result.second, get_face_category(i, dof_handler_index).second);
3030 }
3031
3032 return result;
3033}
3034
3035
3036
3037template <int dim, typename Number, typename VectorizedArrayType>
3038inline unsigned int
3040 const unsigned int cell_batch_index,
3041 const unsigned int dof_handler_index) const
3042{
3043 AssertIndexRange(0, dof_info.size());
3044
3045 const unsigned int dof_handler_index_local =
3046 dof_handler_index == numbers::invalid_unsigned_int ?
3047 first_hp_dof_handler_index :
3048 dof_handler_index;
3049
3051 cell_batch_index,
3052 dof_info[dof_handler_index_local].cell_active_fe_index.size());
3053 if (dof_info[dof_handler_index_local].cell_active_fe_index.empty())
3054 return 0;
3055 else
3056 return dof_info[dof_handler_index_local]
3057 .cell_active_fe_index[cell_batch_index];
3058}
3059
3060
3061
3062template <int dim, typename Number, typename VectorizedArrayType>
3063inline std::pair<unsigned int, unsigned int>
3065 const unsigned int face_batch_index,
3066 const unsigned int dof_handler_index) const
3067{
3068 const unsigned int dof_handler_index_local =
3069 dof_handler_index == numbers::invalid_unsigned_int ?
3070 first_hp_dof_handler_index :
3071 dof_handler_index;
3072
3073 AssertIndexRange(face_batch_index, face_info.faces.size());
3074 if (dof_info[dof_handler_index_local].cell_active_fe_index.empty())
3075 return std::make_pair(0U, 0U);
3076
3077 std::pair<unsigned int, unsigned int> result = std::make_pair(0U, 0U);
3078 for (unsigned int v = 0;
3079 v < VectorizedArrayType::size() &&
3080 face_info.faces[face_batch_index].cells_interior[v] !=
3082 ++v)
3083 result.first =
3084 std::max(result.first,
3085 dof_info[dof_handler_index_local].cell_active_fe_index
3086 [face_info.faces[face_batch_index].cells_interior[v] /
3087 VectorizedArrayType::size()]);
3088 if (face_info.faces[face_batch_index].cells_exterior[0] !=
3090 for (unsigned int v = 0;
3091 v < VectorizedArrayType::size() &&
3092 face_info.faces[face_batch_index].cells_exterior[v] !=
3094 ++v)
3095 result.second =
3096 std::max(result.second,
3097 dof_info[dof_handler_index_local].cell_active_fe_index
3098 [face_info.faces[face_batch_index].cells_exterior[v] /
3099 VectorizedArrayType::size()]);
3100 else
3101 result.second = numbers::invalid_unsigned_int;
3102 return result;
3103}
3104
3105
3106
3107template <int dim, typename Number, typename VectorizedArrayType>
3108inline bool
3110{
3111 return indices_are_initialized;
3112}
3113
3114
3115
3116template <int dim, typename Number, typename VectorizedArrayType>
3117inline bool
3119{
3120 return mapping_is_initialized;
3121}
3122
3123
3124template <int dim, typename Number, typename VectorizedArrayType>
3125inline unsigned int
3127{
3128 return mg_level;
3129}
3130
3131
3132
3133template <int dim, typename Number, typename VectorizedArrayType>
3136{
3137 using list_type =
3138 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3139 list_type &data = scratch_pad.get();
3140 for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
3141 if (it->first == false)
3142 {
3143 it->first = true;
3144 return &it->second;
3145 }
3146 data.emplace_front(true, AlignedVector<VectorizedArrayType>());
3147 return &data.front().second;
3148}
3149
3150
3151
3152template <int dim, typename Number, typename VectorizedArrayType>
3153void
3155 const AlignedVector<VectorizedArrayType> *scratch) const
3156{
3157 using list_type =
3158 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3159 list_type &data = scratch_pad.get();
3160 for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
3161 if (&it->second == scratch)
3162 {
3163 Assert(it->first == true, ExcInternalError());
3164 it->first = false;
3165 return;
3166 }
3167 AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
3168}
3169
3170
3171
3172template <int dim, typename Number, typename VectorizedArrayType>
3176{
3177 for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
3178 scratch_pad_non_threadsafe.begin();
3179 it != scratch_pad_non_threadsafe.end();
3180 ++it)
3181 if (it->first == false)
3182 {
3183 it->first = true;
3184 return &it->second;
3185 }
3186 scratch_pad_non_threadsafe.push_front(
3187 std::make_pair(true, AlignedVector<Number>()));
3188 return &scratch_pad_non_threadsafe.front().second;
3189}
3190
3191
3192
3193template <int dim, typename Number, typename VectorizedArrayType>
3194void
3197 const AlignedVector<Number> *scratch) const
3198{
3199 for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
3200 scratch_pad_non_threadsafe.begin();
3201 it != scratch_pad_non_threadsafe.end();
3202 ++it)
3203 if (&it->second == scratch)
3204 {
3205 Assert(it->first == true, ExcInternalError());
3206 it->first = false;
3207 return;
3208 }
3209 AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
3210}
3211
3212
3213
3214// ------------------------------ reinit functions ---------------------------
3215
3216namespace internal
3217{
3218 namespace MatrixFreeImplementation
3219 {
3220 template <int dim, int spacedim>
3221 inline std::vector<IndexSet>
3222 extract_locally_owned_index_sets(
3223 const std::vector<const DoFHandler<dim, spacedim> *> &dofh,
3224 const unsigned int level)
3225 {
3226 std::vector<IndexSet> locally_owned_set;
3227 locally_owned_set.reserve(dofh.size());
3228 for (unsigned int j = 0; j < dofh.size(); ++j)
3230 locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
3231 else
3232 locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
3233 return locally_owned_set;
3234 }
3235 } // namespace MatrixFreeImplementation
3236} // namespace internal
3237
3238
3239
3240template <int dim, typename Number, typename VectorizedArrayType>
3241template <typename QuadratureType, typename number2, typename MappingType>
3242void
3244 const MappingType &mapping,
3245 const DoFHandler<dim> &dof_handler,
3246 const AffineConstraints<number2> &constraints_in,
3247 const QuadratureType &quad,
3249 &additional_data)
3250{
3251 std::vector<const DoFHandler<dim, dim> *> dof_handlers;
3252 std::vector<const AffineConstraints<number2> *> constraints;
3253
3254 dof_handlers.push_back(&dof_handler);
3255 constraints.push_back(&constraints_in);
3256
3257 std::vector<IndexSet> locally_owned_sets =
3258 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3259 dof_handlers, additional_data.mg_level);
3260
3261 std::vector<hp::QCollection<dim>> quad_hp;
3262 quad_hp.emplace_back(quad);
3263
3264 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3265 dof_handlers,
3266 constraints,
3267 locally_owned_sets,
3268 quad_hp,
3269 additional_data);
3270}
3271
3272
3273
3274template <int dim, typename Number, typename VectorizedArrayType>
3275template <typename QuadratureType, typename number2, typename MappingType>
3276void
3278 const MappingType &mapping,
3279 const std::vector<const DoFHandler<dim> *> &dof_handler,
3280 const std::vector<const AffineConstraints<number2> *> &constraint,
3281 const QuadratureType &quad,
3283 &additional_data)
3284{
3285 std::vector<IndexSet> locally_owned_set =
3286 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3287 dof_handler, additional_data.mg_level);
3288 std::vector<hp::QCollection<dim>> quad_hp;
3289 quad_hp.emplace_back(quad);
3290
3291 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3292 dof_handler,
3293 constraint,
3294 locally_owned_set,
3295 quad_hp,
3296 additional_data);
3297}
3298
3299
3300
3301template <int dim, typename Number, typename VectorizedArrayType>
3302template <typename QuadratureType, typename number2, typename MappingType>
3303void
3305 const MappingType &mapping,
3306 const std::vector<const DoFHandler<dim> *> &dof_handler,
3307 const std::vector<const AffineConstraints<number2> *> &constraint,
3308 const std::vector<QuadratureType> &quad,
3310 &additional_data)
3311{
3312 std::vector<IndexSet> locally_owned_set =
3313 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3314 dof_handler, additional_data.mg_level);
3315 std::vector<hp::QCollection<dim>> quad_hp;
3316 for (unsigned int q = 0; q < quad.size(); ++q)
3317 quad_hp.emplace_back(quad[q]);
3318
3319 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3320 dof_handler,
3321 constraint,
3322 locally_owned_set,
3323 quad_hp,
3324 additional_data);
3325}
3326
3327
3328
3329// ------------------------------ implementation of loops --------------------
3330
3331// internal helper functions that define how to call MPI data exchange
3332// functions: for generic vectors, do nothing at all. For distributed vectors,
3333// call update_ghost_values_start function and so on. If we have collections
3334// of vectors, just do the individual functions of the components. In order to
3335// keep ghost values consistent (whether we are in read or write mode), we
3336// also reset the values at the end. the whole situation is a bit complicated
3337// by the fact that we need to treat block vectors differently, which use some
3338// additional helper functions to select the blocks and template magic.
3339namespace internal
3340{
3344 template <int dim, typename Number, typename VectorizedArrayType>
3345 struct VectorDataExchange
3346 {
3347 // A shift for the MPI messages to reduce the risk for accidental
3348 // interaction with other open communications that a user program might
3349 // set up (parallel vectors support unfinished communication). We let
3350 // the other vectors use the first 20 assigned numbers and start the
3351 // matrix-free communication.
3352 static constexpr unsigned int channel_shift = 20;
3353
3354
3355
3360 VectorDataExchange(
3361 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3362 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3363 DataAccessOnFaces vector_face_access,
3364 const unsigned int n_components)
3365 : matrix_free(matrix_free)
3366 , vector_face_access(
3367 matrix_free.get_task_info().face_partition_data.empty() ?
3368 ::MatrixFree<dim, Number, VectorizedArrayType>::
3369 DataAccessOnFaces::unspecified :
3370 vector_face_access)
3371 , ghosts_were_set(false)
3372# ifdef DEAL_II_WITH_MPI
3373 , tmp_data(n_components)
3374 , requests(n_components)
3375# endif
3376 {
3377 (void)n_components;
3378 if (this->vector_face_access !=
3380 DataAccessOnFaces::unspecified)
3381 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3383 matrix_free.get_dof_info(c).vector_exchanger_face_variants.size(),
3384 5);
3385 }
3386
3387
3388
3392 ~VectorDataExchange() // NOLINT
3393 {
3394# ifdef DEAL_II_WITH_MPI
3395 for (unsigned int i = 0; i < tmp_data.size(); ++i)
3396 if (tmp_data[i] != nullptr)
3397 matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3398# endif
3399 }
3400
3401
3402
3407 template <typename VectorType>
3408 unsigned int
3409 find_vector_in_mf(const VectorType &vec,
3410 const bool check_global_compatibility = true) const
3411 {
3412 // case 1: if there is only one component we don't have anything to do
3413 if (matrix_free.n_components() == 1)
3414 {
3415 Assert(matrix_free.get_dof_info(0).vector_partitioner->is_compatible(
3416 *vec.get_partitioner()),
3417 ExcMessage("Could not find partitioner that fits vector"));
3418 return 0;
3419 }
3420
3421 // case 2: vector was set up with MatrixFree::initialize_dof_vector()
3422 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3423 if (vec.get_partitioner().get() ==
3424 matrix_free.get_dof_info(c).vector_partitioner.get())
3425 return c;
3426
3427 // case 3: user provided own partitioner (compatibility mode)
3428 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3429 if (check_global_compatibility ?
3430 vec.get_partitioner()->is_globally_compatible(
3431 *matrix_free.get_dof_info(c).vector_partitioner) :
3432 vec.get_partitioner()->is_compatible(
3433 *matrix_free.get_dof_info(c).vector_partitioner))
3434 return c;
3435
3436 Assert(false,
3437 ExcNotImplemented("Could not find partitioner that fits vector"));
3438
3440 }
3441
3442
3443
3449 get_partitioner(const unsigned int mf_component) const
3450 {
3451 AssertDimension(matrix_free.get_dof_info(mf_component)
3452 .vector_exchanger_face_variants.size(),
3453 5);
3454 if (vector_face_access ==
3456 DataAccessOnFaces::none)
3457 return *matrix_free.get_dof_info(mf_component)
3458 .vector_exchanger_face_variants[0];
3459 else if (vector_face_access ==
3461 DataAccessOnFaces::values)
3462 return *matrix_free.get_dof_info(mf_component)
3463 .vector_exchanger_face_variants[1];
3464 else if (vector_face_access ==
3466 DataAccessOnFaces::gradients)
3467 return *matrix_free.get_dof_info(mf_component)
3468 .vector_exchanger_face_variants[2];
3469 else if (vector_face_access ==
3471 DataAccessOnFaces::values_all_faces)
3472 return *matrix_free.get_dof_info(mf_component)
3473 .vector_exchanger_face_variants[3];
3474 else if (vector_face_access ==
3476 DataAccessOnFaces::gradients_all_faces)
3477 return *matrix_free.get_dof_info(mf_component)
3478 .vector_exchanger_face_variants[4];
3479 else
3480 return *matrix_free.get_dof_info(mf_component).vector_exchanger.get();
3481 }
3482
3483
3484
3488 template <typename VectorType,
3489 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3490 * = nullptr>
3491 void
3492 update_ghost_values_start(const unsigned int /*component_in_block_vector*/,
3493 const VectorType & /*vec*/)
3494 {}
3495
3496
3501 template <typename VectorType,
3502 std::enable_if_t<!has_update_ghost_values_start<VectorType> &&
3503 !is_not_parallel_vector<VectorType>,
3504 VectorType> * = nullptr>
3505 void
3506 update_ghost_values_start(const unsigned int component_in_block_vector,
3507 const VectorType &vec)
3508 {
3509 (void)component_in_block_vector;
3510 const bool ghosts_set = vec.has_ghost_elements();
3511
3512 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3513 ghosts_set == false,
3515
3516 if (ghosts_set)
3517 {
3518 ghosts_were_set = true;
3519 return;
3520 }
3521
3522 vec.update_ghost_values();
3523 }
3524
3525
3526
3532 template <typename VectorType,
3533 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3534 !has_exchange_on_subset<VectorType>,
3535 VectorType> * = nullptr>
3536 void
3537 update_ghost_values_start(const unsigned int component_in_block_vector,
3538 const VectorType &vec)
3539 {
3540 (void)component_in_block_vector;
3541 const bool ghosts_set = vec.has_ghost_elements();
3542
3543 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3544 ghosts_set == false,
3546
3547 if (ghosts_set)
3548 {
3549 ghosts_were_set = true;
3550 return;
3551 }
3552
3553 vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3554 }
3555
3556
3557
3564 template <typename VectorType,
3565 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3566 has_exchange_on_subset<VectorType>,
3567 VectorType> * = nullptr>
3568 void
3569 update_ghost_values_start(const unsigned int component_in_block_vector,
3570 const VectorType &vec)
3571 {
3572 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3573 "Type mismatch between VectorType and VectorDataExchange");
3574 (void)component_in_block_vector;
3575 const bool ghosts_set = vec.has_ghost_elements();
3576
3577 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3578 ghosts_set == false,
3580
3581 if (ghosts_set)
3582 {
3583 ghosts_were_set = true;
3584 return;
3585 }
3586
3587 if (vec.size() != 0)
3588 {
3589# ifdef DEAL_II_WITH_MPI
3590 const unsigned int mf_component = find_vector_in_mf(vec);
3591
3592 const auto &part = get_partitioner(mf_component);
3593
3594 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3595 part.n_import_sm_procs() == 0)
3596 return;
3597
3598 tmp_data[component_in_block_vector] =
3599 matrix_free.acquire_scratch_data_non_threadsafe();
3600 tmp_data[component_in_block_vector]->resize_fast(
3601 part.n_import_indices());
3602 AssertDimension(requests.size(), tmp_data.size());
3603
3604 part.export_to_ghosted_array_start(
3605 component_in_block_vector * 2 + channel_shift,
3606 ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3607 vec.shared_vector_data(),
3608 ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3609 part.locally_owned_size(),
3610 matrix_free.get_dof_info(mf_component)
3611 .vector_partitioner->n_ghost_indices()),
3612 ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3613 part.n_import_indices()),
3614 this->requests[component_in_block_vector]);
3615# endif
3616 }
3617 }
3618
3619
3620
3625 template <typename VectorType,
3626 std::enable_if_t<!has_update_ghost_values_start<VectorType>,
3627 VectorType> * = nullptr>
3628 void
3629 update_ghost_values_finish(const unsigned int /*component_in_block_vector*/,
3630 const VectorType & /*vec*/)
3631 {}
3632
3633
3634
3640 template <typename VectorType,
3641 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3642 !has_exchange_on_subset<VectorType>,
3643 VectorType> * = nullptr>
3644 void
3645 update_ghost_values_finish(const unsigned int component_in_block_vector,
3646 const VectorType &vec)
3647 {
3648 (void)component_in_block_vector;
3649
3650 if (ghosts_were_set)
3651 return;
3652
3653 vec.update_ghost_values_finish();
3654 }
3655
3656
3657
3664 template <typename VectorType,
3665 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3666 has_exchange_on_subset<VectorType>,
3667 VectorType> * = nullptr>
3668 void
3669 update_ghost_values_finish(const unsigned int component_in_block_vector,
3670 const VectorType &vec)
3671 {
3672 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3673 "Type mismatch between VectorType and VectorDataExchange");
3674 (void)component_in_block_vector;
3675
3676 if (ghosts_were_set)
3677 return;
3678
3679 if (vec.size() != 0)
3680 {
3681# ifdef DEAL_II_WITH_MPI
3682 AssertIndexRange(component_in_block_vector, tmp_data.size());
3683 AssertDimension(requests.size(), tmp_data.size());
3684
3685 const unsigned int mf_component = find_vector_in_mf(vec);
3686
3687 const auto &part = get_partitioner(mf_component);
3688
3689 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3690 part.n_import_sm_procs() != 0)
3691 {
3692 part.export_to_ghosted_array_finish(
3693 ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3694 vec.shared_vector_data(),
3695 ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3696 part.locally_owned_size(),
3697 matrix_free.get_dof_info(mf_component)
3698 .vector_partitioner->n_ghost_indices()),
3699 this->requests[component_in_block_vector]);
3700
3701 matrix_free.release_scratch_data_non_threadsafe(
3702 tmp_data[component_in_block_vector]);
3703 tmp_data[component_in_block_vector] = nullptr;
3704 }
3705# endif
3706 }
3707 // let vector know that ghosts are being updated and we can read from
3708 // them
3709 vec.set_ghost_state(true);
3710 }
3711
3712
3713
3717 template <typename VectorType,
3718 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3719 * = nullptr>
3720 void
3721 compress_start(const unsigned int /*component_in_block_vector*/,
3722 VectorType & /*vec*/)
3723 {}
3724
3725
3726
3731 template <typename VectorType,
3732 std::enable_if_t<!has_compress_start<VectorType> &&
3733 !is_not_parallel_vector<VectorType>,
3734 VectorType> * = nullptr>
3735 void
3736 compress_start(const unsigned int component_in_block_vector,
3737 VectorType &vec)
3738 {
3739 (void)component_in_block_vector;
3740 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3741 vec.compress(VectorOperation::add);
3742 }
3743
3744
3745
3751 template <typename VectorType,
3752 std::enable_if_t<has_compress_start<VectorType> &&
3753 !has_exchange_on_subset<VectorType>,
3754 VectorType> * = nullptr>
3755 void
3756 compress_start(const unsigned int component_in_block_vector,
3757 VectorType &vec)
3758 {
3759 (void)component_in_block_vector;
3760 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3761 vec.compress_start(component_in_block_vector + channel_shift);
3762 }
3763
3764
3765
3772 template <typename VectorType,
3773 std::enable_if_t<has_compress_start<VectorType> &&
3774 has_exchange_on_subset<VectorType>,
3775 VectorType> * = nullptr>
3776 void
3777 compress_start(const unsigned int component_in_block_vector,
3778 VectorType &vec)
3779 {
3780 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3781 "Type mismatch between VectorType and VectorDataExchange");
3782 (void)component_in_block_vector;
3783 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3784
3785 if (vec.size() != 0)
3786 {
3787# ifdef DEAL_II_WITH_MPI
3788 const unsigned int mf_component = find_vector_in_mf(vec);
3789
3790 const auto &part = get_partitioner(mf_component);
3791
3792 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3793 part.n_import_sm_procs() == 0)
3794 return;
3795
3796 tmp_data[component_in_block_vector] =
3797 matrix_free.acquire_scratch_data_non_threadsafe();
3798 tmp_data[component_in_block_vector]->resize_fast(
3799 part.n_import_indices());
3800 AssertDimension(requests.size(), tmp_data.size());
3801
3802 part.import_from_ghosted_array_start(
3804 component_in_block_vector * 2 + channel_shift,
3805 ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3806 vec.shared_vector_data(),
3807 ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3808 matrix_free.get_dof_info(mf_component)
3809 .vector_partitioner->n_ghost_indices()),
3810 ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3811 part.n_import_indices()),
3812 this->requests[component_in_block_vector]);
3813# endif
3814 }
3815 }
3816
3817
3818
3823 template <
3824 typename VectorType,
3825 std::enable_if_t<!has_compress_start<VectorType>, VectorType> * = nullptr>
3826 void
3827 compress_finish(const unsigned int /*component_in_block_vector*/,
3828 VectorType & /*vec*/)
3829 {}
3830
3831
3832
3838 template <typename VectorType,
3839 std::enable_if_t<has_compress_start<VectorType> &&
3840 !has_exchange_on_subset<VectorType>,
3841 VectorType> * = nullptr>
3842 void
3843 compress_finish(const unsigned int component_in_block_vector,
3844 VectorType &vec)
3845 {
3846 (void)component_in_block_vector;
3847 vec.compress_finish(VectorOperation::add);
3848 }
3849
3850
3851
3858 template <typename VectorType,
3859 std::enable_if_t<has_compress_start<VectorType> &&
3860 has_exchange_on_subset<VectorType>,
3861 VectorType> * = nullptr>
3862 void
3863 compress_finish(const unsigned int component_in_block_vector,
3864 VectorType &vec)
3865 {
3866 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3867 "Type mismatch between VectorType and VectorDataExchange");
3868 (void)component_in_block_vector;
3869 if (vec.size() != 0)
3870 {
3871# ifdef DEAL_II_WITH_MPI
3872 AssertIndexRange(component_in_block_vector, tmp_data.size());
3873 AssertDimension(requests.size(), tmp_data.size());
3874
3875 const unsigned int mf_component = find_vector_in_mf(vec);
3876
3877 const auto &part = get_partitioner(mf_component);
3878
3879 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3880 part.n_import_sm_procs() != 0)
3881 {
3882 part.import_from_ghosted_array_finish(
3884 ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3885 vec.shared_vector_data(),
3886 ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3887 matrix_free.get_dof_info(mf_component)
3888 .vector_partitioner->n_ghost_indices()),
3890 tmp_data[component_in_block_vector]->begin(),
3891 part.n_import_indices()),
3892 this->requests[component_in_block_vector]);
3893
3894 matrix_free.release_scratch_data_non_threadsafe(
3895 tmp_data[component_in_block_vector]);
3896 tmp_data[component_in_block_vector] = nullptr;
3897 }
3898
3900 {
3901 const int ierr =
3902 MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3903 AssertThrowMPI(ierr);
3904 }
3905# endif
3906 }
3907 }
3908
3909
3910
3914 template <typename VectorType,
3915 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3916 * = nullptr>
3917 void
3918 reset_ghost_values(const VectorType & /*vec*/) const
3919 {}
3920
3921
3922
3927 template <typename VectorType,
3928 std::enable_if_t<!has_exchange_on_subset<VectorType> &&
3929 !is_not_parallel_vector<VectorType>,
3930 VectorType> * = nullptr>
3931 void
3932 reset_ghost_values(const VectorType &vec) const
3933 {
3934 if (ghosts_were_set == true)
3935 return;
3936
3937 vec.zero_out_ghost_values();
3938 }
3939
3940
3941
3947 template <typename VectorType,
3948 std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3949 * = nullptr>
3950 void
3951 reset_ghost_values(const VectorType &vec) const
3952 {
3953 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3954 "Type mismatch between VectorType and VectorDataExchange");
3955 if (ghosts_were_set == true)
3956 return;
3957
3958 if (vec.size() != 0)
3959 {
3960# ifdef DEAL_II_WITH_MPI
3961 AssertDimension(requests.size(), tmp_data.size());
3962
3963 const unsigned int mf_component = find_vector_in_mf(vec);
3964
3965 const auto &part = get_partitioner(mf_component);
3966
3967 if (part.n_ghost_indices() > 0)
3968 {
3969 part.reset_ghost_values(
3970 ArrayView<Number>(const_cast<VectorType &>(vec).begin() +
3971 part.locally_owned_size(),
3972 matrix_free.get_dof_info(mf_component)
3973 .vector_partitioner->n_ghost_indices()));
3974 }
3975
3976# endif
3977 }
3978 // let vector know that it's not ghosted anymore
3979 vec.set_ghost_state(false);
3980 }
3981
3982
3983
3989 template <typename VectorType>
3990 void
3991 zero_vector_region(const unsigned int range_index, VectorType &vec) const
3992 {
3993 if constexpr (has_exchange_on_subset<VectorType>)
3994 {
3995 static_assert(
3996 std::is_same_v<Number, typename VectorType::value_type>,
3997 "Type mismatch between VectorType and VectorDataExchange");
3998 if (range_index == numbers::invalid_unsigned_int)
3999 vec = Number();
4000 else
4001 {
4002 const unsigned int mf_component = find_vector_in_mf(vec, false);
4004 matrix_free.get_dof_info(mf_component);
4005 Assert(dof_info.vector_zero_range_list_index.empty() == false,
4007
4008 Assert(vec.partitioners_are_compatible(
4009 *dof_info.vector_partitioner),
4011 AssertIndexRange(range_index,
4012 dof_info.vector_zero_range_list_index.size() -
4013 1);
4014 for (unsigned int id =
4015 dof_info.vector_zero_range_list_index[range_index];
4016 id != dof_info.vector_zero_range_list_index[range_index + 1];
4017 ++id)
4018 std::memset(vec.begin() +
4019 dof_info.vector_zero_range_list[id].first,
4020 0,
4021 (dof_info.vector_zero_range_list[id].second -
4022 dof_info.vector_zero_range_list[id].first) *
4023 sizeof(Number));
4024 }
4025 }
4026 else if constexpr (has_assignment_operator<VectorType>)
4027 {
4028 if (range_index == numbers::invalid_unsigned_int || range_index == 0)
4029 {
4030 if constexpr (std::is_same_v<
4032 VectorType>)
4033 {
4034 for (unsigned int i = 0; i < vec.size(); ++i)
4035 vec[i] = typename VectorType::value_type();
4036 }
4037 else
4038 vec = typename VectorType::value_type();
4039 }
4040 }
4041 else
4042 {
4043 Assert(false,
4045 "Zeroing is only implemented for vector types "
4046 "which provide operator=(const VectorType::value_type)"));
4047 }
4048 }
4049
4050
4051
4052 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
4053 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
4054 DataAccessOnFaces vector_face_access;
4055 bool ghosts_were_set;
4056# ifdef DEAL_II_WITH_MPI
4057 // Users typically do not attach many DoFHandlers to a MatrixFree object,
4058 // so make the default buffer size 16 to cover most cases
4059 boost::container::small_vector<AlignedVector<Number> *, 16> tmp_data;
4060 boost::container::small_vector<std::vector<MPI_Request>, 16> requests;
4061# endif
4062 }; // VectorDataExchange
4063
4064 template <typename VectorStruct>
4065 unsigned int
4066 n_components(const VectorStruct &vec);
4067
4068 template <typename VectorStruct>
4069 unsigned int
4070 n_components_block(const VectorStruct &vec, const std::bool_constant<true>)
4071 {
4072 unsigned int components = 0;
4073 for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
4074 components += n_components(vec.block(bl));
4075 return components;
4076 }
4077
4078 template <typename VectorStruct>
4079 unsigned int
4080 n_components_block(const VectorStruct &, const std::bool_constant<false>)
4081 {
4082 return 1;
4083 }
4084
4085 template <typename VectorStruct>
4086 unsigned int
4087 n_components(const VectorStruct &vec)
4088 {
4089 return n_components_block(
4090 vec, std::bool_constant<IsBlockVector<VectorStruct>::value>());
4091 }
4092
4093 template <typename VectorStruct>
4094 inline unsigned int
4095 n_components(const std::vector<VectorStruct> &vec)
4096 {
4097 unsigned int components = 0;
4098 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4099 components += n_components_block(
4100 vec[comp], std::bool_constant<IsBlockVector<VectorStruct>::value>());
4101 return components;
4102 }
4103
4104 template <typename VectorStruct>
4105 inline unsigned int
4106 n_components(const std::vector<VectorStruct *> &vec)
4107 {
4108 unsigned int components = 0;
4109 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4110 components += n_components_block(
4111 *vec[comp], std::bool_constant<IsBlockVector<VectorStruct>::value>());
4112 return components;
4113 }
4114
4115
4116
4117 // A helper function to identify block vectors with many components where we
4118 // should not try to overlap computations and communication because there
4119 // would be too many outstanding communication requests.
4120
4121 // default value for vectors that do not have communication_block_size
4122 template <typename VectorStruct,
4123 std::enable_if_t<!has_communication_block_size<VectorStruct>,
4124 VectorStruct> * = nullptr>
4125 constexpr unsigned int
4126 get_communication_block_size(const VectorStruct &)
4127 {
4129 }
4130
4131
4132
4133 template <typename VectorStruct,
4134 std::enable_if_t<has_communication_block_size<VectorStruct>,
4135 VectorStruct> * = nullptr>
4136 constexpr unsigned int
4137 get_communication_block_size(const VectorStruct &)
4138 {
4139 return VectorStruct::communication_block_size;
4140 }
4141
4142
4143
4144 template <typename VectorType,
4145 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType> * =
4146 nullptr>
4147 bool
4148 has_ghost_elements(const VectorType &vec)
4149 {
4150 (void)vec;
4151 return false;
4152 }
4153
4154
4155
4156 template <typename VectorType,
4157 std::enable_if_t<!is_not_parallel_vector<VectorType>, VectorType>
4158 * = nullptr>
4159 bool
4160 has_ghost_elements(const VectorType &vec)
4161 {
4162 return vec.has_ghost_elements();
4163 }
4164
4165
4166
4167 // --------------------------------------------------------------------------
4168 // below we have wrappers to distinguish between block and non-block vectors.
4169 // --------------------------------------------------------------------------
4170
4171 //
4172 // update_ghost_values_start
4173 //
4174
4175 // update_ghost_values for block vectors
4176 template <int dim,
4177 typename VectorStruct,
4178 typename Number,
4179 typename VectorizedArrayType,
4180 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4181 * = nullptr>
4182 void
4183 update_ghost_values_start(
4184 const VectorStruct &vec,
4185 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4186 const unsigned int channel = 0)
4187 {
4188 if (get_communication_block_size(vec) < vec.n_blocks())
4189 {
4190 const bool ghosts_set = vec.has_ghost_elements();
4191
4192 Assert(exchanger.matrix_free.get_task_info()
4193 .allow_ghosted_vectors_in_loops ||
4194 ghosts_set == false,
4196
4197 if (ghosts_set)
4198 {
4199 exchanger.ghosts_were_set = true;
4200 return;
4201 }
4202
4203 vec.update_ghost_values();
4204 }
4205 else
4206 {
4207 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4208 update_ghost_values_start(vec.block(i), exchanger, channel + i);
4209 }
4210 }
4211
4212
4213
4214 // update_ghost_values for non-block vectors
4215 template <int dim,
4216 typename VectorStruct,
4217 typename Number,
4218 typename VectorizedArrayType,
4219 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4220 * = nullptr>
4221 void
4222 update_ghost_values_start(
4223 const VectorStruct &vec,
4224 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4225 const unsigned int channel = 0)
4226 {
4227 exchanger.update_ghost_values_start(channel, vec);
4228 }
4229
4230
4231
4232 // update_ghost_values_start() for vector of vectors
4233 template <int dim,
4234 typename VectorStruct,
4235 typename Number,
4236 typename VectorizedArrayType>
4237 inline void
4238 update_ghost_values_start(
4239 const std::vector<VectorStruct> &vec,
4240 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4241 {
4242 unsigned int component_index = 0;
4243 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4244 {
4245 update_ghost_values_start(vec[comp], exchanger, component_index);
4246 component_index += n_components(vec[comp]);
4247 }
4248 }
4249
4250
4251
4252 // update_ghost_values_start() for vector of pointers to vectors
4253 template <int dim,
4254 typename VectorStruct,
4255 typename Number,
4256 typename VectorizedArrayType>
4257 inline void
4258 update_ghost_values_start(
4259 const std::vector<VectorStruct *> &vec,
4260 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4261 {
4262 unsigned int component_index = 0;
4263 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4264 {
4265 update_ghost_values_start(*vec[comp], exchanger, component_index);
4266 component_index += n_components(*vec[comp]);
4267 }
4268 }
4269
4270
4271
4272 //
4273 // update_ghost_values_finish
4274 //
4275
4276 // for block vectors
4277 template <int dim,
4278 typename VectorStruct,
4279 typename Number,
4280 typename VectorizedArrayType,
4281 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4282 * = nullptr>
4283 void
4284 update_ghost_values_finish(
4285 const VectorStruct &vec,
4286 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4287 const unsigned int channel = 0)
4288 {
4289 if (get_communication_block_size(vec) < vec.n_blocks())
4290 {
4291 // do nothing, everything has already been completed in the _start()
4292 // call
4293 }
4294 else
4295 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4296 update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4297 }
4298
4299
4300
4301 // for non-block vectors
4302 template <int dim,
4303 typename VectorStruct,
4304 typename Number,
4305 typename VectorizedArrayType,
4306 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4307 * = nullptr>
4308 void
4309 update_ghost_values_finish(
4310 const VectorStruct &vec,
4311 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4312 const unsigned int channel = 0)
4313 {
4314 exchanger.update_ghost_values_finish(channel, vec);
4315 }
4316
4317
4318
4319 // for vector of vectors
4320 template <int dim,
4321 typename VectorStruct,
4322 typename Number,
4323 typename VectorizedArrayType>
4324 inline void
4325 update_ghost_values_finish(
4326 const std::vector<VectorStruct> &vec,
4327 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4328 {
4329 unsigned int component_index = 0;
4330 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4331 {
4332 update_ghost_values_finish(vec[comp], exchanger, component_index);
4333 component_index += n_components(vec[comp]);
4334 }
4335 }
4336
4337
4338
4339 // for vector of pointers to vectors
4340 template <int dim,
4341 typename VectorStruct,
4342 typename Number,
4343 typename VectorizedArrayType>
4344 inline void
4345 update_ghost_values_finish(
4346 const std::vector<VectorStruct *> &vec,
4347 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4348 {
4349 unsigned int component_index = 0;
4350 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4351 {
4352 update_ghost_values_finish(*vec[comp], exchanger, component_index);
4353 component_index += n_components(*vec[comp]);
4354 }
4355 }
4356
4357
4358
4359 //
4360 // compress_start
4361 //
4362
4363 // for block vectors
4364 template <int dim,
4365 typename VectorStruct,
4366 typename Number,
4367 typename VectorizedArrayType,
4368 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4369 * = nullptr>
4370 inline void
4371 compress_start(
4372 VectorStruct &vec,
4373 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4374 const unsigned int channel = 0)
4375 {
4376 if (get_communication_block_size(vec) < vec.n_blocks())
4377 {
4378 // do nothing, everything will be done in the _finish() call
4379 }
4380 else
4381 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4382 compress_start(vec.block(i), exchanger, channel + i);
4383 }
4384
4385
4386
4387 // for non-block vectors
4388 template <int dim,
4389 typename VectorStruct,
4390 typename Number,
4391 typename VectorizedArrayType,
4392 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4393 * = nullptr>
4394 inline void
4395 compress_start(
4396 VectorStruct &vec,
4397 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4398 const unsigned int channel = 0)
4399 {
4400 exchanger.compress_start(channel, vec);
4401 }
4402
4403
4404
4405 // for std::vector of vectors
4406 template <int dim,
4407 typename VectorStruct,
4408 typename Number,
4409 typename VectorizedArrayType>
4410 inline void
4411 compress_start(
4412 std::vector<VectorStruct> &vec,
4413 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4414 {
4415 unsigned int component_index = 0;
4416 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4417 {
4418 compress_start(vec[comp], exchanger, component_index);
4419 component_index += n_components(vec[comp]);
4420 }
4421 }
4422
4423
4424
4425 // for std::vector of pointer to vectors
4426 template <int dim,
4427 typename VectorStruct,
4428 typename Number,
4429 typename VectorizedArrayType>
4430 inline void
4431 compress_start(
4432 std::vector<VectorStruct *> &vec,
4433 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4434 {
4435 unsigned int component_index = 0;
4436 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4437 {
4438 compress_start(*vec[comp], exchanger, component_index);
4439 component_index += n_components(*vec[comp]);
4440 }
4441 }
4442
4443
4444
4445 //
4446 // compress_finish
4447 //
4448
4449 // for block vectors
4450 template <int dim,
4451 typename VectorStruct,
4452 typename Number,
4453 typename VectorizedArrayType,
4454 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4455 * = nullptr>
4456 inline void
4457 compress_finish(
4458 VectorStruct &vec,
4459 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4460 const unsigned int channel = 0)
4461 {
4462 if (get_communication_block_size(vec) < vec.n_blocks())
4463 vec.compress(VectorOperation::add);
4464 else
4465 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4466 compress_finish(vec.block(i), exchanger, channel + i);
4467 }
4468
4469
4470
4471 // for non-block vectors
4472 template <int dim,
4473 typename VectorStruct,
4474 typename Number,
4475 typename VectorizedArrayType,
4476 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4477 * = nullptr>
4478 inline void
4479 compress_finish(
4480 VectorStruct &vec,
4481 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4482 const unsigned int channel = 0)
4483 {
4484 exchanger.compress_finish(channel, vec);
4485 }
4486
4487
4488
4489 // for std::vector of vectors
4490 template <int dim,
4491 typename VectorStruct,
4492 typename Number,
4493 typename VectorizedArrayType>
4494 inline void
4495 compress_finish(
4496 std::vector<VectorStruct> &vec,
4497 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4498 {
4499 unsigned int component_index = 0;
4500 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4501 {
4502 compress_finish(vec[comp], exchanger, component_index);
4503 component_index += n_components(vec[comp]);
4504 }
4505 }
4506
4507
4508
4509 // for std::vector of pointer to vectors
4510 template <int dim,
4511 typename VectorStruct,
4512 typename Number,
4513 typename VectorizedArrayType>
4514 inline void
4515 compress_finish(
4516 std::vector<VectorStruct *> &vec,
4517 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4518 {
4519 unsigned int component_index = 0;
4520 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4521 {
4522 compress_finish(*vec[comp], exchanger, component_index);
4523 component_index += n_components(*vec[comp]);
4524 }
4525 }
4526
4527
4528
4529 //
4530 // reset_ghost_values:
4531 //
4532 // if the input vector did not have ghosts imported, clear them here again
4533 // in order to avoid subsequent operations e.g. in linear solvers to work
4534 // with ghosts all the time
4535 //
4536
4537 // for block vectors
4538 template <int dim,
4539 typename VectorStruct,
4540 typename Number,
4541 typename VectorizedArrayType,
4542 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4543 * = nullptr>
4544 inline void
4545 reset_ghost_values(
4546 const VectorStruct &vec,
4547 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4548 {
4549 // return immediately if there is nothing to do.
4550 if (exchanger.ghosts_were_set == true)
4551 return;
4552
4553 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4554 reset_ghost_values(vec.block(i), exchanger);
4555 }
4556
4557
4558
4559 // for non-block vectors
4560 template <int dim,
4561 typename VectorStruct,
4562 typename Number,
4563 typename VectorizedArrayType,
4564 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4565 * = nullptr>
4566 inline void
4567 reset_ghost_values(
4568 const VectorStruct &vec,
4569 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4570 {
4571 // return immediately if there is nothing to do.
4572 if (exchanger.ghosts_were_set == true)
4573 return;
4574
4575 exchanger.reset_ghost_values(vec);
4576 }
4577
4578
4579
4580 // for std::vector of vectors
4581 template <int dim,
4582 typename VectorStruct,
4583 typename Number,
4584 typename VectorizedArrayType>
4585 inline void
4586 reset_ghost_values(
4587 const std::vector<VectorStruct> &vec,
4588 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4589 {
4590 // return immediately if there is nothing to do.
4591 if (exchanger.ghosts_were_set == true)
4592 return;
4593
4594 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4595 reset_ghost_values(vec[comp], exchanger);
4596 }
4597
4598
4599
4600 // for std::vector of pointer to vectors
4601 template <int dim,
4602 typename VectorStruct,
4603 typename Number,
4604 typename VectorizedArrayType>
4605 inline void
4606 reset_ghost_values(
4607 const std::vector<VectorStruct *> &vec,
4608 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4609 {
4610 // return immediately if there is nothing to do.
4611 if (exchanger.ghosts_were_set == true)
4612 return;
4613
4614 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4615 reset_ghost_values(*vec[comp], exchanger);
4616 }
4617
4618
4619
4620 //
4621 // zero_vector_region
4622 //
4623
4624 // for block vectors
4625 template <int dim,
4626 typename VectorStruct,
4627 typename Number,
4628 typename VectorizedArrayType,
4629 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4630 * = nullptr>
4631 inline void
4632 zero_vector_region(
4633 const unsigned int range_index,
4634 VectorStruct &vec,
4635 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4636 {
4637 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4638 exchanger.zero_vector_region(range_index, vec.block(i));
4639 }
4640
4641
4642
4643 // for non-block vectors
4644 template <int dim,
4645 typename VectorStruct,
4646 typename Number,
4647 typename VectorizedArrayType,
4648 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4649 * = nullptr>
4650 inline void
4651 zero_vector_region(
4652 const unsigned int range_index,
4653 VectorStruct &vec,
4654 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4655 {
4656 exchanger.zero_vector_region(range_index, vec);
4657 }
4658
4659
4660
4661 // for std::vector of vectors
4662 template <int dim,
4663 typename VectorStruct,
4664 typename Number,
4665 typename VectorizedArrayType>
4666 inline void
4667 zero_vector_region(
4668 const unsigned int range_index,
4669 std::vector<VectorStruct> &vec,
4670 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4671 {
4672 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4673 zero_vector_region(range_index, vec[comp], exchanger);
4674 }
4675
4676
4677
4678 // for std::vector of pointers to vectors
4679 template <int dim,
4680 typename VectorStruct,
4681 typename Number,
4682 typename VectorizedArrayType>
4683 inline void
4684 zero_vector_region(
4685 const unsigned int range_index,
4686 std::vector<VectorStruct *> &vec,
4687 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4688 {
4689 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4690 zero_vector_region(range_index, *vec[comp], exchanger);
4691 }
4692
4693
4694
4695 // Apply a unit matrix operation to constrained DoFs: Default cases where we
4696 // cannot detect a LinearAlgebra::distributed::Vector, we do not do
4697 // anything, else we apply the constraints as a unit operation
4698 template <typename VectorStruct1, typename VectorStruct2>
4699 inline void
4700 apply_operation_to_constrained_dofs(const std::vector<unsigned int> &,
4701 const VectorStruct1 &,
4702 VectorStruct2 &)
4703 {}
4704
4705 template <typename Number>
4706 inline void
4707 apply_operation_to_constrained_dofs(
4708 const std::vector<unsigned int> &constrained_dofs,
4711 {
4712 for (const unsigned int i : constrained_dofs)
4713 dst.local_element(i) = src.local_element(i);
4714 }
4715
4716
4717 namespace MatrixFreeFunctions
4718 {
4719 // struct to select between a const interface and a non-const interface
4720 // for MFWorker
4721 template <typename, typename, typename, typename, bool>
4722 struct InterfaceSelector
4723 {};
4724
4725 // Version for constant functions
4726 template <typename MF,
4727 typename InVector,
4728 typename OutVector,
4729 typename Container>
4730 struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4731 {
4732 using function_type = void (Container::*)(
4733 const MF &,
4734 OutVector &,
4735 const InVector &,
4736 const std::pair<unsigned int, unsigned int> &) const;
4737 };
4738
4739 // Version for non-constant functions
4740 template <typename MF,
4741 typename InVector,
4742 typename OutVector,
4743 typename Container>
4744 struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4745 {
4746 using function_type =
4747 void (Container::*)(const MF &,
4748 OutVector &,
4749 const InVector &,
4750 const std::pair<unsigned int, unsigned int> &);
4751 };
4752 } // namespace MatrixFreeFunctions
4753
4754
4755
4756 // A implementation class for the worker object that runs the various
4757 // operations we want to perform during the matrix-free loop
4758 template <typename MF,
4759 typename InVector,
4760 typename OutVector,
4761 typename Container,
4762 bool is_constant>
4763 class MFWorker : public MFWorkerInterface
4764 {
4765 public:
4766 // An alias to make the arguments further down more readable
4767 using function_type = typename MatrixFreeFunctions::
4768 InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4769 function_type;
4770
4771 // constructor, binds all the arguments to this class
4772 MFWorker(const MF &matrix_free,
4773 const InVector &src,
4774 OutVector &dst,
4775 const bool zero_dst_vector_setting,
4776 const Container &container,
4777 function_type cell_function,
4778 function_type face_function,
4779 function_type boundary_function,
4780 const typename MF::DataAccessOnFaces src_vector_face_access =
4781 MF::DataAccessOnFaces::none,
4782 const typename MF::DataAccessOnFaces dst_vector_face_access =
4783 MF::DataAccessOnFaces::none,
4784 const std::function<void(const unsigned int, const unsigned int)>
4785 &operation_before_loop = {},
4786 const std::function<void(const unsigned int, const unsigned int)>
4787 &operation_after_loop = {},
4788 const unsigned int dof_handler_index_pre_post = 0)
4789 : matrix_free(matrix_free)
4790 , container(const_cast<Container &>(container))
4791 , cell_function(cell_function)
4792 , face_function(face_function)
4793 , boundary_function(boundary_function)
4794 , src(src)
4795 , dst(dst)
4796 , src_data_exchanger(matrix_free,
4797 src_vector_face_access,
4798 n_components(src))
4799 , dst_data_exchanger(matrix_free,
4800 dst_vector_face_access,
4801 n_components(dst))
4802 , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4803 , zero_dst_vector_setting(zero_dst_vector_setting &&
4804 !src_and_dst_are_same)
4805 , operation_before_loop(operation_before_loop)
4806 , operation_after_loop(operation_after_loop)
4807 , dof_handler_index_pre_post(dof_handler_index_pre_post)
4808 {
4809 Assert(!has_ghost_elements(dst),
4810 ExcMessage("The destination vector passed to the matrix-free "
4811 "loop is ghosted. This is not allowed."));
4812 }
4813
4814 // Runs the cell work. If no function is given, nothing is done
4815 virtual void
4816 cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4817 {
4818 if (cell_function != nullptr && cell_range.second > cell_range.first)
4819 for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4820 {
4821 const auto cell_subrange =
4822 matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4823
4824 if (cell_subrange.second <= cell_subrange.first)
4825 continue;
4826
4827 (container.*
4828 cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4829 }
4830 }
4831
4832 virtual void
4833 cell(const unsigned int range_index) override
4834 {
4835 process_range(cell_function,
4836 matrix_free.get_task_info().cell_partition_data_hp_ptr,
4837 matrix_free.get_task_info().cell_partition_data_hp,
4838 range_index);
4839 }
4840
4841 virtual void
4842 face(const unsigned int range_index) override
4843 {
4844 process_range(face_function,
4845 matrix_free.get_task_info().face_partition_data_hp_ptr,
4846 matrix_free.get_task_info().face_partition_data_hp,
4847 range_index);
4848 }
4849
4850 virtual void
4851 boundary(const unsigned int range_index) override
4852 {
4853 process_range(boundary_function,
4854 matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4855 matrix_free.get_task_info().boundary_partition_data_hp,
4856 range_index);
4857 }
4858
4859 private:
4860 void
4861 process_range(const function_type &fu,
4862 const std::vector<unsigned int> &ptr,
4863 const std::vector<unsigned int> &data,
4864 const unsigned int range_index)
4865 {
4866 if (fu == nullptr)
4867 return;
4868
4869 AssertIndexRange(range_index + 1, ptr.size());
4870 for (unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4871 {
4872 AssertIndexRange(2 * i + 1, data.size());
4873 (container.*fu)(matrix_free,
4874 this->dst,
4875 this->src,
4876 std::make_pair(data[2 * i], data[2 * i + 1]));
4877 }
4878 }
4879
4880 public:
4881 // Starts the communication for the update ghost values operation. We
4882 // cannot call this update if ghost and destination are the same because
4883 // that would introduce spurious entries in the destination (there is also
4884 // the problem that reading from a vector that we also write to is usually
4885 // not intended in case there is overlap, but this is up to the
4886 // application code to decide and we cannot catch this case here).
4887 virtual void
4888 vector_update_ghosts_start() override
4889 {
4890 if (!src_and_dst_are_same)
4891 internal::update_ghost_values_start(src, src_data_exchanger);
4892 }
4893
4894 // Finishes the communication for the update ghost values operation
4895 virtual void
4896 vector_update_ghosts_finish() override
4897 {
4898 if (!src_and_dst_are_same)
4899 internal::update_ghost_values_finish(src, src_data_exchanger);
4900 }
4901
4902 // Starts the communication for the vector compress operation
4903 virtual void
4904 vector_compress_start() override
4905 {
4906 internal::compress_start(dst, dst_data_exchanger);
4907 }
4908
4909 // Finishes the communication for the vector compress operation
4910 virtual void
4911 vector_compress_finish() override
4912 {
4913 internal::compress_finish(dst, dst_data_exchanger);
4914 if (!src_and_dst_are_same)
4915 internal::reset_ghost_values(src, src_data_exchanger);
4916 }
4917
4918 // Zeros the given input vector
4919 virtual void
4920 zero_dst_vector_range(const unsigned int range_index) override
4921 {
4922 if (zero_dst_vector_setting)
4923 internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4924 }
4925
4926 virtual void
4927 cell_loop_pre_range(const unsigned int range_index) override
4928 {
4929 if (operation_before_loop)
4930 {
4932 matrix_free.get_dof_info(dof_handler_index_pre_post);
4933 if (range_index == numbers::invalid_unsigned_int)
4934 {
4935 // Case with threaded loop -> currently no overlap implemented
4937 0U,
4938 dof_info.vector_partitioner->locally_owned_size(),
4939 operation_before_loop,
4941 }
4942 else
4943 {
4944 AssertIndexRange(range_index,
4945 dof_info.cell_loop_pre_list_index.size() - 1);
4946 for (unsigned int id =
4947 dof_info.cell_loop_pre_list_index[range_index];
4948 id != dof_info.cell_loop_pre_list_index[range_index + 1];
4949 ++id)
4950 operation_before_loop(dof_info.cell_loop_pre_list[id].first,
4951 dof_info.cell_loop_pre_list[id].second);
4952 }
4953 }
4954 }
4955
4956 virtual void
4957 cell_loop_post_range(const unsigned int range_index) override
4958 {
4959 if (operation_after_loop)
4960 {
4961 // Run unit matrix operation on constrained dofs if we are at the
4962 // last range
4963 const std::vector<unsigned int> &partition_row_index =
4964 matrix_free.get_task_info().partition_row_index;
4965 if (range_index ==
4966 partition_row_index[partition_row_index.size() - 2] - 1)
4967 apply_operation_to_constrained_dofs(
4968 matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
4969 src,
4970 dst);
4971
4973 matrix_free.get_dof_info(dof_handler_index_pre_post);
4974 if (range_index == numbers::invalid_unsigned_int)
4975 {
4976 // Case with threaded loop -> currently no overlap implemented
4978 0U,
4979 dof_info.vector_partitioner->locally_owned_size(),
4980 operation_after_loop,
4982 }
4983 else
4984 {
4985 AssertIndexRange(range_index,
4986 dof_info.cell_loop_post_list_index.size() - 1);
4987 for (unsigned int id =
4988 dof_info.cell_loop_post_list_index[range_index];
4989 id != dof_info.cell_loop_post_list_index[range_index + 1];
4990 ++id)
4991 operation_after_loop(dof_info.cell_loop_post_list[id].first,
4992 dof_info.cell_loop_post_list[id].second);
4993 }
4994 }
4995 }
4996
4997 private:
4998 const MF &matrix_free;
4999 Container &container;
5000 function_type cell_function;
5001 function_type face_function;
5002 function_type boundary_function;
5003
5004 const InVector &src;
5005 OutVector &dst;
5006 VectorDataExchange<MF::dimension,
5007 typename MF::value_type,
5008 typename MF::vectorized_value_type>
5009 src_data_exchanger;
5010 VectorDataExchange<MF::dimension,
5011 typename MF::value_type,
5012 typename MF::vectorized_value_type>
5013 dst_data_exchanger;
5014 const bool src_and_dst_are_same;
5015 const bool zero_dst_vector_setting;
5016 const std::function<void(const unsigned int, const unsigned int)>
5017 operation_before_loop;
5018 const std::function<void(const unsigned int, const unsigned int)>
5019 operation_after_loop;
5020 const unsigned int dof_handler_index_pre_post;
5021 };
5022
5023
5024
5029 template <class MF, typename InVector, typename OutVector>
5030 struct MFClassWrapper
5031 {
5032 using function_type =
5033 std::function<void(const MF &,
5034 OutVector &,
5035 const InVector &,
5036 const std::pair<unsigned int, unsigned int> &)>;
5037
5038 MFClassWrapper(const function_type cell,
5039 const function_type face,
5040 const function_type boundary)
5041 : cell(cell)
5042 , face(face)
5043 , boundary(boundary)
5044 {}
5045
5046 void
5047 cell_integrator(const MF &mf,
5048 OutVector &dst,
5049 const InVector &src,
5050 const std::pair<unsigned int, unsigned int> &range) const
5051 {
5052 if (cell)
5053 cell(mf, dst, src, range);
5054 }
5055
5056 void
5057 face_integrator(const MF &mf,
5058 OutVector &dst,
5059 const InVector &src,
5060 const std::pair<unsigned int, unsigned int> &range) const
5061 {
5062 if (face)
5063 face(mf, dst, src, range);
5064 }
5065
5066 void
5067 boundary_integrator(
5068 const MF &mf,
5069 OutVector &dst,
5070 const InVector &src,
5071 const std::pair<unsigned int, unsigned int> &range) const
5072 {
5073 if (boundary)
5074 boundary(mf, dst, src, range);
5075 }
5076
5077 const function_type cell;
5078 const function_type face;
5079 const function_type boundary;
5080 };
5081
5082} // end of namespace internal
5083
5084
5085
5086template <int dim, typename Number, typename VectorizedArrayType>
5087template <typename OutVector, typename InVector>
5088inline void
5090 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5091 OutVector &,
5092 const InVector &,
5093 const std::pair<unsigned int, unsigned int> &)>
5094 &cell_operation,
5095 OutVector &dst,
5096 const InVector &src,
5097 const bool zero_dst_vector) const
5098{
5099 using Wrapper =
5100 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5101 InVector,
5102 OutVector>;
5103 Wrapper wrap(cell_operation, nullptr, nullptr);
5104 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5105 InVector,
5106 OutVector,
5107 Wrapper,
5108 true>
5109 worker(*this,
5110 src,
5111 dst,
5112 zero_dst_vector,
5113 wrap,
5114 &Wrapper::cell_integrator,
5115 &Wrapper::face_integrator,
5116 &Wrapper::boundary_integrator);
5117
5118 task_info.loop(worker);
5119}
5120
5121
5122
5123template <int dim, typename Number, typename VectorizedArrayType>
5124template <typename OutVector, typename InVector>
5125inline void
5127 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5128 OutVector &,
5129 const InVector &,
5130 const std::pair<unsigned int, unsigned int> &)>
5131 &cell_operation,
5132 OutVector &dst,
5133 const InVector &src,
5134 const std::function<void(const unsigned int, const unsigned int)>
5135 &operation_before_loop,
5136 const std::function<void(const unsigned int, const unsigned int)>
5137 &operation_after_loop,
5138 const unsigned int dof_handler_index_pre_post) const
5139{
5140 using Wrapper =
5141 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5142 InVector,
5143 OutVector>;
5144 Wrapper wrap(cell_operation, nullptr, nullptr);
5145 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5146 InVector,
5147 OutVector,
5148 Wrapper,
5149 true>
5150 worker(*this,
5151 src,
5152 dst,
5153 false,
5154 wrap,
5155 &Wrapper::cell_integrator,
5156 &Wrapper::face_integrator,
5157 &Wrapper::boundary_integrator,
5158 DataAccessOnFaces::none,
5159 DataAccessOnFaces::none,
5160 operation_before_loop,
5161 operation_after_loop,
5162 dof_handler_index_pre_post);
5163
5164 task_info.loop(worker);
5165}
5166
5167
5168
5169template <int dim, typename Number, typename VectorizedArrayType>
5170template <typename OutVector, typename InVector>
5171inline void
5173 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5174 OutVector &,
5175 const InVector &,
5176 const std::pair<unsigned int, unsigned int> &)>
5177 &cell_operation,
5178 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5179 OutVector &,
5180 const InVector &,
5181 const std::pair<unsigned int, unsigned int> &)>
5182 &inner_face_operation,
5183 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5184 OutVector &,
5185 const InVector &,
5186 const std::pair<unsigned int, unsigned int> &)>
5187 &boundary_face_operation,
5188 OutVector &dst,
5189 const InVector &src,
5190 const bool zero_dst_vector,
5191 const DataAccessOnFaces dst_vector_face_access,
5192 const DataAccessOnFaces src_vector_face_access) const
5193{
5194 using Wrapper =
5195 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5196 InVector,
5197 OutVector>;
5198 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5199 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5200 InVector,
5201 OutVector,
5202 Wrapper,
5203 true>
5204 worker(*this,
5205 src,
5206 dst,
5207 zero_dst_vector,
5208 wrap,
5209 &Wrapper::cell_integrator,
5210 &Wrapper::face_integrator,
5211 &Wrapper::boundary_integrator,
5212 src_vector_face_access,
5213 dst_vector_face_access);
5214
5215 task_info.loop(worker);
5216}
5217
5218
5219
5220template <int dim, typename Number, typename VectorizedArrayType>
5221template <typename CLASS, typename OutVector, typename InVector>
5222inline void
5224 void (CLASS::*function_pointer)(
5225 const MatrixFree<dim, Number, VectorizedArrayType> &,
5226 OutVector &,
5227 const InVector &,
5228 const std::pair<unsigned int, unsigned int> &) const,
5229 const CLASS *owning_class,
5230 OutVector &dst,
5231 const InVector &src,
5232 const bool zero_dst_vector) const
5233{
5234 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5235 InVector,
5236 OutVector,
5237 CLASS,
5238 true>
5239 worker(*this,
5240 src,
5241 dst,
5242 zero_dst_vector,
5243 *owning_class,
5244 function_pointer,
5245 nullptr,
5246 nullptr);
5247 task_info.loop(worker);
5248}
5249
5250
5251
5252template <int dim, typename Number, typename VectorizedArrayType>
5253template <typename CLASS, typename OutVector, typename InVector>
5254inline void
5256 void (CLASS::*function_pointer)(
5257 const MatrixFree<dim, Number, VectorizedArrayType> &,
5258 OutVector &,
5259 const InVector &,
5260 const std::pair<unsigned int, unsigned int> &) const,
5261 const CLASS *owning_class,
5262 OutVector &dst,
5263 const InVector &src,
5264 const std::function<void(const unsigned int, const unsigned int)>
5265 &operation_before_loop,
5266 const std::function<void(const unsigned int, const unsigned int)>
5267 &operation_after_loop,
5268 const unsigned int dof_handler_index_pre_post) const
5269{
5270 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5271 InVector,
5272 OutVector,
5273 CLASS,
5274 true>
5275 worker(*this,
5276 src,
5277 dst,
5278 false,
5279 *owning_class,
5280 function_pointer,
5281 nullptr,
5282 nullptr,
5285 operation_before_loop,
5286 operation_after_loop,
5287 dof_handler_index_pre_post);
5288 task_info.loop(worker);
5289}
5290
5291
5292
5293template <int dim, typename Number, typename VectorizedArrayType>
5294template <typename CLASS, typename OutVector, typename InVector>
5295inline void
5297 void (CLASS::*cell_operation)(
5298 const MatrixFree<dim, Number, VectorizedArrayType> &,
5299 OutVector &,
5300 const InVector &,
5301 const std::pair<unsigned int, unsigned int> &) const,
5302 void (CLASS::*inner_face_operation)(
5303 const MatrixFree<dim, Number, VectorizedArrayType> &,
5304 OutVector &,
5305 const InVector &,
5306 const std::pair<unsigned int, unsigned int> &) const,
5307 void (CLASS::*boundary_face_operation)(
5308 const MatrixFree<dim, Number, VectorizedArrayType> &,
5309 OutVector &,
5310 const InVector &,
5311 const std::pair<unsigned int, unsigned int> &) const,
5312 const CLASS *owning_class,
5313 OutVector &dst,
5314 const InVector &src,
5315 const bool zero_dst_vector,
5316 const DataAccessOnFaces dst_vector_face_access,
5317 const DataAccessOnFaces src_vector_face_access) const
5318{
5319 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5320 InVector,
5321 OutVector,
5322 CLASS,
5323 true>
5324 worker(*this,
5325 src,
5326 dst,
5327 zero_dst_vector,
5328 *owning_class,
5329 cell_operation,
5330 inner_face_operation,
5331 boundary_face_operation,
5332 src_vector_face_access,
5333 dst_vector_face_access);
5334 task_info.loop(worker);
5335}
5336
5337
5338
5339template <int dim, typename Number, typename VectorizedArrayType>
5340template <typename CLASS, typename OutVector, typename InVector>
5341inline void
5343 void (CLASS::*function_pointer)(
5344 const MatrixFree<dim, Number, VectorizedArrayType> &,
5345 OutVector &,
5346 const InVector &,
5347 const std::pair<unsigned int, unsigned int> &),
5348 CLASS *owning_class,
5349 OutVector &dst,
5350 const InVector &src,
5351 const bool zero_dst_vector) const
5352{
5353 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5354 InVector,
5355 OutVector,
5356 CLASS,
5357 false>
5358 worker(*this,
5359 src,
5360 dst,
5361 zero_dst_vector,
5362 *owning_class,
5363 function_pointer,
5364 nullptr,
5365 nullptr);
5366 task_info.loop(worker);
5367}
5368
5369
5370
5371template <int dim, typename Number, typename VectorizedArrayType>
5372template <typename CLASS, typename OutVector, typename InVector>
5373inline void
5375 void (CLASS::*function_pointer)(
5376 const MatrixFree<dim, Number, VectorizedArrayType> &,
5377 OutVector &,
5378 const InVector &,
5379 const std::pair<unsigned int, unsigned int> &),
5380 CLASS *owning_class,
5381 OutVector &dst,
5382 const InVector &src,
5383 const std::function<void(const unsigned int, const unsigned int)>
5384 &operation_before_loop,
5385 const std::function<void(const unsigned int, const unsigned int)>
5386 &operation_after_loop,
5387 const unsigned int dof_handler_index_pre_post) const
5388{
5389 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5390 InVector,
5391 OutVector,
5392 CLASS,
5393 false>
5394 worker(*this,
5395 src,
5396 dst,
5397 false,
5398 *owning_class,
5399 function_pointer,
5400 nullptr,
5401 nullptr,
5404 operation_before_loop,
5405 operation_after_loop,
5406 dof_handler_index_pre_post);
5407 task_info.loop(worker);
5408}
5409
5410
5411
5412template <int dim, typename Number, typename VectorizedArrayType>
5413template <typename CLASS, typename OutVector, typename InVector>
5414inline void
5416 void (CLASS::*cell_operation)(
5417 const MatrixFree<dim, Number, VectorizedArrayType> &,
5418 OutVector &,
5419 const InVector &,
5420 const std::pair<unsigned int, unsigned int> &),
5421 void (CLASS::*inner_face_operation)(
5422 const MatrixFree<dim, Number, VectorizedArrayType> &,
5423 OutVector &,
5424 const InVector &,
5425 const std::pair<unsigned int, unsigned int> &),
5426 void (CLASS::*boundary_face_operation)(
5427 const MatrixFree<dim, Number, VectorizedArrayType> &,
5428 OutVector &,
5429 const InVector &,
5430 const std::pair<unsigned int, unsigned int> &),
5431 CLASS *owning_class,
5432 OutVector &dst,
5433 const InVector &src,
5434 const bool zero_dst_vector,
5435 const DataAccessOnFaces dst_vector_face_access,
5436 const DataAccessOnFaces src_vector_face_access) const
5437{
5438 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5439 InVector,
5440 OutVector,
5441 CLASS,
5442 false>
5443 worker(*this,
5444 src,
5445 dst,
5446 zero_dst_vector,
5447 *owning_class,
5448 cell_operation,
5449 inner_face_operation,
5450 boundary_face_operation,
5451 src_vector_face_access,
5452 dst_vector_face_access);
5453 task_info.loop(worker);
5454}
5455
5456
5457
5458template <int dim, typename Number, typename VectorizedArrayType>
5459template <typename OutVector, typename InVector>
5460inline void
5462 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5463 OutVector &,
5464 const InVector &,
5465 const std::pair<unsigned int, unsigned int> &)>
5466 &cell_operation,
5467 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5468 OutVector &,
5469 const InVector &,
5470 const std::pair<unsigned int, unsigned int> &)>
5471 &inner_face_operation,
5472 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5473 OutVector &,
5474 const InVector &,
5475 const std::pair<unsigned int, unsigned int> &)>
5476 &boundary_face_operation,
5477 OutVector &dst,
5478 const InVector &src,
5479 const std::function<void(const unsigned int, const unsigned int)>
5480 &operation_before_loop,
5481 const std::function<void(const unsigned int, const unsigned int)>
5482 &operation_after_loop,
5483 const unsigned int dof_handler_index_pre_post,
5484 const DataAccessOnFaces dst_vector_face_access,
5485 const DataAccessOnFaces src_vector_face_access) const
5486{
5487 using Wrapper =
5488 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5489 InVector,
5490 OutVector>;
5491 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5492 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5493 InVector,
5494 OutVector,
5495 Wrapper,
5496 true>
5497 worker(*this,
5498 src,
5499 dst,
5500 false,
5501 wrap,
5502 &Wrapper::cell_integrator,
5503 &Wrapper::face_integrator,
5504 &Wrapper::boundary_integrator,
5505 src_vector_face_access,
5506 dst_vector_face_access,
5507 operation_before_loop,
5508 operation_after_loop,
5509 dof_handler_index_pre_post);
5510
5511 task_info.loop(worker);
5512}
5513
5514
5515
5516template <int dim, typename Number, typename VectorizedArrayType>
5517template <typename CLASS, typename OutVector, typename InVector>
5518inline void
5520 void (CLASS::*cell_operation)(const MatrixFree &,
5521 OutVector &,
5522 const InVector &,
5523 const std::pair<unsigned int, unsigned int> &)
5524 const,
5525 void (CLASS::*inner_face_operation)(
5526 const MatrixFree &,
5527 OutVector &,
5528 const InVector &,
5529 const std::pair<unsigned int, unsigned int> &) const,
5530 void (CLASS::*boundary_face_operation)(
5531 const MatrixFree &,
5532 OutVector &,
5533 const InVector &,
5534 const std::pair<unsigned int, unsigned int> &) const,
5535 const CLASS *owning_class,
5536 OutVector &dst,
5537 const InVector &src,
5538 const std::function<void(const unsigned int, const unsigned int)>
5539 &operation_before_loop,
5540 const std::function<void(const unsigned int, const unsigned int)>
5541 &operation_after_loop,
5542 const unsigned int dof_handler_index_pre_post,
5543 const DataAccessOnFaces dst_vector_face_access,
5544 const DataAccessOnFaces src_vector_face_access) const
5545{
5546 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5547 InVector,
5548 OutVector,
5549 CLASS,
5550 true>
5551 worker(*this,
5552 src,
5553 dst,
5554 false,
5555 *owning_class,
5556 cell_operation,
5557 inner_face_operation,
5558 boundary_face_operation,
5559 src_vector_face_access,
5560 dst_vector_face_access,
5561 operation_before_loop,
5562 operation_after_loop,
5563 dof_handler_index_pre_post);
5564 task_info.loop(worker);
5565}
5566
5567
5568
5569template <int dim, typename Number, typename VectorizedArrayType>
5570template <typename CLASS, typename OutVector, typename InVector>
5571inline void
5573 void (CLASS::*cell_operation)(const MatrixFree &,
5574 OutVector &,
5575 const InVector &,
5576 const std::pair<unsigned int, unsigned int> &),
5577 void (CLASS::*inner_face_operation)(
5578 const MatrixFree &,
5579 OutVector &,
5580 const InVector &,
5581 const std::pair<unsigned int, unsigned int> &),
5582 void (CLASS::*boundary_face_operation)(
5583 const MatrixFree &,
5584 OutVector &,
5585 const InVector &,
5586 const std::pair<unsigned int, unsigned int> &),
5587 const CLASS *owning_class,
5588 OutVector &dst,
5589 const InVector &src,
5590 const std::function<void(const unsigned int, const unsigned int)>
5591 &operation_before_loop,
5592 const std::function<void(const unsigned int, const unsigned int)>
5593 &operation_after_loop,
5594 const unsigned int dof_handler_index_pre_post,
5595 const DataAccessOnFaces dst_vector_face_access,
5596 const DataAccessOnFaces src_vector_face_access) const
5597{
5598 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5599 InVector,
5600 OutVector,
5601 CLASS,
5602 false>
5603 worker(*this,
5604 src,
5605 dst,
5606 false,
5607 *owning_class,
5608 cell_operation,
5609 inner_face_operation,
5610 boundary_face_operation,
5611 src_vector_face_access,
5612 dst_vector_face_access,
5613 operation_before_loop,
5614 operation_after_loop,
5615 dof_handler_index_pre_post);
5616 task_info.loop(worker);
5617}
5618
5619
5620
5621template <int dim, typename Number, typename VectorizedArrayType>
5622template <typename CLASS, typename OutVector, typename InVector>
5623inline void
5625 void (CLASS::*function_pointer)(
5626 const MatrixFree<dim, Number, VectorizedArrayType> &,
5627 OutVector &,
5628 const InVector &,
5629 const std::pair<unsigned int, unsigned int> &) const,
5630 const CLASS *owning_class,
5631 OutVector &dst,
5632 const InVector &src,
5633 const bool zero_dst_vector,
5634 const DataAccessOnFaces src_vector_face_access) const
5635{
5636 auto src_vector_face_access_temp = src_vector_face_access;
5637 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5638 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5639 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5640 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5641
5642 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5643 InVector,
5644 OutVector,
5645 CLASS,
5646 true>
5647 worker(*this,
5648 src,
5649 dst,
5650 zero_dst_vector,
5651 *owning_class,
5652 function_pointer,
5653 nullptr,
5654 nullptr,
5655 src_vector_face_access_temp,
5657 task_info.loop(worker);
5658}
5659
5660
5661
5662template <int dim, typename Number, typename VectorizedArrayType>
5663template <typename CLASS, typename OutVector, typename InVector>
5664inline void
5666 void (CLASS::*function_pointer)(
5667 const MatrixFree<dim, Number, VectorizedArrayType> &,
5668 OutVector &,
5669 const InVector &,
5670 const std::pair<unsigned int, unsigned int> &),
5671 CLASS *owning_class,
5672 OutVector &dst,
5673 const InVector &src,
5674 const bool zero_dst_vector,
5675 const DataAccessOnFaces src_vector_face_access) const
5676{
5677 auto src_vector_face_access_temp = src_vector_face_access;
5678 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5679 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5680 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5681 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5682
5683 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5684 InVector,
5685 OutVector,
5686 CLASS,
5687 false>
5688 worker(*this,
5689 src,
5690 dst,
5691 zero_dst_vector,
5692 *owning_class,
5693 function_pointer,
5694 nullptr,
5695 nullptr,
5696 src_vector_face_access_temp,
5698 task_info.loop(worker);
5699}
5700
5701
5702
5703template <int dim, typename Number, typename VectorizedArrayType>
5704template <typename OutVector, typename InVector>
5705inline void
5707 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5708 OutVector &,
5709 const InVector &,
5710 const std::pair<unsigned int, unsigned int> &)>
5711 &cell_operation,
5712 OutVector &dst,
5713 const InVector &src,
5714 const bool zero_dst_vector,
5715 const DataAccessOnFaces src_vector_face_access) const
5716{
5717 auto src_vector_face_access_temp = src_vector_face_access;
5718 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5719 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5720 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5721 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5722
5723 using Wrapper =
5724 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5725 InVector,
5726 OutVector>;
5727 Wrapper wrap(cell_operation, nullptr, nullptr);
5728
5729 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5730 InVector,
5731 OutVector,
5732 Wrapper,
5733 true>
5734 worker(*this,
5735 src,
5736 dst,
5737 zero_dst_vector,
5738 wrap,
5739 &Wrapper::cell_integrator,
5740 &Wrapper::face_integrator,
5741 &Wrapper::boundary_integrator,
5742 src_vector_face_access_temp,
5743 DataAccessOnFaces::none);
5744 task_info.loop(worker);
5745}
5746
5747
5748#endif // ifndef DOXYGEN
5749
5750
5751
5753
5754#endif
*  *  for(const auto &cell :triangulation.active_cell_iterators())
*  *  iterator begin()
iterator begin()
void resize(const size_type new_size)
void collect_sizes()
BlockType & block(const unsigned int i)
void reinit(const size_type num_blocks, const size_type block_size=0, const bool omit_zeroing_entries=false)
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
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
unsigned int n_ghost_cell_batches() const
void initialize_dof_handlers(const std::vector< const DoFHandler< dim, dim > * > &dof_handlers, const AdditionalData &additional_data)
types::boundary_id get_boundary_id(const unsigned int face_batch_index) const
void print(std::ostream &out) const
void update_mapping(const std::shared_ptr< hp::MappingCollection< dim > > &mapping)
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
unsigned int n_inner_face_batches() const
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
void update_mapping(const Mapping< dim > &mapping)
const Quadrature< dim > & get_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_index) const
unsigned int get_mg_level() const
unsigned int get_cell_category(const unsigned int cell_batch_index, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*boundary_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void print_memory_consumption(StreamType &out) const
std::pair< unsigned int, unsigned int > get_face_range_category(const std::pair< unsigned int, unsigned int > face_batch_range, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) const
bool mapping_initialized() const
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
void clear()
void internal_reinit(const std::shared_ptr< hp::MappingCollection< dim > > &mapping, const std::vector< const DoFHandler< dim, dim > * > &dof_handlers, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const std::vector< hp::QCollection< q_dim > > &quad, const AdditionalData &additional_data)
unsigned int get_n_q_points_face(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &inner_face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_face_operation, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
bool at_irregular_cell(const unsigned int cell_batch_index) const
~MatrixFree() override=default
const AffineConstraints< Number > & get_affine_constraints(const unsigned int dof_handler_index=0) const
void copy_from(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free_base)
void initialize_face_data_vector(AlignedVector< T > &vec) const
unsigned int get_dofs_per_cell(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int n_constraint_pool_entries() const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
const Number * constraint_pool_begin(const unsigned int pool_index) const
AlignedVector< Number > * acquire_scratch_data_non_threadsafe() const
const IndexSet & get_locally_owned_set(const unsigned int dof_handler_index=0) const
unsigned int get_dofs_per_face(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
void initialize_dof_vector(LinearAlgebra::distributed::Vector< Number2, MemorySpace > &vec, const unsigned int dof_handler_index=0) const
internal::MatrixFreeFunctions::TaskInfo task_info
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
std::vector< ObserverPointer< const DoFHandler< dim > > > dof_handlers
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*boundary_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
unsigned int first_hp_dof_handler_index
bool mapping_is_initialized
void loop_cell_centric(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
DoFHandler< dim >::cell_iterator get_cell_iterator(const unsigned int cell_batch_index, const unsigned int lane_index, const unsigned int dof_handler_index=0) const
MatrixFree(const MatrixFree< dim, Number, VectorizedArrayType > &other)
std::pair< unsigned int, unsigned int > create_cell_subrange_hp(const std::pair< unsigned int, unsigned int > &range, const unsigned int fe_degree, const unsigned int dof_handler_index=0) const
std::pair< unsigned int, unsigned int > create_cell_subrange_hp_by_index(const std::pair< unsigned int, unsigned int > &range, const unsigned int fe_index, const unsigned int dof_handler_index=0) const
unsigned int mg_level
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
std::vector< std::pair< unsigned int, unsigned int > > cell_level_index
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner(const unsigned int dof_handler_index=0) const
std::vector< ObserverPointer< const AffineConstraints< Number > > > affine_constraints
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
void initialize_dof_vector(LinearAlgebra::distributed::BlockVector< Number2, MemorySpace > &vec) const
std::vector< Number > constraint_pool_data
VectorizedArrayType vectorized_value_type
unsigned int n_cell_batches() const
const IndexSet & get_ghost_set(const unsigned int dof_handler_index=0) const
std::pair< int, int > get_cell_level_and_index(const unsigned int cell_batch_index, const unsigned int lane_index) const
bool indices_initialized() const
const Quadrature< dim - 1 > & get_face_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int get_face_active_fe_index(const std::pair< unsigned int, unsigned int > range, const bool is_interior_face=true, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) const
unsigned int get_cell_active_fe_index(const std::pair< unsigned int, unsigned int > range, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) const
Threads::ThreadLocalStorage< std::list< std::pair< bool, AlignedVector< VectorizedArrayType > > > > scratch_pad
unsigned int get_matrix_free_cell_index(const typename Triangulation< dim >::cell_iterator &cell) const
Number value_type
std::size_t memory_consumption() const
std::pair< unsigned int, unsigned int > get_face_category(const unsigned int face_batch_index, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) const
unsigned int n_boundary_face_batches() const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*boundary_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
std::list< std::pair< bool, AlignedVector< Number > > > scratch_pad_non_threadsafe
const Number * constraint_pool_end(const unsigned int pool_index) const
void loop_cell_centric(const std::function< void(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &inner_face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_face_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void initialize_indices(const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const AdditionalData &additional_data)
unsigned int n_components() const
unsigned int n_physical_cells() const
void initialize_dof_vector(VectorType &vec, const unsigned int dof_handler_index=0) const
unsigned int get_n_q_points(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
void reinit(const MappingType &mapping, const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
std::vector< unsigned int > mf_cell_indices
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
unsigned int get_cell_range_category(const std::pair< unsigned int, unsigned int > cell_batch_range, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*boundary_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
std::vector< internal::MatrixFreeFunctions::DoFInfo > dof_info
void initialize_cell_data_vector(AlignedVector< T > &vec) const
unsigned int n_ghost_inner_face_batches() const
void release_scratch_data_non_threadsafe(const AlignedVector< Number > *memory) const
void renumber_dofs(std::vector< types::global_dof_index > &renumbering, const unsigned int dof_handler_index=0)
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
bool indices_are_initialized
internal::MatrixFreeFunctions::FaceInfo< VectorizedArrayType::size()> face_info
void reinit(const MappingType &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< number2 > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
void reinit(const MappingType &mapping, const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< QuadratureType > &quad, const AdditionalData &additional_data=AdditionalData())
std::array< types::boundary_id, VectorizedArrayType::size()> get_faces_by_cells_boundary_id(const unsigned int cell_batch_index, const unsigned int face_number) const
std::vector< unsigned int > constraint_pool_row_index
const internal::MatrixFreeFunctions::ShapeInfo< Number > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
const std::vector< unsigned int > & get_constrained_dofs(const unsigned int dof_handler_index=0) const
unsigned int cell_level_index_end_local
unsigned int n_base_elements(const unsigned int dof_handler_index) const
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< Number > > shape_info
static bool is_supported(const FiniteElement< dim, spacedim > &fe)
static constexpr unsigned int dimension
A class that provides a separate storage location on each thread that accesses the object.
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
Point< 2 > second
Definition grid_out.cc:4640
Point< 2 > first
Definition grid_out.cc:4639
unsigned int level
Definition grid_out.cc:4642
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ActiveSelector::cell_iterator cell_iterator
UpdateFlags
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_default
No update.
std::vector< index_type > data
Definition mpi.cc:734
std::size_t size
Definition mpi.cc:733
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
*  *  *  RotationFunction< dim, Number >::RotationFunction Number(dim)
*  *  *  *  std::vector< Number > ThermoPlasticMaterial< dim, ViscoplasticYieldLaw, Number >::get_state_parameters   const
bool job_supports_mpi()
Definition mpi.cc:680
unsigned int minimum_parallel_grain_size
Definition parallel.cc:48
constexpr unsigned int invalid_unsigned_int
Definition types.h:228
constexpr types::boundary_id invalid_boundary_id
Definition types.h:299
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:266
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
unsigned short int fe_index
Definition types.h:70
unsigned int boundary_id
Definition types.h:159
TasksParallelScheme tasks_parallel_scheme
UpdateFlags mapping_update_flags_inner_faces
AdditionalData(const TasksParallelScheme tasks_parallel_scheme=partition_partition, const unsigned int tasks_block_size=0, const UpdateFlags mapping_update_flags=update_gradients|update_JxW_values, const UpdateFlags mapping_update_flags_boundary_faces=update_default, const UpdateFlags mapping_update_flags_inner_faces=update_default, const UpdateFlags mapping_update_flags_faces_by_cells=update_default, const unsigned int mg_level=numbers::invalid_unsigned_int, const bool store_plain_indices=true, const bool initialize_indices=true, const bool initialize_mapping=true, const bool overlap_communication_computation=true, const bool hold_all_faces_to_owned_cells=false, const bool cell_vectorization_categories_strict=false, const bool allow_ghosted_vectors_in_loops=true)
std::vector< unsigned int > cell_vectorization_category
UpdateFlags mapping_update_flags_boundary_faces
UpdateFlags mapping_update_flags_faces_by_cells
AdditionalData & operator=(const AdditionalData &other)=default
AdditionalData(const AdditionalData &other)
std::vector< unsigned int > cell_loop_pre_list_index
Definition dof_info.h:747
std::vector< unsigned int > cell_loop_post_list_index
Definition dof_info.h:760
std::vector< std::pair< unsigned int, unsigned int > > vector_zero_range_list
Definition dof_info.h:740
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition dof_info.h:591
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_pre_list
Definition dof_info.h:753
std::vector< unsigned int > vector_zero_range_list_index
Definition dof_info.h:735
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_post_list
Definition dof_info.h:766
void loop(MFWorkerInterface &worker) const
Definition task_info.cc:348