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.");
587 template <
typename QuadratureType,
typename number2,
typename MappingType>
592 const QuadratureType &quad,
616 template <
typename QuadratureType,
typename number2,
typename MappingType>
621 const std::vector<QuadratureType> &quad,
631 template <
typename QuadratureType,
typename number2,
typename MappingType>
636 const QuadratureType &quad,
793 template <
typename OutVector,
typename InVector>
799 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
802 const bool zero_dst_vector =
false)
const;
850 template <
typename CLASS,
typename OutVector,
typename InVector>
856 const
std::pair<unsigned
int, unsigned
int> &) const,
857 const CLASS *owning_class,
860 const
bool zero_dst_vector = false) const;
865 template <
typename CLASS,
typename OutVector,
typename InVector>
871 const
std::pair<unsigned
int, unsigned
int> &),
875 const
bool zero_dst_vector = false) const;
961 template <
typename CLASS,
typename OutVector,
typename InVector>
967 const
std::pair<unsigned
int, unsigned
int> &) const,
968 const CLASS *owning_class,
971 const
std::function<void(const unsigned
int, const unsigned
int)>
972 &operation_before_loop,
973 const
std::function<void(const unsigned
int, const unsigned
int)>
974 &operation_after_loop,
975 const unsigned
int dof_handler_index_pre_post = 0) const;
980 template <
typename CLASS,
typename OutVector,
typename InVector>
986 const
std::pair<unsigned
int, unsigned
int> &),
990 const
std::function<void(const unsigned
int, const unsigned
int)>
991 &operation_before_loop,
992 const
std::function<void(const unsigned
int, const unsigned
int)>
993 &operation_after_loop,
994 const unsigned
int dof_handler_index_pre_post = 0) const;
1000 template <
typename OutVector,
typename InVector>
1006 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1008 const InVector &src,
1009 const std::function<
void(
const unsigned int,
const unsigned int)>
1010 &operation_before_loop,
1011 const std::function<
void(
const unsigned int,
const unsigned int)>
1012 &operation_after_loop,
1013 const unsigned int dof_handler_index_pre_post = 0)
const;
1090 template <
typename OutVector,
typename InVector>
1096 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1097 const std::function<
1101 const std::pair<unsigned int, unsigned int> &)> &face_operation,
1102 const std::function<
void(
1106 const std::pair<unsigned int, unsigned int> &)> &boundary_operation,
1108 const InVector &src,
1109 const bool zero_dst_vector =
false,
1200 template <
typename CLASS,
typename OutVector,
typename InVector>
1203 void (CLASS::*cell_operation)(const
MatrixFree &,
1206 const
std::pair<unsigned
int, unsigned
int> &)
1208 void (CLASS::*face_operation)(const
MatrixFree &,
1211 const
std::pair<unsigned
int, unsigned
int> &)
1213 void (CLASS::*boundary_operation)(
1217 const
std::pair<unsigned
int, unsigned
int> &) const,
1218 const CLASS *owning_class,
1220 const InVector &src,
1221 const
bool zero_dst_vector = false,
1230 template <
typename CLASS,
typename OutVector,
typename InVector>
1232 loop(
void (CLASS::*cell_operation)(
1236 const
std::pair<unsigned
int, unsigned
int> &),
1237 void (CLASS::*face_operation)(
1241 const
std::pair<unsigned
int, unsigned
int> &),
1242 void (CLASS::*boundary_operation)(
1246 const
std::pair<unsigned
int, unsigned
int> &),
1247 CLASS *owning_class,
1249 const InVector &src,
1250 const
bool zero_dst_vector = false,
1375 template <
typename CLASS,
typename OutVector,
typename InVector>
1378 void (CLASS::*cell_operation)(const
MatrixFree &,
1381 const
std::pair<unsigned
int, unsigned
int> &)
1383 void (CLASS::*face_operation)(const
MatrixFree &,
1386 const
std::pair<unsigned
int, unsigned
int> &)
1388 void (CLASS::*boundary_operation)(
1392 const
std::pair<unsigned
int, unsigned
int> &) const,
1393 const CLASS *owning_class,
1395 const InVector &src,
1396 const
std::function<void(const unsigned
int, const unsigned
int)>
1397 &operation_before_loop,
1398 const
std::function<void(const unsigned
int, const unsigned
int)>
1399 &operation_after_loop,
1400 const unsigned
int dof_handler_index_pre_post = 0,
1409 template <
typename CLASS,
typename OutVector,
typename InVector>
1411 loop(
void (CLASS::*cell_operation)(
1415 const
std::pair<unsigned
int, unsigned
int> &),
1416 void (CLASS::*face_operation)(
1420 const
std::pair<unsigned
int, unsigned
int> &),
1421 void (CLASS::*boundary_operation)(
1425 const
std::pair<unsigned
int, unsigned
int> &),
1426 const CLASS *owning_class,
1428 const InVector &src,
1429 const
std::function<void(const unsigned
int, const unsigned
int)>
1430 &operation_before_loop,
1431 const
std::function<void(const unsigned
int, const unsigned
int)>
1432 &operation_after_loop,
1433 const unsigned
int dof_handler_index_pre_post = 0,
1444 template <
typename OutVector,
typename InVector>
1450 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1451 const std::function<
1455 const std::pair<unsigned int, unsigned int> &)> &face_operation,
1456 const std::function<
void(
1460 const std::pair<unsigned int, unsigned int> &)> &boundary_operation,
1462 const InVector &src,
1463 const std::function<
void(
const unsigned int,
const unsigned int)>
1464 &operation_before_loop,
1465 const std::function<
void(
const unsigned int,
const unsigned int)>
1466 &operation_after_loop,
1467 const unsigned int dof_handler_index_pre_post = 0,
1537 template <
typename CLASS,
typename OutVector,
typename InVector>
1543 const
std::pair<unsigned
int, unsigned
int> &) const,
1544 const CLASS *owning_class,
1546 const InVector &src,
1547 const
bool zero_dst_vector = false,
1554 template <
typename CLASS,
typename OutVector,
typename InVector>
1560 const
std::pair<unsigned
int, unsigned
int> &),
1561 CLASS *owning_class,
1563 const InVector &src,
1564 const
bool zero_dst_vector = false,
1571 template <
typename OutVector,
typename InVector>
1577 const std::pair<unsigned int, unsigned int> &)>
1580 const InVector &src,
1581 const bool zero_dst_vector =
false,
1592 std::pair<unsigned int, unsigned int>
1594 const unsigned int fe_degree,
1595 const unsigned int dof_handler_index = 0)
const;
1603 std::pair<unsigned int, unsigned int>
1605 const std::pair<unsigned int, unsigned int> &range,
1606 const unsigned int fe_index,
1607 const unsigned int dof_handler_index = 0)
const;
1620 const std::pair<unsigned int, unsigned int> range)
const;
1627 const bool is_interior_face =
true)
const;
1640 template <
typename T>
1649 template <
typename T>
1666 template <
typename VectorType>
1669 const unsigned int dof_handler_index = 0)
const;
1688 template <
typename Number2>
1691 const unsigned int dof_handler_index = 0)
const;
1703 const std::shared_ptr<const Utilities::MPI::Partitioner> &
1727 const std::vector<unsigned int> &
1742 const unsigned int dof_handler_index = 0);
1753 template <
int spacedim>
1852 const unsigned int face_number)
const;
1884 const unsigned int lane_index,
1885 const unsigned int dof_handler_index = 0)
const;
1894 const unsigned int lane_index)
const;
1918 std::pair<typename DoFHandler<dim>::cell_iterator,
unsigned int>
1920 const unsigned int lane_index,
1921 const bool interior =
true,
1922 const unsigned int fe_component = 0)
const;
1968 const unsigned int hp_active_fe_index = 0)
const;
1975 const unsigned int hp_active_fe_index = 0)
const;
1983 const unsigned int hp_active_fe_index = 0)
const;
1991 const unsigned int hp_active_fe_index = 0)
const;
1998 const unsigned int hp_active_fe_index = 0)
const;
2005 const unsigned int hp_active_fe_index = 0)
const;
2021 const std::pair<unsigned int, unsigned int> cell_batch_range)
const;
2027 std::pair<unsigned int, unsigned int>
2029 const std::pair<unsigned int, unsigned int> face_batch_range)
const;
2049 std::pair<unsigned int, unsigned int>
2083 template <
typename StreamType>
2111 const internal::MatrixFreeFunctions::
2112 MappingInfo<dim, Number, VectorizedArrayType> &
2147 const unsigned int quad_index = 0,
2148 const unsigned int fe_base_element = 0,
2149 const unsigned int hp_active_fe_index = 0,
2150 const unsigned int hp_active_quad_index = 0)
const;
2156 VectorizedArrayType::size()> &
2216 template <
typename number2,
int q_dim>
2222 const std::vector<IndexSet> &locally_owned_set,
2232 template <
typename number2>
2236 const std::vector<IndexSet> &locally_owned_set,
2264 std::vector<internal::MatrixFreeFunctions::DoFInfo>
dof_info;
2347 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2354 mutable std::list<std::pair<bool, AlignedVector<Number>>>
2371template <
int dim,
typename Number,
typename VectorizedArrayType>
2372template <
typename T>
2377 vec.
resize(this->n_cell_batches() + this->n_ghost_cell_batches());
2382template <
int dim,
typename Number,
typename VectorizedArrayType>
2383template <
typename T>
2388 vec.
resize(this->n_inner_face_batches() + this->n_boundary_face_batches() +
2389 this->n_ghost_inner_face_batches());
2394template <
int dim,
typename Number,
typename VectorizedArrayType>
2395template <
typename VectorType>
2399 const unsigned int comp)
const
2402 "This function is not supported for block vectors.");
2404 Assert(task_info.n_procs == 1,
2405 ExcMessage(
"This function can only be used in serial."));
2408 vec.reinit(dof_info[comp].vector_partitioner->size());
2413template <
int dim,
typename Number,
typename VectorizedArrayType>
2414template <
typename Number2>
2418 const unsigned int comp)
const
2421 vec.
reinit(dof_info[comp].vector_partitioner, task_info.communicator_sm);
2426template <
int dim,
typename Number,
typename VectorizedArrayType>
2427inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2429 const unsigned int comp)
const
2432 return dof_info[comp].vector_partitioner;
2437template <
int dim,
typename Number,
typename VectorizedArrayType>
2438inline const std::vector<unsigned int> &
2440 const unsigned int comp)
const
2443 return dof_info[comp].constrained_dofs;
2448template <
int dim,
typename Number,
typename VectorizedArrayType>
2453 return dof_handlers.size();
2458template <
int dim,
typename Number,
typename VectorizedArrayType>
2461 const unsigned int dof_no)
const
2465 return dof_handlers[dof_no]->get_fe().n_base_elements();
2470template <
int dim,
typename Number,
typename VectorizedArrayType>
2479template <
int dim,
typename Number,
typename VectorizedArrayType>
2483 return task_info.n_active_cells;
2488template <
int dim,
typename Number,
typename VectorizedArrayType>
2492 return *(task_info.cell_partition_data.end() - 2);
2497template <
int dim,
typename Number,
typename VectorizedArrayType>
2501 return *(task_info.cell_partition_data.end() - 1) -
2502 *(task_info.cell_partition_data.end() - 2);
2507template <
int dim,
typename Number,
typename VectorizedArrayType>
2511 if (task_info.face_partition_data.empty())
2513 return task_info.face_partition_data.back();
2518template <
int dim,
typename Number,
typename VectorizedArrayType>
2522 if (task_info.face_partition_data.empty())
2524 return task_info.boundary_partition_data.back() -
2525 task_info.face_partition_data.back();
2530template <
int dim,
typename Number,
typename VectorizedArrayType>
2534 if (task_info.face_partition_data.empty())
2536 return face_info.faces.size() - task_info.boundary_partition_data.back();
2541template <
int dim,
typename Number,
typename VectorizedArrayType>
2544 const unsigned int face_batch_index)
const
2546 Assert(face_batch_index >= task_info.boundary_partition_data[0] &&
2547 face_batch_index < task_info.boundary_partition_data.back(),
2549 task_info.boundary_partition_data[0],
2550 task_info.boundary_partition_data.back()));
2556template <
int dim,
typename Number,
typename VectorizedArrayType>
2559 const unsigned int cell_batch_index,
2560 const unsigned int face_number)
const
2564 Assert(face_info.cell_and_face_boundary_id.size(0) >= n_cell_batches(),
2568 for (
unsigned int v = 0;
2569 v < n_active_entries_per_cell_batch(cell_batch_index);
2572 face_info.cell_and_face_boundary_id(cell_batch_index, face_number, v);
2578template <
int dim,
typename Number,
typename VectorizedArrayType>
2579inline const internal::MatrixFreeFunctions::
2580 MappingInfo<dim, Number, VectorizedArrayType> &
2583 return mapping_info;
2588template <
int dim,
typename Number,
typename VectorizedArrayType>
2591 const unsigned int dof_index)
const
2594 return dof_info[dof_index];
2599template <
int dim,
typename Number,
typename VectorizedArrayType>
2603 return constraint_pool_row_index.size() - 1;
2608template <
int dim,
typename Number,
typename VectorizedArrayType>
2609inline const Number *
2611 const unsigned int row)
const
2614 return constraint_pool_data.empty() ?
2616 constraint_pool_data.data() + constraint_pool_row_index[row];
2621template <
int dim,
typename Number,
typename VectorizedArrayType>
2622inline const Number *
2624 const unsigned int row)
const
2627 return constraint_pool_data.empty() ?
2629 constraint_pool_data.data() + constraint_pool_row_index[row + 1];
2634template <
int dim,
typename Number,
typename VectorizedArrayType>
2635inline std::pair<unsigned int, unsigned int>
2637 const std::pair<unsigned int, unsigned int> &range,
2638 const unsigned int degree,
2639 const unsigned int dof_handler_component)
const
2641 if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2644 dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2646 dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2647 if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2650 return {range.second, range.second};
2654 dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2655 if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2656 return {range.second, range.second};
2658 return create_cell_subrange_hp_by_index(range,
2660 dof_handler_component);
2665template <
int dim,
typename Number,
typename VectorizedArrayType>
2668 const unsigned int cell_batch_index)
const
2671 return VectorizedArrayType::size() > 1 &&
2672 cell_level_index[(cell_batch_index + 1) * VectorizedArrayType::size() -
2673 1] == cell_level_index[(cell_batch_index + 1) *
2674 VectorizedArrayType::size() -
2680template <
int dim,
typename Number,
typename VectorizedArrayType>
2684 return shape_info.size(2);
2688template <
int dim,
typename Number,
typename VectorizedArrayType>
2691 const std::pair<unsigned int, unsigned int> range)
const
2693 const auto &fe_indices = dof_info[0].cell_active_fe_index;
2695 if (fe_indices.empty() ==
true ||
2696 dof_handlers[0]->get_fe_collection().size() == 1)
2699 const auto index = fe_indices[range.first];
2701 for (
unsigned int i = range.first; i < range.second; ++i)
2709template <
int dim,
typename Number,
typename VectorizedArrayType>
2712 const std::pair<unsigned int, unsigned int> range,
2713 const bool is_interior_face)
const
2715 const auto &fe_indices = dof_info[0].cell_active_fe_index;
2717 if (fe_indices.empty() ==
true)
2720 if (is_interior_face)
2722 const unsigned int index =
2723 fe_indices[face_info.faces[range.first].cells_interior[0] /
2724 VectorizedArrayType::size()];
2726 for (
unsigned int i = range.first; i < range.second; ++i)
2728 fe_indices[face_info.faces[i].cells_interior[0] /
2729 VectorizedArrayType::size()]);
2735 const unsigned int index =
2736 fe_indices[face_info.faces[range.first].cells_exterior[0] /
2737 VectorizedArrayType::size()];
2739 for (
unsigned int i = range.first; i < range.second; ++i)
2741 fe_indices[face_info.faces[i].cells_exterior[0] /
2742 VectorizedArrayType::size()]);
2750template <
int dim,
typename Number,
typename VectorizedArrayType>
2753 const unsigned int cell_batch_index)
const
2757 const std::vector<unsigned char> &n_lanes_filled =
2758 dof_info[0].n_vectorization_lanes_filled
2762 return n_lanes_filled[cell_batch_index];
2767template <
int dim,
typename Number,
typename VectorizedArrayType>
2770 const unsigned int face_batch_index)
const
2774 const std::vector<unsigned char> &n_lanes_filled =
2775 dof_info[0].n_vectorization_lanes_filled
2778 return n_lanes_filled[face_batch_index];
2783template <
int dim,
typename Number,
typename VectorizedArrayType>
2786 const unsigned int dof_handler_index,
2787 const unsigned int active_fe_index)
const
2789 return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2794template <
int dim,
typename Number,
typename VectorizedArrayType>
2797 const unsigned int quad_index,
2798 const unsigned int active_fe_index)
const
2801 return mapping_info.cell_data[quad_index]
2802 .descriptor[active_fe_index]
2808template <
int dim,
typename Number,
typename VectorizedArrayType>
2811 const unsigned int dof_handler_index,
2812 const unsigned int active_fe_index)
const
2814 return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2819template <
int dim,
typename Number,
typename VectorizedArrayType>
2822 const unsigned int quad_index,
2823 const unsigned int active_fe_index)
const
2826 return mapping_info.face_data[quad_index]
2827 .descriptor[active_fe_index]
2833template <
int dim,
typename Number,
typename VectorizedArrayType>
2836 const unsigned int dof_handler_index)
const
2838 return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2843template <
int dim,
typename Number,
typename VectorizedArrayType>
2846 const unsigned int dof_handler_index)
const
2848 return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2853template <
int dim,
typename Number,
typename VectorizedArrayType>
2856 const unsigned int dof_handler_index,
2857 const unsigned int index_quad,
2858 const unsigned int index_fe,
2859 const unsigned int active_fe_index,
2860 const unsigned int active_quad_index)
const
2863 const unsigned int ind =
2864 dof_info[dof_handler_index].global_base_element_offset + index_fe;
2869 return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2874template <
int dim,
typename Number,
typename VectorizedArrayType>
2876 VectorizedArrayType::size()> &
2878 const unsigned int face_batch_index)
const
2881 return face_info.faces[face_batch_index];
2886template <
int dim,
typename Number,
typename VectorizedArrayType>
2891 return face_info.cell_and_face_to_plain_faces;
2896template <
int dim,
typename Number,
typename VectorizedArrayType>
2899 const unsigned int quad_index,
2900 const unsigned int active_fe_index)
const
2903 return mapping_info.cell_data[quad_index]
2904 .descriptor[active_fe_index]
2910template <
int dim,
typename Number,
typename VectorizedArrayType>
2913 const unsigned int quad_index,
2914 const unsigned int active_fe_index)
const
2917 return mapping_info.face_data[quad_index]
2918 .descriptor[active_fe_index]
2924template <
int dim,
typename Number,
typename VectorizedArrayType>
2927 const std::pair<unsigned int, unsigned int> range)
const
2929 auto result = get_cell_category(range.first);
2931 for (
unsigned int i = range.first; i < range.second; ++i)
2932 result =
std::max(result, get_cell_category(i));
2939template <
int dim,
typename Number,
typename VectorizedArrayType>
2940inline std::pair<unsigned int, unsigned int>
2942 const std::pair<unsigned int, unsigned int> range)
const
2944 auto result = get_face_category(range.first);
2946 for (
unsigned int i = range.first; i < range.second; ++i)
2948 result.first =
std::max(result.first, get_face_category(i).
first);
2949 result.second =
std::max(result.second, get_face_category(i).
second);
2957template <
int dim,
typename Number,
typename VectorizedArrayType>
2960 const unsigned int cell_batch_index)
const
2963 AssertIndexRange(cell_batch_index, dof_info[0].cell_active_fe_index.size());
2964 if (dof_info[0].cell_active_fe_index.empty())
2967 return dof_info[0].cell_active_fe_index[cell_batch_index];
2972template <
int dim,
typename Number,
typename VectorizedArrayType>
2973inline std::pair<unsigned int, unsigned int>
2975 const unsigned int face_batch_index)
const
2978 if (dof_info[0].cell_active_fe_index.empty())
2979 return std::make_pair(0U, 0U);
2981 std::pair<unsigned int, unsigned int> result = std::make_pair(0U, 0U);
2982 for (
unsigned int v = 0;
2983 v < VectorizedArrayType::size() &&
2984 face_info.faces[face_batch_index].cells_interior[v] !=
2989 dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
2990 .cells_interior[v] /
2991 VectorizedArrayType::size()]);
2992 if (face_info.faces[face_batch_index].cells_exterior[0] !=
2994 for (
unsigned int v = 0;
2995 v < VectorizedArrayType::size() &&
2996 face_info.faces[face_batch_index].cells_exterior[v] !=
3001 dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
3002 .cells_exterior[v] /
3003 VectorizedArrayType::size()]);
3011template <
int dim,
typename Number,
typename VectorizedArrayType>
3015 return indices_are_initialized;
3020template <
int dim,
typename Number,
typename VectorizedArrayType>
3024 return mapping_is_initialized;
3028template <
int dim,
typename Number,
typename VectorizedArrayType>
3037template <
int dim,
typename Number,
typename VectorizedArrayType>
3042 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3043 list_type &data = scratch_pad.get();
3044 for (
typename list_type::iterator it = data.begin();
it != data.end(); ++
it)
3045 if (
it->first ==
false)
3051 return &data.front().second;
3056template <
int dim,
typename Number,
typename VectorizedArrayType>
3062 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3063 list_type &data = scratch_pad.get();
3064 for (
typename list_type::iterator it = data.begin();
it != data.end(); ++
it)
3065 if (&
it->second == scratch)
3076template <
int dim,
typename Number,
typename VectorizedArrayType>
3082 scratch_pad_non_threadsafe.
begin();
3083 it != scratch_pad_non_threadsafe.end();
3085 if (
it->first ==
false)
3090 scratch_pad_non_threadsafe.push_front(
3092 return &scratch_pad_non_threadsafe.front().second;
3097template <
int dim,
typename Number,
typename VectorizedArrayType>
3104 scratch_pad_non_threadsafe.
begin();
3105 it != scratch_pad_non_threadsafe.end();
3107 if (&
it->second == scratch)
3122 namespace MatrixFreeImplementation
3124 template <
int dim,
int spacedim>
3125 inline std::vector<IndexSet>
3126 extract_locally_owned_index_sets(
3128 const unsigned int level)
3130 std::vector<IndexSet> locally_owned_set;
3131 locally_owned_set.reserve(dofh.size());
3132 for (
unsigned int j = 0; j < dofh.size(); ++j)
3134 locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
3136 locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(
level));
3137 return locally_owned_set;
3144template <
int dim,
typename Number,
typename VectorizedArrayType>
3145template <
typename QuadratureType,
typename number2,
typename MappingType>
3148 const MappingType &mapping,
3151 const QuadratureType &quad,
3155 std::vector<const DoFHandler<dim, dim> *> dof_handlers;
3156 std::vector<const AffineConstraints<number2> *> constraints;
3158 dof_handlers.push_back(&dof_handler);
3159 constraints.push_back(&constraints_in);
3161 std::vector<IndexSet> locally_owned_sets =
3162 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3163 dof_handlers, additional_data.
mg_level);
3165 std::vector<hp::QCollection<dim>> quad_hp;
3166 quad_hp.emplace_back(quad);
3178template <
int dim,
typename Number,
typename VectorizedArrayType>
3179template <
typename QuadratureType,
typename number2,
typename MappingType>
3182 const MappingType &mapping,
3185 const QuadratureType &quad,
3189 std::vector<IndexSet> locally_owned_set =
3190 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3191 dof_handler, additional_data.
mg_level);
3192 std::vector<hp::QCollection<dim>> quad_hp;
3193 quad_hp.emplace_back(quad);
3205template <
int dim,
typename Number,
typename VectorizedArrayType>
3206template <
typename QuadratureType,
typename number2,
typename MappingType>
3209 const MappingType &mapping,
3212 const std::vector<QuadratureType> &quad,
3216 std::vector<IndexSet> locally_owned_set =
3217 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3218 dof_handler, additional_data.
mg_level);
3219 std::vector<hp::QCollection<dim>> quad_hp;
3220 for (
unsigned int q = 0; q < quad.size(); ++q)
3221 quad_hp.emplace_back(quad[q]);
3248 template <
int dim,
typename Number,
typename VectorizedArrayType>
3249 struct VectorDataExchange
3256 static constexpr unsigned int channel_shift = 20;
3265 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3266 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3267 DataAccessOnFaces vector_face_access,
3268 const unsigned
int n_components)
3269 : matrix_free(matrix_free)
3270 , vector_face_access(
3271 matrix_free.get_task_info().face_partition_data.empty() ?
3272 ::
MatrixFree<dim, Number, VectorizedArrayType>::
3273 DataAccessOnFaces::unspecified :
3275 , ghosts_were_set(false)
3276# ifdef DEAL_II_WITH_MPI
3277 , tmp_data(n_components)
3278 , requests(n_components)
3282 if (this->vector_face_access !=
3284 DataAccessOnFaces::unspecified)
3285 for (unsigned
int c = 0; c < matrix_free.n_components(); ++c)
3287 matrix_free.get_dof_info(c).vector_exchanger_face_variants.size(),
3296 ~VectorDataExchange()
3298# ifdef DEAL_II_WITH_MPI
3299 for (
unsigned int i = 0; i < tmp_data.size(); ++i)
3300 if (tmp_data[i] !=
nullptr)
3301 matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3311 template <
typename VectorType>
3313 find_vector_in_mf(
const VectorType &vec,
3314 const bool check_global_compatibility =
true)
const
3317 for (
unsigned int c = 0; c < matrix_free.n_components(); ++c)
3318 if (vec.get_partitioner().get() ==
3319 matrix_free.get_dof_info(c).vector_partitioner.get())
3323 for (
unsigned int c = 0; c < matrix_free.n_components(); ++c)
3324 if (check_global_compatibility ?
3325 vec.get_partitioner()->is_globally_compatible(
3326 *matrix_free.get_dof_info(c).vector_partitioner) :
3327 vec.get_partitioner()->is_compatible(
3328 *matrix_free.get_dof_info(c).vector_partitioner))
3344 get_partitioner(
const unsigned int mf_component)
const
3347 .vector_exchanger_face_variants.size(),
3349 if (vector_face_access ==
3351 DataAccessOnFaces::
none)
3352 return *matrix_free.get_dof_info(mf_component)
3353 .vector_exchanger_face_variants[0];
3354 else if (vector_face_access ==
3356 DataAccessOnFaces::
values)
3357 return *matrix_free.get_dof_info(mf_component)
3358 .vector_exchanger_face_variants[1];
3359 else if (vector_face_access ==
3362 return *matrix_free.get_dof_info(mf_component)
3363 .vector_exchanger_face_variants[2];
3364 else if (vector_face_access ==
3366 DataAccessOnFaces::values_all_faces)
3367 return *matrix_free.get_dof_info(mf_component)
3368 .vector_exchanger_face_variants[3];
3369 else if (vector_face_access ==
3371 DataAccessOnFaces::gradients_all_faces)
3372 return *matrix_free.get_dof_info(mf_component)
3373 .vector_exchanger_face_variants[4];
3375 return *matrix_free.get_dof_info(mf_component).vector_exchanger.get();
3383 template <
typename VectorType,
3384 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3387 update_ghost_values_start(
const unsigned int ,
3388 const VectorType & )
3396 template <
typename VectorType,
3397 std::enable_if_t<!has_update_ghost_values_start<VectorType> &&
3398 !is_not_parallel_vector<VectorType>,
3399 VectorType> * =
nullptr>
3401 update_ghost_values_start(
const unsigned int component_in_block_vector,
3402 const VectorType &vec)
3404 (void)component_in_block_vector;
3405 const bool ghosts_set = vec.has_ghost_elements();
3407 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3408 ghosts_set ==
false,
3412 ghosts_were_set =
true;
3414 vec.update_ghost_values();
3424 template <
typename VectorType,
3425 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3426 !has_exchange_on_subset<VectorType>,
3427 VectorType> * =
nullptr>
3429 update_ghost_values_start(
const unsigned int component_in_block_vector,
3430 const VectorType &vec)
3432 (void)component_in_block_vector;
3433 const bool ghosts_set = vec.has_ghost_elements();
3435 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3436 ghosts_set ==
false,
3440 ghosts_were_set =
true;
3442 vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3453 template <
typename VectorType,
3454 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3455 has_exchange_on_subset<VectorType>,
3456 VectorType> * =
nullptr>
3458 update_ghost_values_start(
const unsigned int component_in_block_vector,
3459 const VectorType &vec)
3461 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3462 "Type mismatch between VectorType and VectorDataExchange");
3463 (void)component_in_block_vector;
3464 const bool ghosts_set = vec.has_ghost_elements();
3466 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3467 ghosts_set ==
false,
3471 ghosts_were_set =
true;
3473 if (vec.size() != 0)
3475# ifdef DEAL_II_WITH_MPI
3476 const unsigned int mf_component = find_vector_in_mf(vec);
3478 const auto &part = get_partitioner(mf_component);
3480 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3481 part.n_import_sm_procs() == 0)
3484 tmp_data[component_in_block_vector] =
3485 matrix_free.acquire_scratch_data_non_threadsafe();
3486 tmp_data[component_in_block_vector]->resize_fast(
3487 part.n_import_indices());
3490 part.export_to_ghosted_array_start(
3491 component_in_block_vector * 2 + channel_shift,
3493 vec.shared_vector_data(),
3495 part.locally_owned_size(),
3496 matrix_free.get_dof_info(mf_component)
3497 .vector_partitioner->n_ghost_indices()),
3499 part.n_import_indices()),
3500 this->requests[component_in_block_vector]);
3511 template <
typename VectorType,
3512 std::enable_if_t<!has_update_ghost_values_start<VectorType>,
3513 VectorType> * =
nullptr>
3515 update_ghost_values_finish(
const unsigned int ,
3516 const VectorType & )
3526 template <
typename VectorType,
3527 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3528 !has_exchange_on_subset<VectorType>,
3529 VectorType> * =
nullptr>
3531 update_ghost_values_finish(
const unsigned int component_in_block_vector,
3532 const VectorType &vec)
3534 (void)component_in_block_vector;
3535 vec.update_ghost_values_finish();
3546 template <
typename VectorType,
3547 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3548 has_exchange_on_subset<VectorType>,
3549 VectorType> * =
nullptr>
3551 update_ghost_values_finish(
const unsigned int component_in_block_vector,
3552 const VectorType &vec)
3554 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3555 "Type mismatch between VectorType and VectorDataExchange");
3556 (void)component_in_block_vector;
3558 if (vec.size() != 0)
3560# ifdef DEAL_II_WITH_MPI
3564 const unsigned int mf_component = find_vector_in_mf(vec);
3566 const auto &part = get_partitioner(mf_component);
3568 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3569 part.n_import_sm_procs() != 0)
3571 part.export_to_ghosted_array_finish(
3573 vec.shared_vector_data(),
3575 part.locally_owned_size(),
3576 matrix_free.get_dof_info(mf_component)
3577 .vector_partitioner->n_ghost_indices()),
3578 this->requests[component_in_block_vector]);
3580 matrix_free.release_scratch_data_non_threadsafe(
3581 tmp_data[component_in_block_vector]);
3582 tmp_data[component_in_block_vector] =
nullptr;
3588 vec.set_ghost_state(
true);
3596 template <
typename VectorType,
3597 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3600 compress_start(
const unsigned int ,
3610 template <
typename VectorType,
3611 std::enable_if_t<!has_compress_start<VectorType> &&
3612 !is_not_parallel_vector<VectorType>,
3613 VectorType> * =
nullptr>
3615 compress_start(
const unsigned int component_in_block_vector,
3618 (void)component_in_block_vector;
3630 template <
typename VectorType,
3631 std::enable_if_t<has_compress_start<VectorType> &&
3632 !has_exchange_on_subset<VectorType>,
3633 VectorType> * =
nullptr>
3635 compress_start(
const unsigned int component_in_block_vector,
3638 (void)component_in_block_vector;
3640 vec.compress_start(component_in_block_vector + channel_shift);
3651 template <
typename VectorType,
3652 std::enable_if_t<has_compress_start<VectorType> &&
3653 has_exchange_on_subset<VectorType>,
3654 VectorType> * =
nullptr>
3656 compress_start(
const unsigned int component_in_block_vector,
3659 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3660 "Type mismatch between VectorType and VectorDataExchange");
3661 (void)component_in_block_vector;
3664 if (vec.size() != 0)
3666# ifdef DEAL_II_WITH_MPI
3667 const unsigned int mf_component = find_vector_in_mf(vec);
3669 const auto &part = get_partitioner(mf_component);
3671 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3672 part.n_import_sm_procs() == 0)
3675 tmp_data[component_in_block_vector] =
3676 matrix_free.acquire_scratch_data_non_threadsafe();
3677 tmp_data[component_in_block_vector]->resize_fast(
3678 part.n_import_indices());
3681 part.import_from_ghosted_array_start(
3683 component_in_block_vector * 2 + channel_shift,
3685 vec.shared_vector_data(),
3687 matrix_free.get_dof_info(mf_component)
3688 .vector_partitioner->n_ghost_indices()),
3690 part.n_import_indices()),
3691 this->requests[component_in_block_vector]);
3703 typename VectorType,
3704 std::enable_if_t<!has_compress_start<VectorType>, VectorType> * =
nullptr>
3706 compress_finish(
const unsigned int ,
3717 template <
typename VectorType,
3718 std::enable_if_t<has_compress_start<VectorType> &&
3719 !has_exchange_on_subset<VectorType>,
3720 VectorType> * =
nullptr>
3722 compress_finish(
const unsigned int component_in_block_vector,
3725 (void)component_in_block_vector;
3737 template <
typename VectorType,
3738 std::enable_if_t<has_compress_start<VectorType> &&
3739 has_exchange_on_subset<VectorType>,
3740 VectorType> * =
nullptr>
3742 compress_finish(
const unsigned int component_in_block_vector,
3745 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3746 "Type mismatch between VectorType and VectorDataExchange");
3747 (void)component_in_block_vector;
3748 if (vec.size() != 0)
3750# ifdef DEAL_II_WITH_MPI
3754 const unsigned int mf_component = find_vector_in_mf(vec);
3756 const auto &part = get_partitioner(mf_component);
3758 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3759 part.n_import_sm_procs() != 0)
3761 part.import_from_ghosted_array_finish(
3764 vec.shared_vector_data(),
3766 matrix_free.get_dof_info(mf_component)
3767 .vector_partitioner->n_ghost_indices()),
3769 tmp_data[component_in_block_vector]->begin(),
3770 part.n_import_indices()),
3771 this->requests[component_in_block_vector]);
3773 matrix_free.release_scratch_data_non_threadsafe(
3774 tmp_data[component_in_block_vector]);
3775 tmp_data[component_in_block_vector] =
nullptr;
3781 MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3793 template <
typename VectorType,
3794 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3797 reset_ghost_values(
const VectorType & )
const
3806 template <
typename VectorType,
3807 std::enable_if_t<!has_exchange_on_subset<VectorType> &&
3808 !is_not_parallel_vector<VectorType>,
3809 VectorType> * =
nullptr>
3811 reset_ghost_values(
const VectorType &vec)
const
3813 if (ghosts_were_set ==
true)
3816 vec.zero_out_ghost_values();
3826 template <
typename VectorType,
3827 std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3830 reset_ghost_values(
const VectorType &vec)
const
3832 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3833 "Type mismatch between VectorType and VectorDataExchange");
3834 if (ghosts_were_set ==
true)
3837 if (vec.size() != 0)
3839# ifdef DEAL_II_WITH_MPI
3842 const unsigned int mf_component = find_vector_in_mf(vec);
3844 const auto &part = get_partitioner(mf_component);
3846 if (part.n_ghost_indices() > 0)
3851 part.locally_owned_size(),
3852 matrix_free.get_dof_info(mf_component)
3853 .vector_partitioner->n_ghost_indices()));
3859 vec.set_ghost_state(
false);
3869 template <
typename VectorType,
3870 std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3873 zero_vector_region(
const unsigned int range_index, VectorType &vec)
const
3875 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3876 "Type mismatch between VectorType and VectorDataExchange");
3881 const unsigned int mf_component = find_vector_in_mf(vec,
false);
3883 matrix_free.get_dof_info(mf_component);
3891 for (
unsigned int id =
3910 template <
typename VectorType,
3911 std::enable_if_t<!has_exchange_on_subset<VectorType>, VectorType>
3913 typename VectorType::value_type * =
nullptr>
3915 zero_vector_region(
const unsigned int range_index, VectorType &vec)
const
3918 vec =
typename VectorType::value_type();
3928 zero_vector_region(
const unsigned int, ...)
const
3932 "which provide VectorType::value_type"));
3937 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
3938 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3939 DataAccessOnFaces vector_face_access;
3940 bool ghosts_were_set;
3941# ifdef DEAL_II_WITH_MPI
3942 std::vector<AlignedVector<Number> *> tmp_data;
3943 std::vector<std::vector<MPI_Request>> requests;
3947 template <
typename VectorStruct>
3949 n_components(
const VectorStruct &vec);
3951 template <
typename VectorStruct>
3953 n_components_block(
const VectorStruct &vec,
const std::bool_constant<true>)
3955 unsigned int components = 0;
3956 for (
unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
3957 components += n_components(vec.block(bl));
3961 template <
typename VectorStruct>
3963 n_components_block(
const VectorStruct &,
const std::bool_constant<false>)
3968 template <
typename VectorStruct>
3970 n_components(
const VectorStruct &vec)
3972 return n_components_block(
3976 template <
typename VectorStruct>
3978 n_components(
const std::vector<VectorStruct> &vec)
3980 unsigned int components = 0;
3981 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
3982 components += n_components_block(
3987 template <
typename VectorStruct>
3989 n_components(
const std::vector<VectorStruct *> &vec)
3991 unsigned int components = 0;
3992 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
3993 components += n_components_block(
4005 template <
typename VectorStruct,
4006 std::enable_if_t<!has_communication_block_size<VectorStruct>,
4007 VectorStruct> * =
nullptr>
4008 constexpr unsigned int
4009 get_communication_block_size(
const VectorStruct &)
4016 template <
typename VectorStruct,
4017 std::enable_if_t<has_communication_block_size<VectorStruct>,
4018 VectorStruct> * =
nullptr>
4019 constexpr unsigned int
4020 get_communication_block_size(
const VectorStruct &)
4022 return VectorStruct::communication_block_size;
4027 template <
typename VectorType,
4028 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType> * =
4031 has_ghost_elements(
const VectorType &vec)
4039 template <
typename VectorType,
4040 std::enable_if_t<!is_not_parallel_vector<VectorType>, VectorType>
4043 has_ghost_elements(
const VectorType &vec)
4045 return vec.has_ghost_elements();
4060 typename VectorStruct,
4062 typename VectorizedArrayType,
4063 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4066 update_ghost_values_start(
4067 const VectorStruct &vec,
4068 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4069 const unsigned int channel = 0)
4071 if (get_communication_block_size(vec) < vec.n_blocks())
4073 const bool ghosts_set = vec.has_ghost_elements();
4075 Assert(exchanger.matrix_free.get_task_info()
4076 .allow_ghosted_vectors_in_loops ||
4077 ghosts_set ==
false,
4081 exchanger.ghosts_were_set =
true;
4083 vec.update_ghost_values();
4087 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4088 update_ghost_values_start(vec.block(i), exchanger, channel + i);
4096 typename VectorStruct,
4098 typename VectorizedArrayType,
4099 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4102 update_ghost_values_start(
4103 const VectorStruct &vec,
4104 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4105 const unsigned int channel = 0)
4107 exchanger.update_ghost_values_start(channel, vec);
4114 typename VectorStruct,
4116 typename VectorizedArrayType>
4118 update_ghost_values_start(
4119 const std::vector<VectorStruct> &vec,
4120 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4122 unsigned int component_index = 0;
4123 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4125 update_ghost_values_start(vec[comp], exchanger, component_index);
4126 component_index += n_components(vec[comp]);
4134 typename VectorStruct,
4136 typename VectorizedArrayType>
4138 update_ghost_values_start(
4139 const std::vector<VectorStruct *> &vec,
4140 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4142 unsigned int component_index = 0;
4143 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4145 update_ghost_values_start(*vec[comp], exchanger, component_index);
4146 component_index += n_components(*vec[comp]);
4158 typename VectorStruct,
4160 typename VectorizedArrayType,
4161 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4164 update_ghost_values_finish(
4165 const VectorStruct &vec,
4166 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4167 const unsigned int channel = 0)
4169 if (get_communication_block_size(vec) < vec.n_blocks())
4175 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4176 update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4183 typename VectorStruct,
4185 typename VectorizedArrayType,
4186 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4189 update_ghost_values_finish(
4190 const VectorStruct &vec,
4191 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4192 const unsigned int channel = 0)
4194 exchanger.update_ghost_values_finish(channel, vec);
4201 typename VectorStruct,
4203 typename VectorizedArrayType>
4205 update_ghost_values_finish(
4206 const std::vector<VectorStruct> &vec,
4207 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4209 unsigned int component_index = 0;
4210 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4212 update_ghost_values_finish(vec[comp], exchanger, component_index);
4213 component_index += n_components(vec[comp]);
4221 typename VectorStruct,
4223 typename VectorizedArrayType>
4225 update_ghost_values_finish(
4226 const std::vector<VectorStruct *> &vec,
4227 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4229 unsigned int component_index = 0;
4230 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4232 update_ghost_values_finish(*vec[comp], exchanger, component_index);
4233 component_index += n_components(*vec[comp]);
4245 typename VectorStruct,
4247 typename VectorizedArrayType,
4248 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4253 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4254 const unsigned int channel = 0)
4256 if (get_communication_block_size(vec) < vec.n_blocks())
4259 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4260 compress_start(vec.block(i), exchanger, channel + i);
4267 typename VectorStruct,
4269 typename VectorizedArrayType,
4270 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4275 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4276 const unsigned int channel = 0)
4278 exchanger.compress_start(channel, vec);
4285 typename VectorStruct,
4287 typename VectorizedArrayType>
4290 std::vector<VectorStruct> &vec,
4291 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4293 unsigned int component_index = 0;
4294 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4296 compress_start(vec[comp], exchanger, component_index);
4297 component_index += n_components(vec[comp]);
4305 typename VectorStruct,
4307 typename VectorizedArrayType>
4310 std::vector<VectorStruct *> &vec,
4311 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4313 unsigned int component_index = 0;
4314 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4316 compress_start(*vec[comp], exchanger, component_index);
4317 component_index += n_components(*vec[comp]);
4329 typename VectorStruct,
4331 typename VectorizedArrayType,
4332 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4337 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4338 const unsigned int channel = 0)
4340 if (get_communication_block_size(vec) < vec.n_blocks())
4346 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4347 compress_finish(vec.block(i), exchanger, channel + i);
4354 typename VectorStruct,
4356 typename VectorizedArrayType,
4357 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4362 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4363 const unsigned int channel = 0)
4365 exchanger.compress_finish(channel, vec);
4372 typename VectorStruct,
4374 typename VectorizedArrayType>
4377 std::vector<VectorStruct> &vec,
4378 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4380 unsigned int component_index = 0;
4381 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4383 compress_finish(vec[comp], exchanger, component_index);
4384 component_index += n_components(vec[comp]);
4392 typename VectorStruct,
4394 typename VectorizedArrayType>
4397 std::vector<VectorStruct *> &vec,
4398 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4400 unsigned int component_index = 0;
4401 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4403 compress_finish(*vec[comp], exchanger, component_index);
4404 component_index += n_components(*vec[comp]);
4420 typename VectorStruct,
4422 typename VectorizedArrayType,
4423 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4427 const VectorStruct &vec,
4428 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4431 if (exchanger.ghosts_were_set ==
true)
4434 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4435 reset_ghost_values(vec.block(i), exchanger);
4442 typename VectorStruct,
4444 typename VectorizedArrayType,
4445 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4449 const VectorStruct &vec,
4450 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4452 exchanger.reset_ghost_values(vec);
4459 typename VectorStruct,
4461 typename VectorizedArrayType>
4464 const std::vector<VectorStruct> &vec,
4465 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4468 if (exchanger.ghosts_were_set ==
true)
4471 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4472 reset_ghost_values(vec[comp], exchanger);
4479 typename VectorStruct,
4481 typename VectorizedArrayType>
4484 const std::vector<VectorStruct *> &vec,
4485 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4488 if (exchanger.ghosts_were_set ==
true)
4491 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4492 reset_ghost_values(*vec[comp], exchanger);
4503 typename VectorStruct,
4505 typename VectorizedArrayType,
4506 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4510 const unsigned int range_index,
4512 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4514 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4515 exchanger.zero_vector_region(range_index, vec.block(i));
4522 typename VectorStruct,
4524 typename VectorizedArrayType,
4525 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4529 const unsigned int range_index,
4531 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4533 exchanger.zero_vector_region(range_index, vec);
4540 typename VectorStruct,
4542 typename VectorizedArrayType>
4545 const unsigned int range_index,
4546 std::vector<VectorStruct> &vec,
4547 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4549 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4550 zero_vector_region(range_index, vec[comp], exchanger);
4557 typename VectorStruct,
4559 typename VectorizedArrayType>
4562 const unsigned int range_index,
4563 std::vector<VectorStruct *> &vec,
4564 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4566 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4567 zero_vector_region(range_index, *vec[comp], exchanger);
4575 template <
typename VectorStruct1,
typename VectorStruct2>
4577 apply_operation_to_constrained_dofs(
const std::vector<unsigned int> &,
4578 const VectorStruct1 &,
4582 template <
typename Number>
4584 apply_operation_to_constrained_dofs(
4585 const std::vector<unsigned int> &constrained_dofs,
4589 for (
const unsigned int i : constrained_dofs)
4590 dst.local_element(i) = src.local_element(i);
4594 namespace MatrixFreeFunctions
4598 template <
typename,
typename,
typename,
typename,
bool>
4599 struct InterfaceSelector
4603 template <
typename MF,
4607 struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4609 using function_type = void (Container::*)(
4613 const
std::pair<unsigned
int, unsigned
int> &) const;
4617 template <
typename MF,
4621 struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4623 using function_type =
4624 void (Container::*)(const MF &,
4627 const
std::pair<unsigned
int, unsigned
int> &);
4635 template <
typename MF,
4640 class MFWorker :
public MFWorkerInterface
4644 using function_type =
typename MatrixFreeFunctions::
4645 InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4649 MFWorker(
const MF &matrix_free,
4650 const InVector &src,
4652 const bool zero_dst_vector_setting,
4653 const Container &container,
4654 function_type cell_function,
4655 function_type face_function,
4656 function_type boundary_function,
4657 const typename MF::DataAccessOnFaces src_vector_face_access =
4658 MF::DataAccessOnFaces::none,
4659 const typename MF::DataAccessOnFaces dst_vector_face_access =
4660 MF::DataAccessOnFaces::none,
4661 const std::function<
void(
const unsigned int,
const unsigned int)>
4662 &operation_before_loop = {},
4663 const std::function<void(
const unsigned int,
const unsigned int)>
4664 &operation_after_loop = {},
4665 const unsigned int dof_handler_index_pre_post = 0)
4666 : matrix_free(matrix_free)
4667 , container(const_cast<Container &>(container))
4668 , cell_function(cell_function)
4669 , face_function(face_function)
4670 , boundary_function(boundary_function)
4673 , src_data_exchanger(matrix_free,
4674 src_vector_face_access,
4676 , dst_data_exchanger(matrix_free,
4677 dst_vector_face_access,
4680 , zero_dst_vector_setting(zero_dst_vector_setting &&
4681 !src_and_dst_are_same)
4682 , operation_before_loop(operation_before_loop)
4683 , operation_after_loop(operation_after_loop)
4684 , dof_handler_index_pre_post(dof_handler_index_pre_post)
4686 Assert(!has_ghost_elements(dst),
4687 ExcMessage(
"The destination vector passed to the matrix-free "
4688 "loop is ghosted. This is not allowed."));
4693 cell(
const std::pair<unsigned int, unsigned int> &cell_range)
override
4695 if (cell_function !=
nullptr && cell_range.second > cell_range.first)
4696 for (
unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4698 const auto cell_subrange =
4699 matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4701 if (cell_subrange.second <= cell_subrange.first)
4705 cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4710 cell(
const unsigned int range_index)
override
4712 process_range(cell_function,
4713 matrix_free.get_task_info().cell_partition_data_hp_ptr,
4714 matrix_free.get_task_info().cell_partition_data_hp,
4719 face(
const unsigned int range_index)
override
4721 process_range(face_function,
4722 matrix_free.get_task_info().face_partition_data_hp_ptr,
4723 matrix_free.get_task_info().face_partition_data_hp,
4728 boundary(
const unsigned int range_index)
override
4730 process_range(boundary_function,
4731 matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4732 matrix_free.get_task_info().boundary_partition_data_hp,
4738 process_range(
const function_type &fu,
4739 const std::vector<unsigned int> &ptr,
4740 const std::vector<unsigned int> &data,
4741 const unsigned int range_index)
4747 for (
unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4750 (container.*fu)(matrix_free,
4753 std::make_pair(data[2 * i], data[2 * i + 1]));
4765 vector_update_ghosts_start()
override
4767 if (!src_and_dst_are_same)
4768 internal::update_ghost_values_start(src, src_data_exchanger);
4773 vector_update_ghosts_finish()
override
4775 if (!src_and_dst_are_same)
4776 internal::update_ghost_values_finish(src, src_data_exchanger);
4781 vector_compress_start()
override
4783 internal::compress_start(dst, dst_data_exchanger);
4788 vector_compress_finish()
override
4790 internal::compress_finish(dst, dst_data_exchanger);
4791 if (!src_and_dst_are_same)
4792 internal::reset_ghost_values(src, src_data_exchanger);
4797 zero_dst_vector_range(
const unsigned int range_index)
override
4799 if (zero_dst_vector_setting)
4800 internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4804 cell_loop_pre_range(
const unsigned int range_index)
override
4806 if (operation_before_loop)
4809 matrix_free.get_dof_info(dof_handler_index_pre_post);
4816 operation_before_loop,
4823 for (
unsigned int id =
4834 cell_loop_post_range(
const unsigned int range_index)
override
4836 if (operation_after_loop)
4840 const std::vector<unsigned int> &partition_row_index =
4841 matrix_free.get_task_info().partition_row_index;
4843 partition_row_index[partition_row_index.size() - 2] - 1)
4844 apply_operation_to_constrained_dofs(
4845 matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
4850 matrix_free.get_dof_info(dof_handler_index_pre_post);
4857 operation_after_loop,
4864 for (
unsigned int id =
4875 const MF &matrix_free;
4876 Container &container;
4877 function_type cell_function;
4878 function_type face_function;
4879 function_type boundary_function;
4881 const InVector &src;
4883 VectorDataExchange<MF::dimension,
4884 typename MF::value_type,
4885 typename MF::vectorized_value_type>
4887 VectorDataExchange<MF::dimension,
4888 typename MF::value_type,
4889 typename MF::vectorized_value_type>
4891 const bool src_and_dst_are_same;
4892 const bool zero_dst_vector_setting;
4893 const std::function<void(
const unsigned int,
const unsigned int)>
4894 operation_before_loop;
4895 const std::function<void(
const unsigned int,
const unsigned int)>
4896 operation_after_loop;
4897 const unsigned int dof_handler_index_pre_post;
4906 template <
class MF,
typename InVector,
typename OutVector>
4907 struct MFClassWrapper
4909 using function_type =
4910 std::function<void(
const MF &,
4913 const std::pair<unsigned int, unsigned int> &)>;
4915 MFClassWrapper(
const function_type cell,
4916 const function_type face,
4917 const function_type boundary)
4920 , boundary(boundary)
4924 cell_integrator(
const MF &mf,
4926 const InVector &src,
4927 const std::pair<unsigned int, unsigned int> &range)
const
4930 cell(mf, dst, src, range);
4934 face_integrator(
const MF &mf,
4936 const InVector &src,
4937 const std::pair<unsigned int, unsigned int> &range)
const
4940 face(mf, dst, src, range);
4944 boundary_integrator(
4947 const InVector &src,
4948 const std::pair<unsigned int, unsigned int> &range)
const
4951 boundary(mf, dst, src, range);
4954 const function_type cell;
4955 const function_type face;
4956 const function_type boundary;
4963template <
int dim,
typename Number,
typename VectorizedArrayType>
4964template <
typename OutVector,
typename InVector>
4970 const std::pair<unsigned int, unsigned int> &)>
4973 const InVector &src,
4974 const bool zero_dst_vector)
const
4977 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4980 Wrapper wrap(cell_operation,
nullptr,
nullptr);
4981 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4991 &Wrapper::cell_integrator,
4992 &Wrapper::face_integrator,
4993 &Wrapper::boundary_integrator);
4995 task_info.loop(worker);
5000template <
int dim,
typename Number,
typename VectorizedArrayType>
5001template <
typename OutVector,
typename InVector>
5007 const std::pair<unsigned int, unsigned int> &)>
5010 const InVector &src,
5011 const std::function<
void(
const unsigned int,
const unsigned int)>
5012 &operation_before_loop,
5013 const std::function<
void(
const unsigned int,
const unsigned int)>
5014 &operation_after_loop,
5015 const unsigned int dof_handler_index_pre_post)
const
5018 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5021 Wrapper wrap(cell_operation,
nullptr,
nullptr);
5022 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5032 &Wrapper::cell_integrator,
5033 &Wrapper::face_integrator,
5034 &Wrapper::boundary_integrator,
5035 DataAccessOnFaces::none,
5036 DataAccessOnFaces::none,
5037 operation_before_loop,
5038 operation_after_loop,
5039 dof_handler_index_pre_post);
5041 task_info.loop(worker);
5046template <
int dim,
typename Number,
typename VectorizedArrayType>
5047template <
typename OutVector,
typename InVector>
5053 const std::pair<unsigned int, unsigned int> &)>
5058 const std::pair<unsigned int, unsigned int> &)>
5063 const std::pair<unsigned int, unsigned int> &)>
5064 &boundary_operation,
5066 const InVector &src,
5067 const bool zero_dst_vector,
5068 const DataAccessOnFaces dst_vector_face_access,
5069 const DataAccessOnFaces src_vector_face_access)
const
5072 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5075 Wrapper wrap(cell_operation, face_operation, boundary_operation);
5076 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5086 &Wrapper::cell_integrator,
5087 &Wrapper::face_integrator,
5088 &Wrapper::boundary_integrator,
5089 src_vector_face_access,
5090 dst_vector_face_access);
5092 task_info.loop(worker);
5097template <
int dim,
typename Number,
typename VectorizedArrayType>
5098template <
typename CLASS,
typename OutVector,
typename InVector>
5101 void (CLASS::*function_pointer)(
5102 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5105 const
std::pair<unsigned
int, unsigned
int> &) const,
5106 const CLASS *owning_class,
5108 const InVector &src,
5109 const
bool zero_dst_vector) const
5111 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5129template <
int dim,
typename Number,
typename VectorizedArrayType>
5130template <
typename CLASS,
typename OutVector,
typename InVector>
5133 void (CLASS::*function_pointer)(
5134 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5137 const
std::pair<unsigned
int, unsigned
int> &) const,
5138 const CLASS *owning_class,
5140 const InVector &src,
5142 &operation_before_loop,
5144 &operation_after_loop,
5145 const unsigned
int dof_handler_index_pre_post) const
5147 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5162 operation_before_loop,
5163 operation_after_loop,
5164 dof_handler_index_pre_post);
5170template <
int dim,
typename Number,
typename VectorizedArrayType>
5171template <
typename CLASS,
typename OutVector,
typename InVector>
5174 void (CLASS::*cell_operation)(
5175 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5178 const
std::pair<unsigned
int, unsigned
int> &) const,
5179 void (CLASS::*face_operation)(
5180 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5183 const
std::pair<unsigned
int, unsigned
int> &) const,
5184 void (CLASS::*boundary_operation)(
5185 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5188 const
std::pair<unsigned
int, unsigned
int> &) const,
5189 const CLASS *owning_class,
5191 const InVector &src,
5192 const
bool zero_dst_vector,
5193 const DataAccessOnFaces dst_vector_face_access,
5194 const DataAccessOnFaces src_vector_face_access) const
5196 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5209 src_vector_face_access,
5210 dst_vector_face_access);
5216template <
int dim,
typename Number,
typename VectorizedArrayType>
5217template <
typename CLASS,
typename OutVector,
typename InVector>
5220 void (CLASS::*function_pointer)(
5221 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5224 const
std::pair<unsigned
int, unsigned
int> &),
5225 CLASS *owning_class,
5227 const InVector &src,
5228 const
bool zero_dst_vector) const
5230 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5248template <
int dim,
typename Number,
typename VectorizedArrayType>
5249template <
typename CLASS,
typename OutVector,
typename InVector>
5252 void (CLASS::*function_pointer)(
5253 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5256 const
std::pair<unsigned
int, unsigned
int> &),
5257 CLASS *owning_class,
5259 const InVector &src,
5261 &operation_before_loop,
5263 &operation_after_loop,
5264 const unsigned
int dof_handler_index_pre_post) const
5266 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5281 operation_before_loop,
5282 operation_after_loop,
5283 dof_handler_index_pre_post);
5289template <
int dim,
typename Number,
typename VectorizedArrayType>
5290template <
typename CLASS,
typename OutVector,
typename InVector>
5293 void (CLASS::*cell_operation)(
5294 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5297 const
std::pair<unsigned
int, unsigned
int> &),
5298 void (CLASS::*face_operation)(
5299 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5302 const
std::pair<unsigned
int, unsigned
int> &),
5303 void (CLASS::*boundary_operation)(
5304 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5307 const
std::pair<unsigned
int, unsigned
int> &),
5308 CLASS *owning_class,
5310 const InVector &src,
5311 const
bool zero_dst_vector,
5312 const DataAccessOnFaces dst_vector_face_access,
5313 const DataAccessOnFaces src_vector_face_access) const
5315 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5328 src_vector_face_access,
5329 dst_vector_face_access);
5335template <
int dim,
typename Number,
typename VectorizedArrayType>
5336template <
typename OutVector,
typename InVector>
5342 const std::pair<unsigned int, unsigned int> &)>
5347 const std::pair<unsigned int, unsigned int> &)>
5352 const std::pair<unsigned int, unsigned int> &)>
5353 &boundary_operation,
5355 const InVector &src,
5356 const std::function<
void(
const unsigned int,
const unsigned int)>
5357 &operation_before_loop,
5358 const std::function<
void(
const unsigned int,
const unsigned int)>
5359 &operation_after_loop,
5360 const unsigned int dof_handler_index_pre_post,
5361 const DataAccessOnFaces dst_vector_face_access,
5362 const DataAccessOnFaces src_vector_face_access)
const
5365 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5368 Wrapper wrap(cell_operation, face_operation, boundary_operation);
5369 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5379 &Wrapper::cell_integrator,
5380 &Wrapper::face_integrator,
5381 &Wrapper::boundary_integrator,
5382 src_vector_face_access,
5383 dst_vector_face_access,
5384 operation_before_loop,
5385 operation_after_loop,
5386 dof_handler_index_pre_post);
5388 task_info.loop(worker);
5393template <
int dim,
typename Number,
typename VectorizedArrayType>
5394template <
typename CLASS,
typename OutVector,
typename InVector>
5397 void (CLASS::*cell_operation)(const
MatrixFree &,
5400 const
std::pair<unsigned
int, unsigned
int> &)
5402 void (CLASS::*face_operation)(const
MatrixFree &,
5405 const
std::pair<unsigned
int, unsigned
int> &)
5407 void (CLASS::*boundary_operation)(
5411 const
std::pair<unsigned
int, unsigned
int> &) const,
5412 const CLASS *owning_class,
5414 const InVector &src,
5416 &operation_before_loop,
5418 &operation_after_loop,
5419 const unsigned
int dof_handler_index_pre_post,
5420 const DataAccessOnFaces dst_vector_face_access,
5421 const DataAccessOnFaces src_vector_face_access) const
5423 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5436 src_vector_face_access,
5437 dst_vector_face_access,
5438 operation_before_loop,
5439 operation_after_loop,
5440 dof_handler_index_pre_post);
5446template <
int dim,
typename Number,
typename VectorizedArrayType>
5447template <
typename CLASS,
typename OutVector,
typename InVector>
5450 void (CLASS::*cell_operation)(const
MatrixFree &,
5453 const
std::pair<unsigned
int, unsigned
int> &),
5454 void (CLASS::*face_operation)(const
MatrixFree &,
5457 const
std::pair<unsigned
int, unsigned
int> &),
5458 void (CLASS::*boundary_operation)(
5462 const
std::pair<unsigned
int, unsigned
int> &),
5463 const CLASS *owning_class,
5465 const InVector &src,
5467 &operation_before_loop,
5469 &operation_after_loop,
5470 const unsigned
int dof_handler_index_pre_post,
5471 const DataAccessOnFaces dst_vector_face_access,
5472 const DataAccessOnFaces src_vector_face_access) const
5474 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5487 src_vector_face_access,
5488 dst_vector_face_access,
5489 operation_before_loop,
5490 operation_after_loop,
5491 dof_handler_index_pre_post);
5497template <
int dim,
typename Number,
typename VectorizedArrayType>
5498template <
typename CLASS,
typename OutVector,
typename InVector>
5501 void (CLASS::*function_pointer)(
5502 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5505 const
std::pair<unsigned
int, unsigned
int> &) const,
5506 const CLASS *owning_class,
5508 const InVector &src,
5509 const
bool zero_dst_vector,
5510 const DataAccessOnFaces src_vector_face_access) const
5512 auto src_vector_face_access_temp = src_vector_face_access;
5518 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5531 src_vector_face_access_temp,
5538template <
int dim,
typename Number,
typename VectorizedArrayType>
5539template <
typename CLASS,
typename OutVector,
typename InVector>
5542 void (CLASS::*function_pointer)(
5543 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5546 const
std::pair<unsigned
int, unsigned
int> &),
5547 CLASS *owning_class,
5549 const InVector &src,
5550 const
bool zero_dst_vector,
5551 const DataAccessOnFaces src_vector_face_access) const
5553 auto src_vector_face_access_temp = src_vector_face_access;
5559 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5572 src_vector_face_access_temp,
5579template <
int dim,
typename Number,
typename VectorizedArrayType>
5580template <
typename OutVector,
typename InVector>
5586 const std::pair<unsigned int, unsigned int> &)>
5589 const InVector &src,
5590 const bool zero_dst_vector,
5591 const DataAccessOnFaces src_vector_face_access)
const
5593 auto src_vector_face_access_temp = src_vector_face_access;
5594 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5595 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5596 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5597 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5600 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5603 Wrapper wrap(cell_operation,
nullptr,
nullptr);
5605 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5615 &Wrapper::cell_integrator,
5616 &Wrapper::face_integrator,
5617 &Wrapper::boundary_integrator,
5618 src_vector_face_access_temp,
5619 DataAccessOnFaces::none);
5620 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 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 > &)> &face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_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
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
void initialize_dof_handlers(const std::vector< const DoFHandler< dim, dim > * > &dof_handlers, const AdditionalData &additional_data)
types::boundary_id get_boundary_id(const unsigned int face_batch_index) const
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > > shape_info
void print(std::ostream &out) const
void update_mapping(const std::shared_ptr< hp::MappingCollection< dim > > &mapping)
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
unsigned int n_inner_face_batches() const
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
void update_mapping(const Mapping< dim > &mapping)
const Quadrature< dim > & get_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_index) const
unsigned int get_mg_level() const
void 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
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
bool at_irregular_cell(const unsigned int cell_batch_index) const
~MatrixFree() override=default
const AffineConstraints< Number > & get_affine_constraints(const unsigned int dof_handler_index=0) const
void copy_from(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free_base)
void initialize_face_data_vector(AlignedVector< T > &vec) const
unsigned int get_dofs_per_cell(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int n_constraint_pool_entries() const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
const Number * constraint_pool_begin(const unsigned int pool_index) const
AlignedVector< Number > * acquire_scratch_data_non_threadsafe() const
const IndexSet & get_locally_owned_set(const unsigned int dof_handler_index=0) const
unsigned int get_dofs_per_face(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
internal::MatrixFreeFunctions::TaskInfo task_info
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*boundary_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
bool mapping_is_initialized
void loop_cell_centric(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void initialize_dof_vector(LinearAlgebra::distributed::Vector< Number2 > &vec, const unsigned int dof_handler_index=0) const
DoFHandler< dim >::cell_iterator get_cell_iterator(const unsigned int cell_batch_index, const unsigned int lane_index, const unsigned int dof_handler_index=0) const
MatrixFree(const MatrixFree< dim, Number, VectorizedArrayType > &other)
std::pair< unsigned int, unsigned int > create_cell_subrange_hp(const std::pair< unsigned int, unsigned int > &range, const unsigned int fe_degree, const unsigned int dof_handler_index=0) const
std::pair< unsigned int, unsigned int > create_cell_subrange_hp_by_index(const std::pair< unsigned int, unsigned int > &range, const unsigned int fe_index, const unsigned int dof_handler_index=0) const
std::pair< unsigned int, unsigned int > get_face_range_category(const std::pair< unsigned int, unsigned int > face_batch_range) const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
std::vector< std::pair< unsigned int, unsigned int > > cell_level_index
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner(const unsigned int dof_handler_index=0) const
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
std::vector< Number > constraint_pool_data
VectorizedArrayType vectorized_value_type
unsigned int n_cell_batches() const
const IndexSet & get_ghost_set(const unsigned int dof_handler_index=0) const
std::pair< int, int > get_cell_level_and_index(const unsigned int cell_batch_index, const unsigned int lane_index) const
bool indices_initialized() const
const Quadrature< dim - 1 > & get_face_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int get_cell_category(const unsigned int cell_batch_index) const
Threads::ThreadLocalStorage< std::list< std::pair< bool, AlignedVector< VectorizedArrayType > > > > scratch_pad
unsigned int get_matrix_free_cell_index(const typename Triangulation< dim >::cell_iterator &cell) const
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
std::list< std::pair< bool, AlignedVector< Number > > > scratch_pad_non_threadsafe
const Number * constraint_pool_end(const unsigned int pool_index) const
void loop_cell_centric(const std::function< void(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
std::vector< SmartPointer< const DoFHandler< dim > > > dof_handlers
void initialize_indices(const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const AdditionalData &additional_data)
std::vector< SmartPointer< const AffineConstraints< Number > > > affine_constraints
unsigned int n_components() const
unsigned int n_physical_cells() const
void initialize_dof_vector(VectorType &vec, const unsigned int dof_handler_index=0) const
unsigned int get_n_q_points(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
void reinit(const MappingType &mapping, const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
std::vector< unsigned int > mf_cell_indices
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*boundary_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
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 loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*boundary_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 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
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 > &)> &face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_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 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
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*boundary_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
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.
unsigned int minimum_parallel_grain_size
const types::boundary_id invalid_boundary_id
static const unsigned int invalid_unsigned_int
void apply_to_subranges(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, const Function &f, const unsigned int grainsize)
const InputIterator OutputIterator const Function & function
::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
AdditionalData & operator=(const AdditionalData &other)
bool cell_vectorization_categories_strict
UpdateFlags mapping_update_flags_boundary_faces
UpdateFlags mapping_update_flags
UpdateFlags mapping_update_flags_faces_by_cells
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