Reference documentation for deal.II version GIT relicensing-468-ge3ce94fd06 2024-04-23 15:40:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
matrix_free.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2012 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
16#ifndef dealii_matrix_free_h
17#define dealii_matrix_free_h
18
19#include <deal.II/base/config.h>
20
27
29
30#include <deal.II/fe/fe.h>
31#include <deal.II/fe/mapping.h>
32
35
40
47
48#include <cstdlib>
49#include <limits>
50#include <list>
51#include <memory>
52
53
55
56
57
109template <int dim,
110 typename Number = double,
111 typename VectorizedArrayType = VectorizedArray<Number>>
113{
114 static_assert(
115 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
116 "Type of Number and of VectorizedArrayType do not match.");
117
118public:
123 using value_type = Number;
124 using vectorized_value_type = VectorizedArrayType;
125
129 static constexpr unsigned int dimension = dim;
130
184 {
189
217
223 const unsigned int tasks_block_size = 0,
229 const unsigned int mg_level = numbers::invalid_unsigned_int,
230 const bool store_plain_indices = true,
231 const bool initialize_indices = true,
232 const bool initialize_mapping = true,
233 const bool overlap_communication_computation = true,
234 const bool hold_all_faces_to_owned_cells = false,
235 const bool cell_vectorization_categories_strict = false,
236 const bool allow_ghosted_vectors_in_loops = true)
252 , communicator_sm(MPI_COMM_SELF)
253 {}
254
280
310
348
358 unsigned int tasks_block_size;
359
374
394
414
442
451 unsigned int mg_level;
452
460
471
480
494
503
525 std::vector<unsigned int> cell_vectorization_category;
526
537
549
554 };
555
564
569
573 ~MatrixFree() override = default;
574
587 template <typename QuadratureType, typename number2, typename MappingType>
588 void
589 reinit(const MappingType &mapping,
590 const DoFHandler<dim> &dof_handler,
591 const AffineConstraints<number2> &constraint,
592 const QuadratureType &quad,
593 const AdditionalData &additional_data = AdditionalData());
594
616 template <typename QuadratureType, typename number2, typename MappingType>
617 void
618 reinit(const MappingType &mapping,
619 const std::vector<const DoFHandler<dim> *> &dof_handler,
620 const std::vector<const AffineConstraints<number2> *> &constraint,
621 const std::vector<QuadratureType> &quad,
622 const AdditionalData &additional_data = AdditionalData());
623
631 template <typename QuadratureType, typename number2, typename MappingType>
632 void
633 reinit(const MappingType &mapping,
634 const std::vector<const DoFHandler<dim> *> &dof_handler,
635 const std::vector<const AffineConstraints<number2> *> &constraint,
636 const QuadratureType &quad,
637 const AdditionalData &additional_data = AdditionalData());
638
644 void
646 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free_base);
647
657 void
659
663 void
664 update_mapping(const std::shared_ptr<hp::MappingCollection<dim>> &mapping);
665
670 void
672
687 {
694 none,
695
706 values,
707
717
728 gradients,
729
739
747 };
748
793 template <typename OutVector, typename InVector>
794 void
795 cell_loop(const std::function<void(
797 OutVector &,
798 const InVector &,
799 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
800 OutVector &dst,
801 const InVector &src,
802 const bool zero_dst_vector = false) const;
803
850 template <typename CLASS, typename OutVector, typename InVector>
851 void
852 cell_loop(void (CLASS::*cell_operation)(
853 const MatrixFree &,
854 OutVector &,
855 const InVector &,
856 const std::pair<unsigned int, unsigned int> &) const,
857 const CLASS *owning_class,
858 OutVector &dst,
859 const InVector &src,
860 const bool zero_dst_vector = false) const;
861
865 template <typename CLASS, typename OutVector, typename InVector>
866 void
867 cell_loop(void (CLASS::*cell_operation)(
868 const MatrixFree &,
869 OutVector &,
870 const InVector &,
871 const std::pair<unsigned int, unsigned int> &),
872 CLASS *owning_class,
873 OutVector &dst,
874 const InVector &src,
875 const bool zero_dst_vector = false) const;
876
961 template <typename CLASS, typename OutVector, typename InVector>
962 void
963 cell_loop(void (CLASS::*cell_operation)(
964 const MatrixFree &,
965 OutVector &,
966 const InVector &,
967 const std::pair<unsigned int, unsigned int> &) const,
968 const CLASS *owning_class,
969 OutVector &dst,
970 const InVector &src,
971 const std::function<void(const unsigned int, const unsigned int)>
972 &operation_before_loop,
973 const std::function<void(const unsigned int, const unsigned int)>
974 &operation_after_loop,
975 const unsigned int dof_handler_index_pre_post = 0) const;
976
980 template <typename CLASS, typename OutVector, typename InVector>
981 void
982 cell_loop(void (CLASS::*cell_operation)(
983 const MatrixFree &,
984 OutVector &,
985 const InVector &,
986 const std::pair<unsigned int, unsigned int> &),
987 CLASS *owning_class,
988 OutVector &dst,
989 const InVector &src,
990 const std::function<void(const unsigned int, const unsigned int)>
991 &operation_before_loop,
992 const std::function<void(const unsigned int, const unsigned int)>
993 &operation_after_loop,
994 const unsigned int dof_handler_index_pre_post = 0) const;
995
1000 template <typename OutVector, typename InVector>
1001 void
1002 cell_loop(const std::function<void(
1004 OutVector &,
1005 const InVector &,
1006 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1007 OutVector &dst,
1008 const InVector &src,
1009 const std::function<void(const unsigned int, const unsigned int)>
1010 &operation_before_loop,
1011 const std::function<void(const unsigned int, const unsigned int)>
1012 &operation_after_loop,
1013 const unsigned int dof_handler_index_pre_post = 0) const;
1014
1090 template <typename OutVector, typename InVector>
1091 void
1093 const std::function<
1095 OutVector &,
1096 const InVector &,
1097 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1098 const std::function<void(
1100 OutVector &,
1101 const InVector &,
1102 const std::pair<unsigned int, unsigned int> &)> &inner_face_operation,
1103 const std::function<void(
1105 OutVector &,
1106 const InVector &,
1107 const std::pair<unsigned int, unsigned int> &)> &boundary_face_operation,
1108 OutVector &dst,
1109 const InVector &src,
1110 const bool zero_dst_vector = false,
1111 const DataAccessOnFaces dst_vector_face_access =
1113 const DataAccessOnFaces src_vector_face_access =
1115
1204 template <typename CLASS, typename OutVector, typename InVector>
1205 void
1206 loop(void (CLASS::*cell_operation)(
1207 const MatrixFree &,
1208 OutVector &,
1209 const InVector &,
1210 const std::pair<unsigned int, unsigned int> &) const,
1211 void (CLASS::*inner_face_operation)(
1212 const MatrixFree &,
1213 OutVector &,
1214 const InVector &,
1215 const std::pair<unsigned int, unsigned int> &) const,
1216 void (CLASS::*boundary_face_operation)(
1217 const MatrixFree &,
1218 OutVector &,
1219 const InVector &,
1220 const std::pair<unsigned int, unsigned int> &) const,
1221 const CLASS *owning_class,
1222 OutVector &dst,
1223 const InVector &src,
1224 const bool zero_dst_vector = false,
1225 const DataAccessOnFaces dst_vector_face_access =
1227 const DataAccessOnFaces src_vector_face_access =
1229
1233 template <typename CLASS, typename OutVector, typename InVector>
1234 void
1235 loop(void (CLASS::*cell_operation)(
1236 const MatrixFree &,
1237 OutVector &,
1238 const InVector &,
1239 const std::pair<unsigned int, unsigned int> &),
1240 void (CLASS::*inner_face_operation)(
1241 const MatrixFree &,
1242 OutVector &,
1243 const InVector &,
1244 const std::pair<unsigned int, unsigned int> &),
1245 void (CLASS::*boundary_face_operation)(
1246 const MatrixFree &,
1247 OutVector &,
1248 const InVector &,
1249 const std::pair<unsigned int, unsigned int> &),
1250 CLASS *owning_class,
1251 OutVector &dst,
1252 const InVector &src,
1253 const bool zero_dst_vector = false,
1254 const DataAccessOnFaces dst_vector_face_access =
1256 const DataAccessOnFaces src_vector_face_access =
1258
1380 template <typename CLASS, typename OutVector, typename InVector>
1381 void
1382 loop(void (CLASS::*cell_operation)(
1383 const MatrixFree &,
1384 OutVector &,
1385 const InVector &,
1386 const std::pair<unsigned int, unsigned int> &) const,
1387 void (CLASS::*inner_face_operation)(
1388 const MatrixFree &,
1389 OutVector &,
1390 const InVector &,
1391 const std::pair<unsigned int, unsigned int> &) const,
1392 void (CLASS::*boundary_face_operation)(
1393 const MatrixFree &,
1394 OutVector &,
1395 const InVector &,
1396 const std::pair<unsigned int, unsigned int> &) const,
1397 const CLASS *owning_class,
1398 OutVector &dst,
1399 const InVector &src,
1400 const std::function<void(const unsigned int, const unsigned int)>
1401 &operation_before_loop,
1402 const std::function<void(const unsigned int, const unsigned int)>
1403 &operation_after_loop,
1404 const unsigned int dof_handler_index_pre_post = 0,
1405 const DataAccessOnFaces dst_vector_face_access =
1407 const DataAccessOnFaces src_vector_face_access =
1409
1413 template <typename CLASS, typename OutVector, typename InVector>
1414 void
1415 loop(void (CLASS::*cell_operation)(
1416 const MatrixFree &,
1417 OutVector &,
1418 const InVector &,
1419 const std::pair<unsigned int, unsigned int> &),
1420 void (CLASS::*inner_face_operation)(
1421 const MatrixFree &,
1422 OutVector &,
1423 const InVector &,
1424 const std::pair<unsigned int, unsigned int> &),
1425 void (CLASS::*boundary_face_operation)(
1426 const MatrixFree &,
1427 OutVector &,
1428 const InVector &,
1429 const std::pair<unsigned int, unsigned int> &),
1430 const CLASS *owning_class,
1431 OutVector &dst,
1432 const InVector &src,
1433 const std::function<void(const unsigned int, const unsigned int)>
1434 &operation_before_loop,
1435 const std::function<void(const unsigned int, const unsigned int)>
1436 &operation_after_loop,
1437 const unsigned int dof_handler_index_pre_post = 0,
1438 const DataAccessOnFaces dst_vector_face_access =
1440 const DataAccessOnFaces src_vector_face_access =
1442
1448 template <typename OutVector, typename InVector>
1449 void
1451 const std::function<
1453 OutVector &,
1454 const InVector &,
1455 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1456 const std::function<void(
1458 OutVector &,
1459 const InVector &,
1460 const std::pair<unsigned int, unsigned int> &)> &inner_face_operation,
1461 const std::function<void(
1463 OutVector &,
1464 const InVector &,
1465 const std::pair<unsigned int, unsigned int> &)> &boundary_face_operation,
1466 OutVector &dst,
1467 const InVector &src,
1468 const std::function<void(const unsigned int, const unsigned int)>
1469 &operation_before_loop,
1470 const std::function<void(const unsigned int, const unsigned int)>
1471 &operation_after_loop,
1472 const unsigned int dof_handler_index_pre_post = 0,
1473 const DataAccessOnFaces dst_vector_face_access =
1475 const DataAccessOnFaces src_vector_face_access =
1477
1542 template <typename CLASS, typename OutVector, typename InVector>
1543 void
1544 loop_cell_centric(void (CLASS::*cell_operation)(
1545 const MatrixFree &,
1546 OutVector &,
1547 const InVector &,
1548 const std::pair<unsigned int, unsigned int> &) const,
1549 const CLASS *owning_class,
1550 OutVector &dst,
1551 const InVector &src,
1552 const bool zero_dst_vector = false,
1553 const DataAccessOnFaces src_vector_face_access =
1555
1559 template <typename CLASS, typename OutVector, typename InVector>
1560 void
1561 loop_cell_centric(void (CLASS::*cell_operation)(
1562 const MatrixFree &,
1563 OutVector &,
1564 const InVector &,
1565 const std::pair<unsigned int, unsigned int> &),
1566 CLASS *owning_class,
1567 OutVector &dst,
1568 const InVector &src,
1569 const bool zero_dst_vector = false,
1570 const DataAccessOnFaces src_vector_face_access =
1572
1576 template <typename OutVector, typename InVector>
1577 void
1579 const std::function<void(const MatrixFree &,
1580 OutVector &,
1581 const InVector &,
1582 const std::pair<unsigned int, unsigned int> &)>
1583 &cell_operation,
1584 OutVector &dst,
1585 const InVector &src,
1586 const bool zero_dst_vector = false,
1587 const DataAccessOnFaces src_vector_face_access =
1589
1597 std::pair<unsigned int, unsigned int>
1598 create_cell_subrange_hp(const std::pair<unsigned int, unsigned int> &range,
1599 const unsigned int fe_degree,
1600 const unsigned int dof_handler_index = 0) const;
1601
1608 std::pair<unsigned int, unsigned int>
1610 const std::pair<unsigned int, unsigned int> &range,
1611 const unsigned int fe_index,
1612 const unsigned int dof_handler_index = 0) const;
1613
1617 unsigned int
1619
1623 unsigned int
1625 const std::pair<unsigned int, unsigned int> range) const;
1626
1630 unsigned int
1631 get_face_active_fe_index(const std::pair<unsigned int, unsigned int> range,
1632 const bool is_interior_face = true) const;
1633
1645 template <typename T>
1646 void
1648
1654 template <typename T>
1655 void
1657
1671 template <typename VectorType>
1672 void
1673 initialize_dof_vector(VectorType &vec,
1674 const unsigned int dof_handler_index = 0) const;
1675
1693 template <typename Number2>
1694 void
1696 const unsigned int dof_handler_index = 0) const;
1697
1708 const std::shared_ptr<const Utilities::MPI::Partitioner> &
1709 get_vector_partitioner(const unsigned int dof_handler_index = 0) const;
1710
1714 const IndexSet &
1715 get_locally_owned_set(const unsigned int dof_handler_index = 0) const;
1716
1720 const IndexSet &
1721 get_ghost_set(const unsigned int dof_handler_index = 0) const;
1722
1732 const std::vector<unsigned int> &
1733 get_constrained_dofs(const unsigned int dof_handler_index = 0) const;
1734
1745 void
1746 renumber_dofs(std::vector<types::global_dof_index> &renumbering,
1747 const unsigned int dof_handler_index = 0);
1748
1758 template <int spacedim>
1759 static bool
1761
1765 unsigned int
1767
1772 unsigned int
1773 n_base_elements(const unsigned int dof_handler_index) const;
1774
1782 unsigned int
1784
1794 unsigned int
1796
1803 unsigned int
1805
1815 unsigned int
1817
1828 unsigned int
1830
1835 unsigned int
1837
1849 get_boundary_id(const unsigned int face_batch_index) const;
1850
1855 std::array<types::boundary_id, VectorizedArrayType::size()>
1856 get_faces_by_cells_boundary_id(const unsigned int cell_batch_index,
1857 const unsigned int face_number) const;
1858
1863 const DoFHandler<dim> &
1864 get_dof_handler(const unsigned int dof_handler_index = 0) const;
1865
1873 get_affine_constraints(const unsigned int dof_handler_index = 0) const;
1874
1888 get_cell_iterator(const unsigned int cell_batch_index,
1889 const unsigned int lane_index,
1890 const unsigned int dof_handler_index = 0) const;
1891
1897 std::pair<int, int>
1898 get_cell_level_and_index(const unsigned int cell_batch_index,
1899 const unsigned int lane_index) const;
1900
1907 unsigned int
1909 const typename Triangulation<dim>::cell_iterator &cell) const;
1910
1923 std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>
1924 get_face_iterator(const unsigned int face_batch_index,
1925 const unsigned int lane_index,
1926 const bool interior = true,
1927 const unsigned int fe_component = 0) const;
1928
1941 bool
1942 at_irregular_cell(const unsigned int cell_batch_index) const;
1943
1953 unsigned int
1954 n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const;
1955
1965 unsigned int
1966 n_active_entries_per_face_batch(const unsigned int face_batch_index) const;
1967
1971 unsigned int
1972 get_dofs_per_cell(const unsigned int dof_handler_index = 0,
1973 const unsigned int hp_active_fe_index = 0) const;
1974
1978 unsigned int
1979 get_n_q_points(const unsigned int quad_index = 0,
1980 const unsigned int hp_active_fe_index = 0) const;
1981
1986 unsigned int
1987 get_dofs_per_face(const unsigned int dof_handler_index = 0,
1988 const unsigned int hp_active_fe_index = 0) const;
1989
1994 unsigned int
1995 get_n_q_points_face(const unsigned int quad_index = 0,
1996 const unsigned int hp_active_fe_index = 0) const;
1997
2001 const Quadrature<dim> &
2002 get_quadrature(const unsigned int quad_index = 0,
2003 const unsigned int hp_active_fe_index = 0) const;
2004
2008 const Quadrature<dim - 1> &
2009 get_face_quadrature(const unsigned int quad_index = 0,
2010 const unsigned int hp_active_fe_index = 0) const;
2011
2024 unsigned int
2026 const std::pair<unsigned int, unsigned int> cell_batch_range) const;
2027
2032 std::pair<unsigned int, unsigned int>
2034 const std::pair<unsigned int, unsigned int> face_batch_range) const;
2035
2047 unsigned int
2048 get_cell_category(const unsigned int cell_batch_index) const;
2049
2054 std::pair<unsigned int, unsigned int>
2055 get_face_category(const unsigned int face_batch_index) const;
2056
2060 bool
2062
2067 bool
2069
2074 unsigned int
2076
2081 std::size_t
2083
2088 template <typename StreamType>
2089 void
2090 print_memory_consumption(StreamType &out) const;
2091
2096 void
2097 print(std::ostream &out) const;
2098
2112
2113 /*
2114 * Return geometry-dependent information on the cells.
2115 */
2116 const internal::MatrixFreeFunctions::
2117 MappingInfo<dim, Number, VectorizedArrayType> &
2119
2124 get_dof_info(const unsigned int dof_handler_index_component = 0) const;
2125
2129 unsigned int
2131
2136 const Number *
2137 constraint_pool_begin(const unsigned int pool_index) const;
2138
2144 const Number *
2145 constraint_pool_end(const unsigned int pool_index) const;
2146
2151 get_shape_info(const unsigned int dof_handler_index_component = 0,
2152 const unsigned int quad_index = 0,
2153 const unsigned int fe_base_element = 0,
2154 const unsigned int hp_active_fe_index = 0,
2155 const unsigned int hp_active_quad_index = 0) const;
2156
2161 VectorizedArrayType::size()> &
2162 get_face_info(const unsigned int face_batch_index) const;
2163
2164
2172
2188
2192 void
2194
2206
2210 void
2212 const AlignedVector<Number> *memory) const;
2213
2216private:
2221 template <typename number2, int q_dim>
2222 void
2224 const std::shared_ptr<hp::MappingCollection<dim>> &mapping,
2225 const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
2226 const std::vector<const AffineConstraints<number2> *> &constraint,
2227 const std::vector<IndexSet> &locally_owned_set,
2228 const std::vector<hp::QCollection<q_dim>> &quad,
2229 const AdditionalData &additional_data);
2230
2237 template <typename number2>
2238 void
2240 const std::vector<const AffineConstraints<number2> *> &constraint,
2241 const std::vector<IndexSet> &locally_owned_set,
2242 const AdditionalData &additional_data);
2243
2247 void
2249 const std::vector<const DoFHandler<dim, dim> *> &dof_handlers,
2250 const AdditionalData &additional_data);
2251
2255 std::vector<SmartPointer<const DoFHandler<dim>>> dof_handlers;
2256
2263 std::vector<SmartPointer<const AffineConstraints<Number>>> affine_constraints;
2264
2269 std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
2270
2277 std::vector<Number> constraint_pool_data;
2278
2283 std::vector<unsigned int> constraint_pool_row_index;
2284
2291
2297
2304 std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
2305
2310 std::vector<unsigned int> mf_cell_indices;
2311
2319
2326
2331 internal::MatrixFreeFunctions::FaceInfo<VectorizedArrayType::size()>
2333
2338
2343
2352 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2354
2359 mutable std::list<std::pair<bool, AlignedVector<Number>>>
2361
2365 unsigned int mg_level;
2366};
2367
2368
2369
2370/*----------------------- Inline functions ----------------------------------*/
2371
2372#ifndef DOXYGEN
2373
2374
2375
2376template <int dim, typename Number, typename VectorizedArrayType>
2377template <typename T>
2378inline void
2380 AlignedVector<T> &vec) const
2381{
2382 vec.resize(this->n_cell_batches() + this->n_ghost_cell_batches());
2383}
2384
2385
2386
2387template <int dim, typename Number, typename VectorizedArrayType>
2388template <typename T>
2389inline void
2391 AlignedVector<T> &vec) const
2392{
2393 vec.resize(this->n_inner_face_batches() + this->n_boundary_face_batches() +
2394 this->n_ghost_inner_face_batches());
2395}
2396
2397
2398
2399template <int dim, typename Number, typename VectorizedArrayType>
2400template <typename VectorType>
2401inline void
2403 VectorType &vec,
2404 const unsigned int comp) const
2405{
2406 static_assert(IsBlockVector<VectorType>::value == false,
2407 "This function is not supported for block vectors.");
2408
2409 Assert(task_info.n_procs == 1,
2410 ExcMessage("This function can only be used in serial."));
2411
2412 AssertIndexRange(comp, n_components());
2413 vec.reinit(dof_info[comp].vector_partitioner->size());
2414}
2415
2416
2417
2418template <int dim, typename Number, typename VectorizedArrayType>
2419template <typename Number2>
2420inline void
2423 const unsigned int comp) const
2424{
2425 AssertIndexRange(comp, n_components());
2426 vec.reinit(dof_info[comp].vector_partitioner, task_info.communicator_sm);
2427}
2428
2429
2430
2431template <int dim, typename Number, typename VectorizedArrayType>
2432inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2434 const unsigned int comp) const
2435{
2436 AssertIndexRange(comp, n_components());
2437 return dof_info[comp].vector_partitioner;
2438}
2439
2440
2441
2442template <int dim, typename Number, typename VectorizedArrayType>
2443inline const std::vector<unsigned int> &
2445 const unsigned int comp) const
2446{
2447 AssertIndexRange(comp, n_components());
2448 return dof_info[comp].constrained_dofs;
2449}
2450
2451
2452
2453template <int dim, typename Number, typename VectorizedArrayType>
2454inline unsigned int
2456{
2457 AssertDimension(dof_handlers.size(), dof_info.size());
2458 return dof_handlers.size();
2459}
2460
2461
2462
2463template <int dim, typename Number, typename VectorizedArrayType>
2464inline unsigned int
2466 const unsigned int dof_no) const
2467{
2468 AssertDimension(dof_handlers.size(), dof_info.size());
2469 AssertIndexRange(dof_no, dof_handlers.size());
2470 return dof_handlers[dof_no]->get_fe().n_base_elements();
2471}
2472
2473
2474
2475template <int dim, typename Number, typename VectorizedArrayType>
2478{
2479 return task_info;
2480}
2481
2482
2483
2484template <int dim, typename Number, typename VectorizedArrayType>
2485inline unsigned int
2487{
2488 return task_info.n_active_cells;
2489}
2490
2491
2492
2493template <int dim, typename Number, typename VectorizedArrayType>
2494inline unsigned int
2496{
2497 return *(task_info.cell_partition_data.end() - 2);
2498}
2499
2500
2501
2502template <int dim, typename Number, typename VectorizedArrayType>
2503inline unsigned int
2505{
2506 return *(task_info.cell_partition_data.end() - 1) -
2507 *(task_info.cell_partition_data.end() - 2);
2508}
2509
2510
2511
2512template <int dim, typename Number, typename VectorizedArrayType>
2513inline unsigned int
2515{
2516 if (task_info.face_partition_data.empty())
2517 return 0;
2518 return task_info.face_partition_data.back();
2519}
2520
2521
2522
2523template <int dim, typename Number, typename VectorizedArrayType>
2524inline unsigned int
2526{
2527 if (task_info.face_partition_data.empty())
2528 return 0;
2529 return task_info.boundary_partition_data.back() -
2530 task_info.face_partition_data.back();
2531}
2532
2533
2534
2535template <int dim, typename Number, typename VectorizedArrayType>
2536inline unsigned int
2538{
2539 if (task_info.face_partition_data.empty())
2540 return 0;
2541 return face_info.faces.size() - task_info.boundary_partition_data.back();
2542}
2543
2544
2545
2546template <int dim, typename Number, typename VectorizedArrayType>
2547inline types::boundary_id
2549 const unsigned int face_batch_index) const
2550{
2551 Assert(face_batch_index >= task_info.boundary_partition_data[0] &&
2552 face_batch_index < task_info.boundary_partition_data.back(),
2553 ExcIndexRange(face_batch_index,
2554 task_info.boundary_partition_data[0],
2555 task_info.boundary_partition_data.back()));
2556 return types::boundary_id(face_info.faces[face_batch_index].exterior_face_no);
2557}
2558
2559
2560
2561template <int dim, typename Number, typename VectorizedArrayType>
2562inline std::array<types::boundary_id, VectorizedArrayType::size()>
2564 const unsigned int cell_batch_index,
2565 const unsigned int face_number) const
2566{
2567 AssertIndexRange(cell_batch_index, n_cell_batches());
2569 Assert(face_info.cell_and_face_boundary_id.size(0) >= n_cell_batches(),
2571 std::array<types::boundary_id, VectorizedArrayType::size()> result;
2572 result.fill(numbers::invalid_boundary_id);
2573 for (unsigned int v = 0;
2574 v < n_active_entries_per_cell_batch(cell_batch_index);
2575 ++v)
2576 result[v] =
2577 face_info.cell_and_face_boundary_id(cell_batch_index, face_number, v);
2578 return result;
2579}
2580
2581
2582
2583template <int dim, typename Number, typename VectorizedArrayType>
2584inline const internal::MatrixFreeFunctions::
2585 MappingInfo<dim, Number, VectorizedArrayType> &
2587{
2588 return mapping_info;
2589}
2590
2591
2592
2593template <int dim, typename Number, typename VectorizedArrayType>
2596 const unsigned int dof_index) const
2597{
2598 AssertIndexRange(dof_index, n_components());
2599 return dof_info[dof_index];
2600}
2601
2602
2603
2604template <int dim, typename Number, typename VectorizedArrayType>
2605inline unsigned int
2607{
2608 return constraint_pool_row_index.size() - 1;
2609}
2610
2611
2612
2613template <int dim, typename Number, typename VectorizedArrayType>
2614inline const Number *
2616 const unsigned int row) const
2617{
2618 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2619 return constraint_pool_data.empty() ?
2620 nullptr :
2621 constraint_pool_data.data() + constraint_pool_row_index[row];
2622}
2623
2624
2625
2626template <int dim, typename Number, typename VectorizedArrayType>
2627inline const Number *
2629 const unsigned int row) const
2630{
2631 AssertIndexRange(row, constraint_pool_row_index.size() - 1);
2632 return constraint_pool_data.empty() ?
2633 nullptr :
2634 constraint_pool_data.data() + constraint_pool_row_index[row + 1];
2635}
2636
2637
2638
2639template <int dim, typename Number, typename VectorizedArrayType>
2640inline std::pair<unsigned int, unsigned int>
2642 const std::pair<unsigned int, unsigned int> &range,
2643 const unsigned int degree,
2644 const unsigned int dof_handler_component) const
2645{
2646 if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2647 {
2649 dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2651 dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2652 if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2653 return range;
2654 else
2655 return {range.second, range.second};
2656 }
2657
2658 const unsigned int fe_index =
2659 dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2660 if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2661 return {range.second, range.second};
2662 else
2663 return create_cell_subrange_hp_by_index(range,
2664 fe_index,
2665 dof_handler_component);
2666}
2667
2668
2669
2670template <int dim, typename Number, typename VectorizedArrayType>
2671inline bool
2673 const unsigned int cell_batch_index) const
2674{
2675 AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2676 return VectorizedArrayType::size() > 1 &&
2677 cell_level_index[(cell_batch_index + 1) * VectorizedArrayType::size() -
2678 1] == cell_level_index[(cell_batch_index + 1) *
2679 VectorizedArrayType::size() -
2680 2];
2681}
2682
2683
2684
2685template <int dim, typename Number, typename VectorizedArrayType>
2686unsigned int
2688{
2689 return shape_info.size(2);
2690}
2691
2692
2693template <int dim, typename Number, typename VectorizedArrayType>
2694unsigned int
2696 const std::pair<unsigned int, unsigned int> range) const
2697{
2698 const auto &fe_indices = dof_info[0].cell_active_fe_index;
2699
2700 if (fe_indices.empty() == true ||
2701 dof_handlers[0]->get_fe_collection().size() == 1)
2702 return 0;
2703
2704 const auto index = fe_indices[range.first];
2705
2706 for (unsigned int i = range.first; i < range.second; ++i)
2707 AssertDimension(index, fe_indices[i]);
2708
2709 return index;
2710}
2711
2712
2713
2714template <int dim, typename Number, typename VectorizedArrayType>
2715unsigned int
2717 const std::pair<unsigned int, unsigned int> range,
2718 const bool is_interior_face) const
2719{
2720 const auto &fe_indices = dof_info[0].cell_active_fe_index;
2721
2722 if (fe_indices.empty() == true)
2723 return 0;
2724
2725 if (is_interior_face)
2726 {
2727 const unsigned int index =
2728 fe_indices[face_info.faces[range.first].cells_interior[0] /
2729 VectorizedArrayType::size()];
2730
2731 for (unsigned int i = range.first; i < range.second; ++i)
2732 AssertDimension(index,
2733 fe_indices[face_info.faces[i].cells_interior[0] /
2734 VectorizedArrayType::size()]);
2735
2736 return index;
2737 }
2738 else
2739 {
2740 const unsigned int index =
2741 fe_indices[face_info.faces[range.first].cells_exterior[0] /
2742 VectorizedArrayType::size()];
2743
2744 for (unsigned int i = range.first; i < range.second; ++i)
2745 AssertDimension(index,
2746 fe_indices[face_info.faces[i].cells_exterior[0] /
2747 VectorizedArrayType::size()]);
2748
2749 return index;
2750 }
2751}
2752
2753
2754
2755template <int dim, typename Number, typename VectorizedArrayType>
2756inline unsigned int
2758 const unsigned int cell_batch_index) const
2759{
2760 Assert(!dof_info.empty(), ExcNotInitialized());
2761 AssertIndexRange(cell_batch_index, task_info.cell_partition_data.back());
2762 const std::vector<unsigned char> &n_lanes_filled =
2763 dof_info[0].n_vectorization_lanes_filled
2765 AssertIndexRange(cell_batch_index, n_lanes_filled.size());
2766
2767 return n_lanes_filled[cell_batch_index];
2768}
2769
2770
2771
2772template <int dim, typename Number, typename VectorizedArrayType>
2773inline unsigned int
2775 const unsigned int face_batch_index) const
2776{
2777 AssertIndexRange(face_batch_index, face_info.faces.size());
2778 Assert(!dof_info.empty(), ExcNotInitialized());
2779 const std::vector<unsigned char> &n_lanes_filled =
2780 dof_info[0].n_vectorization_lanes_filled
2782 AssertIndexRange(face_batch_index, n_lanes_filled.size());
2783 return n_lanes_filled[face_batch_index];
2784}
2785
2786
2787
2788template <int dim, typename Number, typename VectorizedArrayType>
2789inline unsigned int
2791 const unsigned int dof_handler_index,
2792 const unsigned int active_fe_index) const
2793{
2794 return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2795}
2796
2797
2798
2799template <int dim, typename Number, typename VectorizedArrayType>
2800inline unsigned int
2802 const unsigned int quad_index,
2803 const unsigned int active_fe_index) const
2804{
2805 AssertIndexRange(quad_index, mapping_info.cell_data.size());
2806 return mapping_info.cell_data[quad_index]
2807 .descriptor[active_fe_index]
2808 .n_q_points;
2809}
2810
2811
2812
2813template <int dim, typename Number, typename VectorizedArrayType>
2814inline unsigned int
2816 const unsigned int dof_handler_index,
2817 const unsigned int active_fe_index) const
2818{
2819 return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2820}
2821
2822
2823
2824template <int dim, typename Number, typename VectorizedArrayType>
2825inline unsigned int
2827 const unsigned int quad_index,
2828 const unsigned int active_fe_index) const
2829{
2830 AssertIndexRange(quad_index, mapping_info.face_data.size());
2831 return mapping_info.face_data[quad_index]
2832 .descriptor[active_fe_index]
2833 .n_q_points;
2834}
2835
2836
2837
2838template <int dim, typename Number, typename VectorizedArrayType>
2839inline const IndexSet &
2841 const unsigned int dof_handler_index) const
2842{
2843 return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2844}
2845
2846
2847
2848template <int dim, typename Number, typename VectorizedArrayType>
2849inline const IndexSet &
2851 const unsigned int dof_handler_index) const
2852{
2853 return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2854}
2855
2856
2857
2858template <int dim, typename Number, typename VectorizedArrayType>
2861 const unsigned int dof_handler_index,
2862 const unsigned int index_quad,
2863 const unsigned int index_fe,
2864 const unsigned int active_fe_index,
2865 const unsigned int active_quad_index) const
2866{
2867 AssertIndexRange(dof_handler_index, dof_info.size());
2868 const unsigned int ind =
2869 dof_info[dof_handler_index].global_base_element_offset + index_fe;
2870 AssertIndexRange(ind, shape_info.size(0));
2871 AssertIndexRange(index_quad, shape_info.size(1));
2872 AssertIndexRange(active_fe_index, shape_info.size(2));
2873 AssertIndexRange(active_quad_index, shape_info.size(3));
2874 return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2875}
2876
2877
2878
2879template <int dim, typename Number, typename VectorizedArrayType>
2881 VectorizedArrayType::size()> &
2883 const unsigned int face_batch_index) const
2884{
2885 AssertIndexRange(face_batch_index, face_info.faces.size());
2886 return face_info.faces[face_batch_index];
2887}
2888
2889
2890
2891template <int dim, typename Number, typename VectorizedArrayType>
2892inline const Table<3, unsigned int> &
2894 const
2895{
2896 return face_info.cell_and_face_to_plain_faces;
2897}
2898
2899
2900
2901template <int dim, typename Number, typename VectorizedArrayType>
2902inline const Quadrature<dim> &
2904 const unsigned int quad_index,
2905 const unsigned int active_fe_index) const
2906{
2907 AssertIndexRange(quad_index, mapping_info.cell_data.size());
2908 return mapping_info.cell_data[quad_index]
2909 .descriptor[active_fe_index]
2910 .quadrature;
2911}
2912
2913
2914
2915template <int dim, typename Number, typename VectorizedArrayType>
2916inline const Quadrature<dim - 1> &
2918 const unsigned int quad_index,
2919 const unsigned int active_fe_index) const
2920{
2921 AssertIndexRange(quad_index, mapping_info.face_data.size());
2922 return mapping_info.face_data[quad_index]
2923 .descriptor[active_fe_index]
2924 .quadrature;
2925}
2926
2927
2928
2929template <int dim, typename Number, typename VectorizedArrayType>
2930inline unsigned int
2932 const std::pair<unsigned int, unsigned int> range) const
2933{
2934 auto result = get_cell_category(range.first);
2935
2936 for (unsigned int i = range.first; i < range.second; ++i)
2937 result = std::max(result, get_cell_category(i));
2938
2939 return result;
2940}
2941
2942
2943
2944template <int dim, typename Number, typename VectorizedArrayType>
2945inline std::pair<unsigned int, unsigned int>
2947 const std::pair<unsigned int, unsigned int> range) const
2948{
2949 auto result = get_face_category(range.first);
2950
2951 for (unsigned int i = range.first; i < range.second; ++i)
2952 {
2953 result.first = std::max(result.first, get_face_category(i).first);
2954 result.second = std::max(result.second, get_face_category(i).second);
2955 }
2956
2957 return result;
2958}
2959
2960
2961
2962template <int dim, typename Number, typename VectorizedArrayType>
2963inline unsigned int
2965 const unsigned int cell_batch_index) const
2966{
2967 AssertIndexRange(0, dof_info.size());
2968 AssertIndexRange(cell_batch_index, dof_info[0].cell_active_fe_index.size());
2969 if (dof_info[0].cell_active_fe_index.empty())
2970 return 0;
2971 else
2972 return dof_info[0].cell_active_fe_index[cell_batch_index];
2973}
2974
2975
2976
2977template <int dim, typename Number, typename VectorizedArrayType>
2978inline std::pair<unsigned int, unsigned int>
2980 const unsigned int face_batch_index) const
2981{
2982 AssertIndexRange(face_batch_index, face_info.faces.size());
2983 if (dof_info[0].cell_active_fe_index.empty())
2984 return std::make_pair(0U, 0U);
2985
2986 std::pair<unsigned int, unsigned int> result = std::make_pair(0U, 0U);
2987 for (unsigned int v = 0;
2988 v < VectorizedArrayType::size() &&
2989 face_info.faces[face_batch_index].cells_interior[v] !=
2991 ++v)
2992 result.first = std::max(
2993 result.first,
2994 dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
2995 .cells_interior[v] /
2996 VectorizedArrayType::size()]);
2997 if (face_info.faces[face_batch_index].cells_exterior[0] !=
2999 for (unsigned int v = 0;
3000 v < VectorizedArrayType::size() &&
3001 face_info.faces[face_batch_index].cells_exterior[v] !=
3003 ++v)
3004 result.second = std::max(
3005 result.second,
3006 dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
3007 .cells_exterior[v] /
3008 VectorizedArrayType::size()]);
3009 else
3010 result.second = numbers::invalid_unsigned_int;
3011 return result;
3012}
3013
3014
3015
3016template <int dim, typename Number, typename VectorizedArrayType>
3017inline bool
3019{
3020 return indices_are_initialized;
3021}
3022
3023
3024
3025template <int dim, typename Number, typename VectorizedArrayType>
3026inline bool
3028{
3029 return mapping_is_initialized;
3030}
3031
3032
3033template <int dim, typename Number, typename VectorizedArrayType>
3034inline unsigned int
3036{
3037 return mg_level;
3038}
3039
3040
3041
3042template <int dim, typename Number, typename VectorizedArrayType>
3045{
3046 using list_type =
3047 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3048 list_type &data = scratch_pad.get();
3049 for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
3050 if (it->first == false)
3051 {
3052 it->first = true;
3053 return &it->second;
3054 }
3055 data.emplace_front(true, AlignedVector<VectorizedArrayType>());
3056 return &data.front().second;
3057}
3058
3059
3060
3061template <int dim, typename Number, typename VectorizedArrayType>
3062void
3064 const AlignedVector<VectorizedArrayType> *scratch) const
3065{
3066 using list_type =
3067 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3068 list_type &data = scratch_pad.get();
3069 for (typename list_type::iterator it = data.begin(); it != data.end(); ++it)
3070 if (&it->second == scratch)
3071 {
3072 Assert(it->first == true, ExcInternalError());
3073 it->first = false;
3074 return;
3075 }
3076 AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
3077}
3078
3079
3080
3081template <int dim, typename Number, typename VectorizedArrayType>
3085{
3086 for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
3087 scratch_pad_non_threadsafe.begin();
3088 it != scratch_pad_non_threadsafe.end();
3089 ++it)
3090 if (it->first == false)
3091 {
3092 it->first = true;
3093 return &it->second;
3094 }
3095 scratch_pad_non_threadsafe.push_front(
3096 std::make_pair(true, AlignedVector<Number>()));
3097 return &scratch_pad_non_threadsafe.front().second;
3098}
3099
3100
3101
3102template <int dim, typename Number, typename VectorizedArrayType>
3103void
3106 const AlignedVector<Number> *scratch) const
3107{
3108 for (typename std::list<std::pair<bool, AlignedVector<Number>>>::iterator it =
3109 scratch_pad_non_threadsafe.begin();
3110 it != scratch_pad_non_threadsafe.end();
3111 ++it)
3112 if (&it->second == scratch)
3113 {
3114 Assert(it->first == true, ExcInternalError());
3115 it->first = false;
3116 return;
3117 }
3118 AssertThrow(false, ExcMessage("Tried to release invalid scratch pad"));
3119}
3120
3121
3122
3123// ------------------------------ reinit functions ---------------------------
3124
3125namespace internal
3126{
3127 namespace MatrixFreeImplementation
3128 {
3129 template <int dim, int spacedim>
3130 inline std::vector<IndexSet>
3131 extract_locally_owned_index_sets(
3132 const std::vector<const DoFHandler<dim, spacedim> *> &dofh,
3133 const unsigned int level)
3134 {
3135 std::vector<IndexSet> locally_owned_set;
3136 locally_owned_set.reserve(dofh.size());
3137 for (unsigned int j = 0; j < dofh.size(); ++j)
3139 locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
3140 else
3141 locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
3142 return locally_owned_set;
3143 }
3144 } // namespace MatrixFreeImplementation
3145} // namespace internal
3146
3147
3148
3149template <int dim, typename Number, typename VectorizedArrayType>
3150template <typename QuadratureType, typename number2, typename MappingType>
3151void
3153 const MappingType &mapping,
3154 const DoFHandler<dim> &dof_handler,
3155 const AffineConstraints<number2> &constraints_in,
3156 const QuadratureType &quad,
3158 &additional_data)
3159{
3160 std::vector<const DoFHandler<dim, dim> *> dof_handlers;
3161 std::vector<const AffineConstraints<number2> *> constraints;
3162
3163 dof_handlers.push_back(&dof_handler);
3164 constraints.push_back(&constraints_in);
3165
3166 std::vector<IndexSet> locally_owned_sets =
3167 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3168 dof_handlers, additional_data.mg_level);
3169
3170 std::vector<hp::QCollection<dim>> quad_hp;
3171 quad_hp.emplace_back(quad);
3172
3173 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3174 dof_handlers,
3175 constraints,
3176 locally_owned_sets,
3177 quad_hp,
3178 additional_data);
3179}
3180
3181
3182
3183template <int dim, typename Number, typename VectorizedArrayType>
3184template <typename QuadratureType, typename number2, typename MappingType>
3185void
3187 const MappingType &mapping,
3188 const std::vector<const DoFHandler<dim> *> &dof_handler,
3189 const std::vector<const AffineConstraints<number2> *> &constraint,
3190 const QuadratureType &quad,
3192 &additional_data)
3193{
3194 std::vector<IndexSet> locally_owned_set =
3195 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3196 dof_handler, additional_data.mg_level);
3197 std::vector<hp::QCollection<dim>> quad_hp;
3198 quad_hp.emplace_back(quad);
3199
3200 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3201 dof_handler,
3202 constraint,
3203 locally_owned_set,
3204 quad_hp,
3205 additional_data);
3206}
3207
3208
3209
3210template <int dim, typename Number, typename VectorizedArrayType>
3211template <typename QuadratureType, typename number2, typename MappingType>
3212void
3214 const MappingType &mapping,
3215 const std::vector<const DoFHandler<dim> *> &dof_handler,
3216 const std::vector<const AffineConstraints<number2> *> &constraint,
3217 const std::vector<QuadratureType> &quad,
3219 &additional_data)
3220{
3221 std::vector<IndexSet> locally_owned_set =
3222 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3223 dof_handler, additional_data.mg_level);
3224 std::vector<hp::QCollection<dim>> quad_hp;
3225 for (unsigned int q = 0; q < quad.size(); ++q)
3226 quad_hp.emplace_back(quad[q]);
3227
3228 internal_reinit(std::make_shared<hp::MappingCollection<dim>>(mapping),
3229 dof_handler,
3230 constraint,
3231 locally_owned_set,
3232 quad_hp,
3233 additional_data);
3234}
3235
3236
3237
3238// ------------------------------ implementation of loops --------------------
3239
3240// internal helper functions that define how to call MPI data exchange
3241// functions: for generic vectors, do nothing at all. For distributed vectors,
3242// call update_ghost_values_start function and so on. If we have collections
3243// of vectors, just do the individual functions of the components. In order to
3244// keep ghost values consistent (whether we are in read or write mode), we
3245// also reset the values at the end. the whole situation is a bit complicated
3246// by the fact that we need to treat block vectors differently, which use some
3247// additional helper functions to select the blocks and template magic.
3248namespace internal
3249{
3253 template <int dim, typename Number, typename VectorizedArrayType>
3254 struct VectorDataExchange
3255 {
3256 // A shift for the MPI messages to reduce the risk for accidental
3257 // interaction with other open communications that a user program might
3258 // set up (parallel vectors support unfinished communication). We let
3259 // the other vectors use the first 20 assigned numbers and start the
3260 // matrix-free communication.
3261 static constexpr unsigned int channel_shift = 20;
3262
3263
3264
3269 VectorDataExchange(
3270 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3271 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3272 DataAccessOnFaces vector_face_access,
3273 const unsigned int n_components)
3274 : matrix_free(matrix_free)
3275 , vector_face_access(
3276 matrix_free.get_task_info().face_partition_data.empty() ?
3277 ::MatrixFree<dim, Number, VectorizedArrayType>::
3278 DataAccessOnFaces::unspecified :
3279 vector_face_access)
3280 , ghosts_were_set(false)
3281# ifdef DEAL_II_WITH_MPI
3282 , tmp_data(n_components)
3283 , requests(n_components)
3284# endif
3285 {
3286 (void)n_components;
3287 if (this->vector_face_access !=
3289 DataAccessOnFaces::unspecified)
3290 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3292 matrix_free.get_dof_info(c).vector_exchanger_face_variants.size(),
3293 5);
3294 }
3295
3296
3297
3301 ~VectorDataExchange() // NOLINT
3302 {
3303# ifdef DEAL_II_WITH_MPI
3304 for (unsigned int i = 0; i < tmp_data.size(); ++i)
3305 if (tmp_data[i] != nullptr)
3306 matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3307# endif
3308 }
3309
3310
3311
3316 template <typename VectorType>
3317 unsigned int
3318 find_vector_in_mf(const VectorType &vec,
3319 const bool check_global_compatibility = true) const
3320 {
3321 // case 1: vector was set up with MatrixFree::initialize_dof_vector()
3322 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3323 if (vec.get_partitioner().get() ==
3324 matrix_free.get_dof_info(c).vector_partitioner.get())
3325 return c;
3326
3327 // case 2: user provided own partitioner (compatibility mode)
3328 for (unsigned int c = 0; c < matrix_free.n_components(); ++c)
3329 if (check_global_compatibility ?
3330 vec.get_partitioner()->is_globally_compatible(
3331 *matrix_free.get_dof_info(c).vector_partitioner) :
3332 vec.get_partitioner()->is_compatible(
3333 *matrix_free.get_dof_info(c).vector_partitioner))
3334 return c;
3335
3336 Assert(false,
3337 ExcNotImplemented("Could not find partitioner that fits vector"));
3338
3340 }
3341
3342
3343
3349 get_partitioner(const unsigned int mf_component) const
3350 {
3351 AssertDimension(matrix_free.get_dof_info(mf_component)
3352 .vector_exchanger_face_variants.size(),
3353 5);
3354 if (vector_face_access ==
3356 DataAccessOnFaces::none)
3357 return *matrix_free.get_dof_info(mf_component)
3358 .vector_exchanger_face_variants[0];
3359 else if (vector_face_access ==
3361 DataAccessOnFaces::values)
3362 return *matrix_free.get_dof_info(mf_component)
3363 .vector_exchanger_face_variants[1];
3364 else if (vector_face_access ==
3366 DataAccessOnFaces::gradients)
3367 return *matrix_free.get_dof_info(mf_component)
3368 .vector_exchanger_face_variants[2];
3369 else if (vector_face_access ==
3371 DataAccessOnFaces::values_all_faces)
3372 return *matrix_free.get_dof_info(mf_component)
3373 .vector_exchanger_face_variants[3];
3374 else if (vector_face_access ==
3376 DataAccessOnFaces::gradients_all_faces)
3377 return *matrix_free.get_dof_info(mf_component)
3378 .vector_exchanger_face_variants[4];
3379 else
3380 return *matrix_free.get_dof_info(mf_component).vector_exchanger.get();
3381 }
3382
3383
3384
3388 template <typename VectorType,
3389 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3390 * = nullptr>
3391 void
3392 update_ghost_values_start(const unsigned int /*component_in_block_vector*/,
3393 const VectorType & /*vec*/)
3394 {}
3395
3396
3401 template <typename VectorType,
3402 std::enable_if_t<!has_update_ghost_values_start<VectorType> &&
3403 !is_not_parallel_vector<VectorType>,
3404 VectorType> * = nullptr>
3405 void
3406 update_ghost_values_start(const unsigned int component_in_block_vector,
3407 const VectorType &vec)
3408 {
3409 (void)component_in_block_vector;
3410 const bool ghosts_set = vec.has_ghost_elements();
3411
3412 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3413 ghosts_set == false,
3415
3416 if (ghosts_set)
3417 ghosts_were_set = true;
3418
3419 vec.update_ghost_values();
3420 }
3421
3422
3423
3429 template <typename VectorType,
3430 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3431 !has_exchange_on_subset<VectorType>,
3432 VectorType> * = nullptr>
3433 void
3434 update_ghost_values_start(const unsigned int component_in_block_vector,
3435 const VectorType &vec)
3436 {
3437 (void)component_in_block_vector;
3438 const bool ghosts_set = vec.has_ghost_elements();
3439
3440 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3441 ghosts_set == false,
3443
3444 if (ghosts_set)
3445 ghosts_were_set = true;
3446
3447 vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3448 }
3449
3450
3451
3458 template <typename VectorType,
3459 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3460 has_exchange_on_subset<VectorType>,
3461 VectorType> * = nullptr>
3462 void
3463 update_ghost_values_start(const unsigned int component_in_block_vector,
3464 const VectorType &vec)
3465 {
3466 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3467 "Type mismatch between VectorType and VectorDataExchange");
3468 (void)component_in_block_vector;
3469 const bool ghosts_set = vec.has_ghost_elements();
3470
3471 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3472 ghosts_set == false,
3474
3475 if (ghosts_set)
3476 ghosts_were_set = true;
3477
3478 if (vec.size() != 0)
3479 {
3480# ifdef DEAL_II_WITH_MPI
3481 const unsigned int mf_component = find_vector_in_mf(vec);
3482
3483 const auto &part = get_partitioner(mf_component);
3484
3485 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3486 part.n_import_sm_procs() == 0)
3487 return;
3488
3489 tmp_data[component_in_block_vector] =
3490 matrix_free.acquire_scratch_data_non_threadsafe();
3491 tmp_data[component_in_block_vector]->resize_fast(
3492 part.n_import_indices());
3493 AssertDimension(requests.size(), tmp_data.size());
3494
3495 part.export_to_ghosted_array_start(
3496 component_in_block_vector * 2 + channel_shift,
3497 ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3498 vec.shared_vector_data(),
3499 ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3500 part.locally_owned_size(),
3501 matrix_free.get_dof_info(mf_component)
3502 .vector_partitioner->n_ghost_indices()),
3503 ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3504 part.n_import_indices()),
3505 this->requests[component_in_block_vector]);
3506# endif
3507 }
3508 }
3509
3510
3511
3516 template <typename VectorType,
3517 std::enable_if_t<!has_update_ghost_values_start<VectorType>,
3518 VectorType> * = nullptr>
3519 void
3520 update_ghost_values_finish(const unsigned int /*component_in_block_vector*/,
3521 const VectorType & /*vec*/)
3522 {}
3523
3524
3525
3531 template <typename VectorType,
3532 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3533 !has_exchange_on_subset<VectorType>,
3534 VectorType> * = nullptr>
3535 void
3536 update_ghost_values_finish(const unsigned int component_in_block_vector,
3537 const VectorType &vec)
3538 {
3539 (void)component_in_block_vector;
3540 vec.update_ghost_values_finish();
3541 }
3542
3543
3544
3551 template <typename VectorType,
3552 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3553 has_exchange_on_subset<VectorType>,
3554 VectorType> * = nullptr>
3555 void
3556 update_ghost_values_finish(const unsigned int component_in_block_vector,
3557 const VectorType &vec)
3558 {
3559 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3560 "Type mismatch between VectorType and VectorDataExchange");
3561 (void)component_in_block_vector;
3562
3563 if (vec.size() != 0)
3564 {
3565# ifdef DEAL_II_WITH_MPI
3566 AssertIndexRange(component_in_block_vector, tmp_data.size());
3567 AssertDimension(requests.size(), tmp_data.size());
3568
3569 const unsigned int mf_component = find_vector_in_mf(vec);
3570
3571 const auto &part = get_partitioner(mf_component);
3572
3573 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3574 part.n_import_sm_procs() != 0)
3575 {
3576 part.export_to_ghosted_array_finish(
3577 ArrayView<const Number>(vec.begin(), part.locally_owned_size()),
3578 vec.shared_vector_data(),
3579 ArrayView<Number>(const_cast<Number *>(vec.begin()) +
3580 part.locally_owned_size(),
3581 matrix_free.get_dof_info(mf_component)
3582 .vector_partitioner->n_ghost_indices()),
3583 this->requests[component_in_block_vector]);
3584
3585 matrix_free.release_scratch_data_non_threadsafe(
3586 tmp_data[component_in_block_vector]);
3587 tmp_data[component_in_block_vector] = nullptr;
3588 }
3589# endif
3590 }
3591 // let vector know that ghosts are being updated and we can read from
3592 // them
3593 vec.set_ghost_state(true);
3594 }
3595
3596
3597
3601 template <typename VectorType,
3602 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3603 * = nullptr>
3604 void
3605 compress_start(const unsigned int /*component_in_block_vector*/,
3606 VectorType & /*vec*/)
3607 {}
3608
3609
3610
3615 template <typename VectorType,
3616 std::enable_if_t<!has_compress_start<VectorType> &&
3617 !is_not_parallel_vector<VectorType>,
3618 VectorType> * = nullptr>
3619 void
3620 compress_start(const unsigned int component_in_block_vector,
3621 VectorType &vec)
3622 {
3623 (void)component_in_block_vector;
3624 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3625 vec.compress(VectorOperation::add);
3626 }
3627
3628
3629
3635 template <typename VectorType,
3636 std::enable_if_t<has_compress_start<VectorType> &&
3637 !has_exchange_on_subset<VectorType>,
3638 VectorType> * = nullptr>
3639 void
3640 compress_start(const unsigned int component_in_block_vector,
3641 VectorType &vec)
3642 {
3643 (void)component_in_block_vector;
3644 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3645 vec.compress_start(component_in_block_vector + channel_shift);
3646 }
3647
3648
3649
3656 template <typename VectorType,
3657 std::enable_if_t<has_compress_start<VectorType> &&
3658 has_exchange_on_subset<VectorType>,
3659 VectorType> * = nullptr>
3660 void
3661 compress_start(const unsigned int component_in_block_vector,
3662 VectorType &vec)
3663 {
3664 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3665 "Type mismatch between VectorType and VectorDataExchange");
3666 (void)component_in_block_vector;
3667 Assert(vec.has_ghost_elements() == false, ExcNotImplemented());
3668
3669 if (vec.size() != 0)
3670 {
3671# ifdef DEAL_II_WITH_MPI
3672 const unsigned int mf_component = find_vector_in_mf(vec);
3673
3674 const auto &part = get_partitioner(mf_component);
3675
3676 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3677 part.n_import_sm_procs() == 0)
3678 return;
3679
3680 tmp_data[component_in_block_vector] =
3681 matrix_free.acquire_scratch_data_non_threadsafe();
3682 tmp_data[component_in_block_vector]->resize_fast(
3683 part.n_import_indices());
3684 AssertDimension(requests.size(), tmp_data.size());
3685
3686 part.import_from_ghosted_array_start(
3688 component_in_block_vector * 2 + channel_shift,
3689 ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3690 vec.shared_vector_data(),
3691 ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3692 matrix_free.get_dof_info(mf_component)
3693 .vector_partitioner->n_ghost_indices()),
3694 ArrayView<Number>(tmp_data[component_in_block_vector]->begin(),
3695 part.n_import_indices()),
3696 this->requests[component_in_block_vector]);
3697# endif
3698 }
3699 }
3700
3701
3702
3707 template <
3708 typename VectorType,
3709 std::enable_if_t<!has_compress_start<VectorType>, VectorType> * = nullptr>
3710 void
3711 compress_finish(const unsigned int /*component_in_block_vector*/,
3712 VectorType & /*vec*/)
3713 {}
3714
3715
3716
3722 template <typename VectorType,
3723 std::enable_if_t<has_compress_start<VectorType> &&
3724 !has_exchange_on_subset<VectorType>,
3725 VectorType> * = nullptr>
3726 void
3727 compress_finish(const unsigned int component_in_block_vector,
3728 VectorType &vec)
3729 {
3730 (void)component_in_block_vector;
3731 vec.compress_finish(VectorOperation::add);
3732 }
3733
3734
3735
3742 template <typename VectorType,
3743 std::enable_if_t<has_compress_start<VectorType> &&
3744 has_exchange_on_subset<VectorType>,
3745 VectorType> * = nullptr>
3746 void
3747 compress_finish(const unsigned int component_in_block_vector,
3748 VectorType &vec)
3749 {
3750 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3751 "Type mismatch between VectorType and VectorDataExchange");
3752 (void)component_in_block_vector;
3753 if (vec.size() != 0)
3754 {
3755# ifdef DEAL_II_WITH_MPI
3756 AssertIndexRange(component_in_block_vector, tmp_data.size());
3757 AssertDimension(requests.size(), tmp_data.size());
3758
3759 const unsigned int mf_component = find_vector_in_mf(vec);
3760
3761 const auto &part = get_partitioner(mf_component);
3762
3763 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3764 part.n_import_sm_procs() != 0)
3765 {
3766 part.import_from_ghosted_array_finish(
3768 ArrayView<Number>(vec.begin(), part.locally_owned_size()),
3769 vec.shared_vector_data(),
3770 ArrayView<Number>(vec.begin() + part.locally_owned_size(),
3771 matrix_free.get_dof_info(mf_component)
3772 .vector_partitioner->n_ghost_indices()),
3774 tmp_data[component_in_block_vector]->begin(),
3775 part.n_import_indices()),
3776 this->requests[component_in_block_vector]);
3777
3778 matrix_free.release_scratch_data_non_threadsafe(
3779 tmp_data[component_in_block_vector]);
3780 tmp_data[component_in_block_vector] = nullptr;
3781 }
3782
3784 {
3785 const int ierr =
3786 MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3787 AssertThrowMPI(ierr);
3788 }
3789# endif
3790 }
3791 }
3792
3793
3794
3798 template <typename VectorType,
3799 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3800 * = nullptr>
3801 void
3802 reset_ghost_values(const VectorType & /*vec*/) const
3803 {}
3804
3805
3806
3811 template <typename VectorType,
3812 std::enable_if_t<!has_exchange_on_subset<VectorType> &&
3813 !is_not_parallel_vector<VectorType>,
3814 VectorType> * = nullptr>
3815 void
3816 reset_ghost_values(const VectorType &vec) const
3817 {
3818 if (ghosts_were_set == true)
3819 return;
3820
3821 vec.zero_out_ghost_values();
3822 }
3823
3824
3825
3831 template <typename VectorType,
3832 std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3833 * = nullptr>
3834 void
3835 reset_ghost_values(const VectorType &vec) const
3836 {
3837 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3838 "Type mismatch between VectorType and VectorDataExchange");
3839 if (ghosts_were_set == true)
3840 return;
3841
3842 if (vec.size() != 0)
3843 {
3844# ifdef DEAL_II_WITH_MPI
3845 AssertDimension(requests.size(), tmp_data.size());
3846
3847 const unsigned int mf_component = find_vector_in_mf(vec);
3848
3849 const auto &part = get_partitioner(mf_component);
3850
3851 if (part.n_ghost_indices() > 0)
3852 {
3853 part.reset_ghost_values(ArrayView<Number>(
3855 .begin() +
3856 part.locally_owned_size(),
3857 matrix_free.get_dof_info(mf_component)
3858 .vector_partitioner->n_ghost_indices()));
3859 }
3860
3861# endif
3862 }
3863 // let vector know that it's not ghosted anymore
3864 vec.set_ghost_state(false);
3865 }
3866
3867
3868
3874 template <typename VectorType,
3875 std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3876 * = nullptr>
3877 void
3878 zero_vector_region(const unsigned int range_index, VectorType &vec) const
3879 {
3880 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3881 "Type mismatch between VectorType and VectorDataExchange");
3882 if (range_index == numbers::invalid_unsigned_int)
3883 vec = Number();
3884 else
3885 {
3886 const unsigned int mf_component = find_vector_in_mf(vec, false);
3888 matrix_free.get_dof_info(mf_component);
3889 Assert(dof_info.vector_zero_range_list_index.empty() == false,
3891
3892 Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3894 AssertIndexRange(range_index,
3895 dof_info.vector_zero_range_list_index.size() - 1);
3896 for (unsigned int id =
3897 dof_info.vector_zero_range_list_index[range_index];
3898 id != dof_info.vector_zero_range_list_index[range_index + 1];
3899 ++id)
3900 std::memset(vec.begin() + dof_info.vector_zero_range_list[id].first,
3901 0,
3902 (dof_info.vector_zero_range_list[id].second -
3903 dof_info.vector_zero_range_list[id].first) *
3904 sizeof(Number));
3905 }
3906 }
3907
3908
3909
3915 template <typename VectorType,
3916 std::enable_if_t<!has_exchange_on_subset<VectorType>, VectorType>
3917 * = nullptr,
3918 typename VectorType::value_type * = nullptr>
3919 void
3920 zero_vector_region(const unsigned int range_index, VectorType &vec) const
3921 {
3922 if (range_index == numbers::invalid_unsigned_int || range_index == 0)
3923 vec = typename VectorType::value_type();
3924 }
3925
3926
3927
3932 void
3933 zero_vector_region(const unsigned int, ...) const
3934 {
3935 Assert(false,
3936 ExcNotImplemented("Zeroing is only implemented for vector types "
3937 "which provide VectorType::value_type"));
3938 }
3939
3940
3941
3942 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
3943 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3944 DataAccessOnFaces vector_face_access;
3945 bool ghosts_were_set;
3946# ifdef DEAL_II_WITH_MPI
3947 std::vector<AlignedVector<Number> *> tmp_data;
3948 std::vector<std::vector<MPI_Request>> requests;
3949# endif
3950 }; // VectorDataExchange
3951
3952 template <typename VectorStruct>
3953 unsigned int
3954 n_components(const VectorStruct &vec);
3955
3956 template <typename VectorStruct>
3957 unsigned int
3958 n_components_block(const VectorStruct &vec, const std::bool_constant<true>)
3959 {
3960 unsigned int components = 0;
3961 for (unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
3962 components += n_components(vec.block(bl));
3963 return components;
3964 }
3965
3966 template <typename VectorStruct>
3967 unsigned int
3968 n_components_block(const VectorStruct &, const std::bool_constant<false>)
3969 {
3970 return 1;
3971 }
3972
3973 template <typename VectorStruct>
3974 unsigned int
3975 n_components(const VectorStruct &vec)
3976 {
3977 return n_components_block(
3978 vec, std::bool_constant<IsBlockVector<VectorStruct>::value>());
3979 }
3980
3981 template <typename VectorStruct>
3982 inline unsigned int
3983 n_components(const std::vector<VectorStruct> &vec)
3984 {
3985 unsigned int components = 0;
3986 for (unsigned int comp = 0; comp < vec.size(); ++comp)
3987 components += n_components_block(
3988 vec[comp], std::bool_constant<IsBlockVector<VectorStruct>::value>());
3989 return components;
3990 }
3991
3992 template <typename VectorStruct>
3993 inline unsigned int
3994 n_components(const std::vector<VectorStruct *> &vec)
3995 {
3996 unsigned int components = 0;
3997 for (unsigned int comp = 0; comp < vec.size(); ++comp)
3998 components += n_components_block(
3999 *vec[comp], std::bool_constant<IsBlockVector<VectorStruct>::value>());
4000 return components;
4001 }
4002
4003
4004
4005 // A helper function to identify block vectors with many components where we
4006 // should not try to overlap computations and communication because there
4007 // would be too many outstanding communication requests.
4008
4009 // default value for vectors that do not have communication_block_size
4010 template <typename VectorStruct,
4011 std::enable_if_t<!has_communication_block_size<VectorStruct>,
4012 VectorStruct> * = nullptr>
4013 constexpr unsigned int
4014 get_communication_block_size(const VectorStruct &)
4015 {
4017 }
4018
4019
4020
4021 template <typename VectorStruct,
4022 std::enable_if_t<has_communication_block_size<VectorStruct>,
4023 VectorStruct> * = nullptr>
4024 constexpr unsigned int
4025 get_communication_block_size(const VectorStruct &)
4026 {
4027 return VectorStruct::communication_block_size;
4028 }
4029
4030
4031
4032 template <typename VectorType,
4033 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType> * =
4034 nullptr>
4035 bool
4036 has_ghost_elements(const VectorType &vec)
4037 {
4038 (void)vec;
4039 return false;
4040 }
4041
4042
4043
4044 template <typename VectorType,
4045 std::enable_if_t<!is_not_parallel_vector<VectorType>, VectorType>
4046 * = nullptr>
4047 bool
4048 has_ghost_elements(const VectorType &vec)
4049 {
4050 return vec.has_ghost_elements();
4051 }
4052
4053
4054
4055 // --------------------------------------------------------------------------
4056 // below we have wrappers to distinguish between block and non-block vectors.
4057 // --------------------------------------------------------------------------
4058
4059 //
4060 // update_ghost_values_start
4061 //
4062
4063 // update_ghost_values for block vectors
4064 template <int dim,
4065 typename VectorStruct,
4066 typename Number,
4067 typename VectorizedArrayType,
4068 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4069 * = nullptr>
4070 void
4071 update_ghost_values_start(
4072 const VectorStruct &vec,
4073 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4074 const unsigned int channel = 0)
4075 {
4076 if (get_communication_block_size(vec) < vec.n_blocks())
4077 {
4078 const bool ghosts_set = vec.has_ghost_elements();
4079
4080 Assert(exchanger.matrix_free.get_task_info()
4081 .allow_ghosted_vectors_in_loops ||
4082 ghosts_set == false,
4084
4085 if (ghosts_set)
4086 exchanger.ghosts_were_set = true;
4087
4088 vec.update_ghost_values();
4089 }
4090 else
4091 {
4092 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4093 update_ghost_values_start(vec.block(i), exchanger, channel + i);
4094 }
4095 }
4096
4097
4098
4099 // update_ghost_values for non-block vectors
4100 template <int dim,
4101 typename VectorStruct,
4102 typename Number,
4103 typename VectorizedArrayType,
4104 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4105 * = nullptr>
4106 void
4107 update_ghost_values_start(
4108 const VectorStruct &vec,
4109 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4110 const unsigned int channel = 0)
4111 {
4112 exchanger.update_ghost_values_start(channel, vec);
4113 }
4114
4115
4116
4117 // update_ghost_values_start() for vector of vectors
4118 template <int dim,
4119 typename VectorStruct,
4120 typename Number,
4121 typename VectorizedArrayType>
4122 inline void
4123 update_ghost_values_start(
4124 const std::vector<VectorStruct> &vec,
4125 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4126 {
4127 unsigned int component_index = 0;
4128 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4129 {
4130 update_ghost_values_start(vec[comp], exchanger, component_index);
4131 component_index += n_components(vec[comp]);
4132 }
4133 }
4134
4135
4136
4137 // update_ghost_values_start() for vector of pointers to vectors
4138 template <int dim,
4139 typename VectorStruct,
4140 typename Number,
4141 typename VectorizedArrayType>
4142 inline void
4143 update_ghost_values_start(
4144 const std::vector<VectorStruct *> &vec,
4145 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4146 {
4147 unsigned int component_index = 0;
4148 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4149 {
4150 update_ghost_values_start(*vec[comp], exchanger, component_index);
4151 component_index += n_components(*vec[comp]);
4152 }
4153 }
4154
4155
4156
4157 //
4158 // update_ghost_values_finish
4159 //
4160
4161 // for block vectors
4162 template <int dim,
4163 typename VectorStruct,
4164 typename Number,
4165 typename VectorizedArrayType,
4166 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4167 * = nullptr>
4168 void
4169 update_ghost_values_finish(
4170 const VectorStruct &vec,
4171 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4172 const unsigned int channel = 0)
4173 {
4174 if (get_communication_block_size(vec) < vec.n_blocks())
4175 {
4176 // do nothing, everything has already been completed in the _start()
4177 // call
4178 }
4179 else
4180 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4181 update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4182 }
4183
4184
4185
4186 // for non-block vectors
4187 template <int dim,
4188 typename VectorStruct,
4189 typename Number,
4190 typename VectorizedArrayType,
4191 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4192 * = nullptr>
4193 void
4194 update_ghost_values_finish(
4195 const VectorStruct &vec,
4196 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4197 const unsigned int channel = 0)
4198 {
4199 exchanger.update_ghost_values_finish(channel, vec);
4200 }
4201
4202
4203
4204 // for vector of vectors
4205 template <int dim,
4206 typename VectorStruct,
4207 typename Number,
4208 typename VectorizedArrayType>
4209 inline void
4210 update_ghost_values_finish(
4211 const std::vector<VectorStruct> &vec,
4212 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4213 {
4214 unsigned int component_index = 0;
4215 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4216 {
4217 update_ghost_values_finish(vec[comp], exchanger, component_index);
4218 component_index += n_components(vec[comp]);
4219 }
4220 }
4221
4222
4223
4224 // for vector of pointers to vectors
4225 template <int dim,
4226 typename VectorStruct,
4227 typename Number,
4228 typename VectorizedArrayType>
4229 inline void
4230 update_ghost_values_finish(
4231 const std::vector<VectorStruct *> &vec,
4232 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4233 {
4234 unsigned int component_index = 0;
4235 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4236 {
4237 update_ghost_values_finish(*vec[comp], exchanger, component_index);
4238 component_index += n_components(*vec[comp]);
4239 }
4240 }
4241
4242
4243
4244 //
4245 // compress_start
4246 //
4247
4248 // for block vectors
4249 template <int dim,
4250 typename VectorStruct,
4251 typename Number,
4252 typename VectorizedArrayType,
4253 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4254 * = nullptr>
4255 inline void
4256 compress_start(
4257 VectorStruct &vec,
4258 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4259 const unsigned int channel = 0)
4260 {
4261 if (get_communication_block_size(vec) < vec.n_blocks())
4262 vec.compress(VectorOperation::add);
4263 else
4264 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4265 compress_start(vec.block(i), exchanger, channel + i);
4266 }
4267
4268
4269
4270 // for non-block vectors
4271 template <int dim,
4272 typename VectorStruct,
4273 typename Number,
4274 typename VectorizedArrayType,
4275 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4276 * = nullptr>
4277 inline void
4278 compress_start(
4279 VectorStruct &vec,
4280 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4281 const unsigned int channel = 0)
4282 {
4283 exchanger.compress_start(channel, vec);
4284 }
4285
4286
4287
4288 // for std::vector of vectors
4289 template <int dim,
4290 typename VectorStruct,
4291 typename Number,
4292 typename VectorizedArrayType>
4293 inline void
4294 compress_start(
4295 std::vector<VectorStruct> &vec,
4296 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4297 {
4298 unsigned int component_index = 0;
4299 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4300 {
4301 compress_start(vec[comp], exchanger, component_index);
4302 component_index += n_components(vec[comp]);
4303 }
4304 }
4305
4306
4307
4308 // for std::vector of pointer to vectors
4309 template <int dim,
4310 typename VectorStruct,
4311 typename Number,
4312 typename VectorizedArrayType>
4313 inline void
4314 compress_start(
4315 std::vector<VectorStruct *> &vec,
4316 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4317 {
4318 unsigned int component_index = 0;
4319 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4320 {
4321 compress_start(*vec[comp], exchanger, component_index);
4322 component_index += n_components(*vec[comp]);
4323 }
4324 }
4325
4326
4327
4328 //
4329 // compress_finish
4330 //
4331
4332 // for block vectors
4333 template <int dim,
4334 typename VectorStruct,
4335 typename Number,
4336 typename VectorizedArrayType,
4337 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4338 * = nullptr>
4339 inline void
4340 compress_finish(
4341 VectorStruct &vec,
4342 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4343 const unsigned int channel = 0)
4344 {
4345 if (get_communication_block_size(vec) < vec.n_blocks())
4346 {
4347 // do nothing, everything has already been completed in the _start()
4348 // call
4349 }
4350 else
4351 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4352 compress_finish(vec.block(i), exchanger, channel + i);
4353 }
4354
4355
4356
4357 // for non-block vectors
4358 template <int dim,
4359 typename VectorStruct,
4360 typename Number,
4361 typename VectorizedArrayType,
4362 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4363 * = nullptr>
4364 inline void
4365 compress_finish(
4366 VectorStruct &vec,
4367 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4368 const unsigned int channel = 0)
4369 {
4370 exchanger.compress_finish(channel, vec);
4371 }
4372
4373
4374
4375 // for std::vector of vectors
4376 template <int dim,
4377 typename VectorStruct,
4378 typename Number,
4379 typename VectorizedArrayType>
4380 inline void
4381 compress_finish(
4382 std::vector<VectorStruct> &vec,
4383 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4384 {
4385 unsigned int component_index = 0;
4386 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4387 {
4388 compress_finish(vec[comp], exchanger, component_index);
4389 component_index += n_components(vec[comp]);
4390 }
4391 }
4392
4393
4394
4395 // for std::vector of pointer to vectors
4396 template <int dim,
4397 typename VectorStruct,
4398 typename Number,
4399 typename VectorizedArrayType>
4400 inline void
4401 compress_finish(
4402 std::vector<VectorStruct *> &vec,
4403 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4404 {
4405 unsigned int component_index = 0;
4406 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4407 {
4408 compress_finish(*vec[comp], exchanger, component_index);
4409 component_index += n_components(*vec[comp]);
4410 }
4411 }
4412
4413
4414
4415 //
4416 // reset_ghost_values:
4417 //
4418 // if the input vector did not have ghosts imported, clear them here again
4419 // in order to avoid subsequent operations e.g. in linear solvers to work
4420 // with ghosts all the time
4421 //
4422
4423 // for block vectors
4424 template <int dim,
4425 typename VectorStruct,
4426 typename Number,
4427 typename VectorizedArrayType,
4428 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4429 * = nullptr>
4430 inline void
4431 reset_ghost_values(
4432 const VectorStruct &vec,
4433 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4434 {
4435 // return immediately if there is nothing to do.
4436 if (exchanger.ghosts_were_set == true)
4437 return;
4438
4439 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4440 reset_ghost_values(vec.block(i), exchanger);
4441 }
4442
4443
4444
4445 // for non-block vectors
4446 template <int dim,
4447 typename VectorStruct,
4448 typename Number,
4449 typename VectorizedArrayType,
4450 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4451 * = nullptr>
4452 inline void
4453 reset_ghost_values(
4454 const VectorStruct &vec,
4455 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4456 {
4457 exchanger.reset_ghost_values(vec);
4458 }
4459
4460
4461
4462 // for std::vector of vectors
4463 template <int dim,
4464 typename VectorStruct,
4465 typename Number,
4466 typename VectorizedArrayType>
4467 inline void
4468 reset_ghost_values(
4469 const std::vector<VectorStruct> &vec,
4470 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4471 {
4472 // return immediately if there is nothing to do.
4473 if (exchanger.ghosts_were_set == true)
4474 return;
4475
4476 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4477 reset_ghost_values(vec[comp], exchanger);
4478 }
4479
4480
4481
4482 // for std::vector of pointer to vectors
4483 template <int dim,
4484 typename VectorStruct,
4485 typename Number,
4486 typename VectorizedArrayType>
4487 inline void
4488 reset_ghost_values(
4489 const std::vector<VectorStruct *> &vec,
4490 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4491 {
4492 // return immediately if there is nothing to do.
4493 if (exchanger.ghosts_were_set == true)
4494 return;
4495
4496 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4497 reset_ghost_values(*vec[comp], exchanger);
4498 }
4499
4500
4501
4502 //
4503 // zero_vector_region
4504 //
4505
4506 // for block vectors
4507 template <int dim,
4508 typename VectorStruct,
4509 typename Number,
4510 typename VectorizedArrayType,
4511 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4512 * = nullptr>
4513 inline void
4514 zero_vector_region(
4515 const unsigned int range_index,
4516 VectorStruct &vec,
4517 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4518 {
4519 for (unsigned int i = 0; i < vec.n_blocks(); ++i)
4520 exchanger.zero_vector_region(range_index, vec.block(i));
4521 }
4522
4523
4524
4525 // for non-block vectors
4526 template <int dim,
4527 typename VectorStruct,
4528 typename Number,
4529 typename VectorizedArrayType,
4530 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4531 * = nullptr>
4532 inline void
4533 zero_vector_region(
4534 const unsigned int range_index,
4535 VectorStruct &vec,
4536 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4537 {
4538 exchanger.zero_vector_region(range_index, vec);
4539 }
4540
4541
4542
4543 // for std::vector of vectors
4544 template <int dim,
4545 typename VectorStruct,
4546 typename Number,
4547 typename VectorizedArrayType>
4548 inline void
4549 zero_vector_region(
4550 const unsigned int range_index,
4551 std::vector<VectorStruct> &vec,
4552 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4553 {
4554 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4555 zero_vector_region(range_index, vec[comp], exchanger);
4556 }
4557
4558
4559
4560 // for std::vector of pointers to vectors
4561 template <int dim,
4562 typename VectorStruct,
4563 typename Number,
4564 typename VectorizedArrayType>
4565 inline void
4566 zero_vector_region(
4567 const unsigned int range_index,
4568 std::vector<VectorStruct *> &vec,
4569 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4570 {
4571 for (unsigned int comp = 0; comp < vec.size(); ++comp)
4572 zero_vector_region(range_index, *vec[comp], exchanger);
4573 }
4574
4575
4576
4577 // Apply a unit matrix operation to constrained DoFs: Default cases where we
4578 // cannot detect a LinearAlgebra::distributed::Vector, we do not do
4579 // anything, else we apply the constraints as a unit operation
4580 template <typename VectorStruct1, typename VectorStruct2>
4581 inline void
4582 apply_operation_to_constrained_dofs(const std::vector<unsigned int> &,
4583 const VectorStruct1 &,
4584 VectorStruct2 &)
4585 {}
4586
4587 template <typename Number>
4588 inline void
4589 apply_operation_to_constrained_dofs(
4590 const std::vector<unsigned int> &constrained_dofs,
4593 {
4594 for (const unsigned int i : constrained_dofs)
4595 dst.local_element(i) = src.local_element(i);
4596 }
4597
4598
4599 namespace MatrixFreeFunctions
4600 {
4601 // struct to select between a const interface and a non-const interface
4602 // for MFWorker
4603 template <typename, typename, typename, typename, bool>
4604 struct InterfaceSelector
4605 {};
4606
4607 // Version for constant functions
4608 template <typename MF,
4609 typename InVector,
4610 typename OutVector,
4611 typename Container>
4612 struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4613 {
4614 using function_type = void (Container::*)(
4615 const MF &,
4616 OutVector &,
4617 const InVector &,
4618 const std::pair<unsigned int, unsigned int> &) const;
4619 };
4620
4621 // Version for non-constant functions
4622 template <typename MF,
4623 typename InVector,
4624 typename OutVector,
4625 typename Container>
4626 struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4627 {
4628 using function_type =
4629 void (Container::*)(const MF &,
4630 OutVector &,
4631 const InVector &,
4632 const std::pair<unsigned int, unsigned int> &);
4633 };
4634 } // namespace MatrixFreeFunctions
4635
4636
4637
4638 // A implementation class for the worker object that runs the various
4639 // operations we want to perform during the matrix-free loop
4640 template <typename MF,
4641 typename InVector,
4642 typename OutVector,
4643 typename Container,
4644 bool is_constant>
4645 class MFWorker : public MFWorkerInterface
4646 {
4647 public:
4648 // An alias to make the arguments further down more readable
4649 using function_type = typename MatrixFreeFunctions::
4650 InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4651 function_type;
4652
4653 // constructor, binds all the arguments to this class
4654 MFWorker(const MF &matrix_free,
4655 const InVector &src,
4656 OutVector &dst,
4657 const bool zero_dst_vector_setting,
4658 const Container &container,
4659 function_type cell_function,
4660 function_type face_function,
4661 function_type boundary_function,
4662 const typename MF::DataAccessOnFaces src_vector_face_access =
4663 MF::DataAccessOnFaces::none,
4664 const typename MF::DataAccessOnFaces dst_vector_face_access =
4665 MF::DataAccessOnFaces::none,
4666 const std::function<void(const unsigned int, const unsigned int)>
4667 &operation_before_loop = {},
4668 const std::function<void(const unsigned int, const unsigned int)>
4669 &operation_after_loop = {},
4670 const unsigned int dof_handler_index_pre_post = 0)
4671 : matrix_free(matrix_free)
4672 , container(const_cast<Container &>(container))
4673 , cell_function(cell_function)
4674 , face_function(face_function)
4675 , boundary_function(boundary_function)
4676 , src(src)
4677 , dst(dst)
4678 , src_data_exchanger(matrix_free,
4679 src_vector_face_access,
4680 n_components(src))
4681 , dst_data_exchanger(matrix_free,
4682 dst_vector_face_access,
4683 n_components(dst))
4684 , src_and_dst_are_same(PointerComparison::equal(&src, &dst))
4685 , zero_dst_vector_setting(zero_dst_vector_setting &&
4686 !src_and_dst_are_same)
4687 , operation_before_loop(operation_before_loop)
4688 , operation_after_loop(operation_after_loop)
4689 , dof_handler_index_pre_post(dof_handler_index_pre_post)
4690 {
4691 Assert(!has_ghost_elements(dst),
4692 ExcMessage("The destination vector passed to the matrix-free "
4693 "loop is ghosted. This is not allowed."));
4694 }
4695
4696 // Runs the cell work. If no function is given, nothing is done
4697 virtual void
4698 cell(const std::pair<unsigned int, unsigned int> &cell_range) override
4699 {
4700 if (cell_function != nullptr && cell_range.second > cell_range.first)
4701 for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4702 {
4703 const auto cell_subrange =
4704 matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4705
4706 if (cell_subrange.second <= cell_subrange.first)
4707 continue;
4708
4709 (container.*
4710 cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4711 }
4712 }
4713
4714 virtual void
4715 cell(const unsigned int range_index) override
4716 {
4717 process_range(cell_function,
4718 matrix_free.get_task_info().cell_partition_data_hp_ptr,
4719 matrix_free.get_task_info().cell_partition_data_hp,
4720 range_index);
4721 }
4722
4723 virtual void
4724 face(const unsigned int range_index) override
4725 {
4726 process_range(face_function,
4727 matrix_free.get_task_info().face_partition_data_hp_ptr,
4728 matrix_free.get_task_info().face_partition_data_hp,
4729 range_index);
4730 }
4731
4732 virtual void
4733 boundary(const unsigned int range_index) override
4734 {
4735 process_range(boundary_function,
4736 matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4737 matrix_free.get_task_info().boundary_partition_data_hp,
4738 range_index);
4739 }
4740
4741 private:
4742 void
4743 process_range(const function_type &fu,
4744 const std::vector<unsigned int> &ptr,
4745 const std::vector<unsigned int> &data,
4746 const unsigned int range_index)
4747 {
4748 if (fu == nullptr)
4749 return;
4750
4751 AssertIndexRange(range_index + 1, ptr.size());
4752 for (unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4753 {
4754 AssertIndexRange(2 * i + 1, data.size());
4755 (container.*fu)(matrix_free,
4756 this->dst,
4757 this->src,
4758 std::make_pair(data[2 * i], data[2 * i + 1]));
4759 }
4760 }
4761
4762 public:
4763 // Starts the communication for the update ghost values operation. We
4764 // cannot call this update if ghost and destination are the same because
4765 // that would introduce spurious entries in the destination (there is also
4766 // the problem that reading from a vector that we also write to is usually
4767 // not intended in case there is overlap, but this is up to the
4768 // application code to decide and we cannot catch this case here).
4769 virtual void
4770 vector_update_ghosts_start() override
4771 {
4772 if (!src_and_dst_are_same)
4773 internal::update_ghost_values_start(src, src_data_exchanger);
4774 }
4775
4776 // Finishes the communication for the update ghost values operation
4777 virtual void
4778 vector_update_ghosts_finish() override
4779 {
4780 if (!src_and_dst_are_same)
4781 internal::update_ghost_values_finish(src, src_data_exchanger);
4782 }
4783
4784 // Starts the communication for the vector compress operation
4785 virtual void
4786 vector_compress_start() override
4787 {
4788 internal::compress_start(dst, dst_data_exchanger);
4789 }
4790
4791 // Finishes the communication for the vector compress operation
4792 virtual void
4793 vector_compress_finish() override
4794 {
4795 internal::compress_finish(dst, dst_data_exchanger);
4796 if (!src_and_dst_are_same)
4797 internal::reset_ghost_values(src, src_data_exchanger);
4798 }
4799
4800 // Zeros the given input vector
4801 virtual void
4802 zero_dst_vector_range(const unsigned int range_index) override
4803 {
4804 if (zero_dst_vector_setting)
4805 internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4806 }
4807
4808 virtual void
4809 cell_loop_pre_range(const unsigned int range_index) override
4810 {
4811 if (operation_before_loop)
4812 {
4814 matrix_free.get_dof_info(dof_handler_index_pre_post);
4815 if (range_index == numbers::invalid_unsigned_int)
4816 {
4817 // Case with threaded loop -> currently no overlap implemented
4819 0U,
4820 dof_info.vector_partitioner->locally_owned_size(),
4821 operation_before_loop,
4823 }
4824 else
4825 {
4826 AssertIndexRange(range_index,
4827 dof_info.cell_loop_pre_list_index.size() - 1);
4828 for (unsigned int id =
4829 dof_info.cell_loop_pre_list_index[range_index];
4830 id != dof_info.cell_loop_pre_list_index[range_index + 1];
4831 ++id)
4832 operation_before_loop(dof_info.cell_loop_pre_list[id].first,
4833 dof_info.cell_loop_pre_list[id].second);
4834 }
4835 }
4836 }
4837
4838 virtual void
4839 cell_loop_post_range(const unsigned int range_index) override
4840 {
4841 if (operation_after_loop)
4842 {
4843 // Run unit matrix operation on constrained dofs if we are at the
4844 // last range
4845 const std::vector<unsigned int> &partition_row_index =
4846 matrix_free.get_task_info().partition_row_index;
4847 if (range_index ==
4848 partition_row_index[partition_row_index.size() - 2] - 1)
4849 apply_operation_to_constrained_dofs(
4850 matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
4851 src,
4852 dst);
4853
4855 matrix_free.get_dof_info(dof_handler_index_pre_post);
4856 if (range_index == numbers::invalid_unsigned_int)
4857 {
4858 // Case with threaded loop -> currently no overlap implemented
4860 0U,
4861 dof_info.vector_partitioner->locally_owned_size(),
4862 operation_after_loop,
4864 }
4865 else
4866 {
4867 AssertIndexRange(range_index,
4868 dof_info.cell_loop_post_list_index.size() - 1);
4869 for (unsigned int id =
4870 dof_info.cell_loop_post_list_index[range_index];
4871 id != dof_info.cell_loop_post_list_index[range_index + 1];
4872 ++id)
4873 operation_after_loop(dof_info.cell_loop_post_list[id].first,
4874 dof_info.cell_loop_post_list[id].second);
4875 }
4876 }
4877 }
4878
4879 private:
4880 const MF &matrix_free;
4881 Container &container;
4882 function_type cell_function;
4883 function_type face_function;
4884 function_type boundary_function;
4885
4886 const InVector &src;
4887 OutVector &dst;
4888 VectorDataExchange<MF::dimension,
4889 typename MF::value_type,
4890 typename MF::vectorized_value_type>
4891 src_data_exchanger;
4892 VectorDataExchange<MF::dimension,
4893 typename MF::value_type,
4894 typename MF::vectorized_value_type>
4895 dst_data_exchanger;
4896 const bool src_and_dst_are_same;
4897 const bool zero_dst_vector_setting;
4898 const std::function<void(const unsigned int, const unsigned int)>
4899 operation_before_loop;
4900 const std::function<void(const unsigned int, const unsigned int)>
4901 operation_after_loop;
4902 const unsigned int dof_handler_index_pre_post;
4903 };
4904
4905
4906
4911 template <class MF, typename InVector, typename OutVector>
4912 struct MFClassWrapper
4913 {
4914 using function_type =
4915 std::function<void(const MF &,
4916 OutVector &,
4917 const InVector &,
4918 const std::pair<unsigned int, unsigned int> &)>;
4919
4920 MFClassWrapper(const function_type cell,
4921 const function_type face,
4922 const function_type boundary)
4923 : cell(cell)
4924 , face(face)
4925 , boundary(boundary)
4926 {}
4927
4928 void
4929 cell_integrator(const MF &mf,
4930 OutVector &dst,
4931 const InVector &src,
4932 const std::pair<unsigned int, unsigned int> &range) const
4933 {
4934 if (cell)
4935 cell(mf, dst, src, range);
4936 }
4937
4938 void
4939 face_integrator(const MF &mf,
4940 OutVector &dst,
4941 const InVector &src,
4942 const std::pair<unsigned int, unsigned int> &range) const
4943 {
4944 if (face)
4945 face(mf, dst, src, range);
4946 }
4947
4948 void
4949 boundary_integrator(
4950 const MF &mf,
4951 OutVector &dst,
4952 const InVector &src,
4953 const std::pair<unsigned int, unsigned int> &range) const
4954 {
4955 if (boundary)
4956 boundary(mf, dst, src, range);
4957 }
4958
4959 const function_type cell;
4960 const function_type face;
4961 const function_type boundary;
4962 };
4963
4964} // end of namespace internal
4965
4966
4967
4968template <int dim, typename Number, typename VectorizedArrayType>
4969template <typename OutVector, typename InVector>
4970inline void
4972 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
4973 OutVector &,
4974 const InVector &,
4975 const std::pair<unsigned int, unsigned int> &)>
4976 &cell_operation,
4977 OutVector &dst,
4978 const InVector &src,
4979 const bool zero_dst_vector) const
4980{
4981 using Wrapper =
4982 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4983 InVector,
4984 OutVector>;
4985 Wrapper wrap(cell_operation, nullptr, nullptr);
4986 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4987 InVector,
4988 OutVector,
4989 Wrapper,
4990 true>
4991 worker(*this,
4992 src,
4993 dst,
4994 zero_dst_vector,
4995 wrap,
4996 &Wrapper::cell_integrator,
4997 &Wrapper::face_integrator,
4998 &Wrapper::boundary_integrator);
4999
5000 task_info.loop(worker);
5001}
5002
5003
5004
5005template <int dim, typename Number, typename VectorizedArrayType>
5006template <typename OutVector, typename InVector>
5007inline void
5009 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5010 OutVector &,
5011 const InVector &,
5012 const std::pair<unsigned int, unsigned int> &)>
5013 &cell_operation,
5014 OutVector &dst,
5015 const InVector &src,
5016 const std::function<void(const unsigned int, const unsigned int)>
5017 &operation_before_loop,
5018 const std::function<void(const unsigned int, const unsigned int)>
5019 &operation_after_loop,
5020 const unsigned int dof_handler_index_pre_post) const
5021{
5022 using Wrapper =
5023 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5024 InVector,
5025 OutVector>;
5026 Wrapper wrap(cell_operation, nullptr, nullptr);
5027 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5028 InVector,
5029 OutVector,
5030 Wrapper,
5031 true>
5032 worker(*this,
5033 src,
5034 dst,
5035 false,
5036 wrap,
5037 &Wrapper::cell_integrator,
5038 &Wrapper::face_integrator,
5039 &Wrapper::boundary_integrator,
5040 DataAccessOnFaces::none,
5041 DataAccessOnFaces::none,
5042 operation_before_loop,
5043 operation_after_loop,
5044 dof_handler_index_pre_post);
5045
5046 task_info.loop(worker);
5047}
5048
5049
5050
5051template <int dim, typename Number, typename VectorizedArrayType>
5052template <typename OutVector, typename InVector>
5053inline void
5055 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5056 OutVector &,
5057 const InVector &,
5058 const std::pair<unsigned int, unsigned int> &)>
5059 &cell_operation,
5060 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5061 OutVector &,
5062 const InVector &,
5063 const std::pair<unsigned int, unsigned int> &)>
5064 &inner_face_operation,
5065 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5066 OutVector &,
5067 const InVector &,
5068 const std::pair<unsigned int, unsigned int> &)>
5069 &boundary_face_operation,
5070 OutVector &dst,
5071 const InVector &src,
5072 const bool zero_dst_vector,
5073 const DataAccessOnFaces dst_vector_face_access,
5074 const DataAccessOnFaces src_vector_face_access) const
5075{
5076 using Wrapper =
5077 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5078 InVector,
5079 OutVector>;
5080 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5081 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5082 InVector,
5083 OutVector,
5084 Wrapper,
5085 true>
5086 worker(*this,
5087 src,
5088 dst,
5089 zero_dst_vector,
5090 wrap,
5091 &Wrapper::cell_integrator,
5092 &Wrapper::face_integrator,
5093 &Wrapper::boundary_integrator,
5094 src_vector_face_access,
5095 dst_vector_face_access);
5096
5097 task_info.loop(worker);
5098}
5099
5100
5101
5102template <int dim, typename Number, typename VectorizedArrayType>
5103template <typename CLASS, typename OutVector, typename InVector>
5104inline void
5106 void (CLASS::*function_pointer)(
5107 const MatrixFree<dim, Number, VectorizedArrayType> &,
5108 OutVector &,
5109 const InVector &,
5110 const std::pair<unsigned int, unsigned int> &) const,
5111 const CLASS *owning_class,
5112 OutVector &dst,
5113 const InVector &src,
5114 const bool zero_dst_vector) const
5115{
5116 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5117 InVector,
5118 OutVector,
5119 CLASS,
5120 true>
5121 worker(*this,
5122 src,
5123 dst,
5124 zero_dst_vector,
5125 *owning_class,
5126 function_pointer,
5127 nullptr,
5128 nullptr);
5129 task_info.loop(worker);
5130}
5131
5132
5133
5134template <int dim, typename Number, typename VectorizedArrayType>
5135template <typename CLASS, typename OutVector, typename InVector>
5136inline void
5138 void (CLASS::*function_pointer)(
5139 const MatrixFree<dim, Number, VectorizedArrayType> &,
5140 OutVector &,
5141 const InVector &,
5142 const std::pair<unsigned int, unsigned int> &) const,
5143 const CLASS *owning_class,
5144 OutVector &dst,
5145 const InVector &src,
5146 const std::function<void(const unsigned int, const unsigned int)>
5147 &operation_before_loop,
5148 const std::function<void(const unsigned int, const unsigned int)>
5149 &operation_after_loop,
5150 const unsigned int dof_handler_index_pre_post) const
5151{
5152 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5153 InVector,
5154 OutVector,
5155 CLASS,
5156 true>
5157 worker(*this,
5158 src,
5159 dst,
5160 false,
5161 *owning_class,
5162 function_pointer,
5163 nullptr,
5164 nullptr,
5167 operation_before_loop,
5168 operation_after_loop,
5169 dof_handler_index_pre_post);
5170 task_info.loop(worker);
5171}
5172
5173
5174
5175template <int dim, typename Number, typename VectorizedArrayType>
5176template <typename CLASS, typename OutVector, typename InVector>
5177inline void
5179 void (CLASS::*cell_operation)(
5180 const MatrixFree<dim, Number, VectorizedArrayType> &,
5181 OutVector &,
5182 const InVector &,
5183 const std::pair<unsigned int, unsigned int> &) const,
5184 void (CLASS::*inner_face_operation)(
5185 const MatrixFree<dim, Number, VectorizedArrayType> &,
5186 OutVector &,
5187 const InVector &,
5188 const std::pair<unsigned int, unsigned int> &) const,
5189 void (CLASS::*boundary_face_operation)(
5190 const MatrixFree<dim, Number, VectorizedArrayType> &,
5191 OutVector &,
5192 const InVector &,
5193 const std::pair<unsigned int, unsigned int> &) const,
5194 const CLASS *owning_class,
5195 OutVector &dst,
5196 const InVector &src,
5197 const bool zero_dst_vector,
5198 const DataAccessOnFaces dst_vector_face_access,
5199 const DataAccessOnFaces src_vector_face_access) const
5200{
5201 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5202 InVector,
5203 OutVector,
5204 CLASS,
5205 true>
5206 worker(*this,
5207 src,
5208 dst,
5209 zero_dst_vector,
5210 *owning_class,
5211 cell_operation,
5212 inner_face_operation,
5213 boundary_face_operation,
5214 src_vector_face_access,
5215 dst_vector_face_access);
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> &),
5230 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 false>
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> &),
5262 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 false>
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> &),
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> &),
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> &),
5313 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 false>
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 OutVector, typename InVector>
5342inline void
5344 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5345 OutVector &,
5346 const InVector &,
5347 const std::pair<unsigned int, unsigned int> &)>
5348 &cell_operation,
5349 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5350 OutVector &,
5351 const InVector &,
5352 const std::pair<unsigned int, unsigned int> &)>
5353 &inner_face_operation,
5354 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5355 OutVector &,
5356 const InVector &,
5357 const std::pair<unsigned int, unsigned int> &)>
5358 &boundary_face_operation,
5359 OutVector &dst,
5360 const InVector &src,
5361 const std::function<void(const unsigned int, const unsigned int)>
5362 &operation_before_loop,
5363 const std::function<void(const unsigned int, const unsigned int)>
5364 &operation_after_loop,
5365 const unsigned int dof_handler_index_pre_post,
5366 const DataAccessOnFaces dst_vector_face_access,
5367 const DataAccessOnFaces src_vector_face_access) const
5368{
5369 using Wrapper =
5370 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5371 InVector,
5372 OutVector>;
5373 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5374 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5375 InVector,
5376 OutVector,
5377 Wrapper,
5378 true>
5379 worker(*this,
5380 src,
5381 dst,
5382 false,
5383 wrap,
5384 &Wrapper::cell_integrator,
5385 &Wrapper::face_integrator,
5386 &Wrapper::boundary_integrator,
5387 src_vector_face_access,
5388 dst_vector_face_access,
5389 operation_before_loop,
5390 operation_after_loop,
5391 dof_handler_index_pre_post);
5392
5393 task_info.loop(worker);
5394}
5395
5396
5397
5398template <int dim, typename Number, typename VectorizedArrayType>
5399template <typename CLASS, typename OutVector, typename InVector>
5400inline void
5402 void (CLASS::*cell_operation)(const MatrixFree &,
5403 OutVector &,
5404 const InVector &,
5405 const std::pair<unsigned int, unsigned int> &)
5406 const,
5407 void (CLASS::*inner_face_operation)(
5408 const MatrixFree &,
5409 OutVector &,
5410 const InVector &,
5411 const std::pair<unsigned int, unsigned int> &) const,
5412 void (CLASS::*boundary_face_operation)(
5413 const MatrixFree &,
5414 OutVector &,
5415 const InVector &,
5416 const std::pair<unsigned int, unsigned int> &) const,
5417 const CLASS *owning_class,
5418 OutVector &dst,
5419 const InVector &src,
5420 const std::function<void(const unsigned int, const unsigned int)>
5421 &operation_before_loop,
5422 const std::function<void(const unsigned int, const unsigned int)>
5423 &operation_after_loop,
5424 const unsigned int dof_handler_index_pre_post,
5425 const DataAccessOnFaces dst_vector_face_access,
5426 const DataAccessOnFaces src_vector_face_access) const
5427{
5428 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5429 InVector,
5430 OutVector,
5431 CLASS,
5432 true>
5433 worker(*this,
5434 src,
5435 dst,
5436 false,
5437 *owning_class,
5438 cell_operation,
5439 inner_face_operation,
5440 boundary_face_operation,
5441 src_vector_face_access,
5442 dst_vector_face_access,
5443 operation_before_loop,
5444 operation_after_loop,
5445 dof_handler_index_pre_post);
5446 task_info.loop(worker);
5447}
5448
5449
5450
5451template <int dim, typename Number, typename VectorizedArrayType>
5452template <typename CLASS, typename OutVector, typename InVector>
5453inline void
5455 void (CLASS::*cell_operation)(const MatrixFree &,
5456 OutVector &,
5457 const InVector &,
5458 const std::pair<unsigned int, unsigned int> &),
5459 void (CLASS::*inner_face_operation)(
5460 const MatrixFree &,
5461 OutVector &,
5462 const InVector &,
5463 const std::pair<unsigned int, unsigned int> &),
5464 void (CLASS::*boundary_face_operation)(
5465 const MatrixFree &,
5466 OutVector &,
5467 const InVector &,
5468 const std::pair<unsigned int, unsigned int> &),
5469 const CLASS *owning_class,
5470 OutVector &dst,
5471 const InVector &src,
5472 const std::function<void(const unsigned int, const unsigned int)>
5473 &operation_before_loop,
5474 const std::function<void(const unsigned int, const unsigned int)>
5475 &operation_after_loop,
5476 const unsigned int dof_handler_index_pre_post,
5477 const DataAccessOnFaces dst_vector_face_access,
5478 const DataAccessOnFaces src_vector_face_access) const
5479{
5480 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5481 InVector,
5482 OutVector,
5483 CLASS,
5484 false>
5485 worker(*this,
5486 src,
5487 dst,
5488 false,
5489 *owning_class,
5490 cell_operation,
5491 inner_face_operation,
5492 boundary_face_operation,
5493 src_vector_face_access,
5494 dst_vector_face_access,
5495 operation_before_loop,
5496 operation_after_loop,
5497 dof_handler_index_pre_post);
5498 task_info.loop(worker);
5499}
5500
5501
5502
5503template <int dim, typename Number, typename VectorizedArrayType>
5504template <typename CLASS, typename OutVector, typename InVector>
5505inline void
5507 void (CLASS::*function_pointer)(
5508 const MatrixFree<dim, Number, VectorizedArrayType> &,
5509 OutVector &,
5510 const InVector &,
5511 const std::pair<unsigned int, unsigned int> &) const,
5512 const CLASS *owning_class,
5513 OutVector &dst,
5514 const InVector &src,
5515 const bool zero_dst_vector,
5516 const DataAccessOnFaces src_vector_face_access) const
5517{
5518 auto src_vector_face_access_temp = src_vector_face_access;
5519 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5520 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5521 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5522 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5523
5524 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5525 InVector,
5526 OutVector,
5527 CLASS,
5528 true>
5529 worker(*this,
5530 src,
5531 dst,
5532 zero_dst_vector,
5533 *owning_class,
5534 function_pointer,
5535 nullptr,
5536 nullptr,
5537 src_vector_face_access_temp,
5539 task_info.loop(worker);
5540}
5541
5542
5543
5544template <int dim, typename Number, typename VectorizedArrayType>
5545template <typename CLASS, typename OutVector, typename InVector>
5546inline void
5548 void (CLASS::*function_pointer)(
5549 const MatrixFree<dim, Number, VectorizedArrayType> &,
5550 OutVector &,
5551 const InVector &,
5552 const std::pair<unsigned int, unsigned int> &),
5553 CLASS *owning_class,
5554 OutVector &dst,
5555 const InVector &src,
5556 const bool zero_dst_vector,
5557 const DataAccessOnFaces src_vector_face_access) const
5558{
5559 auto src_vector_face_access_temp = src_vector_face_access;
5560 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5561 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5562 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5563 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5564
5565 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5566 InVector,
5567 OutVector,
5568 CLASS,
5569 false>
5570 worker(*this,
5571 src,
5572 dst,
5573 zero_dst_vector,
5574 *owning_class,
5575 function_pointer,
5576 nullptr,
5577 nullptr,
5578 src_vector_face_access_temp,
5580 task_info.loop(worker);
5581}
5582
5583
5584
5585template <int dim, typename Number, typename VectorizedArrayType>
5586template <typename OutVector, typename InVector>
5587inline void
5589 const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType> &,
5590 OutVector &,
5591 const InVector &,
5592 const std::pair<unsigned int, unsigned int> &)>
5593 &cell_operation,
5594 OutVector &dst,
5595 const InVector &src,
5596 const bool zero_dst_vector,
5597 const DataAccessOnFaces src_vector_face_access) const
5598{
5599 auto src_vector_face_access_temp = src_vector_face_access;
5600 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5601 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5602 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5603 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5604
5605 using Wrapper =
5606 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5607 InVector,
5608 OutVector>;
5609 Wrapper wrap(cell_operation, nullptr, nullptr);
5610
5611 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5612 InVector,
5613 OutVector,
5614 Wrapper,
5615 true>
5616 worker(*this,
5617 src,
5618 dst,
5619 zero_dst_vector,
5620 wrap,
5621 &Wrapper::cell_integrator,
5622 &Wrapper::face_integrator,
5623 &Wrapper::boundary_integrator,
5624 src_vector_face_access_temp,
5625 DataAccessOnFaces::none);
5626 task_info.loop(worker);
5627}
5628
5629
5630#endif // ifndef DOXYGEN
5631
5632
5633
5635
5636#endif
iterator begin()
void resize(const size_type new_size)
void reinit(const size_type size, const bool omit_zeroing_entries=false)
Abstract base class for mapping classes.
Definition mapping.h:318
unsigned int n_active_fe_indices() const
std::pair< typename DoFHandler< dim >::cell_iterator, unsigned int > get_face_iterator(const unsigned int face_batch_index, const unsigned int lane_index, const bool interior=true, const unsigned int fe_component=0) const
void loop_cell_centric(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > mapping_info
unsigned int get_cell_active_fe_index(const std::pair< unsigned int, unsigned int > range) const
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
unsigned int n_ghost_cell_batches() const
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
void initialize_dof_handlers(const std::vector< const DoFHandler< dim, dim > * > &dof_handlers, const AdditionalData &additional_data)
types::boundary_id get_boundary_id(const unsigned int face_batch_index) const
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > > shape_info
void print(std::ostream &out) const
void update_mapping(const std::shared_ptr< hp::MappingCollection< dim > > &mapping)
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
unsigned int n_inner_face_batches() const
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
void update_mapping(const Mapping< dim > &mapping)
const Quadrature< dim > & get_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_index) const
unsigned int get_mg_level() const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*boundary_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void print_memory_consumption(StreamType &out) const
bool mapping_initialized() const
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
void clear()
void internal_reinit(const std::shared_ptr< hp::MappingCollection< dim > > &mapping, const std::vector< const DoFHandler< dim, dim > * > &dof_handlers, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const std::vector< hp::QCollection< q_dim > > &quad, const AdditionalData &additional_data)
unsigned int get_n_q_points_face(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &inner_face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_face_operation, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
bool at_irregular_cell(const unsigned int cell_batch_index) const
~MatrixFree() override=default
const AffineConstraints< Number > & get_affine_constraints(const unsigned int dof_handler_index=0) const
void copy_from(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free_base)
void initialize_face_data_vector(AlignedVector< T > &vec) const
unsigned int get_dofs_per_cell(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int n_constraint_pool_entries() const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
const Number * constraint_pool_begin(const unsigned int pool_index) const
AlignedVector< Number > * acquire_scratch_data_non_threadsafe() const
const IndexSet & get_locally_owned_set(const unsigned int dof_handler_index=0) const
unsigned int get_dofs_per_face(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
internal::MatrixFreeFunctions::TaskInfo task_info
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
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
bool mapping_is_initialized
void loop_cell_centric(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void initialize_dof_vector(LinearAlgebra::distributed::Vector< Number2 > &vec, const unsigned int dof_handler_index=0) const
DoFHandler< dim >::cell_iterator get_cell_iterator(const unsigned int cell_batch_index, const unsigned int lane_index, const unsigned int dof_handler_index=0) const
MatrixFree(const MatrixFree< dim, Number, VectorizedArrayType > &other)
std::pair< unsigned int, unsigned int > create_cell_subrange_hp(const std::pair< unsigned int, unsigned int > &range, const unsigned int fe_degree, const unsigned int dof_handler_index=0) const
std::pair< unsigned int, unsigned int > create_cell_subrange_hp_by_index(const std::pair< unsigned int, unsigned int > &range, const unsigned int fe_index, const unsigned int dof_handler_index=0) const
unsigned int mg_level
std::pair< unsigned int, unsigned int > get_face_range_category(const std::pair< unsigned int, unsigned int > face_batch_range) const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
std::vector< std::pair< unsigned int, unsigned int > > cell_level_index
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner(const unsigned int dof_handler_index=0) const
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
std::vector< Number > constraint_pool_data
VectorizedArrayType vectorized_value_type
unsigned int n_cell_batches() const
const IndexSet & get_ghost_set(const unsigned int dof_handler_index=0) const
std::pair< int, int > get_cell_level_and_index(const unsigned int cell_batch_index, const unsigned int lane_index) const
bool indices_initialized() const
const Quadrature< dim - 1 > & get_face_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int get_cell_category(const unsigned int cell_batch_index) const
Threads::ThreadLocalStorage< std::list< std::pair< bool, AlignedVector< VectorizedArrayType > > > > scratch_pad
unsigned int get_matrix_free_cell_index(const typename Triangulation< dim >::cell_iterator &cell) const
Number value_type
std::size_t memory_consumption() const
std::pair< unsigned int, unsigned int > get_face_category(const unsigned int face_batch_index) const
unsigned int n_boundary_face_batches() const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
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
std::vector< SmartPointer< const DoFHandler< dim > > > dof_handlers
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)
std::vector< SmartPointer< const AffineConstraints< Number > > > affine_constraints
unsigned int n_components() const
unsigned int n_physical_cells() const
void initialize_dof_vector(VectorType &vec, const unsigned int dof_handler_index=0) const
unsigned int get_n_q_points(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
void reinit(const MappingType &mapping, const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
std::vector< unsigned int > mf_cell_indices
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*boundary_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
std::vector< internal::MatrixFreeFunctions::DoFInfo > dof_info
void initialize_cell_data_vector(AlignedVector< T > &vec) const
unsigned int n_ghost_inner_face_batches() const
void release_scratch_data_non_threadsafe(const AlignedVector< Number > *memory) const
void renumber_dofs(std::vector< types::global_dof_index > &renumbering, const unsigned int dof_handler_index=0)
unsigned int get_face_active_fe_index(const std::pair< unsigned int, unsigned int > range, const bool is_interior_face=true) const
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
bool indices_are_initialized
internal::MatrixFreeFunctions::FaceInfo< VectorizedArrayType::size()> face_info
void reinit(const MappingType &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< number2 > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
unsigned int get_cell_range_category(const std::pair< unsigned int, unsigned int > cell_batch_range) const
void reinit(const MappingType &mapping, const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< QuadratureType > &quad, const AdditionalData &additional_data=AdditionalData())
std::array< types::boundary_id, VectorizedArrayType::size()> get_faces_by_cells_boundary_id(const unsigned int cell_batch_index, const unsigned int face_number) const
std::vector< unsigned int > constraint_pool_row_index
const 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
static bool is_supported(const FiniteElement< dim, spacedim > &fe)
static constexpr unsigned int dimension
A class that provides a separate storage location on each thread that accesses the object.
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
unsigned int level
Definition grid_out.cc:4626
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ActiveSelector::cell_iterator cell_iterator
UpdateFlags
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_default
No update.
for(unsigned int j=best_vertex+1;j< vertices.size();++j) if(vertices_to_use[j])
bool job_supports_mpi()
Definition mpi.cc:697
unsigned int minimum_parallel_grain_size
Definition parallel.cc:33
const types::boundary_id invalid_boundary_id
Definition types.h:292
static const unsigned int invalid_unsigned_int
Definition types.h:220
void apply_to_subranges(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, const Function &f, const unsigned int grainsize)
Definition parallel.h:452
const InputIterator OutputIterator const Function & function
Definition parallel.h:168
const Iterator & begin
Definition parallel.h:609
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
unsigned short int fe_index
Definition types.h:59
unsigned int boundary_id
Definition types.h:144
TasksParallelScheme tasks_parallel_scheme
UpdateFlags mapping_update_flags_inner_faces
AdditionalData(const TasksParallelScheme tasks_parallel_scheme=partition_partition, const unsigned int tasks_block_size=0, const UpdateFlags mapping_update_flags=update_gradients|update_JxW_values, const UpdateFlags mapping_update_flags_boundary_faces=update_default, const UpdateFlags mapping_update_flags_inner_faces=update_default, const UpdateFlags mapping_update_flags_faces_by_cells=update_default, const unsigned int mg_level=numbers::invalid_unsigned_int, const bool store_plain_indices=true, const bool initialize_indices=true, const bool initialize_mapping=true, const bool overlap_communication_computation=true, const bool hold_all_faces_to_owned_cells=false, const bool cell_vectorization_categories_strict=false, const bool allow_ghosted_vectors_in_loops=true)
std::vector< unsigned int > cell_vectorization_category
AdditionalData & operator=(const AdditionalData &other)
UpdateFlags mapping_update_flags_boundary_faces
UpdateFlags mapping_update_flags_faces_by_cells
AdditionalData(const AdditionalData &other)
std::vector< unsigned int > cell_loop_pre_list_index
Definition dof_info.h:748
std::vector< unsigned int > cell_loop_post_list_index
Definition dof_info.h:761
std::vector< std::pair< unsigned int, unsigned int > > vector_zero_range_list
Definition dof_info.h:741
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition dof_info.h:592
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_pre_list
Definition dof_info.h:754
std::vector< unsigned int > vector_zero_range_list_index
Definition dof_info.h:736
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_post_list
Definition dof_info.h:767
void loop(MFWorkerInterface &worker) const
Definition task_info.cc:350