Reference documentation for deal.II version 9.3.2
\(\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 
23 #include <deal.II/base/utilities.h>
24 
26 
27 #include <deal.II/grid/tria.h>
29 
32 
33 #include <fstream>
34 
35 
37 
38 
39 namespace 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  {
124  ghosted,
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,
424  ExcInternalError());
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],
551  ExcInternalError());
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;
602  const types::subdomain_id
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(
738  &triangulation,
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(
767  &triangulation,
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,
857  ExcInternalError());
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(),
886  ExcInternalError());
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  {
1235  AssertIndexRange(f,
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:396
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:397
Point< 2 > second
Definition: grid_out.cc:4588
Point< 2 > first
Definition: grid_out.cc:4587
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
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
unsigned int subdomain_id
Definition: types.h:43
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
std::vector< std::pair< CellId, CellId > > shared_faces
std::vector< FaceToCellTopology< 1 > > inner_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)
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)
void generate_faces(const ::Triangulation< dim > &triangulation, const std::vector< std::pair< unsigned int, unsigned int >> &cell_levels, TaskInfo &task_info)
std::vector< FaceToCellTopology< 1 > > inner_ghost_faces