14#ifndef dealii_matrix_free_h
15#define dealii_matrix_free_h
47#include <boost/container/small_vector.hpp>
116 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
117 "Type of Number and of VectorizedArrayType do not match.");
600 template <
typename QuadratureType,
typename number2,
typename MappingType>
605 const QuadratureType &quad,
629 template <
typename QuadratureType,
typename number2,
typename MappingType>
634 const std::vector<QuadratureType> &quad,
644 template <
typename QuadratureType,
typename number2,
typename MappingType>
649 const QuadratureType &quad,
806 template <
typename OutVector,
typename InVector>
812 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
815 const bool zero_dst_vector =
false)
const;
863 template <
typename CLASS,
typename OutVector,
typename InVector>
869 const
std::pair<unsigned
int, unsigned
int> &) const,
870 const CLASS *owning_class,
873 const
bool zero_dst_vector = false) const;
878 template <
typename CLASS,
typename OutVector,
typename InVector>
884 const
std::pair<unsigned
int, unsigned
int> &),
888 const
bool zero_dst_vector = false) const;
974 template <
typename CLASS,
typename OutVector,
typename InVector>
980 const
std::pair<unsigned
int, unsigned
int> &) const,
981 const CLASS *owning_class,
984 const
std::function<void(const unsigned
int, const unsigned
int)>
985 &operation_before_loop,
986 const
std::function<void(const unsigned
int, const unsigned
int)>
987 &operation_after_loop,
988 const unsigned
int dof_handler_index_pre_post = 0) const;
993 template <
typename CLASS,
typename OutVector,
typename InVector>
999 const
std::pair<unsigned
int, unsigned
int> &),
1000 CLASS *owning_class,
1002 const InVector &src,
1003 const
std::function<void(const unsigned
int, const unsigned
int)>
1004 &operation_before_loop,
1005 const
std::function<void(const unsigned
int, const unsigned
int)>
1006 &operation_after_loop,
1007 const unsigned
int dof_handler_index_pre_post = 0) const;
1013 template <
typename OutVector,
typename InVector>
1019 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1021 const InVector &src,
1022 const std::function<
void(
const unsigned int,
const unsigned int)>
1023 &operation_before_loop,
1024 const std::function<
void(
const unsigned int,
const unsigned int)>
1025 &operation_after_loop,
1026 const unsigned int dof_handler_index_pre_post = 0)
const;
1103 template <
typename OutVector,
typename InVector>
1106 const std::function<
1110 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1111 const std::function<
void(
1115 const std::pair<unsigned int, unsigned int> &)> &inner_face_operation,
1116 const std::function<
void(
1120 const std::pair<unsigned int, unsigned int> &)> &boundary_face_operation,
1122 const InVector &src,
1123 const bool zero_dst_vector =
false,
1217 template <
typename CLASS,
typename OutVector,
typename InVector>
1219 loop(
void (CLASS::*cell_operation)(
1223 const
std::pair<unsigned
int, unsigned
int> &) const,
1224 void (CLASS::*inner_face_operation)(
1228 const
std::pair<unsigned
int, unsigned
int> &) const,
1229 void (CLASS::*boundary_face_operation)(
1233 const
std::pair<unsigned
int, unsigned
int> &) const,
1234 const CLASS *owning_class,
1236 const InVector &src,
1237 const
bool zero_dst_vector = false,
1246 template <
typename CLASS,
typename OutVector,
typename InVector>
1248 loop(
void (CLASS::*cell_operation)(
1252 const
std::pair<unsigned
int, unsigned
int> &),
1253 void (CLASS::*inner_face_operation)(
1257 const
std::pair<unsigned
int, unsigned
int> &),
1258 void (CLASS::*boundary_face_operation)(
1262 const
std::pair<unsigned
int, unsigned
int> &),
1263 CLASS *owning_class,
1265 const InVector &src,
1266 const
bool zero_dst_vector = false,
1393 template <
typename CLASS,
typename OutVector,
typename InVector>
1395 loop(
void (CLASS::*cell_operation)(
1399 const
std::pair<unsigned
int, unsigned
int> &) const,
1400 void (CLASS::*inner_face_operation)(
1404 const
std::pair<unsigned
int, unsigned
int> &) const,
1405 void (CLASS::*boundary_face_operation)(
1409 const
std::pair<unsigned
int, unsigned
int> &) const,
1410 const CLASS *owning_class,
1412 const InVector &src,
1413 const
std::function<void(const unsigned
int, const unsigned
int)>
1414 &operation_before_loop,
1415 const
std::function<void(const unsigned
int, const unsigned
int)>
1416 &operation_after_loop,
1417 const unsigned
int dof_handler_index_pre_post = 0,
1426 template <
typename CLASS,
typename OutVector,
typename InVector>
1428 loop(
void (CLASS::*cell_operation)(
1432 const
std::pair<unsigned
int, unsigned
int> &),
1433 void (CLASS::*inner_face_operation)(
1437 const
std::pair<unsigned
int, unsigned
int> &),
1438 void (CLASS::*boundary_face_operation)(
1442 const
std::pair<unsigned
int, unsigned
int> &),
1443 const CLASS *owning_class,
1445 const InVector &src,
1446 const
std::function<void(const unsigned
int, const unsigned
int)>
1447 &operation_before_loop,
1448 const
std::function<void(const unsigned
int, const unsigned
int)>
1449 &operation_after_loop,
1450 const unsigned
int dof_handler_index_pre_post = 0,
1461 template <
typename OutVector,
typename InVector>
1464 const std::function<
1468 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1469 const std::function<
void(
1473 const std::pair<unsigned int, unsigned int> &)> &inner_face_operation,
1474 const std::function<
void(
1478 const std::pair<unsigned int, unsigned int> &)> &boundary_face_operation,
1480 const InVector &src,
1481 const std::function<
void(
const unsigned int,
const unsigned int)>
1482 &operation_before_loop,
1483 const std::function<
void(
const unsigned int,
const unsigned int)>
1484 &operation_after_loop,
1485 const unsigned int dof_handler_index_pre_post = 0,
1555 template <
typename CLASS,
typename OutVector,
typename InVector>
1561 const
std::pair<unsigned
int, unsigned
int> &) const,
1562 const CLASS *owning_class,
1564 const InVector &src,
1565 const
bool zero_dst_vector = false,
1572 template <
typename CLASS,
typename OutVector,
typename InVector>
1578 const
std::pair<unsigned
int, unsigned
int> &),
1579 CLASS *owning_class,
1581 const InVector &src,
1582 const
bool zero_dst_vector = false,
1589 template <
typename OutVector,
typename InVector>
1595 const std::pair<unsigned int, unsigned int> &)>
1598 const InVector &src,
1599 const bool zero_dst_vector =
false,
1610 std::pair<unsigned int, unsigned int>
1612 const unsigned int fe_degree,
1613 const unsigned int dof_handler_index = 0)
const;
1621 std::pair<unsigned int, unsigned int>
1623 const std::pair<unsigned int, unsigned int> &range,
1624 const unsigned int fe_index,
1625 const unsigned int dof_handler_index = 0)
const;
1638 const std::pair<unsigned int, unsigned int> range,
1646 const std::pair<unsigned int, unsigned int> range,
1647 const bool is_interior_face =
true,
1661 template <
typename T>
1670 template <
typename T>
1687 template <
typename VectorType>
1690 const unsigned int dof_handler_index = 0)
const;
1709 template <
typename Number2,
typename MemorySpace>
1713 const unsigned int dof_handler_index = 0)
const;
1723 template <
typename Number2,
typename MemorySpace>
1738 const std::shared_ptr<const Utilities::MPI::Partitioner> &
1762 const std::vector<unsigned int> &
1777 const unsigned int dof_handler_index = 0);
1788 template <
int spacedim>
1887 const unsigned int face_number)
const;
1919 const unsigned int lane_index,
1920 const unsigned int dof_handler_index = 0)
const;
1929 const unsigned int lane_index)
const;
1953 std::pair<typename DoFHandler<dim>::cell_iterator,
unsigned int>
1955 const unsigned int lane_index,
1956 const bool interior =
true,
1957 const unsigned int fe_component = 0)
const;
2003 const unsigned int hp_active_fe_index = 0)
const;
2010 const unsigned int hp_active_fe_index = 0)
const;
2018 const unsigned int hp_active_fe_index = 0)
const;
2026 const unsigned int hp_active_fe_index = 0)
const;
2033 const unsigned int hp_active_fe_index = 0)
const;
2040 const unsigned int hp_active_fe_index = 0)
const;
2056 const std::pair<unsigned int, unsigned int> cell_batch_range,
2063 std::pair<unsigned int, unsigned int>
2065 const std::pair<unsigned int, unsigned int> face_batch_range,
2081 const unsigned int cell_batch_index,
2088 std::pair<unsigned int, unsigned int>
2090 const unsigned int face_batch_index,
2124 template <
typename StreamType>
2152 const internal::MatrixFreeFunctions::
2153 MappingInfo<dim, Number, VectorizedArrayType> &
2188 const unsigned int quad_index = 0,
2189 const unsigned int fe_base_element = 0,
2190 const unsigned int hp_active_fe_index = 0,
2191 const unsigned int hp_active_quad_index = 0)
const;
2197 VectorizedArrayType::size()> &
2257 template <
typename number2,
int q_dim>
2263 const std::vector<IndexSet> &locally_owned_set,
2273 template <
typename number2>
2277 const std::vector<IndexSet> &locally_owned_set,
2299 std::vector<ObserverPointer<const AffineConstraints<Number>>>
2306 std::vector<internal::MatrixFreeFunctions::DoFInfo>
dof_info;
2388 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2395 mutable std::list<std::pair<bool, AlignedVector<Number>>>
2418template <
int dim,
typename Number,
typename VectorizedArrayType>
2419template <
typename T>
2424 vec.
resize(this->n_cell_batches() + this->n_ghost_cell_batches());
2429template <
int dim,
typename Number,
typename VectorizedArrayType>
2430template <
typename T>
2435 vec.
resize(this->n_inner_face_batches() + this->n_boundary_face_batches() +
2436 this->n_ghost_inner_face_batches());
2441template <
int dim,
typename Number,
typename VectorizedArrayType>
2442template <
typename VectorType>
2446 const unsigned int comp)
const
2449 "This function is not supported for block vectors.");
2451 Assert(task_info.n_procs == 1,
2452 ExcMessage(
"This function can only be used in serial."));
2455 vec.reinit(dof_info[comp].vector_partitioner->size());
2460template <
int dim,
typename Number,
typename VectorizedArrayType>
2461template <
typename Number2,
typename MemorySpace>
2465 const unsigned int comp)
const
2468 vec.
reinit(dof_info[comp].vector_partitioner, task_info.communicator_sm);
2473template <
int dim,
typename Number,
typename VectorizedArrayType>
2474template <
typename Number2,
typename MemorySpace>
2479 vec.
reinit(n_components());
2480 for (
unsigned int c = 0; c < n_components(); ++c)
2481 this->initialize_dof_vector(vec.
block(c), c);
2488template <
int dim,
typename Number,
typename VectorizedArrayType>
2489inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2491 const unsigned int comp)
const
2494 return dof_info[comp].vector_partitioner;
2499template <
int dim,
typename Number,
typename VectorizedArrayType>
2500inline const std::vector<unsigned int> &
2502 const unsigned int comp)
const
2505 return dof_info[comp].constrained_dofs;
2510template <
int dim,
typename Number,
typename VectorizedArrayType>
2515 return dof_handlers.size();
2520template <
int dim,
typename Number,
typename VectorizedArrayType>
2523 const unsigned int dof_handler_index)
const
2527 return dof_handlers[dof_handler_index]->get_fe().n_base_elements();
2532template <
int dim,
typename Number,
typename VectorizedArrayType>
2541template <
int dim,
typename Number,
typename VectorizedArrayType>
2545 return task_info.n_active_cells;
2550template <
int dim,
typename Number,
typename VectorizedArrayType>
2554 return *(task_info.cell_partition_data.end() - 2);
2559template <
int dim,
typename Number,
typename VectorizedArrayType>
2563 return *(task_info.cell_partition_data.end() - 1) -
2564 *(task_info.cell_partition_data.end() - 2);
2569template <
int dim,
typename Number,
typename VectorizedArrayType>
2573 if (task_info.face_partition_data.empty())
2575 return task_info.face_partition_data.back();
2580template <
int dim,
typename Number,
typename VectorizedArrayType>
2584 if (task_info.face_partition_data.empty())
2586 return task_info.boundary_partition_data.back() -
2587 task_info.face_partition_data.back();
2592template <
int dim,
typename Number,
typename VectorizedArrayType>
2596 if (task_info.face_partition_data.empty())
2598 return face_info.faces.size() - task_info.boundary_partition_data.back();
2603template <
int dim,
typename Number,
typename VectorizedArrayType>
2606 const unsigned int face_batch_index)
const
2608 Assert(face_batch_index >= task_info.boundary_partition_data[0] &&
2609 face_batch_index < task_info.boundary_partition_data.back(),
2611 task_info.boundary_partition_data[0],
2612 task_info.boundary_partition_data.back()));
2618template <
int dim,
typename Number,
typename VectorizedArrayType>
2621 const unsigned int cell_batch_index,
2622 const unsigned int face_number)
const
2626 Assert(face_info.cell_and_face_boundary_id.size(0) >= n_cell_batches(),
2630 for (
unsigned int v = 0;
2631 v < n_active_entries_per_cell_batch(cell_batch_index);
2634 face_info.cell_and_face_boundary_id(cell_batch_index, face_number, v);
2640template <
int dim,
typename Number,
typename VectorizedArrayType>
2641inline const internal::MatrixFreeFunctions::
2642 MappingInfo<dim, Number, VectorizedArrayType> &
2645 return mapping_info;
2650template <
int dim,
typename Number,
typename VectorizedArrayType>
2653 const unsigned int dof_index)
const
2656 return dof_info[dof_index];
2661template <
int dim,
typename Number,
typename VectorizedArrayType>
2665 return constraint_pool_row_index.size() - 1;
2670template <
int dim,
typename Number,
typename VectorizedArrayType>
2673 const unsigned int row)
const
2676 return constraint_pool_data.empty() ?
2678 constraint_pool_data.data() + constraint_pool_row_index[row];
2683template <
int dim,
typename Number,
typename VectorizedArrayType>
2686 const unsigned int row)
const
2689 return constraint_pool_data.empty() ?
2691 constraint_pool_data.data() + constraint_pool_row_index[row + 1];
2696template <
int dim,
typename Number,
typename VectorizedArrayType>
2697inline std::pair<unsigned int, unsigned int>
2699 const std::pair<unsigned int, unsigned int> &range,
2700 const unsigned int degree,
2701 const unsigned int dof_handler_component)
const
2703 if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2706 dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2708 dof_info[dof_handler_component].fe_index_conversion[0].
size(), 1);
2709 if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2712 return {range.second, range.second};
2716 dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2717 if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2718 return {range.second, range.second};
2720 return create_cell_subrange_hp_by_index(range,
2722 dof_handler_component);
2727template <
int dim,
typename Number,
typename VectorizedArrayType>
2730 const unsigned int cell_batch_index)
const
2733 return VectorizedArrayType::size() > 1 &&
2734 cell_level_index[(cell_batch_index + 1) * VectorizedArrayType::size() -
2735 1] == cell_level_index[(cell_batch_index + 1) *
2736 VectorizedArrayType::size() -
2742template <
int dim,
typename Number,
typename VectorizedArrayType>
2746 return shape_info.size(2);
2750template <
int dim,
typename Number,
typename VectorizedArrayType>
2753 const std::pair<unsigned int, unsigned int> range,
2754 const unsigned int dof_handler_index)
const
2756 const unsigned int dof_handler_index_local =
2758 first_hp_dof_handler_index :
2761 const auto &fe_indices =
2762 dof_info[dof_handler_index_local].cell_active_fe_index;
2764 if (fe_indices.empty() ==
true ||
2765 dof_handlers[dof_handler_index_local]->get_fe_collection().size() == 1)
2768 const auto index = fe_indices[range.first];
2770 for (
unsigned int i = range.first; i < range.second; ++i)
2778template <
int dim,
typename Number,
typename VectorizedArrayType>
2781 const std::pair<unsigned int, unsigned int> range,
2782 const bool is_interior_face,
2783 const unsigned int dof_handler_index)
const
2785 const unsigned int dof_handler_index_local =
2787 first_hp_dof_handler_index :
2790 const auto &fe_indices =
2791 dof_info[dof_handler_index_local].cell_active_fe_index;
2793 if (fe_indices.empty() ==
true)
2796 if (is_interior_face)
2798 const unsigned int index =
2799 fe_indices[face_info.faces[range.first].cells_interior[0] /
2800 VectorizedArrayType::size()];
2802 for (
unsigned int i = range.first; i < range.second; ++i)
2804 fe_indices[face_info.faces[i].cells_interior[0] /
2805 VectorizedArrayType::size()]);
2811 const unsigned int index =
2812 fe_indices[face_info.faces[range.first].cells_exterior[0] /
2813 VectorizedArrayType::size()];
2815 for (
unsigned int i = range.first; i < range.second; ++i)
2817 fe_indices[face_info.faces[i].cells_exterior[0] /
2818 VectorizedArrayType::size()]);
2826template <
int dim,
typename Number,
typename VectorizedArrayType>
2829 const unsigned int cell_batch_index)
const
2833 const std::vector<unsigned char> &n_lanes_filled =
2834 dof_info[0].n_vectorization_lanes_filled
2838 return n_lanes_filled[cell_batch_index];
2843template <
int dim,
typename Number,
typename VectorizedArrayType>
2846 const unsigned int face_batch_index)
const
2850 const std::vector<unsigned char> &n_lanes_filled =
2851 dof_info[0].n_vectorization_lanes_filled
2854 return n_lanes_filled[face_batch_index];
2859template <
int dim,
typename Number,
typename VectorizedArrayType>
2862 const unsigned int dof_handler_index,
2863 const unsigned int active_fe_index)
const
2865 return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2870template <
int dim,
typename Number,
typename VectorizedArrayType>
2873 const unsigned int quad_index,
2874 const unsigned int active_fe_index)
const
2877 return mapping_info.cell_data[quad_index]
2878 .descriptor[active_fe_index]
2884template <
int dim,
typename Number,
typename VectorizedArrayType>
2887 const unsigned int dof_handler_index,
2888 const unsigned int active_fe_index)
const
2890 return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2895template <
int dim,
typename Number,
typename VectorizedArrayType>
2898 const unsigned int quad_index,
2899 const unsigned int active_fe_index)
const
2902 return mapping_info.face_data[quad_index]
2903 .descriptor[active_fe_index]
2909template <
int dim,
typename Number,
typename VectorizedArrayType>
2912 const unsigned int dof_handler_index)
const
2914 return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2919template <
int dim,
typename Number,
typename VectorizedArrayType>
2922 const unsigned int dof_handler_index)
const
2924 return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2929template <
int dim,
typename Number,
typename VectorizedArrayType>
2932 const unsigned int dof_handler_index,
2933 const unsigned int index_quad,
2934 const unsigned int index_fe,
2935 const unsigned int active_fe_index,
2936 const unsigned int active_quad_index)
const
2939 const unsigned int ind =
2940 dof_info[dof_handler_index].global_base_element_offset + index_fe;
2945 return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2950template <
int dim,
typename Number,
typename VectorizedArrayType>
2952 VectorizedArrayType::size()> &
2954 const unsigned int face_batch_index)
const
2957 return face_info.faces[face_batch_index];
2962template <
int dim,
typename Number,
typename VectorizedArrayType>
2967 return face_info.cell_and_face_to_plain_faces;
2972template <
int dim,
typename Number,
typename VectorizedArrayType>
2975 const unsigned int quad_index,
2976 const unsigned int active_fe_index)
const
2979 return mapping_info.cell_data[quad_index]
2980 .descriptor[active_fe_index]
2986template <
int dim,
typename Number,
typename VectorizedArrayType>
2989 const unsigned int quad_index,
2990 const unsigned int active_fe_index)
const
2993 return mapping_info.face_data[quad_index]
2994 .descriptor[active_fe_index]
3000template <
int dim,
typename Number,
typename VectorizedArrayType>
3003 const std::pair<unsigned int, unsigned int> range,
3004 const unsigned int dof_handler_index)
const
3006 auto result = get_cell_category(range.first, dof_handler_index);
3008 for (
unsigned int i = range.first; i < range.second; ++i)
3009 result =
std::max(result, get_cell_category(i, dof_handler_index));
3016template <
int dim,
typename Number,
typename VectorizedArrayType>
3017inline std::pair<unsigned int, unsigned int>
3019 const std::pair<unsigned int, unsigned int> range,
3020 const unsigned int dof_handler_index)
const
3022 auto result = get_face_category(range.first, dof_handler_index);
3024 for (
unsigned int i = range.first; i < range.second; ++i)
3027 std::max(result.first, get_face_category(i, dof_handler_index).
first);
3029 std::max(result.second, get_face_category(i, dof_handler_index).
second);
3037template <
int dim,
typename Number,
typename VectorizedArrayType>
3040 const unsigned int cell_batch_index,
3041 const unsigned int dof_handler_index)
const
3045 const unsigned int dof_handler_index_local =
3047 first_hp_dof_handler_index :
3052 dof_info[dof_handler_index_local].cell_active_fe_index.size());
3053 if (dof_info[dof_handler_index_local].cell_active_fe_index.empty())
3056 return dof_info[dof_handler_index_local]
3057 .cell_active_fe_index[cell_batch_index];
3062template <
int dim,
typename Number,
typename VectorizedArrayType>
3063inline std::pair<unsigned int, unsigned int>
3065 const unsigned int face_batch_index,
3066 const unsigned int dof_handler_index)
const
3068 const unsigned int dof_handler_index_local =
3070 first_hp_dof_handler_index :
3074 if (dof_info[dof_handler_index_local].cell_active_fe_index.empty())
3075 return std::make_pair(0U, 0U);
3077 std::pair<unsigned int, unsigned int> result = std::make_pair(0U, 0U);
3078 for (
unsigned int v = 0;
3079 v < VectorizedArrayType::size() &&
3080 face_info.faces[face_batch_index].cells_interior[v] !=
3085 dof_info[dof_handler_index_local].cell_active_fe_index
3086 [face_info.faces[face_batch_index].cells_interior[v] /
3087 VectorizedArrayType::size()]);
3088 if (face_info.faces[face_batch_index].cells_exterior[0] !=
3090 for (
unsigned int v = 0;
3091 v < VectorizedArrayType::size() &&
3092 face_info.faces[face_batch_index].cells_exterior[v] !=
3097 dof_info[dof_handler_index_local].cell_active_fe_index
3098 [face_info.faces[face_batch_index].cells_exterior[v] /
3099 VectorizedArrayType::size()]);
3107template <
int dim,
typename Number,
typename VectorizedArrayType>
3111 return indices_are_initialized;
3116template <
int dim,
typename Number,
typename VectorizedArrayType>
3120 return mapping_is_initialized;
3124template <
int dim,
typename Number,
typename VectorizedArrayType>
3133template <
int dim,
typename Number,
typename VectorizedArrayType>
3138 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3139 list_type &
data = scratch_pad.get();
3140 for (
typename list_type::iterator it =
data.begin(); it !=
data.end(); ++it)
3141 if (it->first ==
false)
3147 return &
data.front().second;
3152template <
int dim,
typename Number,
typename VectorizedArrayType>
3158 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3159 list_type &
data = scratch_pad.get();
3160 for (
typename list_type::iterator it =
data.begin(); it !=
data.end(); ++it)
3161 if (&it->second == scratch)
3172template <
int dim,
typename Number,
typename VectorizedArrayType>
3178 scratch_pad_non_threadsafe.
begin();
3179 it != scratch_pad_non_threadsafe.end();
3181 if (it->first ==
false)
3186 scratch_pad_non_threadsafe.push_front(
3188 return &scratch_pad_non_threadsafe.front().second;
3193template <
int dim,
typename Number,
typename VectorizedArrayType>
3200 scratch_pad_non_threadsafe.
begin();
3201 it != scratch_pad_non_threadsafe.end();
3203 if (&it->second == scratch)
3218 namespace MatrixFreeImplementation
3220 template <
int dim,
int spacedim>
3221 inline std::vector<IndexSet>
3222 extract_locally_owned_index_sets(
3224 const unsigned int level)
3226 std::vector<IndexSet> locally_owned_set;
3227 locally_owned_set.reserve(dofh.size());
3228 for (
unsigned int j = 0; j < dofh.size(); ++j)
3230 locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
3232 locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(
level));
3233 return locally_owned_set;
3240template <
int dim,
typename Number,
typename VectorizedArrayType>
3241template <
typename QuadratureType,
typename number2,
typename MappingType>
3244 const MappingType &mapping,
3247 const QuadratureType &quad,
3251 std::vector<const DoFHandler<dim, dim> *> dof_handlers;
3252 std::vector<const AffineConstraints<number2> *> constraints;
3254 dof_handlers.push_back(&dof_handler);
3255 constraints.push_back(&constraints_in);
3257 std::vector<IndexSet> locally_owned_sets =
3258 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3259 dof_handlers, additional_data.
mg_level);
3261 std::vector<hp::QCollection<dim>> quad_hp;
3262 quad_hp.emplace_back(quad);
3274template <
int dim,
typename Number,
typename VectorizedArrayType>
3275template <
typename QuadratureType,
typename number2,
typename MappingType>
3278 const MappingType &mapping,
3281 const QuadratureType &quad,
3285 std::vector<IndexSet> locally_owned_set =
3286 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3287 dof_handler, additional_data.
mg_level);
3288 std::vector<hp::QCollection<dim>> quad_hp;
3289 quad_hp.emplace_back(quad);
3301template <
int dim,
typename Number,
typename VectorizedArrayType>
3302template <
typename QuadratureType,
typename number2,
typename MappingType>
3305 const MappingType &mapping,
3308 const std::vector<QuadratureType> &quad,
3312 std::vector<IndexSet> locally_owned_set =
3313 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3314 dof_handler, additional_data.
mg_level);
3315 std::vector<hp::QCollection<dim>> quad_hp;
3316 for (
unsigned int q = 0; q < quad.size(); ++q)
3317 quad_hp.emplace_back(quad[q]);
3344 template <
int dim,
typename Number,
typename VectorizedArrayType>
3345 struct VectorDataExchange
3352 static constexpr unsigned int channel_shift = 20;
3361 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3362 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3363 DataAccessOnFaces vector_face_access,
3365 : matrix_free(matrix_free)
3366 , vector_face_access(
3367 matrix_free.get_task_info().face_partition_data.empty() ?
3369 DataAccessOnFaces::unspecified :
3371 , ghosts_were_set(false)
3372# ifdef DEAL_II_WITH_MPI
3373 , tmp_data(n_components)
3374 , requests(n_components)
3378 if (this->vector_face_access !=
3380 DataAccessOnFaces::unspecified)
3381 for (unsigned
int c = 0; c < matrix_free.n_components(); ++c)
3383 matrix_free.get_dof_info(c).vector_exchanger_face_variants.size(),
3392 ~VectorDataExchange()
3394# ifdef DEAL_II_WITH_MPI
3395 for (
unsigned int i = 0; i < tmp_data.size(); ++i)
3396 if (tmp_data[i] !=
nullptr)
3397 matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3407 template <
typename VectorType>
3409 find_vector_in_mf(
const VectorType &vec,
3410 const bool check_global_compatibility =
true)
const
3413 if (matrix_free.n_components() == 1)
3415 Assert(matrix_free.get_dof_info(0).vector_partitioner->is_compatible(
3416 *vec.get_partitioner()),
3417 ExcMessage(
"Could not find partitioner that fits vector"));
3422 for (
unsigned int c = 0; c < matrix_free.n_components(); ++c)
3423 if (vec.get_partitioner().get() ==
3424 matrix_free.get_dof_info(c).vector_partitioner.get())
3428 for (
unsigned int c = 0; c < matrix_free.n_components(); ++c)
3429 if (check_global_compatibility ?
3430 vec.get_partitioner()->is_globally_compatible(
3431 *matrix_free.get_dof_info(c).vector_partitioner) :
3432 vec.get_partitioner()->is_compatible(
3433 *matrix_free.get_dof_info(c).vector_partitioner))
3449 get_partitioner(
const unsigned int mf_component)
const
3452 .vector_exchanger_face_variants.size(),
3454 if (vector_face_access ==
3456 DataAccessOnFaces::
none)
3457 return *matrix_free.get_dof_info(mf_component)
3458 .vector_exchanger_face_variants[0];
3459 else if (vector_face_access ==
3461 DataAccessOnFaces::
values)
3462 return *matrix_free.get_dof_info(mf_component)
3463 .vector_exchanger_face_variants[1];
3464 else if (vector_face_access ==
3467 return *matrix_free.get_dof_info(mf_component)
3468 .vector_exchanger_face_variants[2];
3469 else if (vector_face_access ==
3471 DataAccessOnFaces::values_all_faces)
3472 return *matrix_free.get_dof_info(mf_component)
3473 .vector_exchanger_face_variants[3];
3474 else if (vector_face_access ==
3476 DataAccessOnFaces::gradients_all_faces)
3477 return *matrix_free.get_dof_info(mf_component)
3478 .vector_exchanger_face_variants[4];
3480 return *matrix_free.get_dof_info(mf_component).vector_exchanger.get();
3489 std::enable_if_t<is_not_parallel_vector<VectorType>,
VectorType>
3492 update_ghost_values_start(
const unsigned int ,
3493 const VectorType & )
3502 std::enable_if_t<!has_update_ghost_values_start<VectorType> &&
3503 !is_not_parallel_vector<VectorType>,
3506 update_ghost_values_start(
const unsigned int component_in_block_vector,
3507 const VectorType &vec)
3509 (void)component_in_block_vector;
3510 const bool ghosts_set = vec.has_ghost_elements();
3512 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3513 ghosts_set ==
false,
3518 ghosts_were_set =
true;
3522 vec.update_ghost_values();
3533 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3534 !has_exchange_on_subset<VectorType>,
3537 update_ghost_values_start(
const unsigned int component_in_block_vector,
3538 const VectorType &vec)
3540 (void)component_in_block_vector;
3541 const bool ghosts_set = vec.has_ghost_elements();
3543 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3544 ghosts_set ==
false,
3549 ghosts_were_set =
true;
3553 vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3565 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3566 has_exchange_on_subset<VectorType>,
3569 update_ghost_values_start(
const unsigned int component_in_block_vector,
3570 const VectorType &vec)
3572 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3573 "Type mismatch between VectorType and VectorDataExchange");
3574 (void)component_in_block_vector;
3575 const bool ghosts_set = vec.has_ghost_elements();
3577 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3578 ghosts_set ==
false,
3583 ghosts_were_set =
true;
3587 if (vec.size() != 0)
3589# ifdef DEAL_II_WITH_MPI
3590 const unsigned int mf_component = find_vector_in_mf(vec);
3592 const auto &part = get_partitioner(mf_component);
3594 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3595 part.n_import_sm_procs() == 0)
3598 tmp_data[component_in_block_vector] =
3599 matrix_free.acquire_scratch_data_non_threadsafe();
3600 tmp_data[component_in_block_vector]->resize_fast(
3601 part.n_import_indices());
3604 part.export_to_ghosted_array_start(
3605 component_in_block_vector * 2 + channel_shift,
3607 vec.shared_vector_data(),
3609 part.locally_owned_size(),
3610 matrix_free.get_dof_info(mf_component)
3611 .vector_partitioner->n_ghost_indices()),
3613 part.n_import_indices()),
3614 this->requests[component_in_block_vector]);
3626 std::enable_if_t<!has_update_ghost_values_start<VectorType>,
3629 update_ghost_values_finish(
const unsigned int ,
3630 const VectorType & )
3641 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3642 !has_exchange_on_subset<VectorType>,
3645 update_ghost_values_finish(
const unsigned int component_in_block_vector,
3646 const VectorType &vec)
3648 (void)component_in_block_vector;
3650 if (ghosts_were_set)
3653 vec.update_ghost_values_finish();
3665 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3666 has_exchange_on_subset<VectorType>,
3669 update_ghost_values_finish(
const unsigned int component_in_block_vector,
3670 const VectorType &vec)
3672 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3673 "Type mismatch between VectorType and VectorDataExchange");
3674 (void)component_in_block_vector;
3676 if (ghosts_were_set)
3679 if (vec.size() != 0)
3681# ifdef DEAL_II_WITH_MPI
3685 const unsigned int mf_component = find_vector_in_mf(vec);
3687 const auto &part = get_partitioner(mf_component);
3689 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3690 part.n_import_sm_procs() != 0)
3692 part.export_to_ghosted_array_finish(
3694 vec.shared_vector_data(),
3696 part.locally_owned_size(),
3697 matrix_free.get_dof_info(mf_component)
3698 .vector_partitioner->n_ghost_indices()),
3699 this->requests[component_in_block_vector]);
3701 matrix_free.release_scratch_data_non_threadsafe(
3702 tmp_data[component_in_block_vector]);
3703 tmp_data[component_in_block_vector] =
nullptr;
3709 vec.set_ghost_state(
true);
3718 std::enable_if_t<is_not_parallel_vector<VectorType>,
VectorType>
3721 compress_start(
const unsigned int ,
3732 std::enable_if_t<!has_compress_start<VectorType> &&
3733 !is_not_parallel_vector<VectorType>,
3736 compress_start(
const unsigned int component_in_block_vector,
3739 (void)component_in_block_vector;
3752 std::enable_if_t<has_compress_start<VectorType> &&
3753 !has_exchange_on_subset<VectorType>,
3756 compress_start(
const unsigned int component_in_block_vector,
3759 (void)component_in_block_vector;
3761 vec.compress_start(component_in_block_vector + channel_shift);
3773 std::enable_if_t<has_compress_start<VectorType> &&
3774 has_exchange_on_subset<VectorType>,
3777 compress_start(
const unsigned int component_in_block_vector,
3780 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3781 "Type mismatch between VectorType and VectorDataExchange");
3782 (void)component_in_block_vector;
3785 if (vec.size() != 0)
3787# ifdef DEAL_II_WITH_MPI
3788 const unsigned int mf_component = find_vector_in_mf(vec);
3790 const auto &part = get_partitioner(mf_component);
3792 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3793 part.n_import_sm_procs() == 0)
3796 tmp_data[component_in_block_vector] =
3797 matrix_free.acquire_scratch_data_non_threadsafe();
3798 tmp_data[component_in_block_vector]->resize_fast(
3799 part.n_import_indices());
3802 part.import_from_ghosted_array_start(
3804 component_in_block_vector * 2 + channel_shift,
3806 vec.shared_vector_data(),
3808 matrix_free.get_dof_info(mf_component)
3809 .vector_partitioner->n_ghost_indices()),
3811 part.n_import_indices()),
3812 this->requests[component_in_block_vector]);
3825 std::enable_if_t<!has_compress_start<VectorType>,
VectorType> * =
nullptr>
3827 compress_finish(
const unsigned int ,
3839 std::enable_if_t<has_compress_start<VectorType> &&
3840 !has_exchange_on_subset<VectorType>,
3843 compress_finish(
const unsigned int component_in_block_vector,
3846 (void)component_in_block_vector;
3859 std::enable_if_t<has_compress_start<VectorType> &&
3860 has_exchange_on_subset<VectorType>,
3863 compress_finish(
const unsigned int component_in_block_vector,
3866 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3867 "Type mismatch between VectorType and VectorDataExchange");
3868 (void)component_in_block_vector;
3869 if (vec.size() != 0)
3871# ifdef DEAL_II_WITH_MPI
3875 const unsigned int mf_component = find_vector_in_mf(vec);
3877 const auto &part = get_partitioner(mf_component);
3879 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3880 part.n_import_sm_procs() != 0)
3882 part.import_from_ghosted_array_finish(
3885 vec.shared_vector_data(),
3887 matrix_free.get_dof_info(mf_component)
3888 .vector_partitioner->n_ghost_indices()),
3890 tmp_data[component_in_block_vector]->begin(),
3891 part.n_import_indices()),
3892 this->requests[component_in_block_vector]);
3894 matrix_free.release_scratch_data_non_threadsafe(
3895 tmp_data[component_in_block_vector]);
3896 tmp_data[component_in_block_vector] =
nullptr;
3902 MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3915 std::enable_if_t<is_not_parallel_vector<VectorType>,
VectorType>
3918 reset_ghost_values(
const VectorType & )
const
3928 std::enable_if_t<!has_exchange_on_subset<VectorType> &&
3929 !is_not_parallel_vector<VectorType>,
3932 reset_ghost_values(
const VectorType &vec)
const
3934 if (ghosts_were_set ==
true)
3937 vec.zero_out_ghost_values();
3948 std::enable_if_t<has_exchange_on_subset<VectorType>,
VectorType>
3951 reset_ghost_values(
const VectorType &vec)
const
3953 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3954 "Type mismatch between VectorType and VectorDataExchange");
3955 if (ghosts_were_set ==
true)
3958 if (vec.size() != 0)
3960# ifdef DEAL_II_WITH_MPI
3963 const unsigned int mf_component = find_vector_in_mf(vec);
3965 const auto &part = get_partitioner(mf_component);
3967 if (part.n_ghost_indices() > 0)
3969 part.reset_ghost_values(
3971 part.locally_owned_size(),
3972 matrix_free.get_dof_info(mf_component)
3973 .vector_partitioner->n_ghost_indices()));
3979 vec.set_ghost_state(
false);
3989 template <
typename VectorType>
3991 zero_vector_region(
const unsigned int range_index, VectorType &vec)
const
3993 if constexpr (has_exchange_on_subset<VectorType>)
3996 std::is_same_v<Number, typename VectorType::value_type>,
3997 "Type mismatch between VectorType and VectorDataExchange");
4002 const unsigned int mf_component = find_vector_in_mf(vec,
false);
4004 matrix_free.get_dof_info(mf_component);
4008 Assert(vec.partitioners_are_compatible(
4014 for (
unsigned int id =
4018 std::memset(vec.begin() +
4026 else if constexpr (has_assignment_operator<VectorType>)
4030 if constexpr (std::is_same_v<
4034 for (
unsigned int i = 0; i < vec.size(); ++i)
4035 vec[i] =
typename VectorType::value_type();
4038 vec =
typename VectorType::value_type();
4045 "Zeroing is only implemented for vector types "
4046 "which provide operator=(const VectorType::value_type)"));
4052 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
4053 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
4054 DataAccessOnFaces vector_face_access;
4055 bool ghosts_were_set;
4056# ifdef DEAL_II_WITH_MPI
4059 boost::container::small_vector<AlignedVector<Number> *, 16> tmp_data;
4060 boost::container::small_vector<std::vector<MPI_Request>, 16> requests;
4064 template <
typename VectorStruct>
4066 n_components(
const VectorStruct &vec);
4068 template <
typename VectorStruct>
4070 n_components_block(
const VectorStruct &vec,
const std::bool_constant<true>)
4072 unsigned int components = 0;
4073 for (
unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
4074 components += n_components(vec.block(bl));
4078 template <
typename VectorStruct>
4080 n_components_block(
const VectorStruct &,
const std::bool_constant<false>)
4085 template <
typename VectorStruct>
4087 n_components(
const VectorStruct &vec)
4089 return n_components_block(
4093 template <
typename VectorStruct>
4095 n_components(
const std::vector<VectorStruct> &vec)
4097 unsigned int components = 0;
4098 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4099 components += n_components_block(
4104 template <
typename VectorStruct>
4106 n_components(
const std::vector<VectorStruct *> &vec)
4108 unsigned int components = 0;
4109 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4110 components += n_components_block(
4122 template <
typename VectorStruct,
4123 std::enable_if_t<!has_communication_block_size<VectorStruct>,
4124 VectorStruct> * =
nullptr>
4125 constexpr unsigned int
4126 get_communication_block_size(
const VectorStruct &)
4133 template <
typename VectorStruct,
4134 std::enable_if_t<has_communication_block_size<VectorStruct>,
4135 VectorStruct> * =
nullptr>
4136 constexpr unsigned int
4137 get_communication_block_size(
const VectorStruct &)
4139 return VectorStruct::communication_block_size;
4145 std::enable_if_t<is_not_parallel_vector<VectorType>,
VectorType> * =
4148 has_ghost_elements(
const VectorType &vec)
4157 std::enable_if_t<!is_not_parallel_vector<VectorType>,
VectorType>
4160 has_ghost_elements(
const VectorType &vec)
4162 return vec.has_ghost_elements();
4177 typename VectorStruct,
4179 typename VectorizedArrayType,
4180 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4183 update_ghost_values_start(
4184 const VectorStruct &vec,
4185 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4186 const unsigned int channel = 0)
4188 if (get_communication_block_size(vec) < vec.n_blocks())
4190 const bool ghosts_set = vec.has_ghost_elements();
4192 Assert(exchanger.matrix_free.get_task_info()
4193 .allow_ghosted_vectors_in_loops ||
4194 ghosts_set ==
false,
4199 exchanger.ghosts_were_set =
true;
4203 vec.update_ghost_values();
4207 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4208 update_ghost_values_start(vec.block(i), exchanger, channel + i);
4216 typename VectorStruct,
4218 typename VectorizedArrayType,
4219 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4222 update_ghost_values_start(
4223 const VectorStruct &vec,
4224 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4225 const unsigned int channel = 0)
4227 exchanger.update_ghost_values_start(channel, vec);
4234 typename VectorStruct,
4236 typename VectorizedArrayType>
4238 update_ghost_values_start(
4239 const std::vector<VectorStruct> &vec,
4240 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4242 unsigned int component_index = 0;
4243 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4245 update_ghost_values_start(vec[comp], exchanger, component_index);
4246 component_index += n_components(vec[comp]);
4254 typename VectorStruct,
4256 typename VectorizedArrayType>
4258 update_ghost_values_start(
4259 const std::vector<VectorStruct *> &vec,
4260 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4262 unsigned int component_index = 0;
4263 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4265 update_ghost_values_start(*vec[comp], exchanger, component_index);
4266 component_index += n_components(*vec[comp]);
4278 typename VectorStruct,
4280 typename VectorizedArrayType,
4281 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4284 update_ghost_values_finish(
4285 const VectorStruct &vec,
4286 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4287 const unsigned int channel = 0)
4289 if (get_communication_block_size(vec) < vec.n_blocks())
4295 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4296 update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4303 typename VectorStruct,
4305 typename VectorizedArrayType,
4306 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4309 update_ghost_values_finish(
4310 const VectorStruct &vec,
4311 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4312 const unsigned int channel = 0)
4314 exchanger.update_ghost_values_finish(channel, vec);
4321 typename VectorStruct,
4323 typename VectorizedArrayType>
4325 update_ghost_values_finish(
4326 const std::vector<VectorStruct> &vec,
4327 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4329 unsigned int component_index = 0;
4330 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4332 update_ghost_values_finish(vec[comp], exchanger, component_index);
4333 component_index += n_components(vec[comp]);
4341 typename VectorStruct,
4343 typename VectorizedArrayType>
4345 update_ghost_values_finish(
4346 const std::vector<VectorStruct *> &vec,
4347 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4349 unsigned int component_index = 0;
4350 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4352 update_ghost_values_finish(*vec[comp], exchanger, component_index);
4353 component_index += n_components(*vec[comp]);
4365 typename VectorStruct,
4367 typename VectorizedArrayType,
4368 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4373 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4374 const unsigned int channel = 0)
4376 if (get_communication_block_size(vec) < vec.n_blocks())
4379 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4380 compress_start(vec.block(i), exchanger, channel + i);
4387 typename VectorStruct,
4389 typename VectorizedArrayType,
4390 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4395 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4396 const unsigned int channel = 0)
4398 exchanger.compress_start(channel, vec);
4405 typename VectorStruct,
4407 typename VectorizedArrayType>
4410 std::vector<VectorStruct> &vec,
4411 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4413 unsigned int component_index = 0;
4414 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4416 compress_start(vec[comp], exchanger, component_index);
4417 component_index += n_components(vec[comp]);
4425 typename VectorStruct,
4427 typename VectorizedArrayType>
4430 std::vector<VectorStruct *> &vec,
4431 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4433 unsigned int component_index = 0;
4434 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4436 compress_start(*vec[comp], exchanger, component_index);
4437 component_index += n_components(*vec[comp]);
4449 typename VectorStruct,
4451 typename VectorizedArrayType,
4452 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4457 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4458 const unsigned int channel = 0)
4460 if (get_communication_block_size(vec) < vec.n_blocks())
4466 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4467 compress_finish(vec.block(i), exchanger, channel + i);
4474 typename VectorStruct,
4476 typename VectorizedArrayType,
4477 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4482 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4483 const unsigned int channel = 0)
4485 exchanger.compress_finish(channel, vec);
4492 typename VectorStruct,
4494 typename VectorizedArrayType>
4497 std::vector<VectorStruct> &vec,
4498 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4500 unsigned int component_index = 0;
4501 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4503 compress_finish(vec[comp], exchanger, component_index);
4504 component_index += n_components(vec[comp]);
4512 typename VectorStruct,
4514 typename VectorizedArrayType>
4517 std::vector<VectorStruct *> &vec,
4518 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4520 unsigned int component_index = 0;
4521 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4523 compress_finish(*vec[comp], exchanger, component_index);
4524 component_index += n_components(*vec[comp]);
4540 typename VectorStruct,
4542 typename VectorizedArrayType,
4543 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4547 const VectorStruct &vec,
4548 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4551 if (exchanger.ghosts_were_set ==
true)
4554 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4555 reset_ghost_values(vec.block(i), exchanger);
4562 typename VectorStruct,
4564 typename VectorizedArrayType,
4565 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4569 const VectorStruct &vec,
4570 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4573 if (exchanger.ghosts_were_set ==
true)
4576 exchanger.reset_ghost_values(vec);
4583 typename VectorStruct,
4585 typename VectorizedArrayType>
4588 const std::vector<VectorStruct> &vec,
4589 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4592 if (exchanger.ghosts_were_set ==
true)
4595 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4596 reset_ghost_values(vec[comp], exchanger);
4603 typename VectorStruct,
4605 typename VectorizedArrayType>
4608 const std::vector<VectorStruct *> &vec,
4609 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4612 if (exchanger.ghosts_were_set ==
true)
4615 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4616 reset_ghost_values(*vec[comp], exchanger);
4627 typename VectorStruct,
4629 typename VectorizedArrayType,
4630 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4634 const unsigned int range_index,
4636 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4638 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4639 exchanger.zero_vector_region(range_index, vec.block(i));
4646 typename VectorStruct,
4648 typename VectorizedArrayType,
4649 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4653 const unsigned int range_index,
4655 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4657 exchanger.zero_vector_region(range_index, vec);
4664 typename VectorStruct,
4666 typename VectorizedArrayType>
4669 const unsigned int range_index,
4670 std::vector<VectorStruct> &vec,
4671 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4673 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4674 zero_vector_region(range_index, vec[comp], exchanger);
4681 typename VectorStruct,
4683 typename VectorizedArrayType>
4686 const unsigned int range_index,
4687 std::vector<VectorStruct *> &vec,
4688 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4690 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4691 zero_vector_region(range_index, *vec[comp], exchanger);
4699 template <
typename VectorStruct1,
typename VectorStruct2>
4701 apply_operation_to_constrained_dofs(
const std::vector<unsigned int> &,
4702 const VectorStruct1 &,
4706 template <
typename Number>
4708 apply_operation_to_constrained_dofs(
4709 const std::vector<unsigned int> &constrained_dofs,
4713 for (
const unsigned int i : constrained_dofs)
4714 dst.local_element(i) = src.local_element(i);
4718 namespace MatrixFreeFunctions
4722 template <
typename,
typename,
typename,
typename,
bool>
4723 struct InterfaceSelector
4727 template <
typename MF,
4731 struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4733 using function_type = void (Container::*)(
4741 template <
typename MF,
4745 struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4747 using function_type =
4748 void (Container::*)(
const MF &,
4759 template <
typename MF,
4764 class MFWorker :
public MFWorkerInterface
4768 using function_type =
typename MatrixFreeFunctions::
4769 InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4773 MFWorker(
const MF &matrix_free,
4774 const InVector &src,
4776 const bool zero_dst_vector_setting,
4777 const Container &container,
4778 function_type cell_function,
4779 function_type face_function,
4780 function_type boundary_function,
4781 const typename MF::DataAccessOnFaces src_vector_face_access =
4782 MF::DataAccessOnFaces::none,
4783 const typename MF::DataAccessOnFaces dst_vector_face_access =
4784 MF::DataAccessOnFaces::none,
4785 const std::function<
void(
const unsigned int,
const unsigned int)>
4786 &operation_before_loop = {},
4787 const std::function<void(
const unsigned int,
const unsigned int)>
4788 &operation_after_loop = {},
4789 const unsigned int dof_handler_index_pre_post = 0)
4790 : matrix_free(matrix_free)
4791 , container(const_cast<Container &>(container))
4792 , cell_function(cell_function)
4793 , face_function(face_function)
4794 , boundary_function(boundary_function)
4797 , src_data_exchanger(matrix_free,
4798 src_vector_face_access,
4800 , dst_data_exchanger(matrix_free,
4801 dst_vector_face_access,
4804 , zero_dst_vector_setting(zero_dst_vector_setting &&
4805 !src_and_dst_are_same)
4806 , operation_before_loop(operation_before_loop)
4807 , operation_after_loop(operation_after_loop)
4808 , dof_handler_index_pre_post(dof_handler_index_pre_post)
4810 Assert(!has_ghost_elements(dst),
4811 ExcMessage(
"The destination vector passed to the matrix-free "
4812 "loop is ghosted. This is not allowed."));
4817 cell(
const std::pair<unsigned int, unsigned int> &cell_range)
override
4819 if (cell_function !=
nullptr && cell_range.second > cell_range.first)
4820 for (
unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4822 const auto cell_subrange =
4823 matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4825 if (cell_subrange.second <= cell_subrange.first)
4829 cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4834 cell(
const unsigned int range_index)
override
4836 process_range(cell_function,
4837 matrix_free.get_task_info().cell_partition_data_hp_ptr,
4838 matrix_free.get_task_info().cell_partition_data_hp,
4843 face(
const unsigned int range_index)
override
4845 process_range(face_function,
4846 matrix_free.get_task_info().face_partition_data_hp_ptr,
4847 matrix_free.get_task_info().face_partition_data_hp,
4852 boundary(
const unsigned int range_index)
override
4854 process_range(boundary_function,
4855 matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4856 matrix_free.get_task_info().boundary_partition_data_hp,
4862 process_range(
const function_type &fu,
4863 const std::vector<unsigned int> &ptr,
4864 const std::vector<unsigned int> &
data,
4865 const unsigned int range_index)
4871 for (
unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4874 (container.*fu)(matrix_free,
4877 std::make_pair(
data[2 * i],
data[2 * i + 1]));
4889 vector_update_ghosts_start()
override
4891 if (!src_and_dst_are_same)
4892 internal::update_ghost_values_start(src, src_data_exchanger);
4897 vector_update_ghosts_finish()
override
4899 if (!src_and_dst_are_same)
4900 internal::update_ghost_values_finish(src, src_data_exchanger);
4905 vector_compress_start()
override
4907 internal::compress_start(dst, dst_data_exchanger);
4912 vector_compress_finish()
override
4914 internal::compress_finish(dst, dst_data_exchanger);
4915 if (!src_and_dst_are_same)
4916 internal::reset_ghost_values(src, src_data_exchanger);
4921 zero_dst_vector_range(
const unsigned int range_index)
override
4923 if (zero_dst_vector_setting)
4924 internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4928 cell_loop_pre_range(
const unsigned int range_index)
override
4930 if (operation_before_loop)
4933 matrix_free.get_dof_info(dof_handler_index_pre_post);
4940 operation_before_loop,
4947 for (
unsigned int id =
4958 cell_loop_post_range(
const unsigned int range_index)
override
4960 if (operation_after_loop)
4964 const std::vector<unsigned int> &partition_row_index =
4965 matrix_free.get_task_info().partition_row_index;
4967 partition_row_index[partition_row_index.size() - 2] - 1)
4968 apply_operation_to_constrained_dofs(
4969 matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
4974 matrix_free.get_dof_info(dof_handler_index_pre_post);
4981 operation_after_loop,
4988 for (
unsigned int id =
4999 const MF &matrix_free;
5000 Container &container;
5001 function_type cell_function;
5002 function_type face_function;
5003 function_type boundary_function;
5005 const InVector &src;
5007 VectorDataExchange<MF::dimension,
5008 typename MF::value_type,
5009 typename MF::vectorized_value_type>
5011 VectorDataExchange<MF::dimension,
5012 typename MF::value_type,
5013 typename MF::vectorized_value_type>
5015 const bool src_and_dst_are_same;
5016 const bool zero_dst_vector_setting;
5017 const std::function<void(
const unsigned int,
const unsigned int)>
5018 operation_before_loop;
5019 const std::function<void(
const unsigned int,
const unsigned int)>
5020 operation_after_loop;
5021 const unsigned int dof_handler_index_pre_post;
5030 template <
class MF,
typename InVector,
typename OutVector>
5031 struct MFClassWrapper
5033 using function_type =
5034 std::function<void(
const MF &,
5037 const std::pair<unsigned int, unsigned int> &)>;
5039 MFClassWrapper(
const function_type cell,
5040 const function_type face,
5041 const function_type boundary)
5044 , boundary(boundary)
5048 cell_integrator(
const MF &mf,
5050 const InVector &src,
5051 const std::pair<unsigned int, unsigned int> &range)
const
5054 cell(mf, dst, src, range);
5058 face_integrator(
const MF &mf,
5060 const InVector &src,
5061 const std::pair<unsigned int, unsigned int> &range)
const
5064 face(mf, dst, src, range);
5068 boundary_integrator(
5071 const InVector &src,
5072 const std::pair<unsigned int, unsigned int> &range)
const
5075 boundary(mf, dst, src, range);
5078 const function_type cell;
5079 const function_type face;
5080 const function_type boundary;
5087template <
int dim,
typename Number,
typename VectorizedArrayType>
5088template <
typename OutVector,
typename InVector>
5094 const std::pair<unsigned int, unsigned int> &)>
5097 const InVector &src,
5098 const bool zero_dst_vector)
const
5101 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5104 Wrapper wrap(cell_operation,
nullptr,
nullptr);
5105 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5115 &Wrapper::cell_integrator,
5116 &Wrapper::face_integrator,
5117 &Wrapper::boundary_integrator);
5119 task_info.loop(worker);
5124template <
int dim,
typename Number,
typename VectorizedArrayType>
5125template <
typename OutVector,
typename InVector>
5131 const std::pair<unsigned int, unsigned int> &)>
5134 const InVector &src,
5135 const std::function<
void(
const unsigned int,
const unsigned int)>
5136 &operation_before_loop,
5137 const std::function<
void(
const unsigned int,
const unsigned int)>
5138 &operation_after_loop,
5139 const unsigned int dof_handler_index_pre_post)
const
5142 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5145 Wrapper wrap(cell_operation,
nullptr,
nullptr);
5146 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5156 &Wrapper::cell_integrator,
5157 &Wrapper::face_integrator,
5158 &Wrapper::boundary_integrator,
5159 DataAccessOnFaces::none,
5160 DataAccessOnFaces::none,
5161 operation_before_loop,
5162 operation_after_loop,
5163 dof_handler_index_pre_post);
5165 task_info.loop(worker);
5170template <
int dim,
typename Number,
typename VectorizedArrayType>
5171template <
typename OutVector,
typename InVector>
5177 const std::pair<unsigned int, unsigned int> &)>
5182 const std::pair<unsigned int, unsigned int> &)>
5183 &inner_face_operation,
5187 const std::pair<unsigned int, unsigned int> &)>
5188 &boundary_face_operation,
5190 const InVector &src,
5191 const bool zero_dst_vector,
5192 const DataAccessOnFaces dst_vector_face_access,
5193 const DataAccessOnFaces src_vector_face_access)
const
5196 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5199 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5200 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5210 &Wrapper::cell_integrator,
5211 &Wrapper::face_integrator,
5212 &Wrapper::boundary_integrator,
5213 src_vector_face_access,
5214 dst_vector_face_access);
5216 task_info.loop(worker);
5221template <
int dim,
typename Number,
typename VectorizedArrayType>
5222template <
typename CLASS,
typename OutVector,
typename InVector>
5225 void (CLASS::*function_pointer)(
5230 const CLASS *owning_class,
5232 const InVector &src,
5235 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5253template <
int dim,
typename Number,
typename VectorizedArrayType>
5254template <
typename CLASS,
typename OutVector,
typename InVector>
5257 void (CLASS::*function_pointer)(
5262 const CLASS *owning_class,
5264 const InVector &src,
5266 &operation_before_loop,
5268 &operation_after_loop,
5271 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5286 operation_before_loop,
5287 operation_after_loop,
5288 dof_handler_index_pre_post);
5294template <
int dim,
typename Number,
typename VectorizedArrayType>
5295template <
typename CLASS,
typename OutVector,
typename InVector>
5298 void (CLASS::*cell_operation)(
5303 void (CLASS::*inner_face_operation)(
5308 void (CLASS::*boundary_face_operation)(
5313 const CLASS *owning_class,
5315 const InVector &src,
5317 const DataAccessOnFaces dst_vector_face_access,
5318 const DataAccessOnFaces src_vector_face_access)
const
5320 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5331 inner_face_operation,
5332 boundary_face_operation,
5333 src_vector_face_access,
5334 dst_vector_face_access);
5340template <
int dim,
typename Number,
typename VectorizedArrayType>
5341template <
typename CLASS,
typename OutVector,
typename InVector>
5344 void (CLASS::*function_pointer)(
5349 CLASS *owning_class,
5351 const InVector &src,
5354 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5372template <
int dim,
typename Number,
typename VectorizedArrayType>
5373template <
typename CLASS,
typename OutVector,
typename InVector>
5376 void (CLASS::*function_pointer)(
5381 CLASS *owning_class,
5383 const InVector &src,
5385 &operation_before_loop,
5387 &operation_after_loop,
5390 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5405 operation_before_loop,
5406 operation_after_loop,
5407 dof_handler_index_pre_post);
5413template <
int dim,
typename Number,
typename VectorizedArrayType>
5414template <
typename CLASS,
typename OutVector,
typename InVector>
5417 void (CLASS::*cell_operation)(
5422 void (CLASS::*inner_face_operation)(
5427 void (CLASS::*boundary_face_operation)(
5432 CLASS *owning_class,
5434 const InVector &src,
5436 const DataAccessOnFaces dst_vector_face_access,
5437 const DataAccessOnFaces src_vector_face_access)
const
5439 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5450 inner_face_operation,
5451 boundary_face_operation,
5452 src_vector_face_access,
5453 dst_vector_face_access);
5459template <
int dim,
typename Number,
typename VectorizedArrayType>
5460template <
typename OutVector,
typename InVector>
5466 const std::pair<unsigned int, unsigned int> &)>
5471 const std::pair<unsigned int, unsigned int> &)>
5472 &inner_face_operation,
5476 const std::pair<unsigned int, unsigned int> &)>
5477 &boundary_face_operation,
5479 const InVector &src,
5480 const std::function<
void(
const unsigned int,
const unsigned int)>
5481 &operation_before_loop,
5482 const std::function<
void(
const unsigned int,
const unsigned int)>
5483 &operation_after_loop,
5484 const unsigned int dof_handler_index_pre_post,
5485 const DataAccessOnFaces dst_vector_face_access,
5486 const DataAccessOnFaces src_vector_face_access)
const
5489 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5492 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5493 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5503 &Wrapper::cell_integrator,
5504 &Wrapper::face_integrator,
5505 &Wrapper::boundary_integrator,
5506 src_vector_face_access,
5507 dst_vector_face_access,
5508 operation_before_loop,
5509 operation_after_loop,
5510 dof_handler_index_pre_post);
5512 task_info.loop(worker);
5517template <
int dim,
typename Number,
typename VectorizedArrayType>
5518template <
typename CLASS,
typename OutVector,
typename InVector>
5526 void (CLASS::*inner_face_operation)(
5531 void (CLASS::*boundary_face_operation)(
5536 const CLASS *owning_class,
5538 const InVector &src,
5540 &operation_before_loop,
5542 &operation_after_loop,
5543 const unsigned
int dof_handler_index_pre_post,
5544 const DataAccessOnFaces dst_vector_face_access,
5545 const DataAccessOnFaces src_vector_face_access)
const
5547 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5558 inner_face_operation,
5559 boundary_face_operation,
5560 src_vector_face_access,
5561 dst_vector_face_access,
5562 operation_before_loop,
5563 operation_after_loop,
5564 dof_handler_index_pre_post);
5570template <
int dim,
typename Number,
typename VectorizedArrayType>
5571template <
typename CLASS,
typename OutVector,
typename InVector>
5578 void (CLASS::*inner_face_operation)(
5583 void (CLASS::*boundary_face_operation)(
5588 const CLASS *owning_class,
5590 const InVector &src,
5592 &operation_before_loop,
5594 &operation_after_loop,
5595 const unsigned
int dof_handler_index_pre_post,
5596 const DataAccessOnFaces dst_vector_face_access,
5597 const DataAccessOnFaces src_vector_face_access)
const
5599 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5610 inner_face_operation,
5611 boundary_face_operation,
5612 src_vector_face_access,
5613 dst_vector_face_access,
5614 operation_before_loop,
5615 operation_after_loop,
5616 dof_handler_index_pre_post);
5622template <
int dim,
typename Number,
typename VectorizedArrayType>
5623template <
typename CLASS,
typename OutVector,
typename InVector>
5626 void (CLASS::*function_pointer)(
5631 const CLASS *owning_class,
5633 const InVector &src,
5635 const DataAccessOnFaces src_vector_face_access)
const
5637 auto src_vector_face_access_temp = src_vector_face_access;
5643 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5656 src_vector_face_access_temp,
5663template <
int dim,
typename Number,
typename VectorizedArrayType>
5664template <
typename CLASS,
typename OutVector,
typename InVector>
5667 void (CLASS::*function_pointer)(
5672 CLASS *owning_class,
5674 const InVector &src,
5676 const DataAccessOnFaces src_vector_face_access)
const
5678 auto src_vector_face_access_temp = src_vector_face_access;
5684 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5697 src_vector_face_access_temp,
5704template <
int dim,
typename Number,
typename VectorizedArrayType>
5705template <
typename OutVector,
typename InVector>
5711 const std::pair<unsigned int, unsigned int> &)>
5714 const InVector &src,
5715 const bool zero_dst_vector,
5716 const DataAccessOnFaces src_vector_face_access)
const
5718 auto src_vector_face_access_temp = src_vector_face_access;
5719 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5720 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5721 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5722 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5725 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5728 Wrapper wrap(cell_operation,
nullptr,
nullptr);
5730 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5740 &Wrapper::cell_integrator,
5741 &Wrapper::face_integrator,
5742 &Wrapper::boundary_integrator,
5743 src_vector_face_access_temp,
5744 DataAccessOnFaces::none);
5745 task_info.loop(worker);
* * for(const auto &cell :triangulation.active_cell_iterators())
void resize(const size_type new_size)
BlockType & block(const unsigned int i)
void reinit(const size_type num_blocks, const size_type block_size=0, const bool omit_zeroing_entries=false)
void reinit(const size_type size, const bool omit_zeroing_entries=false)
Abstract base class for mapping classes.
unsigned int n_active_fe_indices() const
std::pair< typename DoFHandler< dim >::cell_iterator, unsigned int > get_face_iterator(const unsigned int face_batch_index, const unsigned int lane_index, const bool interior=true, const unsigned int fe_component=0) const
void loop_cell_centric(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > mapping_info
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
unsigned int n_ghost_cell_batches() const
void initialize_dof_handlers(const std::vector< const DoFHandler< dim, dim > * > &dof_handlers, const AdditionalData &additional_data)
types::boundary_id get_boundary_id(const unsigned int face_batch_index) const
void print(std::ostream &out) const
void update_mapping(const std::shared_ptr< hp::MappingCollection< dim > > &mapping)
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
unsigned int n_inner_face_batches() const
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
void update_mapping(const Mapping< dim > &mapping)
const Quadrature< dim > & get_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_index) const
unsigned int get_mg_level() const
unsigned int get_cell_category(const unsigned int cell_batch_index, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*boundary_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void print_memory_consumption(StreamType &out) const
std::pair< unsigned int, unsigned int > get_face_range_category(const std::pair< unsigned int, unsigned int > face_batch_range, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) const
bool mapping_initialized() const
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
void internal_reinit(const std::shared_ptr< hp::MappingCollection< dim > > &mapping, const std::vector< const DoFHandler< dim, dim > * > &dof_handlers, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const std::vector< hp::QCollection< q_dim > > &quad, const AdditionalData &additional_data)
unsigned int get_n_q_points_face(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &inner_face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_face_operation, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
bool at_irregular_cell(const unsigned int cell_batch_index) const
~MatrixFree() override=default
const AffineConstraints< Number > & get_affine_constraints(const unsigned int dof_handler_index=0) const
void copy_from(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free_base)
void initialize_face_data_vector(AlignedVector< T > &vec) const
unsigned int get_dofs_per_cell(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int n_constraint_pool_entries() const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
const Number * constraint_pool_begin(const unsigned int pool_index) const
AlignedVector< Number > * acquire_scratch_data_non_threadsafe() const
const IndexSet & get_locally_owned_set(const unsigned int dof_handler_index=0) const
unsigned int get_dofs_per_face(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
void initialize_dof_vector(LinearAlgebra::distributed::Vector< Number2, MemorySpace > &vec, const unsigned int dof_handler_index=0) const
internal::MatrixFreeFunctions::TaskInfo task_info
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
std::vector< ObserverPointer< const DoFHandler< dim > > > dof_handlers
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*boundary_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
unsigned int first_hp_dof_handler_index
bool mapping_is_initialized
void loop_cell_centric(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
DoFHandler< dim >::cell_iterator get_cell_iterator(const unsigned int cell_batch_index, const unsigned int lane_index, const unsigned int dof_handler_index=0) const
MatrixFree(const MatrixFree< dim, Number, VectorizedArrayType > &other)
std::pair< unsigned int, unsigned int > create_cell_subrange_hp(const std::pair< unsigned int, unsigned int > &range, const unsigned int fe_degree, const unsigned int dof_handler_index=0) const
std::pair< unsigned int, unsigned int > create_cell_subrange_hp_by_index(const std::pair< unsigned int, unsigned int > &range, const unsigned int fe_index, const unsigned int dof_handler_index=0) const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
std::vector< std::pair< unsigned int, unsigned int > > cell_level_index
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner(const unsigned int dof_handler_index=0) const
std::vector< ObserverPointer< const AffineConstraints< Number > > > affine_constraints
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
void initialize_dof_vector(LinearAlgebra::distributed::BlockVector< Number2, MemorySpace > &vec) const
std::vector< Number > constraint_pool_data
VectorizedArrayType vectorized_value_type
unsigned int n_cell_batches() const
const IndexSet & get_ghost_set(const unsigned int dof_handler_index=0) const
std::pair< int, int > get_cell_level_and_index(const unsigned int cell_batch_index, const unsigned int lane_index) const
bool indices_initialized() const
const Quadrature< dim - 1 > & get_face_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int get_face_active_fe_index(const std::pair< unsigned int, unsigned int > range, const bool is_interior_face=true, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) const
unsigned int get_cell_active_fe_index(const std::pair< unsigned int, unsigned int > range, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) const
Threads::ThreadLocalStorage< std::list< std::pair< bool, AlignedVector< VectorizedArrayType > > > > scratch_pad
unsigned int get_matrix_free_cell_index(const typename Triangulation< dim >::cell_iterator &cell) const
std::size_t memory_consumption() const
std::pair< unsigned int, unsigned int > get_face_category(const unsigned int face_batch_index, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) const
unsigned int n_boundary_face_batches() const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*boundary_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
std::list< std::pair< bool, AlignedVector< Number > > > scratch_pad_non_threadsafe
const Number * constraint_pool_end(const unsigned int pool_index) const
void loop_cell_centric(const std::function< void(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &inner_face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_face_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void initialize_indices(const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const AdditionalData &additional_data)
unsigned int n_components() const
unsigned int n_physical_cells() const
void initialize_dof_vector(VectorType &vec, const unsigned int dof_handler_index=0) const
unsigned int get_n_q_points(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
void reinit(const MappingType &mapping, const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
std::vector< unsigned int > mf_cell_indices
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
unsigned int get_cell_range_category(const std::pair< unsigned int, unsigned int > cell_batch_range, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*boundary_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
std::vector< internal::MatrixFreeFunctions::DoFInfo > dof_info
void initialize_cell_data_vector(AlignedVector< T > &vec) const
unsigned int n_ghost_inner_face_batches() const
void release_scratch_data_non_threadsafe(const AlignedVector< Number > *memory) const
void renumber_dofs(std::vector< types::global_dof_index > &renumbering, const unsigned int dof_handler_index=0)
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
bool indices_are_initialized
internal::MatrixFreeFunctions::FaceInfo< VectorizedArrayType::size()> face_info
void reinit(const MappingType &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< number2 > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
void reinit(const MappingType &mapping, const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< QuadratureType > &quad, const AdditionalData &additional_data=AdditionalData())
std::array< types::boundary_id, VectorizedArrayType::size()> get_faces_by_cells_boundary_id(const unsigned int cell_batch_index, const unsigned int face_number) const
std::vector< unsigned int > constraint_pool_row_index
const internal::MatrixFreeFunctions::ShapeInfo< Number > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
const std::vector< unsigned int > & get_constrained_dofs(const unsigned int dof_handler_index=0) const
unsigned int cell_level_index_end_local
unsigned int n_base_elements(const unsigned int dof_handler_index) const
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< Number > > shape_info
static bool is_supported(const FiniteElement< dim, spacedim > &fe)
static constexpr unsigned int dimension
A class that provides a separate storage location on each thread that accesses the object.
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
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
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_default
No update.
std::vector< index_type > data
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
* * * RotationFunction< dim, Number >::RotationFunction Number(dim)
* * * * std::vector< Number > ThermoPlasticMaterial< dim, ViscoplasticYieldLaw, Number >::get_state_parameters const
unsigned int minimum_parallel_grain_size
constexpr unsigned int invalid_unsigned_int
constexpr types::boundary_id invalid_boundary_id
void apply_to_subranges(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, const Function &f, const unsigned int grainsize)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
unsigned short int fe_index
TasksParallelScheme tasks_parallel_scheme
bool hold_all_faces_to_owned_cells
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
bool cell_vectorization_categories_strict
UpdateFlags mapping_update_flags_boundary_faces
UpdateFlags mapping_update_flags
UpdateFlags mapping_update_flags_faces_by_cells
AdditionalData & operator=(const AdditionalData &other)=default
AdditionalData(const AdditionalData &other)
bool allow_ghosted_vectors_in_loops
unsigned int tasks_block_size
bool overlap_communication_computation
@ dof_access_face_interior
std::vector< unsigned int > cell_loop_pre_list_index
std::vector< unsigned int > cell_loop_post_list_index
std::vector< std::pair< unsigned int, unsigned int > > vector_zero_range_list
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_pre_list
std::vector< unsigned int > vector_zero_range_list_index
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_post_list
void loop(MFWorkerInterface &worker) const