39 template <
int dim,
typename Number,
bool is_face_>
52 std::function<std::vector<std::unique_ptr<FEEvalType>>(
53 const std::pair<unsigned int, unsigned int> &)>
55 std::function<void(std::vector<std::unique_ptr<FEEvalType>> &,
58 std::function<void(std::vector<std::unique_ptr<FEEvalType>> &)>
68 template <
int dim,
typename AdditionalData>
71 AdditionalData &additional_data);
91 typename VectorizedArrayType,
96 VectorType &diagonal_global,
102 VectorizedArrayType> &)>
104 const unsigned int dof_no = 0,
105 const unsigned int quad_no = 0,
106 const unsigned int first_selected_component = 0,
107 const unsigned int first_vector_component = 0);
118 typename QuadOperation>
123 const QuadOperation &quad_operation,
126 const unsigned int dof_no = 0,
127 const unsigned int quad_no = 0,
128 const unsigned int first_selected_component = 0,
129 const unsigned int first_vector_component = 0);
134 template <
typename CLASS,
140 typename VectorizedArrayType,
145 VectorType &diagonal_global,
151 VectorizedArrayType> &) const,
152 const CLASS *owning_class,
153 const unsigned
int dof_no = 0,
154 const unsigned
int quad_no = 0,
155 const unsigned
int first_selected_component = 0,
156 const unsigned
int first_vector_component = 0);
175 typename VectorizedArrayType,
180 VectorType &diagonal_global,
186 VectorizedArrayType> &)>
193 VectorizedArrayType> &,
199 VectorizedArrayType> &)>
206 VectorizedArrayType> &)>
208 const unsigned int dof_no = 0,
209 const unsigned int quad_no = 0,
210 const unsigned int first_selected_component = 0,
211 const unsigned int first_vector_component = 0);
218 template <
typename CLASS,
224 typename VectorizedArrayType,
229 VectorType &diagonal_global,
235 VectorizedArrayType> &) const,
241 VectorizedArrayType> &,
247 VectorizedArrayType> &)
254 VectorizedArrayType> &)
256 const CLASS *owning_class,
257 const unsigned
int dof_no = 0,
258 const unsigned
int quad_no = 0,
259 const unsigned
int first_selected_component = 0,
260 const unsigned
int first_vector_component = 0);
277 typename VectorizedArrayType,
289 VectorizedArrayType> &)>
291 const unsigned int dof_no = 0,
292 const unsigned int quad_no = 0,
293 const unsigned int first_selected_component = 0);
300 template <
typename CLASS,
306 typename VectorizedArrayType,
318 VectorizedArrayType> &) const,
319 const CLASS *owning_class,
320 const unsigned
int dof_no = 0,
321 const unsigned
int quad_no = 0,
322 const unsigned
int first_selected_component = 0);
333 typename VectorizedArrayType,
335 typename VectorType2>
345 VectorType &diagonal_global,
346 std::vector<VectorType2 *> &diagonal_global_components);
354 typename VectorizedArrayType,
385 typename VectorizedArrayType,
397 VectorizedArrayType> &)>
404 VectorizedArrayType> &,
410 VectorizedArrayType> &)>
417 VectorizedArrayType> &)>
419 const unsigned int dof_no = 0,
420 const unsigned int quad_no = 0,
421 const unsigned int first_selected_component = 0);
428 template <
typename CLASS,
434 typename VectorizedArrayType,
446 VectorizedArrayType> &) const,
452 VectorizedArrayType> &,
458 VectorizedArrayType> &)
465 VectorizedArrayType> &)
467 const CLASS *owning_class,
468 const unsigned
int dof_no = 0,
469 const unsigned
int quad_no = 0,
470 const unsigned
int first_selected_component = 0);
522 std::vector<unsigned int> valid_fe_indices;
524 const auto &fe_collection =
526 .get_fe_collection();
528 for (
unsigned int i = 0; i < fe_collection.size(); ++i)
529 if (fe_collection[i].n_dofs_per_cell() > 0)
530 valid_fe_indices.push_back(i);
544 template <
typename VectorTypeOut,
typename VectorTypeIn>
549 const VectorTypeIn &,
550 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
552 const VectorTypeIn &src,
553 const bool zero_dst_vector =
false)
const
555 const auto ebd_cell_operation = [&](
const auto &
matrix_free,
559 const auto category =
matrix_free.get_cell_range_category(range);
567 matrix_free->template cell_loop<VectorTypeOut, VectorTypeIn>(
568 ebd_cell_operation, dst, src, zero_dst_vector);
578 template <
typename VectorTypeOut,
typename VectorTypeIn>
583 const VectorTypeIn &,
584 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
588 const VectorTypeIn &,
589 const std::pair<unsigned int, unsigned int> &)> &face_operation,
593 const VectorTypeIn &,
594 const std::pair<unsigned int, unsigned int> &,
595 const bool)> &boundary_operation,
597 const VectorTypeIn &src,
598 const bool zero_dst_vector =
false)
const
600 const auto ebd_cell_operation = [&](
const auto &
matrix_free,
604 const auto category =
matrix_free.get_cell_range_category(range);
612 const auto ebd_internal_or_boundary_face_operation =
617 const auto category =
matrix_free.get_face_range_category(range);
619 const unsigned int type =
633 matrix_free->template loop<VectorTypeOut, VectorTypeIn>(
635 ebd_internal_or_boundary_face_operation,
636 ebd_internal_or_boundary_face_operation,
659 template <
int dim,
typename AdditionalData>
662 AdditionalData &additional_data)
665 const unsigned int level = additional_data.mg_level;
670 additional_data.cell_vectorization_category.assign(
673 additional_data.cell_vectorization_category.assign(tria.
n_active_cells(),
680 std::sort(bids.begin(), bids.end());
683 unsigned int n_bids = bids.size() + 1;
685 for (
unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
686 i++, offset = offset * n_bids)
690 const auto to_category = [&](
const auto &cell) {
691 unsigned int c_num = 0;
692 for (
unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
694 const auto face = cell->face(i);
695 if (face->at_boundary() && !cell->has_periodic_neighbor(i))
697 factors[i] * (1 + std::distance(bids.begin(),
698 std::find(bids.begin(),
700 face->boundary_id())));
709 if (cell->is_locally_owned())
711 .cell_vectorization_category[cell->active_cell_index()] =
719 if (cell->is_locally_owned_on_level())
720 additional_data.cell_vectorization_category[cell->index()] =
726 additional_data.hold_all_faces_to_owned_cells =
true;
727 additional_data.cell_vectorization_categories_strict =
true;
728 additional_data.mapping_update_flags_faces_by_cells =
729 additional_data.mapping_update_flags_inner_faces |
730 additional_data.mapping_update_flags_boundary_faces;
735 template <
typename Number>
738 std::vector<unsigned int> row_lid_to_gid;
739 std::vector<unsigned int> row;
740 std::vector<unsigned int> col;
741 std::vector<Number> val;
743 std::vector<unsigned int> inverse_lookup_rows;
744 std::vector<std::pair<unsigned int, unsigned int>> inverse_lookup_origins;
747 template <
int dim,
typename VectorizedArrayType,
bool is_face>
748 class ComputeDiagonalHelper
751 using FEEvaluationType =
754 using Number =
typename VectorizedArrayType::value_type;
755 static const unsigned int n_lanes = VectorizedArrayType::size();
757 ComputeDiagonalHelper()
759 , matrix_free(nullptr)
760 , dofs_per_component(0)
764 ComputeDiagonalHelper(
const ComputeDiagonalHelper &)
766 , matrix_free(nullptr)
767 , dofs_per_component(0)
773 FEEvaluationType &phi,
775 const unsigned int n_components)
779 if (dofs_per_component !=
780 phi.get_shape_info().dofs_per_component_on_cell)
782 locally_relevant_constraints_hn_map.
clear();
784 phi.get_shape_info().dofs_per_component_on_cell;
786 this->n_components = n_components;
787 this->dofs_per_cell = n_components * dofs_per_component;
789 this->matrix_free = &matrix_free;
793 reinit(
const unsigned int cell)
796 const auto &matrix_free = *this->matrix_free;
804 const unsigned int first_selected_component =
805 n_fe_components == 1 ? 0 : phi->get_first_selected_component();
807 this->n_lanes_filled =
816 std::vector<std::tuple<unsigned int, unsigned int, Number>>
817 locally_relevant_constraints, locally_relevant_constraints_tmp;
818 locally_relevant_constraints.reserve(dofs_per_cell);
819 std::vector<unsigned int> constraint_position;
820 std::vector<unsigned char> is_constrained_hn;
822 const std::array<unsigned int, n_lanes> &cells =
823 this->phi->get_cell_ids();
825 for (
unsigned int v = 0; v < n_lanes_filled; ++v)
830 const unsigned int *dof_indices;
831 unsigned int index_indicators, next_index_indicators;
833 const unsigned int start =
834 cells[v] * n_fe_components + first_selected_component;
836 dof_info.dof_indices.data() + dof_info.row_starts[start].first;
837 index_indicators = dof_info.row_starts[start].second;
838 next_index_indicators = dof_info.row_starts[start + 1].second;
842 locally_relevant_constraints.clear();
844 if (n_components == 1 || n_fe_components == 1)
846 unsigned int ind_local = 0;
847 for (; index_indicators != next_index_indicators;
848 ++index_indicators, ++ind_local)
850 const std::pair<unsigned short, unsigned short> indicator =
851 dof_info.constraint_indicator[index_indicators];
853 for (
unsigned int j = 0; j < indicator.first;
855 locally_relevant_constraints.emplace_back(ind_local,
859 dof_indices += indicator.first;
861 const Number *data_val =
863 const Number *end_pool =
866 for (; data_val != end_pool; ++data_val, ++dof_indices)
867 locally_relevant_constraints.emplace_back(ind_local,
874 for (; ind_local < dofs_per_component;
875 ++dof_indices, ++ind_local)
876 locally_relevant_constraints.emplace_back(ind_local,
887 for (
unsigned int comp = 0; comp < n_components; ++comp)
889 unsigned int ind_local = 0;
893 for (; index_indicators != next_index_indicators;
894 ++index_indicators, ++ind_local)
896 const std::pair<unsigned short, unsigned short>
898 dof_info.constraint_indicator[index_indicators];
901 for (
unsigned int j = 0; j < indicator.first;
903 locally_relevant_constraints.emplace_back(
904 comp * dofs_per_component + ind_local,
907 dof_indices += indicator.first;
909 const Number *data_val =
911 const Number *end_pool =
914 for (; data_val != end_pool; ++data_val, ++dof_indices)
915 locally_relevant_constraints.emplace_back(
916 comp * dofs_per_component + ind_local,
924 for (; ind_local < dofs_per_component;
925 ++dof_indices, ++ind_local)
926 locally_relevant_constraints.emplace_back(
927 comp * dofs_per_component + ind_local,
931 if (comp + 1 < n_components)
932 next_index_indicators =
933 dof_info.row_starts[start + comp + 2].second;
940 for (
unsigned int i = 1; i < locally_relevant_constraints.size();
942 Assert(std::get<0>(locally_relevant_constraints[i]) >=
943 std::get<0>(locally_relevant_constraints[i - 1]),
947 if (dof_info.hanging_node_constraint_masks.size() > 0 &&
948 dof_info.hanging_node_constraint_masks_comp.size() > 0 &&
949 dof_info.hanging_node_constraint_masks_comp
950 [phi->get_active_fe_index()][first_selected_component])
953 dof_info.hanging_node_constraint_masks[cells[v]];
956 if (mask != ::internal::MatrixFreeFunctions::
961 if (locally_relevant_constraints_hn_map.find(mask) ==
962 locally_relevant_constraints_hn_map.end())
963 fill_constraint_type_into_map(mask);
965 const auto &locally_relevant_constraints_hn =
966 locally_relevant_constraints_hn_map[
mask];
968 locally_relevant_constraints_tmp.clear();
969 if (locally_relevant_constraints_tmp.capacity() <
970 locally_relevant_constraints.size())
971 locally_relevant_constraints_tmp.reserve(
972 locally_relevant_constraints.size() +
973 locally_relevant_constraints_hn.size());
978 constraint_position.assign(dofs_per_cell,
980 for (
auto &a : locally_relevant_constraints)
981 if (constraint_position[
std::get<0>(a)] ==
983 constraint_position[
std::get<0>(a)] =
984 std::distance(locally_relevant_constraints.data(),
986 is_constrained_hn.assign(dofs_per_cell,
false);
987 for (
auto &hn : locally_relevant_constraints_hn)
988 is_constrained_hn[
std::get<0>(hn)] = 1;
991 for (
const auto &a : locally_relevant_constraints)
992 if (is_constrained_hn[
std::get<0>(a)] == 0)
993 locally_relevant_constraints_tmp.push_back(a);
997 for (
const auto &hn : locally_relevant_constraints_hn)
998 if (constraint_position[
std::get<1>(hn)] !=
1002 locally_relevant_constraints.size());
1003 auto other = locally_relevant_constraints.begin() +
1004 constraint_position[std::get<1>(hn)];
1007 for (; other != locally_relevant_constraints.end() &&
1008 std::get<0>(*other) == std::get<1>(hn);
1010 locally_relevant_constraints_tmp.emplace_back(
1012 std::get<1>(*other),
1013 std::get<2>(hn) * std::get<2>(*other));
1016 std::swap(locally_relevant_constraints,
1017 locally_relevant_constraints_tmp);
1022 std::sort(locally_relevant_constraints.begin(),
1023 locally_relevant_constraints.end(),
1024 [](
const auto &a,
const auto &b) {
1025 if (std::get<1>(a) < std::get<1>(b))
1027 return (std::get<1>(a) == std::get<1>(b)) &&
1028 (std::get<0>(a) < std::get<0>(b));
1032 auto &c_pool = c_pools[v];
1034 c_pool.row_lid_to_gid.clear();
1036 c_pool.row.push_back(0);
1040 if (locally_relevant_constraints.size() > 0)
1041 c_pool.row_lid_to_gid.emplace_back(
1042 std::get<1>(locally_relevant_constraints.front()));
1043 for (
const auto &j : locally_relevant_constraints)
1045 if (c_pool.row_lid_to_gid.back() != std::get<1>(j))
1047 c_pool.row_lid_to_gid.push_back(std::get<1>(j));
1048 c_pool.row.push_back(c_pool.val.size());
1051 c_pool.col.emplace_back(std::get<0>(j));
1052 c_pool.val.emplace_back(std::get<2>(j));
1055 if (c_pool.val.size() > 0)
1056 c_pool.row.push_back(c_pool.val.size());
1058 c_pool.inverse_lookup_rows.clear();
1059 c_pool.inverse_lookup_rows.resize(1 + dofs_per_cell);
1060 for (
const unsigned int i : c_pool.col)
1061 c_pool.inverse_lookup_rows[1 + i]++;
1063 std::partial_sum(c_pool.inverse_lookup_rows.begin(),
1064 c_pool.inverse_lookup_rows.end(),
1065 c_pool.inverse_lookup_rows.begin());
1069 c_pool.inverse_lookup_origins.resize(c_pool.col.size());
1070 std::vector<unsigned int> inverse_lookup_count(dofs_per_cell);
1071 for (
unsigned int row = 0; row < c_pool.row.size() - 1; ++row)
1072 for (
unsigned int col = c_pool.row[row];
1073 col < c_pool.row[row + 1];
1076 const unsigned int index = c_pool.col[col];
1077 c_pool.inverse_lookup_origins
1078 [c_pool.inverse_lookup_rows[
index] +
1079 inverse_lookup_count[
index]] = std::make_pair(row, col);
1080 ++inverse_lookup_count[
index];
1109 for (
unsigned int v = 0; v < n_lanes_filled; ++v)
1110 diagonals_local_constrained[v].assign(
1111 c_pools[v].row_lid_to_gid.size() *
1112 (n_fe_components == 1 ? n_components : 1),
1116 bool use_fast_path =
true;
1118 for (
unsigned int v = 0; v < n_lanes_filled; ++v)
1120 auto &c_pool = c_pools[v];
1122 for (
unsigned int i = 0; i < c_pool.row.size() - 1; ++i)
1124 if ((c_pool.row[i + 1] - c_pool.row[i]) > 1)
1126 use_fast_path =
false;
1129 else if (((c_pool.row[i + 1] - c_pool.row[i]) == 1) &&
1130 (c_pool.val[c_pool.row[i]] != 1.0))
1132 use_fast_path =
false;
1137 if (use_fast_path ==
false)
1141 this->has_simple_constraints_ = use_fast_path;
1145 fill_constraint_type_into_map(
1146 const ::internal::MatrixFreeFunctions::compressed_constraint_kind
1149 auto &constraints_hn = locally_relevant_constraints_hn_map[
mask];
1154 const unsigned int degree =
1155 phi->get_shape_info().data.front().fe_degree;
1160 values_dofs.resize(dofs_per_component);
1163 VectorizedArrayType::size()>
1165 constraint_mask.fill(::internal::MatrixFreeFunctions::
1167 constraint_mask[0] =
mask;
1169 for (
unsigned int i = 0; i < dofs_per_component; ++i)
1171 for (
unsigned int j = 0; j < dofs_per_component; ++j)
1172 values_dofs[j] = VectorizedArrayType();
1173 values_dofs[i] = Number(1);
1178 VectorizedArrayType>::apply(1,
1180 phi->get_shape_info(),
1183 values_dofs.data());
1185 const Number tolerance =
1187 std::numeric_limits<Number>::epsilon() * 16);
1188 for (
unsigned int j = 0; j < dofs_per_component; ++j)
1189 if (
std::abs(values_dofs[j][0]) > tolerance &&
1191 std::abs(values_dofs[j][0] - Number(1)) > tolerance))
1192 constraints_hn.emplace_back(j, i, values_dofs[j][0]);
1196 const unsigned int n_hn_constraints = constraints_hn.size();
1197 constraints_hn.resize(n_hn_constraints * n_components);
1199 for (
unsigned int c = 1; c < n_components; ++c)
1200 for (
unsigned int i = 0; i < n_hn_constraints; ++i)
1201 constraints_hn[c * n_hn_constraints + i] = std::make_tuple(
1202 std::get<0>(constraints_hn[i]) + c * dofs_per_component,
1203 std::get<1>(constraints_hn[i]) + c * dofs_per_component,
1204 std::get<2>(constraints_hn[i]));
1208 prepare_basis_vector(
const unsigned int i)
1215 VectorizedArrayType *dof_values = phi->begin_dof_values();
1216 for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1217 dof_values[j] = VectorizedArrayType();
1218 dof_values[i] = Number(1);
1224 VectorizedArrayType *dof_values = phi->begin_dof_values();
1225 for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1226 dof_values[j] = VectorizedArrayType();
1235 const unsigned int n_fe_components =
1236 phi->get_dof_info().start_components.back();
1237 const unsigned int comp =
1238 n_fe_components == 1 ? i / dofs_per_component : 0;
1239 const unsigned int i_comp =
1240 n_fe_components == 1 ? (i % dofs_per_component) : i;
1244 for (
unsigned int v = 0; v < n_lanes_filled; ++v)
1246 const auto &c_pool = c_pools[v];
1248 for (
unsigned int jj = c_pool.inverse_lookup_rows[i_comp];
1249 jj < c_pool.inverse_lookup_rows[i_comp + 1];
1252 const unsigned int j = c_pool.inverse_lookup_origins[jj].first;
1255 for (
unsigned int k = c_pool.row[j]; k < c_pool.row[j + 1]; ++k)
1256 temp += c_pool.val[k] *
1257 phi->begin_dof_values()[comp * dofs_per_component +
1261 diagonals_local_constrained
1262 [v][j + comp * c_pools[v].row_lid_to_gid.size()] +=
1263 temp * c_pool.val[c_pool.inverse_lookup_origins[jj].second];
1268 template <
typename VectorType>
1270 distribute_local_to_global(std::vector<VectorType *> &diagonal_global)
1273 const unsigned int n_fe_components =
1274 phi->get_dof_info().start_components.back();
1276 if (n_fe_components == 1)
1279 for (
unsigned int v = 0; v < n_lanes_filled; ++v)
1283 for (
unsigned int j = 0; j < c_pools[v].row.size() - 1; ++j)
1284 for (
unsigned int comp = 0;
1285 comp < (n_fe_components == 1 ?
1286 static_cast<unsigned int>(n_components) :
1290 *diagonal_global[n_fe_components == 1 ? comp : 0],
1291 c_pools[v].row_lid_to_gid[j],
1292 diagonals_local_constrained
1293 [v][j + comp * c_pools[v].row_lid_to_gid.size()]);
1297 has_simple_constraints()
const
1299 return has_simple_constraints_;
1303 FEEvaluationType *phi;
1306 unsigned int dofs_per_component;
1307 unsigned int dofs_per_cell;
1308 unsigned int n_components;
1312 unsigned int n_lanes_filled;
1314 std::array<internal::LocalCSR<Number>, n_lanes> c_pools;
1319 std::array<std::vector<Number>, n_lanes> diagonals_local_constrained;
1323 std::vector<std::tuple<unsigned int, unsigned int, Number>>>
1324 locally_relevant_constraints_hn_map;
1329 bool has_simple_constraints_;
1332 template <
bool is_face,
1335 typename VectorizedArrayType>
1339 const std::pair<unsigned int, unsigned int> &range,
1340 const unsigned int dof_no,
1341 const unsigned int quad_no,
1342 const unsigned int first_selected_component,
1343 const unsigned int fe_degree,
1344 const unsigned int n_q_points_1d,
1345 const bool is_interior_face =
true)
1347 const unsigned int static_n_q_points =
1351 unsigned int active_fe_index = 0;
1354 else if (is_interior_face)
1359 const auto init_data = ::internal::
1360 extract_initialization_data<is_face, dim, Number, VectorizedArrayType>(
1363 first_selected_component,
1371 return init_data.shape_info->dofs_per_component_on_cell == 0;
1377 template <
int dim,
int fe_degree,
typename Number,
typename QuadOperation>
1381 CellAction(
const QuadOperation &quad_operation,
1384 : m_quad_operation(quad_operation)
1385 , m_evaluation_flags(evaluation_flags)
1386 , m_integration_flags(integration_flags)
1389 KOKKOS_FUNCTION
void
1390 operator()(
const unsigned int cell,
1397 gpu_data, shared_data);
1398 m_quad_operation.set_matrix_free_data(*gpu_data);
1399 m_quad_operation.set_cell(cell);
1400 constexpr int dofs_per_cell =
decltype(fe_eval)::tensor_dofs_per_cell;
1401 Number
diagonal[dofs_per_cell] = {};
1402 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1404 Kokkos::parallel_for(
1405 Kokkos::TeamThreadRange(shared_data->
team_member, dofs_per_cell),
1406 [&](
int j) { fe_eval.submit_dof_value(i == j ? 1 : 0, j); });
1408 Portable::internal::
1409 resolve_hanging_nodes<dim, fe_degree, false, Number>(
1413 Kokkos::subview(shared_data->
values, Kokkos::ALL, 0));
1415 fe_eval.evaluate(m_evaluation_flags);
1416 fe_eval.apply_for_each_quad_point(m_quad_operation);
1417 fe_eval.integrate(m_integration_flags);
1419 Portable::internal::
1420 resolve_hanging_nodes<dim, fe_degree, true, Number>(
1424 Kokkos::subview(shared_data->
values, Kokkos::ALL, 0));
1426 Kokkos::single(Kokkos::PerTeam(shared_data->
team_member),
1427 [&] { diagonal[i] = fe_eval.get_dof_value(i); });
1430 Kokkos::single(Kokkos::PerTeam(shared_data->
team_member), [&] {
1431 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1432 fe_eval.submit_dof_value(diagonal[i], i);
1439 Kokkos::parallel_for(Kokkos::TeamThreadRange(shared_data->
team_member,
1442 dst[gpu_data->local_to_global(cell, i)] +=
1443 shared_data->values(i, 0);
1448 Kokkos::parallel_for(
1449 Kokkos::TeamThreadRange(shared_data->
team_member, dofs_per_cell),
1451 Kokkos::atomic_add(&dst[gpu_data->local_to_global(cell, i)],
1452 shared_data->values(i, 0));
1457 static constexpr unsigned int n_local_dofs = QuadOperation::n_local_dofs;
1460 mutable QuadOperation m_quad_operation;
1472 typename QuadOperation>
1477 const QuadOperation &quad_operation,
1480 const unsigned int dof_no,
1481 const unsigned int quad_no,
1482 const unsigned int first_selected_component,
1483 const unsigned int first_vector_component)
1493 CellAction<dim, fe_degree, Number, QuadOperation>
cell_action(
1494 quad_operation, evaluation_flags, integration_flags);
1496 matrix_free.
cell_loop(cell_action, dummy, diagonal_global);
1506 typename VectorizedArrayType,
1511 VectorType &diagonal_global,
1517 VectorizedArrayType> &)>
1519 const unsigned int dof_no,
1520 const unsigned int quad_no,
1521 const unsigned int first_selected_component,
1522 const unsigned int first_vector_component)
1529 VectorizedArrayType,
1537 first_selected_component,
1538 first_vector_component);
1541 template <
typename CLASS,
1547 typename VectorizedArrayType,
1552 VectorType &diagonal_global,
1558 VectorizedArrayType> &) const,
1559 const CLASS *owning_class,
1560 const unsigned
int dof_no,
1561 const unsigned
int quad_no,
1562 const unsigned
int first_selected_component,
1563 const unsigned
int first_vector_component)
1570 VectorizedArrayType,
1574 [&](
auto &phi) { (owning_class->*cell_operation)(phi); },
1577 first_selected_component,
1578 first_vector_component);
1586 typename VectorizedArrayType,
1591 VectorType &diagonal_global,
1597 VectorizedArrayType> &)>
1604 VectorizedArrayType> &,
1610 VectorizedArrayType> &)>
1617 VectorizedArrayType> &)>
1618 &boundary_operation,
1619 const unsigned int dof_no,
1620 const unsigned int quad_no,
1621 const unsigned int first_selected_component,
1622 const unsigned int first_vector_component)
1624 std::vector<typename ::internal::BlockVectorSelector<
1627 diagonal_global_components(n_components);
1629 for (
unsigned int d = 0;
d < n_components; ++
d)
1630 diagonal_global_components[d] = ::internal::
1632 get_vector_component(diagonal_global,
d + first_vector_component);
1634 const auto &dof_info = matrix_free.
get_dof_info(dof_no);
1636 if (dof_info.start_components.back() == 1)
1637 for (
unsigned int comp = 0; comp < n_components; ++comp)
1639 Assert(diagonal_global_components[comp] !=
nullptr,
1640 ExcMessage(
"The finite element underlying this FEEvaluation "
1641 "object is scalar, but you requested " +
1642 std::to_string(n_components) +
1643 " components via the template argument in "
1644 "FEEvaluation. In that case, you must pass an "
1645 "std::vector<VectorType> or a BlockVector to " +
1646 "read_dof_values and distribute_local_to_global."));
1648 *diagonal_global_components[comp], matrix_free, dof_info);
1653 *diagonal_global_components[0], matrix_free, dof_info);
1661 VectorizedArrayType>;
1668 VectorizedArrayType>;
1670 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
1673 data_cell.dof_numbers = {dof_no};
1674 data_cell.quad_numbers = {quad_no};
1675 data_cell.n_components = {n_components};
1676 data_cell.first_selected_components = {first_selected_component};
1677 data_cell.batch_type = {0};
1679 data_cell.op_create =
1680 [&](
const std::pair<unsigned int, unsigned int> &range) {
1682 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, false>>>
1685 if (!internal::is_fe_nothing<false>(matrix_free,
1689 first_selected_component,
1692 phi.emplace_back(std::make_unique<FEEvalType>(
1693 matrix_free, range, dof_no, quad_no, first_selected_component));
1698 data_cell.op_reinit = [](
auto &phi,
const unsigned batch) {
1699 static_cast<FEEvalType &
>(*phi[0]).
reinit(batch);
1703 data_cell.op_compute = [&](
auto &phi) {
1704 cell_operation(
static_cast<FEEvalType &
>(*phi[0]));
1707 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1710 data_face.dof_numbers = {dof_no, dof_no};
1711 data_face.quad_numbers = {quad_no, quad_no};
1712 data_face.n_components = {n_components, n_components};
1713 data_face.first_selected_components = {first_selected_component,
1714 first_selected_component};
1715 data_face.batch_type = {1, 2};
1717 data_face.op_create =
1718 [&](
const std::pair<unsigned int, unsigned int> &range) {
1720 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
1723 if (!internal::is_fe_nothing<true>(matrix_free,
1727 first_selected_component,
1731 !internal::is_fe_nothing<true>(matrix_free,
1735 first_selected_component,
1741 std::make_unique<FEFaceEvalType>(matrix_free,
1746 first_selected_component));
1748 std::make_unique<FEFaceEvalType>(matrix_free,
1753 first_selected_component));
1759 data_face.op_reinit = [](
auto &phi,
const unsigned batch) {
1760 static_cast<FEFaceEvalType &
>(*phi[0]).
reinit(batch);
1761 static_cast<FEFaceEvalType &
>(*phi[1]).
reinit(batch);
1765 data_face.op_compute = [&](
auto &phi) {
1766 face_operation(
static_cast<FEFaceEvalType &
>(*phi[0]),
1767 static_cast<FEFaceEvalType &
>(*phi[1]));
1770 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1773 data_boundary.dof_numbers = {dof_no};
1774 data_boundary.quad_numbers = {quad_no};
1775 data_boundary.n_components = {n_components};
1776 data_boundary.first_selected_components = {first_selected_component};
1777 data_boundary.batch_type = {1};
1780 .op_create = [&](
const std::pair<unsigned int, unsigned int> &range) {
1782 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
1785 if (!internal::is_fe_nothing<true>(matrix_free,
1789 first_selected_component,
1793 phi.emplace_back(std::make_unique<FEFaceEvalType>(
1794 matrix_free, range,
true, dof_no, quad_no, first_selected_component));
1799 data_boundary.op_reinit = [](
auto &phi,
const unsigned batch) {
1800 static_cast<FEFaceEvalType &
>(*phi[0]).
reinit(batch);
1803 if (boundary_operation)
1804 data_boundary.op_compute = [&](
auto &phi) {
1805 boundary_operation(
static_cast<FEFaceEvalType &
>(*phi[0]));
1813 diagonal_global_components);
1820 typename VectorizedArrayType,
1822 typename VectorType2>
1826 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
1828 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1830 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1832 VectorType &diagonal_global,
1833 std::vector<VectorType2 *> &diagonal_global_components)
1840 internal::ComputeDiagonalHelper<dim, VectorizedArrayType, false>;
1843 internal::ComputeDiagonalHelper<dim, VectorizedArrayType, true>;
1847 scratch_data_internal;
1850 const auto batch_operation =
1853 const std::pair<unsigned int, unsigned int> &range) {
1854 if (!data.op_compute)
1857 auto phi = data.op_create(range);
1859 const unsigned int n_blocks = phi.size();
1861 auto &helpers = scratch_data.
get();
1862 helpers.resize(n_blocks);
1865 helpers[b].initialize(*phi[b], matrix_free, data.n_components[b]);
1867 for (
unsigned int batch = range.first; batch < range.second; ++batch)
1869 data.op_reinit(phi, batch);
1872 helpers[b].
reinit(batch);
1876 Assert(std::all_of(helpers.begin(),
1878 [](
const auto &helper) {
1879 return helper.has_simple_constraints();
1886 for (
unsigned int i = 0;
1887 i < phi[
b]->get_shape_info().dofs_per_component_on_cell *
1888 data.n_components[b];
1891 for (
unsigned int bb = 0; bb <
n_blocks; ++bb)
1893 helpers[bb].prepare_basis_vector(i);
1895 helpers[bb].zero_basis_vector();
1897 data.op_compute(phi);
1898 helpers[
b].submit();
1901 helpers[
b].distribute_local_to_global(
1902 diagonal_global_components);
1907 const auto cell_operation_wrapped =
1908 [&](
const auto &,
auto &,
const auto &,
const auto range) {
1909 batch_operation(data_cell, scratch_data, range);
1912 const auto face_operation_wrapped =
1913 [&](
const auto &,
auto &,
const auto &,
const auto range) {
1914 batch_operation(data_face, scratch_data_internal, range);
1917 const auto boundary_operation_wrapped =
1918 [&](
const auto &,
auto &,
const auto &,
const auto range) {
1919 batch_operation(data_boundary, scratch_data_bc, range);
1922 if (data_face.op_compute || data_boundary.op_compute)
1923 matrix_free.template loop<VectorType, int>(cell_operation_wrapped,
1924 face_operation_wrapped,
1925 boundary_operation_wrapped,
1930 matrix_free.template cell_loop<VectorType, int>(cell_operation_wrapped,
1937 template <
typename CLASS,
1943 typename VectorizedArrayType,
1948 VectorType &diagonal_global,
1954 VectorizedArrayType> &) const,
1960 VectorizedArrayType> &,
1966 VectorizedArrayType> &)
1973 VectorizedArrayType> &)
1975 const CLASS *owning_class,
1976 const unsigned
int dof_no,
1977 const unsigned
int quad_no,
1978 const unsigned
int first_selected_component,
1979 const unsigned
int first_vector_component)
1986 VectorizedArrayType,
1990 [&](
auto &phi) { (owning_class->*cell_operation)(phi); },
1991 [&](
auto &phi_m,
auto &phi_p) {
1992 (owning_class->*face_operation)(phi_m, phi_p);
1994 [&](
auto &phi) { (owning_class->*boundary_operation)(phi); },
1997 first_selected_component,
1998 first_vector_component);
2010 std::enable_if_t<std::is_same_v<
2011 std::remove_const_t<
2012 std::remove_reference_t<typename MatrixType::value_type>>,
2013 std::remove_const_t<std::remove_reference_t<Number>>>> * =
nullptr>
2015 create_new_affine_constraints_if_needed(
2031 std::enable_if_t<!std::is_same_v<
2032 std::remove_const_t<
2033 std::remove_reference_t<typename MatrixType::value_type>>,
2034 std::remove_const_t<std::remove_reference_t<Number>>>> * =
nullptr>
2036 create_new_affine_constraints_if_needed(
2043 std::make_unique<AffineConstraints<typename MatrixType::value_type>>();
2044 new_constraints->
copy_from(constraints);
2046 return *new_constraints;
2055 typename VectorizedArrayType,
2067 VectorizedArrayType> &)>
2069 const unsigned int dof_no,
2070 const unsigned int quad_no,
2071 const unsigned int first_selected_component)
2078 VectorizedArrayType,
2087 first_selected_component);
2090 template <
typename CLASS,
2096 typename VectorizedArrayType,
2108 VectorizedArrayType> &) const,
2109 const CLASS *owning_class,
2110 const unsigned
int dof_no,
2111 const unsigned
int quad_no,
2112 const unsigned
int first_selected_component)
2119 VectorizedArrayType,
2124 [&](
auto &phi) { (owning_class->*cell_operation)(phi); },
2127 first_selected_component);
2135 typename VectorizedArrayType,
2147 VectorizedArrayType> &)>
2154 VectorizedArrayType> &,
2160 VectorizedArrayType> &)>
2167 VectorizedArrayType> &)>
2168 &boundary_operation,
2169 const unsigned int dof_no,
2170 const unsigned int quad_no,
2171 const unsigned int first_selected_component)
2178 VectorizedArrayType>;
2185 VectorizedArrayType>;
2187 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
2190 data_cell.dof_numbers = {dof_no};
2191 data_cell.quad_numbers = {quad_no};
2192 data_cell.n_components = {n_components};
2193 data_cell.first_selected_components = {first_selected_component};
2194 data_cell.batch_type = {0};
2196 data_cell.op_create =
2197 [&](
const std::pair<unsigned int, unsigned int> &range) {
2199 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, false>>>
2202 if (!internal::is_fe_nothing<false>(matrix_free,
2206 first_selected_component,
2209 phi.emplace_back(std::make_unique<FEEvalType>(
2210 matrix_free, range, dof_no, quad_no, first_selected_component));
2215 data_cell.op_reinit = [](
auto &phi,
const unsigned batch) {
2216 static_cast<FEEvalType &
>(*phi[0]).
reinit(batch);
2220 data_cell.op_compute = [&](
auto &phi) {
2221 cell_operation(
static_cast<FEEvalType &
>(*phi[0]));
2224 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2227 data_face.dof_numbers = {dof_no, dof_no};
2228 data_face.quad_numbers = {quad_no, quad_no};
2229 data_face.n_components = {n_components, n_components};
2230 data_face.first_selected_components = {first_selected_component,
2231 first_selected_component};
2232 data_face.batch_type = {1, 2};
2234 data_face.op_create =
2235 [&](
const std::pair<unsigned int, unsigned int> &range) {
2237 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
2240 if (!internal::is_fe_nothing<true>(matrix_free,
2244 first_selected_component,
2248 !internal::is_fe_nothing<true>(matrix_free,
2252 first_selected_component,
2258 std::make_unique<FEFaceEvalType>(matrix_free,
2263 first_selected_component));
2265 std::make_unique<FEFaceEvalType>(matrix_free,
2270 first_selected_component));
2276 data_face.op_reinit = [](
auto &phi,
const unsigned batch) {
2277 static_cast<FEFaceEvalType &
>(*phi[0]).
reinit(batch);
2278 static_cast<FEFaceEvalType &
>(*phi[1]).
reinit(batch);
2282 data_face.op_compute = [&](
auto &phi) {
2283 face_operation(
static_cast<FEFaceEvalType &
>(*phi[0]),
2284 static_cast<FEFaceEvalType &
>(*phi[1]));
2287 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2290 data_boundary.dof_numbers = {dof_no};
2291 data_boundary.quad_numbers = {quad_no};
2292 data_boundary.n_components = {n_components};
2293 data_boundary.first_selected_components = {first_selected_component};
2294 data_boundary.batch_type = {1};
2297 .op_create = [&](
const std::pair<unsigned int, unsigned int> &range) {
2299 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
2302 if (!internal::is_fe_nothing<true>(matrix_free,
2306 first_selected_component,
2310 phi.emplace_back(std::make_unique<FEFaceEvalType>(
2311 matrix_free, range,
true, dof_no, quad_no, first_selected_component));
2316 data_boundary.op_reinit = [](
auto &phi,
const unsigned batch) {
2317 static_cast<FEFaceEvalType &
>(*phi[0]).
reinit(batch);
2320 if (boundary_operation)
2321 data_boundary.op_compute = [&](
auto &phi) {
2322 boundary_operation(
static_cast<FEFaceEvalType &
>(*phi[0]));
2326 matrix_free, constraints_in, data_cell, data_face, data_boundary, matrix);
2333 typename VectorizedArrayType,
2339 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
2341 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2343 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2347 std::unique_ptr<AffineConstraints<typename MatrixType::value_type>>
2348 constraints_for_matrix;
2350 internal::create_new_affine_constraints_if_needed(
2351 matrix, constraints_in, constraints_for_matrix);
2353 const auto batch_operation =
2354 [&matrix_free, &constraints, &
matrix](
2355 auto &data,
const std::pair<unsigned int, unsigned int> &range) {
2356 if (!data.op_compute)
2359 auto phi = data.op_create(range);
2361 const unsigned int n_blocks = phi.size();
2370 n_blocks, VectorizedArrayType::size());
2373 std::array<FullMatrix<typename MatrixType::value_type>,
2374 VectorizedArrayType::size()>>
2375 matrices(n_blocks, n_blocks);
2380 .get_fe(phi[b]->get_active_fe_index());
2382 const auto component_base =
2385 const auto component_in_base =
2386 data.first_selected_components[
b] -
2391 data.dof_numbers[b],
2392 data.quad_numbers[b],
2394 phi[b]->get_active_fe_index(),
2395 phi[b]->get_active_quadrature_index());
2398 shape_info.dofs_per_component_on_cell * data.n_components[
b];
2400 dof_indices[
b].resize(fe.n_dofs_per_cell());
2402 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2403 dof_indices_mf[b][v].resize(dofs_per_cell[b]);
2405 lexicographic_numbering[
b].insert(
2406 lexicographic_numbering[b].
begin(),
2407 shape_info.lexicographic_numbering.begin() +
2408 component_in_base * shape_info.dofs_per_component_on_cell,
2409 shape_info.lexicographic_numbering.begin() +
2410 (component_in_base + data.n_components[
b]) *
2411 shape_info.dofs_per_component_on_cell);
2414 for (
unsigned int bj = 0; bj <
n_blocks; ++bj)
2415 for (
unsigned int bi = 0; bi <
n_blocks; ++bi)
2416 std::fill_n(matrices[bi][bj].
begin(),
2417 VectorizedArrayType::size(),
2419 dofs_per_cell[bi], dofs_per_cell[bj]));
2421 for (
auto batch = range.first; batch < range.second; ++batch)
2423 data.op_reinit(phi, batch);
2425 const unsigned int n_filled_lanes =
2430 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
2434 (data.batch_type[
b] == 0) ?
2435 (batch * VectorizedArrayType::size() + v) :
2436 ((data.batch_type[
b] == 1) ?
2437 matrix_free.get_face_info(batch).cells_interior[v] :
2438 matrix_free.get_face_info(batch).cells_exterior[v]);
2443 data.dof_numbers[b]);
2447 cell_iterator->get_mg_dof_indices(dof_indices[b]);
2449 cell_iterator->get_dof_indices(dof_indices[b]);
2451 for (
unsigned int j = 0; j < dofs_per_cell[
b]; ++j)
2452 dof_indices_mf[b][v][j] =
2453 dof_indices[b][lexicographic_numbering[b][j]];
2456 for (
unsigned int bj = 0; bj <
n_blocks; ++bj)
2458 for (
unsigned int j = 0; j < dofs_per_cell[bj]; ++j)
2460 for (
unsigned int bi = 0; bi <
n_blocks; ++bi)
2461 for (
unsigned int i = 0; i < dofs_per_cell[bi]; ++i)
2462 phi[bi]->begin_dof_values()[i] =
2463 (bj == bi) ?
static_cast<Number
>(i == j) : 0.0;
2465 data.op_compute(phi);
2467 for (
unsigned int bi = 0; bi <
n_blocks; ++bi)
2468 for (
unsigned int i = 0; i < dofs_per_cell[bi]; ++i)
2469 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
2470 matrices[bi][bj][v](i, j) =
2471 phi[bi]->begin_dof_values()[i][v];
2474 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
2475 for (
unsigned int bi = 0; bi <
n_blocks; ++bi)
2483 matrices[bi][bi][v], dof_indices_mf[bi][v], matrix);
2486 matrices[bi][bj][v],
2487 dof_indices_mf[bi][v],
2488 dof_indices_mf[bj][v],
2494 const auto cell_operation_wrapped =
2495 [&](
const auto &,
auto &,
const auto &,
const auto range) {
2496 batch_operation(data_cell, range);
2499 const auto face_operation_wrapped =
2500 [&](
const auto &,
auto &,
const auto &,
const auto range) {
2501 batch_operation(data_face, range);
2504 const auto boundary_operation_wrapped =
2505 [&](
const auto &,
auto &,
const auto &,
const auto range) {
2506 batch_operation(data_boundary, range);
2509 if (data_face.op_compute || data_boundary.op_compute)
2511 matrix_free.template loop<MatrixType, MatrixType>(
2512 cell_operation_wrapped,
2513 face_operation_wrapped,
2514 boundary_operation_wrapped,
2519 matrix_free.template cell_loop<MatrixType, MatrixType>(
2520 cell_operation_wrapped, matrix, matrix);
2526 template <
typename CLASS,
2532 typename VectorizedArrayType,
2544 VectorizedArrayType> &) const,
2550 VectorizedArrayType> &,
2556 VectorizedArrayType> &)
2563 VectorizedArrayType> &)
2565 const CLASS *owning_class,
2566 const unsigned
int dof_no,
2567 const unsigned
int quad_no,
2568 const unsigned
int first_selected_component)
2575 VectorizedArrayType,
2580 [&](
auto &phi) { (owning_class->*cell_operation)(phi); },
2581 [&](
auto &phi_m,
auto &phi_p) {
2582 (owning_class->*face_operation)(phi_m, phi_p);
2584 [&](
auto &phi) { (owning_class->*boundary_operation)(phi); },
2587 first_selected_component);