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