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