157 const bool use_vector_data_exchanger_full)
164 const std::size_t n_ghosts =
ghost_dofs.size();
171 std::vector<unsigned int> ghost_numbering(n_ghosts);
175 unsigned int n_unique_ghosts = 0;
179 std::vector<std::pair<types::global_dof_index, unsigned int>>
180 ghost_origin(n_ghosts);
181 for (std::size_t i = 0; i < n_ghosts; ++i)
184 ghost_origin[i].second = i;
186 std::sort(ghost_origin.begin(), ghost_origin.end());
189 ghost_numbering[ghost_origin[0].second] = 0;
190 for (std::size_t i = 1; i < n_ghosts; ++i)
192 if (ghost_origin[i].
first > ghost_origin[i - 1].
first + 1)
194 ghost_indices.
add_range(last_contiguous_start,
195 ghost_origin[i - 1].
first + 1);
196 last_contiguous_start = ghost_origin[i].first;
198 if (ghost_origin[i].
first > ghost_origin[i - 1].
first)
200 ghost_numbering[ghost_origin[i].second] = n_unique_ghosts;
203 ghost_indices.
add_range(last_contiguous_start,
204 ghost_origin.back().first + 1);
211 for (std::size_t i = 0; i < n_ghosts; ++i)
212 Assert(ghost_numbering[i] ==
221 const unsigned int n_boundary_cells = boundary_cells.size();
222 for (
unsigned int i = 0; i < n_boundary_cells; ++i)
224 unsigned int *data_ptr =
227 const unsigned int *row_end =
230 for (; data_ptr != row_end; ++data_ptr)
231 *data_ptr = ((*data_ptr < n_owned) ?
233 n_owned + ghost_numbering[*data_ptr - n_owned]);
238 bool has_hanging_nodes =
false;
240 const unsigned int fe_index =
251 for (
unsigned int comp = 0; comp <
n_components; ++comp)
255 if (has_hanging_nodes ||
260 unsigned int *data_ptr =
263 const unsigned int *row_end =
265 for (; data_ptr != row_end; ++data_ptr)
267 ((*data_ptr < n_owned) ?
269 n_owned + ghost_numbering[*data_ptr - n_owned]);
275 std::vector<types::global_dof_index> empty;
285 if (use_vector_data_exchanger_full ==
false)
291 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
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());
934 compressed_set.
add_indices(ghost_indices.begin(), ghost_indices.end());
936 const bool all_ghosts_equal =
942 std::shared_ptr<const Utilities::MPI::Partitioner> temp_0;
944 if (all_ghosts_equal)
948 temp_0 = std::make_shared<Utilities::MPI::Partitioner>(
954 if (use_vector_data_exchanger_full ==
false)
960 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
961 temp_0, communicator_sm);
965 std::vector<FaceToCellTopology<1>> all_faces(inner_faces);
966 all_faces.insert(all_faces.end(),
967 ghosted_faces.begin(),
968 ghosted_faces.end());
971 2 * shape_info(0, 0).n_dimensions);
973 for (
unsigned int f = 0; f < all_faces.size(); ++f)
975 cell_and_face_to_faces(all_faces[f].cells_interior[0],
976 all_faces[f].interior_face_no) = f;
977 Assert(all_faces[f].cells_exterior[0] !=
980 cell_and_face_to_faces(all_faces[f].cells_exterior[0],
981 all_faces[f].exterior_face_no) = f;
985 const auto loop_over_faces =
986 [&](
const std::function<
987 void(
const unsigned int,
const unsigned int,
const bool)> &fu) {
988 for (
const auto &face : inner_faces)
991 fu(face.cells_exterior[0], face.exterior_face_no,
false );
995 const auto loop_over_all_faces =
996 [&](
const std::function<
997 void(
const unsigned int,
const unsigned int,
const bool)> &fu) {
998 for (
unsigned int c = 0; c < cell_and_face_to_faces.size(0); ++c)
999 for (
unsigned int d = 0; d < cell_and_face_to_faces.size(1); ++d)
1001 const unsigned int f = cell_and_face_to_faces(c, d);
1005 const unsigned int cell_m = all_faces[f].cells_interior[0];
1006 const unsigned int cell_p = all_faces[f].cells_exterior[0];
1008 const bool ext = c == cell_m;
1013 const unsigned int p = ext ? cell_p : cell_m;
1014 const unsigned int face_no = ext ?
1015 all_faces[f].exterior_face_no :
1016 all_faces[f].interior_face_no;
1018 fu(p, face_no,
true);
1022 const auto process_values =
1024 std::shared_ptr<const Utilities::MPI::Partitioner>
1025 &vector_partitioner_values,
1026 const std::function<void(
1027 const std::function<
void(
1028 const unsigned int,
const unsigned int,
const bool)> &)> &loop) {
1029 bool all_nodal_and_tensorial = shape_info.size(1) == 1;
1031 if (all_nodal_and_tensorial)
1036 if (!si.nodal_at_cell_boundaries ||
1039 all_nodal_and_tensorial =
false;
1042 if (all_nodal_and_tensorial ==
false)
1046 bool has_noncontiguous_cell =
false;
1048 loop([&](
const unsigned int cell_no,
1049 const unsigned int face_no,
1051 const unsigned int index =
1056 const unsigned int stride =
1057 dof_indices_interleave_strides[dof_access_cell][cell_no];
1059 for (unsigned int e = 0; e < n_base_elements; ++e)
1060 for (unsigned int c = 0; c < n_components[e]; ++c)
1062 const ShapeInfo<double> &shape =
1063 shape_info(global_base_element_offset + e, 0);
1064 for (unsigned int j = 0;
1065 j < shape.dofs_per_component_on_face;
1067 ghost_indices.push_back(part.local_to_global(
1069 shape.face_to_cell_index_nodal(face_no, j) *
1071 i += shape.dofs_per_component_on_cell * stride;
1076 has_noncontiguous_cell =
true;
1078 has_noncontiguous_cell =
1080 has_noncontiguous_cell),
1083 std::sort(ghost_indices.begin(), ghost_indices.end());
1086 ghost_indices.end());
1088 const bool all_ghosts_equal =
1093 if (all_ghosts_equal || has_noncontiguous_cell)
1097 vector_partitioner_values =
1098 std::make_shared<Utilities::MPI::Partitioner>(
1101 vector_partitioner_values.get())
1108 const auto process_gradients =
1110 const std::shared_ptr<const Utilities::MPI::Partitioner>
1111 &vector_partitoner_values,
1112 std::shared_ptr<const Utilities::MPI::Partitioner>
1113 &vector_partitioner_gradients,
1114 const std::function<void(
1115 const std::function<
void(
1116 const unsigned int,
const unsigned int,
const bool)> &)> &loop) {
1117 bool all_hermite = shape_info.
size(1) == 1;
1120 for (
unsigned int c = 0; c < n_base_elements; ++c)
1121 if (shape_info(global_base_element_offset + c, 0).element_type !=
1123 all_hermite =
false;
1124 if (all_hermite ==
false ||
1125 vector_partitoner_values.get() == vector_partitioner.get())
1126 vector_partitioner_gradients = vector_partitioner;
1129 loop([&](
const unsigned int cell_no,
1130 const unsigned int face_no,
1132 const unsigned int index =
1133 dof_indices_contiguous[dof_access_cell][cell_no];
1135 index >= part.locally_owned_size()))
1137 const unsigned int stride =
1138 dof_indices_interleave_strides[dof_access_cell][cell_no];
1140 for (unsigned int e = 0; e < n_base_elements; ++e)
1141 for (unsigned int c = 0; c < n_components[e]; ++c)
1143 const ShapeInfo<double> &shape =
1144 shape_info(global_base_element_offset + e, 0);
1145 for (unsigned int j = 0;
1146 j < 2 * shape.dofs_per_component_on_face;
1148 ghost_indices.push_back(part.local_to_global(
1150 shape.face_to_cell_index_hermite(face_no, j) *
1152 i += shape.dofs_per_component_on_cell * stride;
1157 std::sort(ghost_indices.begin(), ghost_indices.end());
1158 IndexSet compressed_set(part.size());
1159 compressed_set.add_indices(ghost_indices.begin(),
1160 ghost_indices.end());
1161 compressed_set.subtract_set(part.locally_owned_range());
1162 const bool all_ghosts_equal =
1164 compressed_set.n_elements() ==
1165 part.ghost_indices().n_elements()),
1166 part.get_mpi_communicator()) != 0;
1167 if (all_ghosts_equal)
1168 vector_partitioner_gradients = vector_partitioner;
1171 vector_partitioner_gradients =
1172 std::make_shared<Utilities::MPI::Partitioner>(
1173 part.locally_owned_range(), part.get_mpi_communicator());
1175 vector_partitioner_gradients.get())
1176 ->set_ghost_indices(compressed_set, part.ghost_indices());
1181 std::shared_ptr<const Utilities::MPI::Partitioner> temp_1, temp_2, temp_3,
1185 process_values(temp_1, loop_over_faces);
1188 process_gradients(temp_1, temp_2, loop_over_faces);
1190 if (fill_cell_centric)
1192 ghost_indices.clear();
1194 process_values(temp_3, loop_over_all_faces);
1196 process_gradients(temp_3, temp_4, loop_over_all_faces);
1200 temp_3 = std::make_shared<Utilities::MPI::Partitioner>(
1201 part.locally_owned_range(), part.get_mpi_communicator());
1202 temp_4 = std::make_shared<Utilities::MPI::Partitioner>(
1203 part.locally_owned_range(), part.get_mpi_communicator());
1206 if (use_vector_data_exchanger_full ==
false)
1208 vector_exchanger_face_variants[1] = std::make_shared<
1209 MatrixFreeFunctions::VectorDataExchange::PartitionerWrapper>(
1211 vector_exchanger_face_variants[2] = std::make_shared<
1212 MatrixFreeFunctions::VectorDataExchange::PartitionerWrapper>(
1214 vector_exchanger_face_variants[3] = std::make_shared<
1215 MatrixFreeFunctions::VectorDataExchange::PartitionerWrapper>(
1217 vector_exchanger_face_variants[4] = std::make_shared<
1218 MatrixFreeFunctions::VectorDataExchange::PartitionerWrapper>(
1223 vector_exchanger_face_variants[1] =
1224 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
1225 temp_1, communicator_sm);
1226 vector_exchanger_face_variants[2] =
1227 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
1228 temp_2, communicator_sm);
1229 vector_exchanger_face_variants[3] =
1230 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
1231 temp_3, communicator_sm);
1232 vector_exchanger_face_variants[4] =
1233 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
1234 temp_4, communicator_sm);
1392 DoFInfo::make_connectivity_graph(
1394 const std::vector<unsigned int> &renumbering,
1397 unsigned int n_rows = (vector_partitioner->local_range().second -
1398 vector_partitioner->local_range().first) +
1399 vector_partitioner->ghost_indices().n_elements();
1406 std::vector<unsigned int> row_lengths(n_rows);
1407 std::vector<std::mutex> mutexes(n_rows / internal::bucket_size_threading +
1412 [
this, &mutexes, &row_lengths](
const unsigned int begin,
1413 const unsigned int end) {
1414 internal::compute_row_lengths(
1415 begin, end, *this, mutexes, row_lengths);
1421 for (
unsigned int row = 0; row < n_rows; ++row)
1422 if (row_lengths[row] <= 1)
1423 row_lengths[row] = 0;
1434 [
this, &row_lengths, &mutexes, &connectivity_dof](
1435 const unsigned int begin,
const unsigned int end) {
1436 internal::fill_connectivity_dofs(
1437 begin, end, *this, row_lengths, mutexes, connectivity_dof);
1444 std::vector<unsigned int> reverse_numbering(task_info.
n_active_cells);
1454 [
this, &reverse_numbering, &connectivity_dof, &connectivity](
1455 const unsigned int begin,
const unsigned int end) {
1456 internal::fill_connectivity(begin,
1469 DoFInfo::compute_dof_renumbering(
1470 std::vector<types::global_dof_index> &renumbering)
1472 const unsigned int locally_owned_size =
1473 vector_partitioner->locally_owned_size();
1474 renumbering.resize(0);
1478 const unsigned int n_components = start_components.back();
1479 const unsigned int n_cell_batches =
1480 n_vectorization_lanes_filled[dof_access_cell].size();
1482 (row_starts.size() - 1) / vectorization_length / n_components,
1484 for (
unsigned int cell_no = 0; cell_no < n_cell_batches; ++cell_no)
1487 if (row_starts[cell_no * n_components * vectorization_length]
1489 row_starts[(cell_no + 1) * n_components * vectorization_length]
1492 const unsigned int ndofs =
1493 dofs_per_cell.size() == 1 ?
1495 (dofs_per_cell[cell_active_fe_index.size() > 0 ?
1496 cell_active_fe_index[cell_no] :
1498 const unsigned int *dof_ind =
1499 dof_indices.data() +
1500 row_starts[cell_no * n_components * vectorization_length].first;
1501 for (
unsigned int i = 0; i < ndofs; ++i)
1502 for (
unsigned int j = 0;
1503 j < n_vectorization_lanes_filled[dof_access_cell][cell_no];
1505 if (dof_ind[j * ndofs + i] < locally_owned_size)
1506 if (renumbering[dof_ind[j * ndofs + i]] ==
1508 renumbering[dof_ind[j * ndofs + i]] = counter++;
1515 dof_index = counter++;
1519 dof_index = vector_partitioner->local_to_global(dof_index);