17#ifndef dealii_face_setup_internal_h
18#define dealii_face_setup_internal_h
41 namespace MatrixFreeFunctions
83 const unsigned int mg_level,
84 const bool hold_all_faces_to_owned_cells,
86 std::vector<std::pair<unsigned int, unsigned int>> &
cell_levels);
97 const std::vector<std::pair<unsigned int, unsigned int>> &
cell_levels,
109 const typename ::Triangulation<dim>::cell_iterator &cell,
111 const typename ::Triangulation<dim>::cell_iterator &neighbor,
113 const bool is_mixed_mesh);
143 template <
int vectorization_w
idth>
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,
171 std::vector<std::pair<unsigned int, unsigned int>> &
cell_levels)
177 if (use_active_cells)
180 typename ::Triangulation<dim>::cell_iterator
dcell(
192 FaceCategory::locally_active_done_elsewhere);
196 std::map<types::subdomain_id, FaceIdentifier>
207 typename ::Triangulation<dim>::cell_iterator
dcell(
209 for (
const unsigned int f :
dcell->face_indices())
211 if (
dcell->at_boundary(f) && !
dcell->has_periodic_neighbor(f))
213 typename ::Triangulation<dim>::cell_iterator neighbor =
214 dcell->neighbor_or_periodic_neighbor(f);
221 if (use_active_cells && neighbor->has_children())
222 for (
unsigned int c = 0;
223 c < (
dcell->has_periodic_neighbor(f) ?
224 dcell->periodic_neighbor(f)
225 ->face(
dcell->periodic_neighbor_face_no(f))
227 dcell->face(f)->n_children());
230 typename ::Triangulation<dim>::cell_iterator
232 dcell->at_boundary(f) ?
233 dcell->periodic_neighbor_child_on_subface(f, c) :
234 dcell->neighbor_child_on_subface(f, c);
239 .n_hanging_faces_larger_subdomain++;
242 .n_hanging_faces_smaller_subdomain++;
247 use_active_cells ? neighbor->subdomain_id() :
248 neighbor->level_subdomain_id();
249 if (neighbor->level() <
dcell->level() &&
254 .n_hanging_faces_smaller_subdomain++;
257 .n_hanging_faces_larger_subdomain++;
259 else if (neighbor->level() ==
dcell->level() &&
295# if defined(DEAL_II_WITH_MPI) && defined(DEBUG)
297 if (const ::parallel::TriangulationBase<dim> *
ptria =
298 dynamic_cast<const ::parallel::TriangulationBase<dim>
367 std::vector<std::tuple<CellId, CellId, unsigned int>>
other_range(
371 std::make_tuple(
inner_face.second.shared_faces[i].second,
388 unsigned int count = 0;
389 for (
unsigned int i = 1;
392 if (
inner_face.second.shared_faces[i].first ==
400 for (
unsigned int k = 0;
k <=
count; ++
k)
421 for (
unsigned int k = 0;
k <=
count; ++
k)
424 .shared_faces[std::get<2>(
428 .shared_faces[std::get<2>(
454 inner_face.second.n_hanging_faces_smaller_subdomain;
456 inner_face.second.n_hanging_faces_larger_subdomain;
459 inner_face.second.n_hanging_faces_smaller_subdomain +
460 inner_face.second.n_hanging_faces_larger_subdomain);
472# if defined(DEAL_II_WITH_MPI) && defined(DEBUG)
532 unsigned int i = 0, c = 0;
549 std::vector<std::pair<CellId, CellId>>
check_faces;
577 std::set<std::pair<unsigned int, unsigned int>> ghost_cells;
580 typename ::Triangulation<dim>::cell_iterator
dcell(
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;
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)
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() ||
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())
644 ->face(
dcell->has_periodic_neighbor(f) ?
645 dcell->periodic_neighbor_face_no(f) :
646 dcell->neighbor_face_no(f))
648 FaceCategory::locally_active_done_here;
654 if (use_active_cells)
656 (
dcell->subdomain_id() != neighbor->subdomain_id());
659 neighbor->level_subdomain_id());
661 else if (hold_all_faces_to_owned_cells ==
true)
664 face_is_owned[
dcell->face(f)->index()] =
665 FaceCategory::ghosted;
666 if (use_active_cells)
668 if (neighbor->has_children())
669 for (
unsigned int s = 0;
670 s <
dcell->face(f)->n_children();
672 if (
dcell->at_boundary(f))
675 ->periodic_neighbor_child_on_subface(f,
678 dcell->subdomain_id())
683 if (
dcell->neighbor_child_on_subface(f, s)
685 dcell->subdomain_id())
690 neighbor->subdomain_id());
694 neighbor->level_subdomain_id());
699 if (use_active_cells && neighbor->has_children())
700 for (
unsigned int s = 0;
701 s <
dcell->face(f)->n_children();
704 typename ::Triangulation<dim>::cell_iterator
706 dcell->at_boundary(f) ?
707 dcell->periodic_neighbor_child_on_subface(f,
709 dcell->neighbor_child_on_subface(f, s);
711 dcell->subdomain_id())
713 std::pair<unsigned int, unsigned int>(
719 std::pair<unsigned int, unsigned int>(
720 neighbor->level(), neighbor->index()));
721 at_processor_boundary[i] =
true;
737 FaceSetup<dim>::generate_faces(
739 const std::vector<std::pair<unsigned int, unsigned int>> &
cell_levels,
740 TaskInfo & task_info)
746 std::map<std::pair<unsigned int, unsigned int>,
unsigned int>
751 typename ::Triangulation<dim>::cell_iterator
dcell(
761 const unsigned int vectorization_length = task_info.vectorization_length;
762 task_info.face_partition_data.resize(
763 task_info.cell_partition_data.size() - 1, 0);
764 task_info.boundary_partition_data.resize(
765 task_info.cell_partition_data.size() - 1, 0);
766 std::vector<unsigned char>
face_visited(face_is_owned.size(), 0);
767 for (
unsigned int partition = 0;
768 partition < task_info.cell_partition_data.size() - 2;
773 for (
unsigned int cell = task_info.cell_partition_data[partition] *
774 vectorization_length;
775 cell < task_info.cell_partition_data[
partition + 1] *
776 vectorization_length;
780 typename ::Triangulation<dim>::cell_iterator
dcell(
784 for (
const auto f :
dcell->face_indices())
787 if (face_is_owned[
dcell->face(f)->index()] ==
788 FaceCategory::locally_active_at_boundary)
793 info.cells_interior[0] = cell;
795 info.interior_face_no = f;
796 info.exterior_face_no =
dcell->face(f)->boundary_id();
799 (
dcell->face(f)->reference_cell() !=
804 info.face_orientation = 0;
805 boundary_faces.push_back(
info);
812 typename ::Triangulation<dim>::cell_iterator
813 neighbor =
dcell->neighbor_or_periodic_neighbor(f);
814 if (use_active_cells && neighbor->has_children())
816 for (
unsigned int c = 0;
817 c <
dcell->face(f)->n_children();
820 typename ::Triangulation<
822 dcell->at_boundary(f) ?
823 dcell->periodic_neighbor_child_on_subface(
825 dcell->neighbor_child_on_subface(f, c);
828 const unsigned int neighbor_face_no =
829 dcell->has_periodic_neighbor(f) ?
830 dcell->periodic_neighbor_face_no(f) :
831 dcell->neighbor_face_no(f);
834 [
dcell->face(f)->child(c)->index()] ==
837 std::pair<unsigned int, unsigned int>
841 [
dcell->face(f)->child(c)->index()] ==
842 FaceCategory::locally_active_done_here)
845 inner_faces.push_back(create_face(
853 else if (face_is_owned[
dcell->face(f)
856 FaceCategory::ghosted)
858 inner_ghost_faces.push_back(create_face(
868 face_is_owned[
dcell->face(f)
871 locally_active_done_elsewhere ||
872 face_is_owned[
dcell->face(f)
874 FaceCategory::ghosted,
880 [
dcell->face(f)->child(c)->index()] = 1;
887 use_active_cells ?
dcell->subdomain_id() :
888 dcell->level_subdomain_id();
890 use_active_cells ? neighbor->subdomain_id() :
891 neighbor->level_subdomain_id();
895 std::pair<unsigned int, unsigned int>
898 if (face_is_owned[
dcell->face(f)->index()] ==
899 FaceCategory::locally_active_done_here)
901 Assert(use_active_cells ||
906 inner_faces.push_back(create_face(
914 else if (face_is_owned[
dcell->face(f)
916 FaceCategory::ghosted)
918 inner_ghost_faces.push_back(create_face(
930 if (
dcell->has_periodic_neighbor(f))
934 dcell->periodic_neighbor_face_no(f))
937 if (face_is_owned[
dcell->face(f)->index()] ==
938 FaceCategory::multigrid_refinement_edge)
940 refinement_edge_faces.push_back(
945 refinement_edge_faces.
size(),
952 task_info.face_partition_data[
partition + 1] =
954 task_info.boundary_partition_data[
partition + 1] =
957 task_info.ghost_face_partition_data.resize(2);
958 task_info.ghost_face_partition_data[0] = 0;
959 task_info.ghost_face_partition_data[1] = inner_ghost_faces.
size();
960 task_info.refinement_edge_face_partition_data.resize(2);
961 task_info.refinement_edge_face_partition_data[0] = 0;
962 task_info.refinement_edge_face_partition_data[1] =
963 refinement_edge_faces.size();
970 FaceSetup<dim>::create_face(
972 const typename ::Triangulation<dim>::cell_iterator &cell,
974 const typename ::Triangulation<dim>::cell_iterator &neighbor,
976 const bool is_mixed_mesh)
982 if (cell->has_periodic_neighbor(
face_no))
983 info.exterior_face_no = cell->periodic_neighbor_face_no(
face_no);
985 info.exterior_face_no = cell->neighbor_face_no(
face_no);
987 info.face_type = is_mixed_mesh ?
988 (cell->face(
face_no)->reference_cell() !=
994 if (cell->level() > neighbor->level())
996 if (cell->has_periodic_neighbor(
face_no))
998 cell->periodic_neighbor_of_coarser_periodic_neighbor(
face_no)
1001 info.subface_index =
1002 cell->neighbor_of_coarser_neighbor(
face_no).second;
1006 if (dim == 3 && cell->has_periodic_neighbor(
face_no))
1009 !cell->get_triangulation()
1010 .get_periodic_face_map()
1013 2 * cell->get_triangulation()
1014 .get_periodic_face_map()
1017 4 * cell->get_triangulation()
1018 .get_periodic_face_map()
1027 info.face_orientation = 0;
1029 !cell->face_orientation(
face_no) + 2 * cell->face_flip(
face_no) +
1030 4 * cell->face_rotation(
face_no);
1032 !neighbor->face_orientation(
info.exterior_face_no) +
1033 2 * neighbor->face_flip(
info.exterior_face_no) +
1034 4 * neighbor->face_rotation(
info.exterior_face_no);
1040 "Face seems to be wrongly oriented from both sides"));
1050 ShapeInfo<double>::compute_orientation_table(2);
1052 {0, 1, 2, 3, 6, 5, 4, 7}};
1053 info.subface_index =
1055 [
info.subface_index];
1073 const std::vector<unsigned int> &active_fe_indices,
1074 const unsigned int length)
1076 if (
face1.interior_face_no !=
face2.interior_face_no)
1078 if (
face1.exterior_face_no !=
face2.exterior_face_no)
1080 if (
face1.subface_index !=
face2.subface_index)
1082 if (
face1.face_orientation !=
face2.face_orientation)
1087 if (active_fe_indices.size() > 0)
1089 if (active_fe_indices[
face1.cells_interior[0] / length] !=
1090 active_fe_indices[
face2.cells_interior[0] / length])
1094 if (active_fe_indices[
face1.cells_exterior[0] / length] !=
1095 active_fe_indices[
face2.cells_exterior[0] / length])
1110 template <
int length>
1113 FaceComparator(
const std::vector<unsigned int> &active_fe_indices)
1114 : active_fe_indices(active_fe_indices)
1128 if (active_fe_indices.size() > 0)
1131 if (active_fe_indices[
face1.cells_interior[0] / length] <
1132 active_fe_indices[
face2.cells_interior[0] / length])
1134 else if (active_fe_indices[
face1.cells_interior[0] / length] >
1135 active_fe_indices[
face2.cells_interior[0] / length])
1141 if (active_fe_indices[
face1.cells_exterior[0] / length] <
1142 active_fe_indices[
face2.cells_exterior[0] / length])
1144 else if (active_fe_indices[
face1.cells_exterior[0] / length] >
1145 active_fe_indices[
face2.cells_exterior[0] / length])
1150 for (
unsigned int i = 0; i < length; ++i)
1151 if (
face1.cells_interior[i] <
face2.cells_interior[i])
1153 else if (
face1.cells_interior[i] >
face2.cells_interior[i])
1155 for (
unsigned int i = 0; i < length; ++i)
1156 if (
face1.cells_exterior[i] <
face2.cells_exterior[i])
1158 else if (
face1.cells_exterior[i] >
face2.cells_exterior[i])
1160 if (
face1.interior_face_no <
face2.interior_face_no)
1162 else if (
face1.interior_face_no >
face2.interior_face_no)
1164 if (
face1.exterior_face_no <
face2.exterior_face_no)
1166 else if (
face1.exterior_face_no >
face2.exterior_face_no)
1179 const std::vector<unsigned int> &active_fe_indices;
1184 template <
int vectorization_w
idth>
1189 std::vector<unsigned int> & face_partition_data,
1191 const std::vector<unsigned int> & active_fe_indices)
1194 std::vector<std::vector<unsigned int>>
faces_type;
1196 unsigned int face_start = face_partition_data[0],
1200 for (
unsigned int partition = 0;
1201 partition < face_partition_data.size() - 1;
1225 face_type.push_back(face);
1237 std::set<FaceToCellTopology<vectorization_width>,
1244 faces_in[face_type[0]].interior_face_no;
1246 faces_in[face_type[0]].exterior_face_no;
1249 faces_in[face_type[0]].face_orientation;
1258 for (
unsigned int f = 0; f <
no_faces; ++f)
1259 if (
faces_in[face_type[f]].cells_interior[0] %
1263 bool is_contiguous =
true;
1265 is_contiguous =
false;
1268 if (
faces_in[face_type[f + v]].cells_interior[0] !=
1269 faces_in[face_type[f]].cells_interior[0] + v)
1270 is_contiguous =
false;
1279 faces_in[face_type[f + v]].cells_interior[0];
1281 faces_in[face_type[f + v]].cells_exterior[0];
1292 for (
unsigned int f = 0; f <
no_faces; ++f)
1299 faces_in[face_type[f]].cells_interior[0];
1301 faces_in[face_type[f]].cells_exterior[0];
1313 partition == face_partition_data.
size() - 2)
1329 for (
unsigned int f = 0; f < v; ++f)
1352 for (
unsigned int i = face_partition_data[0];
1353 i < face_partition_data.back();
1364 for (
unsigned int i = face_partition_data[0];
1365 i < face_partition_data.back();
1367 for (
unsigned int v = 0;
#define DEAL_II_NAMESPACE_OPEN
#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)
const types::subdomain_id invalid_subdomain_id
static const unsigned int invalid_unsigned_int
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)