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>
1093 const std::function<
1097 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1098 const std::function<
void(
1102 const std::pair<unsigned int, unsigned int> &)> &inner_face_operation,
1103 const std::function<
void(
1107 const std::pair<unsigned int, unsigned int> &)> &boundary_face_operation,
1109 const InVector &src,
1110 const bool zero_dst_vector =
false,
1204 template <
typename CLASS,
typename OutVector,
typename InVector>
1206 loop(
void (CLASS::*cell_operation)(
1210 const
std::pair<unsigned
int, unsigned
int> &) const,
1211 void (CLASS::*inner_face_operation)(
1215 const
std::pair<unsigned
int, unsigned
int> &) const,
1216 void (CLASS::*boundary_face_operation)(
1220 const
std::pair<unsigned
int, unsigned
int> &) const,
1221 const CLASS *owning_class,
1223 const InVector &src,
1224 const
bool zero_dst_vector = false,
1233 template <
typename CLASS,
typename OutVector,
typename InVector>
1235 loop(
void (CLASS::*cell_operation)(
1239 const
std::pair<unsigned
int, unsigned
int> &),
1240 void (CLASS::*inner_face_operation)(
1244 const
std::pair<unsigned
int, unsigned
int> &),
1245 void (CLASS::*boundary_face_operation)(
1249 const
std::pair<unsigned
int, unsigned
int> &),
1250 CLASS *owning_class,
1252 const InVector &src,
1253 const
bool zero_dst_vector = false,
1380 template <
typename CLASS,
typename OutVector,
typename InVector>
1382 loop(
void (CLASS::*cell_operation)(
1386 const
std::pair<unsigned
int, unsigned
int> &) const,
1387 void (CLASS::*inner_face_operation)(
1391 const
std::pair<unsigned
int, unsigned
int> &) const,
1392 void (CLASS::*boundary_face_operation)(
1396 const
std::pair<unsigned
int, unsigned
int> &) const,
1397 const CLASS *owning_class,
1399 const InVector &src,
1400 const
std::function<void(const unsigned
int, const unsigned
int)>
1401 &operation_before_loop,
1402 const
std::function<void(const unsigned
int, const unsigned
int)>
1403 &operation_after_loop,
1404 const unsigned
int dof_handler_index_pre_post = 0,
1413 template <
typename CLASS,
typename OutVector,
typename InVector>
1415 loop(
void (CLASS::*cell_operation)(
1419 const
std::pair<unsigned
int, unsigned
int> &),
1420 void (CLASS::*inner_face_operation)(
1424 const
std::pair<unsigned
int, unsigned
int> &),
1425 void (CLASS::*boundary_face_operation)(
1429 const
std::pair<unsigned
int, unsigned
int> &),
1430 const CLASS *owning_class,
1432 const InVector &src,
1433 const
std::function<void(const unsigned
int, const unsigned
int)>
1434 &operation_before_loop,
1435 const
std::function<void(const unsigned
int, const unsigned
int)>
1436 &operation_after_loop,
1437 const unsigned
int dof_handler_index_pre_post = 0,
1448 template <
typename OutVector,
typename InVector>
1451 const std::function<
1455 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1456 const std::function<
void(
1460 const std::pair<unsigned int, unsigned int> &)> &inner_face_operation,
1461 const std::function<
void(
1465 const std::pair<unsigned int, unsigned int> &)> &boundary_face_operation,
1467 const InVector &src,
1468 const std::function<
void(
const unsigned int,
const unsigned int)>
1469 &operation_before_loop,
1470 const std::function<
void(
const unsigned int,
const unsigned int)>
1471 &operation_after_loop,
1472 const unsigned int dof_handler_index_pre_post = 0,
1542 template <
typename CLASS,
typename OutVector,
typename InVector>
1548 const
std::pair<unsigned
int, unsigned
int> &) const,
1549 const CLASS *owning_class,
1551 const InVector &src,
1552 const
bool zero_dst_vector = false,
1559 template <
typename CLASS,
typename OutVector,
typename InVector>
1565 const
std::pair<unsigned
int, unsigned
int> &),
1566 CLASS *owning_class,
1568 const InVector &src,
1569 const
bool zero_dst_vector = false,
1576 template <
typename OutVector,
typename InVector>
1582 const std::pair<unsigned int, unsigned int> &)>
1585 const InVector &src,
1586 const bool zero_dst_vector =
false,
1597 std::pair<unsigned int, unsigned int>
1599 const unsigned int fe_degree,
1600 const unsigned int dof_handler_index = 0)
const;
1608 std::pair<unsigned int, unsigned int>
1610 const std::pair<unsigned int, unsigned int> &range,
1611 const unsigned int fe_index,
1612 const unsigned int dof_handler_index = 0)
const;
1625 const std::pair<unsigned int, unsigned int> range)
const;
1632 const bool is_interior_face =
true)
const;
1645 template <
typename T>
1654 template <
typename T>
1671 template <
typename VectorType>
1674 const unsigned int dof_handler_index = 0)
const;
1693 template <
typename Number2>
1696 const unsigned int dof_handler_index = 0)
const;
1708 const std::shared_ptr<const Utilities::MPI::Partitioner> &
1732 const std::vector<unsigned int> &
1747 const unsigned int dof_handler_index = 0);
1758 template <
int spacedim>
1857 const unsigned int face_number)
const;
1889 const unsigned int lane_index,
1890 const unsigned int dof_handler_index = 0)
const;
1899 const unsigned int lane_index)
const;
1923 std::pair<typename DoFHandler<dim>::cell_iterator,
unsigned int>
1925 const unsigned int lane_index,
1926 const bool interior =
true,
1927 const unsigned int fe_component = 0)
const;
1973 const unsigned int hp_active_fe_index = 0)
const;
1980 const unsigned int hp_active_fe_index = 0)
const;
1988 const unsigned int hp_active_fe_index = 0)
const;
1996 const unsigned int hp_active_fe_index = 0)
const;
2003 const unsigned int hp_active_fe_index = 0)
const;
2010 const unsigned int hp_active_fe_index = 0)
const;
2026 const std::pair<unsigned int, unsigned int> cell_batch_range)
const;
2032 std::pair<unsigned int, unsigned int>
2034 const std::pair<unsigned int, unsigned int> face_batch_range)
const;
2054 std::pair<unsigned int, unsigned int>
2088 template <
typename StreamType>
2116 const internal::MatrixFreeFunctions::
2117 MappingInfo<dim, Number, VectorizedArrayType> &
2152 const unsigned int quad_index = 0,
2153 const unsigned int fe_base_element = 0,
2154 const unsigned int hp_active_fe_index = 0,
2155 const unsigned int hp_active_quad_index = 0)
const;
2161 VectorizedArrayType::size()> &
2221 template <
typename number2,
int q_dim>
2227 const std::vector<IndexSet> &locally_owned_set,
2237 template <
typename number2>
2241 const std::vector<IndexSet> &locally_owned_set,
2269 std::vector<internal::MatrixFreeFunctions::DoFInfo>
dof_info;
2352 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2359 mutable std::list<std::pair<bool, AlignedVector<Number>>>
2376template <
int dim,
typename Number,
typename VectorizedArrayType>
2377template <
typename T>
2382 vec.
resize(this->n_cell_batches() + this->n_ghost_cell_batches());
2387template <
int dim,
typename Number,
typename VectorizedArrayType>
2388template <
typename T>
2393 vec.
resize(this->n_inner_face_batches() + this->n_boundary_face_batches() +
2394 this->n_ghost_inner_face_batches());
2399template <
int dim,
typename Number,
typename VectorizedArrayType>
2400template <
typename VectorType>
2404 const unsigned int comp)
const
2407 "This function is not supported for block vectors.");
2409 Assert(task_info.n_procs == 1,
2410 ExcMessage(
"This function can only be used in serial."));
2413 vec.reinit(dof_info[comp].vector_partitioner->size());
2418template <
int dim,
typename Number,
typename VectorizedArrayType>
2419template <
typename Number2>
2423 const unsigned int comp)
const
2426 vec.
reinit(dof_info[comp].vector_partitioner, task_info.communicator_sm);
2431template <
int dim,
typename Number,
typename VectorizedArrayType>
2432inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2434 const unsigned int comp)
const
2437 return dof_info[comp].vector_partitioner;
2442template <
int dim,
typename Number,
typename VectorizedArrayType>
2443inline const std::vector<unsigned int> &
2445 const unsigned int comp)
const
2448 return dof_info[comp].constrained_dofs;
2453template <
int dim,
typename Number,
typename VectorizedArrayType>
2458 return dof_handlers.size();
2463template <
int dim,
typename Number,
typename VectorizedArrayType>
2466 const unsigned int dof_no)
const
2470 return dof_handlers[dof_no]->get_fe().n_base_elements();
2475template <
int dim,
typename Number,
typename VectorizedArrayType>
2484template <
int dim,
typename Number,
typename VectorizedArrayType>
2488 return task_info.n_active_cells;
2493template <
int dim,
typename Number,
typename VectorizedArrayType>
2497 return *(task_info.cell_partition_data.end() - 2);
2502template <
int dim,
typename Number,
typename VectorizedArrayType>
2506 return *(task_info.cell_partition_data.end() - 1) -
2507 *(task_info.cell_partition_data.end() - 2);
2512template <
int dim,
typename Number,
typename VectorizedArrayType>
2516 if (task_info.face_partition_data.empty())
2518 return task_info.face_partition_data.back();
2523template <
int dim,
typename Number,
typename VectorizedArrayType>
2527 if (task_info.face_partition_data.empty())
2529 return task_info.boundary_partition_data.back() -
2530 task_info.face_partition_data.back();
2535template <
int dim,
typename Number,
typename VectorizedArrayType>
2539 if (task_info.face_partition_data.empty())
2541 return face_info.faces.size() - task_info.boundary_partition_data.back();
2546template <
int dim,
typename Number,
typename VectorizedArrayType>
2549 const unsigned int face_batch_index)
const
2551 Assert(face_batch_index >= task_info.boundary_partition_data[0] &&
2552 face_batch_index < task_info.boundary_partition_data.back(),
2554 task_info.boundary_partition_data[0],
2555 task_info.boundary_partition_data.back()));
2561template <
int dim,
typename Number,
typename VectorizedArrayType>
2564 const unsigned int cell_batch_index,
2565 const unsigned int face_number)
const
2569 Assert(face_info.cell_and_face_boundary_id.size(0) >= n_cell_batches(),
2573 for (
unsigned int v = 0;
2574 v < n_active_entries_per_cell_batch(cell_batch_index);
2577 face_info.cell_and_face_boundary_id(cell_batch_index, face_number, v);
2583template <
int dim,
typename Number,
typename VectorizedArrayType>
2584inline const internal::MatrixFreeFunctions::
2585 MappingInfo<dim, Number, VectorizedArrayType> &
2588 return mapping_info;
2593template <
int dim,
typename Number,
typename VectorizedArrayType>
2596 const unsigned int dof_index)
const
2599 return dof_info[dof_index];
2604template <
int dim,
typename Number,
typename VectorizedArrayType>
2608 return constraint_pool_row_index.size() - 1;
2613template <
int dim,
typename Number,
typename VectorizedArrayType>
2614inline const Number *
2616 const unsigned int row)
const
2619 return constraint_pool_data.empty() ?
2621 constraint_pool_data.data() + constraint_pool_row_index[row];
2626template <
int dim,
typename Number,
typename VectorizedArrayType>
2627inline const Number *
2629 const unsigned int row)
const
2632 return constraint_pool_data.empty() ?
2634 constraint_pool_data.data() + constraint_pool_row_index[row + 1];
2639template <
int dim,
typename Number,
typename VectorizedArrayType>
2640inline std::pair<unsigned int, unsigned int>
2642 const std::pair<unsigned int, unsigned int> &range,
2643 const unsigned int degree,
2644 const unsigned int dof_handler_component)
const
2646 if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2649 dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2651 dof_info[dof_handler_component].fe_index_conversion[0].size(), 1);
2652 if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2655 return {range.second, range.second};
2659 dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2660 if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2661 return {range.second, range.second};
2663 return create_cell_subrange_hp_by_index(range,
2665 dof_handler_component);
2670template <
int dim,
typename Number,
typename VectorizedArrayType>
2673 const unsigned int cell_batch_index)
const
2676 return VectorizedArrayType::size() > 1 &&
2677 cell_level_index[(cell_batch_index + 1) * VectorizedArrayType::size() -
2678 1] == cell_level_index[(cell_batch_index + 1) *
2679 VectorizedArrayType::size() -
2685template <
int dim,
typename Number,
typename VectorizedArrayType>
2689 return shape_info.size(2);
2693template <
int dim,
typename Number,
typename VectorizedArrayType>
2696 const std::pair<unsigned int, unsigned int> range)
const
2698 const auto &fe_indices = dof_info[0].cell_active_fe_index;
2700 if (fe_indices.empty() ==
true ||
2701 dof_handlers[0]->get_fe_collection().size() == 1)
2704 const auto index = fe_indices[range.first];
2706 for (
unsigned int i = range.first; i < range.second; ++i)
2714template <
int dim,
typename Number,
typename VectorizedArrayType>
2717 const std::pair<unsigned int, unsigned int> range,
2718 const bool is_interior_face)
const
2720 const auto &fe_indices = dof_info[0].cell_active_fe_index;
2722 if (fe_indices.empty() ==
true)
2725 if (is_interior_face)
2727 const unsigned int index =
2728 fe_indices[face_info.faces[range.first].cells_interior[0] /
2729 VectorizedArrayType::size()];
2731 for (
unsigned int i = range.first; i < range.second; ++i)
2733 fe_indices[face_info.faces[i].cells_interior[0] /
2734 VectorizedArrayType::size()]);
2740 const unsigned int index =
2741 fe_indices[face_info.faces[range.first].cells_exterior[0] /
2742 VectorizedArrayType::size()];
2744 for (
unsigned int i = range.first; i < range.second; ++i)
2746 fe_indices[face_info.faces[i].cells_exterior[0] /
2747 VectorizedArrayType::size()]);
2755template <
int dim,
typename Number,
typename VectorizedArrayType>
2758 const unsigned int cell_batch_index)
const
2762 const std::vector<unsigned char> &n_lanes_filled =
2763 dof_info[0].n_vectorization_lanes_filled
2767 return n_lanes_filled[cell_batch_index];
2772template <
int dim,
typename Number,
typename VectorizedArrayType>
2775 const unsigned int face_batch_index)
const
2779 const std::vector<unsigned char> &n_lanes_filled =
2780 dof_info[0].n_vectorization_lanes_filled
2783 return n_lanes_filled[face_batch_index];
2788template <
int dim,
typename Number,
typename VectorizedArrayType>
2791 const unsigned int dof_handler_index,
2792 const unsigned int active_fe_index)
const
2794 return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2799template <
int dim,
typename Number,
typename VectorizedArrayType>
2802 const unsigned int quad_index,
2803 const unsigned int active_fe_index)
const
2806 return mapping_info.cell_data[quad_index]
2807 .descriptor[active_fe_index]
2813template <
int dim,
typename Number,
typename VectorizedArrayType>
2816 const unsigned int dof_handler_index,
2817 const unsigned int active_fe_index)
const
2819 return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2824template <
int dim,
typename Number,
typename VectorizedArrayType>
2827 const unsigned int quad_index,
2828 const unsigned int active_fe_index)
const
2831 return mapping_info.face_data[quad_index]
2832 .descriptor[active_fe_index]
2838template <
int dim,
typename Number,
typename VectorizedArrayType>
2841 const unsigned int dof_handler_index)
const
2843 return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2848template <
int dim,
typename Number,
typename VectorizedArrayType>
2851 const unsigned int dof_handler_index)
const
2853 return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2858template <
int dim,
typename Number,
typename VectorizedArrayType>
2861 const unsigned int dof_handler_index,
2862 const unsigned int index_quad,
2863 const unsigned int index_fe,
2864 const unsigned int active_fe_index,
2865 const unsigned int active_quad_index)
const
2868 const unsigned int ind =
2869 dof_info[dof_handler_index].global_base_element_offset + index_fe;
2874 return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2879template <
int dim,
typename Number,
typename VectorizedArrayType>
2881 VectorizedArrayType::size()> &
2883 const unsigned int face_batch_index)
const
2886 return face_info.faces[face_batch_index];
2891template <
int dim,
typename Number,
typename VectorizedArrayType>
2896 return face_info.cell_and_face_to_plain_faces;
2901template <
int dim,
typename Number,
typename VectorizedArrayType>
2904 const unsigned int quad_index,
2905 const unsigned int active_fe_index)
const
2908 return mapping_info.cell_data[quad_index]
2909 .descriptor[active_fe_index]
2915template <
int dim,
typename Number,
typename VectorizedArrayType>
2918 const unsigned int quad_index,
2919 const unsigned int active_fe_index)
const
2922 return mapping_info.face_data[quad_index]
2923 .descriptor[active_fe_index]
2929template <
int dim,
typename Number,
typename VectorizedArrayType>
2932 const std::pair<unsigned int, unsigned int> range)
const
2934 auto result = get_cell_category(range.first);
2936 for (
unsigned int i = range.first; i < range.second; ++i)
2937 result =
std::max(result, get_cell_category(i));
2944template <
int dim,
typename Number,
typename VectorizedArrayType>
2945inline std::pair<unsigned int, unsigned int>
2947 const std::pair<unsigned int, unsigned int> range)
const
2949 auto result = get_face_category(range.first);
2951 for (
unsigned int i = range.first; i < range.second; ++i)
2953 result.first =
std::max(result.first, get_face_category(i).
first);
2954 result.second =
std::max(result.second, get_face_category(i).
second);
2962template <
int dim,
typename Number,
typename VectorizedArrayType>
2965 const unsigned int cell_batch_index)
const
2968 AssertIndexRange(cell_batch_index, dof_info[0].cell_active_fe_index.size());
2969 if (dof_info[0].cell_active_fe_index.empty())
2972 return dof_info[0].cell_active_fe_index[cell_batch_index];
2977template <
int dim,
typename Number,
typename VectorizedArrayType>
2978inline std::pair<unsigned int, unsigned int>
2980 const unsigned int face_batch_index)
const
2983 if (dof_info[0].cell_active_fe_index.empty())
2984 return std::make_pair(0U, 0U);
2986 std::pair<unsigned int, unsigned int> result = std::make_pair(0U, 0U);
2987 for (
unsigned int v = 0;
2988 v < VectorizedArrayType::size() &&
2989 face_info.faces[face_batch_index].cells_interior[v] !=
2994 dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
2995 .cells_interior[v] /
2996 VectorizedArrayType::size()]);
2997 if (face_info.faces[face_batch_index].cells_exterior[0] !=
2999 for (
unsigned int v = 0;
3000 v < VectorizedArrayType::size() &&
3001 face_info.faces[face_batch_index].cells_exterior[v] !=
3006 dof_info[0].cell_active_fe_index[face_info.faces[face_batch_index]
3007 .cells_exterior[v] /
3008 VectorizedArrayType::size()]);
3016template <
int dim,
typename Number,
typename VectorizedArrayType>
3020 return indices_are_initialized;
3025template <
int dim,
typename Number,
typename VectorizedArrayType>
3029 return mapping_is_initialized;
3033template <
int dim,
typename Number,
typename VectorizedArrayType>
3042template <
int dim,
typename Number,
typename VectorizedArrayType>
3047 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3048 list_type &data = scratch_pad.get();
3049 for (
typename list_type::iterator it = data.begin();
it != data.end(); ++
it)
3050 if (
it->first ==
false)
3056 return &data.front().second;
3061template <
int dim,
typename Number,
typename VectorizedArrayType>
3067 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3068 list_type &data = scratch_pad.get();
3069 for (
typename list_type::iterator it = data.begin();
it != data.end(); ++
it)
3070 if (&
it->second == scratch)
3081template <
int dim,
typename Number,
typename VectorizedArrayType>
3087 scratch_pad_non_threadsafe.
begin();
3088 it != scratch_pad_non_threadsafe.end();
3090 if (
it->first ==
false)
3095 scratch_pad_non_threadsafe.push_front(
3097 return &scratch_pad_non_threadsafe.front().second;
3102template <
int dim,
typename Number,
typename VectorizedArrayType>
3109 scratch_pad_non_threadsafe.
begin();
3110 it != scratch_pad_non_threadsafe.end();
3112 if (&
it->second == scratch)
3127 namespace MatrixFreeImplementation
3129 template <
int dim,
int spacedim>
3130 inline std::vector<IndexSet>
3131 extract_locally_owned_index_sets(
3133 const unsigned int level)
3135 std::vector<IndexSet> locally_owned_set;
3136 locally_owned_set.reserve(dofh.size());
3137 for (
unsigned int j = 0; j < dofh.size(); ++j)
3139 locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
3141 locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(
level));
3142 return locally_owned_set;
3149template <
int dim,
typename Number,
typename VectorizedArrayType>
3150template <
typename QuadratureType,
typename number2,
typename MappingType>
3153 const MappingType &mapping,
3156 const QuadratureType &quad,
3160 std::vector<const DoFHandler<dim, dim> *> dof_handlers;
3161 std::vector<const AffineConstraints<number2> *> constraints;
3163 dof_handlers.push_back(&dof_handler);
3164 constraints.push_back(&constraints_in);
3166 std::vector<IndexSet> locally_owned_sets =
3167 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3168 dof_handlers, additional_data.
mg_level);
3170 std::vector<hp::QCollection<dim>> quad_hp;
3171 quad_hp.emplace_back(quad);
3183template <
int dim,
typename Number,
typename VectorizedArrayType>
3184template <
typename QuadratureType,
typename number2,
typename MappingType>
3187 const MappingType &mapping,
3190 const QuadratureType &quad,
3194 std::vector<IndexSet> locally_owned_set =
3195 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3196 dof_handler, additional_data.
mg_level);
3197 std::vector<hp::QCollection<dim>> quad_hp;
3198 quad_hp.emplace_back(quad);
3210template <
int dim,
typename Number,
typename VectorizedArrayType>
3211template <
typename QuadratureType,
typename number2,
typename MappingType>
3214 const MappingType &mapping,
3217 const std::vector<QuadratureType> &quad,
3221 std::vector<IndexSet> locally_owned_set =
3222 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3223 dof_handler, additional_data.
mg_level);
3224 std::vector<hp::QCollection<dim>> quad_hp;
3225 for (
unsigned int q = 0; q < quad.size(); ++q)
3226 quad_hp.emplace_back(quad[q]);
3253 template <
int dim,
typename Number,
typename VectorizedArrayType>
3254 struct VectorDataExchange
3261 static constexpr unsigned int channel_shift = 20;
3270 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3271 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3272 DataAccessOnFaces vector_face_access,
3273 const unsigned
int n_components)
3274 : matrix_free(matrix_free)
3275 , vector_face_access(
3276 matrix_free.get_task_info().face_partition_data.empty() ?
3277 ::
MatrixFree<dim, Number, VectorizedArrayType>::
3278 DataAccessOnFaces::unspecified :
3280 , ghosts_were_set(false)
3281# ifdef DEAL_II_WITH_MPI
3282 , tmp_data(n_components)
3283 , requests(n_components)
3287 if (this->vector_face_access !=
3289 DataAccessOnFaces::unspecified)
3290 for (unsigned
int c = 0; c < matrix_free.n_components(); ++c)
3292 matrix_free.get_dof_info(c).vector_exchanger_face_variants.size(),
3301 ~VectorDataExchange()
3303# ifdef DEAL_II_WITH_MPI
3304 for (
unsigned int i = 0; i < tmp_data.size(); ++i)
3305 if (tmp_data[i] !=
nullptr)
3306 matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3316 template <
typename VectorType>
3318 find_vector_in_mf(
const VectorType &vec,
3319 const bool check_global_compatibility =
true)
const
3322 for (
unsigned int c = 0; c < matrix_free.n_components(); ++c)
3323 if (vec.get_partitioner().get() ==
3324 matrix_free.get_dof_info(c).vector_partitioner.get())
3328 for (
unsigned int c = 0; c < matrix_free.n_components(); ++c)
3329 if (check_global_compatibility ?
3330 vec.get_partitioner()->is_globally_compatible(
3331 *matrix_free.get_dof_info(c).vector_partitioner) :
3332 vec.get_partitioner()->is_compatible(
3333 *matrix_free.get_dof_info(c).vector_partitioner))
3349 get_partitioner(
const unsigned int mf_component)
const
3352 .vector_exchanger_face_variants.size(),
3354 if (vector_face_access ==
3356 DataAccessOnFaces::
none)
3357 return *matrix_free.get_dof_info(mf_component)
3358 .vector_exchanger_face_variants[0];
3359 else if (vector_face_access ==
3361 DataAccessOnFaces::
values)
3362 return *matrix_free.get_dof_info(mf_component)
3363 .vector_exchanger_face_variants[1];
3364 else if (vector_face_access ==
3367 return *matrix_free.get_dof_info(mf_component)
3368 .vector_exchanger_face_variants[2];
3369 else if (vector_face_access ==
3371 DataAccessOnFaces::values_all_faces)
3372 return *matrix_free.get_dof_info(mf_component)
3373 .vector_exchanger_face_variants[3];
3374 else if (vector_face_access ==
3376 DataAccessOnFaces::gradients_all_faces)
3377 return *matrix_free.get_dof_info(mf_component)
3378 .vector_exchanger_face_variants[4];
3380 return *matrix_free.get_dof_info(mf_component).vector_exchanger.get();
3388 template <
typename VectorType,
3389 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3392 update_ghost_values_start(
const unsigned int ,
3393 const VectorType & )
3401 template <
typename VectorType,
3402 std::enable_if_t<!has_update_ghost_values_start<VectorType> &&
3403 !is_not_parallel_vector<VectorType>,
3404 VectorType> * =
nullptr>
3406 update_ghost_values_start(
const unsigned int component_in_block_vector,
3407 const VectorType &vec)
3409 (void)component_in_block_vector;
3410 const bool ghosts_set = vec.has_ghost_elements();
3412 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3413 ghosts_set ==
false,
3417 ghosts_were_set =
true;
3419 vec.update_ghost_values();
3429 template <
typename VectorType,
3430 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3431 !has_exchange_on_subset<VectorType>,
3432 VectorType> * =
nullptr>
3434 update_ghost_values_start(
const unsigned int component_in_block_vector,
3435 const VectorType &vec)
3437 (void)component_in_block_vector;
3438 const bool ghosts_set = vec.has_ghost_elements();
3440 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3441 ghosts_set ==
false,
3445 ghosts_were_set =
true;
3447 vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3458 template <
typename VectorType,
3459 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3460 has_exchange_on_subset<VectorType>,
3461 VectorType> * =
nullptr>
3463 update_ghost_values_start(
const unsigned int component_in_block_vector,
3464 const VectorType &vec)
3466 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3467 "Type mismatch between VectorType and VectorDataExchange");
3468 (void)component_in_block_vector;
3469 const bool ghosts_set = vec.has_ghost_elements();
3471 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3472 ghosts_set ==
false,
3476 ghosts_were_set =
true;
3478 if (vec.size() != 0)
3480# ifdef DEAL_II_WITH_MPI
3481 const unsigned int mf_component = find_vector_in_mf(vec);
3483 const auto &part = get_partitioner(mf_component);
3485 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3486 part.n_import_sm_procs() == 0)
3489 tmp_data[component_in_block_vector] =
3490 matrix_free.acquire_scratch_data_non_threadsafe();
3491 tmp_data[component_in_block_vector]->resize_fast(
3492 part.n_import_indices());
3495 part.export_to_ghosted_array_start(
3496 component_in_block_vector * 2 + channel_shift,
3498 vec.shared_vector_data(),
3500 part.locally_owned_size(),
3501 matrix_free.get_dof_info(mf_component)
3502 .vector_partitioner->n_ghost_indices()),
3504 part.n_import_indices()),
3505 this->requests[component_in_block_vector]);
3516 template <
typename VectorType,
3517 std::enable_if_t<!has_update_ghost_values_start<VectorType>,
3518 VectorType> * =
nullptr>
3520 update_ghost_values_finish(
const unsigned int ,
3521 const VectorType & )
3531 template <
typename VectorType,
3532 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3533 !has_exchange_on_subset<VectorType>,
3534 VectorType> * =
nullptr>
3536 update_ghost_values_finish(
const unsigned int component_in_block_vector,
3537 const VectorType &vec)
3539 (void)component_in_block_vector;
3540 vec.update_ghost_values_finish();
3551 template <
typename VectorType,
3552 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3553 has_exchange_on_subset<VectorType>,
3554 VectorType> * =
nullptr>
3556 update_ghost_values_finish(
const unsigned int component_in_block_vector,
3557 const VectorType &vec)
3559 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3560 "Type mismatch between VectorType and VectorDataExchange");
3561 (void)component_in_block_vector;
3563 if (vec.size() != 0)
3565# ifdef DEAL_II_WITH_MPI
3569 const unsigned int mf_component = find_vector_in_mf(vec);
3571 const auto &part = get_partitioner(mf_component);
3573 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3574 part.n_import_sm_procs() != 0)
3576 part.export_to_ghosted_array_finish(
3578 vec.shared_vector_data(),
3580 part.locally_owned_size(),
3581 matrix_free.get_dof_info(mf_component)
3582 .vector_partitioner->n_ghost_indices()),
3583 this->requests[component_in_block_vector]);
3585 matrix_free.release_scratch_data_non_threadsafe(
3586 tmp_data[component_in_block_vector]);
3587 tmp_data[component_in_block_vector] =
nullptr;
3593 vec.set_ghost_state(
true);
3601 template <
typename VectorType,
3602 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3605 compress_start(
const unsigned int ,
3615 template <
typename VectorType,
3616 std::enable_if_t<!has_compress_start<VectorType> &&
3617 !is_not_parallel_vector<VectorType>,
3618 VectorType> * =
nullptr>
3620 compress_start(
const unsigned int component_in_block_vector,
3623 (void)component_in_block_vector;
3635 template <
typename VectorType,
3636 std::enable_if_t<has_compress_start<VectorType> &&
3637 !has_exchange_on_subset<VectorType>,
3638 VectorType> * =
nullptr>
3640 compress_start(
const unsigned int component_in_block_vector,
3643 (void)component_in_block_vector;
3645 vec.compress_start(component_in_block_vector + channel_shift);
3656 template <
typename VectorType,
3657 std::enable_if_t<has_compress_start<VectorType> &&
3658 has_exchange_on_subset<VectorType>,
3659 VectorType> * =
nullptr>
3661 compress_start(
const unsigned int component_in_block_vector,
3664 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3665 "Type mismatch between VectorType and VectorDataExchange");
3666 (void)component_in_block_vector;
3669 if (vec.size() != 0)
3671# ifdef DEAL_II_WITH_MPI
3672 const unsigned int mf_component = find_vector_in_mf(vec);
3674 const auto &part = get_partitioner(mf_component);
3676 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3677 part.n_import_sm_procs() == 0)
3680 tmp_data[component_in_block_vector] =
3681 matrix_free.acquire_scratch_data_non_threadsafe();
3682 tmp_data[component_in_block_vector]->resize_fast(
3683 part.n_import_indices());
3686 part.import_from_ghosted_array_start(
3688 component_in_block_vector * 2 + channel_shift,
3690 vec.shared_vector_data(),
3692 matrix_free.get_dof_info(mf_component)
3693 .vector_partitioner->n_ghost_indices()),
3695 part.n_import_indices()),
3696 this->requests[component_in_block_vector]);
3708 typename VectorType,
3709 std::enable_if_t<!has_compress_start<VectorType>, VectorType> * =
nullptr>
3711 compress_finish(
const unsigned int ,
3722 template <
typename VectorType,
3723 std::enable_if_t<has_compress_start<VectorType> &&
3724 !has_exchange_on_subset<VectorType>,
3725 VectorType> * =
nullptr>
3727 compress_finish(
const unsigned int component_in_block_vector,
3730 (void)component_in_block_vector;
3742 template <
typename VectorType,
3743 std::enable_if_t<has_compress_start<VectorType> &&
3744 has_exchange_on_subset<VectorType>,
3745 VectorType> * =
nullptr>
3747 compress_finish(
const unsigned int component_in_block_vector,
3750 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3751 "Type mismatch between VectorType and VectorDataExchange");
3752 (void)component_in_block_vector;
3753 if (vec.size() != 0)
3755# ifdef DEAL_II_WITH_MPI
3759 const unsigned int mf_component = find_vector_in_mf(vec);
3761 const auto &part = get_partitioner(mf_component);
3763 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3764 part.n_import_sm_procs() != 0)
3766 part.import_from_ghosted_array_finish(
3769 vec.shared_vector_data(),
3771 matrix_free.get_dof_info(mf_component)
3772 .vector_partitioner->n_ghost_indices()),
3774 tmp_data[component_in_block_vector]->begin(),
3775 part.n_import_indices()),
3776 this->requests[component_in_block_vector]);
3778 matrix_free.release_scratch_data_non_threadsafe(
3779 tmp_data[component_in_block_vector]);
3780 tmp_data[component_in_block_vector] =
nullptr;
3786 MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3798 template <
typename VectorType,
3799 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType>
3802 reset_ghost_values(
const VectorType & )
const
3811 template <
typename VectorType,
3812 std::enable_if_t<!has_exchange_on_subset<VectorType> &&
3813 !is_not_parallel_vector<VectorType>,
3814 VectorType> * =
nullptr>
3816 reset_ghost_values(
const VectorType &vec)
const
3818 if (ghosts_were_set ==
true)
3821 vec.zero_out_ghost_values();
3831 template <
typename VectorType,
3832 std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3835 reset_ghost_values(
const VectorType &vec)
const
3837 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3838 "Type mismatch between VectorType and VectorDataExchange");
3839 if (ghosts_were_set ==
true)
3842 if (vec.size() != 0)
3844# ifdef DEAL_II_WITH_MPI
3847 const unsigned int mf_component = find_vector_in_mf(vec);
3849 const auto &part = get_partitioner(mf_component);
3851 if (part.n_ghost_indices() > 0)
3856 part.locally_owned_size(),
3857 matrix_free.get_dof_info(mf_component)
3858 .vector_partitioner->n_ghost_indices()));
3864 vec.set_ghost_state(
false);
3874 template <
typename VectorType,
3875 std::enable_if_t<has_exchange_on_subset<VectorType>, VectorType>
3878 zero_vector_region(
const unsigned int range_index, VectorType &vec)
const
3880 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3881 "Type mismatch between VectorType and VectorDataExchange");
3886 const unsigned int mf_component = find_vector_in_mf(vec,
false);
3888 matrix_free.get_dof_info(mf_component);
3896 for (
unsigned int id =
3915 template <
typename VectorType,
3916 std::enable_if_t<!has_exchange_on_subset<VectorType>, VectorType>
3918 typename VectorType::value_type * =
nullptr>
3920 zero_vector_region(
const unsigned int range_index, VectorType &vec)
const
3923 vec =
typename VectorType::value_type();
3933 zero_vector_region(
const unsigned int, ...)
const
3937 "which provide VectorType::value_type"));
3942 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
3943 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3944 DataAccessOnFaces vector_face_access;
3945 bool ghosts_were_set;
3946# ifdef DEAL_II_WITH_MPI
3947 std::vector<AlignedVector<Number> *> tmp_data;
3948 std::vector<std::vector<MPI_Request>> requests;
3952 template <
typename VectorStruct>
3954 n_components(
const VectorStruct &vec);
3956 template <
typename VectorStruct>
3958 n_components_block(
const VectorStruct &vec,
const std::bool_constant<true>)
3960 unsigned int components = 0;
3961 for (
unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
3962 components += n_components(vec.block(bl));
3966 template <
typename VectorStruct>
3968 n_components_block(
const VectorStruct &,
const std::bool_constant<false>)
3973 template <
typename VectorStruct>
3975 n_components(
const VectorStruct &vec)
3977 return n_components_block(
3981 template <
typename VectorStruct>
3983 n_components(
const std::vector<VectorStruct> &vec)
3985 unsigned int components = 0;
3986 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
3987 components += n_components_block(
3992 template <
typename VectorStruct>
3994 n_components(
const std::vector<VectorStruct *> &vec)
3996 unsigned int components = 0;
3997 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
3998 components += n_components_block(
4010 template <
typename VectorStruct,
4011 std::enable_if_t<!has_communication_block_size<VectorStruct>,
4012 VectorStruct> * =
nullptr>
4013 constexpr unsigned int
4014 get_communication_block_size(
const VectorStruct &)
4021 template <
typename VectorStruct,
4022 std::enable_if_t<has_communication_block_size<VectorStruct>,
4023 VectorStruct> * =
nullptr>
4024 constexpr unsigned int
4025 get_communication_block_size(
const VectorStruct &)
4027 return VectorStruct::communication_block_size;
4032 template <
typename VectorType,
4033 std::enable_if_t<is_not_parallel_vector<VectorType>, VectorType> * =
4036 has_ghost_elements(
const VectorType &vec)
4044 template <
typename VectorType,
4045 std::enable_if_t<!is_not_parallel_vector<VectorType>, VectorType>
4048 has_ghost_elements(
const VectorType &vec)
4050 return vec.has_ghost_elements();
4065 typename VectorStruct,
4067 typename VectorizedArrayType,
4068 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4071 update_ghost_values_start(
4072 const VectorStruct &vec,
4073 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4074 const unsigned int channel = 0)
4076 if (get_communication_block_size(vec) < vec.n_blocks())
4078 const bool ghosts_set = vec.has_ghost_elements();
4080 Assert(exchanger.matrix_free.get_task_info()
4081 .allow_ghosted_vectors_in_loops ||
4082 ghosts_set ==
false,
4086 exchanger.ghosts_were_set =
true;
4088 vec.update_ghost_values();
4092 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4093 update_ghost_values_start(vec.block(i), exchanger, channel + i);
4101 typename VectorStruct,
4103 typename VectorizedArrayType,
4104 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4107 update_ghost_values_start(
4108 const VectorStruct &vec,
4109 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4110 const unsigned int channel = 0)
4112 exchanger.update_ghost_values_start(channel, vec);
4119 typename VectorStruct,
4121 typename VectorizedArrayType>
4123 update_ghost_values_start(
4124 const std::vector<VectorStruct> &vec,
4125 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4127 unsigned int component_index = 0;
4128 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4130 update_ghost_values_start(vec[comp], exchanger, component_index);
4131 component_index += n_components(vec[comp]);
4139 typename VectorStruct,
4141 typename VectorizedArrayType>
4143 update_ghost_values_start(
4144 const std::vector<VectorStruct *> &vec,
4145 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4147 unsigned int component_index = 0;
4148 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4150 update_ghost_values_start(*vec[comp], exchanger, component_index);
4151 component_index += n_components(*vec[comp]);
4163 typename VectorStruct,
4165 typename VectorizedArrayType,
4166 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4169 update_ghost_values_finish(
4170 const VectorStruct &vec,
4171 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4172 const unsigned int channel = 0)
4174 if (get_communication_block_size(vec) < vec.n_blocks())
4180 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4181 update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4188 typename VectorStruct,
4190 typename VectorizedArrayType,
4191 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4194 update_ghost_values_finish(
4195 const VectorStruct &vec,
4196 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4197 const unsigned int channel = 0)
4199 exchanger.update_ghost_values_finish(channel, vec);
4206 typename VectorStruct,
4208 typename VectorizedArrayType>
4210 update_ghost_values_finish(
4211 const std::vector<VectorStruct> &vec,
4212 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4214 unsigned int component_index = 0;
4215 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4217 update_ghost_values_finish(vec[comp], exchanger, component_index);
4218 component_index += n_components(vec[comp]);
4226 typename VectorStruct,
4228 typename VectorizedArrayType>
4230 update_ghost_values_finish(
4231 const std::vector<VectorStruct *> &vec,
4232 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4234 unsigned int component_index = 0;
4235 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4237 update_ghost_values_finish(*vec[comp], exchanger, component_index);
4238 component_index += n_components(*vec[comp]);
4250 typename VectorStruct,
4252 typename VectorizedArrayType,
4253 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4258 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4259 const unsigned int channel = 0)
4261 if (get_communication_block_size(vec) < vec.n_blocks())
4264 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4265 compress_start(vec.block(i), exchanger, channel + i);
4272 typename VectorStruct,
4274 typename VectorizedArrayType,
4275 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4280 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4281 const unsigned int channel = 0)
4283 exchanger.compress_start(channel, vec);
4290 typename VectorStruct,
4292 typename VectorizedArrayType>
4295 std::vector<VectorStruct> &vec,
4296 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4298 unsigned int component_index = 0;
4299 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4301 compress_start(vec[comp], exchanger, component_index);
4302 component_index += n_components(vec[comp]);
4310 typename VectorStruct,
4312 typename VectorizedArrayType>
4315 std::vector<VectorStruct *> &vec,
4316 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4318 unsigned int component_index = 0;
4319 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4321 compress_start(*vec[comp], exchanger, component_index);
4322 component_index += n_components(*vec[comp]);
4334 typename VectorStruct,
4336 typename VectorizedArrayType,
4337 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4342 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4343 const unsigned int channel = 0)
4345 if (get_communication_block_size(vec) < vec.n_blocks())
4351 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4352 compress_finish(vec.block(i), exchanger, channel + i);
4359 typename VectorStruct,
4361 typename VectorizedArrayType,
4362 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4367 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4368 const unsigned int channel = 0)
4370 exchanger.compress_finish(channel, vec);
4377 typename VectorStruct,
4379 typename VectorizedArrayType>
4382 std::vector<VectorStruct> &vec,
4383 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4385 unsigned int component_index = 0;
4386 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4388 compress_finish(vec[comp], exchanger, component_index);
4389 component_index += n_components(vec[comp]);
4397 typename VectorStruct,
4399 typename VectorizedArrayType>
4402 std::vector<VectorStruct *> &vec,
4403 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4405 unsigned int component_index = 0;
4406 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4408 compress_finish(*vec[comp], exchanger, component_index);
4409 component_index += n_components(*vec[comp]);
4425 typename VectorStruct,
4427 typename VectorizedArrayType,
4428 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4432 const VectorStruct &vec,
4433 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4436 if (exchanger.ghosts_were_set ==
true)
4439 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4440 reset_ghost_values(vec.block(i), exchanger);
4447 typename VectorStruct,
4449 typename VectorizedArrayType,
4450 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4454 const VectorStruct &vec,
4455 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4457 exchanger.reset_ghost_values(vec);
4464 typename VectorStruct,
4466 typename VectorizedArrayType>
4469 const std::vector<VectorStruct> &vec,
4470 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4473 if (exchanger.ghosts_were_set ==
true)
4476 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4477 reset_ghost_values(vec[comp], exchanger);
4484 typename VectorStruct,
4486 typename VectorizedArrayType>
4489 const std::vector<VectorStruct *> &vec,
4490 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4493 if (exchanger.ghosts_were_set ==
true)
4496 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4497 reset_ghost_values(*vec[comp], exchanger);
4508 typename VectorStruct,
4510 typename VectorizedArrayType,
4511 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4515 const unsigned int range_index,
4517 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4519 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4520 exchanger.zero_vector_region(range_index, vec.block(i));
4527 typename VectorStruct,
4529 typename VectorizedArrayType,
4530 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4534 const unsigned int range_index,
4536 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4538 exchanger.zero_vector_region(range_index, vec);
4545 typename VectorStruct,
4547 typename VectorizedArrayType>
4550 const unsigned int range_index,
4551 std::vector<VectorStruct> &vec,
4552 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4554 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4555 zero_vector_region(range_index, vec[comp], exchanger);
4562 typename VectorStruct,
4564 typename VectorizedArrayType>
4567 const unsigned int range_index,
4568 std::vector<VectorStruct *> &vec,
4569 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4571 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4572 zero_vector_region(range_index, *vec[comp], exchanger);
4580 template <
typename VectorStruct1,
typename VectorStruct2>
4582 apply_operation_to_constrained_dofs(
const std::vector<unsigned int> &,
4583 const VectorStruct1 &,
4587 template <
typename Number>
4589 apply_operation_to_constrained_dofs(
4590 const std::vector<unsigned int> &constrained_dofs,
4594 for (
const unsigned int i : constrained_dofs)
4595 dst.local_element(i) = src.local_element(i);
4599 namespace MatrixFreeFunctions
4603 template <
typename,
typename,
typename,
typename,
bool>
4604 struct InterfaceSelector
4608 template <
typename MF,
4612 struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4614 using function_type = void (Container::*)(
4618 const
std::pair<unsigned
int, unsigned
int> &) const;
4622 template <
typename MF,
4626 struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4628 using function_type =
4629 void (Container::*)(const MF &,
4632 const
std::pair<unsigned
int, unsigned
int> &);
4640 template <
typename MF,
4645 class MFWorker :
public MFWorkerInterface
4649 using function_type =
typename MatrixFreeFunctions::
4650 InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4654 MFWorker(
const MF &matrix_free,
4655 const InVector &src,
4657 const bool zero_dst_vector_setting,
4658 const Container &container,
4659 function_type cell_function,
4660 function_type face_function,
4661 function_type boundary_function,
4662 const typename MF::DataAccessOnFaces src_vector_face_access =
4663 MF::DataAccessOnFaces::none,
4664 const typename MF::DataAccessOnFaces dst_vector_face_access =
4665 MF::DataAccessOnFaces::none,
4666 const std::function<
void(
const unsigned int,
const unsigned int)>
4667 &operation_before_loop = {},
4668 const std::function<void(
const unsigned int,
const unsigned int)>
4669 &operation_after_loop = {},
4670 const unsigned int dof_handler_index_pre_post = 0)
4671 : matrix_free(matrix_free)
4672 , container(const_cast<Container &>(container))
4673 , cell_function(cell_function)
4674 , face_function(face_function)
4675 , boundary_function(boundary_function)
4678 , src_data_exchanger(matrix_free,
4679 src_vector_face_access,
4681 , dst_data_exchanger(matrix_free,
4682 dst_vector_face_access,
4685 , zero_dst_vector_setting(zero_dst_vector_setting &&
4686 !src_and_dst_are_same)
4687 , operation_before_loop(operation_before_loop)
4688 , operation_after_loop(operation_after_loop)
4689 , dof_handler_index_pre_post(dof_handler_index_pre_post)
4691 Assert(!has_ghost_elements(dst),
4692 ExcMessage(
"The destination vector passed to the matrix-free "
4693 "loop is ghosted. This is not allowed."));
4698 cell(
const std::pair<unsigned int, unsigned int> &cell_range)
override
4700 if (cell_function !=
nullptr && cell_range.second > cell_range.first)
4701 for (
unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4703 const auto cell_subrange =
4704 matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4706 if (cell_subrange.second <= cell_subrange.first)
4710 cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4715 cell(
const unsigned int range_index)
override
4717 process_range(cell_function,
4718 matrix_free.get_task_info().cell_partition_data_hp_ptr,
4719 matrix_free.get_task_info().cell_partition_data_hp,
4724 face(
const unsigned int range_index)
override
4726 process_range(face_function,
4727 matrix_free.get_task_info().face_partition_data_hp_ptr,
4728 matrix_free.get_task_info().face_partition_data_hp,
4733 boundary(
const unsigned int range_index)
override
4735 process_range(boundary_function,
4736 matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4737 matrix_free.get_task_info().boundary_partition_data_hp,
4743 process_range(
const function_type &fu,
4744 const std::vector<unsigned int> &ptr,
4745 const std::vector<unsigned int> &data,
4746 const unsigned int range_index)
4752 for (
unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4755 (container.*fu)(matrix_free,
4758 std::make_pair(data[2 * i], data[2 * i + 1]));
4770 vector_update_ghosts_start()
override
4772 if (!src_and_dst_are_same)
4773 internal::update_ghost_values_start(src, src_data_exchanger);
4778 vector_update_ghosts_finish()
override
4780 if (!src_and_dst_are_same)
4781 internal::update_ghost_values_finish(src, src_data_exchanger);
4786 vector_compress_start()
override
4788 internal::compress_start(dst, dst_data_exchanger);
4793 vector_compress_finish()
override
4795 internal::compress_finish(dst, dst_data_exchanger);
4796 if (!src_and_dst_are_same)
4797 internal::reset_ghost_values(src, src_data_exchanger);
4802 zero_dst_vector_range(
const unsigned int range_index)
override
4804 if (zero_dst_vector_setting)
4805 internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4809 cell_loop_pre_range(
const unsigned int range_index)
override
4811 if (operation_before_loop)
4814 matrix_free.get_dof_info(dof_handler_index_pre_post);
4821 operation_before_loop,
4828 for (
unsigned int id =
4839 cell_loop_post_range(
const unsigned int range_index)
override
4841 if (operation_after_loop)
4845 const std::vector<unsigned int> &partition_row_index =
4846 matrix_free.get_task_info().partition_row_index;
4848 partition_row_index[partition_row_index.size() - 2] - 1)
4849 apply_operation_to_constrained_dofs(
4850 matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
4855 matrix_free.get_dof_info(dof_handler_index_pre_post);
4862 operation_after_loop,
4869 for (
unsigned int id =
4880 const MF &matrix_free;
4881 Container &container;
4882 function_type cell_function;
4883 function_type face_function;
4884 function_type boundary_function;
4886 const InVector &src;
4888 VectorDataExchange<MF::dimension,
4889 typename MF::value_type,
4890 typename MF::vectorized_value_type>
4892 VectorDataExchange<MF::dimension,
4893 typename MF::value_type,
4894 typename MF::vectorized_value_type>
4896 const bool src_and_dst_are_same;
4897 const bool zero_dst_vector_setting;
4898 const std::function<void(
const unsigned int,
const unsigned int)>
4899 operation_before_loop;
4900 const std::function<void(
const unsigned int,
const unsigned int)>
4901 operation_after_loop;
4902 const unsigned int dof_handler_index_pre_post;
4911 template <
class MF,
typename InVector,
typename OutVector>
4912 struct MFClassWrapper
4914 using function_type =
4915 std::function<void(
const MF &,
4918 const std::pair<unsigned int, unsigned int> &)>;
4920 MFClassWrapper(
const function_type cell,
4921 const function_type face,
4922 const function_type boundary)
4925 , boundary(boundary)
4929 cell_integrator(
const MF &mf,
4931 const InVector &src,
4932 const std::pair<unsigned int, unsigned int> &range)
const
4935 cell(mf, dst, src, range);
4939 face_integrator(
const MF &mf,
4941 const InVector &src,
4942 const std::pair<unsigned int, unsigned int> &range)
const
4945 face(mf, dst, src, range);
4949 boundary_integrator(
4952 const InVector &src,
4953 const std::pair<unsigned int, unsigned int> &range)
const
4956 boundary(mf, dst, src, range);
4959 const function_type cell;
4960 const function_type face;
4961 const function_type boundary;
4968template <
int dim,
typename Number,
typename VectorizedArrayType>
4969template <
typename OutVector,
typename InVector>
4975 const std::pair<unsigned int, unsigned int> &)>
4978 const InVector &src,
4979 const bool zero_dst_vector)
const
4982 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
4985 Wrapper wrap(cell_operation,
nullptr,
nullptr);
4986 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
4996 &Wrapper::cell_integrator,
4997 &Wrapper::face_integrator,
4998 &Wrapper::boundary_integrator);
5000 task_info.loop(worker);
5005template <
int dim,
typename Number,
typename VectorizedArrayType>
5006template <
typename OutVector,
typename InVector>
5012 const std::pair<unsigned int, unsigned int> &)>
5015 const InVector &src,
5016 const std::function<
void(
const unsigned int,
const unsigned int)>
5017 &operation_before_loop,
5018 const std::function<
void(
const unsigned int,
const unsigned int)>
5019 &operation_after_loop,
5020 const unsigned int dof_handler_index_pre_post)
const
5023 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5026 Wrapper wrap(cell_operation,
nullptr,
nullptr);
5027 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5037 &Wrapper::cell_integrator,
5038 &Wrapper::face_integrator,
5039 &Wrapper::boundary_integrator,
5040 DataAccessOnFaces::none,
5041 DataAccessOnFaces::none,
5042 operation_before_loop,
5043 operation_after_loop,
5044 dof_handler_index_pre_post);
5046 task_info.loop(worker);
5051template <
int dim,
typename Number,
typename VectorizedArrayType>
5052template <
typename OutVector,
typename InVector>
5058 const std::pair<unsigned int, unsigned int> &)>
5063 const std::pair<unsigned int, unsigned int> &)>
5064 &inner_face_operation,
5068 const std::pair<unsigned int, unsigned int> &)>
5069 &boundary_face_operation,
5071 const InVector &src,
5072 const bool zero_dst_vector,
5073 const DataAccessOnFaces dst_vector_face_access,
5074 const DataAccessOnFaces src_vector_face_access)
const
5077 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5080 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5081 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5091 &Wrapper::cell_integrator,
5092 &Wrapper::face_integrator,
5093 &Wrapper::boundary_integrator,
5094 src_vector_face_access,
5095 dst_vector_face_access);
5097 task_info.loop(worker);
5102template <
int dim,
typename Number,
typename VectorizedArrayType>
5103template <
typename CLASS,
typename OutVector,
typename InVector>
5106 void (CLASS::*function_pointer)(
5107 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5110 const
std::pair<unsigned
int, unsigned
int> &) const,
5111 const CLASS *owning_class,
5113 const InVector &src,
5114 const
bool zero_dst_vector) const
5116 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5134template <
int dim,
typename Number,
typename VectorizedArrayType>
5135template <
typename CLASS,
typename OutVector,
typename InVector>
5138 void (CLASS::*function_pointer)(
5139 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5142 const
std::pair<unsigned
int, unsigned
int> &) const,
5143 const CLASS *owning_class,
5145 const InVector &src,
5147 &operation_before_loop,
5149 &operation_after_loop,
5150 const unsigned
int dof_handler_index_pre_post) const
5152 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5167 operation_before_loop,
5168 operation_after_loop,
5169 dof_handler_index_pre_post);
5175template <
int dim,
typename Number,
typename VectorizedArrayType>
5176template <
typename CLASS,
typename OutVector,
typename InVector>
5179 void (CLASS::*cell_operation)(
5180 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5183 const
std::pair<unsigned
int, unsigned
int> &) const,
5184 void (CLASS::*inner_face_operation)(
5185 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5188 const
std::pair<unsigned
int, unsigned
int> &) const,
5189 void (CLASS::*boundary_face_operation)(
5190 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5193 const
std::pair<unsigned
int, unsigned
int> &) const,
5194 const CLASS *owning_class,
5196 const InVector &src,
5197 const
bool zero_dst_vector,
5198 const DataAccessOnFaces dst_vector_face_access,
5199 const DataAccessOnFaces src_vector_face_access) const
5201 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5212 inner_face_operation,
5213 boundary_face_operation,
5214 src_vector_face_access,
5215 dst_vector_face_access);
5221template <
int dim,
typename Number,
typename VectorizedArrayType>
5222template <
typename CLASS,
typename OutVector,
typename InVector>
5225 void (CLASS::*function_pointer)(
5226 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5229 const
std::pair<unsigned
int, unsigned
int> &),
5230 CLASS *owning_class,
5232 const InVector &src,
5233 const
bool zero_dst_vector) const
5235 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5253template <
int dim,
typename Number,
typename VectorizedArrayType>
5254template <
typename CLASS,
typename OutVector,
typename InVector>
5257 void (CLASS::*function_pointer)(
5258 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5261 const
std::pair<unsigned
int, unsigned
int> &),
5262 CLASS *owning_class,
5264 const InVector &src,
5266 &operation_before_loop,
5268 &operation_after_loop,
5269 const unsigned
int dof_handler_index_pre_post) const
5271 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5286 operation_before_loop,
5287 operation_after_loop,
5288 dof_handler_index_pre_post);
5294template <
int dim,
typename Number,
typename VectorizedArrayType>
5295template <
typename CLASS,
typename OutVector,
typename InVector>
5298 void (CLASS::*cell_operation)(
5299 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5302 const
std::pair<unsigned
int, unsigned
int> &),
5303 void (CLASS::*inner_face_operation)(
5304 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5307 const
std::pair<unsigned
int, unsigned
int> &),
5308 void (CLASS::*boundary_face_operation)(
5309 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5312 const
std::pair<unsigned
int, unsigned
int> &),
5313 CLASS *owning_class,
5315 const InVector &src,
5316 const
bool zero_dst_vector,
5317 const DataAccessOnFaces dst_vector_face_access,
5318 const DataAccessOnFaces src_vector_face_access) const
5320 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5331 inner_face_operation,
5332 boundary_face_operation,
5333 src_vector_face_access,
5334 dst_vector_face_access);
5340template <
int dim,
typename Number,
typename VectorizedArrayType>
5341template <
typename OutVector,
typename InVector>
5347 const std::pair<unsigned int, unsigned int> &)>
5352 const std::pair<unsigned int, unsigned int> &)>
5353 &inner_face_operation,
5357 const std::pair<unsigned int, unsigned int> &)>
5358 &boundary_face_operation,
5360 const InVector &src,
5361 const std::function<
void(
const unsigned int,
const unsigned int)>
5362 &operation_before_loop,
5363 const std::function<
void(
const unsigned int,
const unsigned int)>
5364 &operation_after_loop,
5365 const unsigned int dof_handler_index_pre_post,
5366 const DataAccessOnFaces dst_vector_face_access,
5367 const DataAccessOnFaces src_vector_face_access)
const
5370 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5373 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5374 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5384 &Wrapper::cell_integrator,
5385 &Wrapper::face_integrator,
5386 &Wrapper::boundary_integrator,
5387 src_vector_face_access,
5388 dst_vector_face_access,
5389 operation_before_loop,
5390 operation_after_loop,
5391 dof_handler_index_pre_post);
5393 task_info.loop(worker);
5398template <
int dim,
typename Number,
typename VectorizedArrayType>
5399template <
typename CLASS,
typename OutVector,
typename InVector>
5402 void (CLASS::*cell_operation)(const
MatrixFree &,
5405 const
std::pair<unsigned
int, unsigned
int> &)
5407 void (CLASS::*inner_face_operation)(
5411 const
std::pair<unsigned
int, unsigned
int> &) const,
5412 void (CLASS::*boundary_face_operation)(
5416 const
std::pair<unsigned
int, unsigned
int> &) const,
5417 const CLASS *owning_class,
5419 const InVector &src,
5421 &operation_before_loop,
5423 &operation_after_loop,
5424 const unsigned
int dof_handler_index_pre_post,
5425 const DataAccessOnFaces dst_vector_face_access,
5426 const DataAccessOnFaces src_vector_face_access) const
5428 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5439 inner_face_operation,
5440 boundary_face_operation,
5441 src_vector_face_access,
5442 dst_vector_face_access,
5443 operation_before_loop,
5444 operation_after_loop,
5445 dof_handler_index_pre_post);
5451template <
int dim,
typename Number,
typename VectorizedArrayType>
5452template <
typename CLASS,
typename OutVector,
typename InVector>
5455 void (CLASS::*cell_operation)(const
MatrixFree &,
5458 const
std::pair<unsigned
int, unsigned
int> &),
5459 void (CLASS::*inner_face_operation)(
5463 const
std::pair<unsigned
int, unsigned
int> &),
5464 void (CLASS::*boundary_face_operation)(
5468 const
std::pair<unsigned
int, unsigned
int> &),
5469 const CLASS *owning_class,
5471 const InVector &src,
5473 &operation_before_loop,
5475 &operation_after_loop,
5476 const unsigned
int dof_handler_index_pre_post,
5477 const DataAccessOnFaces dst_vector_face_access,
5478 const DataAccessOnFaces src_vector_face_access) const
5480 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5491 inner_face_operation,
5492 boundary_face_operation,
5493 src_vector_face_access,
5494 dst_vector_face_access,
5495 operation_before_loop,
5496 operation_after_loop,
5497 dof_handler_index_pre_post);
5503template <
int dim,
typename Number,
typename VectorizedArrayType>
5504template <
typename CLASS,
typename OutVector,
typename InVector>
5507 void (CLASS::*function_pointer)(
5508 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5511 const
std::pair<unsigned
int, unsigned
int> &) const,
5512 const CLASS *owning_class,
5514 const InVector &src,
5515 const
bool zero_dst_vector,
5516 const DataAccessOnFaces src_vector_face_access) const
5518 auto src_vector_face_access_temp = src_vector_face_access;
5524 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5537 src_vector_face_access_temp,
5544template <
int dim,
typename Number,
typename VectorizedArrayType>
5545template <
typename CLASS,
typename OutVector,
typename InVector>
5548 void (CLASS::*function_pointer)(
5549 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5552 const
std::pair<unsigned
int, unsigned
int> &),
5553 CLASS *owning_class,
5555 const InVector &src,
5556 const
bool zero_dst_vector,
5557 const DataAccessOnFaces src_vector_face_access) const
5559 auto src_vector_face_access_temp = src_vector_face_access;
5565 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5578 src_vector_face_access_temp,
5585template <
int dim,
typename Number,
typename VectorizedArrayType>
5586template <
typename OutVector,
typename InVector>
5592 const std::pair<unsigned int, unsigned int> &)>
5595 const InVector &src,
5596 const bool zero_dst_vector,
5597 const DataAccessOnFaces src_vector_face_access)
const
5599 auto src_vector_face_access_temp = src_vector_face_access;
5600 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5601 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5602 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5603 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5606 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5609 Wrapper wrap(cell_operation,
nullptr,
nullptr);
5611 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5621 &Wrapper::cell_integrator,
5622 &Wrapper::face_integrator,
5623 &Wrapper::boundary_integrator,
5624 src_vector_face_access_temp,
5625 DataAccessOnFaces::none);
5626 task_info.loop(worker);
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
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
void initialize_dof_handlers(const std::vector< const DoFHandler< dim, dim > * > &dof_handlers, const AdditionalData &additional_data)
types::boundary_id get_boundary_id(const unsigned int face_batch_index) const
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > > shape_info
void print(std::ostream &out) const
void update_mapping(const std::shared_ptr< hp::MappingCollection< dim > > &mapping)
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
unsigned int n_inner_face_batches() const
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
void update_mapping(const Mapping< dim > &mapping)
const Quadrature< dim > & get_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_index) const
unsigned int get_mg_level() const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*boundary_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void print_memory_consumption(StreamType &out) const
bool mapping_initialized() const
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
void internal_reinit(const std::shared_ptr< hp::MappingCollection< dim > > &mapping, const std::vector< const DoFHandler< dim, dim > * > &dof_handlers, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const std::vector< hp::QCollection< q_dim > > &quad, const AdditionalData &additional_data)
unsigned int get_n_q_points_face(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &inner_face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_face_operation, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
bool at_irregular_cell(const unsigned int cell_batch_index) const
~MatrixFree() override=default
const AffineConstraints< Number > & get_affine_constraints(const unsigned int dof_handler_index=0) const
void copy_from(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free_base)
void initialize_face_data_vector(AlignedVector< T > &vec) const
unsigned int get_dofs_per_cell(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int n_constraint_pool_entries() const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
const Number * constraint_pool_begin(const unsigned int pool_index) const
AlignedVector< Number > * acquire_scratch_data_non_threadsafe() const
const IndexSet & get_locally_owned_set(const unsigned int dof_handler_index=0) const
unsigned int get_dofs_per_face(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
internal::MatrixFreeFunctions::TaskInfo task_info
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*boundary_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
bool mapping_is_initialized
void loop_cell_centric(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void initialize_dof_vector(LinearAlgebra::distributed::Vector< Number2 > &vec, const unsigned int dof_handler_index=0) const
DoFHandler< dim >::cell_iterator get_cell_iterator(const unsigned int cell_batch_index, const unsigned int lane_index, const unsigned int dof_handler_index=0) const
MatrixFree(const MatrixFree< dim, Number, VectorizedArrayType > &other)
std::pair< unsigned int, unsigned int > create_cell_subrange_hp(const std::pair< unsigned int, unsigned int > &range, const unsigned int fe_degree, const unsigned int dof_handler_index=0) const
std::pair< unsigned int, unsigned int > create_cell_subrange_hp_by_index(const std::pair< unsigned int, unsigned int > &range, const unsigned int fe_index, const unsigned int dof_handler_index=0) const
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
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*boundary_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
std::list< std::pair< bool, AlignedVector< Number > > > scratch_pad_non_threadsafe
const Number * constraint_pool_end(const unsigned int pool_index) const
void loop_cell_centric(const std::function< void(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
std::vector< SmartPointer< const DoFHandler< dim > > > dof_handlers
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &inner_face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_face_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void initialize_indices(const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const AdditionalData &additional_data)
std::vector< SmartPointer< const AffineConstraints< Number > > > affine_constraints
unsigned int n_components() const
unsigned int n_physical_cells() const
void initialize_dof_vector(VectorType &vec, const unsigned int dof_handler_index=0) const
unsigned int get_n_q_points(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
void reinit(const MappingType &mapping, const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
std::vector< unsigned int > mf_cell_indices
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*boundary_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
std::vector< internal::MatrixFreeFunctions::DoFInfo > dof_info
void initialize_cell_data_vector(AlignedVector< T > &vec) const
unsigned int n_ghost_inner_face_batches() const
void release_scratch_data_non_threadsafe(const AlignedVector< Number > *memory) const
void renumber_dofs(std::vector< types::global_dof_index > &renumbering, const unsigned int dof_handler_index=0)
unsigned int get_face_active_fe_index(const std::pair< unsigned int, unsigned int > range, const bool is_interior_face=true) const
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
bool indices_are_initialized
internal::MatrixFreeFunctions::FaceInfo< VectorizedArrayType::size()> face_info
void reinit(const MappingType &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< number2 > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
unsigned int get_cell_range_category(const std::pair< unsigned int, unsigned int > cell_batch_range) const
void reinit(const MappingType &mapping, const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< QuadratureType > &quad, const AdditionalData &additional_data=AdditionalData())
std::array< types::boundary_id, VectorizedArrayType::size()> get_faces_by_cells_boundary_id(const unsigned int cell_batch_index, const unsigned int face_number) const
std::vector< unsigned int > constraint_pool_row_index
const std::vector< unsigned int > & get_constrained_dofs(const unsigned int dof_handler_index=0) const
unsigned int cell_level_index_end_local
unsigned int n_base_elements(const unsigned int dof_handler_index) const
static bool is_supported(const FiniteElement< dim, spacedim > &fe)
static constexpr unsigned int dimension
A class that provides a separate storage location on each thread that accesses the object.
#define DEAL_II_NAMESPACE_OPEN
#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