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