39 template <
int dim,
typename Number,
bool is_face_>
52 std::function<std::vector<std::unique_ptr<FEEvalType>>(
53 const std::pair<unsigned int, unsigned int> &)>
55 std::function<
void(std::vector<std::unique_ptr<FEEvalType>> &,
58 std::function<
void(std::vector<std::unique_ptr<FEEvalType>> &)>
68 template <
int dim,
typename AdditionalData>
71 AdditionalData &additional_data);
91 typename VectorizedArrayType,
102 VectorizedArrayType> &)>
104 const unsigned int dof_no = 0,
105 const unsigned int quad_no = 0,
106 const unsigned int first_selected_component = 0,
107 const unsigned int first_vector_component = 0);
126 const unsigned int dof_no = 0,
127 const unsigned int quad_no = 0,
128 const unsigned int first_selected_component = 0,
129 const unsigned int first_vector_component = 0);
134 template <
typename CLASS,
140 typename VectorizedArrayType,
151 VectorizedArrayType> &)
const,
155 const unsigned
int first_selected_component = 0,
156 const unsigned
int first_vector_component = 0);
175 typename VectorizedArrayType,
186 VectorizedArrayType> &)>
193 VectorizedArrayType> &,
199 VectorizedArrayType> &)>
206 VectorizedArrayType> &)>
208 const unsigned int dof_no = 0,
209 const unsigned int quad_no = 0,
210 const unsigned int first_selected_component = 0,
211 const unsigned int first_vector_component = 0);
218 template <
typename CLASS,
224 typename VectorizedArrayType,
235 VectorizedArrayType> &)
const,
241 VectorizedArrayType> &,
247 VectorizedArrayType> &)
254 VectorizedArrayType> &)
259 const unsigned
int first_selected_component = 0,
260 const unsigned
int first_vector_component = 0);
277 typename VectorizedArrayType,
289 VectorizedArrayType> &)>
291 const unsigned int dof_no = 0,
292 const unsigned int quad_no = 0,
293 const unsigned int first_selected_component = 0);
300 template <
typename CLASS,
306 typename VectorizedArrayType,
318 VectorizedArrayType> &)
const,
322 const unsigned
int first_selected_component = 0);
333 typename VectorizedArrayType,
354 typename VectorizedArrayType,
385 typename VectorizedArrayType,
397 VectorizedArrayType> &)>
404 VectorizedArrayType> &,
410 VectorizedArrayType> &)>
417 VectorizedArrayType> &)>
419 const unsigned int dof_no = 0,
420 const unsigned int quad_no = 0,
421 const unsigned int first_selected_component = 0);
428 template <
typename CLASS,
434 typename VectorizedArrayType,
446 VectorizedArrayType> &)
const,
452 VectorizedArrayType> &,
458 VectorizedArrayType> &)
465 VectorizedArrayType> &)
470 const unsigned
int first_selected_component = 0);
524 const auto &fe_collection =
525 matrix_free.get_dof_handler(additional_data.dof_index)
526 .get_fe_collection();
528 for (
unsigned int i = 0; i < fe_collection.size(); ++i)
529 if (fe_collection[i].n_dofs_per_cell() > 0)
544 template <
typename VectorTypeOut,
typename VectorTypeIn>
578 template <
typename VectorTypeOut,
typename VectorTypeIn>
594 const std::pair<unsigned int, unsigned int> &,
619 const unsigned int type =
659 template <
int dim,
typename AdditionalData>
662 AdditionalData &additional_data)
665 const unsigned int level = additional_data.mg_level;
670 additional_data.cell_vectorization_category.assign(
671 std::distance(tria.begin(
level), tria.end(
level)), 0);
673 additional_data.cell_vectorization_category.assign(tria.n_active_cells(),
679 auto bids = tria.get_boundary_ids();
680 std::sort(
bids.begin(),
bids.end());
686 i++, offset = offset *
n_bids)
691 unsigned int c_num = 0;
694 const auto face = cell->face(i);
695 if (face->at_boundary() && !cell->has_periodic_neighbor(i))
697 factors[i] * (1 + std::distance(
bids.begin(),
698 std::find(
bids.begin(),
700 face->boundary_id())));
707 for (
auto cell = tria.begin_active(); cell != tria.end(); ++cell)
709 if (cell->is_locally_owned())
711 .cell_vectorization_category[cell->active_cell_index()] =
717 for (
auto cell = tria.begin(
level); cell != tria.end(
level); ++cell)
719 if (cell->is_locally_owned_on_level())
720 additional_data.cell_vectorization_category[cell->index()] =
726 additional_data.hold_all_faces_to_owned_cells =
true;
727 additional_data.cell_vectorization_categories_strict =
true;
728 additional_data.mapping_update_flags_faces_by_cells =
729 additional_data.mapping_update_flags_inner_faces |
730 additional_data.mapping_update_flags_boundary_faces;
735 template <
typename Number>
739 std::vector<unsigned int> row;
740 std::vector<unsigned int> col;
741 std::vector<Number> val;
747 template <
int dim,
typename VectorizedArrayType,
bool is_face>
754 using Number =
typename VectorizedArrayType::value_type;
755 static const unsigned int n_lanes = VectorizedArrayType::size();
760 , dofs_per_component(0)
767 , dofs_per_component(0)
775 const unsigned int n_components)
779 if (dofs_per_component !=
780 phi.get_shape_info().dofs_per_component_on_cell)
784 phi.get_shape_info().dofs_per_component_on_cell;
786 this->n_components = n_components;
787 this->dofs_per_cell = n_components * dofs_per_component;
789 this->matrix_free = &matrix_free;
793 reinit(
const unsigned int cell)
796 const auto &matrix_free = *this->matrix_free;
797 const auto &dof_info =
phi->get_dof_info();
798 const unsigned int n_fe_components = dof_info.start_components.back();
804 const unsigned int first_selected_component =
805 n_fe_components == 1 ? 0 :
phi->get_first_selected_component();
808 is_face ? matrix_free.n_active_entries_per_face_batch(cell) :
809 matrix_free.n_active_entries_per_cell_batch(cell);
816 std::vector<std::tuple<unsigned int, unsigned int, Number>>
822 const std::array<unsigned int, n_lanes> &cells =
823 this->phi->get_cell_ids();
830 const unsigned int *dof_indices;
833 const unsigned int start =
834 cells[v] * n_fe_components + first_selected_component;
836 dof_info.dof_indices.data() + dof_info.row_starts[start].first;
844 if (n_components == 1 || n_fe_components == 1)
850 const std::pair<unsigned short, unsigned short>
indicator =
862 matrix_free.constraint_pool_begin(
indicator.second);
864 matrix_free.constraint_pool_end(
indicator.second);
887 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
896 const std::pair<unsigned short, unsigned short>
910 matrix_free.constraint_pool_begin(
indicator.second);
912 matrix_free.constraint_pool_end(
indicator.second);
931 if (
comp + 1 < n_components)
933 dof_info.row_starts[start +
comp + 2].second;
947 if (dof_info.hanging_node_constraint_masks.size() > 0 &&
948 dof_info.hanging_node_constraint_masks_comp.size() > 0 &&
949 dof_info.hanging_node_constraint_masks_comp
950 [
phi->get_active_fe_index()][first_selected_component])
953 dof_info.hanging_node_constraint_masks[cells[v]];
956 if (mask != ::internal::MatrixFreeFunctions::
1008 std::get<0>(*
other) == std::get<1>(
hn);
1012 std::get<1>(*
other),
1013 std::get<2>(
hn) * std::get<2>(*
other));
1024 [](
const auto &a,
const auto &b) {
1025 if (std::get<1>(a) < std::get<1>(b))
1027 return (std::get<1>(a) == std::get<1>(b)) &&
1028 (std::get<0>(a) < std::get<0>(b));
1041 c_pool.row_lid_to_gid.emplace_back(
1045 if (
c_pool.row_lid_to_gid.back() != std::get<1>(
j))
1047 c_pool.row_lid_to_gid.push_back(std::get<1>(
j));
1051 c_pool.col.emplace_back(std::get<0>(
j));
1052 c_pool.val.emplace_back(std::get<2>(
j));
1055 if (
c_pool.val.size() > 0)
1059 c_pool.inverse_lookup_rows.resize(1 + dofs_per_cell);
1060 for (
const unsigned int i :
c_pool.col)
1063 std::partial_sum(
c_pool.inverse_lookup_rows.begin(),
1064 c_pool.inverse_lookup_rows.end(),
1065 c_pool.inverse_lookup_rows.begin());
1069 c_pool.inverse_lookup_origins.resize(
c_pool.col.size());
1071 for (
unsigned int row = 0; row <
c_pool.row.size() - 1; ++row)
1072 for (
unsigned int col =
c_pool.row[row];
1073 col <
c_pool.row[row + 1];
1077 c_pool.inverse_lookup_origins
1112 (n_fe_components == 1 ? n_components : 1),
1122 for (
unsigned int i = 0; i <
c_pool.row.size() - 1; ++i)
1146 const ::internal::MatrixFreeFunctions::compressed_constraint_kind
1154 const unsigned int degree =
1155 phi->get_shape_info().data.front().fe_degree;
1160 values_dofs.resize(dofs_per_component);
1163 VectorizedArrayType::size()>
1165 constraint_mask.fill(::internal::MatrixFreeFunctions::
1167 constraint_mask[0] =
mask;
1169 for (
unsigned int i = 0; i < dofs_per_component; ++i)
1171 for (
unsigned int j = 0;
j < dofs_per_component; ++
j)
1172 values_dofs[
j] = VectorizedArrayType();
1173 values_dofs[i] = Number(1);
1178 VectorizedArrayType>::apply(1,
1180 phi->get_shape_info(),
1183 values_dofs.data());
1185 const Number tolerance =
1187 std::numeric_limits<Number>::epsilon() * 16);
1188 for (
unsigned int j = 0;
j < dofs_per_component; ++
j)
1189 if (
std::abs(values_dofs[
j][0]) > tolerance &&
1191 std::abs(values_dofs[
j][0] - Number(1)) > tolerance))
1199 for (
unsigned int c = 1; c < n_components; ++c)
1216 for (
unsigned int j = 0;
j < dofs_per_cell; ++
j)
1225 for (
unsigned int j = 0;
j < dofs_per_cell; ++
j)
1235 const unsigned int n_fe_components =
1236 phi->get_dof_info().start_components.back();
1237 const unsigned int comp =
1238 n_fe_components == 1 ? i / dofs_per_component : 0;
1239 const unsigned int i_comp =
1240 n_fe_components == 1 ? (i % dofs_per_component) : i;
1252 const unsigned int j =
c_pool.inverse_lookup_origins[
jj].first;
1257 phi->begin_dof_values()[
comp * dofs_per_component +
1268 template <
typename VectorType>
1270 distribute_local_to_global(std::vector<VectorType *> &
diagonal_global)
1273 const unsigned int n_fe_components =
1274 phi->get_dof_info().start_components.back();
1276 if (n_fe_components == 1)
1283 for (
unsigned int j = 0;
j <
c_pools[v].row.size() - 1; ++
j)
1284 for (
unsigned int comp = 0;
1285 comp < (n_fe_components == 1 ?
1286 static_cast<unsigned int>(n_components) :
1306 unsigned int dofs_per_component;
1307 unsigned int dofs_per_cell;
1308 unsigned int n_components;
1314 std::array<internal::LocalCSR<Number>, n_lanes>
c_pools;
1323 std::vector<std::tuple<unsigned int, unsigned int, Number>>>
1332 template <
bool is_face,
1335 typename VectorizedArrayType>
1339 const std::pair<unsigned int, unsigned int> &
range,
1340 const unsigned int dof_no,
1341 const unsigned int quad_no,
1342 const unsigned int first_selected_component,
1343 const unsigned int fe_degree,
1344 const unsigned int n_q_points_1d,
1345 const bool is_interior_face =
true)
1347 const unsigned int static_n_q_points =
1351 unsigned int active_fe_index = 0;
1353 active_fe_index = matrix_free.get_cell_active_fe_index(
range);
1354 else if (is_interior_face)
1355 active_fe_index = matrix_free.get_face_range_category(
range).first;
1357 active_fe_index = matrix_free.get_face_range_category(
range).second;
1360 extract_initialization_data<is_face, dim, Number, VectorizedArrayType>(
1363 first_selected_component,
1371 return init_data.shape_info->dofs_per_component_on_cell == 0;
1395 operator()(
const unsigned int cell,
1402 FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number>
1406 constexpr int dofs_per_cell =
decltype(fe_eval)::tensor_dofs_per_cell;
1407 typename decltype(fe_eval)::value_type
1408 diagonal[dofs_per_cell / n_components] = {};
1409 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1411 const auto c = i % n_components;
1413 Kokkos::parallel_for(
1415 dofs_per_cell / n_components),
1417 typename decltype(fe_eval)::value_type val = {};
1419 if constexpr (n_components == 1)
1421 val = (i ==
j) ? 1 : 0;
1425 val[c] = (i / n_components ==
j) ? 1 : 0;
1428 fe_eval.submit_dof_value(val,
j);
1433 Portable::internal::
1434 resolve_hanging_nodes<dim, fe_degree, false, Number>(
1437 gpu_data->constraint_mask(cell * n_components + c),
1438 Kokkos::subview(shared_data->
values, Kokkos::ALL, c));
1444 Portable::internal::
1445 resolve_hanging_nodes<dim, fe_degree, true, Number>(
1448 gpu_data->constraint_mask(cell * n_components + c),
1449 Kokkos::subview(shared_data->
values, Kokkos::ALL, c));
1451 Kokkos::single(Kokkos::PerTeam(shared_data->
team_member), [&] {
1452 if constexpr (n_components == 1)
1453 diagonal[i] = fe_eval.get_dof_value(i);
1455 diagonal[i / n_components][i % n_components] =
1456 fe_eval.get_dof_value(i / n_components)[i % n_components];
1462 Kokkos::single(Kokkos::PerTeam(shared_data->
team_member), [&] {
1463 for (unsigned int i = 0; i < dofs_per_cell / n_components; ++i)
1464 fe_eval.submit_dof_value(diagonal[i], i);
1473 Kokkos::parallel_for(
1474 Kokkos::TeamThreadRange(shared_data->
team_member, dofs_per_cell),
1476 dst[gpu_data->local_to_global(cell, i)] +=
1477 shared_data->values(i % (dofs_per_cell / n_components),
1478 i / (dofs_per_cell / n_components));
1483 Kokkos::parallel_for(
1484 Kokkos::TeamThreadRange(shared_data->
team_member, dofs_per_cell),
1487 &dst[gpu_data->local_to_global(cell, i)],
1488 shared_data->values(i % (dofs_per_cell / n_components),
1489 i / (dofs_per_cell / n_components)));
1496 static constexpr unsigned int n_q_points =
1520 const unsigned int dof_no,
1521 const unsigned int quad_no,
1522 const unsigned int first_selected_component,
1523 const unsigned int first_vector_component)
1551 typename VectorizedArrayType,
1562 VectorizedArrayType> &)>
1564 const unsigned int dof_no,
1565 const unsigned int quad_no,
1566 const unsigned int first_selected_component,
1567 const unsigned int first_vector_component)
1574 VectorizedArrayType,
1582 first_selected_component,
1583 first_vector_component);
1586 template <
typename CLASS,
1592 typename VectorizedArrayType,
1603 VectorizedArrayType> &)
const,
1607 const unsigned
int first_selected_component,
1608 const unsigned
int first_vector_component)
1615 VectorizedArrayType,
1622 first_selected_component,
1623 first_vector_component);
1631 typename VectorizedArrayType,
1642 VectorizedArrayType> &)>
1649 VectorizedArrayType> &,
1655 VectorizedArrayType> &)>
1662 VectorizedArrayType> &)>
1664 const unsigned int dof_no,
1665 const unsigned int quad_no,
1666 const unsigned int first_selected_component,
1667 const unsigned int first_vector_component)
1669 std::vector<typename ::internal::BlockVectorSelector<
1674 for (
unsigned int d = 0;
d < n_components; ++
d)
1679 const auto &dof_info = matrix_free.get_dof_info(
dof_no);
1681 if (dof_info.start_components.back() == 1)
1682 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
1685 ExcMessage(
"The finite element underlying this FEEvaluation "
1686 "object is scalar, but you requested " +
1687 std::to_string(n_components) +
1688 " components via the template argument in "
1689 "FEEvaluation. In that case, you must pass an "
1690 "std::vector<VectorType> or a BlockVector to " +
1691 "read_dof_values and distribute_local_to_global."));
1706 VectorizedArrayType>;
1713 VectorizedArrayType>;
1715 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
1720 data_cell.n_components = {n_components};
1721 data_cell.first_selected_components = {first_selected_component};
1725 [&](
const std::pair<unsigned int, unsigned int> &
range) {
1727 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, false>>>
1730 if (!internal::is_fe_nothing<false>(matrix_free,
1734 first_selected_component,
1737 phi.emplace_back(std::make_unique<FEEvalType>(
1738 matrix_free,
range,
dof_no, quad_no, first_selected_component));
1744 if (
phi.size() == 1)
1745 static_cast<FEEvalType &
>(*
phi[0]).reinit(
batch);
1753 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1757 data_face.quad_numbers = {quad_no, quad_no};
1758 data_face.n_components = {n_components, n_components};
1759 data_face.first_selected_components = {first_selected_component,
1760 first_selected_component};
1764 [&](
const std::pair<unsigned int, unsigned int> &
range) {
1766 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
1769 if (!internal::is_fe_nothing<true>(matrix_free,
1773 first_selected_component,
1777 !internal::is_fe_nothing<true>(matrix_free,
1781 first_selected_component,
1787 std::make_unique<FEFaceEvalType>(matrix_free,
1792 first_selected_component));
1794 std::make_unique<FEFaceEvalType>(matrix_free,
1799 first_selected_component));
1806 if (
phi.size() == 2)
1819 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1825 data_boundary.first_selected_components = {first_selected_component};
1829 .op_create = [&](
const std::pair<unsigned int, unsigned int> &
range) {
1831 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
1834 if (!internal::is_fe_nothing<true>(matrix_free,
1838 first_selected_component,
1842 phi.emplace_back(std::make_unique<FEFaceEvalType>(
1843 matrix_free,
range,
true,
dof_no, quad_no, first_selected_component));
1849 if (
phi.size() == 1)
1858 internal::compute_diagonal(matrix_free,
1870 typename VectorizedArrayType,
1876 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
1878 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1880 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1890 internal::ComputeDiagonalHelper<dim, VectorizedArrayType, false>;
1893 internal::ComputeDiagonalHelper<dim, VectorizedArrayType, true>;
1903 const std::pair<unsigned int, unsigned int> &
range) {
1904 if (!
data.op_compute)
1915 helpers[b].initialize(*
phi[b], matrix_free,
data.n_components[b]);
1929 return helper.has_simple_constraints();
1936 for (
unsigned int i = 0;
1937 i <
phi[
b]->get_shape_info().dofs_per_component_on_cell *
1938 data.n_components[b];
1951 helpers[
b].distribute_local_to_global(
1958 [&](
const auto &,
auto &,
const auto &,
const auto range) {
1963 [&](
const auto &,
auto &,
const auto &,
const auto range) {
1968 [&](
const auto &,
auto &,
const auto &,
const auto range) {
1987 template <
typename CLASS,
1993 typename VectorizedArrayType,
2004 VectorizedArrayType> &)
const,
2010 VectorizedArrayType> &,
2016 VectorizedArrayType> &)
2023 VectorizedArrayType> &)
2028 const unsigned
int first_selected_component,
2029 const unsigned
int first_vector_component)
2036 VectorizedArrayType,
2047 first_selected_component,
2048 first_vector_component);
2060 std::enable_if_t<std::is_same_v<
2061 std::remove_const_t<
2062 std::remove_reference_t<typename MatrixType::value_type>>,
2063 std::remove_const_t<std::remove_reference_t<Number>>>> * =
nullptr>
2081 std::enable_if_t<!std::is_same_v<
2082 std::remove_const_t<
2083 std::remove_reference_t<typename MatrixType::value_type>>,
2084 std::remove_const_t<std::remove_reference_t<Number>>>> * =
nullptr>
2093 std::make_unique<AffineConstraints<typename MatrixType::value_type>>();
2105 typename VectorizedArrayType,
2117 VectorizedArrayType> &)>
2119 const unsigned int dof_no,
2120 const unsigned int quad_no,
2121 const unsigned int first_selected_component)
2128 VectorizedArrayType,
2137 first_selected_component);
2140 template <
typename CLASS,
2146 typename VectorizedArrayType,
2158 VectorizedArrayType> &)
const,
2162 const unsigned
int first_selected_component)
2169 VectorizedArrayType,
2177 first_selected_component);
2185 typename VectorizedArrayType,
2197 VectorizedArrayType> &)>
2204 VectorizedArrayType> &,
2210 VectorizedArrayType> &)>
2217 VectorizedArrayType> &)>
2219 const unsigned int dof_no,
2220 const unsigned int quad_no,
2221 const unsigned int first_selected_component)
2228 VectorizedArrayType>;
2235 VectorizedArrayType>;
2237 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
2242 data_cell.n_components = {n_components};
2243 data_cell.first_selected_components = {first_selected_component};
2247 [&](
const std::pair<unsigned int, unsigned int> &
range) {
2249 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, false>>>
2252 if (!internal::is_fe_nothing<false>(matrix_free,
2256 first_selected_component,
2259 phi.emplace_back(std::make_unique<FEEvalType>(
2260 matrix_free,
range,
dof_no, quad_no, first_selected_component));
2266 if (
phi.size() == 1)
2267 static_cast<FEEvalType &
>(*
phi[0]).reinit(
batch);
2275 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2279 data_face.quad_numbers = {quad_no, quad_no};
2280 data_face.n_components = {n_components, n_components};
2281 data_face.first_selected_components = {first_selected_component,
2282 first_selected_component};
2286 [&](
const std::pair<unsigned int, unsigned int> &
range) {
2288 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
2291 if (!internal::is_fe_nothing<true>(matrix_free,
2295 first_selected_component,
2299 !internal::is_fe_nothing<true>(matrix_free,
2303 first_selected_component,
2309 std::make_unique<FEFaceEvalType>(matrix_free,
2314 first_selected_component));
2316 std::make_unique<FEFaceEvalType>(matrix_free,
2321 first_selected_component));
2328 if (
phi.size() == 2)
2341 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2347 data_boundary.first_selected_components = {first_selected_component};
2351 .op_create = [&](
const std::pair<unsigned int, unsigned int> &
range) {
2353 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
2356 if (!internal::is_fe_nothing<true>(matrix_free,
2360 first_selected_component,
2364 phi.emplace_back(std::make_unique<FEFaceEvalType>(
2365 matrix_free,
range,
true,
dof_no, quad_no, first_selected_component));
2371 if (
phi.size() == 1)
2380 internal::compute_matrix(
2388 typename VectorizedArrayType,
2394 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
2396 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2398 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2402 std::unique_ptr<AffineConstraints<typename MatrixType::value_type>>
2405 internal::create_new_affine_constraints_if_needed(
2409 [&matrix_free, &constraints, &
matrix](
2410 auto &
data,
const std::pair<unsigned int, unsigned int> &
range) {
2411 if (!
data.op_compute)
2425 n_blocks, VectorizedArrayType::size());
2428 std::array<FullMatrix<typename MatrixType::value_type>,
2429 VectorizedArrayType::size()>>
2430 matrices(n_blocks, n_blocks);
2434 const auto &fe = matrix_free.get_dof_handler(
data.dof_numbers[b])
2435 .get_fe(
phi[b]->get_active_fe_index());
2438 matrix_free.get_dof_info(
data.dof_numbers[b])
2439 .component_to_base_index[
data.first_selected_components[
b]];
2441 data.first_selected_components[
b] -
2442 matrix_free.get_dof_info(
data.dof_numbers[b])
2445 const auto &shape_info = matrix_free.get_shape_info(
2446 data.dof_numbers[b],
2447 data.quad_numbers[b],
2449 phi[b]->get_active_fe_index(),
2450 phi[b]->get_active_quadrature_index());
2453 shape_info.dofs_per_component_on_cell *
data.n_components[
b];
2455 dof_indices[
b].resize(fe.n_dofs_per_cell());
2457 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2460 lexicographic_numbering[
b].insert(
2461 lexicographic_numbering[b].
begin(),
2462 shape_info.lexicographic_numbering.begin() +
2464 shape_info.lexicographic_numbering.begin() +
2466 shape_info.dofs_per_component_on_cell);
2472 VectorizedArrayType::size(),
2474 dofs_per_cell[
bi], dofs_per_cell[
bj]));
2482 matrix_free.n_active_entries_per_face_batch(
batch) :
2483 matrix_free.n_active_entries_per_cell_batch(
batch);
2489 (
data.batch_type[
b] == 0) ?
2490 (
batch * VectorizedArrayType::size() + v) :
2491 ((
data.batch_type[
b] == 1) ?
2492 matrix_free.get_face_info(
batch).cells_interior[v] :
2493 matrix_free.get_face_info(
batch).cells_exterior[v]);
2495 const auto cell_iterator = matrix_free.get_cell_iterator(
2498 data.dof_numbers[b]);
2500 if (matrix_free.get_mg_level() !=
2502 cell_iterator->get_mg_dof_indices(dof_indices[b]);
2504 cell_iterator->get_dof_indices(dof_indices[b]);
2506 for (
unsigned int j = 0;
j < dofs_per_cell[
b]; ++
j)
2508 dof_indices[b][lexicographic_numbering[b][
j]];
2513 for (
unsigned int j = 0;
j < dofs_per_cell[
bj]; ++
j)
2516 for (
unsigned int i = 0; i < dofs_per_cell[
bi]; ++i)
2517 phi[
bi]->begin_dof_values()[i] =
2518 (
bj ==
bi) ?
static_cast<Number
>(i ==
j) : 0.0;
2523 for (
unsigned int i = 0; i < dofs_per_cell[
bi]; ++i)
2525 matrices[
bi][
bj][v](i,
j) =
2526 phi[
bi]->begin_dof_values()[i][v];
2537 constraints.distribute_local_to_global(
2540 constraints.distribute_local_to_global(
2541 matrices[
bi][
bj][v],
2550 [&](
const auto &,
auto &,
const auto &,
const auto range) {
2555 [&](
const auto &,
auto &,
const auto &,
const auto range) {
2560 [&](
const auto &,
auto &,
const auto &,
const auto range) {
2581 template <
typename CLASS,
2587 typename VectorizedArrayType,
2599 VectorizedArrayType> &)
const,
2605 VectorizedArrayType> &,
2611 VectorizedArrayType> &)
2618 VectorizedArrayType> &)
2623 const unsigned
int first_selected_component)
2630 VectorizedArrayType,
2642 first_selected_component);