14#ifndef dealii_face_setup_internal_h
15#define dealii_face_setup_internal_h
39 namespace MatrixFreeFunctions
80 const ::Triangulation<dim> &triangulation,
81 const unsigned int mg_level,
82 const bool hold_all_faces_to_owned_cells,
83 const bool build_inner_faces,
84 std::vector<std::pair<unsigned int, unsigned int>> &cell_levels);
94 const ::Triangulation<dim> &triangulation,
95 const std::vector<std::pair<unsigned int, unsigned int>> &cell_levels,
106 const unsigned int face_no,
107 const typename ::Triangulation<dim>::cell_iterator &cell,
108 const unsigned int number_cell_interior,
109 const typename ::Triangulation<dim>::cell_iterator &neighbor,
110 const unsigned int number_cell_exterior,
111 const bool is_mixed_mesh);
141 template <
int vectorization_w
idth>
145 const std::vector<bool> &hard_vectorization_boundary,
146 std::vector<unsigned int> &face_partition_data,
157 : use_active_cells(true)
164 FaceSetup<dim>::initialize(
165 const ::Triangulation<dim> &triangulation,
166 const unsigned int mg_level,
167 const bool hold_all_faces_to_owned_cells,
168 const bool build_inner_faces,
169 std::vector<std::pair<unsigned int, unsigned int>> &cell_levels)
176 if (use_active_cells)
177 for (
const auto &cell_level : cell_levels)
179 typename ::Triangulation<dim>::cell_iterator dcell(
180 &triangulation, cell_level.first, cell_level.second);
188 at_processor_boundary.resize(cell_levels.size(),
false);
189 face_is_owned.resize(dim > 1 ? triangulation.n_raw_faces() :
190 triangulation.n_vertices(),
191 FaceCategory::locally_active_done_elsewhere);
195 std::map<types::subdomain_id, FaceIdentifier>
196 inner_faces_at_proc_boundary;
197 if (triangulation.locally_owned_subdomain() !=
201 triangulation.locally_owned_subdomain();
202 for (
unsigned int i = 0; i < cell_levels.size(); ++i)
204 if (i > 0 && cell_levels[i] == cell_levels[i - 1])
206 typename ::Triangulation<dim>::cell_iterator dcell(
207 &triangulation, cell_levels[i].
first, cell_levels[i].
second);
208 for (
const unsigned int f : dcell->face_indices())
210 if (dcell->at_boundary(f) && !dcell->has_periodic_neighbor(f))
212 typename ::Triangulation<dim>::cell_iterator neighbor =
213 dcell->neighbor_or_periodic_neighbor(f);
219 const CellId id_mine = dcell->id();
220 if (use_active_cells && neighbor->has_children())
221 for (
unsigned int c = 0;
222 c < (dcell->has_periodic_neighbor(f) ?
223 dcell->periodic_neighbor(f)
224 ->face(dcell->periodic_neighbor_face_no(f))
226 dcell->face(f)->n_children());
229 typename ::Triangulation<dim>::cell_iterator
231 dcell->at_boundary(f) ?
232 dcell->periodic_neighbor_child_on_subface(f, c) :
233 dcell->neighbor_child_on_subface(f, c);
235 neighbor_c->subdomain_id();
236 if (my_domain < neigh_domain)
237 inner_faces_at_proc_boundary[neigh_domain]
238 .n_hanging_faces_larger_subdomain++;
239 else if (my_domain > neigh_domain)
240 inner_faces_at_proc_boundary[neigh_domain]
241 .n_hanging_faces_smaller_subdomain++;
246 use_active_cells ? neighbor->subdomain_id() :
247 neighbor->level_subdomain_id();
248 if (neighbor->level() < dcell->level() &&
251 if (my_domain < neigh_domain)
252 inner_faces_at_proc_boundary[neigh_domain]
253 .n_hanging_faces_smaller_subdomain++;
254 else if (my_domain > neigh_domain)
255 inner_faces_at_proc_boundary[neigh_domain]
256 .n_hanging_faces_larger_subdomain++;
258 else if (neighbor->level() == dcell->level() &&
259 my_domain != neigh_domain)
265 const CellId id_neigh = neighbor->id();
266 if (my_domain < neigh_domain)
267 inner_faces_at_proc_boundary[neigh_domain]
268 .shared_faces.emplace_back(id_mine, id_neigh);
270 inner_faces_at_proc_boundary[neigh_domain]
271 .shared_faces.emplace_back(id_neigh, id_mine);
280 for (
auto &inner_face : inner_faces_at_proc_boundary)
282 Assert(inner_face.first != my_domain,
284 std::sort(inner_face.second.shared_faces.begin(),
285 inner_face.second.shared_faces.end());
286 inner_face.second.shared_faces.erase(
287 std::unique(inner_face.second.shared_faces.begin(),
288 inner_face.second.shared_faces.end()),
289 inner_face.second.shared_faces.end());
294# if defined(DEAL_II_WITH_MPI) && defined(DEBUG)
296 if (const ::parallel::TriangulationBase<dim> *ptria =
297 dynamic_cast<const ::parallel::TriangulationBase<dim>
299 comm = ptria->get_mpi_communicator();
302 unsigned int mysize = inner_face.second.shared_faces.size();
305 int ierr = MPI_Sendrecv(&mysize,
314 600 + inner_face.first,
319 mysize = inner_face.second.n_hanging_faces_smaller_subdomain;
320 ierr = MPI_Sendrecv(&mysize,
329 700 + inner_face.first,
334 mysize = inner_face.second.n_hanging_faces_larger_subdomain;
335 ierr = MPI_Sendrecv(&mysize,
344 800 + inner_face.first,
366 std::vector<std::tuple<CellId, CellId, unsigned int>> other_range(
367 inner_face.second.shared_faces.size());
368 for (
unsigned int i = 0; i < other_range.size(); ++i)
370 std::make_tuple(inner_face.second.shared_faces[i].second,
371 inner_face.second.shared_faces[i].first,
373 std::sort(other_range.begin(), other_range.end());
381 unsigned int n_faces_lower_proc = 0, n_faces_higher_proc = 0;
382 std::vector<signed char> assignment(other_range.size(), 0);
383 if (inner_face.second.shared_faces.size() > 0)
387 unsigned int count = 0;
388 for (
unsigned int i = 1;
389 i < inner_face.second.shared_faces.size();
391 if (inner_face.second.shared_faces[i].first ==
392 inner_face.second.shared_faces[i - 1 - count].first)
399 for (
unsigned int k = 0; k <= count; ++k)
400 assignment[i - 1 - k] = 1;
401 n_faces_higher_proc += count + 1;
411 for (
unsigned int i = 1; i < other_range.size(); ++i)
412 if (std::get<0>(other_range[i]) ==
413 std::get<0>(other_range[i - 1 - count]))
420 for (
unsigned int k = 0; k <= count; ++k)
423 .shared_faces[std::get<2>(
427 .shared_faces[std::get<2>(
428 other_range[i - 1 - k])]
433 if (assignment[std::get<2>(
434 other_range[i - 1 - k])] == 0)
436 assignment[std::get<2>(
437 other_range[i - 1 - k])] = -1;
438 ++n_faces_lower_proc;
452 n_faces_lower_proc +=
453 inner_face.second.n_hanging_faces_smaller_subdomain;
454 n_faces_higher_proc +=
455 inner_face.second.n_hanging_faces_larger_subdomain;
456 const unsigned int n_total_faces_at_proc_boundary =
457 (inner_face.second.shared_faces.size() +
458 inner_face.second.n_hanging_faces_smaller_subdomain +
459 inner_face.second.n_hanging_faces_larger_subdomain);
460 unsigned int split_index = n_total_faces_at_proc_boundary / 2;
461 if (split_index < n_faces_lower_proc)
463 else if (split_index <
464 n_total_faces_at_proc_boundary - n_faces_higher_proc)
465 split_index -= n_faces_lower_proc;
467 split_index = n_total_faces_at_proc_boundary -
468 n_faces_higher_proc - n_faces_lower_proc;
471# if defined(DEAL_II_WITH_MPI) && defined(DEBUG)
472 ierr = MPI_Sendrecv(&split_index,
481 900 + inner_face.first,
486 ierr = MPI_Sendrecv(&n_faces_lower_proc,
495 1000 + inner_face.first,
500 ierr = MPI_Sendrecv(&n_faces_higher_proc,
509 1100 + inner_face.first,
517 std::vector<std::pair<CellId, CellId>> owned_faces_lower,
519 for (
unsigned int i = 0; i < assignment.size(); ++i)
520 if (assignment[i] < 0)
521 owned_faces_lower.push_back(
522 inner_face.second.shared_faces[i]);
523 else if (assignment[i] > 0)
524 owned_faces_higher.push_back(
525 inner_face.second.shared_faces[i]);
527 inner_face.second.shared_faces.size() + 1 -
528 owned_faces_lower.size() -
529 owned_faces_higher.size());
531 unsigned int i = 0, c = 0;
532 for (; i < assignment.size() && c < split_index; ++i)
533 if (assignment[i] == 0)
535 owned_faces_lower.push_back(
536 inner_face.second.shared_faces[i]);
539 for (; i < assignment.size(); ++i)
540 if (assignment[i] == 0)
542 owned_faces_higher.push_back(
543 inner_face.second.shared_faces[i]);
549 std::vector<std::pair<CellId, CellId>> check_faces;
550 check_faces.insert(check_faces.end(),
551 owned_faces_lower.begin(),
552 owned_faces_lower.end());
553 check_faces.insert(check_faces.end(),
554 owned_faces_higher.begin(),
555 owned_faces_higher.end());
556 std::sort(check_faces.begin(), check_faces.end());
558 inner_face.second.shared_faces.size());
559 for (
unsigned int i = 0; i < check_faces.size(); ++i)
560 Assert(check_faces[i] == inner_face.second.shared_faces[i],
565 if (my_domain < inner_face.first)
566 inner_face.second.shared_faces.swap(owned_faces_lower);
568 inner_face.second.shared_faces.swap(owned_faces_higher);
570 std::sort(inner_face.second.shared_faces.begin(),
571 inner_face.second.shared_faces.end());
577 std::set<std::pair<unsigned int, unsigned int>> ghost_cells;
578 for (
unsigned int i = 0; i < cell_levels.size(); ++i)
580 typename ::Triangulation<dim>::cell_iterator dcell(
581 &triangulation, cell_levels[i].
first, cell_levels[i].
second);
582 if (use_active_cells)
584 for (
const auto f : dcell->face_indices())
586 if (dcell->at_boundary(f) && !dcell->has_periodic_neighbor(f))
587 face_is_owned[dcell->face(f)->index()] =
588 FaceCategory::locally_active_at_boundary;
589 else if (!build_inner_faces)
594 else if ((dcell->at_boundary(f) ==
false ||
595 dcell->has_periodic_neighbor(f)) &&
597 dcell->neighbor_or_periodic_neighbor(f)->level() <
600 face_is_owned[dcell->face(f)->index()] =
601 FaceCategory::multigrid_refinement_edge;
605 typename ::Triangulation<dim>::cell_iterator neighbor =
606 dcell->neighbor_or_periodic_neighbor(f);
609 if (use_active_cells && neighbor->has_children() &&
610 hold_all_faces_to_owned_cells ==
false)
613 bool add_to_ghost =
false;
615 id1 = use_active_cells ? dcell->subdomain_id() :
616 dcell->level_subdomain_id(),
617 id2 = use_active_cells ?
618 (neighbor->has_children() ?
619 dcell->neighbor_child_on_subface(f, 0)
621 neighbor->subdomain_id()) :
622 neighbor->level_subdomain_id();
630 (use_active_cells ==
false || neighbor->is_active())) ||
631 dcell->level() > neighbor->level() ||
633 inner_faces_at_proc_boundary[id2].shared_faces.begin(),
634 inner_faces_at_proc_boundary[id2].shared_faces.end(),
635 std::make_pair(id1 < id2 ? dcell->
id() : neighbor->id(),
636 id1 < id2 ? neighbor->id() :
639 face_is_owned[dcell->face(f)->index()] =
640 FaceCategory::locally_active_done_here;
641 if (dcell->level() == neighbor->level() ||
642 dcell->has_periodic_neighbor(f))
645 ->face(dcell->has_periodic_neighbor(f) ?
646 dcell->periodic_neighbor_face_no(f) :
647 dcell->neighbor_face_no(f))
649 FaceCategory::locally_active_done_here;
655 if (use_active_cells)
657 (dcell->subdomain_id() != neighbor->subdomain_id());
659 add_to_ghost = (dcell->level_subdomain_id() !=
660 neighbor->level_subdomain_id());
662 else if (hold_all_faces_to_owned_cells ==
true)
665 face_is_owned[dcell->face(f)->index()] =
666 FaceCategory::ghosted;
667 if (use_active_cells)
669 if (neighbor->has_children())
670 for (
unsigned int s = 0;
671 s < dcell->face(f)->n_children();
673 if (dcell->at_boundary(f))
676 ->periodic_neighbor_child_on_subface(f,
679 dcell->subdomain_id())
684 if (dcell->neighbor_child_on_subface(f, s)
686 dcell->subdomain_id())
690 add_to_ghost = (dcell->subdomain_id() !=
691 neighbor->subdomain_id());
694 add_to_ghost = (dcell->level_subdomain_id() !=
695 neighbor->level_subdomain_id());
700 if (use_active_cells && neighbor->has_children())
701 for (
unsigned int s = 0;
702 s < dcell->face(f)->n_children();
705 typename ::Triangulation<dim>::cell_iterator
707 dcell->at_boundary(f) ?
708 dcell->periodic_neighbor_child_on_subface(f,
710 dcell->neighbor_child_on_subface(f, s);
711 if (neighbor_child->subdomain_id() !=
712 dcell->subdomain_id())
714 std::pair<unsigned int, unsigned int>(
715 neighbor_child->level(),
716 neighbor_child->index()));
720 std::pair<unsigned int, unsigned int>(
721 neighbor->level(), neighbor->index()));
722 at_processor_boundary[i] =
true;
730 for (
const auto &ghost_cell : ghost_cells)
738 FaceSetup<dim>::generate_faces(
739 const ::Triangulation<dim> &triangulation,
740 const std::vector<std::pair<unsigned int, unsigned int>> &cell_levels,
743 const bool is_mixed_mesh = triangulation.is_mixed_mesh();
747 std::map<std::pair<unsigned int, unsigned int>,
unsigned int>
749 for (
unsigned int cell = 0; cell < cell_levels.size(); ++cell)
750 if (cell == 0 || cell_levels[cell] != cell_levels[cell - 1])
752 typename ::Triangulation<dim>::cell_iterator dcell(
754 cell_levels[cell].
first,
755 cell_levels[cell].
second);
756 std::pair<unsigned int, unsigned int> level_index(dcell->level(),
758 map_to_vectorized[level_index] = cell;
762 const unsigned int vectorization_length = task_info.vectorization_length;
763 task_info.face_partition_data.resize(
764 task_info.cell_partition_data.size() - 1, 0);
765 task_info.boundary_partition_data.resize(
766 task_info.cell_partition_data.size() - 1, 0);
767 std::vector<unsigned char> face_visited(face_is_owned.size(), 0);
768 for (
unsigned int partition = 0;
769 partition < task_info.cell_partition_data.size() - 2;
772 unsigned int boundary_counter = 0;
773 unsigned int inner_counter = 0;
774 for (
unsigned int cell = task_info.cell_partition_data[partition] *
775 vectorization_length;
776 cell < task_info.cell_partition_data[
partition + 1] *
777 vectorization_length;
779 if (cell == 0 || cell_levels[cell] != cell_levels[cell - 1])
781 typename ::Triangulation<dim>::cell_iterator dcell(
783 cell_levels[cell].
first,
784 cell_levels[cell].
second);
785 for (
const auto f : dcell->face_indices())
788 if (face_is_owned[dcell->face(f)->index()] ==
789 FaceCategory::locally_active_at_boundary)
793 FaceToCellTopology<1> info;
794 info.cells_interior[0] = cell;
796 info.interior_face_no = f;
797 info.exterior_face_no = dcell->face(f)->boundary_id();
800 (dcell->face(f)->reference_cell() !=
805 info.face_orientation = 0;
806 boundary_faces.push_back(info);
808 face_visited[dcell->face(f)->index()]++;
813 typename ::Triangulation<dim>::cell_iterator
814 neighbor = dcell->neighbor_or_periodic_neighbor(f);
815 if (use_active_cells && neighbor->has_children())
817 const unsigned int n_children =
818 (dim == 1) ? 1 : dcell->face(f)->n_children();
819 for (
unsigned int c = 0; c < n_children; ++c)
821 typename ::Triangulation<
822 dim>::cell_iterator neighbor_c;
825 (dcell->at_boundary(f) ?
827 ->periodic_neighbor_child_on_subface(
829 dcell->neighbor_child_on_subface(f, c));
834 neighbor_c = neighbor->child(1 - f);
835 while (!neighbor_c->is_active())
836 neighbor_c = neighbor_c->child(1 - f);
839 neighbor_c->subdomain_id();
840 const unsigned int neighbor_face_no =
841 dcell->has_periodic_neighbor(f) ?
842 dcell->periodic_neighbor_face_no(f) :
843 dcell->neighbor_face_no(f);
844 const unsigned int child_face_index =
845 dim > 1 ? dcell->face(f)->child(c)->index() :
846 dcell->face(f)->index();
847 if (neigh_domain != dcell->subdomain_id() ||
848 face_visited[child_face_index] == 1)
850 std::pair<unsigned int, unsigned int>
851 level_index(neighbor_c->level(),
852 neighbor_c->index());
853 if (face_is_owned[child_face_index] ==
854 FaceCategory::locally_active_done_here)
857 inner_faces.push_back(create_face(
860 map_to_vectorized[level_index],
865 else if (face_is_owned[child_face_index] ==
866 FaceCategory::ghosted ||
867 face_is_owned[dcell->face(f)
869 FaceCategory::ghosted)
871 inner_ghost_faces.push_back(create_face(
874 map_to_vectorized[level_index],
880 Assert(face_is_owned[dcell->face(f)
883 locally_active_done_elsewhere,
888 face_visited[child_face_index] = 1;
889 if (dcell->has_periodic_neighbor(f))
890 face_visited[dcell->face(f)->index()] = 1;
897 use_active_cells ? dcell->subdomain_id() :
898 dcell->level_subdomain_id();
900 use_active_cells ? neighbor->subdomain_id() :
901 neighbor->level_subdomain_id();
902 const unsigned int face_index =
903 dcell->face(f)->index();
904 if (neigh_domain != my_domain ||
905 face_visited[face_index] == 1)
907 std::pair<unsigned int, unsigned int>
908 level_index(neighbor->level(),
910 if (face_is_owned[dcell->face(f)->index()] ==
911 FaceCategory::locally_active_done_here)
913 Assert(use_active_cells ||
918 inner_faces.push_back(create_face(
923 map_to_vectorized[level_index],
926 else if (face_is_owned[face_index] ==
927 FaceCategory::ghosted)
929 inner_ghost_faces.push_back(create_face(
934 map_to_vectorized[level_index],
940 face_visited[face_index] = 1;
941 if (dcell->has_periodic_neighbor(f))
945 dcell->periodic_neighbor_face_no(f))
948 if (face_is_owned[face_index] ==
949 FaceCategory::multigrid_refinement_edge)
951 refinement_edge_faces.push_back(
956 refinement_edge_faces.size(),
963 task_info.face_partition_data[
partition + 1] =
964 task_info.face_partition_data[
partition] + inner_counter;
965 task_info.boundary_partition_data[
partition + 1] =
966 task_info.boundary_partition_data[
partition] + boundary_counter;
968 task_info.ghost_face_partition_data.resize(2);
969 task_info.ghost_face_partition_data[0] = 0;
970 task_info.ghost_face_partition_data[1] = inner_ghost_faces.size();
971 task_info.refinement_edge_face_partition_data.resize(2);
972 task_info.refinement_edge_face_partition_data[0] = 0;
973 task_info.refinement_edge_face_partition_data[1] =
974 refinement_edge_faces.size();
980 FaceToCellTopology<1>
981 FaceSetup<dim>::create_face(
982 const unsigned int face_no,
983 const typename ::Triangulation<dim>::cell_iterator &cell,
984 const unsigned int number_cell_interior,
985 const typename ::Triangulation<dim>::cell_iterator &neighbor,
986 const unsigned int number_cell_exterior,
987 const bool is_mixed_mesh)
989 FaceToCellTopology<1> info;
990 info.cells_interior[0] = number_cell_interior;
991 info.cells_exterior[0] = number_cell_exterior;
992 info.interior_face_no = face_no;
993 if (cell->has_periodic_neighbor(face_no))
994 info.exterior_face_no = cell->periodic_neighbor_face_no(face_no);
996 info.exterior_face_no = cell->neighbor_face_no(face_no);
998 info.face_type = is_mixed_mesh ?
999 (cell->face(face_no)->reference_cell() !=
1007 if (dim > 1 && cell->level() > neighbor->level())
1009 if (cell->has_periodic_neighbor(face_no))
1010 info.subface_index =
1011 cell->periodic_neighbor_of_coarser_periodic_neighbor(face_no)
1014 info.subface_index =
1015 cell->neighbor_of_coarser_neighbor(face_no).second;
1019 if (cell->has_periodic_neighbor(face_no))
1021 info.face_orientation = cell->get_triangulation()
1022 .get_periodic_face_map()
1023 .at({cell, face_no})
1028 const auto interior_face_orientation =
1029 cell->combined_face_orientation(face_no);
1030 const auto exterior_face_orientation =
1031 neighbor->combined_face_orientation(info.exterior_face_no);
1032 if (interior_face_orientation !=
1035 info.face_orientation = 8 + interior_face_orientation;
1036 Assert(exterior_face_orientation ==
1039 "Face seems to be wrongly oriented from both sides"));
1042 info.face_orientation = exterior_face_orientation;
1046 if (cell->level() > neighbor->level() &&
1047 exterior_face_orientation > 0)
1050 ShapeInfo<double>::compute_orientation_table(2);
1051 const auto face_reference_cell =
1052 cell->face(face_no)->reference_cell();
1053 info.subface_index = orientation(
1054 face_reference_cell.get_inverse_combined_orientation(
1055 exterior_face_orientation),
1056 info.subface_index);
1072 compare_faces_for_vectorization(
1073 const FaceToCellTopology<1> &face1,
1074 const FaceToCellTopology<1> &face2,
1075 const std::vector<unsigned int> &active_fe_indices,
1076 const unsigned int length)
1078 if (face1.interior_face_no != face2.interior_face_no)
1080 if (face1.exterior_face_no != face2.exterior_face_no)
1082 if (face1.subface_index != face2.subface_index)
1084 if (face1.face_orientation != face2.face_orientation)
1086 if (face1.face_type != face2.face_type)
1089 if (active_fe_indices.size() > 0)
1091 if (active_fe_indices[face1.cells_interior[0] / length] !=
1092 active_fe_indices[face2.cells_interior[0] / length])
1096 if (active_fe_indices[face1.cells_exterior[0] / length] !=
1097 active_fe_indices[face2.cells_exterior[0] / length])
1112 template <
int length>
1113 struct FaceComparator
1115 FaceComparator(
const std::vector<unsigned int> &active_fe_indices)
1116 : active_fe_indices(active_fe_indices)
1120 operator()(
const FaceToCellTopology<length> &face1,
1121 const FaceToCellTopology<length> &face2)
const
1124 if (face1.face_type < face2.face_type)
1126 else if (face1.face_type > face2.face_type)
1130 if (active_fe_indices.size() > 0)
1133 if (active_fe_indices[face1.cells_interior[0] / length] <
1134 active_fe_indices[face2.cells_interior[0] / length])
1136 else if (active_fe_indices[face1.cells_interior[0] / length] >
1137 active_fe_indices[face2.cells_interior[0] / length])
1143 if (active_fe_indices[face1.cells_exterior[0] / length] <
1144 active_fe_indices[face2.cells_exterior[0] / length])
1146 else if (active_fe_indices[face1.cells_exterior[0] / length] >
1147 active_fe_indices[face2.cells_exterior[0] / length])
1152 for (
unsigned int i = 0; i < length; ++i)
1153 if (face1.cells_interior[i] < face2.cells_interior[i])
1155 else if (face1.cells_interior[i] > face2.cells_interior[i])
1157 for (
unsigned int i = 0; i < length; ++i)
1158 if (face1.cells_exterior[i] < face2.cells_exterior[i])
1160 else if (face1.cells_exterior[i] > face2.cells_exterior[i])
1162 if (face1.interior_face_no < face2.interior_face_no)
1164 else if (face1.interior_face_no > face2.interior_face_no)
1166 if (face1.exterior_face_no < face2.exterior_face_no)
1168 else if (face1.exterior_face_no > face2.exterior_face_no)
1181 const std::vector<unsigned int> &active_fe_indices;
1186 template <
int vectorization_w
idth>
1189 const std::vector<FaceToCellTopology<1>> &faces_in,
1190 const std::vector<bool> &hard_vectorization_boundary,
1191 std::vector<unsigned int> &face_partition_data,
1192 std::vector<FaceToCellTopology<vectorization_width>> &faces_out,
1193 const std::vector<unsigned int> &active_fe_indices)
1195 FaceToCellTopology<vectorization_width> face_batch;
1196 std::vector<std::vector<unsigned int>> faces_type;
1198 unsigned int face_start = face_partition_data[0],
1199 face_end = face_partition_data[0];
1201 face_partition_data[0] = faces_out.size();
1202 for (
unsigned int partition = 0;
1203 partition < face_partition_data.size() - 1;
1206 std::vector<std::vector<unsigned int>> new_faces_type;
1209 face_start = face_end;
1210 face_end = face_partition_data[
partition + 1];
1217 for (
unsigned int face = face_start; face < face_end; ++face)
1219 for (
auto &face_type : faces_type)
1222 if (compare_faces_for_vectorization(faces_in[face],
1223 faces_in[face_type[0]],
1225 vectorization_width))
1227 face_type.push_back(face);
1231 faces_type.emplace_back(1, face);
1237 FaceComparator<vectorization_width> face_comparator(
1239 std::set<FaceToCellTopology<vectorization_width>,
1240 FaceComparator<vectorization_width>>
1241 new_faces(face_comparator);
1242 for (
const auto &face_type : faces_type)
1244 face_batch.face_type = faces_in[face_type[0]].face_type;
1245 face_batch.interior_face_no =
1246 faces_in[face_type[0]].interior_face_no;
1247 face_batch.exterior_face_no =
1248 faces_in[face_type[0]].exterior_face_no;
1249 face_batch.subface_index = faces_in[face_type[0]].subface_index;
1250 face_batch.face_orientation =
1251 faces_in[face_type[0]].face_orientation;
1252 unsigned int no_faces = face_type.size();
1253 std::vector<unsigned char> touched(no_faces, 0);
1259 unsigned int n_vectorized = 0;
1260 for (
unsigned int f = 0; f < no_faces; ++f)
1261 if (faces_in[face_type[f]].cells_interior[0] %
1262 vectorization_width ==
1265 bool is_contiguous =
true;
1266 if (f + vectorization_width > no_faces)
1267 is_contiguous =
false;
1269 for (
unsigned int v = 1; v < vectorization_width; ++v)
1270 if (faces_in[face_type[f + v]].cells_interior[0] !=
1271 faces_in[face_type[f]].cells_interior[0] + v)
1272 is_contiguous =
false;
1277 vectorization_width + 1);
1278 for (
unsigned int v = 0; v < vectorization_width; ++v)
1280 face_batch.cells_interior[v] =
1281 faces_in[face_type[f + v]].cells_interior[0];
1282 face_batch.cells_exterior[v] =
1283 faces_in[face_type[f + v]].cells_exterior[0];
1286 new_faces.insert(face_batch);
1287 f += vectorization_width - 1;
1288 n_vectorized += vectorization_width;
1292 std::vector<unsigned int> untouched;
1293 untouched.reserve(no_faces - n_vectorized);
1294 for (
unsigned int f = 0; f < no_faces; ++f)
1295 if (touched[f] == 0)
1296 untouched.push_back(f);
1298 for (
const auto f : untouched)
1300 face_batch.cells_interior[v] =
1301 faces_in[face_type[f]].cells_interior[0];
1302 face_batch.cells_exterior[v] =
1303 faces_in[face_type[f]].cells_exterior[0];
1305 if (v == vectorization_width)
1307 new_faces.insert(face_batch);
1311 if (v > 0 && v < vectorization_width)
1314 if (hard_vectorization_boundary[partition + 1] ||
1315 partition == face_partition_data.size() - 2)
1317 for (; v < vectorization_width; ++v)
1320 face_batch.cells_interior[v] =
1322 face_batch.cells_exterior[v] =
1325 new_faces.insert(face_batch);
1330 std::vector<unsigned int> untreated(v);
1331 for (
unsigned int f = 0; f < v; ++f)
1332 untreated[f] = face_type[*(untouched.end() - 1 - f)];
1333 new_faces_type.push_back(untreated);
1339 for (
auto it = new_faces.begin(); it != new_faces.end(); ++it)
1340 faces_out.push_back(*it);
1341 face_partition_data[
partition + 1] += new_faces.size();
1344 faces_type = std::move(new_faces_type);
1350 for (
const auto &face_type : faces_type)
1354 unsigned int nfaces = 0;
1355 for (
unsigned int i = face_partition_data[0];
1356 i < face_partition_data.back();
1358 for (
unsigned int v = 0; v < vectorization_width; ++v)
1359 nfaces += (faces_out[i].cells_interior[v] !=
1363 std::vector<std::pair<unsigned int, unsigned int>> in_faces,
1365 for (
const auto &face_in : faces_in)
1366 in_faces.emplace_back(face_in.cells_interior[0],
1367 face_in.cells_exterior[0]);
1368 for (
unsigned int i = face_partition_data[0];
1369 i < face_partition_data.back();
1371 for (
unsigned int v = 0;
1372 v < vectorization_width && faces_out[i].cells_interior[v] !=
1375 out_faces.emplace_back(faces_out[i].cells_interior[v],
1376 faces_out[i].cells_exterior[v]);
1377 std::sort(in_faces.begin(), in_faces.end());
1378 std::sort(out_faces.begin(), out_faces.end());
1380 for (
unsigned int i = 0; i < in_faces.size(); ++i)
* * Point< dim > operator()(const Point< dim > &p) const *
***mech_lbc_system increment_interpolation_handlers push_back(scale_z_handler)
#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< dim > & 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)
constexpr unsigned int invalid_unsigned_int
constexpr types::subdomain_id invalid_subdomain_id
constexpr types::geometric_orientation default_geometric_orientation
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)