Reference documentation for deal.II version 9.3.3
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
face_setup_internal.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2018 - 2021 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
17#ifndef dealii_face_setup_internal_h
18#define dealii_face_setup_internal_h
19
20#include <deal.II/base/config.h>
21
24
26
27#include <deal.II/grid/tria.h>
29
32
33#include <fstream>
34
35
37
38
39namespace internal
40{
41 namespace MatrixFreeFunctions
42 {
48 {
52 {}
53
54 std::vector<std::pair<CellId, CellId>> shared_faces;
57 };
58
59
60
69 template <int dim>
70 struct FaceSetup
71 {
73
80 void
82 const ::Triangulation<dim> &triangulation,
83 const unsigned int mg_level,
84 const bool hold_all_faces_to_owned_cells,
85 std::vector<std::pair<unsigned int, unsigned int>> &cell_levels);
86
93 void
95 const ::Triangulation<dim> & triangulation,
96 const std::vector<std::pair<unsigned int, unsigned int>> &cell_levels,
97 TaskInfo & task_info);
98
107 const unsigned int face_no,
108 const typename ::Triangulation<dim>::cell_iterator &cell,
109 const unsigned int number_cell_interior,
110 const typename ::Triangulation<dim>::cell_iterator &neighbor,
111 const unsigned int number_cell_exterior);
112
114
119 enum class FaceCategory : char
120 {
121 locally_active_at_boundary,
122 locally_active_done_here,
123 locally_active_done_elsewhere,
124 ghosted,
125 multigrid_refinement_edge
126 };
127
128 std::vector<FaceCategory> face_is_owned;
129 std::vector<bool> at_processor_boundary;
130 std::vector<FaceToCellTopology<1>> inner_faces;
131 std::vector<FaceToCellTopology<1>> boundary_faces;
132 std::vector<FaceToCellTopology<1>> inner_ghost_faces;
133 std::vector<FaceToCellTopology<1>> refinement_edge_faces;
134 };
135
136
137
141 template <int vectorization_width>
142 void
144 const std::vector<FaceToCellTopology<1>> &faces_in,
145 const std::vector<bool> & hard_vectorization_boundary,
146 std::vector<unsigned int> & face_partition_data,
147 std::vector<FaceToCellTopology<vectorization_width>> &faces_out);
148
149
150
151 /* -------------------------------------------------------------------- */
152
153#ifndef DOXYGEN
154
155 template <int dim>
157 : use_active_cells(true)
158 {}
159
160
161
162 template <int dim>
163 void
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 std::vector<std::pair<unsigned int, unsigned int>> &cell_levels)
169 {
170 use_active_cells = mg_level == numbers::invalid_unsigned_int;
171
172# ifdef DEBUG
173 // safety check
174 if (use_active_cells)
175 for (const auto &cell_level : cell_levels)
176 {
177 typename ::Triangulation<dim>::cell_iterator dcell(
178 &triangulation, cell_level.first, cell_level.second);
179 Assert(dcell->is_active(), ExcInternalError());
180 }
181# endif
182
183 // step 1: add ghost cells for those cells that we identify as
184 // interesting
185
186 at_processor_boundary.resize(cell_levels.size(), false);
187 face_is_owned.resize(dim > 1 ? triangulation.n_raw_faces() :
188 triangulation.n_vertices(),
189 FaceCategory::locally_active_done_elsewhere);
190
191 // go through the mesh and divide the faces on the processor
192 // boundaries as evenly as possible between the processors
193 std::map<types::subdomain_id, FaceIdentifier>
194 inner_faces_at_proc_boundary;
195 if (triangulation.locally_owned_subdomain() !=
197 {
198 const types::subdomain_id my_domain =
199 triangulation.locally_owned_subdomain();
200 for (unsigned int i = 0; i < cell_levels.size(); ++i)
201 {
202 if (i > 0 && cell_levels[i] == cell_levels[i - 1])
203 continue;
204 typename ::Triangulation<dim>::cell_iterator dcell(
205 &triangulation, cell_levels[i].first, cell_levels[i].second);
206 for (const unsigned int f : dcell->face_indices())
207 {
208 if (dcell->at_boundary(f) && !dcell->has_periodic_neighbor(f))
209 continue;
210 typename ::Triangulation<dim>::cell_iterator neighbor =
211 dcell->neighbor_or_periodic_neighbor(f);
212
213 // faces at hanging nodes are always treated by the processor
214 // who owns the element on the fine side. but we need to count
215 // the number of inner faces in order to balance the remaining
216 // faces properly
217 const CellId id_mine = dcell->id();
218 if (use_active_cells && neighbor->has_children())
219 for (unsigned int c = 0;
220 c < (dcell->has_periodic_neighbor(f) ?
221 dcell->periodic_neighbor(f)
222 ->face(dcell->periodic_neighbor_face_no(f))
223 ->n_children() :
224 dcell->face(f)->n_children());
225 ++c)
226 {
227 typename ::Triangulation<dim>::cell_iterator
228 neighbor_c =
229 dcell->at_boundary(f) ?
230 dcell->periodic_neighbor_child_on_subface(f, c) :
231 dcell->neighbor_child_on_subface(f, c);
232 const types::subdomain_id neigh_domain =
233 neighbor_c->subdomain_id();
234 if (my_domain < neigh_domain)
235 inner_faces_at_proc_boundary[neigh_domain]
236 .n_hanging_faces_larger_subdomain++;
237 else if (my_domain > neigh_domain)
238 inner_faces_at_proc_boundary[neigh_domain]
239 .n_hanging_faces_smaller_subdomain++;
240 }
241 else
242 {
243 const types::subdomain_id neigh_domain =
244 use_active_cells ? neighbor->subdomain_id() :
245 neighbor->level_subdomain_id();
246 if (neighbor->level() < dcell->level() &&
247 use_active_cells)
248 {
249 if (my_domain < neigh_domain)
250 inner_faces_at_proc_boundary[neigh_domain]
251 .n_hanging_faces_smaller_subdomain++;
252 else if (my_domain > neigh_domain)
253 inner_faces_at_proc_boundary[neigh_domain]
254 .n_hanging_faces_larger_subdomain++;
255 }
256 else if (neighbor->level() == dcell->level() &&
257 my_domain != neigh_domain)
258 {
259 // always list the cell whose owner has the lower
260 // subdomain id first. this applies to both processors
261 // involved, so both processors will generate the same
262 // list that we will later order
263 const CellId id_neigh = neighbor->id();
264 if (my_domain < neigh_domain)
265 inner_faces_at_proc_boundary[neigh_domain]
266 .shared_faces.emplace_back(id_mine, id_neigh);
267 else
268 inner_faces_at_proc_boundary[neigh_domain]
269 .shared_faces.emplace_back(id_neigh, id_mine);
270 }
271 }
272 }
273 }
274
275 // sort the cell ids related to each neighboring processor. This
276 // algorithm is symmetric so every processor combination should
277 // arrive here and no deadlock should be possible
278 for (auto &inner_face : inner_faces_at_proc_boundary)
279 {
280 Assert(inner_face.first != my_domain,
281 ExcInternalError("Should not send info to myself"));
282 std::sort(inner_face.second.shared_faces.begin(),
283 inner_face.second.shared_faces.end());
284 inner_face.second.shared_faces.erase(
285 std::unique(inner_face.second.shared_faces.begin(),
286 inner_face.second.shared_faces.end()),
287 inner_face.second.shared_faces.end());
288
289 // safety check: both involved processors should see the same list
290 // because the pattern of ghosting is symmetric. We test this by
291 // looking at the length of the lists of faces
292# if defined(DEAL_II_WITH_MPI) && defined(DEBUG)
293 MPI_Comm comm = MPI_COMM_SELF;
294 if (const ::parallel::TriangulationBase<dim> *ptria =
295 dynamic_cast<const ::parallel::TriangulationBase<dim>
296 *>(&triangulation))
297 comm = ptria->get_communicator();
298
299 MPI_Status status;
300 unsigned int mysize = inner_face.second.shared_faces.size();
301 unsigned int othersize = numbers::invalid_unsigned_int;
302 MPI_Sendrecv(&mysize,
303 1,
304 MPI_UNSIGNED,
305 inner_face.first,
306 600 + my_domain,
307 &othersize,
308 1,
309 MPI_UNSIGNED,
310 inner_face.first,
311 600 + inner_face.first,
312 comm,
313 &status);
314 AssertDimension(mysize, othersize);
315 mysize = inner_face.second.n_hanging_faces_smaller_subdomain;
316 MPI_Sendrecv(&mysize,
317 1,
318 MPI_UNSIGNED,
319 inner_face.first,
320 700 + my_domain,
321 &othersize,
322 1,
323 MPI_UNSIGNED,
324 inner_face.first,
325 700 + inner_face.first,
326 comm,
327 &status);
328 AssertDimension(mysize, othersize);
329 mysize = inner_face.second.n_hanging_faces_larger_subdomain;
330 MPI_Sendrecv(&mysize,
331 1,
332 MPI_UNSIGNED,
333 inner_face.first,
334 800 + my_domain,
335 &othersize,
336 1,
337 MPI_UNSIGNED,
338 inner_face.first,
339 800 + inner_face.first,
340 comm,
341 &status);
342 AssertDimension(mysize, othersize);
343# endif
344
345 // Arrange the face "ownership" such that cells that are access
346 // by more than one face (think of a cell in a corner) get
347 // ghosted. This arrangement has the advantage that we need to
348 // send less data because the same data is used twice. The
349 // strategy applied here is to ensure the same order of face
350 // pairs on both processors that share some faces, and make the
351 // same decision on both sides.
352
353 // Create a vector with cell ids sorted over the processor with
354 // the larger rank. In the code below we need to be able to
355 // identify the same cell once for the processor with higher
356 // rank and once for the processor with the lower rank. The
357 // format for the processor with the higher rank is already
358 // contained in `shared_faces`, whereas we need a copy that we
359 // sort differently for the other way around.
360 std::vector<std::tuple<CellId, CellId, unsigned int>> other_range(
361 inner_face.second.shared_faces.size());
362 for (unsigned int i = 0; i < other_range.size(); ++i)
363 other_range[i] =
364 std::make_tuple(inner_face.second.shared_faces[i].second,
365 inner_face.second.shared_faces[i].first,
366 i);
367 std::sort(other_range.begin(), other_range.end());
368
369 // the vector 'assignment' sets whether a particular cell
370 // appears more often and acts as a pre-selection of the rank. A
371 // value of 1 means that the process with the higher rank gets
372 // those faces, a value -1 means that the process with the lower
373 // rank gets it, whereas a value 0 means that the decision can
374 // be made in an arbitrary way.
375 unsigned int n_faces_lower_proc = 0, n_faces_higher_proc = 0;
376 std::vector<char> assignment(other_range.size(), 0);
377 if (inner_face.second.shared_faces.size() > 0)
378 {
379 // identify faces that go to the processor with the higher
380 // rank
381 unsigned int count = 0;
382 for (unsigned int i = 1;
383 i < inner_face.second.shared_faces.size();
384 ++i)
385 if (inner_face.second.shared_faces[i].first ==
386 inner_face.second.shared_faces[i - 1 - count].first)
387 ++count;
388 else
389 {
390 AssertThrow(count < 2 * dim, ExcInternalError());
391 if (count > 0)
392 {
393 for (unsigned int k = 0; k <= count; ++k)
394 assignment[i - 1 - k] = 1;
395 n_faces_higher_proc += count + 1;
396 }
397 count = 0;
398 }
399
400 // identify faces that definitely go to the processor with
401 // the lower rank - this must use the sorting of CellId
402 // variables from the processor with the higher rank, i.e.,
403 // other_range rather than `shared_faces`.
404 count = 0;
405 for (unsigned int i = 1; i < other_range.size(); ++i)
406 if (std::get<0>(other_range[i]) ==
407 std::get<0>(other_range[i - 1 - count]))
408 ++count;
409 else
410 {
411 AssertThrow(count < 2 * dim, ExcInternalError());
412 if (count > 0)
413 {
414 for (unsigned int k = 0; k <= count; ++k)
415 {
416 Assert(inner_face.second
417 .shared_faces[std::get<2>(
418 other_range[i - 1])]
419 .second ==
420 inner_face.second
421 .shared_faces[std::get<2>(
422 other_range[i - 1 - k])]
423 .second,
425 // only assign to -1 if higher rank was not
426 // yet set
427 if (assignment[std::get<2>(
428 other_range[i - 1 - k])] == 0)
429 {
430 assignment[std::get<2>(
431 other_range[i - 1 - k])] = -1;
432 ++n_faces_lower_proc;
433 }
434 }
435 }
436 count = 0;
437 }
438 }
439
440
441 // divide the faces evenly between the two processors. the
442 // processor with small rank takes the first half, the processor
443 // with larger rank the second half. Adjust for the hanging
444 // faces that always get assigned to one side, and the faces we
445 // have already assigned due to the criterion above
446 n_faces_lower_proc +=
447 inner_face.second.n_hanging_faces_smaller_subdomain;
448 n_faces_higher_proc +=
449 inner_face.second.n_hanging_faces_larger_subdomain;
450 const unsigned int n_total_faces_at_proc_boundary =
451 (inner_face.second.shared_faces.size() +
452 inner_face.second.n_hanging_faces_smaller_subdomain +
453 inner_face.second.n_hanging_faces_larger_subdomain);
454 unsigned int split_index = n_total_faces_at_proc_boundary / 2;
455 if (split_index < n_faces_lower_proc)
456 split_index = 0;
457 else if (split_index <
458 n_total_faces_at_proc_boundary - n_faces_higher_proc)
459 split_index -= n_faces_lower_proc;
460 else
461 split_index = n_total_faces_at_proc_boundary -
462 n_faces_higher_proc - n_faces_lower_proc;
463
464 // make sure the splitting is consistent between both sides
465# if defined(DEAL_II_WITH_MPI) && defined(DEBUG)
466 MPI_Sendrecv(&split_index,
467 1,
468 MPI_UNSIGNED,
469 inner_face.first,
470 900 + my_domain,
471 &othersize,
472 1,
473 MPI_UNSIGNED,
474 inner_face.first,
475 900 + inner_face.first,
476 comm,
477 &status);
478 AssertDimension(split_index, othersize);
479 MPI_Sendrecv(&n_faces_lower_proc,
480 1,
481 MPI_UNSIGNED,
482 inner_face.first,
483 1000 + my_domain,
484 &othersize,
485 1,
486 MPI_UNSIGNED,
487 inner_face.first,
488 1000 + inner_face.first,
489 comm,
490 &status);
491 AssertDimension(n_faces_lower_proc, othersize);
492 MPI_Sendrecv(&n_faces_higher_proc,
493 1,
494 MPI_UNSIGNED,
495 inner_face.first,
496 1100 + my_domain,
497 &othersize,
498 1,
499 MPI_UNSIGNED,
500 inner_face.first,
501 1100 + inner_face.first,
502 comm,
503 &status);
504 AssertDimension(n_faces_higher_proc, othersize);
505# endif
506
507 // collect the faces on both sides
508 std::vector<std::pair<CellId, CellId>> owned_faces_lower,
509 owned_faces_higher;
510 for (unsigned int i = 0; i < assignment.size(); ++i)
511 if (assignment[i] < 0)
512 owned_faces_lower.push_back(
513 inner_face.second.shared_faces[i]);
514 else if (assignment[i] > 0)
515 owned_faces_higher.push_back(
516 inner_face.second.shared_faces[i]);
517 AssertIndexRange(split_index,
518 inner_face.second.shared_faces.size() + 1 -
519 owned_faces_lower.size() -
520 owned_faces_higher.size());
521
522 unsigned int i = 0, c = 0;
523 for (; i < assignment.size() && c < split_index; ++i)
524 if (assignment[i] == 0)
525 {
526 owned_faces_lower.push_back(
527 inner_face.second.shared_faces[i]);
528 ++c;
529 }
530 for (; i < assignment.size(); ++i)
531 if (assignment[i] == 0)
532 {
533 owned_faces_higher.push_back(
534 inner_face.second.shared_faces[i]);
535 }
536
537# ifdef DEBUG
538 // check consistency of faces on both sides
539 std::vector<std::pair<CellId, CellId>> check_faces;
540 check_faces.insert(check_faces.end(),
541 owned_faces_lower.begin(),
542 owned_faces_lower.end());
543 check_faces.insert(check_faces.end(),
544 owned_faces_higher.begin(),
545 owned_faces_higher.end());
546 std::sort(check_faces.begin(), check_faces.end());
547 AssertDimension(check_faces.size(),
548 inner_face.second.shared_faces.size());
549 for (unsigned int i = 0; i < check_faces.size(); ++i)
550 Assert(check_faces[i] == inner_face.second.shared_faces[i],
552# endif
553
554 // now only set half of the faces as the ones to keep
555 if (my_domain < inner_face.first)
556 inner_face.second.shared_faces.swap(owned_faces_lower);
557 else
558 inner_face.second.shared_faces.swap(owned_faces_higher);
559
560 std::sort(inner_face.second.shared_faces.begin(),
561 inner_face.second.shared_faces.end());
562 }
563 }
564
565 // fill in the additional cells that we need access to via ghosting to
566 // cell_levels
567 std::set<std::pair<unsigned int, unsigned int>> ghost_cells;
568 for (unsigned int i = 0; i < cell_levels.size(); ++i)
569 {
570 typename ::Triangulation<dim>::cell_iterator dcell(
571 &triangulation, cell_levels[i].first, cell_levels[i].second);
572 if (use_active_cells)
573 Assert(dcell->is_active(), ExcNotImplemented());
574 for (const auto f : dcell->face_indices())
575 {
576 if (dcell->at_boundary(f) && !dcell->has_periodic_neighbor(f))
577 face_is_owned[dcell->face(f)->index()] =
578 FaceCategory::locally_active_at_boundary;
579
580 // treat boundaries of cells of different refinement level
581 // inside the domain in case of multigrid separately
582 else if ((dcell->at_boundary(f) == false ||
583 dcell->has_periodic_neighbor(f)) &&
584 mg_level != numbers::invalid_unsigned_int &&
585 dcell->neighbor_or_periodic_neighbor(f)->level() <
586 dcell->level())
587 {
588 face_is_owned[dcell->face(f)->index()] =
589 FaceCategory::multigrid_refinement_edge;
590 }
591 else
592 {
593 typename ::Triangulation<dim>::cell_iterator neighbor =
594 dcell->neighbor_or_periodic_neighbor(f);
595
596 // neighbor is refined -> face will be treated by neighbor
597 if (use_active_cells && neighbor->has_children() &&
598 hold_all_faces_to_owned_cells == false)
599 continue;
600
601 bool add_to_ghost = false;
603 id1 = use_active_cells ? dcell->subdomain_id() :
604 dcell->level_subdomain_id(),
605 id2 = use_active_cells ?
606 (neighbor->has_children() ?
607 dcell->neighbor_child_on_subface(f, 0)
608 ->subdomain_id() :
609 neighbor->subdomain_id()) :
610 neighbor->level_subdomain_id();
611
612 // Check whether the current face should be processed
613 // locally (instead of being processed from the other
614 // side). We process a face locally when we are more refined
615 // (in the active cell case) or when the face is listed in
616 // the `shared_faces` data structure that we built above.
617 if ((id1 == id2 &&
618 (use_active_cells == false || neighbor->is_active())) ||
619 dcell->level() > neighbor->level() ||
620 std::binary_search(
621 inner_faces_at_proc_boundary[id2].shared_faces.begin(),
622 inner_faces_at_proc_boundary[id2].shared_faces.end(),
623 std::make_pair(id1 < id2 ? dcell->id() : neighbor->id(),
624 id1 < id2 ? neighbor->id() :
625 dcell->id())))
626 {
627 face_is_owned[dcell->face(f)->index()] =
628 FaceCategory::locally_active_done_here;
629 if (dcell->level() == neighbor->level())
630 face_is_owned
631 [neighbor
632 ->face(dcell->has_periodic_neighbor(f) ?
633 dcell->periodic_neighbor_face_no(f) :
634 dcell->neighbor_face_no(f))
635 ->index()] =
636 FaceCategory::locally_active_done_here;
637
638 // If neighbor is a ghost element (i.e.
639 // dcell->subdomain_id !
640 // dcell->neighbor(f)->subdomain_id()), we need to add its
641 // index into cell level list.
642 if (use_active_cells)
643 add_to_ghost =
644 (dcell->subdomain_id() != neighbor->subdomain_id());
645 else
646 add_to_ghost = (dcell->level_subdomain_id() !=
647 neighbor->level_subdomain_id());
648 }
649 else if (hold_all_faces_to_owned_cells == true)
650 {
651 // add all cells to ghost layer...
652 face_is_owned[dcell->face(f)->index()] =
653 FaceCategory::ghosted;
654 if (use_active_cells)
655 {
656 if (neighbor->has_children())
657 for (unsigned int s = 0;
658 s < dcell->face(f)->n_children();
659 ++s)
660 if (dcell->at_boundary(f))
661 {
662 if (dcell
663 ->periodic_neighbor_child_on_subface(f,
664 s)
665 ->subdomain_id() !=
666 dcell->subdomain_id())
667 add_to_ghost = true;
668 }
669 else
670 {
671 if (dcell->neighbor_child_on_subface(f, s)
672 ->subdomain_id() !=
673 dcell->subdomain_id())
674 add_to_ghost = true;
675 }
676 else
677 add_to_ghost = (dcell->subdomain_id() !=
678 neighbor->subdomain_id());
679 }
680 else
681 add_to_ghost = (dcell->level_subdomain_id() !=
682 neighbor->level_subdomain_id());
683 }
684
685 if (add_to_ghost)
686 {
687 if (use_active_cells && neighbor->has_children())
688 for (unsigned int s = 0;
689 s < dcell->face(f)->n_children();
690 ++s)
691 {
692 typename ::Triangulation<dim>::cell_iterator
693 neighbor_child =
694 dcell->at_boundary(f) ?
695 dcell->periodic_neighbor_child_on_subface(f,
696 s) :
697 dcell->neighbor_child_on_subface(f, s);
698 if (neighbor_child->subdomain_id() !=
699 dcell->subdomain_id())
700 ghost_cells.insert(
701 std::pair<unsigned int, unsigned int>(
702 neighbor_child->level(),
703 neighbor_child->index()));
704 }
705 else
706 ghost_cells.insert(
707 std::pair<unsigned int, unsigned int>(
708 neighbor->level(), neighbor->index()));
709 at_processor_boundary[i] = true;
710 }
711 }
712 }
713 }
714
715 // step 2: append the ghost cells at the end of the locally owned
716 // cells
717 for (const auto &ghost_cell : ghost_cells)
718 cell_levels.push_back(ghost_cell);
719 }
720
721
722
723 template <int dim>
724 void
725 FaceSetup<dim>::generate_faces(
726 const ::Triangulation<dim> & triangulation,
727 const std::vector<std::pair<unsigned int, unsigned int>> &cell_levels,
728 TaskInfo & task_info)
729 {
730 // step 1: create the inverse map between cell iterators and the
731 // cell_level_index field
732 std::map<std::pair<unsigned int, unsigned int>, unsigned int>
733 map_to_vectorized;
734 for (unsigned int cell = 0; cell < cell_levels.size(); ++cell)
735 if (cell == 0 || cell_levels[cell] != cell_levels[cell - 1])
736 {
737 typename ::Triangulation<dim>::cell_iterator dcell(
739 cell_levels[cell].first,
740 cell_levels[cell].second);
741 std::pair<unsigned int, unsigned int> level_index(dcell->level(),
742 dcell->index());
743 map_to_vectorized[level_index] = cell;
744 }
745
746 // step 2: fill the information about inner faces and boundary faces
747 const unsigned int vectorization_length = task_info.vectorization_length;
748 task_info.face_partition_data.resize(
749 task_info.cell_partition_data.size() - 1, 0);
750 task_info.boundary_partition_data.resize(
751 task_info.cell_partition_data.size() - 1, 0);
752 std::vector<unsigned char> face_visited(face_is_owned.size(), 0);
753 for (unsigned int partition = 0;
754 partition < task_info.cell_partition_data.size() - 2;
755 ++partition)
756 {
757 unsigned int boundary_counter = 0;
758 unsigned int inner_counter = 0;
759 for (unsigned int cell = task_info.cell_partition_data[partition] *
760 vectorization_length;
761 cell < task_info.cell_partition_data[partition + 1] *
762 vectorization_length;
763 ++cell)
764 if (cell == 0 || cell_levels[cell] != cell_levels[cell - 1])
765 {
766 typename ::Triangulation<dim>::cell_iterator dcell(
768 cell_levels[cell].first,
769 cell_levels[cell].second);
770 for (const auto f : dcell->face_indices())
771 {
772 // boundary face
773 if (face_is_owned[dcell->face(f)->index()] ==
774 FaceCategory::locally_active_at_boundary)
775 {
776 Assert(dcell->at_boundary(f), ExcInternalError());
777 ++boundary_counter;
778 FaceToCellTopology<1> info;
779 info.cells_interior[0] = cell;
780 info.cells_exterior[0] = numbers::invalid_unsigned_int;
781 info.interior_face_no = f;
782 info.exterior_face_no = dcell->face(f)->boundary_id();
783 info.face_type =
784 dcell->face(f)->reference_cell() !=
786 info.subface_index =
788 info.face_orientation = 0;
789 boundary_faces.push_back(info);
790
791 face_visited[dcell->face(f)->index()]++;
792 }
793 // interior face, including faces over periodic boundaries
794 else
795 {
796 typename ::Triangulation<dim>::cell_iterator
797 neighbor = dcell->neighbor_or_periodic_neighbor(f);
798 if (use_active_cells && neighbor->has_children())
799 {
800 for (unsigned int c = 0;
801 c < dcell->face(f)->n_children();
802 ++c)
803 {
805 dim>::cell_iterator neighbor_c =
806 dcell->at_boundary(f) ?
807 dcell->periodic_neighbor_child_on_subface(
808 f, c) :
809 dcell->neighbor_child_on_subface(f, c);
810 const types::subdomain_id neigh_domain =
811 neighbor_c->subdomain_id();
812 const unsigned int neighbor_face_no =
813 dcell->has_periodic_neighbor(f) ?
814 dcell->periodic_neighbor_face_no(f) :
815 dcell->neighbor_face_no(f);
816 if (neigh_domain != dcell->subdomain_id() ||
817 face_visited
818 [dcell->face(f)->child(c)->index()] ==
819 1)
820 {
821 std::pair<unsigned int, unsigned int>
822 level_index(neighbor_c->level(),
823 neighbor_c->index());
824 if (face_is_owned
825 [dcell->face(f)->child(c)->index()] ==
826 FaceCategory::locally_active_done_here)
827 {
828 ++inner_counter;
829 inner_faces.push_back(create_face(
830 neighbor_face_no,
831 neighbor_c,
832 map_to_vectorized[level_index],
833 dcell,
834 cell));
835 }
836 else if (face_is_owned[dcell->face(f)
837 ->child(c)
838 ->index()] ==
839 FaceCategory::ghosted)
840 {
841 inner_ghost_faces.push_back(create_face(
842 neighbor_face_no,
843 neighbor_c,
844 map_to_vectorized[level_index],
845 dcell,
846 cell));
847 }
848 else
849 Assert(
850 face_is_owned[dcell->face(f)
851 ->index()] ==
852 FaceCategory::
853 locally_active_done_elsewhere ||
854 face_is_owned[dcell->face(f)
855 ->index()] ==
856 FaceCategory::ghosted,
858 }
859 else
860 {
861 face_visited
862 [dcell->face(f)->child(c)->index()] = 1;
863 }
864 }
865 }
866 else
867 {
868 const types::subdomain_id my_domain =
869 use_active_cells ? dcell->subdomain_id() :
870 dcell->level_subdomain_id();
871 const types::subdomain_id neigh_domain =
872 use_active_cells ? neighbor->subdomain_id() :
873 neighbor->level_subdomain_id();
874 if (neigh_domain != my_domain ||
875 face_visited[dcell->face(f)->index()] == 1)
876 {
877 std::pair<unsigned int, unsigned int>
878 level_index(neighbor->level(),
879 neighbor->index());
880 if (face_is_owned[dcell->face(f)->index()] ==
881 FaceCategory::locally_active_done_here)
882 {
883 Assert(use_active_cells ||
884 dcell->level() ==
885 neighbor->level(),
887 ++inner_counter;
888 inner_faces.push_back(create_face(
889 f,
890 dcell,
891 cell,
892 neighbor,
893 map_to_vectorized[level_index]));
894 }
895 else if (face_is_owned[dcell->face(f)
896 ->index()] ==
897 FaceCategory::ghosted)
898 {
899 inner_ghost_faces.push_back(create_face(
900 f,
901 dcell,
902 cell,
903 neighbor,
904 map_to_vectorized[level_index]));
905 }
906 }
907 else
908 {
909 face_visited[dcell->face(f)->index()] = 1;
910 if (dcell->has_periodic_neighbor(f))
911 face_visited
912 [neighbor
913 ->face(
914 dcell->periodic_neighbor_face_no(f))
915 ->index()] = 1;
916 }
917 if (face_is_owned[dcell->face(f)->index()] ==
918 FaceCategory::multigrid_refinement_edge)
919 {
920 refinement_edge_faces.push_back(
921 create_face(f,
922 dcell,
923 cell,
924 neighbor,
925 refinement_edge_faces.size()));
926 }
927 }
928 }
929 }
930 }
931 task_info.face_partition_data[partition + 1] =
932 task_info.face_partition_data[partition] + inner_counter;
933 task_info.boundary_partition_data[partition + 1] =
934 task_info.boundary_partition_data[partition] + boundary_counter;
935 }
936 task_info.ghost_face_partition_data.resize(2);
937 task_info.ghost_face_partition_data[0] = 0;
938 task_info.ghost_face_partition_data[1] = inner_ghost_faces.size();
939 task_info.refinement_edge_face_partition_data.resize(2);
940 task_info.refinement_edge_face_partition_data[0] = 0;
941 task_info.refinement_edge_face_partition_data[1] =
942 refinement_edge_faces.size();
943 }
944
945
946
947 template <int dim>
948 FaceToCellTopology<1>
949 FaceSetup<dim>::create_face(
950 const unsigned int face_no,
951 const typename ::Triangulation<dim>::cell_iterator &cell,
952 const unsigned int number_cell_interior,
953 const typename ::Triangulation<dim>::cell_iterator &neighbor,
954 const unsigned int number_cell_exterior)
955 {
956 FaceToCellTopology<1> info;
957 info.cells_interior[0] = number_cell_interior;
958 info.cells_exterior[0] = number_cell_exterior;
959 info.interior_face_no = face_no;
960 if (cell->has_periodic_neighbor(face_no))
961 info.exterior_face_no = cell->periodic_neighbor_face_no(face_no);
962 else
963 info.exterior_face_no = cell->neighbor_face_no(face_no);
964
965 info.face_type = cell->face(face_no)->reference_cell() !=
967
968 info.subface_index = GeometryInfo<dim>::max_children_per_cell;
969 Assert(neighbor->level() <= cell->level(), ExcInternalError());
970 if (cell->level() > neighbor->level())
971 {
972 if (cell->has_periodic_neighbor(face_no))
973 info.subface_index =
974 cell->periodic_neighbor_of_coarser_periodic_neighbor(face_no)
975 .second;
976 else
977 info.subface_index =
978 cell->neighbor_of_coarser_neighbor(face_no).second;
979 }
980
981 // special treatment of periodic boundaries
982 if (dim == 3 && cell->has_periodic_neighbor(face_no))
983 {
984 const unsigned int exterior_face_orientation =
985 !cell->get_triangulation()
986 .get_periodic_face_map()
987 .at({cell, face_no})
988 .second[0] +
989 2 * cell->get_triangulation()
990 .get_periodic_face_map()
991 .at({cell, face_no})
992 .second[1] +
993 4 * cell->get_triangulation()
994 .get_periodic_face_map()
995 .at({cell, face_no})
996 .second[2];
997
998 info.face_orientation = exterior_face_orientation;
999
1000 return info;
1001 }
1002
1003 info.face_orientation = 0;
1004 const unsigned int interior_face_orientation =
1005 !cell->face_orientation(face_no) + 2 * cell->face_flip(face_no) +
1006 4 * cell->face_rotation(face_no);
1007 const unsigned int exterior_face_orientation =
1008 !neighbor->face_orientation(info.exterior_face_no) +
1009 2 * neighbor->face_flip(info.exterior_face_no) +
1010 4 * neighbor->face_rotation(info.exterior_face_no);
1011 if (interior_face_orientation != 0)
1012 {
1013 info.face_orientation = 8 + interior_face_orientation;
1014 Assert(exterior_face_orientation == 0,
1015 ExcMessage(
1016 "Face seems to be wrongly oriented from both sides"));
1017 }
1018 else
1019 info.face_orientation = exterior_face_orientation;
1020 return info;
1021 }
1022
1023
1024
1031 inline bool
1032 compare_faces_for_vectorization(
1033 const FaceToCellTopology<1> & face1,
1034 const FaceToCellTopology<1> & face2,
1035 const std::vector<unsigned int> &active_fe_indices,
1036 const unsigned int length)
1037 {
1038 if (face1.interior_face_no != face2.interior_face_no)
1039 return false;
1040 if (face1.exterior_face_no != face2.exterior_face_no)
1041 return false;
1042 if (face1.subface_index != face2.subface_index)
1043 return false;
1044 if (face1.face_orientation != face2.face_orientation)
1045 return false;
1046 if (face1.face_type != face2.face_type)
1047 return false;
1048
1049 if (active_fe_indices.size() > 0)
1050 {
1051 if (active_fe_indices[face1.cells_interior[0] / length] !=
1052 active_fe_indices[face2.cells_interior[0] / length])
1053 return false;
1054
1055 if (face2.cells_exterior[0] != numbers::invalid_unsigned_int)
1056 if (active_fe_indices[face1.cells_exterior[0] / length] !=
1057 active_fe_indices[face2.cells_exterior[0] / length])
1058 return false;
1059 }
1060
1061 return true;
1062 }
1063
1064
1065
1072 template <int length>
1073 struct FaceComparator
1074 {
1075 FaceComparator(const std::vector<unsigned int> &active_fe_indices)
1076 : active_fe_indices(active_fe_indices)
1077 {}
1078
1079 bool
1080 operator()(const FaceToCellTopology<length> &face1,
1081 const FaceToCellTopology<length> &face2) const
1082 {
1083 // check if active FE indices match
1084 if (face1.face_type < face2.face_type)
1085 return true;
1086 else if (face1.face_type > face2.face_type)
1087 return false;
1088
1089 // check if active FE indices match
1090 if (active_fe_indices.size() > 0)
1091 {
1092 // ... for interior faces
1093 if (active_fe_indices[face1.cells_interior[0] / length] <
1094 active_fe_indices[face2.cells_interior[0] / length])
1095 return true;
1096 else if (active_fe_indices[face1.cells_interior[0] / length] >
1097 active_fe_indices[face2.cells_interior[0] / length])
1098 return false;
1099
1100 // ... for exterior faces
1101 if (face2.cells_exterior[0] != numbers::invalid_unsigned_int)
1102 {
1103 if (active_fe_indices[face1.cells_exterior[0] / length] <
1104 active_fe_indices[face2.cells_exterior[0] / length])
1105 return true;
1106 else if (active_fe_indices[face1.cells_exterior[0] / length] >
1107 active_fe_indices[face2.cells_exterior[0] / length])
1108 return false;
1109 }
1110 }
1111
1112 for (unsigned int i = 0; i < length; ++i)
1113 if (face1.cells_interior[i] < face2.cells_interior[i])
1114 return true;
1115 else if (face1.cells_interior[i] > face2.cells_interior[i])
1116 return false;
1117 for (unsigned int i = 0; i < length; ++i)
1118 if (face1.cells_exterior[i] < face2.cells_exterior[i])
1119 return true;
1120 else if (face1.cells_exterior[i] > face2.cells_exterior[i])
1121 return false;
1122 if (face1.interior_face_no < face2.interior_face_no)
1123 return true;
1124 else if (face1.interior_face_no > face2.interior_face_no)
1125 return false;
1126 if (face1.exterior_face_no < face2.exterior_face_no)
1127 return true;
1128 else if (face1.exterior_face_no > face2.exterior_face_no)
1129 return false;
1130
1131 // we do not need to check for subface_index and orientation because
1132 // those cannot be different if when all the other values are the
1133 // same.
1134 AssertDimension(face1.subface_index, face2.subface_index);
1135 AssertDimension(face1.face_orientation, face2.face_orientation);
1136
1137 return false;
1138 }
1139
1140 private:
1141 const std::vector<unsigned int> &active_fe_indices;
1142 };
1143
1144
1145
1146 template <int vectorization_width>
1147 void
1149 const std::vector<FaceToCellTopology<1>> &faces_in,
1150 const std::vector<bool> & hard_vectorization_boundary,
1151 std::vector<unsigned int> & face_partition_data,
1152 std::vector<FaceToCellTopology<vectorization_width>> &faces_out,
1153 const std::vector<unsigned int> & active_fe_indices)
1154 {
1155 FaceToCellTopology<vectorization_width> macro_face;
1156 std::vector<std::vector<unsigned int>> faces_type;
1157
1158 unsigned int face_start = face_partition_data[0],
1159 face_end = face_partition_data[0];
1160
1161 face_partition_data[0] = faces_out.size();
1162 for (unsigned int partition = 0;
1163 partition < face_partition_data.size() - 1;
1164 ++partition)
1165 {
1166 std::vector<std::vector<unsigned int>> new_faces_type;
1167
1168 // start with the end point for the last partition
1169 face_start = face_end;
1170 face_end = face_partition_data[partition + 1];
1171
1172 // set the partitioner to the new vectorized lengths
1173 face_partition_data[partition + 1] = face_partition_data[partition];
1174
1175 // loop over the faces in the current partition and reorder according
1176 // to the face type
1177 for (unsigned int face = face_start; face < face_end; ++face)
1178 {
1179 for (auto &face_type : faces_type)
1180 {
1181 // Compare current face with first face of type type
1182 if (compare_faces_for_vectorization(faces_in[face],
1183 faces_in[face_type[0]],
1184 active_fe_indices,
1185 vectorization_width))
1186 {
1187 face_type.push_back(face);
1188 goto face_found;
1189 }
1190 }
1191 faces_type.emplace_back(1, face);
1192 face_found:
1193 {}
1194 }
1195
1196 // insert new faces in sorted list to get good data locality
1197 FaceComparator<vectorization_width> face_comparator(
1198 active_fe_indices);
1199 std::set<FaceToCellTopology<vectorization_width>,
1200 FaceComparator<vectorization_width>>
1201 new_faces(face_comparator);
1202 for (const auto &face_type : faces_type)
1203 {
1204 macro_face.face_type = faces_in[face_type[0]].face_type;
1205 macro_face.interior_face_no =
1206 faces_in[face_type[0]].interior_face_no;
1207 macro_face.exterior_face_no =
1208 faces_in[face_type[0]].exterior_face_no;
1209 macro_face.subface_index = faces_in[face_type[0]].subface_index;
1210 macro_face.face_orientation =
1211 faces_in[face_type[0]].face_orientation;
1212 unsigned int no_faces = face_type.size();
1213 std::vector<unsigned char> touched(no_faces, 0);
1214
1215 // do two passes through the data. The first is to identify
1216 // similar faces within the same index range as the cells which
1217 // will allow for vectorized read operations, the second picks up
1218 // all the rest
1219 unsigned int n_vectorized = 0;
1220 for (unsigned int f = 0; f < no_faces; ++f)
1221 if (faces_in[face_type[f]].cells_interior[0] %
1222 vectorization_width ==
1223 0)
1224 {
1225 bool is_contiguous = true;
1226 if (f + vectorization_width > no_faces)
1227 is_contiguous = false;
1228 else
1229 for (unsigned int v = 1; v < vectorization_width; ++v)
1230 if (faces_in[face_type[f + v]].cells_interior[0] !=
1231 faces_in[face_type[f]].cells_interior[0] + v)
1232 is_contiguous = false;
1233 if (is_contiguous)
1234 {
1236 face_type.size() -
1237 vectorization_width + 1);
1238 for (unsigned int v = 0; v < vectorization_width; ++v)
1239 {
1240 macro_face.cells_interior[v] =
1241 faces_in[face_type[f + v]].cells_interior[0];
1242 macro_face.cells_exterior[v] =
1243 faces_in[face_type[f + v]].cells_exterior[0];
1244 touched[f + v] = 1;
1245 }
1246 new_faces.insert(macro_face);
1247 f += vectorization_width - 1;
1248 n_vectorized += vectorization_width;
1249 }
1250 }
1251
1252 std::vector<unsigned int> untouched;
1253 untouched.reserve(no_faces - n_vectorized);
1254 for (unsigned int f = 0; f < no_faces; ++f)
1255 if (touched[f] == 0)
1256 untouched.push_back(f);
1257 unsigned int v = 0;
1258 for (const auto f : untouched)
1259 {
1260 macro_face.cells_interior[v] =
1261 faces_in[face_type[f]].cells_interior[0];
1262 macro_face.cells_exterior[v] =
1263 faces_in[face_type[f]].cells_exterior[0];
1264 ++v;
1265 if (v == vectorization_width)
1266 {
1267 new_faces.insert(macro_face);
1268 v = 0;
1269 }
1270 }
1271 if (v > 0 && v < vectorization_width)
1272 {
1273 // must add non-filled face
1274 if (hard_vectorization_boundary[partition + 1] ||
1275 partition == face_partition_data.size() - 2)
1276 {
1277 for (; v < vectorization_width; ++v)
1278 {
1279 // Dummy cell, not used
1280 macro_face.cells_interior[v] =
1282 macro_face.cells_exterior[v] =
1284 }
1285 new_faces.insert(macro_face);
1286 }
1287 else
1288 {
1289 // postpone to the next partition
1290 std::vector<unsigned int> untreated(v);
1291 for (unsigned int f = 0; f < v; ++f)
1292 untreated[f] = face_type[*(untouched.end() - 1 - f)];
1293 new_faces_type.push_back(untreated);
1294 }
1295 }
1296 }
1297
1298 // insert sorted list to vector of faces
1299 for (auto it = new_faces.begin(); it != new_faces.end(); ++it)
1300 faces_out.push_back(*it);
1301 face_partition_data[partition + 1] += new_faces.size();
1302
1303 // set the faces that were left over to faces_type for the next round
1304 faces_type = std::move(new_faces_type);
1305 }
1306
1307# ifdef DEBUG
1308 // final safety checks
1309 for (const auto &face_type : faces_type)
1310 AssertDimension(face_type.size(), 0U);
1311
1312 AssertDimension(faces_out.size(), face_partition_data.back());
1313 unsigned int nfaces = 0;
1314 for (unsigned int i = face_partition_data[0];
1315 i < face_partition_data.back();
1316 ++i)
1317 for (unsigned int v = 0; v < vectorization_width; ++v)
1318 nfaces +=
1319 (faces_out[i].cells_interior[v] != numbers::invalid_unsigned_int);
1320 AssertDimension(nfaces, faces_in.size());
1321
1322 std::vector<std::pair<unsigned int, unsigned int>> in_faces, out_faces;
1323 for (const auto &face_in : faces_in)
1324 in_faces.emplace_back(face_in.cells_interior[0],
1325 face_in.cells_exterior[0]);
1326 for (unsigned int i = face_partition_data[0];
1327 i < face_partition_data.back();
1328 ++i)
1329 for (unsigned int v = 0;
1330 v < vectorization_width &&
1331 faces_out[i].cells_interior[v] != numbers::invalid_unsigned_int;
1332 ++v)
1333 out_faces.emplace_back(faces_out[i].cells_interior[v],
1334 faces_out[i].cells_exterior[v]);
1335 std::sort(in_faces.begin(), in_faces.end());
1336 std::sort(out_faces.begin(), out_faces.end());
1337 AssertDimension(in_faces.size(), out_faces.size());
1338 for (unsigned int i = 0; i < in_faces.size(); ++i)
1339 {
1340 AssertDimension(in_faces[i].first, out_faces[i].first);
1341 AssertDimension(in_faces[i].second, out_faces[i].second);
1342 }
1343# endif
1344 }
1345
1346#endif // ifndef DOXYGEN
1347
1348 } // namespace MatrixFreeFunctions
1349} // namespace internal
1350
1351
1353
1354#endif
Definition: cell_id.h:71
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
Point< 2 > second
Definition: grid_out.cc:4588
Point< 2 > first
Definition: grid_out.cc:4587
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
static const char U
constexpr const ReferenceCell & get_hypercube()
void partition(const SparsityPattern &sparsity_pattern, const unsigned int n_partitions, std::vector< unsigned int > &partition_indices, const Partitioner partitioner=Partitioner::metis)
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
Definition: types.h:276
static const unsigned int invalid_unsigned_int
Definition: types.h:196
TriangulationBase< dim, spacedim > Triangulation
Definition: tria_base.h:396
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
*braid_SplitCommworld & comm
std::vector< std::pair< CellId, CellId > > shared_faces
std::vector< FaceToCellTopology< 1 > > inner_faces
std::vector< FaceToCellTopology< 1 > > boundary_faces
std::vector< FaceToCellTopology< 1 > > refinement_edge_faces
void initialize(const ::Triangulation< dim > &triangulation, const unsigned int mg_level, const bool hold_all_faces_to_owned_cells, std::vector< std::pair< unsigned int, unsigned int > > &cell_levels)
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)
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)