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;
1382 typename QuadOperation>
1386 CellAction(
const QuadOperation &quad_operation,
1389 : m_quad_operation(quad_operation)
1390 , m_evaluation_flags(evaluation_flags)
1391 , m_integration_flags(integration_flags)
1394 KOKKOS_FUNCTION
void
1395 operator()(
const unsigned int cell,
1402 FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number>
1403 fe_eval(gpu_data, shared_data);
1404 m_quad_operation.set_matrix_free_data(*gpu_data);
1405 m_quad_operation.set_cell(cell);
1406 constexpr int dofs_per_cell =
decltype(fe_eval)::tensor_dofs_per_cell;
1407 typename decltype(fe_eval)::value_type
1408 diagonal[dofs_per_cell / n_components] = {};
1409 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1411 const auto c = i % n_components;
1413 Kokkos::parallel_for(
1415 dofs_per_cell / n_components),
1417 typename decltype(fe_eval)::value_type val = {};
1419 if constexpr (n_components == 1)
1421 val = (i == j) ? 1 : 0;
1425 val[c] = (i / n_components == j) ? 1 : 0;
1428 fe_eval.submit_dof_value(val, j);
1433 Portable::internal::
1434 resolve_hanging_nodes<dim, fe_degree, false, Number>(
1438 Kokkos::subview(shared_data->
values, Kokkos::ALL, c));
1440 fe_eval.evaluate(m_evaluation_flags);
1441 fe_eval.apply_for_each_quad_point(m_quad_operation);
1442 fe_eval.integrate(m_integration_flags);
1444 Portable::internal::
1445 resolve_hanging_nodes<dim, fe_degree, true, Number>(
1449 Kokkos::subview(shared_data->
values, Kokkos::ALL, c));
1451 Kokkos::single(Kokkos::PerTeam(shared_data->
team_member), [&] {
1452 if constexpr (n_components == 1)
1453 diagonal[i] = fe_eval.get_dof_value(i);
1455 diagonal[i / n_components][i % n_components] =
1456 fe_eval.get_dof_value(i / n_components)[i % n_components];
1462 Kokkos::single(Kokkos::PerTeam(shared_data->
team_member), [&] {
1463 for (unsigned int i = 0; i < dofs_per_cell / n_components; ++i)
1464 fe_eval.submit_dof_value(diagonal[i], i);
1473 Kokkos::parallel_for(
1474 Kokkos::TeamThreadRange(shared_data->
team_member, dofs_per_cell),
1476 dst[gpu_data->local_to_global(cell, i)] +=
1477 shared_data->values(i % (dofs_per_cell / n_components),
1478 i / (dofs_per_cell / n_components));
1483 Kokkos::parallel_for(
1484 Kokkos::TeamThreadRange(shared_data->
team_member, dofs_per_cell),
1487 &dst[gpu_data->local_to_global(cell, i)],
1488 shared_data->values(i % (dofs_per_cell / n_components),
1489 i / (dofs_per_cell / n_components)));
1494 static constexpr unsigned int n_local_dofs =
1496 static constexpr unsigned int n_q_points =
1500 mutable QuadOperation m_quad_operation;
1512 typename QuadOperation>
1517 const QuadOperation &quad_operation,
1520 const unsigned int dof_no,
1521 const unsigned int quad_no,
1522 const unsigned int first_selected_component,
1523 const unsigned int first_vector_component)
1539 cell_action(quad_operation, evaluation_flags, integration_flags);
1541 matrix_free.
cell_loop(cell_action, dummy, diagonal_global);
1551 typename VectorizedArrayType,
1556 VectorType &diagonal_global,
1562 VectorizedArrayType> &)>
1564 const unsigned int dof_no,
1565 const unsigned int quad_no,
1566 const unsigned int first_selected_component,
1567 const unsigned int first_vector_component)
1574 VectorizedArrayType,
1582 first_selected_component,
1583 first_vector_component);
1586 template <
typename CLASS,
1592 typename VectorizedArrayType,
1597 VectorType &diagonal_global,
1603 VectorizedArrayType> &) const,
1604 const CLASS *owning_class,
1605 const unsigned
int dof_no,
1606 const unsigned
int quad_no,
1607 const unsigned
int first_selected_component,
1608 const unsigned
int first_vector_component)
1615 VectorizedArrayType,
1619 [&](
auto &phi) { (owning_class->*cell_operation)(phi); },
1622 first_selected_component,
1623 first_vector_component);
1631 typename VectorizedArrayType,
1636 VectorType &diagonal_global,
1642 VectorizedArrayType> &)>
1649 VectorizedArrayType> &,
1655 VectorizedArrayType> &)>
1662 VectorizedArrayType> &)>
1663 &boundary_operation,
1664 const unsigned int dof_no,
1665 const unsigned int quad_no,
1666 const unsigned int first_selected_component,
1667 const unsigned int first_vector_component)
1669 std::vector<typename ::internal::BlockVectorSelector<
1672 diagonal_global_components(n_components);
1674 for (
unsigned int d = 0;
d < n_components; ++
d)
1675 diagonal_global_components[d] = ::internal::
1677 get_vector_component(diagonal_global,
d + first_vector_component);
1679 const auto &dof_info = matrix_free.
get_dof_info(dof_no);
1681 if (dof_info.start_components.back() == 1)
1682 for (
unsigned int comp = 0; comp < n_components; ++comp)
1684 Assert(diagonal_global_components[comp] !=
nullptr,
1685 ExcMessage(
"The finite element underlying this FEEvaluation "
1686 "object is scalar, but you requested " +
1687 std::to_string(n_components) +
1688 " components via the template argument in "
1689 "FEEvaluation. In that case, you must pass an "
1690 "std::vector<VectorType> or a BlockVector to " +
1691 "read_dof_values and distribute_local_to_global."));
1693 *diagonal_global_components[comp], matrix_free, dof_info);
1698 *diagonal_global_components[0], matrix_free, dof_info);
1706 VectorizedArrayType>;
1713 VectorizedArrayType>;
1715 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
1718 data_cell.dof_numbers = {dof_no};
1719 data_cell.quad_numbers = {quad_no};
1720 data_cell.n_components = {n_components};
1721 data_cell.first_selected_components = {first_selected_component};
1722 data_cell.batch_type = {0};
1724 data_cell.op_create =
1725 [&](
const std::pair<unsigned int, unsigned int> &range) {
1727 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, false>>>
1730 if (!internal::is_fe_nothing<false>(matrix_free,
1734 first_selected_component,
1737 phi.emplace_back(std::make_unique<FEEvalType>(
1738 matrix_free, range, dof_no, quad_no, first_selected_component));
1743 data_cell.op_reinit = [](
auto &phi,
const unsigned batch) {
1744 if (phi.size() == 1)
1745 static_cast<FEEvalType &
>(*phi[0]).reinit(batch);
1749 data_cell.op_compute = [&](
auto &phi) {
1750 cell_operation(
static_cast<FEEvalType &
>(*phi[0]));
1753 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1756 data_face.dof_numbers = {dof_no, dof_no};
1757 data_face.quad_numbers = {quad_no, quad_no};
1758 data_face.n_components = {n_components, n_components};
1759 data_face.first_selected_components = {first_selected_component,
1760 first_selected_component};
1761 data_face.batch_type = {1, 2};
1763 data_face.op_create =
1764 [&](
const std::pair<unsigned int, unsigned int> &range) {
1766 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
1769 if (!internal::is_fe_nothing<true>(matrix_free,
1773 first_selected_component,
1777 !internal::is_fe_nothing<true>(matrix_free,
1781 first_selected_component,
1787 std::make_unique<FEFaceEvalType>(matrix_free,
1792 first_selected_component));
1794 std::make_unique<FEFaceEvalType>(matrix_free,
1799 first_selected_component));
1805 data_face.op_reinit = [](
auto &phi,
const unsigned batch) {
1806 if (phi.size() == 2)
1808 static_cast<FEFaceEvalType &
>(*phi[0]).
reinit(batch);
1809 static_cast<FEFaceEvalType &
>(*phi[1]).
reinit(batch);
1814 data_face.op_compute = [&](
auto &phi) {
1815 face_operation(
static_cast<FEFaceEvalType &
>(*phi[0]),
1816 static_cast<FEFaceEvalType &
>(*phi[1]));
1819 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1822 data_boundary.dof_numbers = {dof_no};
1823 data_boundary.quad_numbers = {quad_no};
1824 data_boundary.n_components = {n_components};
1825 data_boundary.first_selected_components = {first_selected_component};
1826 data_boundary.batch_type = {1};
1829 .op_create = [&](
const std::pair<unsigned int, unsigned int> &range) {
1831 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
1834 if (!internal::is_fe_nothing<true>(matrix_free,
1838 first_selected_component,
1842 phi.emplace_back(std::make_unique<FEFaceEvalType>(
1843 matrix_free, range,
true, dof_no, quad_no, first_selected_component));
1848 data_boundary.op_reinit = [](
auto &phi,
const unsigned batch) {
1849 if (phi.size() == 1)
1850 static_cast<FEFaceEvalType &
>(*phi[0]).reinit(batch);
1853 if (boundary_operation)
1854 data_boundary.op_compute = [&](
auto &phi) {
1855 boundary_operation(
static_cast<FEFaceEvalType &
>(*phi[0]));
1858 internal::compute_diagonal(matrix_free,
1863 diagonal_global_components);
1870 typename VectorizedArrayType,
1872 typename VectorType2>
1876 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
1878 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1880 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1882 VectorType &diagonal_global,
1883 std::vector<VectorType2 *> &diagonal_global_components)
1890 internal::ComputeDiagonalHelper<dim, VectorizedArrayType, false>;
1893 internal::ComputeDiagonalHelper<dim, VectorizedArrayType, true>;
1897 scratch_data_internal;
1900 const auto batch_operation =
1903 const std::pair<unsigned int, unsigned int> &range) {
1904 if (!
data.op_compute)
1907 auto phi =
data.op_create(range);
1909 const unsigned int n_blocks = phi.size();
1911 auto &helpers = scratch_data.
get();
1912 helpers.resize(n_blocks);
1915 helpers[b].initialize(*phi[b], matrix_free,
data.n_components[b]);
1917 for (
unsigned int batch = range.first; batch < range.second; ++batch)
1919 data.op_reinit(phi, batch);
1922 helpers[b].
reinit(batch);
1926 Assert(std::all_of(helpers.begin(),
1928 [](
const auto &helper) {
1929 return helper.has_simple_constraints();
1936 for (
unsigned int i = 0;
1937 i < phi[
b]->get_shape_info().dofs_per_component_on_cell *
1938 data.n_components[b];
1941 for (
unsigned int bb = 0; bb <
n_blocks; ++bb)
1943 helpers[bb].prepare_basis_vector(i);
1945 helpers[bb].zero_basis_vector();
1947 data.op_compute(phi);
1948 helpers[
b].submit();
1951 helpers[
b].distribute_local_to_global(
1952 diagonal_global_components);
1957 const auto cell_operation_wrapped =
1958 [&](
const auto &,
auto &,
const auto &,
const auto range) {
1959 batch_operation(data_cell, scratch_data, range);
1962 const auto face_operation_wrapped =
1963 [&](
const auto &,
auto &,
const auto &,
const auto range) {
1964 batch_operation(data_face, scratch_data_internal, range);
1967 const auto boundary_operation_wrapped =
1968 [&](
const auto &,
auto &,
const auto &,
const auto range) {
1969 batch_operation(data_boundary, scratch_data_bc, range);
1972 if (data_face.op_compute || data_boundary.op_compute)
1973 matrix_free.template loop<VectorType, int>(cell_operation_wrapped,
1974 face_operation_wrapped,
1975 boundary_operation_wrapped,
1980 matrix_free.template cell_loop<VectorType, int>(cell_operation_wrapped,
1987 template <
typename CLASS,
1993 typename VectorizedArrayType,
1998 VectorType &diagonal_global,
2004 VectorizedArrayType> &) const,
2010 VectorizedArrayType> &,
2016 VectorizedArrayType> &)
2023 VectorizedArrayType> &)
2025 const CLASS *owning_class,
2026 const unsigned
int dof_no,
2027 const unsigned
int quad_no,
2028 const unsigned
int first_selected_component,
2029 const unsigned
int first_vector_component)
2036 VectorizedArrayType,
2040 [&](
auto &phi) { (owning_class->*cell_operation)(phi); },
2041 [&](
auto &phi_m,
auto &phi_p) {
2042 (owning_class->*face_operation)(phi_m, phi_p);
2044 [&](
auto &phi) { (owning_class->*boundary_operation)(phi); },
2047 first_selected_component,
2048 first_vector_component);
2060 std::enable_if_t<std::is_same_v<
2061 std::remove_const_t<
2062 std::remove_reference_t<typename MatrixType::value_type>>,
2063 std::remove_const_t<std::remove_reference_t<Number>>>> * =
nullptr>
2065 create_new_affine_constraints_if_needed(
2081 std::enable_if_t<!std::is_same_v<
2082 std::remove_const_t<
2083 std::remove_reference_t<typename MatrixType::value_type>>,
2084 std::remove_const_t<std::remove_reference_t<Number>>>> * =
nullptr>
2086 create_new_affine_constraints_if_needed(
2093 std::make_unique<AffineConstraints<typename MatrixType::value_type>>();
2094 new_constraints->
copy_from(constraints);
2096 return *new_constraints;
2105 typename VectorizedArrayType,
2117 VectorizedArrayType> &)>
2119 const unsigned int dof_no,
2120 const unsigned int quad_no,
2121 const unsigned int first_selected_component)
2128 VectorizedArrayType,
2137 first_selected_component);
2140 template <
typename CLASS,
2146 typename VectorizedArrayType,
2158 VectorizedArrayType> &) const,
2159 const CLASS *owning_class,
2160 const unsigned
int dof_no,
2161 const unsigned
int quad_no,
2162 const unsigned
int first_selected_component)
2169 VectorizedArrayType,
2174 [&](
auto &phi) { (owning_class->*cell_operation)(phi); },
2177 first_selected_component);
2185 typename VectorizedArrayType,
2197 VectorizedArrayType> &)>
2204 VectorizedArrayType> &,
2210 VectorizedArrayType> &)>
2217 VectorizedArrayType> &)>
2218 &boundary_operation,
2219 const unsigned int dof_no,
2220 const unsigned int quad_no,
2221 const unsigned int first_selected_component)
2228 VectorizedArrayType>;
2235 VectorizedArrayType>;
2237 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
2240 data_cell.dof_numbers = {dof_no};
2241 data_cell.quad_numbers = {quad_no};
2242 data_cell.n_components = {n_components};
2243 data_cell.first_selected_components = {first_selected_component};
2244 data_cell.batch_type = {0};
2246 data_cell.op_create =
2247 [&](
const std::pair<unsigned int, unsigned int> &range) {
2249 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, false>>>
2252 if (!internal::is_fe_nothing<false>(matrix_free,
2256 first_selected_component,
2259 phi.emplace_back(std::make_unique<FEEvalType>(
2260 matrix_free, range, dof_no, quad_no, first_selected_component));
2265 data_cell.op_reinit = [](
auto &phi,
const unsigned batch) {
2266 if (phi.size() == 1)
2267 static_cast<FEEvalType &
>(*phi[0]).reinit(batch);
2271 data_cell.op_compute = [&](
auto &phi) {
2272 cell_operation(
static_cast<FEEvalType &
>(*phi[0]));
2275 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2278 data_face.dof_numbers = {dof_no, dof_no};
2279 data_face.quad_numbers = {quad_no, quad_no};
2280 data_face.n_components = {n_components, n_components};
2281 data_face.first_selected_components = {first_selected_component,
2282 first_selected_component};
2283 data_face.batch_type = {1, 2};
2285 data_face.op_create =
2286 [&](
const std::pair<unsigned int, unsigned int> &range) {
2288 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
2291 if (!internal::is_fe_nothing<true>(matrix_free,
2295 first_selected_component,
2299 !internal::is_fe_nothing<true>(matrix_free,
2303 first_selected_component,
2309 std::make_unique<FEFaceEvalType>(matrix_free,
2314 first_selected_component));
2316 std::make_unique<FEFaceEvalType>(matrix_free,
2321 first_selected_component));
2327 data_face.op_reinit = [](
auto &phi,
const unsigned batch) {
2328 if (phi.size() == 2)
2330 static_cast<FEFaceEvalType &
>(*phi[0]).
reinit(batch);
2331 static_cast<FEFaceEvalType &
>(*phi[1]).
reinit(batch);
2336 data_face.op_compute = [&](
auto &phi) {
2337 face_operation(
static_cast<FEFaceEvalType &
>(*phi[0]),
2338 static_cast<FEFaceEvalType &
>(*phi[1]));
2341 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2344 data_boundary.dof_numbers = {dof_no};
2345 data_boundary.quad_numbers = {quad_no};
2346 data_boundary.n_components = {n_components};
2347 data_boundary.first_selected_components = {first_selected_component};
2348 data_boundary.batch_type = {1};
2351 .op_create = [&](
const std::pair<unsigned int, unsigned int> &range) {
2353 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
2356 if (!internal::is_fe_nothing<true>(matrix_free,
2360 first_selected_component,
2364 phi.emplace_back(std::make_unique<FEFaceEvalType>(
2365 matrix_free, range,
true, dof_no, quad_no, first_selected_component));
2370 data_boundary.op_reinit = [](
auto &phi,
const unsigned batch) {
2371 if (phi.size() == 1)
2372 static_cast<FEFaceEvalType &
>(*phi[0]).reinit(batch);
2375 if (boundary_operation)
2376 data_boundary.op_compute = [&](
auto &phi) {
2377 boundary_operation(
static_cast<FEFaceEvalType &
>(*phi[0]));
2380 internal::compute_matrix(
2381 matrix_free, constraints_in, data_cell, data_face, data_boundary, matrix);
2388 typename VectorizedArrayType,
2394 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
2396 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2398 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2402 std::unique_ptr<AffineConstraints<typename MatrixType::value_type>>
2403 constraints_for_matrix;
2405 internal::create_new_affine_constraints_if_needed(
2406 matrix, constraints_in, constraints_for_matrix);
2408 const auto batch_operation =
2409 [&matrix_free, &constraints, &
matrix](
2410 auto &
data,
const std::pair<unsigned int, unsigned int> &range) {
2411 if (!
data.op_compute)
2414 auto phi =
data.op_create(range);
2416 const unsigned int n_blocks = phi.size();
2425 n_blocks, VectorizedArrayType::size());
2428 std::array<FullMatrix<typename MatrixType::value_type>,
2429 VectorizedArrayType::size()>>
2430 matrices(n_blocks, n_blocks);
2435 .get_fe(phi[b]->get_active_fe_index());
2437 const auto component_base =
2440 const auto component_in_base =
2441 data.first_selected_components[
b] -
2446 data.dof_numbers[b],
2447 data.quad_numbers[b],
2449 phi[b]->get_active_fe_index(),
2450 phi[b]->get_active_quadrature_index());
2453 shape_info.dofs_per_component_on_cell *
data.n_components[
b];
2455 dof_indices[
b].resize(fe.n_dofs_per_cell());
2457 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2458 dof_indices_mf[b][v].resize(dofs_per_cell[b]);
2460 lexicographic_numbering[
b].insert(
2461 lexicographic_numbering[b].
begin(),
2462 shape_info.lexicographic_numbering.begin() +
2463 component_in_base * shape_info.dofs_per_component_on_cell,
2464 shape_info.lexicographic_numbering.begin() +
2465 (component_in_base +
data.n_components[
b]) *
2466 shape_info.dofs_per_component_on_cell);
2469 for (
unsigned int bj = 0; bj <
n_blocks; ++bj)
2470 for (
unsigned int bi = 0; bi <
n_blocks; ++bi)
2471 std::fill_n(matrices[bi][bj].
begin(),
2472 VectorizedArrayType::size(),
2474 dofs_per_cell[bi], dofs_per_cell[bj]));
2476 for (
auto batch = range.first; batch < range.second; ++batch)
2478 data.op_reinit(phi, batch);
2480 const unsigned int n_filled_lanes =
2485 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
2489 (
data.batch_type[
b] == 0) ?
2490 (batch * VectorizedArrayType::size() + v) :
2491 ((
data.batch_type[
b] == 1) ?
2492 matrix_free.get_face_info(batch).cells_interior[v] :
2493 matrix_free.get_face_info(batch).cells_exterior[v]);
2498 data.dof_numbers[b]);
2502 cell_iterator->get_mg_dof_indices(dof_indices[b]);
2504 cell_iterator->get_dof_indices(dof_indices[b]);
2506 for (
unsigned int j = 0; j < dofs_per_cell[
b]; ++j)
2507 dof_indices_mf[b][v][j] =
2508 dof_indices[b][lexicographic_numbering[b][j]];
2511 for (
unsigned int bj = 0; bj <
n_blocks; ++bj)
2513 for (
unsigned int j = 0; j < dofs_per_cell[bj]; ++j)
2515 for (
unsigned int bi = 0; bi <
n_blocks; ++bi)
2516 for (
unsigned int i = 0; i < dofs_per_cell[bi]; ++i)
2517 phi[bi]->begin_dof_values()[i] =
2518 (bj == bi) ?
static_cast<Number
>(i == j) : 0.0;
2520 data.op_compute(phi);
2522 for (
unsigned int bi = 0; bi <
n_blocks; ++bi)
2523 for (
unsigned int i = 0; i < dofs_per_cell[bi]; ++i)
2524 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
2525 matrices[bi][bj][v](i, j) =
2526 phi[bi]->begin_dof_values()[i][v];
2529 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
2530 for (
unsigned int bi = 0; bi <
n_blocks; ++bi)
2538 matrices[bi][bi][v], dof_indices_mf[bi][v], matrix);
2541 matrices[bi][bj][v],
2542 dof_indices_mf[bi][v],
2543 dof_indices_mf[bj][v],
2549 const auto cell_operation_wrapped =
2550 [&](
const auto &,
auto &,
const auto &,
const auto range) {
2551 batch_operation(data_cell, range);
2554 const auto face_operation_wrapped =
2555 [&](
const auto &,
auto &,
const auto &,
const auto range) {
2556 batch_operation(data_face, range);
2559 const auto boundary_operation_wrapped =
2560 [&](
const auto &,
auto &,
const auto &,
const auto range) {
2561 batch_operation(data_boundary, range);
2564 if (data_face.op_compute || data_boundary.op_compute)
2566 matrix_free.template loop<MatrixType, MatrixType>(
2567 cell_operation_wrapped,
2568 face_operation_wrapped,
2569 boundary_operation_wrapped,
2574 matrix_free.template cell_loop<MatrixType, MatrixType>(
2575 cell_operation_wrapped, matrix, matrix);
2581 template <
typename CLASS,
2587 typename VectorizedArrayType,
2599 VectorizedArrayType> &) const,
2605 VectorizedArrayType> &,
2611 VectorizedArrayType> &)
2618 VectorizedArrayType> &)
2620 const CLASS *owning_class,
2621 const unsigned
int dof_no,
2622 const unsigned
int quad_no,
2623 const unsigned
int first_selected_component)
2630 VectorizedArrayType,
2635 [&](
auto &phi) { (owning_class->*cell_operation)(phi); },
2636 [&](
auto &phi_m,
auto &phi_p) {
2637 (owning_class->*face_operation)(phi_m, phi_p);
2639 [&](
auto &phi) { (owning_class->*boundary_operation)(phi); },
2642 first_selected_component);