16#ifndef dealii_matrix_free_h
17#define dealii_matrix_free_h
110 typename Number = double,
115 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
116 "Type of Number and of VectorizedArrayType do not match.");
572 template <
typename QuadratureType,
typename number2,
typename MappingType>
577 const QuadratureType &quad,
601 template <
typename QuadratureType,
typename number2,
typename MappingType>
606 const std::vector<QuadratureType> &quad,
616 template <
typename QuadratureType,
typename number2,
typename MappingType>
621 const QuadratureType &quad,
778 template <
typename OutVector,
typename InVector>
784 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
787 const bool zero_dst_vector =
false)
const;
835 template <
typename CLASS,
typename OutVector,
typename InVector>
841 const
std::pair<unsigned
int, unsigned
int> &) const,
842 const CLASS *owning_class,
845 const
bool zero_dst_vector = false) const;
850 template <
typename CLASS,
typename OutVector,
typename InVector>
856 const
std::pair<unsigned
int, unsigned
int> &),
860 const
bool zero_dst_vector = false) const;
946 template <
typename CLASS,
typename OutVector,
typename InVector>
952 const
std::pair<unsigned
int, unsigned
int> &) const,
953 const CLASS *owning_class,
956 const
std::function<void(const unsigned
int, const unsigned
int)>
957 &operation_before_loop,
958 const
std::function<void(const unsigned
int, const unsigned
int)>
959 &operation_after_loop,
960 const unsigned
int dof_handler_index_pre_post = 0) const;
965 template <
typename CLASS,
typename OutVector,
typename InVector>
971 const
std::pair<unsigned
int, unsigned
int> &),
975 const
std::function<void(const unsigned
int, const unsigned
int)>
976 &operation_before_loop,
977 const
std::function<void(const unsigned
int, const unsigned
int)>
978 &operation_after_loop,
979 const unsigned
int dof_handler_index_pre_post = 0) const;
985 template <
typename OutVector,
typename InVector>
991 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
994 const std::function<
void(
const unsigned int,
const unsigned int)>
995 &operation_before_loop,
996 const std::function<
void(
const unsigned int,
const unsigned int)>
997 &operation_after_loop,
998 const unsigned int dof_handler_index_pre_post = 0)
const;
1075 template <
typename OutVector,
typename InVector>
1078 const std::function<
1082 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1083 const std::function<
void(
1087 const std::pair<unsigned int, unsigned int> &)> &inner_face_operation,
1088 const std::function<
void(
1092 const std::pair<unsigned int, unsigned int> &)> &boundary_face_operation,
1094 const InVector &src,
1095 const bool zero_dst_vector =
false,
1189 template <
typename CLASS,
typename OutVector,
typename InVector>
1191 loop(
void (CLASS::*cell_operation)(
1195 const
std::pair<unsigned
int, unsigned
int> &) const,
1196 void (CLASS::*inner_face_operation)(
1200 const
std::pair<unsigned
int, unsigned
int> &) const,
1201 void (CLASS::*boundary_face_operation)(
1205 const
std::pair<unsigned
int, unsigned
int> &) const,
1206 const CLASS *owning_class,
1208 const InVector &src,
1209 const
bool zero_dst_vector = false,
1218 template <
typename CLASS,
typename OutVector,
typename InVector>
1220 loop(
void (CLASS::*cell_operation)(
1224 const
std::pair<unsigned
int, unsigned
int> &),
1225 void (CLASS::*inner_face_operation)(
1229 const
std::pair<unsigned
int, unsigned
int> &),
1230 void (CLASS::*boundary_face_operation)(
1234 const
std::pair<unsigned
int, unsigned
int> &),
1235 CLASS *owning_class,
1237 const InVector &src,
1238 const
bool zero_dst_vector = false,
1365 template <
typename CLASS,
typename OutVector,
typename InVector>
1367 loop(
void (CLASS::*cell_operation)(
1371 const
std::pair<unsigned
int, unsigned
int> &) const,
1372 void (CLASS::*inner_face_operation)(
1376 const
std::pair<unsigned
int, unsigned
int> &) const,
1377 void (CLASS::*boundary_face_operation)(
1381 const
std::pair<unsigned
int, unsigned
int> &) const,
1382 const CLASS *owning_class,
1384 const InVector &src,
1385 const
std::function<void(const unsigned
int, const unsigned
int)>
1386 &operation_before_loop,
1387 const
std::function<void(const unsigned
int, const unsigned
int)>
1388 &operation_after_loop,
1389 const unsigned
int dof_handler_index_pre_post = 0,
1398 template <
typename CLASS,
typename OutVector,
typename InVector>
1400 loop(
void (CLASS::*cell_operation)(
1404 const
std::pair<unsigned
int, unsigned
int> &),
1405 void (CLASS::*inner_face_operation)(
1409 const
std::pair<unsigned
int, unsigned
int> &),
1410 void (CLASS::*boundary_face_operation)(
1414 const
std::pair<unsigned
int, unsigned
int> &),
1415 const CLASS *owning_class,
1417 const InVector &src,
1418 const
std::function<void(const unsigned
int, const unsigned
int)>
1419 &operation_before_loop,
1420 const
std::function<void(const unsigned
int, const unsigned
int)>
1421 &operation_after_loop,
1422 const unsigned
int dof_handler_index_pre_post = 0,
1433 template <
typename OutVector,
typename InVector>
1436 const std::function<
1440 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1441 const std::function<
void(
1445 const std::pair<unsigned int, unsigned int> &)> &inner_face_operation,
1446 const std::function<
void(
1450 const std::pair<unsigned int, unsigned int> &)> &boundary_face_operation,
1452 const InVector &src,
1453 const std::function<
void(
const unsigned int,
const unsigned int)>
1454 &operation_before_loop,
1455 const std::function<
void(
const unsigned int,
const unsigned int)>
1456 &operation_after_loop,
1457 const unsigned int dof_handler_index_pre_post = 0,
1527 template <
typename CLASS,
typename OutVector,
typename InVector>
1533 const
std::pair<unsigned
int, unsigned
int> &) const,
1534 const CLASS *owning_class,
1536 const InVector &src,
1537 const
bool zero_dst_vector = false,
1544 template <
typename CLASS,
typename OutVector,
typename InVector>
1550 const
std::pair<unsigned
int, unsigned
int> &),
1551 CLASS *owning_class,
1553 const InVector &src,
1554 const
bool zero_dst_vector = false,
1561 template <
typename OutVector,
typename InVector>
1567 const std::pair<unsigned int, unsigned int> &)>
1570 const InVector &src,
1571 const bool zero_dst_vector =
false,
1582 std::pair<unsigned int, unsigned int>
1584 const unsigned int fe_degree,
1585 const unsigned int dof_handler_index = 0)
const;
1593 std::pair<unsigned int, unsigned int>
1595 const std::pair<unsigned int, unsigned int> &range,
1596 const unsigned int fe_index,
1597 const unsigned int dof_handler_index = 0)
const;
1610 const std::pair<unsigned int, unsigned int> range)
const;
1617 const bool is_interior_face =
true)
const;
1630 template <
typename T>
1639 template <
typename T>
1656 template <
typename VectorType>
1659 const unsigned int dof_handler_index = 0)
const;
1678 template <
typename Number2>
1681 const unsigned int dof_handler_index = 0)
const;
1693 const std::shared_ptr<const Utilities::MPI::Partitioner> &
1717 const std::vector<unsigned int> &
1732 const unsigned int dof_handler_index = 0);
1743 template <
int spacedim>
1842 const unsigned int face_number)
const;
1874 const unsigned int lane_index,
1875 const unsigned int dof_handler_index = 0)
const;
1884 const unsigned int lane_index)
const;
1908 std::pair<typename DoFHandler<dim>::cell_iterator,
unsigned int>
1910 const unsigned int lane_index,
1911 const bool interior =
true,
1912 const unsigned int fe_component = 0)
const;
1958 const unsigned int hp_active_fe_index = 0)
const;
1965 const unsigned int hp_active_fe_index = 0)
const;
1973 const unsigned int hp_active_fe_index = 0)
const;
1981 const unsigned int hp_active_fe_index = 0)
const;
1988 const unsigned int hp_active_fe_index = 0)
const;
1995 const unsigned int hp_active_fe_index = 0)
const;
2011 const std::pair<unsigned int, unsigned int> cell_batch_range)
const;
2017 std::pair<unsigned int, unsigned int>
2019 const std::pair<unsigned int, unsigned int> face_batch_range)
const;
2039 std::pair<unsigned int, unsigned int>
2073 template <
typename StreamType>
2101 const internal::MatrixFreeFunctions::
2102 MappingInfo<dim, Number, VectorizedArrayType> &
2137 const unsigned int quad_index = 0,
2138 const unsigned int fe_base_element = 0,
2139 const unsigned int hp_active_fe_index = 0,
2140 const unsigned int hp_active_quad_index = 0)
const;
2146 VectorizedArrayType::size()> &
2206 template <
typename number2,
int q_dim>
2212 const std::vector<IndexSet> &locally_owned_set,
2222 template <
typename number2>
2226 const std::vector<IndexSet> &locally_owned_set,
2248 std::vector<ObserverPointer<const AffineConstraints<Number>>>
2255 std::vector<internal::MatrixFreeFunctions::DoFInfo>
dof_info;
2337 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2344 mutable std::list<std::pair<bool, AlignedVector<Number>>>
2367template <
int dim,
typename Number,
typename VectorizedArrayType>
2368template <
typename T>
2373 vec.
resize(this->n_cell_batches() + this->n_ghost_cell_batches());
2378template <
int dim,
typename Number,
typename VectorizedArrayType>
2379template <
typename T>
2384 vec.
resize(this->n_inner_face_batches() + this->n_boundary_face_batches() +
2385 this->n_ghost_inner_face_batches());
2390template <
int dim,
typename Number,
typename VectorizedArrayType>
2391template <
typename VectorType>
2395 const unsigned int comp)
const
2398 "This function is not supported for block vectors.");
2400 Assert(task_info.n_procs == 1,
2401 ExcMessage(
"This function can only be used in serial."));
2404 vec.reinit(dof_info[comp].vector_partitioner->size());
2409template <
int dim,
typename Number,
typename VectorizedArrayType>
2410template <
typename Number2>
2414 const unsigned int comp)
const
2417 vec.
reinit(dof_info[comp].vector_partitioner, task_info.communicator_sm);
2422template <
int dim,
typename Number,
typename VectorizedArrayType>
2423inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2425 const unsigned int comp)
const
2428 return dof_info[comp].vector_partitioner;
2433template <
int dim,
typename Number,
typename VectorizedArrayType>
2434inline const std::vector<unsigned int> &
2436 const unsigned int comp)
const
2439 return dof_info[comp].constrained_dofs;
2444template <
int dim,
typename Number,
typename VectorizedArrayType>
2449 return dof_handlers.size();
2454template <
int dim,
typename Number,
typename VectorizedArrayType>
2457 const unsigned int dof_no)
const
2461 return dof_handlers[dof_no]->get_fe().n_base_elements();
2466template <
int dim,
typename Number,
typename VectorizedArrayType>
2475template <
int dim,
typename Number,
typename VectorizedArrayType>
2479 return task_info.n_active_cells;
2484template <
int dim,
typename Number,
typename VectorizedArrayType>
2488 return *(task_info.cell_partition_data.end() - 2);
2493template <
int dim,
typename Number,
typename VectorizedArrayType>
2497 return *(task_info.cell_partition_data.end() - 1) -
2498 *(task_info.cell_partition_data.end() - 2);
2503template <
int dim,
typename Number,
typename VectorizedArrayType>
2507 if (task_info.face_partition_data.empty())
2509 return task_info.face_partition_data.back();
2514template <
int dim,
typename Number,
typename VectorizedArrayType>
2518 if (task_info.face_partition_data.empty())
2520 return task_info.boundary_partition_data.back() -
2521 task_info.face_partition_data.back();
2526template <
int dim,
typename Number,
typename VectorizedArrayType>
2530 if (task_info.face_partition_data.empty())
2532 return face_info.faces.size() - task_info.boundary_partition_data.back();
2537template <
int dim,
typename Number,
typename VectorizedArrayType>
2540 const unsigned int face_batch_index)
const
2542 Assert(face_batch_index >= task_info.boundary_partition_data[0] &&
2543 face_batch_index < task_info.boundary_partition_data.back(),
2545 task_info.boundary_partition_data[0],
2546 task_info.boundary_partition_data.back()));
2552template <
int dim,
typename Number,
typename VectorizedArrayType>
2555 const unsigned int cell_batch_index,
2556 const unsigned int face_number)
const
2560 Assert(face_info.cell_and_face_boundary_id.size(0) >= n_cell_batches(),
2564 for (
unsigned int v = 0;
2565 v < n_active_entries_per_cell_batch(cell_batch_index);
2568 face_info.cell_and_face_boundary_id(cell_batch_index, face_number, v);
2574template <
int dim,
typename Number,
typename VectorizedArrayType>
2575inline const internal::MatrixFreeFunctions::
2576 MappingInfo<dim, Number, VectorizedArrayType> &
2579 return mapping_info;
2584template <
int dim,
typename Number,
typename VectorizedArrayType>
2587 const unsigned int dof_index)
const
2590 return dof_info[dof_index];
2595template <
int dim,
typename Number,
typename VectorizedArrayType>
2599 return constraint_pool_row_index.size() - 1;
2604template <
int dim,
typename Number,
typename VectorizedArrayType>
2605inline const Number *
2607 const unsigned int row)
const
2610 return constraint_pool_data.empty() ?
2612 constraint_pool_data.data() + constraint_pool_row_index[row];
2617template <
int dim,
typename Number,
typename VectorizedArrayType>
2618inline const Number *
2620 const unsigned int row)
const
2623 return constraint_pool_data.empty() ?
2625 constraint_pool_data.data() + constraint_pool_row_index[row + 1];
2630template <
int dim,
typename Number,
typename VectorizedArrayType>
2631inline std::pair<unsigned int, unsigned int>
2633 const std::pair<unsigned int, unsigned int> &range,
2634 const unsigned int degree,
2635 const unsigned int dof_handler_component)
const
2637 if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2640 dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2642 dof_info[dof_handler_component].fe_index_conversion[0].
size(), 1);
2643 if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2646 return {range.second, range.second};
2650 dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2651 if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2652 return {range.second, range.second};
2654 return create_cell_subrange_hp_by_index(range,
2656 dof_handler_component);
2661template <
int dim,
typename Number,
typename VectorizedArrayType>
2664 const unsigned int cell_batch_index)
const
2667 return VectorizedArrayType::size() > 1 &&
2668 cell_level_index[(cell_batch_index + 1) * VectorizedArrayType::size() -
2669 1] == cell_level_index[(cell_batch_index + 1) *
2670 VectorizedArrayType::size() -
2676template <
int dim,
typename Number,
typename VectorizedArrayType>
2680 return shape_info.size(2);
2684template <
int dim,
typename Number,
typename VectorizedArrayType>
2687 const std::pair<unsigned int, unsigned int> range)
const
2689 const auto &fe_indices =
2690 dof_info[first_hp_dof_handler_index].cell_active_fe_index;
2692 if (fe_indices.empty() ==
true ||
2693 dof_handlers[first_hp_dof_handler_index]->get_fe_collection().size() == 1)
2696 const auto index = fe_indices[range.first];
2698 for (
unsigned int i = range.first; i < range.second; ++i)
2706template <
int dim,
typename Number,
typename VectorizedArrayType>
2709 const std::pair<unsigned int, unsigned int> range,
2710 const bool is_interior_face)
const
2712 const auto &fe_indices =
2713 dof_info[first_hp_dof_handler_index].cell_active_fe_index;
2715 if (fe_indices.empty() ==
true)
2718 if (is_interior_face)
2720 const unsigned int index =
2721 fe_indices[face_info.faces[range.first].cells_interior[0] /
2722 VectorizedArrayType::size()];
2724 for (
unsigned int i = range.first; i < range.second; ++i)
2726 fe_indices[face_info.faces[i].cells_interior[0] /
2727 VectorizedArrayType::size()]);
2733 const unsigned int index =
2734 fe_indices[face_info.faces[range.first].cells_exterior[0] /
2735 VectorizedArrayType::size()];
2737 for (
unsigned int i = range.first; i < range.second; ++i)
2739 fe_indices[face_info.faces[i].cells_exterior[0] /
2740 VectorizedArrayType::size()]);
2748template <
int dim,
typename Number,
typename VectorizedArrayType>
2751 const unsigned int cell_batch_index)
const
2755 const std::vector<unsigned char> &n_lanes_filled =
2756 dof_info[0].n_vectorization_lanes_filled
2760 return n_lanes_filled[cell_batch_index];
2765template <
int dim,
typename Number,
typename VectorizedArrayType>
2768 const unsigned int face_batch_index)
const
2772 const std::vector<unsigned char> &n_lanes_filled =
2773 dof_info[0].n_vectorization_lanes_filled
2776 return n_lanes_filled[face_batch_index];
2781template <
int dim,
typename Number,
typename VectorizedArrayType>
2784 const unsigned int dof_handler_index,
2785 const unsigned int active_fe_index)
const
2787 return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2792template <
int dim,
typename Number,
typename VectorizedArrayType>
2795 const unsigned int quad_index,
2796 const unsigned int active_fe_index)
const
2799 return mapping_info.cell_data[quad_index]
2800 .descriptor[active_fe_index]
2806template <
int dim,
typename Number,
typename VectorizedArrayType>
2809 const unsigned int dof_handler_index,
2810 const unsigned int active_fe_index)
const
2812 return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2817template <
int dim,
typename Number,
typename VectorizedArrayType>
2820 const unsigned int quad_index,
2821 const unsigned int active_fe_index)
const
2824 return mapping_info.face_data[quad_index]
2825 .descriptor[active_fe_index]
2831template <
int dim,
typename Number,
typename VectorizedArrayType>
2834 const unsigned int dof_handler_index)
const
2836 return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2841template <
int dim,
typename Number,
typename VectorizedArrayType>
2844 const unsigned int dof_handler_index)
const
2846 return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2851template <
int dim,
typename Number,
typename VectorizedArrayType>
2854 const unsigned int dof_handler_index,
2855 const unsigned int index_quad,
2856 const unsigned int index_fe,
2857 const unsigned int active_fe_index,
2858 const unsigned int active_quad_index)
const
2861 const unsigned int ind =
2862 dof_info[dof_handler_index].global_base_element_offset + index_fe;
2867 return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2872template <
int dim,
typename Number,
typename VectorizedArrayType>
2874 VectorizedArrayType::size()> &
2876 const unsigned int face_batch_index)
const
2879 return face_info.faces[face_batch_index];
2884template <
int dim,
typename Number,
typename VectorizedArrayType>
2889 return face_info.cell_and_face_to_plain_faces;
2894template <
int dim,
typename Number,
typename VectorizedArrayType>
2897 const unsigned int quad_index,
2898 const unsigned int active_fe_index)
const
2901 return mapping_info.cell_data[quad_index]
2902 .descriptor[active_fe_index]
2908template <
int dim,
typename Number,
typename VectorizedArrayType>
2911 const unsigned int quad_index,
2912 const unsigned int active_fe_index)
const
2915 return mapping_info.face_data[quad_index]
2916 .descriptor[active_fe_index]
2922template <
int dim,
typename Number,
typename VectorizedArrayType>
2925 const std::pair<unsigned int, unsigned int> range)
const
2927 auto result = get_cell_category(range.first);
2929 for (
unsigned int i = range.first; i < range.second; ++i)
2930 result =
std::max(result, get_cell_category(i));
2937template <
int dim,
typename Number,
typename VectorizedArrayType>
2938inline std::pair<unsigned int, unsigned int>
2940 const std::pair<unsigned int, unsigned int> range)
const
2942 auto result = get_face_category(range.first);
2944 for (
unsigned int i = range.first; i < range.second; ++i)
2946 result.first =
std::max(result.first, get_face_category(i).
first);
2947 result.second =
std::max(result.second, get_face_category(i).
second);
2955template <
int dim,
typename Number,
typename VectorizedArrayType>
2958 const unsigned int cell_batch_index)
const
2963 dof_info[first_hp_dof_handler_index].cell_active_fe_index.size());
2964 if (dof_info[first_hp_dof_handler_index].cell_active_fe_index.empty())
2967 return dof_info[first_hp_dof_handler_index]
2968 .cell_active_fe_index[cell_batch_index];
2973template <
int dim,
typename Number,
typename VectorizedArrayType>
2974inline std::pair<unsigned int, unsigned int>
2976 const unsigned int face_batch_index)
const
2979 if (dof_info[first_hp_dof_handler_index].cell_active_fe_index.empty())
2980 return std::make_pair(0U, 0U);
2982 std::pair<unsigned int, unsigned int> result = std::make_pair(0U, 0U);
2983 for (
unsigned int v = 0;
2984 v < VectorizedArrayType::size() &&
2985 face_info.faces[face_batch_index].cells_interior[v] !=
2990 dof_info[first_hp_dof_handler_index].cell_active_fe_index
2991 [face_info.faces[face_batch_index].cells_interior[v] /
2992 VectorizedArrayType::size()]);
2993 if (face_info.faces[face_batch_index].cells_exterior[0] !=
2995 for (
unsigned int v = 0;
2996 v < VectorizedArrayType::size() &&
2997 face_info.faces[face_batch_index].cells_exterior[v] !=
3002 dof_info[first_hp_dof_handler_index].cell_active_fe_index
3003 [face_info.faces[face_batch_index].cells_exterior[v] /
3004 VectorizedArrayType::size()]);
3012template <
int dim,
typename Number,
typename VectorizedArrayType>
3016 return indices_are_initialized;
3021template <
int dim,
typename Number,
typename VectorizedArrayType>
3025 return mapping_is_initialized;
3029template <
int dim,
typename Number,
typename VectorizedArrayType>
3038template <
int dim,
typename Number,
typename VectorizedArrayType>
3043 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3044 list_type &
data = scratch_pad.get();
3045 for (
typename list_type::iterator it =
data.begin(); it !=
data.end(); ++it)
3046 if (it->first ==
false)
3052 return &
data.front().second;
3057template <
int dim,
typename Number,
typename VectorizedArrayType>
3063 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3064 list_type &
data = scratch_pad.get();
3065 for (
typename list_type::iterator it =
data.begin(); it !=
data.end(); ++it)
3066 if (&it->second == scratch)
3077template <
int dim,
typename Number,
typename VectorizedArrayType>
3083 scratch_pad_non_threadsafe.
begin();
3084 it != scratch_pad_non_threadsafe.end();
3086 if (it->first ==
false)
3091 scratch_pad_non_threadsafe.push_front(
3093 return &scratch_pad_non_threadsafe.front().second;
3098template <
int dim,
typename Number,
typename VectorizedArrayType>
3105 scratch_pad_non_threadsafe.
begin();
3106 it != scratch_pad_non_threadsafe.end();
3108 if (&it->second == scratch)
3123 namespace MatrixFreeImplementation
3125 template <
int dim,
int spacedim>
3126 inline std::vector<IndexSet>
3127 extract_locally_owned_index_sets(
3129 const unsigned int level)
3131 std::vector<IndexSet> locally_owned_set;
3132 locally_owned_set.reserve(dofh.size());
3133 for (
unsigned int j = 0; j < dofh.size(); ++j)
3135 locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
3137 locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(
level));
3138 return locally_owned_set;
3145template <
int dim,
typename Number,
typename VectorizedArrayType>
3146template <
typename QuadratureType,
typename number2,
typename MappingType>
3149 const MappingType &mapping,
3152 const QuadratureType &quad,
3156 std::vector<const DoFHandler<dim, dim> *> dof_handlers;
3157 std::vector<const AffineConstraints<number2> *> constraints;
3159 dof_handlers.push_back(&dof_handler);
3160 constraints.push_back(&constraints_in);
3162 std::vector<IndexSet> locally_owned_sets =
3163 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3164 dof_handlers, additional_data.
mg_level);
3166 std::vector<hp::QCollection<dim>> quad_hp;
3167 quad_hp.emplace_back(quad);
3179template <
int dim,
typename Number,
typename VectorizedArrayType>
3180template <
typename QuadratureType,
typename number2,
typename MappingType>
3183 const MappingType &mapping,
3186 const QuadratureType &quad,
3190 std::vector<IndexSet> locally_owned_set =
3191 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3192 dof_handler, additional_data.
mg_level);
3193 std::vector<hp::QCollection<dim>> quad_hp;
3194 quad_hp.emplace_back(quad);
3206template <
int dim,
typename Number,
typename VectorizedArrayType>
3207template <
typename QuadratureType,
typename number2,
typename MappingType>
3210 const MappingType &mapping,
3213 const std::vector<QuadratureType> &quad,
3217 std::vector<IndexSet> locally_owned_set =
3218 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3219 dof_handler, additional_data.
mg_level);
3220 std::vector<hp::QCollection<dim>> quad_hp;
3221 for (
unsigned int q = 0; q < quad.size(); ++q)
3222 quad_hp.emplace_back(quad[q]);
3249 template <
int dim,
typename Number,
typename VectorizedArrayType>
3250 struct VectorDataExchange
3257 static constexpr unsigned int channel_shift = 20;
3266 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3267 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3268 DataAccessOnFaces vector_face_access,
3269 const unsigned
int n_components)
3270 : matrix_free(matrix_free)
3271 , vector_face_access(
3272 matrix_free.get_task_info().face_partition_data.empty() ?
3273 ::
MatrixFree<dim, Number, VectorizedArrayType>::
3274 DataAccessOnFaces::unspecified :
3276 , ghosts_were_set(false)
3277# ifdef DEAL_II_WITH_MPI
3278 , tmp_data(n_components)
3279 , requests(n_components)
3283 if (this->vector_face_access !=
3285 DataAccessOnFaces::unspecified)
3286 for (unsigned
int c = 0; c < matrix_free.n_components(); ++c)
3288 matrix_free.get_dof_info(c).vector_exchanger_face_variants.size(),
3297 ~VectorDataExchange()
3299# ifdef DEAL_II_WITH_MPI
3300 for (
unsigned int i = 0; i < tmp_data.size(); ++i)
3301 if (tmp_data[i] !=
nullptr)
3302 matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3312 template <
typename VectorType>
3314 find_vector_in_mf(
const VectorType &vec,
3315 const bool check_global_compatibility =
true)
const
3318 for (
unsigned int c = 0; c < matrix_free.n_components(); ++c)
3319 if (vec.get_partitioner().get() ==
3320 matrix_free.get_dof_info(c).vector_partitioner.get())
3324 for (
unsigned int c = 0; c < matrix_free.n_components(); ++c)
3325 if (check_global_compatibility ?
3326 vec.get_partitioner()->is_globally_compatible(
3327 *matrix_free.get_dof_info(c).vector_partitioner) :
3328 vec.get_partitioner()->is_compatible(
3329 *matrix_free.get_dof_info(c).vector_partitioner))
3345 get_partitioner(
const unsigned int mf_component)
const
3348 .vector_exchanger_face_variants.size(),
3350 if (vector_face_access ==
3352 DataAccessOnFaces::
none)
3353 return *matrix_free.get_dof_info(mf_component)
3354 .vector_exchanger_face_variants[0];
3355 else if (vector_face_access ==
3357 DataAccessOnFaces::
values)
3358 return *matrix_free.get_dof_info(mf_component)
3359 .vector_exchanger_face_variants[1];
3360 else if (vector_face_access ==
3363 return *matrix_free.get_dof_info(mf_component)
3364 .vector_exchanger_face_variants[2];
3365 else if (vector_face_access ==
3367 DataAccessOnFaces::values_all_faces)
3368 return *matrix_free.get_dof_info(mf_component)
3369 .vector_exchanger_face_variants[3];
3370 else if (vector_face_access ==
3372 DataAccessOnFaces::gradients_all_faces)
3373 return *matrix_free.get_dof_info(mf_component)
3374 .vector_exchanger_face_variants[4];
3376 return *matrix_free.get_dof_info(mf_component).vector_exchanger.get();
3385 std::enable_if_t<is_not_parallel_vector<VectorType>,
VectorType>
3388 update_ghost_values_start(
const unsigned int ,
3389 const VectorType & )
3398 std::enable_if_t<!has_update_ghost_values_start<VectorType> &&
3399 !is_not_parallel_vector<VectorType>,
3402 update_ghost_values_start(
const unsigned int component_in_block_vector,
3403 const VectorType &vec)
3405 (void)component_in_block_vector;
3406 const bool ghosts_set = vec.has_ghost_elements();
3408 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3409 ghosts_set ==
false,
3414 ghosts_were_set =
true;
3418 vec.update_ghost_values();
3429 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3430 !has_exchange_on_subset<VectorType>,
3433 update_ghost_values_start(
const unsigned int component_in_block_vector,
3434 const VectorType &vec)
3436 (void)component_in_block_vector;
3437 const bool ghosts_set = vec.has_ghost_elements();
3439 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3440 ghosts_set ==
false,
3445 ghosts_were_set =
true;
3449 vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3461 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3462 has_exchange_on_subset<VectorType>,
3465 update_ghost_values_start(
const unsigned int component_in_block_vector,
3466 const VectorType &vec)
3468 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3469 "Type mismatch between VectorType and VectorDataExchange");
3470 (void)component_in_block_vector;
3471 const bool ghosts_set = vec.has_ghost_elements();
3473 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3474 ghosts_set ==
false,
3479 ghosts_were_set =
true;
3483 if (vec.size() != 0)
3485# ifdef DEAL_II_WITH_MPI
3486 const unsigned int mf_component = find_vector_in_mf(vec);
3488 const auto &part = get_partitioner(mf_component);
3490 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3491 part.n_import_sm_procs() == 0)
3494 tmp_data[component_in_block_vector] =
3495 matrix_free.acquire_scratch_data_non_threadsafe();
3496 tmp_data[component_in_block_vector]->resize_fast(
3497 part.n_import_indices());
3500 part.export_to_ghosted_array_start(
3501 component_in_block_vector * 2 + channel_shift,
3503 vec.shared_vector_data(),
3505 part.locally_owned_size(),
3506 matrix_free.get_dof_info(mf_component)
3507 .vector_partitioner->n_ghost_indices()),
3509 part.n_import_indices()),
3510 this->requests[component_in_block_vector]);
3522 std::enable_if_t<!has_update_ghost_values_start<VectorType>,
3525 update_ghost_values_finish(
const unsigned int ,
3526 const VectorType & )
3537 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3538 !has_exchange_on_subset<VectorType>,
3541 update_ghost_values_finish(
const unsigned int component_in_block_vector,
3542 const VectorType &vec)
3544 (void)component_in_block_vector;
3546 if (ghosts_were_set)
3549 vec.update_ghost_values_finish();
3561 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3562 has_exchange_on_subset<VectorType>,
3565 update_ghost_values_finish(
const unsigned int component_in_block_vector,
3566 const VectorType &vec)
3568 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3569 "Type mismatch between VectorType and VectorDataExchange");
3570 (void)component_in_block_vector;
3572 if (ghosts_were_set)
3575 if (vec.size() != 0)
3577# ifdef DEAL_II_WITH_MPI
3581 const unsigned int mf_component = find_vector_in_mf(vec);
3583 const auto &part = get_partitioner(mf_component);
3585 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3586 part.n_import_sm_procs() != 0)
3588 part.export_to_ghosted_array_finish(
3590 vec.shared_vector_data(),
3592 part.locally_owned_size(),
3593 matrix_free.get_dof_info(mf_component)
3594 .vector_partitioner->n_ghost_indices()),
3595 this->requests[component_in_block_vector]);
3597 matrix_free.release_scratch_data_non_threadsafe(
3598 tmp_data[component_in_block_vector]);
3599 tmp_data[component_in_block_vector] =
nullptr;
3605 vec.set_ghost_state(
true);
3614 std::enable_if_t<is_not_parallel_vector<VectorType>,
VectorType>
3617 compress_start(
const unsigned int ,
3628 std::enable_if_t<!has_compress_start<VectorType> &&
3629 !is_not_parallel_vector<VectorType>,
3632 compress_start(
const unsigned int component_in_block_vector,
3635 (void)component_in_block_vector;
3648 std::enable_if_t<has_compress_start<VectorType> &&
3649 !has_exchange_on_subset<VectorType>,
3652 compress_start(
const unsigned int component_in_block_vector,
3655 (void)component_in_block_vector;
3657 vec.compress_start(component_in_block_vector + channel_shift);
3669 std::enable_if_t<has_compress_start<VectorType> &&
3670 has_exchange_on_subset<VectorType>,
3673 compress_start(
const unsigned int component_in_block_vector,
3676 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3677 "Type mismatch between VectorType and VectorDataExchange");
3678 (void)component_in_block_vector;
3681 if (vec.size() != 0)
3683# ifdef DEAL_II_WITH_MPI
3684 const unsigned int mf_component = find_vector_in_mf(vec);
3686 const auto &part = get_partitioner(mf_component);
3688 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3689 part.n_import_sm_procs() == 0)
3692 tmp_data[component_in_block_vector] =
3693 matrix_free.acquire_scratch_data_non_threadsafe();
3694 tmp_data[component_in_block_vector]->resize_fast(
3695 part.n_import_indices());
3698 part.import_from_ghosted_array_start(
3700 component_in_block_vector * 2 + channel_shift,
3702 vec.shared_vector_data(),
3704 matrix_free.get_dof_info(mf_component)
3705 .vector_partitioner->n_ghost_indices()),
3707 part.n_import_indices()),
3708 this->requests[component_in_block_vector]);
3721 std::enable_if_t<!has_compress_start<VectorType>,
VectorType> * =
nullptr>
3723 compress_finish(
const unsigned int ,
3735 std::enable_if_t<has_compress_start<VectorType> &&
3736 !has_exchange_on_subset<VectorType>,
3739 compress_finish(
const unsigned int component_in_block_vector,
3742 (void)component_in_block_vector;
3755 std::enable_if_t<has_compress_start<VectorType> &&
3756 has_exchange_on_subset<VectorType>,
3759 compress_finish(
const unsigned int component_in_block_vector,
3762 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3763 "Type mismatch between VectorType and VectorDataExchange");
3764 (void)component_in_block_vector;
3765 if (vec.size() != 0)
3767# ifdef DEAL_II_WITH_MPI
3771 const unsigned int mf_component = find_vector_in_mf(vec);
3773 const auto &part = get_partitioner(mf_component);
3775 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3776 part.n_import_sm_procs() != 0)
3778 part.import_from_ghosted_array_finish(
3781 vec.shared_vector_data(),
3783 matrix_free.get_dof_info(mf_component)
3784 .vector_partitioner->n_ghost_indices()),
3786 tmp_data[component_in_block_vector]->begin(),
3787 part.n_import_indices()),
3788 this->requests[component_in_block_vector]);
3790 matrix_free.release_scratch_data_non_threadsafe(
3791 tmp_data[component_in_block_vector]);
3792 tmp_data[component_in_block_vector] =
nullptr;
3798 MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3811 std::enable_if_t<is_not_parallel_vector<VectorType>,
VectorType>
3814 reset_ghost_values(
const VectorType & )
const
3824 std::enable_if_t<!has_exchange_on_subset<VectorType> &&
3825 !is_not_parallel_vector<VectorType>,
3828 reset_ghost_values(
const VectorType &vec)
const
3830 if (ghosts_were_set ==
true)
3833 vec.zero_out_ghost_values();
3844 std::enable_if_t<has_exchange_on_subset<VectorType>,
VectorType>
3847 reset_ghost_values(
const VectorType &vec)
const
3849 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3850 "Type mismatch between VectorType and VectorDataExchange");
3851 if (ghosts_were_set ==
true)
3854 if (vec.size() != 0)
3856# ifdef DEAL_II_WITH_MPI
3859 const unsigned int mf_component = find_vector_in_mf(vec);
3861 const auto &part = get_partitioner(mf_component);
3863 if (part.n_ghost_indices() > 0)
3865 part.reset_ghost_values(
3867 part.locally_owned_size(),
3868 matrix_free.get_dof_info(mf_component)
3869 .vector_partitioner->n_ghost_indices()));
3875 vec.set_ghost_state(
false);
3886 std::enable_if_t<has_exchange_on_subset<VectorType>,
VectorType>
3889 zero_vector_region(
const unsigned int range_index, VectorType &vec)
const
3891 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3892 "Type mismatch between VectorType and VectorDataExchange");
3897 const unsigned int mf_component = find_vector_in_mf(vec,
false);
3899 matrix_free.get_dof_info(mf_component);
3907 for (
unsigned int id =
3927 std::enable_if_t<!has_exchange_on_subset<VectorType>,
VectorType>
3929 std::enable_if_t<has_assignment_operator<VectorType>,
VectorType>
3932 zero_vector_region(
const unsigned int range_index, VectorType &vec)
const
3936 if constexpr (std::is_same_v<
3940 for (
unsigned int i = 0; i < vec.size(); ++i)
3941 vec[i] =
typename VectorType::value_type();
3944 vec =
typename VectorType::value_type();
3955 zero_vector_region(
const unsigned int, ...)
const
3959 "Zeroing is only implemented for vector types "
3960 "which provide operator=(const VectorType::value_type)"));
3965 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
3966 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3967 DataAccessOnFaces vector_face_access;
3968 bool ghosts_were_set;
3969# ifdef DEAL_II_WITH_MPI
3970 std::vector<AlignedVector<Number> *> tmp_data;
3971 std::vector<std::vector<MPI_Request>> requests;
3975 template <
typename VectorStruct>
3977 n_components(
const VectorStruct &vec);
3979 template <
typename VectorStruct>
3981 n_components_block(
const VectorStruct &vec,
const std::bool_constant<true>)
3983 unsigned int components = 0;
3984 for (
unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
3985 components += n_components(vec.block(bl));
3989 template <
typename VectorStruct>
3991 n_components_block(
const VectorStruct &,
const std::bool_constant<false>)
3996 template <
typename VectorStruct>
3998 n_components(
const VectorStruct &vec)
4000 return n_components_block(
4004 template <
typename VectorStruct>
4006 n_components(
const std::vector<VectorStruct> &vec)
4008 unsigned int components = 0;
4009 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4010 components += n_components_block(
4015 template <
typename VectorStruct>
4017 n_components(
const std::vector<VectorStruct *> &vec)
4019 unsigned int components = 0;
4020 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4021 components += n_components_block(
4033 template <
typename VectorStruct,
4034 std::enable_if_t<!has_communication_block_size<VectorStruct>,
4035 VectorStruct> * =
nullptr>
4036 constexpr unsigned int
4037 get_communication_block_size(
const VectorStruct &)
4044 template <
typename VectorStruct,
4045 std::enable_if_t<has_communication_block_size<VectorStruct>,
4046 VectorStruct> * =
nullptr>
4047 constexpr unsigned int
4048 get_communication_block_size(
const VectorStruct &)
4050 return VectorStruct::communication_block_size;
4056 std::enable_if_t<is_not_parallel_vector<VectorType>,
VectorType> * =
4059 has_ghost_elements(
const VectorType &vec)
4068 std::enable_if_t<!is_not_parallel_vector<VectorType>,
VectorType>
4071 has_ghost_elements(
const VectorType &vec)
4073 return vec.has_ghost_elements();
4088 typename VectorStruct,
4090 typename VectorizedArrayType,
4091 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4094 update_ghost_values_start(
4095 const VectorStruct &vec,
4096 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4097 const unsigned int channel = 0)
4099 if (get_communication_block_size(vec) < vec.n_blocks())
4101 const bool ghosts_set = vec.has_ghost_elements();
4103 Assert(exchanger.matrix_free.get_task_info()
4104 .allow_ghosted_vectors_in_loops ||
4105 ghosts_set ==
false,
4110 exchanger.ghosts_were_set =
true;
4114 vec.update_ghost_values();
4118 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4119 update_ghost_values_start(vec.block(i), exchanger, channel + i);
4127 typename VectorStruct,
4129 typename VectorizedArrayType,
4130 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4133 update_ghost_values_start(
4134 const VectorStruct &vec,
4135 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4136 const unsigned int channel = 0)
4138 exchanger.update_ghost_values_start(channel, vec);
4145 typename VectorStruct,
4147 typename VectorizedArrayType>
4149 update_ghost_values_start(
4150 const std::vector<VectorStruct> &vec,
4151 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4153 unsigned int component_index = 0;
4154 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4156 update_ghost_values_start(vec[comp], exchanger, component_index);
4157 component_index += n_components(vec[comp]);
4165 typename VectorStruct,
4167 typename VectorizedArrayType>
4169 update_ghost_values_start(
4170 const std::vector<VectorStruct *> &vec,
4171 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4173 unsigned int component_index = 0;
4174 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4176 update_ghost_values_start(*vec[comp], exchanger, component_index);
4177 component_index += n_components(*vec[comp]);
4189 typename VectorStruct,
4191 typename VectorizedArrayType,
4192 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4195 update_ghost_values_finish(
4196 const VectorStruct &vec,
4197 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4198 const unsigned int channel = 0)
4200 if (get_communication_block_size(vec) < vec.n_blocks())
4206 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4207 update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4214 typename VectorStruct,
4216 typename VectorizedArrayType,
4217 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4220 update_ghost_values_finish(
4221 const VectorStruct &vec,
4222 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4223 const unsigned int channel = 0)
4225 exchanger.update_ghost_values_finish(channel, vec);
4232 typename VectorStruct,
4234 typename VectorizedArrayType>
4236 update_ghost_values_finish(
4237 const std::vector<VectorStruct> &vec,
4238 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4240 unsigned int component_index = 0;
4241 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4243 update_ghost_values_finish(vec[comp], exchanger, component_index);
4244 component_index += n_components(vec[comp]);
4252 typename VectorStruct,
4254 typename VectorizedArrayType>
4256 update_ghost_values_finish(
4257 const std::vector<VectorStruct *> &vec,
4258 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4260 unsigned int component_index = 0;
4261 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4263 update_ghost_values_finish(*vec[comp], exchanger, component_index);
4264 component_index += n_components(*vec[comp]);
4276 typename VectorStruct,
4278 typename VectorizedArrayType,
4279 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4284 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4285 const unsigned int channel = 0)
4287 if (get_communication_block_size(vec) < vec.n_blocks())
4290 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4291 compress_start(vec.block(i), exchanger, channel + i);
4298 typename VectorStruct,
4300 typename VectorizedArrayType,
4301 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4306 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4307 const unsigned int channel = 0)
4309 exchanger.compress_start(channel, vec);
4316 typename VectorStruct,
4318 typename VectorizedArrayType>
4321 std::vector<VectorStruct> &vec,
4322 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4324 unsigned int component_index = 0;
4325 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4327 compress_start(vec[comp], exchanger, component_index);
4328 component_index += n_components(vec[comp]);
4336 typename VectorStruct,
4338 typename VectorizedArrayType>
4341 std::vector<VectorStruct *> &vec,
4342 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4344 unsigned int component_index = 0;
4345 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4347 compress_start(*vec[comp], exchanger, component_index);
4348 component_index += n_components(*vec[comp]);
4360 typename VectorStruct,
4362 typename VectorizedArrayType,
4363 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4368 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4369 const unsigned int channel = 0)
4371 if (get_communication_block_size(vec) < vec.n_blocks())
4377 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4378 compress_finish(vec.block(i), exchanger, channel + i);
4385 typename VectorStruct,
4387 typename VectorizedArrayType,
4388 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4393 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4394 const unsigned int channel = 0)
4396 exchanger.compress_finish(channel, vec);
4403 typename VectorStruct,
4405 typename VectorizedArrayType>
4408 std::vector<VectorStruct> &vec,
4409 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4411 unsigned int component_index = 0;
4412 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4414 compress_finish(vec[comp], exchanger, component_index);
4415 component_index += n_components(vec[comp]);
4423 typename VectorStruct,
4425 typename VectorizedArrayType>
4428 std::vector<VectorStruct *> &vec,
4429 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4431 unsigned int component_index = 0;
4432 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4434 compress_finish(*vec[comp], exchanger, component_index);
4435 component_index += n_components(*vec[comp]);
4451 typename VectorStruct,
4453 typename VectorizedArrayType,
4454 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4458 const VectorStruct &vec,
4459 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4462 if (exchanger.ghosts_were_set ==
true)
4465 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4466 reset_ghost_values(vec.block(i), exchanger);
4473 typename VectorStruct,
4475 typename VectorizedArrayType,
4476 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4480 const VectorStruct &vec,
4481 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4484 if (exchanger.ghosts_were_set ==
true)
4487 exchanger.reset_ghost_values(vec);
4494 typename VectorStruct,
4496 typename VectorizedArrayType>
4499 const std::vector<VectorStruct> &vec,
4500 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4503 if (exchanger.ghosts_were_set ==
true)
4506 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4507 reset_ghost_values(vec[comp], exchanger);
4514 typename VectorStruct,
4516 typename VectorizedArrayType>
4519 const std::vector<VectorStruct *> &vec,
4520 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4523 if (exchanger.ghosts_were_set ==
true)
4526 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4527 reset_ghost_values(*vec[comp], exchanger);
4538 typename VectorStruct,
4540 typename VectorizedArrayType,
4541 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4545 const unsigned int range_index,
4547 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4549 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4550 exchanger.zero_vector_region(range_index, vec.block(i));
4557 typename VectorStruct,
4559 typename VectorizedArrayType,
4560 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4564 const unsigned int range_index,
4566 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4568 exchanger.zero_vector_region(range_index, vec);
4575 typename VectorStruct,
4577 typename VectorizedArrayType>
4580 const unsigned int range_index,
4581 std::vector<VectorStruct> &vec,
4582 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4584 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4585 zero_vector_region(range_index, vec[comp], exchanger);
4592 typename VectorStruct,
4594 typename VectorizedArrayType>
4597 const unsigned int range_index,
4598 std::vector<VectorStruct *> &vec,
4599 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4601 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4602 zero_vector_region(range_index, *vec[comp], exchanger);
4610 template <
typename VectorStruct1,
typename VectorStruct2>
4612 apply_operation_to_constrained_dofs(
const std::vector<unsigned int> &,
4613 const VectorStruct1 &,
4617 template <
typename Number>
4619 apply_operation_to_constrained_dofs(
4620 const std::vector<unsigned int> &constrained_dofs,
4624 for (
const unsigned int i : constrained_dofs)
4625 dst.local_element(i) = src.local_element(i);
4629 namespace MatrixFreeFunctions
4633 template <
typename,
typename,
typename,
typename,
bool>
4634 struct InterfaceSelector
4638 template <
typename MF,
4642 struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4644 using function_type = void (Container::*)(
4648 const
std::pair<unsigned
int, unsigned
int> &) const;
4652 template <
typename MF,
4656 struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4658 using function_type =
4659 void (Container::*)(const MF &,
4662 const
std::pair<unsigned
int, unsigned
int> &);
4670 template <
typename MF,
4675 class MFWorker :
public MFWorkerInterface
4679 using function_type =
typename MatrixFreeFunctions::
4680 InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4684 MFWorker(
const MF &matrix_free,
4685 const InVector &src,
4687 const bool zero_dst_vector_setting,
4688 const Container &container,
4689 function_type cell_function,
4690 function_type face_function,
4691 function_type boundary_function,
4692 const typename MF::DataAccessOnFaces src_vector_face_access =
4693 MF::DataAccessOnFaces::none,
4694 const typename MF::DataAccessOnFaces dst_vector_face_access =
4695 MF::DataAccessOnFaces::none,
4696 const std::function<
void(
const unsigned int,
const unsigned int)>
4697 &operation_before_loop = {},
4698 const std::function<void(
const unsigned int,
const unsigned int)>
4699 &operation_after_loop = {},
4700 const unsigned int dof_handler_index_pre_post = 0)
4701 : matrix_free(matrix_free)
4702 , container(const_cast<Container &>(container))
4703 , cell_function(cell_function)
4704 , face_function(face_function)
4705 , boundary_function(boundary_function)
4708 , src_data_exchanger(matrix_free,
4709 src_vector_face_access,
4711 , dst_data_exchanger(matrix_free,
4712 dst_vector_face_access,
4715 , zero_dst_vector_setting(zero_dst_vector_setting &&
4716 !src_and_dst_are_same)
4717 , operation_before_loop(operation_before_loop)
4718 , operation_after_loop(operation_after_loop)
4719 , dof_handler_index_pre_post(dof_handler_index_pre_post)
4721 Assert(!has_ghost_elements(dst),
4722 ExcMessage(
"The destination vector passed to the matrix-free "
4723 "loop is ghosted. This is not allowed."));
4728 cell(
const std::pair<unsigned int, unsigned int> &cell_range)
override
4730 if (cell_function !=
nullptr && cell_range.second > cell_range.first)
4731 for (
unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4733 const auto cell_subrange =
4734 matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4736 if (cell_subrange.second <= cell_subrange.first)
4740 cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4745 cell(
const unsigned int range_index)
override
4747 process_range(cell_function,
4748 matrix_free.get_task_info().cell_partition_data_hp_ptr,
4749 matrix_free.get_task_info().cell_partition_data_hp,
4754 face(
const unsigned int range_index)
override
4756 process_range(face_function,
4757 matrix_free.get_task_info().face_partition_data_hp_ptr,
4758 matrix_free.get_task_info().face_partition_data_hp,
4763 boundary(
const unsigned int range_index)
override
4765 process_range(boundary_function,
4766 matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4767 matrix_free.get_task_info().boundary_partition_data_hp,
4773 process_range(
const function_type &fu,
4774 const std::vector<unsigned int> &ptr,
4775 const std::vector<unsigned int> &
data,
4776 const unsigned int range_index)
4782 for (
unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4785 (container.*fu)(matrix_free,
4788 std::make_pair(
data[2 * i],
data[2 * i + 1]));
4800 vector_update_ghosts_start()
override
4802 if (!src_and_dst_are_same)
4803 internal::update_ghost_values_start(src, src_data_exchanger);
4808 vector_update_ghosts_finish()
override
4810 if (!src_and_dst_are_same)
4811 internal::update_ghost_values_finish(src, src_data_exchanger);
4816 vector_compress_start()
override
4818 internal::compress_start(dst, dst_data_exchanger);
4823 vector_compress_finish()
override
4825 internal::compress_finish(dst, dst_data_exchanger);
4826 if (!src_and_dst_are_same)
4827 internal::reset_ghost_values(src, src_data_exchanger);
4832 zero_dst_vector_range(
const unsigned int range_index)
override
4834 if (zero_dst_vector_setting)
4835 internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4839 cell_loop_pre_range(
const unsigned int range_index)
override
4841 if (operation_before_loop)
4844 matrix_free.get_dof_info(dof_handler_index_pre_post);
4851 operation_before_loop,
4858 for (
unsigned int id =
4869 cell_loop_post_range(
const unsigned int range_index)
override
4871 if (operation_after_loop)
4875 const std::vector<unsigned int> &partition_row_index =
4876 matrix_free.get_task_info().partition_row_index;
4878 partition_row_index[partition_row_index.size() - 2] - 1)
4879 apply_operation_to_constrained_dofs(
4880 matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
4885 matrix_free.get_dof_info(dof_handler_index_pre_post);
4892 operation_after_loop,
4899 for (
unsigned int id =
4910 const MF &matrix_free;
4911 Container &container;
4912 function_type cell_function;
4913 function_type face_function;
4914 function_type boundary_function;
4916 const InVector &src;
4918 VectorDataExchange<MF::dimension,
4919 typename MF::value_type,
4920 typename MF::vectorized_value_type>
4922 VectorDataExchange<MF::dimension,
4923 typename MF::value_type,
4924 typename MF::vectorized_value_type>
4926 const bool src_and_dst_are_same;
4927 const bool zero_dst_vector_setting;
4928 const std::function<void(
const unsigned int,
const unsigned int)>
4929 operation_before_loop;
4930 const std::function<void(
const unsigned int,
const unsigned int)>
4931 operation_after_loop;
4932 const unsigned int dof_handler_index_pre_post;
4941 template <
class MF,
typename InVector,
typename OutVector>
4942 struct MFClassWrapper
4944 using function_type =
4945 std::function<void(
const MF &,
4948 const std::pair<unsigned int, unsigned int> &)>;
4950 MFClassWrapper(
const function_type cell,
4951 const function_type face,
4952 const function_type boundary)
4955 , boundary(boundary)
4959 cell_integrator(
const MF &mf,
4961 const InVector &src,
4962 const std::pair<unsigned int, unsigned int> &range)
const
4965 cell(mf, dst, src, range);
4969 face_integrator(
const MF &mf,
4971 const InVector &src,
4972 const std::pair<unsigned int, unsigned int> &range)
const
4975 face(mf, dst, src, range);
4979 boundary_integrator(
4982 const InVector &src,
4983 const std::pair<unsigned int, unsigned int> &range)
const
4986 boundary(mf, dst, src, range);
4989 const function_type cell;
4990 const function_type face;
4991 const function_type boundary;
4998template <
int dim,
typename Number,
typename VectorizedArrayType>
4999template <
typename OutVector,
typename InVector>
5005 const std::pair<unsigned int, unsigned int> &)>
5008 const InVector &src,
5009 const bool zero_dst_vector)
const
5012 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5015 Wrapper wrap(cell_operation,
nullptr,
nullptr);
5016 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5026 &Wrapper::cell_integrator,
5027 &Wrapper::face_integrator,
5028 &Wrapper::boundary_integrator);
5030 task_info.loop(worker);
5035template <
int dim,
typename Number,
typename VectorizedArrayType>
5036template <
typename OutVector,
typename InVector>
5042 const std::pair<unsigned int, unsigned int> &)>
5045 const InVector &src,
5046 const std::function<
void(
const unsigned int,
const unsigned int)>
5047 &operation_before_loop,
5048 const std::function<
void(
const unsigned int,
const unsigned int)>
5049 &operation_after_loop,
5050 const unsigned int dof_handler_index_pre_post)
const
5053 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5056 Wrapper wrap(cell_operation,
nullptr,
nullptr);
5057 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5067 &Wrapper::cell_integrator,
5068 &Wrapper::face_integrator,
5069 &Wrapper::boundary_integrator,
5070 DataAccessOnFaces::none,
5071 DataAccessOnFaces::none,
5072 operation_before_loop,
5073 operation_after_loop,
5074 dof_handler_index_pre_post);
5076 task_info.loop(worker);
5081template <
int dim,
typename Number,
typename VectorizedArrayType>
5082template <
typename OutVector,
typename InVector>
5088 const std::pair<unsigned int, unsigned int> &)>
5093 const std::pair<unsigned int, unsigned int> &)>
5094 &inner_face_operation,
5098 const std::pair<unsigned int, unsigned int> &)>
5099 &boundary_face_operation,
5101 const InVector &src,
5102 const bool zero_dst_vector,
5103 const DataAccessOnFaces dst_vector_face_access,
5104 const DataAccessOnFaces src_vector_face_access)
const
5107 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5110 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5111 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5121 &Wrapper::cell_integrator,
5122 &Wrapper::face_integrator,
5123 &Wrapper::boundary_integrator,
5124 src_vector_face_access,
5125 dst_vector_face_access);
5127 task_info.loop(worker);
5132template <
int dim,
typename Number,
typename VectorizedArrayType>
5133template <
typename CLASS,
typename OutVector,
typename InVector>
5136 void (CLASS::*function_pointer)(
5137 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5140 const
std::pair<unsigned
int, unsigned
int> &) const,
5141 const CLASS *owning_class,
5143 const InVector &src,
5144 const
bool zero_dst_vector) const
5146 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5164template <
int dim,
typename Number,
typename VectorizedArrayType>
5165template <
typename CLASS,
typename OutVector,
typename InVector>
5168 void (CLASS::*function_pointer)(
5169 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5172 const
std::pair<unsigned
int, unsigned
int> &) const,
5173 const CLASS *owning_class,
5175 const InVector &src,
5176 const
std::function<void(const unsigned
int, const unsigned
int)>
5177 &operation_before_loop,
5178 const
std::function<void(const unsigned
int, const unsigned
int)>
5179 &operation_after_loop,
5180 const unsigned
int dof_handler_index_pre_post) const
5182 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5197 operation_before_loop,
5198 operation_after_loop,
5199 dof_handler_index_pre_post);
5205template <
int dim,
typename Number,
typename VectorizedArrayType>
5206template <
typename CLASS,
typename OutVector,
typename InVector>
5209 void (CLASS::*cell_operation)(
5210 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5213 const
std::pair<unsigned
int, unsigned
int> &) const,
5214 void (CLASS::*inner_face_operation)(
5215 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5218 const
std::pair<unsigned
int, unsigned
int> &) const,
5219 void (CLASS::*boundary_face_operation)(
5220 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5223 const
std::pair<unsigned
int, unsigned
int> &) const,
5224 const CLASS *owning_class,
5226 const InVector &src,
5227 const
bool zero_dst_vector,
5228 const DataAccessOnFaces dst_vector_face_access,
5229 const DataAccessOnFaces src_vector_face_access) const
5231 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5242 inner_face_operation,
5243 boundary_face_operation,
5244 src_vector_face_access,
5245 dst_vector_face_access);
5251template <
int dim,
typename Number,
typename VectorizedArrayType>
5252template <
typename CLASS,
typename OutVector,
typename InVector>
5255 void (CLASS::*function_pointer)(
5256 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5259 const
std::pair<unsigned
int, unsigned
int> &),
5260 CLASS *owning_class,
5262 const InVector &src,
5263 const
bool zero_dst_vector) const
5265 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5283template <
int dim,
typename Number,
typename VectorizedArrayType>
5284template <
typename CLASS,
typename OutVector,
typename InVector>
5287 void (CLASS::*function_pointer)(
5288 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5291 const
std::pair<unsigned
int, unsigned
int> &),
5292 CLASS *owning_class,
5294 const InVector &src,
5295 const
std::function<void(const unsigned
int, const unsigned
int)>
5296 &operation_before_loop,
5297 const
std::function<void(const unsigned
int, const unsigned
int)>
5298 &operation_after_loop,
5299 const unsigned
int dof_handler_index_pre_post) const
5301 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5316 operation_before_loop,
5317 operation_after_loop,
5318 dof_handler_index_pre_post);
5324template <
int dim,
typename Number,
typename VectorizedArrayType>
5325template <
typename CLASS,
typename OutVector,
typename InVector>
5328 void (CLASS::*cell_operation)(
5329 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5332 const
std::pair<unsigned
int, unsigned
int> &),
5333 void (CLASS::*inner_face_operation)(
5334 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5337 const
std::pair<unsigned
int, unsigned
int> &),
5338 void (CLASS::*boundary_face_operation)(
5339 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5342 const
std::pair<unsigned
int, unsigned
int> &),
5343 CLASS *owning_class,
5345 const InVector &src,
5346 const
bool zero_dst_vector,
5347 const DataAccessOnFaces dst_vector_face_access,
5348 const DataAccessOnFaces src_vector_face_access) const
5350 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5361 inner_face_operation,
5362 boundary_face_operation,
5363 src_vector_face_access,
5364 dst_vector_face_access);
5370template <
int dim,
typename Number,
typename VectorizedArrayType>
5371template <
typename OutVector,
typename InVector>
5377 const std::pair<unsigned int, unsigned int> &)>
5382 const std::pair<unsigned int, unsigned int> &)>
5383 &inner_face_operation,
5387 const std::pair<unsigned int, unsigned int> &)>
5388 &boundary_face_operation,
5390 const InVector &src,
5391 const std::function<
void(
const unsigned int,
const unsigned int)>
5392 &operation_before_loop,
5393 const std::function<
void(
const unsigned int,
const unsigned int)>
5394 &operation_after_loop,
5395 const unsigned int dof_handler_index_pre_post,
5396 const DataAccessOnFaces dst_vector_face_access,
5397 const DataAccessOnFaces src_vector_face_access)
const
5400 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5403 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5404 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5414 &Wrapper::cell_integrator,
5415 &Wrapper::face_integrator,
5416 &Wrapper::boundary_integrator,
5417 src_vector_face_access,
5418 dst_vector_face_access,
5419 operation_before_loop,
5420 operation_after_loop,
5421 dof_handler_index_pre_post);
5423 task_info.loop(worker);
5428template <
int dim,
typename Number,
typename VectorizedArrayType>
5429template <
typename CLASS,
typename OutVector,
typename InVector>
5432 void (CLASS::*cell_operation)(const
MatrixFree &,
5435 const
std::pair<unsigned
int, unsigned
int> &)
5437 void (CLASS::*inner_face_operation)(
5441 const
std::pair<unsigned
int, unsigned
int> &) const,
5442 void (CLASS::*boundary_face_operation)(
5446 const
std::pair<unsigned
int, unsigned
int> &) const,
5447 const CLASS *owning_class,
5449 const InVector &src,
5450 const
std::function<void(const unsigned
int, const unsigned
int)>
5451 &operation_before_loop,
5452 const
std::function<void(const unsigned
int, const unsigned
int)>
5453 &operation_after_loop,
5454 const unsigned
int dof_handler_index_pre_post,
5455 const DataAccessOnFaces dst_vector_face_access,
5456 const DataAccessOnFaces src_vector_face_access) const
5458 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5469 inner_face_operation,
5470 boundary_face_operation,
5471 src_vector_face_access,
5472 dst_vector_face_access,
5473 operation_before_loop,
5474 operation_after_loop,
5475 dof_handler_index_pre_post);
5481template <
int dim,
typename Number,
typename VectorizedArrayType>
5482template <
typename CLASS,
typename OutVector,
typename InVector>
5485 void (CLASS::*cell_operation)(const
MatrixFree &,
5488 const
std::pair<unsigned
int, unsigned
int> &),
5489 void (CLASS::*inner_face_operation)(
5493 const
std::pair<unsigned
int, unsigned
int> &),
5494 void (CLASS::*boundary_face_operation)(
5498 const
std::pair<unsigned
int, unsigned
int> &),
5499 const CLASS *owning_class,
5501 const InVector &src,
5502 const
std::function<void(const unsigned
int, const unsigned
int)>
5503 &operation_before_loop,
5504 const
std::function<void(const unsigned
int, const unsigned
int)>
5505 &operation_after_loop,
5506 const unsigned
int dof_handler_index_pre_post,
5507 const DataAccessOnFaces dst_vector_face_access,
5508 const DataAccessOnFaces src_vector_face_access) const
5510 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5521 inner_face_operation,
5522 boundary_face_operation,
5523 src_vector_face_access,
5524 dst_vector_face_access,
5525 operation_before_loop,
5526 operation_after_loop,
5527 dof_handler_index_pre_post);
5533template <
int dim,
typename Number,
typename VectorizedArrayType>
5534template <
typename CLASS,
typename OutVector,
typename InVector>
5537 void (CLASS::*function_pointer)(
5538 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5541 const
std::pair<unsigned
int, unsigned
int> &) const,
5542 const CLASS *owning_class,
5544 const InVector &src,
5545 const
bool zero_dst_vector,
5546 const DataAccessOnFaces src_vector_face_access) const
5548 auto src_vector_face_access_temp = src_vector_face_access;
5554 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5567 src_vector_face_access_temp,
5574template <
int dim,
typename Number,
typename VectorizedArrayType>
5575template <
typename CLASS,
typename OutVector,
typename InVector>
5578 void (CLASS::*function_pointer)(
5579 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5582 const
std::pair<unsigned
int, unsigned
int> &),
5583 CLASS *owning_class,
5585 const InVector &src,
5586 const
bool zero_dst_vector,
5587 const DataAccessOnFaces src_vector_face_access) const
5589 auto src_vector_face_access_temp = src_vector_face_access;
5595 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5608 src_vector_face_access_temp,
5615template <
int dim,
typename Number,
typename VectorizedArrayType>
5616template <
typename OutVector,
typename InVector>
5622 const std::pair<unsigned int, unsigned int> &)>
5625 const InVector &src,
5626 const bool zero_dst_vector,
5627 const DataAccessOnFaces src_vector_face_access)
const
5629 auto src_vector_face_access_temp = src_vector_face_access;
5630 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5631 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5632 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5633 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5636 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5639 Wrapper wrap(cell_operation,
nullptr,
nullptr);
5641 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5651 &Wrapper::cell_integrator,
5652 &Wrapper::face_integrator,
5653 &Wrapper::boundary_integrator,
5654 src_vector_face_access_temp,
5655 DataAccessOnFaces::none);
5656 task_info.loop(worker);
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.
unsigned int n_active_fe_indices() const
std::pair< typename DoFHandler< dim >::cell_iterator, unsigned int > get_face_iterator(const unsigned int face_batch_index, const unsigned int lane_index, const bool interior=true, const unsigned int fe_component=0) const
void loop_cell_centric(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > mapping_info
unsigned int get_cell_active_fe_index(const std::pair< unsigned int, unsigned int > range) const
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
unsigned int n_ghost_cell_batches() const
void 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
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 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
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
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
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
std::vector< ObserverPointer< const AffineConstraints< Number > > > affine_constraints
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
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
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
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 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
VectorType::value_type * begin(VectorType &V)
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