deal.II version GIT relicensing-6809-ge913b9bb34 2026-09-25 17:20:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
matrix_free.h
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: 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 vec.compress(VectorOperation::add);
4378 else
4379 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4380 compress_start(vec.block(i), exchanger, channel + i);
4381 }
4382
4383
4384
4385 // for non-block vectors
4386 template <int dim,
4387 typename VectorStruct,
4388 typename Number,
4389 typename VectorizedArrayType,
4390 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4391 * = nullptr>
4392 inline void
4393 compress_start(
4394 VectorStruct &vec,
4395 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4396 const unsigned int channel = 0)
4397 {
4398 exchanger.compress_start(channel, vec);
4399 }
4400
4401
4402
4403 // for std::vector of vectors
4404 template <int dim,
4405 typename VectorStruct,
4406 typename Number,
4407 typename VectorizedArrayType>
4408 inline void
4409 compress_start(
4410 std::vector<VectorStruct> &vec,
4411 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4412 {
4413 unsigned int component_index = 0;
4414 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4415 {
4416 compress_start(vec[comp], exchanger, component_index);
4417 component_index += n_components(vec[comp]);
4418 }
4419 }
4420
4421
4422
4423 // for std::vector of pointer to vectors
4424 template <int dim,
4425 typename VectorStruct,
4426 typename Number,
4427 typename VectorizedArrayType>
4428 inline void
4429 compress_start(
4430 std::vector<VectorStruct *> &vec,
4431 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4432 {
4433 unsigned int component_index = 0;
4434 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4435 {
4436 compress_start(*vec[comp], exchanger, component_index);
4437 component_index += n_components(*vec[comp]);
4438 }
4439 }
4440
4441
4442
4443 //
4444 // compress_finish
4445 //
4446
4447 // for block vectors
4448 template <int dim,
4449 typename VectorStruct,
4450 typename Number,
4451 typename VectorizedArrayType,
4452 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4453 * = nullptr>
4454 inline void
4455 compress_finish(
4456 VectorStruct &vec,
4457 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4458 const unsigned int channel = 0)
4459 {
4460 if (get_communication_block_size(vec) < vec.n_blocks())
4461 {
4462 // do nothing, everything has already been completed in the _start()
4463 // call
4464 }
4465 else
4466 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4467 compress_finish(vec.block(i), exchanger, channel + i);
4468 }
4469
4470
4471
4472 // for non-block vectors
4473 template <int dim,
4474 typename VectorStruct,
4475 typename Number,
4476 typename VectorizedArrayType,
4477 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4478 * = nullptr>
4479 inline void
4480 compress_finish(
4481 VectorStruct &vec,
4482 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4483 const unsigned int channel = 0)
4484 {
4485 exchanger.compress_finish(channel, vec);
4486 }
4487
4488
4489
4490 // for std::vector of vectors
4491 template <int dim,
4492 typename VectorStruct,
4493 typename Number,
4494 typename VectorizedArrayType>
4495 inline void
4496 compress_finish(
4497 std::vector<VectorStruct> &vec,
4498 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4499 {
4500 unsigned int component_index = 0;
4501 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4502 {
4503 compress_finish(vec[comp], exchanger, component_index);
4504 component_index += n_components(vec[comp]);
4505 }
4506 }
4507
4508
4509
4510 // for std::vector of pointer to vectors
4511 template <int dim,
4512 typename VectorStruct,
4513 typename Number,
4514 typename VectorizedArrayType>
4515 inline void
4516 compress_finish(
4517 std::vector<VectorStruct *> &vec,
4518 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4519 {
4520 unsigned int component_index = 0;
4521 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4522 {
4523 compress_finish(*vec[comp], exchanger, component_index);
4524 component_index += n_components(*vec[comp]);
4525 }
4526 }
4527
4528
4529
4530 //
4531 // reset_ghost_values:
4532 //
4533 // if the input vector did not have ghosts imported, clear them here again
4534 // in order to avoid subsequent operations e.g. in linear solvers to work
4535 // with ghosts all the time
4536 //
4537
4538 // for block vectors
4539 template <int dim,
4540 typename VectorStruct,
4541 typename Number,
4542 typename VectorizedArrayType,
4543 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4544 * = nullptr>
4545 inline void
4546 reset_ghost_values(
4547 const VectorStruct &vec,
4548 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4549 {
4550 // return immediately if there is nothing to do.
4551 if (exchanger.ghosts_were_set == true)
4552 return;
4553
4554 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4555 reset_ghost_values(vec.block(i), exchanger);
4556 }
4557
4558
4559
4560 // for non-block vectors
4561 template <int dim,
4562 typename VectorStruct,
4563 typename Number,
4564 typename VectorizedArrayType,
4565 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4566 * = nullptr>
4567 inline void
4568 reset_ghost_values(
4569 const VectorStruct &vec,
4570 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4571 {
4572 // return immediately if there is nothing to do.
4573 if (exchanger.ghosts_were_set == true)
4574 return;
4575
4576 exchanger.reset_ghost_values(vec);
4577 }
4578
4579
4580
4581 // for std::vector of vectors
4582 template <int dim,
4583 typename VectorStruct,
4584 typename Number,
4585 typename VectorizedArrayType>
4586 inline void
4587 reset_ghost_values(
4588 const std::vector<VectorStruct> &vec,
4589 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4590 {
4591 // return immediately if there is nothing to do.
4592 if (exchanger.ghosts_were_set == true)
4593 return;
4594
4595 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4596 reset_ghost_values(vec[comp], exchanger);
4597 }
4598
4599
4600
4601 // for std::vector of pointer to vectors
4602 template <int dim,
4603 typename VectorStruct,
4604 typename Number,
4605 typename VectorizedArrayType>
4606 inline void
4607 reset_ghost_values(
4608 const std::vector<VectorStruct *> &vec,
4609 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4610 {
4611 // return immediately if there is nothing to do.
4612 if (exchanger.ghosts_were_set == true)
4613 return;
4614
4615 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4616 reset_ghost_values(*vec[comp], exchanger);
4617 }
4618
4619
4620
4621 //
4622 // zero_vector_region
4623 //
4624
4625 // for block vectors
4626 template <int dim,
4627 typename VectorStruct,
4628 typename Number,
4629 typename VectorizedArrayType,
4630 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4631 * = nullptr>
4632 inline void
4633 zero_vector_region(
4634 const unsigned int range_index,
4635 VectorStruct &vec,
4636 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4637 {
4638 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4639 exchanger.zero_vector_region(range_index, vec.block(i));
4640 }
4641
4642
4643
4644 // for non-block vectors
4645 template <int dim,
4646 typename VectorStruct,
4647 typename Number,
4648 typename VectorizedArrayType,
4649 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4650 * = nullptr>
4651 inline void
4652 zero_vector_region(
4653 const unsigned int range_index,
4654 VectorStruct &vec,
4655 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4656 {
4657 exchanger.zero_vector_region(range_index, vec);
4658 }
4659
4660
4661
4662 // for std::vector of vectors
4663 template <int dim,
4664 typename VectorStruct,
4665 typename Number,
4666 typename VectorizedArrayType>
4667 inline void
4668 zero_vector_region(
4669 const unsigned int range_index,
4670 std::vector<VectorStruct> &vec,
4671 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4672 {
4673 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4674 zero_vector_region(range_index, vec[comp], exchanger);
4675 }
4676
4677
4678
4679 // for std::vector of pointers to vectors
4680 template <int dim,
4681 typename VectorStruct,
4682 typename Number,
4683 typename VectorizedArrayType>
4684 inline void
4685 zero_vector_region(
4686 const unsigned int range_index,
4687 std::vector<VectorStruct *> &vec,
4688 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4689 {
4690 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4691 zero_vector_region(range_index, *vec[comp], exchanger);
4692 }
4693
4694
4695
4696 // Apply a unit matrix operation to constrained DoFs: Default cases where we
4697 // cannot detect a LinearAlgebra::distributed::Vector, we do not do
4698 // anything, else we apply the constraints as a unit operation
4699 template <typename VectorStruct1, typename VectorStruct2>
4700 inline void
4701 apply_operation_to_constrained_dofs(const std::vector<unsigned int> &,
4702 const VectorStruct1 &,
4703 VectorStruct2 &)
4704 {}
4705
4706 template <typename Number>
4707 inline void
4708 apply_operation_to_constrained_dofs(
4709 const std::vector<unsigned int> &constrained_dofs,
4712 {
4713 for (const unsigned int i : constrained_dofs)
4714 dst.local_element(i) = src.local_element(i);
4715 }
4716
4717
4718 namespace MatrixFreeFunctions
4719 {
4720 // struct to select between a const interface and a non-const interface
4721 // for MFWorker
4722 template <typename, typename, typename, typename, bool>
4723 struct InterfaceSelector
4724 {};
4725
4726 // Version for constant functions
4727 template <typename MF,
4728 typename InVector,
4729 typename OutVector,
4730 typename Container>
4731 struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4732 {
4733 using function_type = void (Container::*)(
4734 const MF &,
4735 OutVector &,
4736 const InVector &,
4737 const std::pair<unsigned int, unsigned int> &) const;
4738 };
4739
4740 // Version for non-constant functions
4741 template <typename MF,
4742 typename InVector,
4743 typename OutVector,
4744 typename Container>
4745 struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4746 {
4747 using function_type =
4748 void (Container::*)(const MF &,
4749 OutVector &,
4750 const InVector &,
4751 const std::pair<unsigned int, unsigned int> &);
4752 };
4753 } // namespace MatrixFreeFunctions
4754
4755
4756
4757 // A implementation class for the worker object that runs the various
4758 // operations we want to perform during the matrix-free loop
4759 template <typename MF,
4760 typename InVector,
4761 typename OutVector,
4762 typename Container,
4763 bool is_constant>
4764 class MFWorker : public MFWorkerInterface
4765 {
4766 public:
4767 // An alias to make the arguments further down more readable
4768 using function_type = typename MatrixFreeFunctions::
4769 InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4770 function_type;
4771
4772 // constructor, binds all the arguments to this class
4773 MFWorker(const MF &matrix_free,
4774 const InVector &src,
4775 OutVector &dst,
4776 const bool zero_dst_vector_setting,
4777 const Container &container,
4778 function_type cell_function,
4779 function_type face_function,
4780 function_type boundary_function,
4781 const typename MF::DataAccessOnFaces src_vector_face_access =
4782 MF::DataAccessOnFaces::none,
4783 const typename MF::DataAccessOnFaces dst_vector_face_access =
4784 MF::DataAccessOnFaces::none,
4785 const std::function<void(const unsigned int, const unsigned int)>
4786 &operation_before_loop = {},
4787 const std::function<void(const unsigned int, const unsigned int)>
4788 &operation_after_loop = {},
4789 const unsigned int dof_handler_index_pre_post = 0)
4790 : matrix_free(matrix_free)
4791 , container(const_cast<Container &>(container))
4792 , cell_function(cell_function)
4793 , face_function(face_function)
4794 , boundary_function(boundary_function)
4795 , src(src)
4796 , dst(dst)
4797 , src_data_exchanger(matrix_free,
4798 src_vector_face_access,
4799 n_components(src))
4800 , dst_data_exchanger(matrix_free,
4801 dst_vector_face_access,
4802 n_components(dst))
4803 , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4804 , zero_dst_vector_setting(zero_dst_vector_setting &&
4805 !src_and_dst_are_same)
4806 , operation_before_loop(operation_before_loop)
4807 , operation_after_loop(operation_after_loop)
4808 , dof_handler_index_pre_post(dof_handler_index_pre_post)
4809 {
4810 Assert(!has_ghost_elements(dst),
4811 ExcMessage("The destination vector passed to the matrix-free "
4812 "loop is ghosted. This is not allowed."));
4813 }
4814
4815 // Runs the cell work. If no function is given, nothing is done
4816 virtual void
4817 cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4818 {
4819 if (cell_function != nullptr && cell_range.second > cell_range.first)
4820 for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4821 {
4822 const auto cell_subrange =
4823 matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4824
4825 if (cell_subrange.second <= cell_subrange.first)
4826 continue;
4827
4828 (container.*
4829 cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4830 }
4831 }
4832
4833 virtual void
4834 cell(const unsigned int range_index) override
4835 {
4836 process_range(cell_function,
4837 matrix_free.get_task_info().cell_partition_data_hp_ptr,
4838 matrix_free.get_task_info().cell_partition_data_hp,
4839 range_index);
4840 }
4841
4842 virtual void
4843 face(const unsigned int range_index) override
4844 {
4845 process_range(face_function,
4846 matrix_free.get_task_info().face_partition_data_hp_ptr,
4847 matrix_free.get_task_info().face_partition_data_hp,
4848 range_index);
4849 }
4850
4851 virtual void
4852 boundary(const unsigned int range_index) override
4853 {
4854 process_range(boundary_function,
4855 matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4856 matrix_free.get_task_info().boundary_partition_data_hp,
4857 range_index);
4858 }
4859
4860 private:
4861 void
4862 process_range(const function_type &fu,
4863 const std::vector<unsigned int> &ptr,
4864 const std::vector<unsigned int> &data,
4865 const unsigned int range_index)
4866 {
4867 if (fu == nullptr)
4868 return;
4869
4870 AssertIndexRange(range_index + 1, ptr.size());
4871 for (unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4872 {
4873 AssertIndexRange(2 * i + 1, data.size());
4874 (container.*fu)(matrix_free,
4875 this->dst,
4876 this->src,
4877 std::make_pair(data[2 * i], data[2 * i + 1]));
4878 }
4879 }
4880
4881 public:
4882 // Starts the communication for the update ghost values operation. We
4883 // cannot call this update if ghost and destination are the same because
4884 // that would introduce spurious entries in the destination (there is also
4885 // the problem that reading from a vector that we also write to is usually
4886 // not intended in case there is overlap, but this is up to the
4887 // application code to decide and we cannot catch this case here).
4888 virtual void
4889 vector_update_ghosts_start() override
4890 {
4891 if (!src_and_dst_are_same)
4892 internal::update_ghost_values_start(src, src_data_exchanger);
4893 }
4894
4895 // Finishes the communication for the update ghost values operation
4896 virtual void
4897 vector_update_ghosts_finish() override
4898 {
4899 if (!src_and_dst_are_same)
4900 internal::update_ghost_values_finish(src, src_data_exchanger);
4901 }
4902
4903 // Starts the communication for the vector compress operation
4904 virtual void
4905 vector_compress_start() override
4906 {
4907 internal::compress_start(dst, dst_data_exchanger);
4908 }
4909
4910 // Finishes the communication for the vector compress operation
4911 virtual void
4912 vector_compress_finish() override
4913 {
4914 internal::compress_finish(dst, dst_data_exchanger);
4915 if (!src_and_dst_are_same)
4916 internal::reset_ghost_values(src, src_data_exchanger);
4917 }
4918
4919 // Zeros the given input vector
4920 virtual void
4921 zero_dst_vector_range(const unsigned int range_index) override
4922 {
4923 if (zero_dst_vector_setting)
4924 internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4925 }
4926
4927 virtual void
4928 cell_loop_pre_range(const unsigned int range_index) override
4929 {
4930 if (operation_before_loop)
4931 {
4933 matrix_free.get_dof_info(dof_handler_index_pre_post);
4934 if (range_index == numbers::invalid_unsigned_int)
4935 {
4936 // Case with threaded loop -> currently no overlap implemented
4938 0U,
4939 dof_info.vector_partitioner->locally_owned_size(),
4940 operation_before_loop,
4942 }
4943 else
4944 {
4945 AssertIndexRange(range_index,
4946 dof_info.cell_loop_pre_list_index.size() - 1);
4947 for (unsigned int id =
4948 dof_info.cell_loop_pre_list_index[range_index];
4949 id != dof_info.cell_loop_pre_list_index[range_index + 1];
4950 ++id)
4951 operation_before_loop(dof_info.cell_loop_pre_list[id].first,
4952 dof_info.cell_loop_pre_list[id].second);
4953 }
4954 }
4955 }
4956
4957 virtual void
4958 cell_loop_post_range(const unsigned int range_index) override
4959 {
4960 if (operation_after_loop)
4961 {
4962 // Run unit matrix operation on constrained dofs if we are at the
4963 // last range
4964 const std::vector<unsigned int> &partition_row_index =
4965 matrix_free.get_task_info().partition_row_index;
4966 if (range_index ==
4967 partition_row_index[partition_row_index.size() - 2] - 1)
4968 apply_operation_to_constrained_dofs(
4969 matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
4970 src,
4971 dst);
4972
4974 matrix_free.get_dof_info(dof_handler_index_pre_post);
4975 if (range_index == numbers::invalid_unsigned_int)
4976 {
4977 // Case with threaded loop -> currently no overlap implemented
4979 0U,
4980 dof_info.vector_partitioner->locally_owned_size(),
4981 operation_after_loop,
4983 }
4984 else
4985 {
4986 AssertIndexRange(range_index,
4987 dof_info.cell_loop_post_list_index.size() - 1);
4988 for (unsigned int id =
4989 dof_info.cell_loop_post_list_index[range_index];
4990 id != dof_info.cell_loop_post_list_index[range_index + 1];
4991 ++id)
4992 operation_after_loop(dof_info.cell_loop_post_list[id].first,
4993 dof_info.cell_loop_post_list[id].second);
4994 }
4995 }
4996 }
4997
4998 private:
4999 const MF &matrix_free;
5000 Container &container;
5001 function_type cell_function;
5002 function_type face_function;
5003 function_type boundary_function;
5004
5005 const InVector &src;
5006 OutVector &dst;
5007 VectorDataExchange<MF::dimension,
5008 typename MF::value_type,
5009 typename MF::vectorized_value_type>
5010 src_data_exchanger;
5011 VectorDataExchange<MF::dimension,
5012 typename MF::value_type,
5013 typename MF::vectorized_value_type>
5014 dst_data_exchanger;
5015 const bool src_and_dst_are_same;
5016 const bool zero_dst_vector_setting;
5017 const std::function<void(const unsigned int, const unsigned int)>
5018 operation_before_loop;
5019 const std::function<void(const unsigned int, const unsigned int)>
5020 operation_after_loop;
5021 const unsigned int dof_handler_index_pre_post;
5022 };
5023
5024
5025
5030 template <class MF, typename InVector, typename OutVector>
5031 struct MFClassWrapper
5032 {
5033 using function_type =
5034 std::function<void(const MF &,
5035 OutVector &,
5036 const InVector &,
5037 const std::pair<unsigned int, unsigned int> &)>;
5038
5039 MFClassWrapper(const function_type cell,
5040 const function_type face,
5041 const function_type boundary)
5042 : cell(cell)
5043 , face(face)
5044 , boundary(boundary)
5045 {}
5046
5047 void
5048 cell_integrator(const MF &mf,
5049 OutVector &dst,
5050 const InVector &src,
5051 const std::pair<unsigned int, unsigned int> &range) const
5052 {
5053 if (cell)
5054 cell(mf, dst, src, range);
5055 }
5056
5057 void
5058 face_integrator(const MF &mf,
5059 OutVector &dst,
5060 const InVector &src,
5061 const std::pair<unsigned int, unsigned int> &range) const
5062 {
5063 if (face)
5064 face(mf, dst, src, range);
5065 }
5066
5067 void
5068 boundary_integrator(
5069 const MF &mf,
5070 OutVector &dst,
5071 const InVector &src,
5072 const std::pair<unsigned int, unsigned int> &range) const
5073 {
5074 if (boundary)
5075 boundary(mf, dst, src, range);
5076 }
5077
5078 const function_type cell;
5079 const function_type face;
5080 const function_type boundary;
5081 };
5082
5083} // end of namespace internal
5084
5085
5086
5087template <int dim, typename Number, typename VectorizedArrayType>
5088template <typename OutVector, typename InVector>
5089inline void
5091 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5092 OutVector &,
5093 const InVector &,
5094 const std::pair<unsigned int, unsigned int> &)>
5095 &cell_operation,
5096 OutVector &dst,
5097 const InVector &src,
5098 const bool zero_dst_vector) const
5099{
5100 using Wrapper =
5101 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5102 InVector,
5103 OutVector>;
5104 Wrapper wrap(cell_operation, nullptr, nullptr);
5105 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5106 InVector,
5107 OutVector,
5108 Wrapper,
5109 true>
5110 worker(*this,
5111 src,
5112 dst,
5113 zero_dst_vector,
5114 wrap,
5115 &Wrapper::cell_integrator,
5116 &Wrapper::face_integrator,
5117 &Wrapper::boundary_integrator);
5118
5119 task_info.loop(worker);
5120}
5121
5122
5123
5124template <int dim, typename Number, typename VectorizedArrayType>
5125template <typename OutVector, typename InVector>
5126inline void
5128 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5129 OutVector &,
5130 const InVector &,
5131 const std::pair<unsigned int, unsigned int> &)>
5132 &cell_operation,
5133 OutVector &dst,
5134 const InVector &src,
5135 const std::function<void(const unsigned int, const unsigned int)>
5136 &operation_before_loop,
5137 const std::function<void(const unsigned int, const unsigned int)>
5138 &operation_after_loop,
5139 const unsigned int dof_handler_index_pre_post) const
5140{
5141 using Wrapper =
5142 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5143 InVector,
5144 OutVector>;
5145 Wrapper wrap(cell_operation, nullptr, nullptr);
5146 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5147 InVector,
5148 OutVector,
5149 Wrapper,
5150 true>
5151 worker(*this,
5152 src,
5153 dst,
5154 false,
5155 wrap,
5156 &Wrapper::cell_integrator,
5157 &Wrapper::face_integrator,
5158 &Wrapper::boundary_integrator,
5159 DataAccessOnFaces::none,
5160 DataAccessOnFaces::none,
5161 operation_before_loop,
5162 operation_after_loop,
5163 dof_handler_index_pre_post);
5164
5165 task_info.loop(worker);
5166}
5167
5168
5169
5170template <int dim, typename Number, typename VectorizedArrayType>
5171template <typename OutVector, typename InVector>
5172inline void
5174 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5175 OutVector &,
5176 const InVector &,
5177 const std::pair<unsigned int, unsigned int> &)>
5178 &cell_operation,
5179 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5180 OutVector &,
5181 const InVector &,
5182 const std::pair<unsigned int, unsigned int> &)>
5183 &inner_face_operation,
5184 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5185 OutVector &,
5186 const InVector &,
5187 const std::pair<unsigned int, unsigned int> &)>
5188 &boundary_face_operation,
5189 OutVector &dst,
5190 const InVector &src,
5191 const bool zero_dst_vector,
5192 const DataAccessOnFaces dst_vector_face_access,
5193 const DataAccessOnFaces src_vector_face_access) const
5194{
5195 using Wrapper =
5196 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5197 InVector,
5198 OutVector>;
5199 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5200 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5201 InVector,
5202 OutVector,
5203 Wrapper,
5204 true>
5205 worker(*this,
5206 src,
5207 dst,
5208 zero_dst_vector,
5209 wrap,
5210 &Wrapper::cell_integrator,
5211 &Wrapper::face_integrator,
5212 &Wrapper::boundary_integrator,
5213 src_vector_face_access,
5214 dst_vector_face_access);
5215
5216 task_info.loop(worker);
5217}
5218
5219
5220
5221template <int dim, typename Number, typename VectorizedArrayType>
5222template <typename CLASS, typename OutVector, typename InVector>
5223inline void
5225 void (CLASS::*function_pointer)(
5226 const MatrixFree<dim, Number, VectorizedArrayType> &,
5227 OutVector &,
5228 const InVector &,
5229 const std::pair<unsigned int, unsigned int> &) const,
5230 const CLASS *owning_class,
5231 OutVector &dst,
5232 const InVector &src,
5233 const bool zero_dst_vector) const
5234{
5235 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5236 InVector,
5237 OutVector,
5238 CLASS,
5239 true>
5240 worker(*this,
5241 src,
5242 dst,
5243 zero_dst_vector,
5244 *owning_class,
5245 function_pointer,
5246 nullptr,
5247 nullptr);
5248 task_info.loop(worker);
5249}
5250
5251
5252
5253template <int dim, typename Number, typename VectorizedArrayType>
5254template <typename CLASS, typename OutVector, typename InVector>
5255inline void
5257 void (CLASS::*function_pointer)(
5258 const MatrixFree<dim, Number, VectorizedArrayType> &,
5259 OutVector &,
5260 const InVector &,
5261 const std::pair<unsigned int, unsigned int> &) const,
5262 const CLASS *owning_class,
5263 OutVector &dst,
5264 const InVector &src,
5265 const std::function<void(const unsigned int, const unsigned int)>
5266 &operation_before_loop,
5267 const std::function<void(const unsigned int, const unsigned int)>
5268 &operation_after_loop,
5269 const unsigned int dof_handler_index_pre_post) const
5270{
5271 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5272 InVector,
5273 OutVector,
5274 CLASS,
5275 true>
5276 worker(*this,
5277 src,
5278 dst,
5279 false,
5280 *owning_class,
5281 function_pointer,
5282 nullptr,
5283 nullptr,
5286 operation_before_loop,
5287 operation_after_loop,
5288 dof_handler_index_pre_post);
5289 task_info.loop(worker);
5290}
5291
5292
5293
5294template <int dim, typename Number, typename VectorizedArrayType>
5295template <typename CLASS, typename OutVector, typename InVector>
5296inline void
5298 void (CLASS::*cell_operation)(
5299 const MatrixFree<dim, Number, VectorizedArrayType> &,
5300 OutVector &,
5301 const InVector &,
5302 const std::pair<unsigned int, unsigned int> &) const,
5303 void (CLASS::*inner_face_operation)(
5304 const MatrixFree<dim, Number, VectorizedArrayType> &,
5305 OutVector &,
5306 const InVector &,
5307 const std::pair<unsigned int, unsigned int> &) const,
5308 void (CLASS::*boundary_face_operation)(
5309 const MatrixFree<dim, Number, VectorizedArrayType> &,
5310 OutVector &,
5311 const InVector &,
5312 const std::pair<unsigned int, unsigned int> &) const,
5313 const CLASS *owning_class,
5314 OutVector &dst,
5315 const InVector &src,
5316 const bool zero_dst_vector,
5317 const DataAccessOnFaces dst_vector_face_access,
5318 const DataAccessOnFaces src_vector_face_access) const
5319{
5320 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5321 InVector,
5322 OutVector,
5323 CLASS,
5324 true>
5325 worker(*this,
5326 src,
5327 dst,
5328 zero_dst_vector,
5329 *owning_class,
5330 cell_operation,
5331 inner_face_operation,
5332 boundary_face_operation,
5333 src_vector_face_access,
5334 dst_vector_face_access);
5335 task_info.loop(worker);
5336}
5337
5338
5339
5340template <int dim, typename Number, typename VectorizedArrayType>
5341template <typename CLASS, typename OutVector, typename InVector>
5342inline void
5344 void (CLASS::*function_pointer)(
5345 const MatrixFree<dim, Number, VectorizedArrayType> &,
5346 OutVector &,
5347 const InVector &,
5348 const std::pair<unsigned int, unsigned int> &),
5349 CLASS *owning_class,
5350 OutVector &dst,
5351 const InVector &src,
5352 const bool zero_dst_vector) const
5353{
5354 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5355 InVector,
5356 OutVector,
5357 CLASS,
5358 false>
5359 worker(*this,
5360 src,
5361 dst,
5362 zero_dst_vector,
5363 *owning_class,
5364 function_pointer,
5365 nullptr,
5366 nullptr);
5367 task_info.loop(worker);
5368}
5369
5370
5371
5372template <int dim, typename Number, typename VectorizedArrayType>
5373template <typename CLASS, typename OutVector, typename InVector>
5374inline void
5376 void (CLASS::*function_pointer)(
5377 const MatrixFree<dim, Number, VectorizedArrayType> &,
5378 OutVector &,
5379 const InVector &,
5380 const std::pair<unsigned int, unsigned int> &),
5381 CLASS *owning_class,
5382 OutVector &dst,
5383 const InVector &src,
5384 const std::function<void(const unsigned int, const unsigned int)>
5385 &operation_before_loop,
5386 const std::function<void(const unsigned int, const unsigned int)>
5387 &operation_after_loop,
5388 const unsigned int dof_handler_index_pre_post) const
5389{
5390 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5391 InVector,
5392 OutVector,
5393 CLASS,
5394 false>
5395 worker(*this,
5396 src,
5397 dst,
5398 false,
5399 *owning_class,
5400 function_pointer,
5401 nullptr,
5402 nullptr,
5405 operation_before_loop,
5406 operation_after_loop,
5407 dof_handler_index_pre_post);
5408 task_info.loop(worker);
5409}
5410
5411
5412
5413template <int dim, typename Number, typename VectorizedArrayType>
5414template <typename CLASS, typename OutVector, typename InVector>
5415inline void
5417 void (CLASS::*cell_operation)(
5418 const MatrixFree<dim, Number, VectorizedArrayType> &,
5419 OutVector &,
5420 const InVector &,
5421 const std::pair<unsigned int, unsigned int> &),
5422 void (CLASS::*inner_face_operation)(
5423 const MatrixFree<dim, Number, VectorizedArrayType> &,
5424 OutVector &,
5425 const InVector &,
5426 const std::pair<unsigned int, unsigned int> &),
5427 void (CLASS::*boundary_face_operation)(
5428 const MatrixFree<dim, Number, VectorizedArrayType> &,
5429 OutVector &,
5430 const InVector &,
5431 const std::pair<unsigned int, unsigned int> &),
5432 CLASS *owning_class,
5433 OutVector &dst,
5434 const InVector &src,
5435 const bool zero_dst_vector,
5436 const DataAccessOnFaces dst_vector_face_access,
5437 const DataAccessOnFaces src_vector_face_access) const
5438{
5439 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5440 InVector,
5441 OutVector,
5442 CLASS,
5443 false>
5444 worker(*this,
5445 src,
5446 dst,
5447 zero_dst_vector,
5448 *owning_class,
5449 cell_operation,
5450 inner_face_operation,
5451 boundary_face_operation,
5452 src_vector_face_access,
5453 dst_vector_face_access);
5454 task_info.loop(worker);
5455}
5456
5457
5458
5459template <int dim, typename Number, typename VectorizedArrayType>
5460template <typename OutVector, typename InVector>
5461inline void
5463 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5464 OutVector &,
5465 const InVector &,
5466 const std::pair<unsigned int, unsigned int> &)>
5467 &cell_operation,
5468 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5469 OutVector &,
5470 const InVector &,
5471 const std::pair<unsigned int, unsigned int> &)>
5472 &inner_face_operation,
5473 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5474 OutVector &,
5475 const InVector &,
5476 const std::pair<unsigned int, unsigned int> &)>
5477 &boundary_face_operation,
5478 OutVector &dst,
5479 const InVector &src,
5480 const std::function<void(const unsigned int, const unsigned int)>
5481 &operation_before_loop,
5482 const std::function<void(const unsigned int, const unsigned int)>
5483 &operation_after_loop,
5484 const unsigned int dof_handler_index_pre_post,
5485 const DataAccessOnFaces dst_vector_face_access,
5486 const DataAccessOnFaces src_vector_face_access) const
5487{
5488 using Wrapper =
5489 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5490 InVector,
5491 OutVector>;
5492 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5493 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5494 InVector,
5495 OutVector,
5496 Wrapper,
5497 true>
5498 worker(*this,
5499 src,
5500 dst,
5501 false,
5502 wrap,
5503 &Wrapper::cell_integrator,
5504 &Wrapper::face_integrator,
5505 &Wrapper::boundary_integrator,
5506 src_vector_face_access,
5507 dst_vector_face_access,
5508 operation_before_loop,
5509 operation_after_loop,
5510 dof_handler_index_pre_post);
5511
5512 task_info.loop(worker);
5513}
5514
5515
5516
5517template <int dim, typename Number, typename VectorizedArrayType>
5518template <typename CLASS, typename OutVector, typename InVector>
5519inline void
5521 void (CLASS::*cell_operation)(const MatrixFree &,
5522 OutVector &,
5523 const InVector &,
5524 const std::pair<unsigned int, unsigned int> &)
5525 const,
5526 void (CLASS::*inner_face_operation)(
5527 const MatrixFree &,
5528 OutVector &,
5529 const InVector &,
5530 const std::pair<unsigned int, unsigned int> &) const,
5531 void (CLASS::*boundary_face_operation)(
5532 const MatrixFree &,
5533 OutVector &,
5534 const InVector &,
5535 const std::pair<unsigned int, unsigned int> &) const,
5536 const CLASS *owning_class,
5537 OutVector &dst,
5538 const InVector &src,
5539 const std::function<void(const unsigned int, const unsigned int)>
5540 &operation_before_loop,
5541 const std::function<void(const unsigned int, const unsigned int)>
5542 &operation_after_loop,
5543 const unsigned int dof_handler_index_pre_post,
5544 const DataAccessOnFaces dst_vector_face_access,
5545 const DataAccessOnFaces src_vector_face_access) const
5546{
5547 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5548 InVector,
5549 OutVector,
5550 CLASS,
5551 true>
5552 worker(*this,
5553 src,
5554 dst,
5555 false,
5556 *owning_class,
5557 cell_operation,
5558 inner_face_operation,
5559 boundary_face_operation,
5560 src_vector_face_access,
5561 dst_vector_face_access,
5562 operation_before_loop,
5563 operation_after_loop,
5564 dof_handler_index_pre_post);
5565 task_info.loop(worker);
5566}
5567
5568
5569
5570template <int dim, typename Number, typename VectorizedArrayType>
5571template <typename CLASS, typename OutVector, typename InVector>
5572inline void
5574 void (CLASS::*cell_operation)(const MatrixFree &,
5575 OutVector &,
5576 const InVector &,
5577 const std::pair<unsigned int, unsigned int> &),
5578 void (CLASS::*inner_face_operation)(
5579 const MatrixFree &,
5580 OutVector &,
5581 const InVector &,
5582 const std::pair<unsigned int, unsigned int> &),
5583 void (CLASS::*boundary_face_operation)(
5584 const MatrixFree &,
5585 OutVector &,
5586 const InVector &,
5587 const std::pair<unsigned int, unsigned int> &),
5588 const CLASS *owning_class,
5589 OutVector &dst,
5590 const InVector &src,
5591 const std::function<void(const unsigned int, const unsigned int)>
5592 &operation_before_loop,
5593 const std::function<void(const unsigned int, const unsigned int)>
5594 &operation_after_loop,
5595 const unsigned int dof_handler_index_pre_post,
5596 const DataAccessOnFaces dst_vector_face_access,
5597 const DataAccessOnFaces src_vector_face_access) const
5598{
5599 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5600 InVector,
5601 OutVector,
5602 CLASS,
5603 false>
5604 worker(*this,
5605 src,
5606 dst,
5607 false,
5608 *owning_class,
5609 cell_operation,
5610 inner_face_operation,
5611 boundary_face_operation,
5612 src_vector_face_access,
5613 dst_vector_face_access,
5614 operation_before_loop,
5615 operation_after_loop,
5616 dof_handler_index_pre_post);
5617 task_info.loop(worker);
5618}
5619
5620
5621
5622template <int dim, typename Number, typename VectorizedArrayType>
5623template <typename CLASS, typename OutVector, typename InVector>
5624inline void
5626 void (CLASS::*function_pointer)(
5627 const MatrixFree<dim, Number, VectorizedArrayType> &,
5628 OutVector &,
5629 const InVector &,
5630 const std::pair<unsigned int, unsigned int> &) const,
5631 const CLASS *owning_class,
5632 OutVector &dst,
5633 const InVector &src,
5634 const bool zero_dst_vector,
5635 const DataAccessOnFaces src_vector_face_access) const
5636{
5637 auto src_vector_face_access_temp = src_vector_face_access;
5638 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5639 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5640 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5641 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5642
5643 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5644 InVector,
5645 OutVector,
5646 CLASS,
5647 true>
5648 worker(*this,
5649 src,
5650 dst,
5651 zero_dst_vector,
5652 *owning_class,
5653 function_pointer,
5654 nullptr,
5655 nullptr,
5656 src_vector_face_access_temp,
5658 task_info.loop(worker);
5659}
5660
5661
5662
5663template <int dim, typename Number, typename VectorizedArrayType>
5664template <typename CLASS, typename OutVector, typename InVector>
5665inline void
5667 void (CLASS::*function_pointer)(
5668 const MatrixFree<dim, Number, VectorizedArrayType> &,
5669 OutVector &,
5670 const InVector &,
5671 const std::pair<unsigned int, unsigned int> &),
5672 CLASS *owning_class,
5673 OutVector &dst,
5674 const InVector &src,
5675 const bool zero_dst_vector,
5676 const DataAccessOnFaces src_vector_face_access) const
5677{
5678 auto src_vector_face_access_temp = src_vector_face_access;
5679 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5680 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5681 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5682 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5683
5684 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5685 InVector,
5686 OutVector,
5687 CLASS,
5688 false>
5689 worker(*this,
5690 src,
5691 dst,
5692 zero_dst_vector,
5693 *owning_class,
5694 function_pointer,
5695 nullptr,
5696 nullptr,
5697 src_vector_face_access_temp,
5699 task_info.loop(worker);
5700}
5701
5702
5703
5704template <int dim, typename Number, typename VectorizedArrayType>
5705template <typename OutVector, typename InVector>
5706inline void
5708 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5709 OutVector &,
5710 const InVector &,
5711 const std::pair<unsigned int, unsigned int> &)>
5712 &cell_operation,
5713 OutVector &dst,
5714 const InVector &src,
5715 const bool zero_dst_vector,
5716 const DataAccessOnFaces src_vector_face_access) const
5717{
5718 auto src_vector_face_access_temp = src_vector_face_access;
5719 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5720 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5721 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5722 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5723
5724 using Wrapper =
5725 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5726 InVector,
5727 OutVector>;
5728 Wrapper wrap(cell_operation, nullptr, nullptr);
5729
5730 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5731 InVector,
5732 OutVector,
5733 Wrapper,
5734 true>
5735 worker(*this,
5736 src,
5737 dst,
5738 zero_dst_vector,
5739 wrap,
5740 &Wrapper::cell_integrator,
5741 &Wrapper::face_integrator,
5742 &Wrapper::boundary_integrator,
5743 src_vector_face_access_temp,
5744 DataAccessOnFaces::none);
5745 task_info.loop(worker);
5746}
5747
5748
5749#endif // ifndef DOXYGEN
5750
5751
5752
5754
5755#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