16#ifndef dealii_matrix_free_h
17#define dealii_matrix_free_h
111 typename Number = double,
116 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
117 "Type of Number and of VectorizedArrayType do not match.");
573 template <
typename QuadratureType,
typename number2,
typename MappingType>
578 const QuadratureType &quad,
602 template <
typename QuadratureType,
typename number2,
typename MappingType>
607 const std::vector<QuadratureType> &quad,
617 template <
typename QuadratureType,
typename number2,
typename MappingType>
622 const QuadratureType &quad,
779 template <
typename OutVector,
typename InVector>
785 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
788 const bool zero_dst_vector =
false)
const;
836 template <
typename CLASS,
typename OutVector,
typename InVector>
842 const
std::pair<unsigned
int, unsigned
int> &) const,
843 const CLASS *owning_class,
846 const
bool zero_dst_vector = false) const;
851 template <
typename CLASS,
typename OutVector,
typename InVector>
857 const
std::pair<unsigned
int, unsigned
int> &),
861 const
bool zero_dst_vector = false) const;
947 template <
typename CLASS,
typename OutVector,
typename InVector>
953 const
std::pair<unsigned
int, unsigned
int> &) const,
954 const CLASS *owning_class,
957 const
std::function<void(const unsigned
int, const unsigned
int)>
958 &operation_before_loop,
959 const
std::function<void(const unsigned
int, const unsigned
int)>
960 &operation_after_loop,
961 const unsigned
int dof_handler_index_pre_post = 0) const;
966 template <
typename CLASS,
typename OutVector,
typename InVector>
972 const
std::pair<unsigned
int, unsigned
int> &),
976 const
std::function<void(const unsigned
int, const unsigned
int)>
977 &operation_before_loop,
978 const
std::function<void(const unsigned
int, const unsigned
int)>
979 &operation_after_loop,
980 const unsigned
int dof_handler_index_pre_post = 0) const;
986 template <
typename OutVector,
typename InVector>
992 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
995 const std::function<
void(
const unsigned int,
const unsigned int)>
996 &operation_before_loop,
997 const std::function<
void(
const unsigned int,
const unsigned int)>
998 &operation_after_loop,
999 const unsigned int dof_handler_index_pre_post = 0)
const;
1076 template <
typename OutVector,
typename InVector>
1079 const std::function<
1083 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1084 const std::function<
void(
1088 const std::pair<unsigned int, unsigned int> &)> &inner_face_operation,
1089 const std::function<
void(
1093 const std::pair<unsigned int, unsigned int> &)> &boundary_face_operation,
1095 const InVector &src,
1096 const bool zero_dst_vector =
false,
1190 template <
typename CLASS,
typename OutVector,
typename InVector>
1192 loop(
void (CLASS::*cell_operation)(
1196 const
std::pair<unsigned
int, unsigned
int> &) const,
1197 void (CLASS::*inner_face_operation)(
1201 const
std::pair<unsigned
int, unsigned
int> &) const,
1202 void (CLASS::*boundary_face_operation)(
1206 const
std::pair<unsigned
int, unsigned
int> &) const,
1207 const CLASS *owning_class,
1209 const InVector &src,
1210 const
bool zero_dst_vector = false,
1219 template <
typename CLASS,
typename OutVector,
typename InVector>
1221 loop(
void (CLASS::*cell_operation)(
1225 const
std::pair<unsigned
int, unsigned
int> &),
1226 void (CLASS::*inner_face_operation)(
1230 const
std::pair<unsigned
int, unsigned
int> &),
1231 void (CLASS::*boundary_face_operation)(
1235 const
std::pair<unsigned
int, unsigned
int> &),
1236 CLASS *owning_class,
1238 const InVector &src,
1239 const
bool zero_dst_vector = false,
1366 template <
typename CLASS,
typename OutVector,
typename InVector>
1368 loop(
void (CLASS::*cell_operation)(
1372 const
std::pair<unsigned
int, unsigned
int> &) const,
1373 void (CLASS::*inner_face_operation)(
1377 const
std::pair<unsigned
int, unsigned
int> &) const,
1378 void (CLASS::*boundary_face_operation)(
1382 const
std::pair<unsigned
int, unsigned
int> &) const,
1383 const CLASS *owning_class,
1385 const InVector &src,
1386 const
std::function<void(const unsigned
int, const unsigned
int)>
1387 &operation_before_loop,
1388 const
std::function<void(const unsigned
int, const unsigned
int)>
1389 &operation_after_loop,
1390 const unsigned
int dof_handler_index_pre_post = 0,
1399 template <
typename CLASS,
typename OutVector,
typename InVector>
1401 loop(
void (CLASS::*cell_operation)(
1405 const
std::pair<unsigned
int, unsigned
int> &),
1406 void (CLASS::*inner_face_operation)(
1410 const
std::pair<unsigned
int, unsigned
int> &),
1411 void (CLASS::*boundary_face_operation)(
1415 const
std::pair<unsigned
int, unsigned
int> &),
1416 const CLASS *owning_class,
1418 const InVector &src,
1419 const
std::function<void(const unsigned
int, const unsigned
int)>
1420 &operation_before_loop,
1421 const
std::function<void(const unsigned
int, const unsigned
int)>
1422 &operation_after_loop,
1423 const unsigned
int dof_handler_index_pre_post = 0,
1434 template <
typename OutVector,
typename InVector>
1437 const std::function<
1441 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
1442 const std::function<
void(
1446 const std::pair<unsigned int, unsigned int> &)> &inner_face_operation,
1447 const std::function<
void(
1451 const std::pair<unsigned int, unsigned int> &)> &boundary_face_operation,
1453 const InVector &src,
1454 const std::function<
void(
const unsigned int,
const unsigned int)>
1455 &operation_before_loop,
1456 const std::function<
void(
const unsigned int,
const unsigned int)>
1457 &operation_after_loop,
1458 const unsigned int dof_handler_index_pre_post = 0,
1528 template <
typename CLASS,
typename OutVector,
typename InVector>
1534 const
std::pair<unsigned
int, unsigned
int> &) const,
1535 const CLASS *owning_class,
1537 const InVector &src,
1538 const
bool zero_dst_vector = false,
1545 template <
typename CLASS,
typename OutVector,
typename InVector>
1551 const
std::pair<unsigned
int, unsigned
int> &),
1552 CLASS *owning_class,
1554 const InVector &src,
1555 const
bool zero_dst_vector = false,
1562 template <
typename OutVector,
typename InVector>
1568 const std::pair<unsigned int, unsigned int> &)>
1571 const InVector &src,
1572 const bool zero_dst_vector =
false,
1583 std::pair<unsigned int, unsigned int>
1585 const unsigned int fe_degree,
1586 const unsigned int dof_handler_index = 0)
const;
1594 std::pair<unsigned int, unsigned int>
1596 const std::pair<unsigned int, unsigned int> &range,
1597 const unsigned int fe_index,
1598 const unsigned int dof_handler_index = 0)
const;
1611 const std::pair<unsigned int, unsigned int> range,
1619 const std::pair<unsigned int, unsigned int> range,
1620 const bool is_interior_face =
true,
1634 template <
typename T>
1643 template <
typename T>
1660 template <
typename VectorType>
1663 const unsigned int dof_handler_index = 0)
const;
1682 template <
typename Number2>
1685 const unsigned int dof_handler_index = 0)
const;
1697 const std::shared_ptr<const Utilities::MPI::Partitioner> &
1721 const std::vector<unsigned int> &
1736 const unsigned int dof_handler_index = 0);
1747 template <
int spacedim>
1846 const unsigned int face_number)
const;
1878 const unsigned int lane_index,
1879 const unsigned int dof_handler_index = 0)
const;
1888 const unsigned int lane_index)
const;
1912 std::pair<typename DoFHandler<dim>::cell_iterator,
unsigned int>
1914 const unsigned int lane_index,
1915 const bool interior =
true,
1916 const unsigned int fe_component = 0)
const;
1962 const unsigned int hp_active_fe_index = 0)
const;
1969 const unsigned int hp_active_fe_index = 0)
const;
1977 const unsigned int hp_active_fe_index = 0)
const;
1985 const unsigned int hp_active_fe_index = 0)
const;
1992 const unsigned int hp_active_fe_index = 0)
const;
1999 const unsigned int hp_active_fe_index = 0)
const;
2015 const std::pair<unsigned int, unsigned int> cell_batch_range,
2022 std::pair<unsigned int, unsigned int>
2024 const std::pair<unsigned int, unsigned int> face_batch_range,
2040 const unsigned int cell_batch_index,
2047 std::pair<unsigned int, unsigned int>
2049 const unsigned int face_batch_index,
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,
2258 std::vector<ObserverPointer<const AffineConstraints<Number>>>
2265 std::vector<internal::MatrixFreeFunctions::DoFInfo>
dof_info;
2347 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>>
2354 mutable std::list<std::pair<bool, AlignedVector<Number>>>
2377template <
int dim,
typename Number,
typename VectorizedArrayType>
2378template <
typename T>
2383 vec.
resize(this->n_cell_batches() + this->n_ghost_cell_batches());
2388template <
int dim,
typename Number,
typename VectorizedArrayType>
2389template <
typename T>
2394 vec.
resize(this->n_inner_face_batches() + this->n_boundary_face_batches() +
2395 this->n_ghost_inner_face_batches());
2400template <
int dim,
typename Number,
typename VectorizedArrayType>
2401template <
typename VectorType>
2405 const unsigned int comp)
const
2408 "This function is not supported for block vectors.");
2410 Assert(task_info.n_procs == 1,
2411 ExcMessage(
"This function can only be used in serial."));
2414 vec.reinit(dof_info[comp].vector_partitioner->size());
2419template <
int dim,
typename Number,
typename VectorizedArrayType>
2420template <
typename Number2>
2424 const unsigned int comp)
const
2427 vec.
reinit(dof_info[comp].vector_partitioner, task_info.communicator_sm);
2432template <
int dim,
typename Number,
typename VectorizedArrayType>
2433inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
2435 const unsigned int comp)
const
2438 return dof_info[comp].vector_partitioner;
2443template <
int dim,
typename Number,
typename VectorizedArrayType>
2444inline const std::vector<unsigned int> &
2446 const unsigned int comp)
const
2449 return dof_info[comp].constrained_dofs;
2454template <
int dim,
typename Number,
typename VectorizedArrayType>
2459 return dof_handlers.size();
2464template <
int dim,
typename Number,
typename VectorizedArrayType>
2467 const unsigned int dof_no)
const
2471 return dof_handlers[dof_no]->get_fe().n_base_elements();
2476template <
int dim,
typename Number,
typename VectorizedArrayType>
2485template <
int dim,
typename Number,
typename VectorizedArrayType>
2489 return task_info.n_active_cells;
2494template <
int dim,
typename Number,
typename VectorizedArrayType>
2498 return *(task_info.cell_partition_data.end() - 2);
2503template <
int dim,
typename Number,
typename VectorizedArrayType>
2507 return *(task_info.cell_partition_data.end() - 1) -
2508 *(task_info.cell_partition_data.end() - 2);
2513template <
int dim,
typename Number,
typename VectorizedArrayType>
2517 if (task_info.face_partition_data.empty())
2519 return task_info.face_partition_data.back();
2524template <
int dim,
typename Number,
typename VectorizedArrayType>
2528 if (task_info.face_partition_data.empty())
2530 return task_info.boundary_partition_data.back() -
2531 task_info.face_partition_data.back();
2536template <
int dim,
typename Number,
typename VectorizedArrayType>
2540 if (task_info.face_partition_data.empty())
2542 return face_info.faces.size() - task_info.boundary_partition_data.back();
2547template <
int dim,
typename Number,
typename VectorizedArrayType>
2550 const unsigned int face_batch_index)
const
2552 Assert(face_batch_index >= task_info.boundary_partition_data[0] &&
2553 face_batch_index < task_info.boundary_partition_data.back(),
2555 task_info.boundary_partition_data[0],
2556 task_info.boundary_partition_data.back()));
2562template <
int dim,
typename Number,
typename VectorizedArrayType>
2565 const unsigned int cell_batch_index,
2566 const unsigned int face_number)
const
2570 Assert(face_info.cell_and_face_boundary_id.size(0) >= n_cell_batches(),
2574 for (
unsigned int v = 0;
2575 v < n_active_entries_per_cell_batch(cell_batch_index);
2578 face_info.cell_and_face_boundary_id(cell_batch_index, face_number, v);
2584template <
int dim,
typename Number,
typename VectorizedArrayType>
2585inline const internal::MatrixFreeFunctions::
2586 MappingInfo<dim, Number, VectorizedArrayType> &
2589 return mapping_info;
2594template <
int dim,
typename Number,
typename VectorizedArrayType>
2597 const unsigned int dof_index)
const
2600 return dof_info[dof_index];
2605template <
int dim,
typename Number,
typename VectorizedArrayType>
2609 return constraint_pool_row_index.size() - 1;
2614template <
int dim,
typename Number,
typename VectorizedArrayType>
2615inline const Number *
2617 const unsigned int row)
const
2620 return constraint_pool_data.empty() ?
2622 constraint_pool_data.data() + constraint_pool_row_index[row];
2627template <
int dim,
typename Number,
typename VectorizedArrayType>
2628inline const Number *
2630 const unsigned int row)
const
2633 return constraint_pool_data.empty() ?
2635 constraint_pool_data.data() + constraint_pool_row_index[row + 1];
2640template <
int dim,
typename Number,
typename VectorizedArrayType>
2641inline std::pair<unsigned int, unsigned int>
2643 const std::pair<unsigned int, unsigned int> &range,
2644 const unsigned int degree,
2645 const unsigned int dof_handler_component)
const
2647 if (dof_info[dof_handler_component].cell_active_fe_index.empty())
2650 dof_info[dof_handler_component].fe_index_conversion.size(), 1);
2652 dof_info[dof_handler_component].fe_index_conversion[0].
size(), 1);
2653 if (dof_info[dof_handler_component].fe_index_conversion[0][0] == degree)
2656 return {range.second, range.second};
2660 dof_info[dof_handler_component].fe_index_from_degree(0, degree);
2661 if (fe_index >= dof_info[dof_handler_component].max_fe_index)
2662 return {range.second, range.second};
2664 return create_cell_subrange_hp_by_index(range,
2666 dof_handler_component);
2671template <
int dim,
typename Number,
typename VectorizedArrayType>
2674 const unsigned int cell_batch_index)
const
2677 return VectorizedArrayType::size() > 1 &&
2678 cell_level_index[(cell_batch_index + 1) * VectorizedArrayType::size() -
2679 1] == cell_level_index[(cell_batch_index + 1) *
2680 VectorizedArrayType::size() -
2686template <
int dim,
typename Number,
typename VectorizedArrayType>
2690 return shape_info.size(2);
2694template <
int dim,
typename Number,
typename VectorizedArrayType>
2697 const std::pair<unsigned int, unsigned int> range,
2698 const unsigned int dof_handler_index)
const
2700 const unsigned int dof_no =
2702 first_hp_dof_handler_index :
2705 const auto &fe_indices = dof_info[dof_no].cell_active_fe_index;
2707 if (fe_indices.empty() ==
true ||
2708 dof_handlers[dof_no]->get_fe_collection().size() == 1)
2711 const auto index = fe_indices[range.first];
2713 for (
unsigned int i = range.first; i < range.second; ++i)
2721template <
int dim,
typename Number,
typename VectorizedArrayType>
2724 const std::pair<unsigned int, unsigned int> range,
2725 const bool is_interior_face,
2726 const unsigned int dof_handler_index)
const
2728 const unsigned int dof_no =
2730 first_hp_dof_handler_index :
2733 const auto &fe_indices = dof_info[dof_no].cell_active_fe_index;
2735 if (fe_indices.empty() ==
true)
2738 if (is_interior_face)
2740 const unsigned int index =
2741 fe_indices[face_info.faces[range.first].cells_interior[0] /
2742 VectorizedArrayType::size()];
2744 for (
unsigned int i = range.first; i < range.second; ++i)
2746 fe_indices[face_info.faces[i].cells_interior[0] /
2747 VectorizedArrayType::size()]);
2753 const unsigned int index =
2754 fe_indices[face_info.faces[range.first].cells_exterior[0] /
2755 VectorizedArrayType::size()];
2757 for (
unsigned int i = range.first; i < range.second; ++i)
2759 fe_indices[face_info.faces[i].cells_exterior[0] /
2760 VectorizedArrayType::size()]);
2768template <
int dim,
typename Number,
typename VectorizedArrayType>
2771 const unsigned int cell_batch_index)
const
2775 const std::vector<unsigned char> &n_lanes_filled =
2776 dof_info[0].n_vectorization_lanes_filled
2780 return n_lanes_filled[cell_batch_index];
2785template <
int dim,
typename Number,
typename VectorizedArrayType>
2788 const unsigned int face_batch_index)
const
2792 const std::vector<unsigned char> &n_lanes_filled =
2793 dof_info[0].n_vectorization_lanes_filled
2796 return n_lanes_filled[face_batch_index];
2801template <
int dim,
typename Number,
typename VectorizedArrayType>
2804 const unsigned int dof_handler_index,
2805 const unsigned int active_fe_index)
const
2807 return dof_info[dof_handler_index].dofs_per_cell[active_fe_index];
2812template <
int dim,
typename Number,
typename VectorizedArrayType>
2815 const unsigned int quad_index,
2816 const unsigned int active_fe_index)
const
2819 return mapping_info.cell_data[quad_index]
2820 .descriptor[active_fe_index]
2826template <
int dim,
typename Number,
typename VectorizedArrayType>
2829 const unsigned int dof_handler_index,
2830 const unsigned int active_fe_index)
const
2832 return dof_info[dof_handler_index].dofs_per_face[active_fe_index];
2837template <
int dim,
typename Number,
typename VectorizedArrayType>
2840 const unsigned int quad_index,
2841 const unsigned int active_fe_index)
const
2844 return mapping_info.face_data[quad_index]
2845 .descriptor[active_fe_index]
2851template <
int dim,
typename Number,
typename VectorizedArrayType>
2854 const unsigned int dof_handler_index)
const
2856 return dof_info[dof_handler_index].vector_partitioner->locally_owned_range();
2861template <
int dim,
typename Number,
typename VectorizedArrayType>
2864 const unsigned int dof_handler_index)
const
2866 return dof_info[dof_handler_index].vector_partitioner->ghost_indices();
2871template <
int dim,
typename Number,
typename VectorizedArrayType>
2874 const unsigned int dof_handler_index,
2875 const unsigned int index_quad,
2876 const unsigned int index_fe,
2877 const unsigned int active_fe_index,
2878 const unsigned int active_quad_index)
const
2881 const unsigned int ind =
2882 dof_info[dof_handler_index].global_base_element_offset + index_fe;
2887 return shape_info(ind, index_quad, active_fe_index, active_quad_index);
2892template <
int dim,
typename Number,
typename VectorizedArrayType>
2894 VectorizedArrayType::size()> &
2896 const unsigned int face_batch_index)
const
2899 return face_info.faces[face_batch_index];
2904template <
int dim,
typename Number,
typename VectorizedArrayType>
2909 return face_info.cell_and_face_to_plain_faces;
2914template <
int dim,
typename Number,
typename VectorizedArrayType>
2917 const unsigned int quad_index,
2918 const unsigned int active_fe_index)
const
2921 return mapping_info.cell_data[quad_index]
2922 .descriptor[active_fe_index]
2928template <
int dim,
typename Number,
typename VectorizedArrayType>
2931 const unsigned int quad_index,
2932 const unsigned int active_fe_index)
const
2935 return mapping_info.face_data[quad_index]
2936 .descriptor[active_fe_index]
2942template <
int dim,
typename Number,
typename VectorizedArrayType>
2945 const std::pair<unsigned int, unsigned int> range,
2946 const unsigned int dof_handler_index)
const
2948 auto result = get_cell_category(range.first, dof_handler_index);
2950 for (
unsigned int i = range.first; i < range.second; ++i)
2951 result =
std::max(result, get_cell_category(i, dof_handler_index));
2958template <
int dim,
typename Number,
typename VectorizedArrayType>
2959inline std::pair<unsigned int, unsigned int>
2961 const std::pair<unsigned int, unsigned int> range,
2962 const unsigned int dof_handler_index)
const
2964 auto result = get_face_category(range.first, dof_handler_index);
2966 for (
unsigned int i = range.first; i < range.second; ++i)
2969 std::max(result.first, get_face_category(i, dof_handler_index).
first);
2971 std::max(result.second, get_face_category(i, dof_handler_index).
second);
2979template <
int dim,
typename Number,
typename VectorizedArrayType>
2982 const unsigned int cell_batch_index,
2983 const unsigned int dof_handler_index)
const
2987 const unsigned int dof_no =
2989 first_hp_dof_handler_index :
2993 dof_info[dof_no].cell_active_fe_index.size());
2994 if (dof_info[dof_no].cell_active_fe_index.empty())
2997 return dof_info[dof_no].cell_active_fe_index[cell_batch_index];
3002template <
int dim,
typename Number,
typename VectorizedArrayType>
3003inline std::pair<unsigned int, unsigned int>
3005 const unsigned int face_batch_index,
3006 const unsigned int dof_handler_index)
const
3008 const unsigned int dof_no =
3010 first_hp_dof_handler_index :
3014 if (dof_info[dof_no].cell_active_fe_index.empty())
3015 return std::make_pair(0U, 0U);
3017 std::pair<unsigned int, unsigned int> result = std::make_pair(0U, 0U);
3018 for (
unsigned int v = 0;
3019 v < VectorizedArrayType::size() &&
3020 face_info.faces[face_batch_index].cells_interior[v] !=
3025 dof_info[dof_no].cell_active_fe_index[face_info.faces[face_batch_index]
3026 .cells_interior[v] /
3027 VectorizedArrayType::size()]);
3028 if (face_info.faces[face_batch_index].cells_exterior[0] !=
3030 for (
unsigned int v = 0;
3031 v < VectorizedArrayType::size() &&
3032 face_info.faces[face_batch_index].cells_exterior[v] !=
3037 dof_info[dof_no].cell_active_fe_index[face_info.faces[face_batch_index]
3038 .cells_exterior[v] /
3039 VectorizedArrayType::size()]);
3047template <
int dim,
typename Number,
typename VectorizedArrayType>
3051 return indices_are_initialized;
3056template <
int dim,
typename Number,
typename VectorizedArrayType>
3060 return mapping_is_initialized;
3064template <
int dim,
typename Number,
typename VectorizedArrayType>
3073template <
int dim,
typename Number,
typename VectorizedArrayType>
3078 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3079 list_type &
data = scratch_pad.get();
3080 for (
typename list_type::iterator it =
data.begin(); it !=
data.end(); ++it)
3081 if (it->first ==
false)
3087 return &
data.front().second;
3092template <
int dim,
typename Number,
typename VectorizedArrayType>
3098 std::list<std::pair<bool, AlignedVector<VectorizedArrayType>>>;
3099 list_type &
data = scratch_pad.get();
3100 for (
typename list_type::iterator it =
data.begin(); it !=
data.end(); ++it)
3101 if (&it->second == scratch)
3112template <
int dim,
typename Number,
typename VectorizedArrayType>
3118 scratch_pad_non_threadsafe.
begin();
3119 it != scratch_pad_non_threadsafe.end();
3121 if (it->first ==
false)
3126 scratch_pad_non_threadsafe.push_front(
3128 return &scratch_pad_non_threadsafe.front().second;
3133template <
int dim,
typename Number,
typename VectorizedArrayType>
3140 scratch_pad_non_threadsafe.
begin();
3141 it != scratch_pad_non_threadsafe.end();
3143 if (&it->second == scratch)
3158 namespace MatrixFreeImplementation
3160 template <
int dim,
int spacedim>
3161 inline std::vector<IndexSet>
3162 extract_locally_owned_index_sets(
3164 const unsigned int level)
3166 std::vector<IndexSet> locally_owned_set;
3167 locally_owned_set.reserve(dofh.size());
3168 for (
unsigned int j = 0; j < dofh.size(); ++j)
3170 locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
3172 locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(
level));
3173 return locally_owned_set;
3180template <
int dim,
typename Number,
typename VectorizedArrayType>
3181template <
typename QuadratureType,
typename number2,
typename MappingType>
3184 const MappingType &mapping,
3187 const QuadratureType &quad,
3191 std::vector<const DoFHandler<dim, dim> *> dof_handlers;
3192 std::vector<const AffineConstraints<number2> *> constraints;
3194 dof_handlers.push_back(&dof_handler);
3195 constraints.push_back(&constraints_in);
3197 std::vector<IndexSet> locally_owned_sets =
3198 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3199 dof_handlers, additional_data.
mg_level);
3201 std::vector<hp::QCollection<dim>> quad_hp;
3202 quad_hp.emplace_back(quad);
3214template <
int dim,
typename Number,
typename VectorizedArrayType>
3215template <
typename QuadratureType,
typename number2,
typename MappingType>
3218 const MappingType &mapping,
3221 const QuadratureType &quad,
3225 std::vector<IndexSet> locally_owned_set =
3226 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3227 dof_handler, additional_data.
mg_level);
3228 std::vector<hp::QCollection<dim>> quad_hp;
3229 quad_hp.emplace_back(quad);
3241template <
int dim,
typename Number,
typename VectorizedArrayType>
3242template <
typename QuadratureType,
typename number2,
typename MappingType>
3245 const MappingType &mapping,
3248 const std::vector<QuadratureType> &quad,
3252 std::vector<IndexSet> locally_owned_set =
3253 internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
3254 dof_handler, additional_data.
mg_level);
3255 std::vector<hp::QCollection<dim>> quad_hp;
3256 for (
unsigned int q = 0; q < quad.size(); ++q)
3257 quad_hp.emplace_back(quad[q]);
3284 template <
int dim,
typename Number,
typename VectorizedArrayType>
3285 struct VectorDataExchange
3292 static constexpr unsigned int channel_shift = 20;
3301 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
3302 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
3303 DataAccessOnFaces vector_face_access,
3304 const unsigned
int n_components)
3305 : matrix_free(matrix_free)
3306 , vector_face_access(
3307 matrix_free.get_task_info().face_partition_data.empty() ?
3308 ::
MatrixFree<dim, Number, VectorizedArrayType>::
3309 DataAccessOnFaces::unspecified :
3311 , ghosts_were_set(false)
3312# ifdef DEAL_II_WITH_MPI
3313 , tmp_data(n_components)
3314 , requests(n_components)
3318 if (this->vector_face_access !=
3320 DataAccessOnFaces::unspecified)
3321 for (unsigned
int c = 0; c < matrix_free.n_components(); ++c)
3323 matrix_free.get_dof_info(c).vector_exchanger_face_variants.size(),
3332 ~VectorDataExchange()
3334# ifdef DEAL_II_WITH_MPI
3335 for (
unsigned int i = 0; i < tmp_data.size(); ++i)
3336 if (tmp_data[i] !=
nullptr)
3337 matrix_free.release_scratch_data_non_threadsafe(tmp_data[i]);
3347 template <
typename VectorType>
3349 find_vector_in_mf(
const VectorType &vec,
3350 const bool check_global_compatibility =
true)
const
3353 for (
unsigned int c = 0; c < matrix_free.n_components(); ++c)
3354 if (vec.get_partitioner().get() ==
3355 matrix_free.get_dof_info(c).vector_partitioner.get())
3359 for (
unsigned int c = 0; c < matrix_free.n_components(); ++c)
3360 if (check_global_compatibility ?
3361 vec.get_partitioner()->is_globally_compatible(
3362 *matrix_free.get_dof_info(c).vector_partitioner) :
3363 vec.get_partitioner()->is_compatible(
3364 *matrix_free.get_dof_info(c).vector_partitioner))
3380 get_partitioner(
const unsigned int mf_component)
const
3383 .vector_exchanger_face_variants.size(),
3385 if (vector_face_access ==
3387 DataAccessOnFaces::
none)
3388 return *matrix_free.get_dof_info(mf_component)
3389 .vector_exchanger_face_variants[0];
3390 else if (vector_face_access ==
3392 DataAccessOnFaces::
values)
3393 return *matrix_free.get_dof_info(mf_component)
3394 .vector_exchanger_face_variants[1];
3395 else if (vector_face_access ==
3398 return *matrix_free.get_dof_info(mf_component)
3399 .vector_exchanger_face_variants[2];
3400 else if (vector_face_access ==
3402 DataAccessOnFaces::values_all_faces)
3403 return *matrix_free.get_dof_info(mf_component)
3404 .vector_exchanger_face_variants[3];
3405 else if (vector_face_access ==
3407 DataAccessOnFaces::gradients_all_faces)
3408 return *matrix_free.get_dof_info(mf_component)
3409 .vector_exchanger_face_variants[4];
3411 return *matrix_free.get_dof_info(mf_component).vector_exchanger.get();
3420 std::enable_if_t<is_not_parallel_vector<VectorType>,
VectorType>
3423 update_ghost_values_start(
const unsigned int ,
3424 const VectorType & )
3433 std::enable_if_t<!has_update_ghost_values_start<VectorType> &&
3434 !is_not_parallel_vector<VectorType>,
3437 update_ghost_values_start(
const unsigned int component_in_block_vector,
3438 const VectorType &vec)
3440 (void)component_in_block_vector;
3441 const bool ghosts_set = vec.has_ghost_elements();
3443 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3444 ghosts_set ==
false,
3449 ghosts_were_set =
true;
3453 vec.update_ghost_values();
3464 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3465 !has_exchange_on_subset<VectorType>,
3468 update_ghost_values_start(
const unsigned int component_in_block_vector,
3469 const VectorType &vec)
3471 (void)component_in_block_vector;
3472 const bool ghosts_set = vec.has_ghost_elements();
3474 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3475 ghosts_set ==
false,
3480 ghosts_were_set =
true;
3484 vec.update_ghost_values_start(component_in_block_vector + channel_shift);
3496 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3497 has_exchange_on_subset<VectorType>,
3500 update_ghost_values_start(
const unsigned int component_in_block_vector,
3501 const VectorType &vec)
3503 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3504 "Type mismatch between VectorType and VectorDataExchange");
3505 (void)component_in_block_vector;
3506 const bool ghosts_set = vec.has_ghost_elements();
3508 Assert(matrix_free.get_task_info().allow_ghosted_vectors_in_loops ||
3509 ghosts_set ==
false,
3514 ghosts_were_set =
true;
3518 if (vec.size() != 0)
3520# ifdef DEAL_II_WITH_MPI
3521 const unsigned int mf_component = find_vector_in_mf(vec);
3523 const auto &part = get_partitioner(mf_component);
3525 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3526 part.n_import_sm_procs() == 0)
3529 tmp_data[component_in_block_vector] =
3530 matrix_free.acquire_scratch_data_non_threadsafe();
3531 tmp_data[component_in_block_vector]->resize_fast(
3532 part.n_import_indices());
3535 part.export_to_ghosted_array_start(
3536 component_in_block_vector * 2 + channel_shift,
3538 vec.shared_vector_data(),
3540 part.locally_owned_size(),
3541 matrix_free.get_dof_info(mf_component)
3542 .vector_partitioner->n_ghost_indices()),
3544 part.n_import_indices()),
3545 this->requests[component_in_block_vector]);
3557 std::enable_if_t<!has_update_ghost_values_start<VectorType>,
3560 update_ghost_values_finish(
const unsigned int ,
3561 const VectorType & )
3572 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3573 !has_exchange_on_subset<VectorType>,
3576 update_ghost_values_finish(
const unsigned int component_in_block_vector,
3577 const VectorType &vec)
3579 (void)component_in_block_vector;
3581 if (ghosts_were_set)
3584 vec.update_ghost_values_finish();
3596 std::enable_if_t<has_update_ghost_values_start<VectorType> &&
3597 has_exchange_on_subset<VectorType>,
3600 update_ghost_values_finish(
const unsigned int component_in_block_vector,
3601 const VectorType &vec)
3603 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3604 "Type mismatch between VectorType and VectorDataExchange");
3605 (void)component_in_block_vector;
3607 if (ghosts_were_set)
3610 if (vec.size() != 0)
3612# ifdef DEAL_II_WITH_MPI
3616 const unsigned int mf_component = find_vector_in_mf(vec);
3618 const auto &part = get_partitioner(mf_component);
3620 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3621 part.n_import_sm_procs() != 0)
3623 part.export_to_ghosted_array_finish(
3625 vec.shared_vector_data(),
3627 part.locally_owned_size(),
3628 matrix_free.get_dof_info(mf_component)
3629 .vector_partitioner->n_ghost_indices()),
3630 this->requests[component_in_block_vector]);
3632 matrix_free.release_scratch_data_non_threadsafe(
3633 tmp_data[component_in_block_vector]);
3634 tmp_data[component_in_block_vector] =
nullptr;
3640 vec.set_ghost_state(
true);
3649 std::enable_if_t<is_not_parallel_vector<VectorType>,
VectorType>
3652 compress_start(
const unsigned int ,
3663 std::enable_if_t<!has_compress_start<VectorType> &&
3664 !is_not_parallel_vector<VectorType>,
3667 compress_start(
const unsigned int component_in_block_vector,
3670 (void)component_in_block_vector;
3683 std::enable_if_t<has_compress_start<VectorType> &&
3684 !has_exchange_on_subset<VectorType>,
3687 compress_start(
const unsigned int component_in_block_vector,
3690 (void)component_in_block_vector;
3692 vec.compress_start(component_in_block_vector + channel_shift);
3704 std::enable_if_t<has_compress_start<VectorType> &&
3705 has_exchange_on_subset<VectorType>,
3708 compress_start(
const unsigned int component_in_block_vector,
3711 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3712 "Type mismatch between VectorType and VectorDataExchange");
3713 (void)component_in_block_vector;
3716 if (vec.size() != 0)
3718# ifdef DEAL_II_WITH_MPI
3719 const unsigned int mf_component = find_vector_in_mf(vec);
3721 const auto &part = get_partitioner(mf_component);
3723 if (part.n_ghost_indices() == 0 && part.n_import_indices() == 0 &&
3724 part.n_import_sm_procs() == 0)
3727 tmp_data[component_in_block_vector] =
3728 matrix_free.acquire_scratch_data_non_threadsafe();
3729 tmp_data[component_in_block_vector]->resize_fast(
3730 part.n_import_indices());
3733 part.import_from_ghosted_array_start(
3735 component_in_block_vector * 2 + channel_shift,
3737 vec.shared_vector_data(),
3739 matrix_free.get_dof_info(mf_component)
3740 .vector_partitioner->n_ghost_indices()),
3742 part.n_import_indices()),
3743 this->requests[component_in_block_vector]);
3756 std::enable_if_t<!has_compress_start<VectorType>,
VectorType> * =
nullptr>
3758 compress_finish(
const unsigned int ,
3770 std::enable_if_t<has_compress_start<VectorType> &&
3771 !has_exchange_on_subset<VectorType>,
3774 compress_finish(
const unsigned int component_in_block_vector,
3777 (void)component_in_block_vector;
3790 std::enable_if_t<has_compress_start<VectorType> &&
3791 has_exchange_on_subset<VectorType>,
3794 compress_finish(
const unsigned int component_in_block_vector,
3797 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3798 "Type mismatch between VectorType and VectorDataExchange");
3799 (void)component_in_block_vector;
3800 if (vec.size() != 0)
3802# ifdef DEAL_II_WITH_MPI
3806 const unsigned int mf_component = find_vector_in_mf(vec);
3808 const auto &part = get_partitioner(mf_component);
3810 if (part.n_ghost_indices() != 0 || part.n_import_indices() != 0 ||
3811 part.n_import_sm_procs() != 0)
3813 part.import_from_ghosted_array_finish(
3816 vec.shared_vector_data(),
3818 matrix_free.get_dof_info(mf_component)
3819 .vector_partitioner->n_ghost_indices()),
3821 tmp_data[component_in_block_vector]->begin(),
3822 part.n_import_indices()),
3823 this->requests[component_in_block_vector]);
3825 matrix_free.release_scratch_data_non_threadsafe(
3826 tmp_data[component_in_block_vector]);
3827 tmp_data[component_in_block_vector] =
nullptr;
3833 MPI_Barrier(matrix_free.get_task_info().communicator_sm);
3846 std::enable_if_t<is_not_parallel_vector<VectorType>,
VectorType>
3849 reset_ghost_values(
const VectorType & )
const
3859 std::enable_if_t<!has_exchange_on_subset<VectorType> &&
3860 !is_not_parallel_vector<VectorType>,
3863 reset_ghost_values(
const VectorType &vec)
const
3865 if (ghosts_were_set ==
true)
3868 vec.zero_out_ghost_values();
3879 std::enable_if_t<has_exchange_on_subset<VectorType>,
VectorType>
3882 reset_ghost_values(
const VectorType &vec)
const
3884 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3885 "Type mismatch between VectorType and VectorDataExchange");
3886 if (ghosts_were_set ==
true)
3889 if (vec.size() != 0)
3891# ifdef DEAL_II_WITH_MPI
3894 const unsigned int mf_component = find_vector_in_mf(vec);
3896 const auto &part = get_partitioner(mf_component);
3898 if (part.n_ghost_indices() > 0)
3900 part.reset_ghost_values(
3902 part.locally_owned_size(),
3903 matrix_free.get_dof_info(mf_component)
3904 .vector_partitioner->n_ghost_indices()));
3910 vec.set_ghost_state(
false);
3921 std::enable_if_t<has_exchange_on_subset<VectorType>,
VectorType>
3924 zero_vector_region(
const unsigned int range_index, VectorType &vec)
const
3926 static_assert(std::is_same_v<Number, typename VectorType::value_type>,
3927 "Type mismatch between VectorType and VectorDataExchange");
3932 const unsigned int mf_component = find_vector_in_mf(vec,
false);
3934 matrix_free.get_dof_info(mf_component);
3942 for (
unsigned int id =
3962 std::enable_if_t<!has_exchange_on_subset<VectorType>,
VectorType>
3964 std::enable_if_t<has_assignment_operator<VectorType>,
VectorType>
3967 zero_vector_region(
const unsigned int range_index, VectorType &vec)
const
3971 if constexpr (std::is_same_v<
3975 for (
unsigned int i = 0; i < vec.size(); ++i)
3976 vec[i] =
typename VectorType::value_type();
3979 vec =
typename VectorType::value_type();
3990 zero_vector_region(
const unsigned int, ...)
const
3994 "Zeroing is only implemented for vector types "
3995 "which provide operator=(const VectorType::value_type)"));
4000 const ::MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
4001 const typename ::MatrixFree<dim, Number, VectorizedArrayType>::
4002 DataAccessOnFaces vector_face_access;
4003 bool ghosts_were_set;
4004# ifdef DEAL_II_WITH_MPI
4005 std::vector<AlignedVector<Number> *> tmp_data;
4006 std::vector<std::vector<MPI_Request>> requests;
4010 template <
typename VectorStruct>
4012 n_components(
const VectorStruct &vec);
4014 template <
typename VectorStruct>
4016 n_components_block(
const VectorStruct &vec,
const std::bool_constant<true>)
4018 unsigned int components = 0;
4019 for (
unsigned int bl = 0; bl < vec.n_blocks(); ++bl)
4020 components += n_components(vec.block(bl));
4024 template <
typename VectorStruct>
4026 n_components_block(
const VectorStruct &,
const std::bool_constant<false>)
4031 template <
typename VectorStruct>
4033 n_components(
const VectorStruct &vec)
4035 return n_components_block(
4039 template <
typename VectorStruct>
4041 n_components(
const std::vector<VectorStruct> &vec)
4043 unsigned int components = 0;
4044 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4045 components += n_components_block(
4050 template <
typename VectorStruct>
4052 n_components(
const std::vector<VectorStruct *> &vec)
4054 unsigned int components = 0;
4055 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4056 components += n_components_block(
4068 template <
typename VectorStruct,
4069 std::enable_if_t<!has_communication_block_size<VectorStruct>,
4070 VectorStruct> * =
nullptr>
4071 constexpr unsigned int
4072 get_communication_block_size(
const VectorStruct &)
4079 template <
typename VectorStruct,
4080 std::enable_if_t<has_communication_block_size<VectorStruct>,
4081 VectorStruct> * =
nullptr>
4082 constexpr unsigned int
4083 get_communication_block_size(
const VectorStruct &)
4085 return VectorStruct::communication_block_size;
4091 std::enable_if_t<is_not_parallel_vector<VectorType>,
VectorType> * =
4094 has_ghost_elements(
const VectorType &vec)
4103 std::enable_if_t<!is_not_parallel_vector<VectorType>,
VectorType>
4106 has_ghost_elements(
const VectorType &vec)
4108 return vec.has_ghost_elements();
4123 typename VectorStruct,
4125 typename VectorizedArrayType,
4126 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4129 update_ghost_values_start(
4130 const VectorStruct &vec,
4131 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4132 const unsigned int channel = 0)
4134 if (get_communication_block_size(vec) < vec.n_blocks())
4136 const bool ghosts_set = vec.has_ghost_elements();
4138 Assert(exchanger.matrix_free.get_task_info()
4139 .allow_ghosted_vectors_in_loops ||
4140 ghosts_set ==
false,
4145 exchanger.ghosts_were_set =
true;
4149 vec.update_ghost_values();
4153 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4154 update_ghost_values_start(vec.block(i), exchanger, channel + i);
4162 typename VectorStruct,
4164 typename VectorizedArrayType,
4165 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4168 update_ghost_values_start(
4169 const VectorStruct &vec,
4170 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4171 const unsigned int channel = 0)
4173 exchanger.update_ghost_values_start(channel, vec);
4180 typename VectorStruct,
4182 typename VectorizedArrayType>
4184 update_ghost_values_start(
4185 const std::vector<VectorStruct> &vec,
4186 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4188 unsigned int component_index = 0;
4189 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4191 update_ghost_values_start(vec[comp], exchanger, component_index);
4192 component_index += n_components(vec[comp]);
4200 typename VectorStruct,
4202 typename VectorizedArrayType>
4204 update_ghost_values_start(
4205 const std::vector<VectorStruct *> &vec,
4206 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4208 unsigned int component_index = 0;
4209 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4211 update_ghost_values_start(*vec[comp], exchanger, component_index);
4212 component_index += n_components(*vec[comp]);
4224 typename VectorStruct,
4226 typename VectorizedArrayType,
4227 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4230 update_ghost_values_finish(
4231 const VectorStruct &vec,
4232 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4233 const unsigned int channel = 0)
4235 if (get_communication_block_size(vec) < vec.n_blocks())
4241 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4242 update_ghost_values_finish(vec.block(i), exchanger, channel + i);
4249 typename VectorStruct,
4251 typename VectorizedArrayType,
4252 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4255 update_ghost_values_finish(
4256 const VectorStruct &vec,
4257 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4258 const unsigned int channel = 0)
4260 exchanger.update_ghost_values_finish(channel, vec);
4267 typename VectorStruct,
4269 typename VectorizedArrayType>
4271 update_ghost_values_finish(
4272 const std::vector<VectorStruct> &vec,
4273 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4275 unsigned int component_index = 0;
4276 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4278 update_ghost_values_finish(vec[comp], exchanger, component_index);
4279 component_index += n_components(vec[comp]);
4287 typename VectorStruct,
4289 typename VectorizedArrayType>
4291 update_ghost_values_finish(
4292 const std::vector<VectorStruct *> &vec,
4293 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4295 unsigned int component_index = 0;
4296 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4298 update_ghost_values_finish(*vec[comp], exchanger, component_index);
4299 component_index += n_components(*vec[comp]);
4311 typename VectorStruct,
4313 typename VectorizedArrayType,
4314 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4319 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4320 const unsigned int channel = 0)
4322 if (get_communication_block_size(vec) < vec.n_blocks())
4325 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4326 compress_start(vec.block(i), exchanger, channel + i);
4333 typename VectorStruct,
4335 typename VectorizedArrayType,
4336 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4341 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4342 const unsigned int channel = 0)
4344 exchanger.compress_start(channel, vec);
4351 typename VectorStruct,
4353 typename VectorizedArrayType>
4356 std::vector<VectorStruct> &vec,
4357 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4359 unsigned int component_index = 0;
4360 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4362 compress_start(vec[comp], exchanger, component_index);
4363 component_index += n_components(vec[comp]);
4371 typename VectorStruct,
4373 typename VectorizedArrayType>
4376 std::vector<VectorStruct *> &vec,
4377 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4379 unsigned int component_index = 0;
4380 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4382 compress_start(*vec[comp], exchanger, component_index);
4383 component_index += n_components(*vec[comp]);
4395 typename VectorStruct,
4397 typename VectorizedArrayType,
4398 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4403 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4404 const unsigned int channel = 0)
4406 if (get_communication_block_size(vec) < vec.n_blocks())
4412 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4413 compress_finish(vec.block(i), exchanger, channel + i);
4420 typename VectorStruct,
4422 typename VectorizedArrayType,
4423 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4428 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger,
4429 const unsigned int channel = 0)
4431 exchanger.compress_finish(channel, vec);
4438 typename VectorStruct,
4440 typename VectorizedArrayType>
4443 std::vector<VectorStruct> &vec,
4444 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4446 unsigned int component_index = 0;
4447 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4449 compress_finish(vec[comp], exchanger, component_index);
4450 component_index += n_components(vec[comp]);
4458 typename VectorStruct,
4460 typename VectorizedArrayType>
4463 std::vector<VectorStruct *> &vec,
4464 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4466 unsigned int component_index = 0;
4467 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4469 compress_finish(*vec[comp], exchanger, component_index);
4470 component_index += n_components(*vec[comp]);
4486 typename VectorStruct,
4488 typename VectorizedArrayType,
4489 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4493 const VectorStruct &vec,
4494 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4497 if (exchanger.ghosts_were_set ==
true)
4500 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4501 reset_ghost_values(vec.block(i), exchanger);
4508 typename VectorStruct,
4510 typename VectorizedArrayType,
4511 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4515 const VectorStruct &vec,
4516 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4519 if (exchanger.ghosts_were_set ==
true)
4522 exchanger.reset_ghost_values(vec);
4529 typename VectorStruct,
4531 typename VectorizedArrayType>
4534 const std::vector<VectorStruct> &vec,
4535 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4538 if (exchanger.ghosts_were_set ==
true)
4541 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4542 reset_ghost_values(vec[comp], exchanger);
4549 typename VectorStruct,
4551 typename VectorizedArrayType>
4554 const std::vector<VectorStruct *> &vec,
4555 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4558 if (exchanger.ghosts_were_set ==
true)
4561 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4562 reset_ghost_values(*vec[comp], exchanger);
4573 typename VectorStruct,
4575 typename VectorizedArrayType,
4576 std::enable_if_t<IsBlockVector<VectorStruct>::value, VectorStruct>
4580 const unsigned int range_index,
4582 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4584 for (
unsigned int i = 0; i < vec.n_blocks(); ++i)
4585 exchanger.zero_vector_region(range_index, vec.block(i));
4592 typename VectorStruct,
4594 typename VectorizedArrayType,
4595 std::enable_if_t<!IsBlockVector<VectorStruct>::value, VectorStruct>
4599 const unsigned int range_index,
4601 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4603 exchanger.zero_vector_region(range_index, vec);
4610 typename VectorStruct,
4612 typename VectorizedArrayType>
4615 const unsigned int range_index,
4616 std::vector<VectorStruct> &vec,
4617 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4619 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4620 zero_vector_region(range_index, vec[comp], exchanger);
4627 typename VectorStruct,
4629 typename VectorizedArrayType>
4632 const unsigned int range_index,
4633 std::vector<VectorStruct *> &vec,
4634 VectorDataExchange<dim, Number, VectorizedArrayType> &exchanger)
4636 for (
unsigned int comp = 0; comp < vec.size(); ++comp)
4637 zero_vector_region(range_index, *vec[comp], exchanger);
4645 template <
typename VectorStruct1,
typename VectorStruct2>
4647 apply_operation_to_constrained_dofs(
const std::vector<unsigned int> &,
4648 const VectorStruct1 &,
4652 template <
typename Number>
4654 apply_operation_to_constrained_dofs(
4655 const std::vector<unsigned int> &constrained_dofs,
4659 for (
const unsigned int i : constrained_dofs)
4660 dst.local_element(i) = src.local_element(i);
4664 namespace MatrixFreeFunctions
4668 template <
typename,
typename,
typename,
typename,
bool>
4669 struct InterfaceSelector
4673 template <
typename MF,
4677 struct InterfaceSelector<MF, InVector, OutVector, Container, true>
4679 using function_type = void (Container::*)(
4683 const
std::pair<unsigned
int, unsigned
int> &) const;
4687 template <
typename MF,
4691 struct InterfaceSelector<MF, InVector, OutVector, Container, false>
4693 using function_type =
4694 void (Container::*)(const MF &,
4697 const
std::pair<unsigned
int, unsigned
int> &);
4705 template <
typename MF,
4710 class MFWorker :
public MFWorkerInterface
4714 using function_type =
typename MatrixFreeFunctions::
4715 InterfaceSelector<MF, InVector, OutVector, Container, is_constant>::
4719 MFWorker(
const MF &matrix_free,
4720 const InVector &src,
4722 const bool zero_dst_vector_setting,
4723 const Container &container,
4724 function_type cell_function,
4725 function_type face_function,
4726 function_type boundary_function,
4727 const typename MF::DataAccessOnFaces src_vector_face_access =
4728 MF::DataAccessOnFaces::none,
4729 const typename MF::DataAccessOnFaces dst_vector_face_access =
4730 MF::DataAccessOnFaces::none,
4731 const std::function<
void(
const unsigned int,
const unsigned int)>
4732 &operation_before_loop = {},
4733 const std::function<void(
const unsigned int,
const unsigned int)>
4734 &operation_after_loop = {},
4735 const unsigned int dof_handler_index_pre_post = 0)
4736 : matrix_free(matrix_free)
4737 , container(const_cast<Container &>(container))
4738 , cell_function(cell_function)
4739 , face_function(face_function)
4740 , boundary_function(boundary_function)
4743 , src_data_exchanger(matrix_free,
4744 src_vector_face_access,
4746 , dst_data_exchanger(matrix_free,
4747 dst_vector_face_access,
4750 , zero_dst_vector_setting(zero_dst_vector_setting &&
4751 !src_and_dst_are_same)
4752 , operation_before_loop(operation_before_loop)
4753 , operation_after_loop(operation_after_loop)
4754 , dof_handler_index_pre_post(dof_handler_index_pre_post)
4756 Assert(!has_ghost_elements(dst),
4757 ExcMessage(
"The destination vector passed to the matrix-free "
4758 "loop is ghosted. This is not allowed."));
4763 cell(
const std::pair<unsigned int, unsigned int> &cell_range)
override
4765 if (cell_function !=
nullptr && cell_range.second > cell_range.first)
4766 for (
unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
4768 const auto cell_subrange =
4769 matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
4771 if (cell_subrange.second <= cell_subrange.first)
4775 cell_function)(matrix_free, this->dst, this->src, cell_subrange);
4780 cell(
const unsigned int range_index)
override
4782 process_range(cell_function,
4783 matrix_free.get_task_info().cell_partition_data_hp_ptr,
4784 matrix_free.get_task_info().cell_partition_data_hp,
4789 face(
const unsigned int range_index)
override
4791 process_range(face_function,
4792 matrix_free.get_task_info().face_partition_data_hp_ptr,
4793 matrix_free.get_task_info().face_partition_data_hp,
4798 boundary(
const unsigned int range_index)
override
4800 process_range(boundary_function,
4801 matrix_free.get_task_info().boundary_partition_data_hp_ptr,
4802 matrix_free.get_task_info().boundary_partition_data_hp,
4808 process_range(
const function_type &fu,
4809 const std::vector<unsigned int> &ptr,
4810 const std::vector<unsigned int> &
data,
4811 const unsigned int range_index)
4817 for (
unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i)
4820 (container.*fu)(matrix_free,
4823 std::make_pair(
data[2 * i],
data[2 * i + 1]));
4835 vector_update_ghosts_start()
override
4837 if (!src_and_dst_are_same)
4838 internal::update_ghost_values_start(src, src_data_exchanger);
4843 vector_update_ghosts_finish()
override
4845 if (!src_and_dst_are_same)
4846 internal::update_ghost_values_finish(src, src_data_exchanger);
4851 vector_compress_start()
override
4853 internal::compress_start(dst, dst_data_exchanger);
4858 vector_compress_finish()
override
4860 internal::compress_finish(dst, dst_data_exchanger);
4861 if (!src_and_dst_are_same)
4862 internal::reset_ghost_values(src, src_data_exchanger);
4867 zero_dst_vector_range(
const unsigned int range_index)
override
4869 if (zero_dst_vector_setting)
4870 internal::zero_vector_region(range_index, dst, dst_data_exchanger);
4874 cell_loop_pre_range(
const unsigned int range_index)
override
4876 if (operation_before_loop)
4879 matrix_free.get_dof_info(dof_handler_index_pre_post);
4886 operation_before_loop,
4893 for (
unsigned int id =
4904 cell_loop_post_range(
const unsigned int range_index)
override
4906 if (operation_after_loop)
4910 const std::vector<unsigned int> &partition_row_index =
4911 matrix_free.get_task_info().partition_row_index;
4913 partition_row_index[partition_row_index.size() - 2] - 1)
4914 apply_operation_to_constrained_dofs(
4915 matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
4920 matrix_free.get_dof_info(dof_handler_index_pre_post);
4927 operation_after_loop,
4934 for (
unsigned int id =
4945 const MF &matrix_free;
4946 Container &container;
4947 function_type cell_function;
4948 function_type face_function;
4949 function_type boundary_function;
4951 const InVector &src;
4953 VectorDataExchange<MF::dimension,
4954 typename MF::value_type,
4955 typename MF::vectorized_value_type>
4957 VectorDataExchange<MF::dimension,
4958 typename MF::value_type,
4959 typename MF::vectorized_value_type>
4961 const bool src_and_dst_are_same;
4962 const bool zero_dst_vector_setting;
4963 const std::function<void(
const unsigned int,
const unsigned int)>
4964 operation_before_loop;
4965 const std::function<void(
const unsigned int,
const unsigned int)>
4966 operation_after_loop;
4967 const unsigned int dof_handler_index_pre_post;
4976 template <
class MF,
typename InVector,
typename OutVector>
4977 struct MFClassWrapper
4979 using function_type =
4980 std::function<void(
const MF &,
4983 const std::pair<unsigned int, unsigned int> &)>;
4985 MFClassWrapper(
const function_type cell,
4986 const function_type face,
4987 const function_type boundary)
4990 , boundary(boundary)
4994 cell_integrator(
const MF &mf,
4996 const InVector &src,
4997 const std::pair<unsigned int, unsigned int> &range)
const
5000 cell(mf, dst, src, range);
5004 face_integrator(
const MF &mf,
5006 const InVector &src,
5007 const std::pair<unsigned int, unsigned int> &range)
const
5010 face(mf, dst, src, range);
5014 boundary_integrator(
5017 const InVector &src,
5018 const std::pair<unsigned int, unsigned int> &range)
const
5021 boundary(mf, dst, src, range);
5024 const function_type cell;
5025 const function_type face;
5026 const function_type boundary;
5033template <
int dim,
typename Number,
typename VectorizedArrayType>
5034template <
typename OutVector,
typename InVector>
5040 const std::pair<unsigned int, unsigned int> &)>
5043 const InVector &src,
5044 const bool zero_dst_vector)
const
5047 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5050 Wrapper wrap(cell_operation,
nullptr,
nullptr);
5051 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5061 &Wrapper::cell_integrator,
5062 &Wrapper::face_integrator,
5063 &Wrapper::boundary_integrator);
5065 task_info.loop(worker);
5070template <
int dim,
typename Number,
typename VectorizedArrayType>
5071template <
typename OutVector,
typename InVector>
5077 const std::pair<unsigned int, unsigned int> &)>
5080 const InVector &src,
5081 const std::function<
void(
const unsigned int,
const unsigned int)>
5082 &operation_before_loop,
5083 const std::function<
void(
const unsigned int,
const unsigned int)>
5084 &operation_after_loop,
5085 const unsigned int dof_handler_index_pre_post)
const
5088 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5091 Wrapper wrap(cell_operation,
nullptr,
nullptr);
5092 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5102 &Wrapper::cell_integrator,
5103 &Wrapper::face_integrator,
5104 &Wrapper::boundary_integrator,
5105 DataAccessOnFaces::none,
5106 DataAccessOnFaces::none,
5107 operation_before_loop,
5108 operation_after_loop,
5109 dof_handler_index_pre_post);
5111 task_info.loop(worker);
5116template <
int dim,
typename Number,
typename VectorizedArrayType>
5117template <
typename OutVector,
typename InVector>
5123 const std::pair<unsigned int, unsigned int> &)>
5128 const std::pair<unsigned int, unsigned int> &)>
5129 &inner_face_operation,
5133 const std::pair<unsigned int, unsigned int> &)>
5134 &boundary_face_operation,
5136 const InVector &src,
5137 const bool zero_dst_vector,
5138 const DataAccessOnFaces dst_vector_face_access,
5139 const DataAccessOnFaces src_vector_face_access)
const
5142 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5145 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5146 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5156 &Wrapper::cell_integrator,
5157 &Wrapper::face_integrator,
5158 &Wrapper::boundary_integrator,
5159 src_vector_face_access,
5160 dst_vector_face_access);
5162 task_info.loop(worker);
5167template <
int dim,
typename Number,
typename VectorizedArrayType>
5168template <
typename CLASS,
typename OutVector,
typename InVector>
5171 void (CLASS::*function_pointer)(
5172 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5175 const
std::pair<unsigned
int, unsigned
int> &) const,
5176 const CLASS *owning_class,
5178 const InVector &src,
5179 const
bool zero_dst_vector) const
5181 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5199template <
int dim,
typename Number,
typename VectorizedArrayType>
5200template <
typename CLASS,
typename OutVector,
typename InVector>
5203 void (CLASS::*function_pointer)(
5204 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5207 const
std::pair<unsigned
int, unsigned
int> &) const,
5208 const CLASS *owning_class,
5210 const InVector &src,
5211 const
std::function<void(const unsigned
int, const unsigned
int)>
5212 &operation_before_loop,
5213 const
std::function<void(const unsigned
int, const unsigned
int)>
5214 &operation_after_loop,
5215 const unsigned
int dof_handler_index_pre_post) const
5217 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5232 operation_before_loop,
5233 operation_after_loop,
5234 dof_handler_index_pre_post);
5240template <
int dim,
typename Number,
typename VectorizedArrayType>
5241template <
typename CLASS,
typename OutVector,
typename InVector>
5244 void (CLASS::*cell_operation)(
5245 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5248 const
std::pair<unsigned
int, unsigned
int> &) const,
5249 void (CLASS::*inner_face_operation)(
5250 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5253 const
std::pair<unsigned
int, unsigned
int> &) const,
5254 void (CLASS::*boundary_face_operation)(
5255 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5258 const
std::pair<unsigned
int, unsigned
int> &) const,
5259 const CLASS *owning_class,
5261 const InVector &src,
5262 const
bool zero_dst_vector,
5263 const DataAccessOnFaces dst_vector_face_access,
5264 const DataAccessOnFaces src_vector_face_access) const
5266 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5277 inner_face_operation,
5278 boundary_face_operation,
5279 src_vector_face_access,
5280 dst_vector_face_access);
5286template <
int dim,
typename Number,
typename VectorizedArrayType>
5287template <
typename CLASS,
typename OutVector,
typename InVector>
5290 void (CLASS::*function_pointer)(
5291 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5294 const
std::pair<unsigned
int, unsigned
int> &),
5295 CLASS *owning_class,
5297 const InVector &src,
5298 const
bool zero_dst_vector) const
5300 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5318template <
int dim,
typename Number,
typename VectorizedArrayType>
5319template <
typename CLASS,
typename OutVector,
typename InVector>
5322 void (CLASS::*function_pointer)(
5323 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5326 const
std::pair<unsigned
int, unsigned
int> &),
5327 CLASS *owning_class,
5329 const InVector &src,
5330 const
std::function<void(const unsigned
int, const unsigned
int)>
5331 &operation_before_loop,
5332 const
std::function<void(const unsigned
int, const unsigned
int)>
5333 &operation_after_loop,
5334 const unsigned
int dof_handler_index_pre_post) const
5336 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5351 operation_before_loop,
5352 operation_after_loop,
5353 dof_handler_index_pre_post);
5359template <
int dim,
typename Number,
typename VectorizedArrayType>
5360template <
typename CLASS,
typename OutVector,
typename InVector>
5363 void (CLASS::*cell_operation)(
5364 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5367 const
std::pair<unsigned
int, unsigned
int> &),
5368 void (CLASS::*inner_face_operation)(
5369 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5372 const
std::pair<unsigned
int, unsigned
int> &),
5373 void (CLASS::*boundary_face_operation)(
5374 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5377 const
std::pair<unsigned
int, unsigned
int> &),
5378 CLASS *owning_class,
5380 const InVector &src,
5381 const
bool zero_dst_vector,
5382 const DataAccessOnFaces dst_vector_face_access,
5383 const DataAccessOnFaces src_vector_face_access) const
5385 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5396 inner_face_operation,
5397 boundary_face_operation,
5398 src_vector_face_access,
5399 dst_vector_face_access);
5405template <
int dim,
typename Number,
typename VectorizedArrayType>
5406template <
typename OutVector,
typename InVector>
5412 const std::pair<unsigned int, unsigned int> &)>
5417 const std::pair<unsigned int, unsigned int> &)>
5418 &inner_face_operation,
5422 const std::pair<unsigned int, unsigned int> &)>
5423 &boundary_face_operation,
5425 const InVector &src,
5426 const std::function<
void(
const unsigned int,
const unsigned int)>
5427 &operation_before_loop,
5428 const std::function<
void(
const unsigned int,
const unsigned int)>
5429 &operation_after_loop,
5430 const unsigned int dof_handler_index_pre_post,
5431 const DataAccessOnFaces dst_vector_face_access,
5432 const DataAccessOnFaces src_vector_face_access)
const
5435 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5438 Wrapper wrap(cell_operation, inner_face_operation, boundary_face_operation);
5439 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5449 &Wrapper::cell_integrator,
5450 &Wrapper::face_integrator,
5451 &Wrapper::boundary_integrator,
5452 src_vector_face_access,
5453 dst_vector_face_access,
5454 operation_before_loop,
5455 operation_after_loop,
5456 dof_handler_index_pre_post);
5458 task_info.loop(worker);
5463template <
int dim,
typename Number,
typename VectorizedArrayType>
5464template <
typename CLASS,
typename OutVector,
typename InVector>
5467 void (CLASS::*cell_operation)(const
MatrixFree &,
5470 const
std::pair<unsigned
int, unsigned
int> &)
5472 void (CLASS::*inner_face_operation)(
5476 const
std::pair<unsigned
int, unsigned
int> &) const,
5477 void (CLASS::*boundary_face_operation)(
5481 const
std::pair<unsigned
int, unsigned
int> &) const,
5482 const CLASS *owning_class,
5484 const InVector &src,
5485 const
std::function<void(const unsigned
int, const unsigned
int)>
5486 &operation_before_loop,
5487 const
std::function<void(const unsigned
int, const unsigned
int)>
5488 &operation_after_loop,
5489 const unsigned
int dof_handler_index_pre_post,
5490 const DataAccessOnFaces dst_vector_face_access,
5491 const DataAccessOnFaces src_vector_face_access) const
5493 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5504 inner_face_operation,
5505 boundary_face_operation,
5506 src_vector_face_access,
5507 dst_vector_face_access,
5508 operation_before_loop,
5509 operation_after_loop,
5510 dof_handler_index_pre_post);
5516template <
int dim,
typename Number,
typename VectorizedArrayType>
5517template <
typename CLASS,
typename OutVector,
typename InVector>
5520 void (CLASS::*cell_operation)(const
MatrixFree &,
5523 const
std::pair<unsigned
int, unsigned
int> &),
5524 void (CLASS::*inner_face_operation)(
5528 const
std::pair<unsigned
int, unsigned
int> &),
5529 void (CLASS::*boundary_face_operation)(
5533 const
std::pair<unsigned
int, unsigned
int> &),
5534 const CLASS *owning_class,
5536 const InVector &src,
5537 const
std::function<void(const unsigned
int, const unsigned
int)>
5538 &operation_before_loop,
5539 const
std::function<void(const unsigned
int, const unsigned
int)>
5540 &operation_after_loop,
5541 const unsigned
int dof_handler_index_pre_post,
5542 const DataAccessOnFaces dst_vector_face_access,
5543 const DataAccessOnFaces src_vector_face_access) const
5545 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5556 inner_face_operation,
5557 boundary_face_operation,
5558 src_vector_face_access,
5559 dst_vector_face_access,
5560 operation_before_loop,
5561 operation_after_loop,
5562 dof_handler_index_pre_post);
5568template <
int dim,
typename Number,
typename VectorizedArrayType>
5569template <
typename CLASS,
typename OutVector,
typename InVector>
5572 void (CLASS::*function_pointer)(
5573 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5576 const
std::pair<unsigned
int, unsigned
int> &) const,
5577 const CLASS *owning_class,
5579 const InVector &src,
5580 const
bool zero_dst_vector,
5581 const DataAccessOnFaces src_vector_face_access) const
5583 auto src_vector_face_access_temp = src_vector_face_access;
5589 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5602 src_vector_face_access_temp,
5609template <
int dim,
typename Number,
typename VectorizedArrayType>
5610template <
typename CLASS,
typename OutVector,
typename InVector>
5613 void (CLASS::*function_pointer)(
5614 const
MatrixFree<dim, Number, VectorizedArrayType> &,
5617 const
std::pair<unsigned
int, unsigned
int> &),
5618 CLASS *owning_class,
5620 const InVector &src,
5621 const
bool zero_dst_vector,
5622 const DataAccessOnFaces src_vector_face_access) const
5624 auto src_vector_face_access_temp = src_vector_face_access;
5630 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5643 src_vector_face_access_temp,
5650template <
int dim,
typename Number,
typename VectorizedArrayType>
5651template <
typename OutVector,
typename InVector>
5657 const std::pair<unsigned int, unsigned int> &)>
5660 const InVector &src,
5661 const bool zero_dst_vector,
5662 const DataAccessOnFaces src_vector_face_access)
const
5664 auto src_vector_face_access_temp = src_vector_face_access;
5665 if (DataAccessOnFaces::gradients == src_vector_face_access_temp)
5666 src_vector_face_access_temp = DataAccessOnFaces::gradients_all_faces;
5667 else if (DataAccessOnFaces::values == src_vector_face_access_temp)
5668 src_vector_face_access_temp = DataAccessOnFaces::values_all_faces;
5671 internal::MFClassWrapper<MatrixFree<dim, Number, VectorizedArrayType>,
5674 Wrapper wrap(cell_operation,
nullptr,
nullptr);
5676 internal::MFWorker<MatrixFree<dim, Number, VectorizedArrayType>,
5686 &Wrapper::cell_integrator,
5687 &Wrapper::face_integrator,
5688 &Wrapper::boundary_integrator,
5689 src_vector_face_access_temp,
5690 DataAccessOnFaces::none);
5691 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
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
unsigned int n_ghost_cell_batches() const
void initialize_dof_handlers(const std::vector< const DoFHandler< dim, dim > * > &dof_handlers, const AdditionalData &additional_data)
types::boundary_id get_boundary_id(const unsigned int face_batch_index) const
void print(std::ostream &out) const
void update_mapping(const std::shared_ptr< hp::MappingCollection< dim > > &mapping)
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
unsigned int n_inner_face_batches() const
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
void update_mapping(const Mapping< dim > &mapping)
const Quadrature< dim > & get_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_index) const
unsigned int get_mg_level() const
unsigned int get_cell_category(const unsigned int cell_batch_index, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*boundary_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void print_memory_consumption(StreamType &out) const
std::pair< unsigned int, unsigned int > get_face_range_category(const std::pair< unsigned int, unsigned int > face_batch_range, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) const
bool mapping_initialized() const
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
void internal_reinit(const std::shared_ptr< hp::MappingCollection< dim > > &mapping, const std::vector< const DoFHandler< dim, dim > * > &dof_handlers, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const std::vector< hp::QCollection< q_dim > > &quad, const AdditionalData &additional_data)
unsigned int get_n_q_points_face(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &inner_face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_face_operation, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
bool at_irregular_cell(const unsigned int cell_batch_index) const
~MatrixFree() override=default
const AffineConstraints< Number > & get_affine_constraints(const unsigned int dof_handler_index=0) const
void copy_from(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free_base)
void initialize_face_data_vector(AlignedVector< T > &vec) const
unsigned int get_dofs_per_cell(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int n_constraint_pool_entries() const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
const Number * constraint_pool_begin(const unsigned int pool_index) const
AlignedVector< Number > * acquire_scratch_data_non_threadsafe() const
const IndexSet & get_locally_owned_set(const unsigned int dof_handler_index=0) const
unsigned int get_dofs_per_face(const unsigned int dof_handler_index=0, const unsigned int hp_active_fe_index=0) const
internal::MatrixFreeFunctions::TaskInfo task_info
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
std::vector< ObserverPointer< const DoFHandler< dim > > > dof_handlers
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, void(CLASS::*boundary_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
unsigned int first_hp_dof_handler_index
bool mapping_is_initialized
void loop_cell_centric(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void initialize_dof_vector(LinearAlgebra::distributed::Vector< Number2 > &vec, const unsigned int dof_handler_index=0) const
DoFHandler< dim >::cell_iterator get_cell_iterator(const unsigned int cell_batch_index, const unsigned int lane_index, const unsigned int dof_handler_index=0) const
MatrixFree(const MatrixFree< dim, Number, VectorizedArrayType > &other)
std::pair< unsigned int, unsigned int > create_cell_subrange_hp(const std::pair< unsigned int, unsigned int > &range, const unsigned int fe_degree, const unsigned int dof_handler_index=0) const
std::pair< unsigned int, unsigned int > create_cell_subrange_hp_by_index(const std::pair< unsigned int, unsigned int > &range, const unsigned int fe_index, const unsigned int dof_handler_index=0) const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0) const
std::vector< std::pair< unsigned int, unsigned int > > cell_level_index
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner(const unsigned int dof_handler_index=0) const
std::vector< ObserverPointer< const AffineConstraints< Number > > > affine_constraints
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
std::vector< Number > constraint_pool_data
VectorizedArrayType vectorized_value_type
unsigned int n_cell_batches() const
const IndexSet & get_ghost_set(const unsigned int dof_handler_index=0) const
std::pair< int, int > get_cell_level_and_index(const unsigned int cell_batch_index, const unsigned int lane_index) const
bool indices_initialized() const
const Quadrature< dim - 1 > & get_face_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
unsigned int get_face_active_fe_index(const std::pair< unsigned int, unsigned int > range, const bool is_interior_face=true, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) const
unsigned int get_cell_active_fe_index(const std::pair< unsigned int, unsigned int > range, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) const
Threads::ThreadLocalStorage< std::list< std::pair< bool, AlignedVector< VectorizedArrayType > > > > scratch_pad
unsigned int get_matrix_free_cell_index(const typename Triangulation< dim >::cell_iterator &cell) const
std::size_t memory_consumption() const
std::pair< unsigned int, unsigned int > get_face_category(const unsigned int face_batch_index, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) const
unsigned int n_boundary_face_batches() const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*boundary_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), const CLASS *owning_class, OutVector &dst, const InVector &src, const std::function< void(const unsigned int, const unsigned int)> &operation_before_loop, const std::function< void(const unsigned int, const unsigned int)> &operation_after_loop, const unsigned int dof_handler_index_pre_post=0, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
std::list< std::pair< bool, AlignedVector< Number > > > scratch_pad_non_threadsafe
const Number * constraint_pool_end(const unsigned int pool_index) const
void loop_cell_centric(const std::function< void(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &inner_face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_face_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void initialize_indices(const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< IndexSet > &locally_owned_set, const AdditionalData &additional_data)
unsigned int n_components() const
unsigned int n_physical_cells() const
void initialize_dof_vector(VectorType &vec, const unsigned int dof_handler_index=0) const
unsigned int get_n_q_points(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
void reinit(const MappingType &mapping, const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
std::vector< unsigned int > mf_cell_indices
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
unsigned int get_cell_range_category(const std::pair< unsigned int, unsigned int > cell_batch_range, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) const
void loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*inner_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), void(CLASS::*boundary_face_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void cell_loop(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &), CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
std::vector< internal::MatrixFreeFunctions::DoFInfo > dof_info
void initialize_cell_data_vector(AlignedVector< T > &vec) const
unsigned int n_ghost_inner_face_batches() const
void release_scratch_data_non_threadsafe(const AlignedVector< Number > *memory) const
void renumber_dofs(std::vector< types::global_dof_index > &renumbering, const unsigned int dof_handler_index=0)
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
bool indices_are_initialized
internal::MatrixFreeFunctions::FaceInfo< VectorizedArrayType::size()> face_info
void reinit(const MappingType &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< number2 > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
void reinit(const MappingType &mapping, const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const AffineConstraints< number2 > * > &constraint, const std::vector< QuadratureType > &quad, const AdditionalData &additional_data=AdditionalData())
std::array< types::boundary_id, VectorizedArrayType::size()> get_faces_by_cells_boundary_id(const unsigned int cell_batch_index, const unsigned int face_number) const
std::vector< unsigned int > constraint_pool_row_index
const internal::MatrixFreeFunctions::ShapeInfo< Number > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
const std::vector< unsigned int > & get_constrained_dofs(const unsigned int dof_handler_index=0) const
unsigned int cell_level_index_end_local
unsigned int n_base_elements(const unsigned int dof_handler_index) const
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< Number > > shape_info
static bool is_supported(const FiniteElement< dim, spacedim > &fe)
static constexpr unsigned int dimension
A class that provides a separate storage location on each thread that accesses the object.
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ActiveSelector::cell_iterator cell_iterator
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_default
No update.
std::vector< index_type > data
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
VectorType::value_type * begin(VectorType &V)
unsigned int minimum_parallel_grain_size
constexpr unsigned int invalid_unsigned_int
constexpr types::boundary_id invalid_boundary_id
void apply_to_subranges(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, const Function &f, const unsigned int grainsize)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
unsigned short int fe_index
TasksParallelScheme tasks_parallel_scheme
bool hold_all_faces_to_owned_cells
UpdateFlags mapping_update_flags_inner_faces
AdditionalData(const TasksParallelScheme tasks_parallel_scheme=partition_partition, const unsigned int tasks_block_size=0, const UpdateFlags mapping_update_flags=update_gradients|update_JxW_values, const UpdateFlags mapping_update_flags_boundary_faces=update_default, const UpdateFlags mapping_update_flags_inner_faces=update_default, const UpdateFlags mapping_update_flags_faces_by_cells=update_default, const unsigned int mg_level=numbers::invalid_unsigned_int, const bool store_plain_indices=true, const bool initialize_indices=true, const bool initialize_mapping=true, const bool overlap_communication_computation=true, const bool hold_all_faces_to_owned_cells=false, const bool cell_vectorization_categories_strict=false, const bool allow_ghosted_vectors_in_loops=true)
std::vector< unsigned int > cell_vectorization_category
bool cell_vectorization_categories_strict
UpdateFlags mapping_update_flags_boundary_faces
UpdateFlags mapping_update_flags
UpdateFlags mapping_update_flags_faces_by_cells
AdditionalData & operator=(const AdditionalData &other)=default
AdditionalData(const AdditionalData &other)
bool allow_ghosted_vectors_in_loops
unsigned int tasks_block_size
bool overlap_communication_computation
@ dof_access_face_interior
std::vector< unsigned int > cell_loop_pre_list_index
std::vector< unsigned int > cell_loop_post_list_index
std::vector< std::pair< unsigned int, unsigned int > > vector_zero_range_list
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_pre_list
std::vector< unsigned int > vector_zero_range_list_index
std::vector< std::pair< unsigned int, unsigned int > > cell_loop_post_list
void loop(MFWorkerInterface &worker) const