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