41 template <
int dim,
typename Number,
bool is_face_>
54 std::function<std::vector<std::unique_ptr<FEEvalType>>(
55 const std::pair<unsigned int, unsigned int> &)>
57 std::function<void(std::vector<std::unique_ptr<FEEvalType>> &,
60 std::function<void(std::vector<std::unique_ptr<FEEvalType>> &)>
70 template <
int dim,
typename AdditionalData>
73 AdditionalData &additional_data);
93 typename VectorizedArrayType,
98 VectorType &diagonal_global,
104 VectorizedArrayType> &)>
106 const unsigned int dof_no = 0,
107 const unsigned int quad_no = 0,
108 const unsigned int first_selected_component = 0,
109 const unsigned int first_vector_component = 0);
120 typename QuadOperation>
125 const QuadOperation &quad_operation,
128 const unsigned int dof_no = 0,
129 const unsigned int quad_no = 0,
130 const unsigned int first_selected_component = 0,
131 const unsigned int first_vector_component = 0);
136 template <
typename CLASS,
142 typename VectorizedArrayType,
147 VectorType &diagonal_global,
153 VectorizedArrayType> &) const,
154 const CLASS *owning_class,
155 const unsigned
int dof_no = 0,
156 const unsigned
int quad_no = 0,
157 const unsigned
int first_selected_component = 0,
158 const unsigned
int first_vector_component = 0);
177 typename VectorizedArrayType,
182 VectorType &diagonal_global,
188 VectorizedArrayType> &)>
195 VectorizedArrayType> &,
201 VectorizedArrayType> &)>
208 VectorizedArrayType> &)>
210 const unsigned int dof_no = 0,
211 const unsigned int quad_no = 0,
212 const unsigned int first_selected_component = 0,
213 const unsigned int first_vector_component = 0);
220 template <
typename CLASS,
226 typename VectorizedArrayType,
231 VectorType &diagonal_global,
237 VectorizedArrayType> &) const,
243 VectorizedArrayType> &,
249 VectorizedArrayType> &)
256 VectorizedArrayType> &)
258 const CLASS *owning_class,
259 const unsigned
int dof_no = 0,
260 const unsigned
int quad_no = 0,
261 const unsigned
int first_selected_component = 0,
262 const unsigned
int first_vector_component = 0);
279 typename VectorizedArrayType,
291 VectorizedArrayType> &)>
293 const unsigned int dof_no = 0,
294 const unsigned int quad_no = 0,
295 const unsigned int first_selected_component = 0);
302 template <
typename CLASS,
308 typename VectorizedArrayType,
320 VectorizedArrayType> &) const,
321 const CLASS *owning_class,
322 const unsigned
int dof_no = 0,
323 const unsigned
int quad_no = 0,
324 const unsigned
int first_selected_component = 0);
335 typename VectorizedArrayType,
337 typename VectorType2>
347 VectorType &diagonal_global,
348 std::vector<VectorType2 *> &diagonal_global_components);
356 typename VectorizedArrayType,
387 typename VectorizedArrayType,
399 VectorizedArrayType> &)>
406 VectorizedArrayType> &,
412 VectorizedArrayType> &)>
419 VectorizedArrayType> &)>
421 const unsigned int dof_no = 0,
422 const unsigned int quad_no = 0,
423 const unsigned int first_selected_component = 0);
430 template <
typename CLASS,
436 typename VectorizedArrayType,
448 VectorizedArrayType> &) const,
454 VectorizedArrayType> &,
460 VectorizedArrayType> &)
467 VectorizedArrayType> &)
469 const CLASS *owning_class,
470 const unsigned
int dof_no = 0,
471 const unsigned
int quad_no = 0,
472 const unsigned
int first_selected_component = 0);
524 std::vector<unsigned int> valid_fe_indices;
526 const auto &fe_collection =
528 .get_fe_collection();
530 for (
unsigned int i = 0; i < fe_collection.size(); ++i)
531 if (fe_collection[i].n_dofs_per_cell() > 0)
532 valid_fe_indices.push_back(i);
546 template <
typename VectorTypeOut,
typename VectorTypeIn>
551 const VectorTypeIn &,
552 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
554 const VectorTypeIn &src,
555 const bool zero_dst_vector =
false)
const
557 const auto ebd_cell_operation = [&](
const auto &
matrix_free,
561 const auto category =
matrix_free.get_cell_range_category(range);
569 matrix_free->template cell_loop<VectorTypeOut, VectorTypeIn>(
570 ebd_cell_operation, dst, src, zero_dst_vector);
580 template <
typename VectorTypeOut,
typename VectorTypeIn>
585 const VectorTypeIn &,
586 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
590 const VectorTypeIn &,
591 const std::pair<unsigned int, unsigned int> &)> &face_operation,
595 const VectorTypeIn &,
596 const std::pair<unsigned int, unsigned int> &,
597 const bool)> &boundary_operation,
599 const VectorTypeIn &src,
600 const bool zero_dst_vector =
false)
const
602 const auto ebd_cell_operation = [&](
const auto &
matrix_free,
606 const auto category =
matrix_free.get_cell_range_category(range);
614 const auto ebd_internal_or_boundary_face_operation =
619 const auto category =
matrix_free.get_face_range_category(range);
621 const unsigned int type =
635 matrix_free->template loop<VectorTypeOut, VectorTypeIn>(
637 ebd_internal_or_boundary_face_operation,
638 ebd_internal_or_boundary_face_operation,
661 template <
int dim,
typename AdditionalData>
664 AdditionalData &additional_data)
667 const unsigned int level = additional_data.mg_level;
672 additional_data.cell_vectorization_category.assign(
675 additional_data.cell_vectorization_category.assign(tria.
n_active_cells(),
682 std::sort(bids.begin(), bids.end());
685 unsigned int n_bids = bids.size() + 1;
687 for (
unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
688 i++, offset = offset * n_bids)
692 const auto to_category = [&](
const auto &cell) {
693 unsigned int c_num = 0;
694 for (
unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
696 const auto face = cell->face(i);
697 if (face->at_boundary() && !cell->has_periodic_neighbor(i))
699 factors[i] * (1 + std::distance(bids.begin(),
700 std::find(bids.begin(),
702 face->boundary_id())));
711 if (cell->is_locally_owned())
713 .cell_vectorization_category[cell->active_cell_index()] =
721 if (cell->is_locally_owned_on_level())
722 additional_data.cell_vectorization_category[cell->index()] =
728 additional_data.hold_all_faces_to_owned_cells =
true;
729 additional_data.cell_vectorization_categories_strict =
true;
730 additional_data.mapping_update_flags_faces_by_cells =
731 additional_data.mapping_update_flags_inner_faces |
732 additional_data.mapping_update_flags_boundary_faces;
737 template <
typename Number>
740 std::vector<unsigned int> row_lid_to_gid;
741 std::vector<unsigned int> row;
742 std::vector<unsigned int> col;
743 std::vector<Number> val;
745 std::vector<unsigned int> inverse_lookup_rows;
746 std::vector<std::pair<unsigned int, unsigned int>> inverse_lookup_origins;
749 template <
int dim,
typename VectorizedArrayType,
bool is_face>
750 class ComputeDiagonalHelper
753 using FEEvaluationType =
756 using Number =
typename VectorizedArrayType::value_type;
757 static const unsigned int n_lanes = VectorizedArrayType::size();
759 ComputeDiagonalHelper()
761 , matrix_free(nullptr)
762 , dofs_per_component(0)
766 ComputeDiagonalHelper(
const ComputeDiagonalHelper &)
768 , matrix_free(nullptr)
769 , dofs_per_component(0)
775 FEEvaluationType &phi,
777 const unsigned int n_components)
781 if (dofs_per_component !=
782 phi.get_shape_info().dofs_per_component_on_cell)
784 locally_relevant_constraints_hn_map.
clear();
786 phi.get_shape_info().dofs_per_component_on_cell;
788 this->n_components = n_components;
789 this->dofs_per_cell = n_components * dofs_per_component;
791 this->matrix_free = &matrix_free;
795 reinit(
const unsigned int cell)
798 const auto &matrix_free = *this->matrix_free;
806 const unsigned int first_selected_component =
807 n_fe_components == 1 ? 0 : phi->get_first_selected_component();
809 this->n_lanes_filled =
818 std::vector<std::tuple<unsigned int, unsigned int, Number>>
819 locally_relevant_constraints, locally_relevant_constraints_tmp;
820 locally_relevant_constraints.reserve(dofs_per_cell);
821 std::vector<unsigned int> constraint_position;
822 std::vector<unsigned char> is_constrained_hn;
824 const std::array<unsigned int, n_lanes> &cells =
825 this->phi->get_cell_ids();
827 for (
unsigned int v = 0; v < n_lanes_filled; ++v)
832 const unsigned int *dof_indices;
833 unsigned int index_indicators, next_index_indicators;
835 const unsigned int start =
836 cells[v] * n_fe_components + first_selected_component;
838 dof_info.dof_indices.data() + dof_info.row_starts[start].first;
839 index_indicators = dof_info.row_starts[start].second;
840 next_index_indicators = dof_info.row_starts[start + 1].second;
844 locally_relevant_constraints.clear();
846 if (n_components == 1 || n_fe_components == 1)
848 unsigned int ind_local = 0;
849 for (; index_indicators != next_index_indicators;
850 ++index_indicators, ++ind_local)
852 const std::pair<unsigned short, unsigned short> indicator =
853 dof_info.constraint_indicator[index_indicators];
855 for (
unsigned int j = 0; j < indicator.first;
857 locally_relevant_constraints.emplace_back(ind_local,
861 dof_indices += indicator.first;
863 const Number *data_val =
865 const Number *end_pool =
868 for (; data_val != end_pool; ++data_val, ++dof_indices)
869 locally_relevant_constraints.emplace_back(ind_local,
876 for (; ind_local < dofs_per_component;
877 ++dof_indices, ++ind_local)
878 locally_relevant_constraints.emplace_back(ind_local,
889 for (
unsigned int comp = 0; comp < n_components; ++comp)
891 unsigned int ind_local = 0;
895 for (; index_indicators != next_index_indicators;
896 ++index_indicators, ++ind_local)
898 const std::pair<unsigned short, unsigned short>
900 dof_info.constraint_indicator[index_indicators];
903 for (
unsigned int j = 0; j < indicator.first;
905 locally_relevant_constraints.emplace_back(
906 comp * dofs_per_component + ind_local,
909 dof_indices += indicator.first;
911 const Number *data_val =
913 const Number *end_pool =
916 for (; data_val != end_pool; ++data_val, ++dof_indices)
917 locally_relevant_constraints.emplace_back(
918 comp * dofs_per_component + ind_local,
926 for (; ind_local < dofs_per_component;
927 ++dof_indices, ++ind_local)
928 locally_relevant_constraints.emplace_back(
929 comp * dofs_per_component + ind_local,
933 if (comp + 1 < n_components)
934 next_index_indicators =
935 dof_info.row_starts[start + comp + 2].second;
942 for (
unsigned int i = 1; i < locally_relevant_constraints.size();
944 Assert(std::get<0>(locally_relevant_constraints[i]) >=
945 std::get<0>(locally_relevant_constraints[i - 1]),
949 if (dof_info.hanging_node_constraint_masks.size() > 0 &&
950 dof_info.hanging_node_constraint_masks_comp.size() > 0 &&
951 dof_info.hanging_node_constraint_masks_comp
952 [phi->get_active_fe_index()][first_selected_component])
955 dof_info.hanging_node_constraint_masks[cells[v]];
958 if (mask != ::internal::MatrixFreeFunctions::
963 if (locally_relevant_constraints_hn_map.find(mask) ==
964 locally_relevant_constraints_hn_map.end())
965 fill_constraint_type_into_map(mask);
967 const auto &locally_relevant_constraints_hn =
968 locally_relevant_constraints_hn_map[
mask];
970 locally_relevant_constraints_tmp.clear();
971 if (locally_relevant_constraints_tmp.capacity() <
972 locally_relevant_constraints.size())
973 locally_relevant_constraints_tmp.reserve(
974 locally_relevant_constraints.size() +
975 locally_relevant_constraints_hn.size());
980 constraint_position.assign(dofs_per_cell,
982 for (
auto &a : locally_relevant_constraints)
983 if (constraint_position[
std::get<0>(a)] ==
985 constraint_position[
std::get<0>(a)] =
986 std::distance(locally_relevant_constraints.
data(),
988 is_constrained_hn.assign(dofs_per_cell,
false);
989 for (
auto &hn : locally_relevant_constraints_hn)
990 is_constrained_hn[
std::get<0>(hn)] = 1;
993 for (
const auto &a : locally_relevant_constraints)
994 if (is_constrained_hn[
std::get<0>(a)] == 0)
995 locally_relevant_constraints_tmp.push_back(a);
999 for (
const auto &hn : locally_relevant_constraints_hn)
1000 if (constraint_position[
std::get<1>(hn)] !=
1004 locally_relevant_constraints.size());
1005 auto other = locally_relevant_constraints.begin() +
1006 constraint_position[std::get<1>(hn)];
1009 for (; other != locally_relevant_constraints.end() &&
1010 std::get<0>(*other) == std::get<1>(hn);
1012 locally_relevant_constraints_tmp.emplace_back(
1014 std::get<1>(*other),
1015 std::get<2>(hn) * std::get<2>(*other));
1018 std::swap(locally_relevant_constraints,
1019 locally_relevant_constraints_tmp);
1024 std::sort(locally_relevant_constraints.begin(),
1025 locally_relevant_constraints.end(),
1026 [](
const auto &a,
const auto &b) {
1027 if (std::get<1>(a) < std::get<1>(b))
1029 return (std::get<1>(a) == std::get<1>(b)) &&
1030 (std::get<0>(a) < std::get<0>(b));
1034 auto &c_pool = c_pools[v];
1036 c_pool.row_lid_to_gid.clear();
1038 c_pool.row.push_back(0);
1042 if (locally_relevant_constraints.size() > 0)
1043 c_pool.row_lid_to_gid.emplace_back(
1044 std::get<1>(locally_relevant_constraints.front()));
1045 for (
const auto &j : locally_relevant_constraints)
1047 if (c_pool.row_lid_to_gid.back() != std::get<1>(j))
1049 c_pool.row_lid_to_gid.push_back(std::get<1>(j));
1050 c_pool.row.push_back(c_pool.val.size());
1053 c_pool.col.emplace_back(std::get<0>(j));
1054 c_pool.val.emplace_back(std::get<2>(j));
1057 if (c_pool.val.size() > 0)
1058 c_pool.row.push_back(c_pool.val.size());
1060 c_pool.inverse_lookup_rows.clear();
1061 c_pool.inverse_lookup_rows.resize(1 + dofs_per_cell);
1062 for (
const unsigned int i : c_pool.col)
1063 c_pool.inverse_lookup_rows[1 + i]++;
1065 std::partial_sum(c_pool.inverse_lookup_rows.begin(),
1066 c_pool.inverse_lookup_rows.end(),
1067 c_pool.inverse_lookup_rows.begin());
1071 c_pool.inverse_lookup_origins.resize(c_pool.col.size());
1072 std::vector<unsigned int> inverse_lookup_count(dofs_per_cell);
1073 for (
unsigned int row = 0; row < c_pool.row.size() - 1; ++row)
1074 for (
unsigned int col = c_pool.row[row];
1075 col < c_pool.row[row + 1];
1078 const unsigned int index = c_pool.col[col];
1079 c_pool.inverse_lookup_origins
1080 [c_pool.inverse_lookup_rows[
index] +
1081 inverse_lookup_count[
index]] = std::make_pair(row, col);
1082 ++inverse_lookup_count[
index];
1111 for (
unsigned int v = 0; v < n_lanes_filled; ++v)
1112 diagonals_local_constrained[v].assign(
1113 c_pools[v].row_lid_to_gid.size() *
1114 (n_fe_components == 1 ? n_components : 1),
1118 bool use_fast_path =
true;
1120 for (
unsigned int v = 0; v < n_lanes_filled; ++v)
1122 auto &c_pool = c_pools[v];
1124 for (
unsigned int i = 0; i < c_pool.row.size() - 1; ++i)
1126 if ((c_pool.row[i + 1] - c_pool.row[i]) > 1)
1128 use_fast_path =
false;
1131 else if (((c_pool.row[i + 1] - c_pool.row[i]) == 1) &&
1132 (c_pool.val[c_pool.row[i]] != 1.0))
1134 use_fast_path =
false;
1139 if (use_fast_path ==
false)
1143 this->has_simple_constraints_ = use_fast_path;
1147 fill_constraint_type_into_map(
1148 const ::internal::MatrixFreeFunctions::compressed_constraint_kind
1151 auto &constraints_hn = locally_relevant_constraints_hn_map[
mask];
1156 const unsigned int degree =
1157 phi->get_shape_info().data.front().fe_degree;
1162 values_dofs.resize(dofs_per_component);
1165 VectorizedArrayType::size()>
1167 constraint_mask.fill(::internal::MatrixFreeFunctions::
1169 constraint_mask[0] =
mask;
1171 for (
unsigned int i = 0; i < dofs_per_component; ++i)
1173 for (
unsigned int j = 0; j < dofs_per_component; ++j)
1174 values_dofs[j] = VectorizedArrayType();
1175 values_dofs[i] = Number(1);
1180 VectorizedArrayType>::apply(1,
1182 phi->get_shape_info(),
1185 values_dofs.data());
1187 const Number tolerance =
1189 std::numeric_limits<Number>::epsilon() * 16);
1190 for (
unsigned int j = 0; j < dofs_per_component; ++j)
1191 if (
std::abs(values_dofs[j][0]) > tolerance &&
1193 std::abs(values_dofs[j][0] - Number(1)) > tolerance))
1194 constraints_hn.emplace_back(j, i, values_dofs[j][0]);
1198 const unsigned int n_hn_constraints = constraints_hn.size();
1199 constraints_hn.resize(n_hn_constraints * n_components);
1201 for (
unsigned int c = 1; c < n_components; ++c)
1202 for (
unsigned int i = 0; i < n_hn_constraints; ++i)
1203 constraints_hn[c * n_hn_constraints + i] = std::make_tuple(
1204 std::get<0>(constraints_hn[i]) + c * dofs_per_component,
1205 std::get<1>(constraints_hn[i]) + c * dofs_per_component,
1206 std::get<2>(constraints_hn[i]));
1210 prepare_basis_vector(
const unsigned int i)
1217 VectorizedArrayType *dof_values = phi->begin_dof_values();
1218 for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1219 dof_values[j] = VectorizedArrayType();
1220 dof_values[i] = Number(1);
1226 VectorizedArrayType *dof_values = phi->begin_dof_values();
1227 for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1228 dof_values[j] = VectorizedArrayType();
1237 const unsigned int n_fe_components =
1238 phi->get_dof_info().start_components.back();
1239 const unsigned int comp =
1240 n_fe_components == 1 ? i / dofs_per_component : 0;
1241 const unsigned int i_comp =
1242 n_fe_components == 1 ? (i % dofs_per_component) : i;
1246 for (
unsigned int v = 0; v < n_lanes_filled; ++v)
1248 const auto &c_pool = c_pools[v];
1250 for (
unsigned int jj = c_pool.inverse_lookup_rows[i_comp];
1251 jj < c_pool.inverse_lookup_rows[i_comp + 1];
1254 const unsigned int j = c_pool.inverse_lookup_origins[jj].first;
1257 for (
unsigned int k = c_pool.row[j]; k < c_pool.row[j + 1]; ++k)
1258 temp += c_pool.val[k] *
1259 phi->begin_dof_values()[comp * dofs_per_component +
1263 diagonals_local_constrained
1264 [v][j + comp * c_pools[v].row_lid_to_gid.size()] +=
1265 temp * c_pool.val[c_pool.inverse_lookup_origins[jj].second];
1270 template <
typename VectorType>
1272 distribute_local_to_global(std::vector<VectorType *> &diagonal_global)
1275 const unsigned int n_fe_components =
1276 phi->get_dof_info().start_components.back();
1278 if (n_fe_components == 1)
1281 for (
unsigned int v = 0; v < n_lanes_filled; ++v)
1285 for (
unsigned int j = 0; j < c_pools[v].row.size() - 1; ++j)
1286 for (
unsigned int comp = 0;
1287 comp < (n_fe_components == 1 ?
1288 static_cast<unsigned int>(n_components) :
1292 *diagonal_global[n_fe_components == 1 ? comp : 0],
1293 c_pools[v].row_lid_to_gid[j],
1294 diagonals_local_constrained
1295 [v][j + comp * c_pools[v].row_lid_to_gid.
size()]);
1299 has_simple_constraints()
const
1301 return has_simple_constraints_;
1305 FEEvaluationType *phi;
1308 unsigned int dofs_per_component;
1309 unsigned int dofs_per_cell;
1310 unsigned int n_components;
1314 unsigned int n_lanes_filled;
1316 std::array<internal::LocalCSR<Number>, n_lanes> c_pools;
1321 std::array<std::vector<Number>, n_lanes> diagonals_local_constrained;
1325 std::vector<std::tuple<unsigned int, unsigned int, Number>>>
1326 locally_relevant_constraints_hn_map;
1331 bool has_simple_constraints_;
1334 template <
bool is_face,
1337 typename VectorizedArrayType>
1341 const std::pair<unsigned int, unsigned int> &range,
1342 const unsigned int dof_no,
1343 const unsigned int quad_no,
1344 const unsigned int first_selected_component,
1345 const unsigned int fe_degree,
1346 const unsigned int n_q_points_1d,
1347 const bool is_interior_face =
true)
1349 const unsigned int static_n_q_points =
1353 unsigned int active_fe_index = 0;
1356 else if (is_interior_face)
1361 const auto init_data = ::internal::
1362 extract_initialization_data<is_face, dim, Number, VectorizedArrayType>(
1365 first_selected_component,
1373 return init_data.shape_info->dofs_per_component_on_cell == 0;
1386 typename QuadOperation>
1387 class ComputeDiagonalCellAction
1390 ComputeDiagonalCellAction(
1391 const QuadOperation &quad_operation,
1394 : m_quad_operation(quad_operation)
1395 , m_evaluation_flags(evaluation_flags)
1396 , m_integration_flags(integration_flags)
1399 KOKKOS_FUNCTION
void
1405 FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number>
1408 *gpu_data =
data->precomputed_data;
1409 const int cell =
data->cell_index;
1411 constexpr int dofs_per_cell =
decltype(fe_eval)::tensor_dofs_per_cell;
1412 typename decltype(fe_eval)::value_type
1413 diagonal[dofs_per_cell / n_components] = {};
1414 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1416 const auto c = i % n_components;
1418 Kokkos::parallel_for(
1419 Kokkos::TeamThreadRange(
data->team_member,
1420 dofs_per_cell / n_components),
1422 typename decltype(fe_eval)::value_type val = {};
1424 if constexpr (n_components == 1)
1426 val = (i == j) ? 1 : 0;
1430 val[c] = (i / n_components == j) ? 1 : 0;
1433 fe_eval.submit_dof_value(val, j);
1436 data->team_member.team_barrier();
1438 Portable::internal::
1439 resolve_hanging_nodes<dim, fe_degree, false, Number>(
1443 Kokkos::subview(
data->shared_data->values, Kokkos::ALL, c));
1445 fe_eval.evaluate(m_evaluation_flags);
1446 fe_eval.apply_for_each_quad_point(m_quad_operation);
1447 fe_eval.integrate(m_integration_flags);
1449 Portable::internal::
1450 resolve_hanging_nodes<dim, fe_degree, true, Number>(
1454 Kokkos::subview(
data->shared_data->values, Kokkos::ALL, c));
1456 Kokkos::single(Kokkos::PerTeam(
data->team_member), [&] {
1457 if constexpr (n_components == 1)
1458 diagonal[i] = fe_eval.get_dof_value(i);
1460 diagonal[i / n_components][i % n_components] =
1461 fe_eval.get_dof_value(i / n_components)[i % n_components];
1464 data->team_member.team_barrier();
1467 Kokkos::single(Kokkos::PerTeam(
data->team_member), [&] {
1468 for (unsigned int i = 0; i < dofs_per_cell / n_components; ++i)
1469 fe_eval.submit_dof_value(diagonal[i], i);
1472 data->team_member.team_barrier();
1478 Kokkos::parallel_for(
1479 Kokkos::TeamThreadRange(
data->team_member, dofs_per_cell),
1481 dst[gpu_data->local_to_global(cell, i)] +=
1482 data->shared_data->values(i % (dofs_per_cell / n_components),
1483 i / (dofs_per_cell / n_components));
1488 Kokkos::parallel_for(
1489 Kokkos::TeamThreadRange(
data->team_member, dofs_per_cell),
1491 Kokkos::atomic_add(&dst[gpu_data->local_to_global(cell, i)],
1492 data->shared_data->values(
1493 i % (dofs_per_cell / n_components),
1494 i / (dofs_per_cell / n_components)));
1499 static constexpr unsigned int n_local_dofs =
1501 static constexpr unsigned int n_q_points =
1505 const QuadOperation m_quad_operation;
1520 typename QuadOperation>
1525 const QuadOperation &quad_operation,
1528 const unsigned int dof_no,
1529 const unsigned int quad_no,
1530 const unsigned int first_selected_component,
1531 const unsigned int first_vector_component)
1541 internal::ComputeDiagonalCellAction<dim,
1547 cell_action(quad_operation, evaluation_flags, integration_flags);
1549 matrix_free.
cell_loop(cell_action, dummy, diagonal_global);
1559 typename VectorizedArrayType,
1564 VectorType &diagonal_global,
1570 VectorizedArrayType> &)>
1572 const unsigned int dof_no,
1573 const unsigned int quad_no,
1574 const unsigned int first_selected_component,
1575 const unsigned int first_vector_component)
1582 VectorizedArrayType,
1590 first_selected_component,
1591 first_vector_component);
1594 template <
typename CLASS,
1600 typename VectorizedArrayType,
1605 VectorType &diagonal_global,
1611 VectorizedArrayType> &) const,
1612 const CLASS *owning_class,
1613 const unsigned
int dof_no,
1614 const unsigned
int quad_no,
1615 const unsigned
int first_selected_component,
1616 const unsigned
int first_vector_component)
1623 VectorizedArrayType,
1627 [&](
auto &phi) { (owning_class->*cell_operation)(phi); },
1630 first_selected_component,
1631 first_vector_component);
1639 typename VectorizedArrayType,
1644 VectorType &diagonal_global,
1650 VectorizedArrayType> &)>
1657 VectorizedArrayType> &,
1663 VectorizedArrayType> &)>
1670 VectorizedArrayType> &)>
1671 &boundary_operation,
1672 const unsigned int dof_no,
1673 const unsigned int quad_no,
1674 const unsigned int first_selected_component,
1675 const unsigned int first_vector_component)
1677 std::vector<typename ::internal::BlockVectorSelector<
1680 diagonal_global_components(n_components);
1682 for (
unsigned int d = 0;
d < n_components; ++
d)
1683 diagonal_global_components[d] = ::internal::
1685 get_vector_component(diagonal_global,
d + first_vector_component);
1687 const auto &dof_info = matrix_free.
get_dof_info(dof_no);
1689 if (dof_info.start_components.back() == 1)
1690 for (
unsigned int comp = 0; comp < n_components; ++comp)
1692 Assert(diagonal_global_components[comp] !=
nullptr,
1693 ExcMessage(
"The finite element underlying this FEEvaluation "
1694 "object is scalar, but you requested " +
1695 std::to_string(n_components) +
1696 " components via the template argument in "
1697 "FEEvaluation. In that case, you must pass an "
1698 "std::vector<VectorType> or a BlockVector to " +
1699 "read_dof_values and distribute_local_to_global."));
1701 *diagonal_global_components[comp], matrix_free, dof_info);
1706 *diagonal_global_components[0], matrix_free, dof_info);
1714 VectorizedArrayType>;
1721 VectorizedArrayType>;
1723 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
1726 data_cell.dof_numbers = {dof_no};
1727 data_cell.quad_numbers = {quad_no};
1728 data_cell.n_components = {n_components};
1729 data_cell.first_selected_components = {first_selected_component};
1730 data_cell.batch_type = {0};
1732 data_cell.op_create =
1733 [&](
const std::pair<unsigned int, unsigned int> &range) {
1735 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, false>>>
1738 if (!internal::is_fe_nothing<false>(matrix_free,
1742 first_selected_component,
1745 phi.emplace_back(std::make_unique<FEEvalType>(
1746 matrix_free, range, dof_no, quad_no, first_selected_component));
1751 data_cell.op_reinit = [](
auto &phi,
const unsigned batch) {
1752 if (phi.size() == 1)
1753 static_cast<FEEvalType &
>(*phi[0]).reinit(batch);
1757 data_cell.op_compute = [&](
auto &phi) {
1758 cell_operation(
static_cast<FEEvalType &
>(*phi[0]));
1761 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1764 data_face.dof_numbers = {dof_no, dof_no};
1765 data_face.quad_numbers = {quad_no, quad_no};
1766 data_face.n_components = {n_components, n_components};
1767 data_face.first_selected_components = {first_selected_component,
1768 first_selected_component};
1769 data_face.batch_type = {1, 2};
1771 data_face.op_create =
1772 [&](
const std::pair<unsigned int, unsigned int> &range) {
1774 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
1777 if (!internal::is_fe_nothing<true>(matrix_free,
1781 first_selected_component,
1785 !internal::is_fe_nothing<true>(matrix_free,
1789 first_selected_component,
1795 std::make_unique<FEFaceEvalType>(matrix_free,
1800 first_selected_component));
1802 std::make_unique<FEFaceEvalType>(matrix_free,
1807 first_selected_component));
1813 data_face.op_reinit = [](
auto &phi,
const unsigned batch) {
1814 if (phi.size() == 2)
1816 static_cast<FEFaceEvalType &
>(*phi[0]).
reinit(batch);
1817 static_cast<FEFaceEvalType &
>(*phi[1]).
reinit(batch);
1822 data_face.op_compute = [&](
auto &phi) {
1823 face_operation(
static_cast<FEFaceEvalType &
>(*phi[0]),
1824 static_cast<FEFaceEvalType &
>(*phi[1]));
1827 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1830 data_boundary.dof_numbers = {dof_no};
1831 data_boundary.quad_numbers = {quad_no};
1832 data_boundary.n_components = {n_components};
1833 data_boundary.first_selected_components = {first_selected_component};
1834 data_boundary.batch_type = {1};
1837 .op_create = [&](
const std::pair<unsigned int, unsigned int> &range) {
1839 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
1842 if (!internal::is_fe_nothing<true>(matrix_free,
1846 first_selected_component,
1850 phi.emplace_back(std::make_unique<FEFaceEvalType>(
1851 matrix_free, range,
true, dof_no, quad_no, first_selected_component));
1856 data_boundary.op_reinit = [](
auto &phi,
const unsigned batch) {
1857 if (phi.size() == 1)
1858 static_cast<FEFaceEvalType &
>(*phi[0]).reinit(batch);
1861 if (boundary_operation)
1862 data_boundary.op_compute = [&](
auto &phi) {
1863 boundary_operation(
static_cast<FEFaceEvalType &
>(*phi[0]));
1866 internal::compute_diagonal(matrix_free,
1871 diagonal_global_components);
1878 typename VectorizedArrayType,
1880 typename VectorType2>
1884 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
1886 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1888 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1890 VectorType &diagonal_global,
1891 std::vector<VectorType2 *> &diagonal_global_components)
1898 internal::ComputeDiagonalHelper<dim, VectorizedArrayType, false>;
1901 internal::ComputeDiagonalHelper<dim, VectorizedArrayType, true>;
1905 scratch_data_internal;
1908 const auto batch_operation =
1911 const std::pair<unsigned int, unsigned int> &range) {
1912 if (!
data.op_compute)
1915 auto phi =
data.op_create(range);
1917 const unsigned int n_blocks = phi.size();
1919 auto &helpers = scratch_data.
get();
1920 helpers.resize(n_blocks);
1923 helpers[b].initialize(*phi[b], matrix_free,
data.n_components[b]);
1925 for (
unsigned int batch = range.first; batch < range.second; ++batch)
1927 data.op_reinit(phi, batch);
1930 helpers[b].
reinit(batch);
1934 Assert(std::all_of(helpers.begin(),
1936 [](
const auto &helper) {
1937 return helper.has_simple_constraints();
1944 for (
unsigned int i = 0;
1945 i < phi[
b]->get_shape_info().dofs_per_component_on_cell *
1946 data.n_components[b];
1949 for (
unsigned int bb = 0; bb <
n_blocks; ++bb)
1951 helpers[bb].prepare_basis_vector(i);
1953 helpers[bb].zero_basis_vector();
1955 data.op_compute(phi);
1956 helpers[
b].submit();
1959 helpers[
b].distribute_local_to_global(
1960 diagonal_global_components);
1965 const auto cell_operation_wrapped =
1966 [&](
const auto &,
auto &,
const auto &,
const auto range) {
1967 batch_operation(data_cell, scratch_data, range);
1970 const auto face_operation_wrapped =
1971 [&](
const auto &,
auto &,
const auto &,
const auto range) {
1972 batch_operation(data_face, scratch_data_internal, range);
1975 const auto boundary_operation_wrapped =
1976 [&](
const auto &,
auto &,
const auto &,
const auto range) {
1977 batch_operation(data_boundary, scratch_data_bc, range);
1980 if (data_face.op_compute || data_boundary.op_compute)
1981 matrix_free.template loop<VectorType, int>(cell_operation_wrapped,
1982 face_operation_wrapped,
1983 boundary_operation_wrapped,
1988 matrix_free.template cell_loop<VectorType, int>(cell_operation_wrapped,
1995 template <
typename CLASS,
2001 typename VectorizedArrayType,
2006 VectorType &diagonal_global,
2012 VectorizedArrayType> &) const,
2018 VectorizedArrayType> &,
2024 VectorizedArrayType> &)
2031 VectorizedArrayType> &)
2033 const CLASS *owning_class,
2034 const unsigned
int dof_no,
2035 const unsigned
int quad_no,
2036 const unsigned
int first_selected_component,
2037 const unsigned
int first_vector_component)
2044 VectorizedArrayType,
2048 [&](
auto &phi) { (owning_class->*cell_operation)(phi); },
2049 [&](
auto &phi_m,
auto &phi_p) {
2050 (owning_class->*face_operation)(phi_m, phi_p);
2052 [&](
auto &phi) { (owning_class->*boundary_operation)(phi); },
2055 first_selected_component,
2056 first_vector_component);
2068 std::enable_if_t<std::is_same_v<
2069 std::remove_const_t<
2070 std::remove_reference_t<typename MatrixType::value_type>>,
2071 std::remove_const_t<std::remove_reference_t<Number>>>> * =
nullptr>
2073 create_new_affine_constraints_if_needed(
2089 std::enable_if_t<!std::is_same_v<
2090 std::remove_const_t<
2091 std::remove_reference_t<typename MatrixType::value_type>>,
2092 std::remove_const_t<std::remove_reference_t<Number>>>> * =
nullptr>
2094 create_new_affine_constraints_if_needed(
2101 std::make_unique<AffineConstraints<typename MatrixType::value_type>>();
2102 new_constraints->
copy_from(constraints);
2104 return *new_constraints;
2113 typename VectorizedArrayType,
2125 VectorizedArrayType> &)>
2127 const unsigned int dof_no,
2128 const unsigned int quad_no,
2129 const unsigned int first_selected_component)
2136 VectorizedArrayType,
2145 first_selected_component);
2148 template <
typename CLASS,
2154 typename VectorizedArrayType,
2166 VectorizedArrayType> &) const,
2167 const CLASS *owning_class,
2168 const unsigned
int dof_no,
2169 const unsigned
int quad_no,
2170 const unsigned
int first_selected_component)
2177 VectorizedArrayType,
2182 [&](
auto &phi) { (owning_class->*cell_operation)(phi); },
2185 first_selected_component);
2193 typename VectorizedArrayType,
2205 VectorizedArrayType> &)>
2212 VectorizedArrayType> &,
2218 VectorizedArrayType> &)>
2225 VectorizedArrayType> &)>
2226 &boundary_operation,
2227 const unsigned int dof_no,
2228 const unsigned int quad_no,
2229 const unsigned int first_selected_component)
2236 VectorizedArrayType>;
2243 VectorizedArrayType>;
2245 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
2248 data_cell.dof_numbers = {dof_no};
2249 data_cell.quad_numbers = {quad_no};
2250 data_cell.n_components = {n_components};
2251 data_cell.first_selected_components = {first_selected_component};
2252 data_cell.batch_type = {0};
2254 data_cell.op_create =
2255 [&](
const std::pair<unsigned int, unsigned int> &range) {
2257 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, false>>>
2260 if (!internal::is_fe_nothing<false>(matrix_free,
2264 first_selected_component,
2267 phi.emplace_back(std::make_unique<FEEvalType>(
2268 matrix_free, range, dof_no, quad_no, first_selected_component));
2273 data_cell.op_reinit = [](
auto &phi,
const unsigned batch) {
2274 if (phi.size() == 1)
2275 static_cast<FEEvalType &
>(*phi[0]).reinit(batch);
2279 data_cell.op_compute = [&](
auto &phi) {
2280 cell_operation(
static_cast<FEEvalType &
>(*phi[0]));
2283 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2286 data_face.dof_numbers = {dof_no, dof_no};
2287 data_face.quad_numbers = {quad_no, quad_no};
2288 data_face.n_components = {n_components, n_components};
2289 data_face.first_selected_components = {first_selected_component,
2290 first_selected_component};
2291 data_face.batch_type = {1, 2};
2293 data_face.op_create =
2294 [&](
const std::pair<unsigned int, unsigned int> &range) {
2296 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
2299 if (!internal::is_fe_nothing<true>(matrix_free,
2303 first_selected_component,
2307 !internal::is_fe_nothing<true>(matrix_free,
2311 first_selected_component,
2317 std::make_unique<FEFaceEvalType>(matrix_free,
2322 first_selected_component));
2324 std::make_unique<FEFaceEvalType>(matrix_free,
2329 first_selected_component));
2335 data_face.op_reinit = [](
auto &phi,
const unsigned batch) {
2336 if (phi.size() == 2)
2338 static_cast<FEFaceEvalType &
>(*phi[0]).
reinit(batch);
2339 static_cast<FEFaceEvalType &
>(*phi[1]).
reinit(batch);
2344 data_face.op_compute = [&](
auto &phi) {
2345 face_operation(
static_cast<FEFaceEvalType &
>(*phi[0]),
2346 static_cast<FEFaceEvalType &
>(*phi[1]));
2349 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2352 data_boundary.dof_numbers = {dof_no};
2353 data_boundary.quad_numbers = {quad_no};
2354 data_boundary.n_components = {n_components};
2355 data_boundary.first_selected_components = {first_selected_component};
2356 data_boundary.batch_type = {1};
2359 .op_create = [&](
const std::pair<unsigned int, unsigned int> &range) {
2361 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
2364 if (!internal::is_fe_nothing<true>(matrix_free,
2368 first_selected_component,
2372 phi.emplace_back(std::make_unique<FEFaceEvalType>(
2373 matrix_free, range,
true, dof_no, quad_no, first_selected_component));
2378 data_boundary.op_reinit = [](
auto &phi,
const unsigned batch) {
2379 if (phi.size() == 1)
2380 static_cast<FEFaceEvalType &
>(*phi[0]).reinit(batch);
2383 if (boundary_operation)
2384 data_boundary.op_compute = [&](
auto &phi) {
2385 boundary_operation(
static_cast<FEFaceEvalType &
>(*phi[0]));
2388 internal::compute_matrix(
2389 matrix_free, constraints_in, data_cell, data_face, data_boundary, matrix);
2396 typename VectorizedArrayType,
2402 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
2404 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2406 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2410 std::unique_ptr<AffineConstraints<typename MatrixType::value_type>>
2411 constraints_for_matrix;
2413 internal::create_new_affine_constraints_if_needed(
2414 matrix, constraints_in, constraints_for_matrix);
2416 const auto batch_operation =
2417 [&matrix_free, &constraints, &
matrix](
2418 auto &
data,
const std::pair<unsigned int, unsigned int> &range) {
2419 if (!
data.op_compute)
2422 auto phi =
data.op_create(range);
2424 const unsigned int n_blocks = phi.size();
2433 n_blocks, VectorizedArrayType::size());
2436 std::array<FullMatrix<typename MatrixType::value_type>,
2437 VectorizedArrayType::size()>>
2438 matrices(n_blocks, n_blocks);
2443 .get_fe(phi[b]->get_active_fe_index());
2445 const auto component_base =
2448 const auto component_in_base =
2449 data.first_selected_components[
b] -
2454 data.dof_numbers[b],
2455 data.quad_numbers[b],
2457 phi[b]->get_active_fe_index(),
2458 phi[b]->get_active_quadrature_index());
2461 shape_info.dofs_per_component_on_cell *
data.n_components[
b];
2463 dof_indices[
b].resize(fe.n_dofs_per_cell());
2465 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2466 dof_indices_mf[b][v].resize(dofs_per_cell[b]);
2468 lexicographic_numbering[
b].insert(
2469 lexicographic_numbering[b].
begin(),
2470 shape_info.lexicographic_numbering.begin() +
2471 component_in_base * shape_info.dofs_per_component_on_cell,
2472 shape_info.lexicographic_numbering.begin() +
2473 (component_in_base +
data.n_components[
b]) *
2474 shape_info.dofs_per_component_on_cell);
2477 for (
unsigned int bj = 0; bj <
n_blocks; ++bj)
2478 for (
unsigned int bi = 0; bi <
n_blocks; ++bi)
2479 std::fill_n(matrices[bi][bj].
begin(),
2480 VectorizedArrayType::size(),
2482 dofs_per_cell[bi], dofs_per_cell[bj]));
2484 for (
auto batch = range.first; batch < range.second; ++batch)
2486 data.op_reinit(phi, batch);
2488 const unsigned int n_filled_lanes =
2493 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
2497 (
data.batch_type[
b] == 0) ?
2498 (batch * VectorizedArrayType::size() + v) :
2499 ((
data.batch_type[
b] == 1) ?
2500 matrix_free.get_face_info(batch).cells_interior[v] :
2501 matrix_free.get_face_info(batch).cells_exterior[v]);
2506 data.dof_numbers[b]);
2510 cell_iterator->get_mg_dof_indices(dof_indices[b]);
2512 cell_iterator->get_dof_indices(dof_indices[b]);
2514 for (
unsigned int j = 0; j < dofs_per_cell[
b]; ++j)
2515 dof_indices_mf[b][v][j] =
2516 dof_indices[b][lexicographic_numbering[b][j]];
2519 for (
unsigned int bj = 0; bj <
n_blocks; ++bj)
2521 for (
unsigned int j = 0; j < dofs_per_cell[bj]; ++j)
2523 for (
unsigned int bi = 0; bi <
n_blocks; ++bi)
2524 for (
unsigned int i = 0; i < dofs_per_cell[bi]; ++i)
2525 phi[bi]->begin_dof_values()[i] =
2526 (bj == bi) ?
static_cast<Number
>(i == j) : 0.0;
2528 data.op_compute(phi);
2530 for (
unsigned int bi = 0; bi <
n_blocks; ++bi)
2531 for (
unsigned int i = 0; i < dofs_per_cell[bi]; ++i)
2532 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
2533 matrices[bi][bj][v](i, j) =
2534 phi[bi]->begin_dof_values()[i][v];
2537 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
2538 for (
unsigned int bi = 0; bi <
n_blocks; ++bi)
2546 matrices[bi][bi][v], dof_indices_mf[bi][v], matrix);
2549 matrices[bi][bj][v],
2550 dof_indices_mf[bi][v],
2551 dof_indices_mf[bj][v],
2557 const auto cell_operation_wrapped =
2558 [&](
const auto &,
auto &,
const auto &,
const auto range) {
2559 batch_operation(data_cell, range);
2562 const auto face_operation_wrapped =
2563 [&](
const auto &,
auto &,
const auto &,
const auto range) {
2564 batch_operation(data_face, range);
2567 const auto boundary_operation_wrapped =
2568 [&](
const auto &,
auto &,
const auto &,
const auto range) {
2569 batch_operation(data_boundary, range);
2572 if (data_face.op_compute || data_boundary.op_compute)
2574 matrix_free.template loop<MatrixType, MatrixType>(
2575 cell_operation_wrapped,
2576 face_operation_wrapped,
2577 boundary_operation_wrapped,
2582 matrix_free.template cell_loop<MatrixType, MatrixType>(
2583 cell_operation_wrapped, matrix, matrix);
2589 template <
typename CLASS,
2595 typename VectorizedArrayType,
2607 VectorizedArrayType> &) const,
2613 VectorizedArrayType> &,
2619 VectorizedArrayType> &)
2626 VectorizedArrayType> &)
2628 const CLASS *owning_class,
2629 const unsigned
int dof_no,
2630 const unsigned
int quad_no,
2631 const unsigned
int first_selected_component)
2638 VectorizedArrayType,
2643 [&](
auto &phi) { (owning_class->*cell_operation)(phi); },
2644 [&](
auto &phi_m,
auto &phi_p) {
2645 (owning_class->*face_operation)(phi_m, phi_p);
2647 [&](
auto &phi) { (owning_class->*boundary_operation)(phi); },
2650 first_selected_component);