16#ifndef dealii_face_setup_internal_h
17#define dealii_face_setup_internal_h
41 namespace MatrixFreeFunctions
83 const unsigned int mg_level,
84 const bool hold_all_faces_to_owned_cells,
85 const bool build_inner_faces,
86 std::vector<std::pair<unsigned int, unsigned int>> &cell_levels);
97 const std::vector<std::pair<unsigned int, unsigned int>> &cell_levels,
108 const unsigned int face_no,
109 const typename ::Triangulation<dim>::cell_iterator &cell,
110 const unsigned int number_cell_interior,
111 const typename ::Triangulation<dim>::cell_iterator &neighbor,
112 const unsigned int number_cell_exterior,
113 const bool is_mixed_mesh);
143 template <
int vectorization_w
idth>
147 const std::vector<bool> &hard_vectorization_boundary,
148 std::vector<unsigned int> &face_partition_data,
159 : use_active_cells(true)
166 FaceSetup<dim>::initialize(
168 const unsigned int mg_level,
169 const bool hold_all_faces_to_owned_cells,
170 const bool build_inner_faces,
171 std::vector<std::pair<unsigned int, unsigned int>> &cell_levels)
178 if (use_active_cells)
179 for (
const auto &cell_level : cell_levels)
181 typename ::Triangulation<dim>::cell_iterator dcell(
190 at_processor_boundary.resize(cell_levels.size(),
false);
193 FaceCategory::locally_active_done_elsewhere);
197 std::map<types::subdomain_id, FaceIdentifier>
198 inner_faces_at_proc_boundary;
204 for (
unsigned int i = 0; i < cell_levels.size(); ++i)
206 if (i > 0 && cell_levels[i] == cell_levels[i - 1])
208 typename ::Triangulation<dim>::cell_iterator dcell(
210 for (
const unsigned int f : dcell->face_indices())
212 if (dcell->at_boundary(f) && !dcell->has_periodic_neighbor(f))
214 typename ::Triangulation<dim>::cell_iterator neighbor =
215 dcell->neighbor_or_periodic_neighbor(f);
221 const CellId id_mine = dcell->id();
222 if (use_active_cells && neighbor->has_children())
223 for (
unsigned int c = 0;
224 c < (dcell->has_periodic_neighbor(f) ?
225 dcell->periodic_neighbor(f)
226 ->face(dcell->periodic_neighbor_face_no(f))
228 dcell->face(f)->n_children());
231 typename ::Triangulation<dim>::cell_iterator
233 dcell->at_boundary(f) ?
234 dcell->periodic_neighbor_child_on_subface(f, c) :
235 dcell->neighbor_child_on_subface(f, c);
237 neighbor_c->subdomain_id();
238 if (my_domain < neigh_domain)
239 inner_faces_at_proc_boundary[neigh_domain]
240 .n_hanging_faces_larger_subdomain++;
241 else if (my_domain > neigh_domain)
242 inner_faces_at_proc_boundary[neigh_domain]
243 .n_hanging_faces_smaller_subdomain++;
248 use_active_cells ? neighbor->subdomain_id() :
249 neighbor->level_subdomain_id();
250 if (neighbor->level() < dcell->level() &&
253 if (my_domain < neigh_domain)
254 inner_faces_at_proc_boundary[neigh_domain]
255 .n_hanging_faces_smaller_subdomain++;
256 else if (my_domain > neigh_domain)
257 inner_faces_at_proc_boundary[neigh_domain]
258 .n_hanging_faces_larger_subdomain++;
260 else if (neighbor->level() == dcell->level() &&
261 my_domain != neigh_domain)
267 const CellId id_neigh = neighbor->id();
268 if (my_domain < neigh_domain)
269 inner_faces_at_proc_boundary[neigh_domain]
270 .shared_faces.emplace_back(id_mine, id_neigh);
272 inner_faces_at_proc_boundary[neigh_domain]
273 .shared_faces.emplace_back(id_neigh, id_mine);
282 for (
auto &inner_face : inner_faces_at_proc_boundary)
284 Assert(inner_face.first != my_domain,
286 std::sort(inner_face.second.shared_faces.begin(),
287 inner_face.second.shared_faces.end());
288 inner_face.second.shared_faces.erase(
289 std::unique(inner_face.second.shared_faces.begin(),
290 inner_face.second.shared_faces.end()),
291 inner_face.second.shared_faces.end());
296# if defined(DEAL_II_WITH_MPI) && defined(DEBUG)
298 if (const ::parallel::TriangulationBase<dim> *ptria =
299 dynamic_cast<const ::parallel::TriangulationBase<dim>
301 comm = ptria->get_mpi_communicator();
304 unsigned int mysize = inner_face.second.shared_faces.size();
307 int ierr = MPI_Sendrecv(&mysize,
316 600 + inner_face.first,
321 mysize = inner_face.second.n_hanging_faces_smaller_subdomain;
322 ierr = MPI_Sendrecv(&mysize,
331 700 + inner_face.first,
336 mysize = inner_face.second.n_hanging_faces_larger_subdomain;
337 ierr = MPI_Sendrecv(&mysize,
346 800 + inner_face.first,
368 std::vector<std::tuple<CellId, CellId, unsigned int>> other_range(
369 inner_face.second.shared_faces.size());
370 for (
unsigned int i = 0; i < other_range.size(); ++i)
372 std::make_tuple(inner_face.second.shared_faces[i].second,
373 inner_face.second.shared_faces[i].first,
375 std::sort(other_range.begin(), other_range.end());
383 unsigned int n_faces_lower_proc = 0, n_faces_higher_proc = 0;
384 std::vector<signed char> assignment(other_range.size(), 0);
385 if (inner_face.second.shared_faces.size() > 0)
389 unsigned int count = 0;
390 for (
unsigned int i = 1;
391 i < inner_face.second.shared_faces.size();
393 if (inner_face.second.shared_faces[i].first ==
394 inner_face.second.shared_faces[i - 1 - count].first)
401 for (
unsigned int k = 0; k <= count; ++k)
402 assignment[i - 1 - k] = 1;
403 n_faces_higher_proc += count + 1;
413 for (
unsigned int i = 1; i < other_range.size(); ++i)
414 if (std::get<0>(other_range[i]) ==
415 std::get<0>(other_range[i - 1 - count]))
422 for (
unsigned int k = 0; k <= count; ++k)
425 .shared_faces[std::get<2>(
429 .shared_faces[std::get<2>(
430 other_range[i - 1 - k])]
435 if (assignment[std::get<2>(
436 other_range[i - 1 - k])] == 0)
438 assignment[std::get<2>(
439 other_range[i - 1 - k])] = -1;
440 ++n_faces_lower_proc;
454 n_faces_lower_proc +=
455 inner_face.second.n_hanging_faces_smaller_subdomain;
456 n_faces_higher_proc +=
457 inner_face.second.n_hanging_faces_larger_subdomain;
458 const unsigned int n_total_faces_at_proc_boundary =
459 (inner_face.second.shared_faces.size() +
460 inner_face.second.n_hanging_faces_smaller_subdomain +
461 inner_face.second.n_hanging_faces_larger_subdomain);
462 unsigned int split_index = n_total_faces_at_proc_boundary / 2;
463 if (split_index < n_faces_lower_proc)
465 else if (split_index <
466 n_total_faces_at_proc_boundary - n_faces_higher_proc)
467 split_index -= n_faces_lower_proc;
469 split_index = n_total_faces_at_proc_boundary -
470 n_faces_higher_proc - n_faces_lower_proc;
473# if defined(DEAL_II_WITH_MPI) && defined(DEBUG)
474 ierr = MPI_Sendrecv(&split_index,
483 900 + inner_face.first,
488 ierr = MPI_Sendrecv(&n_faces_lower_proc,
497 1000 + inner_face.first,
502 ierr = MPI_Sendrecv(&n_faces_higher_proc,
511 1100 + inner_face.first,
519 std::vector<std::pair<CellId, CellId>> owned_faces_lower,
521 for (
unsigned int i = 0; i < assignment.size(); ++i)
522 if (assignment[i] < 0)
523 owned_faces_lower.push_back(
524 inner_face.second.shared_faces[i]);
525 else if (assignment[i] > 0)
526 owned_faces_higher.push_back(
527 inner_face.second.shared_faces[i]);
529 inner_face.second.shared_faces.size() + 1 -
530 owned_faces_lower.size() -
531 owned_faces_higher.size());
533 unsigned int i = 0, c = 0;
534 for (; i < assignment.size() && c < split_index; ++i)
535 if (assignment[i] == 0)
537 owned_faces_lower.push_back(
538 inner_face.second.shared_faces[i]);
541 for (; i < assignment.size(); ++i)
542 if (assignment[i] == 0)
544 owned_faces_higher.push_back(
545 inner_face.second.shared_faces[i]);
551 std::vector<std::pair<CellId, CellId>> check_faces;
552 check_faces.insert(check_faces.end(),
553 owned_faces_lower.begin(),
554 owned_faces_lower.end());
555 check_faces.insert(check_faces.end(),
556 owned_faces_higher.begin(),
557 owned_faces_higher.end());
558 std::sort(check_faces.begin(), check_faces.end());
560 inner_face.second.shared_faces.size());
561 for (
unsigned int i = 0; i < check_faces.size(); ++i)
562 Assert(check_faces[i] == inner_face.second.shared_faces[i],
567 if (my_domain < inner_face.first)
568 inner_face.second.shared_faces.swap(owned_faces_lower);
570 inner_face.second.shared_faces.swap(owned_faces_higher);
572 std::sort(inner_face.second.shared_faces.begin(),
573 inner_face.second.shared_faces.end());
579 std::set<std::pair<unsigned int, unsigned int>> ghost_cells;
580 for (
unsigned int i = 0; i < cell_levels.size(); ++i)
582 typename ::Triangulation<dim>::cell_iterator dcell(
584 if (use_active_cells)
586 for (
const auto f : dcell->face_indices())
588 if (dcell->at_boundary(f) && !dcell->has_periodic_neighbor(f))
589 face_is_owned[dcell->face(f)->index()] =
590 FaceCategory::locally_active_at_boundary;
591 else if (!build_inner_faces)
596 else if ((dcell->at_boundary(f) ==
false ||
597 dcell->has_periodic_neighbor(f)) &&
599 dcell->neighbor_or_periodic_neighbor(f)->level() <
602 face_is_owned[dcell->face(f)->index()] =
603 FaceCategory::multigrid_refinement_edge;
607 typename ::Triangulation<dim>::cell_iterator neighbor =
608 dcell->neighbor_or_periodic_neighbor(f);
611 if (use_active_cells && neighbor->has_children() &&
612 hold_all_faces_to_owned_cells ==
false)
615 bool add_to_ghost =
false;
617 id1 = use_active_cells ? dcell->subdomain_id() :
618 dcell->level_subdomain_id(),
619 id2 = use_active_cells ?
620 (neighbor->has_children() ?
621 dcell->neighbor_child_on_subface(f, 0)
623 neighbor->subdomain_id()) :
624 neighbor->level_subdomain_id();
632 (use_active_cells ==
false || neighbor->is_active())) ||
633 dcell->level() > neighbor->level() ||
635 inner_faces_at_proc_boundary[id2].shared_faces.begin(),
636 inner_faces_at_proc_boundary[id2].shared_faces.end(),
637 std::make_pair(id1 < id2 ? dcell->
id() : neighbor->id(),
638 id1 < id2 ? neighbor->id() :
641 face_is_owned[dcell->face(f)->index()] =
642 FaceCategory::locally_active_done_here;
643 if (dcell->level() == neighbor->level())
646 ->face(dcell->has_periodic_neighbor(f) ?
647 dcell->periodic_neighbor_face_no(f) :
648 dcell->neighbor_face_no(f))
650 FaceCategory::locally_active_done_here;
656 if (use_active_cells)
658 (dcell->subdomain_id() != neighbor->subdomain_id());
660 add_to_ghost = (dcell->level_subdomain_id() !=
661 neighbor->level_subdomain_id());
663 else if (hold_all_faces_to_owned_cells ==
true)
666 face_is_owned[dcell->face(f)->index()] =
667 FaceCategory::ghosted;
668 if (use_active_cells)
670 if (neighbor->has_children())
671 for (
unsigned int s = 0;
672 s < dcell->face(f)->n_children();
674 if (dcell->at_boundary(f))
677 ->periodic_neighbor_child_on_subface(f,
680 dcell->subdomain_id())
685 if (dcell->neighbor_child_on_subface(f, s)
687 dcell->subdomain_id())
691 add_to_ghost = (dcell->subdomain_id() !=
692 neighbor->subdomain_id());
695 add_to_ghost = (dcell->level_subdomain_id() !=
696 neighbor->level_subdomain_id());
701 if (use_active_cells && neighbor->has_children())
702 for (
unsigned int s = 0;
703 s < dcell->face(f)->n_children();
706 typename ::Triangulation<dim>::cell_iterator
708 dcell->at_boundary(f) ?
709 dcell->periodic_neighbor_child_on_subface(f,
711 dcell->neighbor_child_on_subface(f, s);
712 if (neighbor_child->subdomain_id() !=
713 dcell->subdomain_id())
715 std::pair<unsigned int, unsigned int>(
716 neighbor_child->level(),
717 neighbor_child->index()));
721 std::pair<unsigned int, unsigned int>(
722 neighbor->level(), neighbor->index()));
723 at_processor_boundary[i] =
true;
731 for (
const auto &ghost_cell : ghost_cells)
732 cell_levels.push_back(ghost_cell);
739 FaceSetup<dim>::generate_faces(
741 const std::vector<std::pair<unsigned int, unsigned int>> &cell_levels,
748 std::map<std::pair<unsigned int, unsigned int>,
unsigned int>
750 for (
unsigned int cell = 0; cell < cell_levels.size(); ++cell)
751 if (cell == 0 || cell_levels[cell] != cell_levels[cell - 1])
753 typename ::Triangulation<dim>::cell_iterator dcell(
755 cell_levels[cell].
first,
756 cell_levels[cell].
second);
757 std::pair<unsigned int, unsigned int> level_index(dcell->level(),
759 map_to_vectorized[level_index] = cell;
763 const unsigned int vectorization_length = task_info.vectorization_length;
764 task_info.face_partition_data.resize(
765 task_info.cell_partition_data.size() - 1, 0);
766 task_info.boundary_partition_data.resize(
767 task_info.cell_partition_data.size() - 1, 0);
768 std::vector<unsigned char> face_visited(face_is_owned.size(), 0);
769 for (
unsigned int partition = 0;
770 partition < task_info.cell_partition_data.size() - 2;
773 unsigned int boundary_counter = 0;
774 unsigned int inner_counter = 0;
775 for (
unsigned int cell = task_info.cell_partition_data[partition] *
776 vectorization_length;
777 cell < task_info.cell_partition_data[
partition + 1] *
778 vectorization_length;
780 if (cell == 0 || cell_levels[cell] != cell_levels[cell - 1])
782 typename ::Triangulation<dim>::cell_iterator dcell(
784 cell_levels[cell].
first,
785 cell_levels[cell].
second);
786 for (
const auto f : dcell->face_indices())
789 if (face_is_owned[dcell->face(f)->index()] ==
790 FaceCategory::locally_active_at_boundary)
794 FaceToCellTopology<1> info;
795 info.cells_interior[0] = cell;
797 info.interior_face_no = f;
798 info.exterior_face_no = dcell->face(f)->boundary_id();
801 (dcell->face(f)->reference_cell() !=
806 info.face_orientation = 0;
807 boundary_faces.push_back(info);
809 face_visited[dcell->face(f)->index()]++;
814 typename ::Triangulation<dim>::cell_iterator
815 neighbor = dcell->neighbor_or_periodic_neighbor(f);
816 if (use_active_cells && neighbor->has_children() &&
821 for (
unsigned int c = 0;
822 c < dcell->face(f)->n_children();
825 typename ::Triangulation<
826 dim>::cell_iterator neighbor_c =
827 dcell->at_boundary(f) ?
828 dcell->periodic_neighbor_child_on_subface(
830 dcell->neighbor_child_on_subface(f, c);
832 neighbor_c->subdomain_id();
833 const unsigned int neighbor_face_no =
834 dcell->has_periodic_neighbor(f) ?
835 dcell->periodic_neighbor_face_no(f) :
836 dcell->neighbor_face_no(f);
837 if (neigh_domain != dcell->subdomain_id() ||
839 [dcell->face(f)->child(c)->index()] ==
842 std::pair<unsigned int, unsigned int>
843 level_index(neighbor_c->level(),
844 neighbor_c->index());
846 [dcell->face(f)->child(c)->index()] ==
847 FaceCategory::locally_active_done_here)
850 inner_faces.push_back(create_face(
853 map_to_vectorized[level_index],
858 else if (face_is_owned[dcell->face(f)
861 FaceCategory::ghosted ||
862 face_is_owned[dcell->face(f)
864 FaceCategory::ghosted)
866 inner_ghost_faces.push_back(create_face(
869 map_to_vectorized[level_index],
875 Assert(face_is_owned[dcell->face(f)
878 locally_active_done_elsewhere,
884 [dcell->face(f)->child(c)->index()] = 1;
895 if (face_visited[dcell->face(f)->index()] == 0 &&
896 !(neighbor->has_children()))
898 std::pair<unsigned int, unsigned int>
899 level_index(neighbor->level(),
901 if (face_is_owned[dcell->face(f)->index()] ==
902 FaceCategory::locally_active_done_here)
904 Assert(use_active_cells ||
909 inner_faces.push_back(create_face(
914 map_to_vectorized[level_index],
916 face_visited[dcell->face(f)->index()] = 1;
918 else if (face_is_owned[dcell->face(f)
920 FaceCategory::ghosted)
922 inner_ghost_faces.push_back(create_face(
927 map_to_vectorized[level_index],
929 face_visited[dcell->face(f)->index()] = 1;
932 if (face_is_owned[dcell->face(f)->index()] ==
933 FaceCategory::multigrid_refinement_edge)
935 refinement_edge_faces.push_back(
940 refinement_edge_faces.size(),
947 use_active_cells ? dcell->subdomain_id() :
948 dcell->level_subdomain_id();
950 use_active_cells ? neighbor->subdomain_id() :
951 neighbor->level_subdomain_id();
952 if (neigh_domain != my_domain ||
953 face_visited[dcell->face(f)->index()] == 1)
955 std::pair<unsigned int, unsigned int>
956 level_index(neighbor->level(),
958 if (face_is_owned[dcell->face(f)->index()] ==
959 FaceCategory::locally_active_done_here)
961 Assert(use_active_cells ||
966 inner_faces.push_back(create_face(
971 map_to_vectorized[level_index],
974 else if (face_is_owned[dcell->face(f)
976 FaceCategory::ghosted)
978 inner_ghost_faces.push_back(create_face(
983 map_to_vectorized[level_index],
989 face_visited[dcell->face(f)->index()] = 1;
990 if (dcell->has_periodic_neighbor(f))
994 dcell->periodic_neighbor_face_no(f))
997 if (face_is_owned[dcell->face(f)->index()] ==
998 FaceCategory::multigrid_refinement_edge)
1000 refinement_edge_faces.push_back(
1005 refinement_edge_faces.size(),
1012 task_info.face_partition_data[
partition + 1] =
1013 task_info.face_partition_data[
partition] + inner_counter;
1014 task_info.boundary_partition_data[
partition + 1] =
1015 task_info.boundary_partition_data[
partition] + boundary_counter;
1017 task_info.ghost_face_partition_data.resize(2);
1018 task_info.ghost_face_partition_data[0] = 0;
1019 task_info.ghost_face_partition_data[1] = inner_ghost_faces.size();
1020 task_info.refinement_edge_face_partition_data.resize(2);
1021 task_info.refinement_edge_face_partition_data[0] = 0;
1022 task_info.refinement_edge_face_partition_data[1] =
1023 refinement_edge_faces.size();
1029 FaceToCellTopology<1>
1030 FaceSetup<dim>::create_face(
1031 const unsigned int face_no,
1032 const typename ::Triangulation<dim>::cell_iterator &cell,
1033 const unsigned int number_cell_interior,
1034 const typename ::Triangulation<dim>::cell_iterator &neighbor,
1035 const unsigned int number_cell_exterior,
1036 const bool is_mixed_mesh)
1038 FaceToCellTopology<1> info;
1039 info.cells_interior[0] = number_cell_interior;
1040 info.cells_exterior[0] = number_cell_exterior;
1041 info.interior_face_no = face_no;
1042 if (cell->has_periodic_neighbor(face_no))
1043 info.exterior_face_no = cell->periodic_neighbor_face_no(face_no);
1045 info.exterior_face_no = cell->neighbor_face_no(face_no);
1047 info.face_type = is_mixed_mesh ?
1048 (cell->face(face_no)->reference_cell() !=
1056 if (dim > 1 && cell->level() > neighbor->level())
1058 if (cell->has_periodic_neighbor(face_no))
1059 info.subface_index =
1060 cell->periodic_neighbor_of_coarser_periodic_neighbor(face_no)
1063 info.subface_index =
1064 cell->neighbor_of_coarser_neighbor(face_no).second;
1068 if (dim == 3 && cell->has_periodic_neighbor(face_no))
1071 cell->get_triangulation()
1072 .get_periodic_face_map()
1073 .at({cell, face_no})
1075 const auto [orientation, rotation, flip] =
1077 exterior_face_orientation);
1079 info.face_orientation =
1080 (orientation ? 0u : 1u) + 2 * rotation + 4 * flip;
1085 info.face_orientation = 0;
1086 const unsigned int interior_face_orientation =
1087 !cell->face_orientation(face_no) + 2 * cell->face_rotation(face_no) +
1088 4 * cell->face_flip(face_no);
1089 const unsigned int exterior_face_orientation =
1090 !neighbor->face_orientation(info.exterior_face_no) +
1091 2 * neighbor->face_rotation(info.exterior_face_no) +
1092 4 * neighbor->face_flip(info.exterior_face_no);
1093 if (interior_face_orientation != 0)
1095 info.face_orientation = 8 + interior_face_orientation;
1096 Assert(exterior_face_orientation == 0,
1098 "Face seems to be wrongly oriented from both sides"));
1101 info.face_orientation = exterior_face_orientation;
1105 if (cell->level() > neighbor->level() && exterior_face_orientation > 0)
1108 ShapeInfo<double>::compute_orientation_table(2);
1109 const std::array<unsigned int, 8> inverted_orientations{
1110 {0, 1, 6, 3, 4, 5, 2, 7}};
1111 info.subface_index =
1112 orientation[inverted_orientations[exterior_face_orientation]]
1113 [info.subface_index];
1128 compare_faces_for_vectorization(
1129 const FaceToCellTopology<1> &face1,
1130 const FaceToCellTopology<1> &face2,
1131 const std::vector<unsigned int> &active_fe_indices,
1132 const unsigned int length)
1134 if (face1.interior_face_no != face2.interior_face_no)
1136 if (face1.exterior_face_no != face2.exterior_face_no)
1138 if (face1.subface_index != face2.subface_index)
1140 if (face1.face_orientation != face2.face_orientation)
1142 if (face1.face_type != face2.face_type)
1145 if (active_fe_indices.size() > 0)
1147 if (active_fe_indices[face1.cells_interior[0] / length] !=
1148 active_fe_indices[face2.cells_interior[0] / length])
1152 if (active_fe_indices[face1.cells_exterior[0] / length] !=
1153 active_fe_indices[face2.cells_exterior[0] / length])
1168 template <
int length>
1169 struct FaceComparator
1171 FaceComparator(
const std::vector<unsigned int> &active_fe_indices)
1172 : active_fe_indices(active_fe_indices)
1176 operator()(
const FaceToCellTopology<length> &face1,
1177 const FaceToCellTopology<length> &face2)
const
1180 if (face1.face_type < face2.face_type)
1182 else if (face1.face_type > face2.face_type)
1186 if (active_fe_indices.size() > 0)
1189 if (active_fe_indices[face1.cells_interior[0] / length] <
1190 active_fe_indices[face2.cells_interior[0] / length])
1192 else if (active_fe_indices[face1.cells_interior[0] / length] >
1193 active_fe_indices[face2.cells_interior[0] / length])
1199 if (active_fe_indices[face1.cells_exterior[0] / length] <
1200 active_fe_indices[face2.cells_exterior[0] / length])
1202 else if (active_fe_indices[face1.cells_exterior[0] / length] >
1203 active_fe_indices[face2.cells_exterior[0] / length])
1208 for (
unsigned int i = 0; i < length; ++i)
1209 if (face1.cells_interior[i] < face2.cells_interior[i])
1211 else if (face1.cells_interior[i] > face2.cells_interior[i])
1213 for (
unsigned int i = 0; i < length; ++i)
1214 if (face1.cells_exterior[i] < face2.cells_exterior[i])
1216 else if (face1.cells_exterior[i] > face2.cells_exterior[i])
1218 if (face1.interior_face_no < face2.interior_face_no)
1220 else if (face1.interior_face_no > face2.interior_face_no)
1222 if (face1.exterior_face_no < face2.exterior_face_no)
1224 else if (face1.exterior_face_no > face2.exterior_face_no)
1237 const std::vector<unsigned int> &active_fe_indices;
1242 template <
int vectorization_w
idth>
1245 const std::vector<FaceToCellTopology<1>> &faces_in,
1246 const std::vector<bool> &hard_vectorization_boundary,
1247 std::vector<unsigned int> &face_partition_data,
1248 std::vector<FaceToCellTopology<vectorization_width>> &faces_out,
1249 const std::vector<unsigned int> &active_fe_indices)
1251 FaceToCellTopology<vectorization_width> face_batch;
1252 std::vector<std::vector<unsigned int>> faces_type;
1254 unsigned int face_start = face_partition_data[0],
1255 face_end = face_partition_data[0];
1257 face_partition_data[0] = faces_out.size();
1258 for (
unsigned int partition = 0;
1259 partition < face_partition_data.size() - 1;
1262 std::vector<std::vector<unsigned int>> new_faces_type;
1265 face_start = face_end;
1266 face_end = face_partition_data[
partition + 1];
1273 for (
unsigned int face = face_start; face < face_end; ++face)
1275 for (
auto &face_type : faces_type)
1278 if (compare_faces_for_vectorization(faces_in[face],
1279 faces_in[face_type[0]],
1281 vectorization_width))
1283 face_type.push_back(face);
1287 faces_type.emplace_back(1, face);
1293 FaceComparator<vectorization_width> face_comparator(
1295 std::set<FaceToCellTopology<vectorization_width>,
1296 FaceComparator<vectorization_width>>
1297 new_faces(face_comparator);
1298 for (
const auto &face_type : faces_type)
1300 face_batch.face_type = faces_in[face_type[0]].face_type;
1301 face_batch.interior_face_no =
1302 faces_in[face_type[0]].interior_face_no;
1303 face_batch.exterior_face_no =
1304 faces_in[face_type[0]].exterior_face_no;
1305 face_batch.subface_index = faces_in[face_type[0]].subface_index;
1306 face_batch.face_orientation =
1307 faces_in[face_type[0]].face_orientation;
1308 unsigned int no_faces = face_type.size();
1309 std::vector<unsigned char> touched(no_faces, 0);
1315 unsigned int n_vectorized = 0;
1316 for (
unsigned int f = 0; f < no_faces; ++f)
1317 if (faces_in[face_type[f]].cells_interior[0] %
1318 vectorization_width ==
1321 bool is_contiguous =
true;
1322 if (f + vectorization_width > no_faces)
1323 is_contiguous =
false;
1325 for (
unsigned int v = 1; v < vectorization_width; ++v)
1326 if (faces_in[face_type[f + v]].cells_interior[0] !=
1327 faces_in[face_type[f]].cells_interior[0] + v)
1328 is_contiguous =
false;
1333 vectorization_width + 1);
1334 for (
unsigned int v = 0; v < vectorization_width; ++v)
1336 face_batch.cells_interior[v] =
1337 faces_in[face_type[f + v]].cells_interior[0];
1338 face_batch.cells_exterior[v] =
1339 faces_in[face_type[f + v]].cells_exterior[0];
1342 new_faces.insert(face_batch);
1343 f += vectorization_width - 1;
1344 n_vectorized += vectorization_width;
1348 std::vector<unsigned int> untouched;
1349 untouched.reserve(no_faces - n_vectorized);
1350 for (
unsigned int f = 0; f < no_faces; ++f)
1351 if (touched[f] == 0)
1352 untouched.push_back(f);
1354 for (
const auto f : untouched)
1356 face_batch.cells_interior[v] =
1357 faces_in[face_type[f]].cells_interior[0];
1358 face_batch.cells_exterior[v] =
1359 faces_in[face_type[f]].cells_exterior[0];
1361 if (v == vectorization_width)
1363 new_faces.insert(face_batch);
1367 if (v > 0 && v < vectorization_width)
1370 if (hard_vectorization_boundary[partition + 1] ||
1371 partition == face_partition_data.size() - 2)
1373 for (; v < vectorization_width; ++v)
1376 face_batch.cells_interior[v] =
1378 face_batch.cells_exterior[v] =
1381 new_faces.insert(face_batch);
1386 std::vector<unsigned int> untreated(v);
1387 for (
unsigned int f = 0; f < v; ++f)
1388 untreated[f] = face_type[*(untouched.end() - 1 - f)];
1389 new_faces_type.push_back(untreated);
1395 for (
auto it = new_faces.begin(); it != new_faces.end(); ++it)
1396 faces_out.push_back(*it);
1397 face_partition_data[
partition + 1] += new_faces.size();
1400 faces_type = std::move(new_faces_type);
1406 for (
const auto &face_type : faces_type)
1410 unsigned int nfaces = 0;
1411 for (
unsigned int i = face_partition_data[0];
1412 i < face_partition_data.back();
1414 for (
unsigned int v = 0; v < vectorization_width; ++v)
1415 nfaces += (faces_out[i].cells_interior[v] !=
1419 std::vector<std::pair<unsigned int, unsigned int>> in_faces,
1421 for (
const auto &face_in : faces_in)
1422 in_faces.emplace_back(face_in.cells_interior[0],
1423 face_in.cells_exterior[0]);
1424 for (
unsigned int i = face_partition_data[0];
1425 i < face_partition_data.back();
1427 for (
unsigned int v = 0;
1428 v < vectorization_width && faces_out[i].cells_interior[v] !=
1431 out_faces.emplace_back(faces_out[i].cells_interior[v],
1432 faces_out[i].cells_exterior[v]);
1433 std::sort(in_faces.begin(), in_faces.end());
1434 std::sort(out_faces.begin(), out_faces.end());
1436 for (
unsigned int i = 0; i < in_faces.size(); ++i)
#define DEAL_II_NAMESPACE_OPEN
constexpr bool running_in_debug_mode()
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
constexpr const ReferenceCell & get_hypercube()
void collect_faces_vectorization(const std::vector< FaceToCellTopology< 1 > > &faces_in, const std::vector< bool > &hard_vectorization_boundary, std::vector< unsigned int > &face_partition_data, std::vector< FaceToCellTopology< vectorization_width > > &faces_out)
std::tuple< bool, bool, bool > split_face_orientation(const types::geometric_orientation combined_face_orientation)
constexpr unsigned int invalid_unsigned_int
constexpr types::subdomain_id invalid_subdomain_id
unsigned char geometric_orientation
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
unsigned int n_hanging_faces_smaller_subdomain
unsigned int n_hanging_faces_larger_subdomain
std::vector< std::pair< CellId, CellId > > shared_faces
std::vector< FaceToCellTopology< 1 > > inner_faces
std::vector< bool > at_processor_boundary
std::vector< FaceToCellTopology< 1 > > boundary_faces
std::vector< FaceToCellTopology< 1 > > refinement_edge_faces
FaceToCellTopology< 1 > create_face(const unsigned int face_no, const typename ::Triangulation< dim >::cell_iterator &cell, const unsigned int number_cell_interior, const typename ::Triangulation< dim >::cell_iterator &neighbor, const unsigned int number_cell_exterior, const bool is_mixed_mesh)
std::vector< FaceCategory > face_is_owned
@ locally_active_done_here
@ multigrid_refinement_edge
@ locally_active_at_boundary
@ locally_active_done_elsewhere
void initialize(const ::Triangulation< dim > &triangulation, const unsigned int mg_level, const bool hold_all_faces_to_owned_cells, const bool build_inner_faces, std::vector< std::pair< unsigned int, unsigned int > > &cell_levels)
std::vector< FaceToCellTopology< 1 > > inner_ghost_faces
void generate_faces(const ::Triangulation< dim > &triangulation, const std::vector< std::pair< unsigned int, unsigned int > > &cell_levels, TaskInfo &task_info)