159 const bool use_vector_data_exchanger_full)
166 const std::size_t n_ghosts =
ghost_dofs.size();
174 std::vector<unsigned int> ghost_numbering(n_ghosts);
178 unsigned int n_unique_ghosts = 0;
182 std::vector<std::pair<types::global_dof_index, unsigned int>>
183 ghost_origin(n_ghosts);
184 for (std::size_t i = 0; i < n_ghosts; ++i)
187 ghost_origin[i].second = i;
189 std::sort(ghost_origin.begin(), ghost_origin.end());
192 ghost_numbering[ghost_origin[0].second] = 0;
193 for (std::size_t i = 1; i < n_ghosts; ++i)
195 if (ghost_origin[i].
first > ghost_origin[i - 1].
first + 1)
197 ghost_indices.
add_range(last_contiguous_start,
198 ghost_origin[i - 1].
first + 1);
199 last_contiguous_start = ghost_origin[i].first;
201 if (ghost_origin[i].
first > ghost_origin[i - 1].
first)
203 ghost_numbering[ghost_origin[i].second] = n_unique_ghosts;
206 ghost_indices.
add_range(last_contiguous_start,
207 ghost_origin.back().first + 1);
214 for (std::size_t i = 0; i < n_ghosts; ++i)
215 Assert(ghost_numbering[i] ==
224 const unsigned int n_boundary_cells = boundary_cells.size();
225 for (
unsigned int i = 0; i < n_boundary_cells; ++i)
227 unsigned int *data_ptr =
230 const unsigned int *row_end =
233 for (; data_ptr != row_end; ++data_ptr)
234 *data_ptr = ((*data_ptr < n_owned) ?
236 n_owned + ghost_numbering[*data_ptr - n_owned]);
241 bool has_hanging_nodes =
false;
243 const unsigned int fe_index =
254 for (
unsigned int comp = 0; comp <
n_components; ++comp)
258 if (has_hanging_nodes ||
263 unsigned int *data_ptr =
266 const unsigned int *row_end =
268 for (; data_ptr != row_end; ++data_ptr)
270 ((*data_ptr < n_owned) ?
272 n_owned + ghost_numbering[*data_ptr - n_owned]);
278 std::vector<types::global_dof_index> empty;
288 if (use_vector_data_exchanger_full ==
false)
294 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
303 const std::vector<unsigned int> &renumbering,
304 const std::vector<unsigned int> &constraint_pool_row_index,
305 const std::vector<unsigned char> &irregular_cells)
307 (void)constraint_pool_row_index;
313 std::vector<unsigned int> new_active_fe_index;
315 unsigned int position_cell = 0;
316 for (
unsigned int cell = 0;
320 const unsigned int n_comp =
321 (irregular_cells[cell] > 0 ? irregular_cells[cell] :
326 unsigned int fe_index =
328 for (
unsigned int j = 1; j < n_comp; ++j)
333 new_active_fe_index.push_back(fe_index);
334 position_cell += n_comp;
344 std::vector<std::pair<unsigned int, unsigned int>> new_row_starts(
348 std::vector<unsigned int> new_dof_indices;
349 std::vector<std::pair<unsigned short, unsigned short>>
350 new_constraint_indicator;
351 std::vector<unsigned int> new_plain_indices, new_rowstart_plain;
352 unsigned int position_cell = 0;
356 std::vector<compressed_constraint_kind> new_hanging_node_constraint_masks;
357 new_hanging_node_constraint_masks.reserve(
377 const unsigned int n_vect =
378 (irregular_cells[i] > 0 ? irregular_cells[i] :
382 this->dofs_per_cell[0];
384 for (
unsigned int j = 0; j < n_vect; ++j)
386 const unsigned int cell_no =
389 bool has_hanging_nodes =
false;
397 new_hanging_node_constraint_masks.push_back(
mask);
400 for (
unsigned int comp = 0; comp <
n_components; ++comp)
405 for (
unsigned int comp = 0; comp <
n_components; ++comp)
409 .
first = new_dof_indices.size();
412 .
second = new_constraint_indicator.size();
414 new_dof_indices.insert(
415 new_dof_indices.end(),
421 new_constraint_indicator.push_back(
430 new_plain_indices.size();
431 new_plain_indices.insert(
432 new_plain_indices.end(),
441 for (
unsigned int comp = 0; comp <
n_components; ++comp)
445 .
first = new_dof_indices.size();
448 .
second = new_constraint_indicator.size();
453 new_hanging_node_constraint_masks.push_back(
456 position_cell += n_vect;
463 .first = new_dof_indices.size();
466 .second = new_constraint_indicator.size();
469 new_constraint_indicator.size());
482 const unsigned int index_range =
496 const unsigned int row_length_ind =
504 const std::pair<unsigned short, unsigned short>
512 for (; con_it != end_con; ++con_it)
516 constraint_pool_row_index.size() - 1);
521 unsigned int n_active_cells = 0;
522 for (
unsigned int c = 0;
525 if (irregular_cells[c] > 0)
526 n_active_cells += irregular_cells[c];
539 const std::vector<unsigned char> &irregular_cells)
549 irregular_cells.size());
553 irregular_cells.size());
554 for (
unsigned int i = 0; i < irregular_cells.size(); ++i)
555 if (irregular_cells[i] > 0)
570 std::vector<unsigned int> index_kinds(
571 static_cast<unsigned int>(
575 for (
unsigned int i = 0; i < irregular_cells.size(); ++i)
577 const unsigned int ndofs =
579 const unsigned int n_comp =
583 bool has_constraints =
false;
584 for (
unsigned int j = 0; j < n_comp; ++j)
590 has_constraints =
true;
599 bool indices_are_contiguous = (ndofs > 0);
600 for (
unsigned int j = 0; j < n_comp; ++j)
604 this->dof_indices.data() +
610 for (
unsigned int i = 1; i < ndofs; ++i)
613 indices_are_contiguous =
false;
618 bool indices_are_interleaved_and_contiguous =
623 this->dof_indices.data() +
625 for (
unsigned int k = 0;
626 k < ndofs && indices_are_interleaved_and_contiguous;
628 for (
unsigned int j = 0; j < n_comp; ++j)
632 indices_are_interleaved_and_contiguous =
false;
637 if (indices_are_contiguous ||
638 indices_are_interleaved_and_contiguous)
640 for (
unsigned int j = 0; j < n_comp; ++j)
642 const unsigned int start_index =
655 if (indices_are_interleaved_and_contiguous)
660 for (
unsigned int j = 0; j < n_comp; ++j)
664 else if (indices_are_contiguous)
668 for (
unsigned int j = 0; j < n_comp; ++j)
674 int indices_are_interleaved_and_mixed = 2;
679 for (
unsigned int j = 0; j < n_comp; ++j)
682 for (
unsigned int k = 0;
683 k < ndofs && indices_are_interleaved_and_mixed != 0;
685 for (
unsigned int j = 0; j < n_comp; ++j)
692 indices_are_interleaved_and_mixed = 0;
695 if (indices_are_interleaved_and_mixed == 2)
697 for (
unsigned int j = 0; j < n_comp; ++j)
701 for (
unsigned int j = 0; j < n_comp; ++j)
705 for (
unsigned int j = 0; j < n_comp; ++j)
708 indices_are_interleaved_and_mixed = 1;
711 if (indices_are_interleaved_and_mixed == 1 ||
723 this->dof_indices.data() +
740 bool is_sorted =
true;
741 unsigned int previous = indices[0];
742 for (
unsigned int l = 1; l < n_comp; ++l)
744 const unsigned int current = indices[l * ndofs];
745 if (current <= previous)
752 for (
unsigned int j = 0; j < l; ++j)
753 if (indices[j * ndofs] == current)
769 index_kinds[
static_cast<unsigned int>(
779 auto fix_single_interleaved_indices =
781 if (index_kinds[
static_cast<unsigned int>(
784 index_kinds[
static_cast<unsigned int>(variant)] > 0)
785 for (
unsigned int i = 0; i < irregular_cells.size(); ++i)
797 index_kinds[
static_cast<unsigned int>(
800 index_kinds[
static_cast<unsigned int>(variant)]++;
809 unsigned int n_interleaved =
810 index_kinds[
static_cast<unsigned int>(
812 index_kinds[
static_cast<unsigned int>(
814 index_kinds[
static_cast<unsigned int>(
819 if (n_interleaved > 0 && index_kinds[
static_cast<unsigned int>(
821 for (
unsigned int i = 0; i < irregular_cells.size(); ++i)
827 index_kinds[
static_cast<unsigned int>(
829 index_kinds[
static_cast<unsigned int>(
835 if (n_interleaved > 0 &&
837 index_kinds[
static_cast<unsigned int>(
840 for (
unsigned int i = 0; i < irregular_cells.size(); ++i)
844 index_kinds[
static_cast<unsigned int>(
853 index_kinds[
static_cast<unsigned int>(
858 for (
unsigned int i = 0; i < irregular_cells.size(); ++i)
869 const unsigned int ndofs =
874 unsigned int *interleaved_dof_indices =
878 this->dof_indices_interleaved.size());
888 for (
unsigned int k = 0; k < ndofs; ++k)
890 const unsigned int *my_dof_indices =
dof_indices + k;
891 const unsigned int *end =
893 for (; interleaved_dof_indices != end;
894 ++interleaved_dof_indices, my_dof_indices += ndofs)
895 *interleaved_dof_indices = *my_dof_indices;
905 const unsigned int n_owned_cells,
906 const unsigned int n_lanes,
909 const bool fill_cell_centric,
911 const bool use_vector_data_exchanger_full)
917 std::vector<types::global_dof_index> ghost_indices;
920 for (
unsigned int cell = 0; cell < n_owned_cells; ++cell)
928 const unsigned int fe_index =
937 ghost_indices.push_back(
940 std::sort(ghost_indices.begin(), ghost_indices.end());
942 compressed_set.
add_indices(ghost_indices.begin(), ghost_indices.end());
944 const bool all_ghosts_equal =
945 Utilities::MPI::min<int>(
static_cast<int>(
950 std::shared_ptr<const Utilities::MPI::Partitioner> temp_0;
952 if (all_ghosts_equal)
956 temp_0 = std::make_shared<Utilities::MPI::Partitioner>(
962 if (use_vector_data_exchanger_full ==
false)
968 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
969 temp_0, communicator_sm);
973 std::vector<FaceToCellTopology<1>> all_faces(inner_faces);
974 all_faces.insert(all_faces.end(),
975 ghosted_faces.begin(),
976 ghosted_faces.end());
979 2 * shape_info(0, 0).n_dimensions);
981 for (
unsigned int f = 0; f < all_faces.size(); ++f)
983 cell_and_face_to_faces(all_faces[f].cells_interior[0],
984 all_faces[f].interior_face_no) = f;
985 Assert(all_faces[f].cells_exterior[0] !=
988 cell_and_face_to_faces(all_faces[f].cells_exterior[0],
989 all_faces[f].exterior_face_no) = f;
993 const auto loop_over_faces =
994 [&](
const std::function<
995 void(
const unsigned int,
const unsigned int,
const bool)> &fu) {
996 for (
const auto &face : inner_faces)
999 fu(face.cells_exterior[0], face.exterior_face_no,
false );
1003 const auto loop_over_all_faces =
1004 [&](
const std::function<
1005 void(
const unsigned int,
const unsigned int,
const bool)> &fu) {
1006 for (
unsigned int c = 0; c < cell_and_face_to_faces.size(0); ++c)
1007 for (
unsigned int d = 0; d < cell_and_face_to_faces.size(1); ++d)
1009 const unsigned int f = cell_and_face_to_faces(c, d);
1013 const unsigned int cell_m = all_faces[f].cells_interior[0];
1014 const unsigned int cell_p = all_faces[f].cells_exterior[0];
1016 const bool ext = c == cell_m;
1021 const unsigned int p = ext ? cell_p : cell_m;
1022 const unsigned int face_no = ext ?
1023 all_faces[f].exterior_face_no :
1024 all_faces[f].interior_face_no;
1026 fu(p, face_no,
true);
1030 const auto process_values =
1032 std::shared_ptr<const Utilities::MPI::Partitioner>
1033 &vector_partitioner_values,
1034 const std::function<void(
1035 const std::function<
void(
1036 const unsigned int,
const unsigned int,
const bool)> &)> &loop) {
1037 bool all_nodal_and_tensorial = shape_info.size(1) == 1;
1039 if (all_nodal_and_tensorial)
1044 if (!si.nodal_at_cell_boundaries ||
1047 all_nodal_and_tensorial =
false;
1050 if (all_nodal_and_tensorial ==
false)
1054 bool has_noncontiguous_cell =
false;
1056 loop([&](
const unsigned int cell_no,
1057 const unsigned int face_no,
1059 const unsigned int index =
1064 const unsigned int stride =
1065 dof_indices_interleave_strides[dof_access_cell][cell_no];
1067 for (unsigned int e = 0; e < n_base_elements; ++e)
1068 for (unsigned int c = 0; c < n_components[e]; ++c)
1070 const ShapeInfo<double> &shape =
1071 shape_info(global_base_element_offset + e, 0);
1072 for (unsigned int j = 0;
1073 j < shape.dofs_per_component_on_face;
1075 ghost_indices.push_back(part.local_to_global(
1077 shape.face_to_cell_index_nodal(face_no, j) *
1079 i += shape.dofs_per_component_on_cell * stride;
1084 has_noncontiguous_cell =
true;
1086 has_noncontiguous_cell =
1087 Utilities::MPI::min<int>(
static_cast<int>(
1088 has_noncontiguous_cell),
1091 std::sort(ghost_indices.begin(), ghost_indices.end());
1094 ghost_indices.end());
1096 const bool all_ghosts_equal =
1097 Utilities::MPI::min<int>(
static_cast<int>(
1101 if (all_ghosts_equal || has_noncontiguous_cell)
1105 vector_partitioner_values =
1106 std::make_shared<Utilities::MPI::Partitioner>(
1109 vector_partitioner_values.get())
1116 const auto process_gradients =
1118 const std::shared_ptr<const Utilities::MPI::Partitioner>
1119 &vector_partitoner_values,
1120 std::shared_ptr<const Utilities::MPI::Partitioner>
1121 &vector_partitioner_gradients,
1122 const std::function<void(
1123 const std::function<
void(
1124 const unsigned int,
const unsigned int,
const bool)> &)> &loop) {
1125 bool all_hermite = shape_info.
size(1) == 1;
1128 for (
unsigned int c = 0; c < n_base_elements; ++c)
1129 if (shape_info(global_base_element_offset + c, 0).element_type !=
1131 all_hermite =
false;
1132 if (all_hermite ==
false ||
1133 vector_partitoner_values.get() == vector_partitioner.get())
1134 vector_partitioner_gradients = vector_partitioner;
1137 loop([&](
const unsigned int cell_no,
1138 const unsigned int face_no,
1140 const unsigned int index =
1141 dof_indices_contiguous[dof_access_cell][cell_no];
1143 index >= part.locally_owned_size()))
1145 const unsigned int stride =
1146 dof_indices_interleave_strides[dof_access_cell][cell_no];
1148 for (unsigned int e = 0; e < n_base_elements; ++e)
1149 for (unsigned int c = 0; c < n_components[e]; ++c)
1151 const ShapeInfo<double> &shape =
1152 shape_info(global_base_element_offset + e, 0);
1153 for (unsigned int j = 0;
1154 j < 2 * shape.dofs_per_component_on_face;
1156 ghost_indices.push_back(part.local_to_global(
1158 shape.face_to_cell_index_hermite(face_no, j) *
1160 i += shape.dofs_per_component_on_cell * stride;
1165 std::sort(ghost_indices.begin(), ghost_indices.end());
1166 IndexSet compressed_set(part.size());
1167 compressed_set.add_indices(ghost_indices.begin(),
1168 ghost_indices.end());
1169 compressed_set.subtract_set(part.locally_owned_range());
1170 const bool all_ghosts_equal =
1171 Utilities::MPI::min<int>(
static_cast<int>(
1172 compressed_set.n_elements() ==
1173 part.ghost_indices().n_elements()),
1174 part.get_mpi_communicator()) != 0;
1175 if (all_ghosts_equal)
1176 vector_partitioner_gradients = vector_partitioner;
1179 vector_partitioner_gradients =
1180 std::make_shared<Utilities::MPI::Partitioner>(
1181 part.locally_owned_range(), part.get_mpi_communicator());
1183 vector_partitioner_gradients.get())
1184 ->set_ghost_indices(compressed_set, part.ghost_indices());
1189 std::shared_ptr<const Utilities::MPI::Partitioner> temp_1, temp_2, temp_3,
1193 process_values(temp_1, loop_over_faces);
1196 process_gradients(temp_1, temp_2, loop_over_faces);
1198 if (fill_cell_centric)
1200 ghost_indices.clear();
1202 process_values(temp_3, loop_over_all_faces);
1204 process_gradients(temp_3, temp_4, loop_over_all_faces);
1208 temp_3 = std::make_shared<Utilities::MPI::Partitioner>(
1209 part.locally_owned_range(), part.get_mpi_communicator());
1210 temp_4 = std::make_shared<Utilities::MPI::Partitioner>(
1211 part.locally_owned_range(), part.get_mpi_communicator());
1214 if (use_vector_data_exchanger_full ==
false)
1216 vector_exchanger_face_variants[1] = std::make_shared<
1217 MatrixFreeFunctions::VectorDataExchange::PartitionerWrapper>(
1219 vector_exchanger_face_variants[2] = std::make_shared<
1220 MatrixFreeFunctions::VectorDataExchange::PartitionerWrapper>(
1222 vector_exchanger_face_variants[3] = std::make_shared<
1223 MatrixFreeFunctions::VectorDataExchange::PartitionerWrapper>(
1225 vector_exchanger_face_variants[4] = std::make_shared<
1226 MatrixFreeFunctions::VectorDataExchange::PartitionerWrapper>(
1231 vector_exchanger_face_variants[1] =
1232 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
1233 temp_1, communicator_sm);
1234 vector_exchanger_face_variants[2] =
1235 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
1236 temp_2, communicator_sm);
1237 vector_exchanger_face_variants[3] =
1238 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
1239 temp_3, communicator_sm);
1240 vector_exchanger_face_variants[4] =
1241 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
1242 temp_4, communicator_sm);