300 const std::vector<unsigned int> & renumbering,
301 const std::vector<unsigned int> & constraint_pool_row_index,
302 const std::vector<unsigned char> &irregular_cells)
304 (void)constraint_pool_row_index;
310 std::vector<unsigned int> new_active_fe_index;
312 unsigned int position_cell = 0;
313 for (
unsigned int cell = 0;
317 const unsigned int n_comp =
318 (irregular_cells[cell] > 0 ? irregular_cells[cell] :
323 unsigned int fe_index =
325 for (
unsigned int j = 1; j < n_comp; ++j)
330 new_active_fe_index.push_back(fe_index);
331 position_cell += n_comp;
341 std::vector<std::pair<unsigned int, unsigned int>> new_row_starts(
345 std::vector<unsigned int> new_dof_indices;
346 std::vector<std::pair<unsigned short, unsigned short>>
347 new_constraint_indicator;
348 std::vector<unsigned int> new_plain_indices, new_rowstart_plain;
349 unsigned int position_cell = 0;
353 std::vector<compressed_constraint_kind> new_hanging_node_constraint_masks;
354 new_hanging_node_constraint_masks.reserve(
374 const unsigned int n_vect =
375 (irregular_cells[i] > 0 ? irregular_cells[i] :
379 this->dofs_per_cell[0];
381 for (
unsigned int j = 0; j < n_vect; ++j)
383 const unsigned int cell_no =
386 bool has_hanging_nodes =
false;
394 new_hanging_node_constraint_masks.push_back(
mask);
397 for (
unsigned int comp = 0; comp <
n_components; ++comp)
402 for (
unsigned int comp = 0; comp <
n_components; ++comp)
406 .
first = new_dof_indices.size();
409 .
second = new_constraint_indicator.size();
411 new_dof_indices.insert(
412 new_dof_indices.end(),
418 new_constraint_indicator.push_back(
427 new_plain_indices.size();
428 new_plain_indices.insert(
429 new_plain_indices.end(),
438 for (
unsigned int comp = 0; comp <
n_components; ++comp)
442 .
first = new_dof_indices.size();
445 .
second = new_constraint_indicator.size();
450 new_hanging_node_constraint_masks.push_back(
453 position_cell += n_vect;
460 .first = new_dof_indices.size();
463 .second = new_constraint_indicator.size();
466 new_constraint_indicator.size());
478 const unsigned int index_range =
492 const unsigned int row_length_ind =
498 const std::pair<unsigned short, unsigned short> *
505 for (; con_it != end_con; ++con_it)
509 constraint_pool_row_index.size() - 1);
514 unsigned int n_active_cells = 0;
517 if (irregular_cells[c] > 0)
518 n_active_cells += irregular_cells[c];
531 const std::vector<unsigned char> &irregular_cells)
541 irregular_cells.size());
545 irregular_cells.size());
546 for (
unsigned int i = 0; i < irregular_cells.size(); ++i)
547 if (irregular_cells[i] > 0)
562 std::vector<unsigned int> index_kinds(
563 static_cast<unsigned int>(
567 for (
unsigned int i = 0; i < irregular_cells.size(); ++i)
569 const unsigned int ndofs =
571 const unsigned int n_comp =
575 bool has_constraints =
false;
576 for (
unsigned int j = 0; j < n_comp; ++j)
582 has_constraints =
true;
591 bool indices_are_contiguous = (ndofs > 0);
592 for (
unsigned int j = 0; j < n_comp; ++j)
596 this->dof_indices.data() +
602 for (
unsigned int i = 1; i < ndofs; ++i)
605 indices_are_contiguous =
false;
610 bool indices_are_interleaved_and_contiguous =
615 this->dof_indices.data() +
617 for (
unsigned int k = 0;
618 k < ndofs && indices_are_interleaved_and_contiguous;
620 for (
unsigned int j = 0; j < n_comp; ++j)
624 indices_are_interleaved_and_contiguous =
false;
629 if (indices_are_contiguous ||
630 indices_are_interleaved_and_contiguous)
632 for (
unsigned int j = 0; j < n_comp; ++j)
634 const unsigned int start_index =
647 if (indices_are_interleaved_and_contiguous)
652 for (
unsigned int j = 0; j < n_comp; ++j)
656 else if (indices_are_contiguous)
660 for (
unsigned int j = 0; j < n_comp; ++j)
666 int indices_are_interleaved_and_mixed = 2;
671 for (
unsigned int j = 0; j < n_comp; ++j)
674 for (
unsigned int k = 0;
675 k < ndofs && indices_are_interleaved_and_mixed != 0;
677 for (
unsigned int j = 0; j < n_comp; ++j)
684 indices_are_interleaved_and_mixed = 0;
687 if (indices_are_interleaved_and_mixed == 2)
689 for (
unsigned int j = 0; j < n_comp; ++j)
693 for (
unsigned int j = 0; j < n_comp; ++j)
697 for (
unsigned int j = 0; j < n_comp; ++j)
700 indices_are_interleaved_and_mixed = 1;
703 if (indices_are_interleaved_and_mixed == 1 ||
715 this->dof_indices.data() +
732 bool is_sorted =
true;
733 unsigned int previous = indices[0];
734 for (
unsigned int l = 1; l < n_comp; ++l)
736 const unsigned int current = indices[l * ndofs];
737 if (current <= previous)
744 for (
unsigned int j = 0; j < l; ++j)
745 if (indices[j * ndofs] == current)
761 index_kinds[
static_cast<unsigned int>(
771 auto fix_single_interleaved_indices =
773 if (index_kinds[
static_cast<unsigned int>(
776 index_kinds[
static_cast<unsigned int>(variant)] > 0)
777 for (
unsigned int i = 0; i < irregular_cells.size(); ++i)
789 index_kinds[
static_cast<unsigned int>(
792 index_kinds[
static_cast<unsigned int>(variant)]++;
801 unsigned int n_interleaved =
802 index_kinds[
static_cast<unsigned int>(
804 index_kinds[
static_cast<unsigned int>(
806 index_kinds[
static_cast<unsigned int>(
811 if (n_interleaved > 0 && index_kinds[
static_cast<unsigned int>(
813 for (
unsigned int i = 0; i < irregular_cells.size(); ++i)
819 index_kinds[
static_cast<unsigned int>(
821 index_kinds[
static_cast<unsigned int>(
827 if (n_interleaved > 0 &&
829 index_kinds[
static_cast<unsigned int>(
832 for (
unsigned int i = 0; i < irregular_cells.size(); ++i)
836 index_kinds[
static_cast<unsigned int>(
845 index_kinds[
static_cast<unsigned int>(
850 for (
unsigned int i = 0; i < irregular_cells.size(); ++i)
861 const unsigned int ndofs =
866 unsigned int *interleaved_dof_indices =
870 this->dof_indices_interleaved.size());
880 for (
unsigned int k = 0; k < ndofs; ++k)
882 const unsigned int *my_dof_indices =
dof_indices + k;
883 const unsigned int *end =
885 for (; interleaved_dof_indices != end;
886 ++interleaved_dof_indices, my_dof_indices += ndofs)
887 *interleaved_dof_indices = *my_dof_indices;
897 const unsigned int n_owned_cells,
898 const unsigned int n_lanes,
901 const bool fill_cell_centric,
903 const bool use_vector_data_exchanger_full)
909 std::vector<types::global_dof_index> ghost_indices;
912 for (
unsigned int cell = 0; cell < n_owned_cells; ++cell)
920 const unsigned int fe_index =
929 ghost_indices.push_back(
932 std::sort(ghost_indices.begin(), ghost_indices.end());
933 ghost_indices.erase(std::unique(ghost_indices.begin(),
934 ghost_indices.end()),
935 ghost_indices.end());
937 compressed_set.
add_indices(ghost_indices.begin(), ghost_indices.end());
939 const bool all_ghosts_equal =
940 Utilities::MPI::min<int>(
static_cast<int>(
945 std::shared_ptr<const Utilities::MPI::Partitioner> temp_0;
947 if (all_ghosts_equal)
951 temp_0 = std::make_shared<Utilities::MPI::Partitioner>(
957 if (use_vector_data_exchanger_full ==
false)
963 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
964 temp_0, communicator_sm);
968 std::vector<FaceToCellTopology<1>> all_faces(inner_faces);
969 all_faces.insert(all_faces.end(),
970 ghosted_faces.begin(),
971 ghosted_faces.end());
974 2 * shape_info(0, 0).n_dimensions);
976 for (
unsigned int f = 0; f < all_faces.size(); ++f)
978 cell_and_face_to_faces(all_faces[f].cells_interior[0],
979 all_faces[f].interior_face_no) = f;
980 Assert(all_faces[f].cells_exterior[0] !=
983 cell_and_face_to_faces(all_faces[f].cells_exterior[0],
984 all_faces[f].exterior_face_no) = f;
988 const auto loop_over_faces =
989 [&](
const std::function<
990 void(
const unsigned int,
const unsigned int,
const bool)> &fu) {
991 for (
const auto &face : inner_faces)
994 fu(face.cells_exterior[0], face.exterior_face_no,
false );
998 const auto loop_over_all_faces =
999 [&](
const std::function<
1000 void(
const unsigned int,
const unsigned int,
const bool)> &fu) {
1001 for (
unsigned int c = 0; c < cell_and_face_to_faces.size(0); ++c)
1002 for (
unsigned int d = 0; d < cell_and_face_to_faces.size(1); ++d)
1004 const unsigned int f = cell_and_face_to_faces(c, d);
1008 const unsigned int cell_m = all_faces[f].cells_interior[0];
1009 const unsigned int cell_p = all_faces[f].cells_exterior[0];
1011 const bool ext = c == cell_m;
1016 const unsigned int p = ext ? cell_p : cell_m;
1017 const unsigned int face_no = ext ?
1018 all_faces[f].exterior_face_no :
1019 all_faces[f].interior_face_no;
1021 fu(p, face_no,
true);
1025 const auto process_values =
1027 std::shared_ptr<const Utilities::MPI::Partitioner>
1028 &vector_partitioner_values,
1029 const std::function<void(
1030 const std::function<
void(
1031 const unsigned int,
const unsigned int,
const bool)> &)> &loop) {
1032 bool all_nodal_and_tensorial = shape_info.size(1) == 1;
1034 if (all_nodal_and_tensorial)
1039 if (!si.nodal_at_cell_boundaries ||
1042 all_nodal_and_tensorial =
false;
1045 if (all_nodal_and_tensorial ==
false)
1049 bool has_noncontiguous_cell =
false;
1051 loop([&](
const unsigned int cell_no,
1052 const unsigned int face_no,
1054 const unsigned int index =
1059 const unsigned int stride =
1060 dof_indices_interleave_strides[dof_access_cell][cell_no];
1062 for (unsigned int e = 0; e < n_base_elements; ++e)
1063 for (unsigned int c = 0; c < n_components[e]; ++c)
1065 const ShapeInfo<double> &shape =
1066 shape_info(global_base_element_offset + e, 0);
1067 for (unsigned int j = 0;
1068 j < shape.dofs_per_component_on_face;
1070 ghost_indices.push_back(part.local_to_global(
1072 shape.face_to_cell_index_nodal(face_no, j) *
1074 i += shape.dofs_per_component_on_cell * stride;
1079 has_noncontiguous_cell =
true;
1081 has_noncontiguous_cell =
1082 Utilities::MPI::min<int>(
static_cast<int>(
1083 has_noncontiguous_cell),
1086 std::sort(ghost_indices.begin(), ghost_indices.end());
1087 ghost_indices.erase(std::unique(ghost_indices.begin(),
1088 ghost_indices.end()),
1089 ghost_indices.end());
1092 ghost_indices.end());
1094 const bool all_ghosts_equal =
1095 Utilities::MPI::min<int>(
static_cast<int>(
1099 if (all_ghosts_equal || has_noncontiguous_cell)
1103 vector_partitioner_values =
1104 std::make_shared<Utilities::MPI::Partitioner>(
1107 vector_partitioner_values.get())
1114 const auto process_gradients =
1116 const std::shared_ptr<const Utilities::MPI::Partitioner>
1117 &vector_partitoner_values,
1118 std::shared_ptr<const Utilities::MPI::Partitioner>
1119 &vector_partitioner_gradients,
1120 const std::function<void(
1121 const std::function<
void(
1122 const unsigned int,
const unsigned int,
const bool)> &)> &loop) {
1123 bool all_hermite = shape_info.
size(1) == 1;
1126 for (
unsigned int c = 0; c < n_base_elements; ++c)
1127 if (shape_info(global_base_element_offset + c, 0).element_type !=
1129 all_hermite =
false;
1130 if (all_hermite ==
false ||
1131 vector_partitoner_values.get() == vector_partitioner.get())
1132 vector_partitioner_gradients = vector_partitioner;
1135 loop([&](
const unsigned int cell_no,
1136 const unsigned int face_no,
1138 const unsigned int index =
1139 dof_indices_contiguous[dof_access_cell][cell_no];
1141 index >= part.locally_owned_size()))
1143 const unsigned int stride =
1144 dof_indices_interleave_strides[dof_access_cell][cell_no];
1146 for (unsigned int e = 0; e < n_base_elements; ++e)
1147 for (unsigned int c = 0; c < n_components[e]; ++c)
1149 const ShapeInfo<double> &shape =
1150 shape_info(global_base_element_offset + e, 0);
1151 for (unsigned int j = 0;
1152 j < 2 * shape.dofs_per_component_on_face;
1154 ghost_indices.push_back(part.local_to_global(
1156 shape.face_to_cell_index_hermite(face_no, j) *
1158 i += shape.dofs_per_component_on_cell * stride;
1163 std::sort(ghost_indices.begin(), ghost_indices.end());
1164 ghost_indices.erase(std::unique(ghost_indices.begin(),
1165 ghost_indices.end()),
1166 ghost_indices.end());
1167 IndexSet compressed_set(part.size());
1168 compressed_set.add_indices(ghost_indices.begin(),
1169 ghost_indices.end());
1170 compressed_set.subtract_set(part.locally_owned_range());
1171 const bool all_ghosts_equal =
1172 Utilities::MPI::min<int>(
static_cast<int>(
1173 compressed_set.n_elements() ==
1174 part.ghost_indices().n_elements()),
1175 part.get_mpi_communicator()) != 0;
1176 if (all_ghosts_equal)
1177 vector_partitioner_gradients = vector_partitioner;
1180 vector_partitioner_gradients =
1181 std::make_shared<Utilities::MPI::Partitioner>(
1182 part.locally_owned_range(), part.get_mpi_communicator());
1184 vector_partitioner_gradients.get())
1185 ->set_ghost_indices(compressed_set, part.ghost_indices());
1190 std::shared_ptr<const Utilities::MPI::Partitioner> temp_1, temp_2, temp_3,
1194 process_values(temp_1, loop_over_faces);
1197 process_gradients(temp_1, temp_2, loop_over_faces);
1199 if (fill_cell_centric)
1201 ghost_indices.clear();
1203 process_values(temp_3, loop_over_all_faces);
1205 process_gradients(temp_3, temp_4, loop_over_all_faces);
1209 temp_3 = std::make_shared<Utilities::MPI::Partitioner>(
1210 part.locally_owned_range(), part.get_mpi_communicator());
1211 temp_4 = std::make_shared<Utilities::MPI::Partitioner>(
1212 part.locally_owned_range(), part.get_mpi_communicator());
1215 if (use_vector_data_exchanger_full ==
false)
1217 vector_exchanger_face_variants[1] = std::make_shared<
1218 MatrixFreeFunctions::VectorDataExchange::PartitionerWrapper>(
1220 vector_exchanger_face_variants[2] = std::make_shared<
1221 MatrixFreeFunctions::VectorDataExchange::PartitionerWrapper>(
1223 vector_exchanger_face_variants[3] = std::make_shared<
1224 MatrixFreeFunctions::VectorDataExchange::PartitionerWrapper>(
1226 vector_exchanger_face_variants[4] = std::make_shared<
1227 MatrixFreeFunctions::VectorDataExchange::PartitionerWrapper>(
1232 vector_exchanger_face_variants[1] =
1233 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
1234 temp_1, communicator_sm);
1235 vector_exchanger_face_variants[2] =
1236 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
1237 temp_2, communicator_sm);
1238 vector_exchanger_face_variants[3] =
1239 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
1240 temp_3, communicator_sm);
1241 vector_exchanger_face_variants[4] =
1242 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
1243 temp_4, communicator_sm);