40 template <
int dim,
typename AdditionalData>
43 AdditionalData &additional_data);
58 typename VectorizedArrayType,
63 VectorType &diagonal_global,
69 VectorizedArrayType> &)> &local_vmult,
70 const unsigned int dof_no = 0,
71 const unsigned int quad_no = 0,
72 const unsigned int first_selected_component = 0);
77 template <
typename CLASS,
83 typename VectorizedArrayType,
88 VectorType &diagonal_global,
94 VectorizedArrayType> &) const,
95 const CLASS *owning_class,
96 const unsigned
int dof_no = 0,
97 const unsigned
int quad_no = 0,
98 const unsigned
int first_selected_component = 0);
114 typename VectorizedArrayType,
126 VectorizedArrayType> &)> &local_vmult,
127 const unsigned int dof_no = 0,
128 const unsigned int quad_no = 0,
129 const unsigned int first_selected_component = 0);
134 template <
typename CLASS,
140 typename VectorizedArrayType,
152 VectorizedArrayType> &) const,
153 const CLASS *owning_class,
154 const unsigned
int dof_no = 0,
155 const unsigned
int quad_no = 0,
156 const unsigned
int first_selected_component = 0);
208 std::vector<unsigned int> valid_fe_indices;
210 const auto &fe_collection =
212 .get_fe_collection();
214 for (
unsigned int i = 0; i < fe_collection.size(); ++i)
215 if (fe_collection[i].n_dofs_per_cell() > 0)
216 valid_fe_indices.push_back(i);
230 template <
typename VectorTypeOut,
typename VectorTypeIn>
235 const VectorTypeIn &,
236 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
238 const VectorTypeIn &src,
239 const bool zero_dst_vector =
false)
const
241 const auto ebd_cell_operation = [&](
const auto &
matrix_free,
245 const auto category =
matrix_free.get_cell_range_category(range);
253 matrix_free->template cell_loop<VectorTypeOut, VectorTypeIn>(
254 ebd_cell_operation, dst, src, zero_dst_vector);
264 template <
typename VectorTypeOut,
typename VectorTypeIn>
269 const VectorTypeIn &,
270 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
274 const VectorTypeIn &,
275 const std::pair<unsigned int, unsigned int> &)> &face_operation,
279 const VectorTypeIn &,
280 const std::pair<unsigned int, unsigned int> &,
281 const bool)> &boundary_operation,
283 const VectorTypeIn &src,
284 const bool zero_dst_vector =
false)
const
286 const auto ebd_cell_operation = [&](
const auto &
matrix_free,
290 const auto category =
matrix_free.get_cell_range_category(range);
298 const auto ebd_internal_or_boundary_face_operation =
303 const auto category =
matrix_free.get_face_range_category(range);
305 const unsigned int type =
319 matrix_free->template loop<VectorTypeOut, VectorTypeIn>(
321 ebd_internal_or_boundary_face_operation,
322 ebd_internal_or_boundary_face_operation,
345 template <
int dim,
typename AdditionalData>
348 AdditionalData &additional_data)
351 const unsigned int level = additional_data.mg_level;
356 additional_data.cell_vectorization_category.assign(
359 additional_data.cell_vectorization_category.assign(tria.
n_active_cells(),
366 std::sort(bids.begin(), bids.end());
369 unsigned int n_bids = bids.size() + 1;
371 for (
unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
372 i++, offset = offset * n_bids)
376 const auto to_category = [&](
const auto &cell) {
377 unsigned int c_num = 0;
378 for (
unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
380 const auto face = cell->face(i);
381 if (face->at_boundary() && !cell->has_periodic_neighbor(i))
383 factors[i] * (1 + std::distance(bids.begin(),
384 std::find(bids.begin(),
386 face->boundary_id())));
395 if (cell->is_locally_owned())
397 .cell_vectorization_category[cell->active_cell_index()] =
405 if (cell->is_locally_owned_on_level())
406 additional_data.cell_vectorization_category[cell->index()] =
412 additional_data.hold_all_faces_to_owned_cells =
true;
413 additional_data.cell_vectorization_categories_strict =
true;
414 additional_data.mapping_update_flags_faces_by_cells =
415 additional_data.mapping_update_flags_inner_faces |
416 additional_data.mapping_update_flags_boundary_faces;
421 template <
typename Number>
424 std::vector<unsigned int> row_lid_to_gid;
425 std::vector<unsigned int> row;
426 std::vector<unsigned int> col;
427 std::vector<Number> val;
429 std::vector<unsigned int> inverse_lookup_rows;
430 std::vector<std::pair<unsigned int, unsigned int>> inverse_lookup_origins;
438 typename VectorizedArrayType>
439 class ComputeDiagonalHelper
442 static const unsigned int n_lanes = VectorizedArrayType::size();
444 ComputeDiagonalHelper()
446 , dofs_per_component(0)
449 ComputeDiagonalHelper(
const ComputeDiagonalHelper &)
451 , dofs_per_component(0)
460 VectorizedArrayType> &phi)
464 if (dofs_per_component != phi.dofs_per_component)
466 locally_relevant_constraints_hn_map.clear();
467 dofs_per_component = phi.dofs_per_component;
473 reinit(
const unsigned int cell)
475 this->phi->reinit(cell);
478 const auto &dof_info = phi->get_dof_info();
479 const unsigned int n_fe_components = dof_info.start_components.back();
480 const auto &matrix_free = phi->get_matrix_free();
486 const unsigned int first_selected_component =
487 n_fe_components == 1 ? 0 : phi->get_first_selected_component();
489 const unsigned int n_lanes_filled =
497 std::vector<std::tuple<unsigned int, unsigned int, Number>>
498 locally_relevant_constraints, locally_relevant_constraints_tmp;
499 locally_relevant_constraints.reserve(phi->dofs_per_cell);
500 std::vector<unsigned int> constraint_position;
501 std::vector<unsigned char> is_constrained_hn;
503 for (
unsigned int v = 0; v < n_lanes_filled; ++v)
505 const unsigned int *dof_indices;
506 unsigned int index_indicators, next_index_indicators;
508 const unsigned int start =
509 (cell * n_lanes + v) * n_fe_components + first_selected_component;
511 dof_info.dof_indices.data() + dof_info.row_starts[start].first;
512 index_indicators = dof_info.row_starts[start].second;
513 next_index_indicators = dof_info.row_starts[start + 1].second;
517 locally_relevant_constraints.clear();
519 if (n_components == 1 || n_fe_components == 1)
521 unsigned int ind_local = 0;
522 for (; index_indicators != next_index_indicators;
523 ++index_indicators, ++ind_local)
525 const std::pair<unsigned short, unsigned short> indicator =
526 dof_info.constraint_indicator[index_indicators];
528 for (
unsigned int j = 0; j < indicator.first;
530 locally_relevant_constraints.emplace_back(ind_local,
534 dof_indices += indicator.first;
536 const Number *data_val =
538 const Number *end_pool =
541 for (; data_val != end_pool; ++data_val, ++dof_indices)
542 locally_relevant_constraints.emplace_back(ind_local,
549 for (; ind_local < dofs_per_component;
550 ++dof_indices, ++ind_local)
551 locally_relevant_constraints.emplace_back(ind_local,
562 for (
unsigned int comp = 0; comp < n_components; ++comp)
564 unsigned int ind_local = 0;
568 for (; index_indicators != next_index_indicators;
569 ++index_indicators, ++ind_local)
571 const std::pair<unsigned short, unsigned short>
573 dof_info.constraint_indicator[index_indicators];
576 for (
unsigned int j = 0; j < indicator.first;
578 locally_relevant_constraints.emplace_back(
579 comp * dofs_per_component + ind_local,
582 dof_indices += indicator.first;
584 const Number *data_val =
586 const Number *end_pool =
589 for (; data_val != end_pool; ++data_val, ++dof_indices)
590 locally_relevant_constraints.emplace_back(
591 comp * dofs_per_component + ind_local,
599 for (; ind_local < dofs_per_component;
600 ++dof_indices, ++ind_local)
601 locally_relevant_constraints.emplace_back(
602 comp * dofs_per_component + ind_local,
606 if (comp + 1 < n_components)
607 next_index_indicators =
608 dof_info.row_starts[start + comp + 2].second;
615 for (
unsigned int i = 1; i < locally_relevant_constraints.size();
617 Assert(std::get<0>(locally_relevant_constraints[i]) >=
618 std::get<0>(locally_relevant_constraints[i - 1]),
622 if (dof_info.hanging_node_constraint_masks.size() > 0 &&
623 dof_info.hanging_node_constraint_masks_comp.size() > 0 &&
624 dof_info.hanging_node_constraint_masks_comp
625 [phi->get_active_fe_index()][first_selected_component])
628 dof_info.hanging_node_constraint_masks[cell * n_lanes + v];
631 if (mask != ::internal::MatrixFreeFunctions::
636 if (locally_relevant_constraints_hn_map.find(mask) ==
637 locally_relevant_constraints_hn_map.end())
638 fill_constraint_type_into_map(mask);
640 const auto &locally_relevant_constraints_hn =
641 locally_relevant_constraints_hn_map[
mask];
643 locally_relevant_constraints_tmp.clear();
644 if (locally_relevant_constraints_tmp.capacity() <
645 locally_relevant_constraints.size())
646 locally_relevant_constraints_tmp.reserve(
647 locally_relevant_constraints.size() +
648 locally_relevant_constraints_hn.size());
653 constraint_position.assign(phi->dofs_per_cell,
655 for (
auto &a : locally_relevant_constraints)
656 if (constraint_position[
std::get<0>(a)] ==
658 constraint_position[
std::get<0>(a)] =
659 std::distance(locally_relevant_constraints.data(),
661 is_constrained_hn.assign(phi->dofs_per_cell,
false);
662 for (
auto &hn : locally_relevant_constraints_hn)
663 is_constrained_hn[
std::get<0>(hn)] = 1;
666 for (
const auto &a : locally_relevant_constraints)
667 if (is_constrained_hn[
std::get<0>(a)] == 0)
668 locally_relevant_constraints_tmp.push_back(a);
672 for (
const auto &hn : locally_relevant_constraints_hn)
673 if (constraint_position[
std::get<1>(hn)] !=
677 locally_relevant_constraints.size());
678 auto other = locally_relevant_constraints.begin() +
679 constraint_position[std::get<1>(hn)];
682 for (; other != locally_relevant_constraints.end() &&
683 std::get<0>(*other) == std::get<1>(hn);
685 locally_relevant_constraints_tmp.emplace_back(
688 std::get<2>(hn) * std::get<2>(*other));
691 std::swap(locally_relevant_constraints,
692 locally_relevant_constraints_tmp);
697 std::sort(locally_relevant_constraints.begin(),
698 locally_relevant_constraints.end(),
699 [](
const auto &a,
const auto &b) {
700 if (std::get<1>(a) < std::get<1>(b))
702 return (std::get<1>(a) == std::get<1>(b)) &&
703 (std::get<0>(a) < std::get<0>(b));
707 auto &c_pool = c_pools[v];
709 c_pool.row_lid_to_gid.clear();
711 c_pool.row.push_back(0);
715 if (locally_relevant_constraints.size() > 0)
716 c_pool.row_lid_to_gid.emplace_back(
717 std::get<1>(locally_relevant_constraints.front()));
718 for (
const auto &j : locally_relevant_constraints)
720 if (c_pool.row_lid_to_gid.back() != std::get<1>(j))
722 c_pool.row_lid_to_gid.push_back(std::get<1>(j));
723 c_pool.row.push_back(c_pool.val.size());
726 c_pool.col.emplace_back(std::get<0>(j));
727 c_pool.val.emplace_back(std::get<2>(j));
730 if (c_pool.val.size() > 0)
731 c_pool.row.push_back(c_pool.val.size());
733 c_pool.inverse_lookup_rows.clear();
734 c_pool.inverse_lookup_rows.resize(1 + phi->dofs_per_cell);
735 for (
const unsigned int i : c_pool.col)
736 c_pool.inverse_lookup_rows[1 + i]++;
738 std::partial_sum(c_pool.inverse_lookup_rows.begin(),
739 c_pool.inverse_lookup_rows.end(),
740 c_pool.inverse_lookup_rows.begin());
744 c_pool.inverse_lookup_origins.resize(c_pool.col.size());
745 std::vector<unsigned int> inverse_lookup_count(
747 for (
unsigned int row = 0; row < c_pool.row.size() - 1; ++row)
748 for (
unsigned int col = c_pool.row[row];
749 col < c_pool.row[row + 1];
752 const unsigned int index = c_pool.col[col];
753 c_pool.inverse_lookup_origins
754 [c_pool.inverse_lookup_rows[
index] +
755 inverse_lookup_count[
index]] = std::make_pair(row, col);
756 ++inverse_lookup_count[
index];
785 for (
unsigned int v = 0; v < n_lanes_filled; ++v)
786 diagonals_local_constrained[v].assign(
787 c_pools[v].row_lid_to_gid.size() *
788 (n_fe_components == 1 ? n_components : 1),
793 fill_constraint_type_into_map(
794 const ::internal::MatrixFreeFunctions::compressed_constraint_kind
797 auto &constraints_hn = locally_relevant_constraints_hn_map[
mask];
802 const unsigned int degree =
803 phi->get_shape_info().data.front().fe_degree;
808 values_dofs.resize(dofs_per_component);
811 VectorizedArrayType::size()>
813 constraint_mask.fill(::internal::MatrixFreeFunctions::
815 constraint_mask[0] =
mask;
817 for (
unsigned int i = 0; i < dofs_per_component; ++i)
819 for (
unsigned int j = 0; j < dofs_per_component; ++j)
820 values_dofs[j] = VectorizedArrayType();
821 values_dofs[i] = Number(1);
826 VectorizedArrayType>::apply(1,
828 phi->get_shape_info(),
833 const Number tolerance =
835 std::numeric_limits<Number>::epsilon() * 16);
836 for (
unsigned int j = 0; j < dofs_per_component; ++j)
837 if (
std::abs(values_dofs[j][0]) > tolerance &&
839 std::abs(values_dofs[j][0] - Number(1)) > tolerance))
840 constraints_hn.emplace_back(j, i, values_dofs[j][0]);
844 const unsigned int n_hn_constraints = constraints_hn.size();
845 constraints_hn.resize(n_hn_constraints * n_components);
847 for (
unsigned int c = 1; c < n_components; ++c)
848 for (
unsigned int i = 0; i < n_hn_constraints; ++i)
849 constraints_hn[c * n_hn_constraints + i] = std::make_tuple(
850 std::get<0>(constraints_hn[i]) + c * dofs_per_component,
851 std::get<1>(constraints_hn[i]) + c * dofs_per_component,
852 std::get<2>(constraints_hn[i]));
856 prepare_basis_vector(
const unsigned int i)
863 VectorizedArrayType *dof_values = phi->begin_dof_values();
864 for (
const unsigned int j : phi->dof_indices())
865 dof_values[j] = VectorizedArrayType();
866 dof_values[i] = Number(1);
875 const unsigned int n_fe_components =
876 phi->get_dof_info().start_components.back();
877 const unsigned int comp =
878 n_fe_components == 1 ? i / dofs_per_component : 0;
879 const unsigned int i_comp =
880 n_fe_components == 1 ? (i % dofs_per_component) : i;
884 for (
unsigned int v = 0;
885 v < phi->get_matrix_free().n_active_entries_per_cell_batch(
886 phi->get_current_cell_index());
889 const auto &c_pool = c_pools[v];
891 for (
unsigned int jj = c_pool.inverse_lookup_rows[i_comp];
892 jj < c_pool.inverse_lookup_rows[i_comp + 1];
895 const unsigned int j = c_pool.inverse_lookup_origins[jj].first;
898 for (
unsigned int k = c_pool.row[j]; k < c_pool.row[j + 1]; ++k)
899 temp += c_pool.val[k] *
900 phi->begin_dof_values()[comp * dofs_per_component +
904 diagonals_local_constrained
905 [v][j + comp * c_pools[v].row_lid_to_gid.size()] +=
906 temp * c_pool.val[c_pool.inverse_lookup_origins[jj].second];
911 template <
typename VectorType>
913 distribute_local_to_global(
914 std::array<VectorType *, n_components> &diagonal_global)
917 const unsigned int n_fe_components =
918 phi->get_dof_info().start_components.back();
920 for (
unsigned int v = 0;
921 v < phi->get_matrix_free().n_active_entries_per_cell_batch(
922 phi->get_current_cell_index());
927 for (
unsigned int j = 0; j < c_pools[v].row.size() - 1; ++j)
928 for (
unsigned int comp = 0;
929 comp < (n_fe_components == 1 ?
930 static_cast<unsigned int>(n_components) :
934 *diagonal_global[n_fe_components == 1 ? comp : 0],
935 c_pools[v].row_lid_to_gid[j],
936 diagonals_local_constrained
937 [v][j + comp * c_pools[v].row_lid_to_gid.size()]);
946 VectorizedArrayType> *phi;
948 unsigned int dofs_per_component;
952 std::array<internal::LocalCSR<Number>, n_lanes> c_pools;
957 std::array<std::vector<Number>, n_lanes> diagonals_local_constrained;
961 std::vector<std::tuple<unsigned int, unsigned int, Number>>>
962 locally_relevant_constraints_hn_map;
975 typename VectorizedArrayType,
980 VectorType &diagonal_global,
986 VectorizedArrayType> &)> &local_vmult,
987 const unsigned int dof_no,
988 const unsigned int quad_no,
989 const unsigned int first_selected_component)
993 std::array<typename ::internal::BlockVectorSelector<
997 diagonal_global_components;
999 for (
unsigned int d = 0;
d < n_components; ++
d)
1000 diagonal_global_components[d] = ::internal::
1002 get_vector_component(diagonal_global,
d + first_selected_component);
1004 const auto &dof_info = matrix_free.
get_dof_info(dof_no);
1006 if (dof_info.start_components.back() == 1)
1007 for (
unsigned int comp = 0; comp < n_components; ++comp)
1009 Assert(diagonal_global_components[comp] !=
nullptr,
1010 ExcMessage(
"The finite element underlying this FEEvaluation "
1011 "object is scalar, but you requested " +
1012 std::to_string(n_components) +
1013 " components via the template argument in "
1014 "FEEvaluation. In that case, you must pass an "
1015 "std::vector<VectorType> or a BlockVector to " +
1016 "read_dof_values and distribute_local_to_global."));
1018 *diagonal_global_components[comp], matrix_free, dof_info);
1023 *diagonal_global_components[0], matrix_free, dof_info);
1026 using Helper = internal::ComputeDiagonalHelper<dim,
1031 VectorizedArrayType>;
1034 matrix_free.template cell_loop<VectorType, int>(
1038 const std::pair<unsigned int, unsigned int> &range)
mutable {
1039 Helper &helper = scratch_data.
get();
1046 VectorizedArrayType>
1047 phi(matrix_free, range, dof_no, quad_no, first_selected_component);
1048 helper.initialize(phi);
1050 for (
unsigned int cell = range.first; cell < range.second; ++cell)
1052 helper.reinit(cell);
1054 for (
unsigned int i = 0; i < phi.dofs_per_cell; ++i)
1056 helper.prepare_basis_vector(i);
1061 helper.distribute_local_to_global(diagonal_global_components);
1069 template <
typename CLASS,
1075 typename VectorizedArrayType,
1076 typename VectorType>
1080 VectorType &diagonal_global,
1086 VectorizedArrayType> &) const,
1087 const CLASS *owning_class,
1088 const unsigned
int dof_no,
1089 const unsigned
int quad_no,
1090 const unsigned
int first_selected_component)
1097 VectorizedArrayType,
1101 [&](
auto &feeval) { (owning_class->*cell_operation)(feeval); },
1104 first_selected_component);
1113 template <
typename MatrixType,
1115 std::enable_if_t<std::is_same_v<
1116 typename std::remove_const<
typename std::remove_reference<
1117 typename MatrixType::value_type>::type>::type,
1118 typename std::remove_const<
typename std::remove_reference<
1119 Number>::type>::type>> * =
nullptr>
1121 create_new_affine_constraints_if_needed(
1134 template <
typename MatrixType,
1136 std::enable_if_t<!std::is_same_v<
1137 typename std::remove_const<
typename std::remove_reference<
1138 typename MatrixType::value_type>::type>::type,
1139 typename std::remove_const<
typename std::remove_reference<
1140 Number>::type>::type>> * =
nullptr>
1142 create_new_affine_constraints_if_needed(
1149 std::make_unique<AffineConstraints<typename MatrixType::value_type>>();
1150 new_constraints->
copy_from(constraints);
1152 return *new_constraints;
1161 typename VectorizedArrayType,
1162 typename MatrixType>
1173 VectorizedArrayType> &)> &local_vmult,
1174 const unsigned int dof_no,
1175 const unsigned int quad_no,
1176 const unsigned int first_selected_component)
1178 std::unique_ptr<AffineConstraints<typename MatrixType::value_type>>
1179 constraints_for_matrix;
1181 internal::create_new_affine_constraints_if_needed(matrix,
1183 constraints_for_matrix);
1185 matrix_free.template cell_loop<MatrixType, MatrixType>(
1186 [&](
const auto &,
auto &dst,
const auto &,
const auto range) {
1192 VectorizedArrayType>
1194 matrix_free, range, dof_no, quad_no, first_selected_component);
1196 const unsigned int dofs_per_cell = integrator.dofs_per_cell;
1198 std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
1199 std::vector<types::global_dof_index> dof_indices_mf(dofs_per_cell);
1201 std::array<FullMatrix<typename MatrixType::value_type>,
1202 VectorizedArrayType::size()>
1205 std::fill_n(matrices.begin(),
1206 VectorizedArrayType::size(),
1210 const auto lexicographic_numbering =
1214 first_selected_component,
1215 integrator.get_active_fe_index(),
1216 integrator.get_active_quadrature_index())
1217 .lexicographic_numbering;
1219 for (
auto cell = range.first; cell < range.second; ++cell)
1221 integrator.reinit(cell);
1223 const unsigned int n_filled_lanes =
1226 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
1229 for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1231 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1232 integrator.begin_dof_values()[i] =
1233 static_cast<Number
>(i == j);
1235 local_vmult(integrator);
1237 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1238 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
1239 matrices[v](i, j) = integrator.begin_dof_values()[i][v];
1242 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
1248 cell_v->get_mg_dof_indices(dof_indices);
1250 cell_v->get_dof_indices(dof_indices);
1252 for (
unsigned int j = 0; j < dof_indices.size(); ++j)
1253 dof_indices_mf[j] = dof_indices[lexicographic_numbering[j]];
1267 template <
typename CLASS,
1273 typename VectorizedArrayType,
1274 typename MatrixType>
1285 VectorizedArrayType> &) const,
1286 const CLASS *owning_class,
1287 const unsigned
int dof_no,
1288 const unsigned
int quad_no,
1289 const unsigned
int first_selected_component)
1296 VectorizedArrayType,
1301 [&](
auto &feeval) { (owning_class->*cell_operation)(feeval); },
1304 first_selected_component);