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)
1363 const auto init_data = ::internal::
1364 extract_initialization_data<is_face, dim, Number, VectorizedArrayType>(
1367 first_selected_component,
1375 return init_data.shape_info->dofs_per_component_on_cell == 0;
1388 typename QuadOperation>
1389 class ComputeDiagonalCellAction
1392 ComputeDiagonalCellAction(
1393 const QuadOperation &quad_operation,
1396 : m_quad_operation(quad_operation)
1397 , m_evaluation_flags(evaluation_flags)
1398 , m_integration_flags(integration_flags)
1401 KOKKOS_FUNCTION
void
1407 FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number>
1410 *gpu_data =
data->precomputed_data;
1411 const int cell =
data->cell_index;
1413 constexpr int dofs_per_cell =
decltype(fe_eval)::tensor_dofs_per_cell;
1414 typename decltype(fe_eval)::value_type
1415 diagonal[dofs_per_cell / n_components] = {};
1416 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1418 const auto c = i % n_components;
1420 Kokkos::parallel_for(
1421 Kokkos::TeamThreadRange(
data->team_member,
1422 dofs_per_cell / n_components),
1424 typename decltype(fe_eval)::value_type val = {};
1426 if constexpr (n_components == 1)
1428 val = (i == j) ? 1 : 0;
1432 val[c] = (i / n_components == j) ? 1 : 0;
1435 fe_eval.submit_dof_value(val, j);
1438 data->team_member.team_barrier();
1440 Portable::internal::
1441 resolve_hanging_nodes<dim, fe_degree, false, Number>(
1445 Kokkos::subview(
data->shared_data->values, Kokkos::ALL, c));
1447 fe_eval.evaluate(m_evaluation_flags);
1448 fe_eval.apply_for_each_quad_point(m_quad_operation);
1449 fe_eval.integrate(m_integration_flags);
1451 Portable::internal::
1452 resolve_hanging_nodes<dim, fe_degree, true, Number>(
1456 Kokkos::subview(
data->shared_data->values, Kokkos::ALL, c));
1458 Kokkos::single(Kokkos::PerTeam(
data->team_member), [&] {
1459 if constexpr (n_components == 1)
1460 diagonal[i] = fe_eval.get_dof_value(i);
1462 diagonal[i / n_components][i % n_components] =
1463 fe_eval.get_dof_value(i / n_components)[i % n_components];
1466 data->team_member.team_barrier();
1469 Kokkos::single(Kokkos::PerTeam(
data->team_member), [&] {
1470 for (unsigned int i = 0; i < dofs_per_cell / n_components; ++i)
1471 fe_eval.submit_dof_value(diagonal[i], i);
1474 data->team_member.team_barrier();
1480 Kokkos::parallel_for(
1481 Kokkos::TeamThreadRange(
data->team_member, dofs_per_cell),
1483 dst[gpu_data->local_to_global(i, cell)] +=
1484 data->shared_data->values(i % (dofs_per_cell / n_components),
1485 i / (dofs_per_cell / n_components));
1490 Kokkos::parallel_for(
1491 Kokkos::TeamThreadRange(
data->team_member, dofs_per_cell),
1493 Kokkos::atomic_add(&dst[gpu_data->local_to_global(i, cell)],
1494 data->shared_data->values(
1495 i % (dofs_per_cell / n_components),
1496 i / (dofs_per_cell / n_components)));
1501 static constexpr unsigned int n_local_dofs =
1503 static constexpr unsigned int n_q_points =
1507 const QuadOperation m_quad_operation;
1522 typename QuadOperation>
1527 const QuadOperation &quad_operation,
1530 const unsigned int dof_no,
1531 const unsigned int quad_no,
1532 const unsigned int first_selected_component,
1533 const unsigned int first_vector_component)
1543 internal::ComputeDiagonalCellAction<dim,
1549 cell_action(quad_operation, evaluation_flags, integration_flags);
1551 matrix_free.
cell_loop(cell_action, dummy, diagonal_global);
1561 typename VectorizedArrayType,
1566 VectorType &diagonal_global,
1572 VectorizedArrayType> &)>
1574 const unsigned int dof_no,
1575 const unsigned int quad_no,
1576 const unsigned int first_selected_component,
1577 const unsigned int first_vector_component)
1584 VectorizedArrayType,
1592 first_selected_component,
1593 first_vector_component);
1596 template <
typename CLASS,
1602 typename VectorizedArrayType,
1607 VectorType &diagonal_global,
1613 VectorizedArrayType> &) const,
1614 const CLASS *owning_class,
1615 const unsigned
int dof_no,
1616 const unsigned
int quad_no,
1617 const unsigned
int first_selected_component,
1618 const unsigned
int first_vector_component)
1625 VectorizedArrayType,
1629 [&](
auto &phi) { (owning_class->*cell_operation)(phi); },
1632 first_selected_component,
1633 first_vector_component);
1641 typename VectorizedArrayType,
1646 VectorType &diagonal_global,
1652 VectorizedArrayType> &)>
1659 VectorizedArrayType> &,
1665 VectorizedArrayType> &)>
1672 VectorizedArrayType> &)>
1673 &boundary_operation,
1674 const unsigned int dof_no,
1675 const unsigned int quad_no,
1676 const unsigned int first_selected_component,
1677 const unsigned int first_vector_component)
1679 std::vector<typename ::internal::BlockVectorSelector<
1682 diagonal_global_components(n_components);
1684 for (
unsigned int d = 0;
d < n_components; ++
d)
1685 diagonal_global_components[d] = ::internal::
1687 get_vector_component(diagonal_global,
d + first_vector_component);
1689 const auto &dof_info = matrix_free.
get_dof_info(dof_no);
1691 if (dof_info.start_components.back() == 1)
1692 for (
unsigned int comp = 0; comp < n_components; ++comp)
1694 Assert(diagonal_global_components[comp] !=
nullptr,
1695 ExcMessage(
"The finite element underlying this FEEvaluation "
1696 "object is scalar, but you requested " +
1697 std::to_string(n_components) +
1698 " components via the template argument in "
1699 "FEEvaluation. In that case, you must pass an "
1700 "std::vector<VectorType> or a BlockVector to " +
1701 "read_dof_values and distribute_local_to_global."));
1703 *diagonal_global_components[comp], matrix_free, dof_info);
1708 *diagonal_global_components[0], matrix_free, dof_info);
1716 VectorizedArrayType>;
1723 VectorizedArrayType>;
1725 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
1728 data_cell.dof_numbers = {dof_no};
1729 data_cell.quad_numbers = {quad_no};
1730 data_cell.n_components = {n_components};
1731 data_cell.first_selected_components = {first_selected_component};
1732 data_cell.batch_type = {0};
1734 data_cell.op_create =
1735 [&](
const std::pair<unsigned int, unsigned int> &range) {
1737 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, false>>>
1740 if (!internal::is_fe_nothing<false>(matrix_free,
1744 first_selected_component,
1747 phi.emplace_back(std::make_unique<FEEvalType>(
1748 matrix_free, range, dof_no, quad_no, first_selected_component));
1753 data_cell.op_reinit = [](
auto &phi,
const unsigned batch) {
1754 if (phi.size() == 1)
1755 static_cast<FEEvalType &
>(*phi[0]).reinit(batch);
1759 data_cell.op_compute = [&](
auto &phi) {
1760 cell_operation(
static_cast<FEEvalType &
>(*phi[0]));
1763 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1766 data_face.dof_numbers = {dof_no, dof_no};
1767 data_face.quad_numbers = {quad_no, quad_no};
1768 data_face.n_components = {n_components, n_components};
1769 data_face.first_selected_components = {first_selected_component,
1770 first_selected_component};
1771 data_face.batch_type = {1, 2};
1773 data_face.op_create =
1774 [&](
const std::pair<unsigned int, unsigned int> &range) {
1776 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
1779 if (!internal::is_fe_nothing<true>(matrix_free,
1783 first_selected_component,
1787 !internal::is_fe_nothing<true>(matrix_free,
1791 first_selected_component,
1797 std::make_unique<FEFaceEvalType>(matrix_free,
1802 first_selected_component));
1804 std::make_unique<FEFaceEvalType>(matrix_free,
1809 first_selected_component));
1815 data_face.op_reinit = [](
auto &phi,
const unsigned batch) {
1816 if (phi.size() == 2)
1818 static_cast<FEFaceEvalType &
>(*phi[0]).
reinit(batch);
1819 static_cast<FEFaceEvalType &
>(*phi[1]).
reinit(batch);
1824 data_face.op_compute = [&](
auto &phi) {
1825 face_operation(
static_cast<FEFaceEvalType &
>(*phi[0]),
1826 static_cast<FEFaceEvalType &
>(*phi[1]));
1829 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1832 data_boundary.dof_numbers = {dof_no};
1833 data_boundary.quad_numbers = {quad_no};
1834 data_boundary.n_components = {n_components};
1835 data_boundary.first_selected_components = {first_selected_component};
1836 data_boundary.batch_type = {1};
1839 .op_create = [&](
const std::pair<unsigned int, unsigned int> &range) {
1841 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
1844 if (!internal::is_fe_nothing<true>(matrix_free,
1848 first_selected_component,
1852 phi.emplace_back(std::make_unique<FEFaceEvalType>(
1853 matrix_free, range,
true, dof_no, quad_no, first_selected_component));
1858 data_boundary.op_reinit = [](
auto &phi,
const unsigned batch) {
1859 if (phi.size() == 1)
1860 static_cast<FEFaceEvalType &
>(*phi[0]).reinit(batch);
1863 if (boundary_operation)
1864 data_boundary.op_compute = [&](
auto &phi) {
1865 boundary_operation(
static_cast<FEFaceEvalType &
>(*phi[0]));
1868 internal::compute_diagonal(matrix_free,
1873 diagonal_global_components);
1880 typename VectorizedArrayType,
1882 typename VectorType2>
1886 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
1888 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1890 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1892 VectorType &diagonal_global,
1893 std::vector<VectorType2 *> &diagonal_global_components)
1900 internal::ComputeDiagonalHelper<dim, VectorizedArrayType, false>;
1903 internal::ComputeDiagonalHelper<dim, VectorizedArrayType, true>;
1907 scratch_data_internal;
1910 const auto batch_operation =
1913 const std::pair<unsigned int, unsigned int> &range) {
1914 if (!
data.op_compute)
1917 auto phi =
data.op_create(range);
1919 const unsigned int n_blocks = phi.size();
1921 auto &helpers = scratch_data.
get();
1922 helpers.resize(n_blocks);
1925 helpers[b].initialize(*phi[b], matrix_free,
data.n_components[b]);
1927 for (
unsigned int batch = range.first; batch < range.second; ++batch)
1929 data.op_reinit(phi, batch);
1932 helpers[b].
reinit(batch);
1936 Assert(std::all_of(helpers.begin(),
1938 [](
const auto &helper) {
1939 return helper.has_simple_constraints();
1946 for (
unsigned int i = 0;
1947 i < phi[
b]->get_shape_info().dofs_per_component_on_cell *
1948 data.n_components[b];
1951 for (
unsigned int bb = 0; bb <
n_blocks; ++bb)
1953 helpers[bb].prepare_basis_vector(i);
1955 helpers[bb].zero_basis_vector();
1957 data.op_compute(phi);
1958 helpers[
b].submit();
1961 helpers[
b].distribute_local_to_global(
1962 diagonal_global_components);
1967 const auto cell_operation_wrapped =
1968 [&](
const auto &,
auto &,
const auto &,
const auto range) {
1969 batch_operation(data_cell, scratch_data, range);
1972 const auto face_operation_wrapped =
1973 [&](
const auto &,
auto &,
const auto &,
const auto range) {
1974 batch_operation(data_face, scratch_data_internal, range);
1977 const auto boundary_operation_wrapped =
1978 [&](
const auto &,
auto &,
const auto &,
const auto range) {
1979 batch_operation(data_boundary, scratch_data_bc, range);
1982 if (data_face.op_compute || data_boundary.op_compute)
1983 matrix_free.template loop<VectorType, int>(cell_operation_wrapped,
1984 face_operation_wrapped,
1985 boundary_operation_wrapped,
1990 matrix_free.template cell_loop<VectorType, int>(cell_operation_wrapped,
1997 template <
typename CLASS,
2003 typename VectorizedArrayType,
2008 VectorType &diagonal_global,
2014 VectorizedArrayType> &) const,
2020 VectorizedArrayType> &,
2026 VectorizedArrayType> &)
2033 VectorizedArrayType> &)
2035 const CLASS *owning_class,
2036 const unsigned
int dof_no,
2037 const unsigned
int quad_no,
2038 const unsigned
int first_selected_component,
2039 const unsigned
int first_vector_component)
2046 VectorizedArrayType,
2050 [&](
auto &phi) { (owning_class->*cell_operation)(phi); },
2051 [&](
auto &phi_m,
auto &phi_p) {
2052 (owning_class->*face_operation)(phi_m, phi_p);
2054 [&](
auto &phi) { (owning_class->*boundary_operation)(phi); },
2057 first_selected_component,
2058 first_vector_component);
2070 std::enable_if_t<std::is_same_v<
2071 std::remove_const_t<
2072 std::remove_reference_t<typename MatrixType::value_type>>,
2073 std::remove_const_t<std::remove_reference_t<Number>>>> * =
nullptr>
2075 create_new_affine_constraints_if_needed(
2091 std::enable_if_t<!std::is_same_v<
2092 std::remove_const_t<
2093 std::remove_reference_t<typename MatrixType::value_type>>,
2094 std::remove_const_t<std::remove_reference_t<Number>>>> * =
nullptr>
2096 create_new_affine_constraints_if_needed(
2103 std::make_unique<AffineConstraints<typename MatrixType::value_type>>();
2104 new_constraints->
copy_from(constraints);
2106 return *new_constraints;
2115 typename VectorizedArrayType,
2127 VectorizedArrayType> &)>
2129 const unsigned int dof_no,
2130 const unsigned int quad_no,
2131 const unsigned int first_selected_component)
2138 VectorizedArrayType,
2147 first_selected_component);
2150 template <
typename CLASS,
2156 typename VectorizedArrayType,
2168 VectorizedArrayType> &) const,
2169 const CLASS *owning_class,
2170 const unsigned
int dof_no,
2171 const unsigned
int quad_no,
2172 const unsigned
int first_selected_component)
2179 VectorizedArrayType,
2184 [&](
auto &phi) { (owning_class->*cell_operation)(phi); },
2187 first_selected_component);
2195 typename VectorizedArrayType,
2207 VectorizedArrayType> &)>
2214 VectorizedArrayType> &,
2220 VectorizedArrayType> &)>
2227 VectorizedArrayType> &)>
2228 &boundary_operation,
2229 const unsigned int dof_no,
2230 const unsigned int quad_no,
2231 const unsigned int first_selected_component)
2238 VectorizedArrayType>;
2245 VectorizedArrayType>;
2247 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
2250 data_cell.dof_numbers = {dof_no};
2251 data_cell.quad_numbers = {quad_no};
2252 data_cell.n_components = {n_components};
2253 data_cell.first_selected_components = {first_selected_component};
2254 data_cell.batch_type = {0};
2256 data_cell.op_create =
2257 [&](
const std::pair<unsigned int, unsigned int> &range) {
2259 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, false>>>
2262 if (!internal::is_fe_nothing<false>(matrix_free,
2266 first_selected_component,
2269 phi.emplace_back(std::make_unique<FEEvalType>(
2270 matrix_free, range, dof_no, quad_no, first_selected_component));
2275 data_cell.op_reinit = [](
auto &phi,
const unsigned batch) {
2276 if (phi.size() == 1)
2277 static_cast<FEEvalType &
>(*phi[0]).reinit(batch);
2281 data_cell.op_compute = [&](
auto &phi) {
2282 cell_operation(
static_cast<FEEvalType &
>(*phi[0]));
2285 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2288 data_face.dof_numbers = {dof_no, dof_no};
2289 data_face.quad_numbers = {quad_no, quad_no};
2290 data_face.n_components = {n_components, n_components};
2291 data_face.first_selected_components = {first_selected_component,
2292 first_selected_component};
2293 data_face.batch_type = {1, 2};
2295 data_face.op_create =
2296 [&](
const std::pair<unsigned int, unsigned int> &range) {
2298 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
2301 if (!internal::is_fe_nothing<true>(matrix_free,
2305 first_selected_component,
2309 !internal::is_fe_nothing<true>(matrix_free,
2313 first_selected_component,
2319 std::make_unique<FEFaceEvalType>(matrix_free,
2324 first_selected_component));
2326 std::make_unique<FEFaceEvalType>(matrix_free,
2331 first_selected_component));
2337 data_face.op_reinit = [](
auto &phi,
const unsigned batch) {
2338 if (phi.size() == 2)
2340 static_cast<FEFaceEvalType &
>(*phi[0]).
reinit(batch);
2341 static_cast<FEFaceEvalType &
>(*phi[1]).
reinit(batch);
2346 data_face.op_compute = [&](
auto &phi) {
2347 face_operation(
static_cast<FEFaceEvalType &
>(*phi[0]),
2348 static_cast<FEFaceEvalType &
>(*phi[1]));
2351 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2354 data_boundary.dof_numbers = {dof_no};
2355 data_boundary.quad_numbers = {quad_no};
2356 data_boundary.n_components = {n_components};
2357 data_boundary.first_selected_components = {first_selected_component};
2358 data_boundary.batch_type = {1};
2361 .op_create = [&](
const std::pair<unsigned int, unsigned int> &range) {
2363 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
2366 if (!internal::is_fe_nothing<true>(matrix_free,
2370 first_selected_component,
2374 phi.emplace_back(std::make_unique<FEFaceEvalType>(
2375 matrix_free, range,
true, dof_no, quad_no, first_selected_component));
2380 data_boundary.op_reinit = [](
auto &phi,
const unsigned batch) {
2381 if (phi.size() == 1)
2382 static_cast<FEFaceEvalType &
>(*phi[0]).reinit(batch);
2385 if (boundary_operation)
2386 data_boundary.op_compute = [&](
auto &phi) {
2387 boundary_operation(
static_cast<FEFaceEvalType &
>(*phi[0]));
2390 internal::compute_matrix(
2391 matrix_free, constraints_in, data_cell, data_face, data_boundary, matrix);
2398 typename VectorizedArrayType,
2404 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
2406 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2408 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2412 std::unique_ptr<AffineConstraints<typename MatrixType::value_type>>
2413 constraints_for_matrix;
2415 internal::create_new_affine_constraints_if_needed(
2416 matrix, constraints_in, constraints_for_matrix);
2418 const auto batch_operation =
2419 [&matrix_free, &constraints, &
matrix](
2420 auto &
data,
const std::pair<unsigned int, unsigned int> &range) {
2421 if (!
data.op_compute)
2424 auto phi =
data.op_create(range);
2426 const unsigned int n_blocks = phi.size();
2435 n_blocks, VectorizedArrayType::size());
2438 std::array<FullMatrix<typename MatrixType::value_type>,
2439 VectorizedArrayType::size()>>
2440 matrices(n_blocks, n_blocks);
2445 .get_fe(phi[b]->get_active_fe_index());
2447 const auto component_base =
2450 const auto component_in_base =
2451 data.first_selected_components[
b] -
2456 data.dof_numbers[b],
2457 data.quad_numbers[b],
2459 phi[b]->get_active_fe_index(),
2460 phi[b]->get_active_quadrature_index());
2463 shape_info.dofs_per_component_on_cell *
data.n_components[
b];
2465 dof_indices[
b].resize(fe.n_dofs_per_cell());
2467 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2468 dof_indices_mf[b][v].resize(dofs_per_cell[b]);
2470 lexicographic_numbering[
b].insert(
2471 lexicographic_numbering[b].
begin(),
2472 shape_info.lexicographic_numbering.begin() +
2473 component_in_base * shape_info.dofs_per_component_on_cell,
2474 shape_info.lexicographic_numbering.begin() +
2475 (component_in_base +
data.n_components[
b]) *
2476 shape_info.dofs_per_component_on_cell);
2479 for (
unsigned int bj = 0; bj <
n_blocks; ++bj)
2480 for (
unsigned int bi = 0; bi <
n_blocks; ++bi)
2481 std::fill_n(matrices[bi][bj].
begin(),
2482 VectorizedArrayType::size(),
2484 dofs_per_cell[bi], dofs_per_cell[bj]));
2486 for (
auto batch = range.first; batch < range.second; ++batch)
2488 data.op_reinit(phi, batch);
2490 const unsigned int n_filled_lanes =
2495 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
2499 (
data.batch_type[
b] == 0) ?
2500 (batch * VectorizedArrayType::size() + v) :
2501 ((
data.batch_type[
b] == 1) ?
2502 matrix_free.get_face_info(batch).cells_interior[v] :
2503 matrix_free.get_face_info(batch).cells_exterior[v]);
2508 data.dof_numbers[b]);
2512 cell_iterator->get_mg_dof_indices(dof_indices[b]);
2514 cell_iterator->get_dof_indices(dof_indices[b]);
2516 for (
unsigned int j = 0; j < dofs_per_cell[
b]; ++j)
2517 dof_indices_mf[b][v][j] =
2518 dof_indices[b][lexicographic_numbering[b][j]];
2521 for (
unsigned int bj = 0; bj <
n_blocks; ++bj)
2523 for (
unsigned int j = 0; j < dofs_per_cell[bj]; ++j)
2525 for (
unsigned int bi = 0; bi <
n_blocks; ++bi)
2526 for (
unsigned int i = 0; i < dofs_per_cell[bi]; ++i)
2527 phi[bi]->begin_dof_values()[i] =
2528 (bj == bi) ?
static_cast<Number
>(i == j) : 0.0;
2530 data.op_compute(phi);
2532 for (
unsigned int bi = 0; bi <
n_blocks; ++bi)
2533 for (
unsigned int i = 0; i < dofs_per_cell[bi]; ++i)
2534 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
2535 matrices[bi][bj][v](i, j) =
2536 phi[bi]->begin_dof_values()[i][v];
2539 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
2540 for (
unsigned int bi = 0; bi <
n_blocks; ++bi)
2548 matrices[bi][bi][v], dof_indices_mf[bi][v], matrix);
2551 matrices[bi][bj][v],
2552 dof_indices_mf[bi][v],
2553 dof_indices_mf[bj][v],
2559 const auto cell_operation_wrapped =
2560 [&](
const auto &,
auto &,
const auto &,
const auto range) {
2561 batch_operation(data_cell, range);
2564 const auto face_operation_wrapped =
2565 [&](
const auto &,
auto &,
const auto &,
const auto range) {
2566 batch_operation(data_face, range);
2569 const auto boundary_operation_wrapped =
2570 [&](
const auto &,
auto &,
const auto &,
const auto range) {
2571 batch_operation(data_boundary, range);
2574 if (data_face.op_compute || data_boundary.op_compute)
2576 matrix_free.template loop<MatrixType, MatrixType>(
2577 cell_operation_wrapped,
2578 face_operation_wrapped,
2579 boundary_operation_wrapped,
2584 matrix_free.template cell_loop<MatrixType, MatrixType>(
2585 cell_operation_wrapped, matrix, matrix);
2591 template <
typename CLASS,
2597 typename VectorizedArrayType,
2609 VectorizedArrayType> &) const,
2615 VectorizedArrayType> &,
2621 VectorizedArrayType> &)
2628 VectorizedArrayType> &)
2630 const CLASS *owning_class,
2631 const unsigned
int dof_no,
2632 const unsigned
int quad_no,
2633 const unsigned
int first_selected_component)
2640 VectorizedArrayType,
2645 [&](
auto &phi) { (owning_class->*cell_operation)(phi); },
2646 [&](
auto &phi_m,
auto &phi_p) {
2647 (owning_class->*face_operation)(phi_m, phi_p);
2649 [&](
auto &phi) { (owning_class->*boundary_operation)(phi); },
2652 first_selected_component);