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