Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
face_setup_internal.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2018 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii_face_setup_internal_h
18 #define dealii_face_setup_internal_h
19 
20 #include <deal.II/base/config.h>
21 
23 #include <deal.II/base/utilities.h>
24 
26 
27 #include <deal.II/grid/tria.h>
29 
32 
33 #include <fstream>
34 
35 
37 
38 
39 namespace internal
40 {
41  namespace MatrixFreeFunctions
42  {
50  {
54  {}
55 
56  std::vector<std::pair<CellId, CellId>> shared_faces;
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
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 
118 
123  enum class FaceCategory : char
124  {
128  ghosted,
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
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;
785  FaceToCellTopology<1> info;
786  info.cells_interior[0] = cell;
787  info.cells_exterior[0] = numbers::invalid_unsigned_int;
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  {
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  task_info.ghost_face_partition_data.resize(2);
941  task_info.ghost_face_partition_data[0] = 0;
942  task_info.ghost_face_partition_data[1] = inner_ghost_faces.size();
943  task_info.refinement_edge_face_partition_data.resize(2);
944  task_info.refinement_edge_face_partition_data[0] = 0;
945  task_info.refinement_edge_face_partition_data[1] =
946  refinement_edge_faces.size();
947  }
948 
949 
950 
951  template <int dim>
952  FaceToCellTopology<1>
954  const unsigned int face_no,
955  const typename ::Triangulation<dim>::cell_iterator &cell,
956  const unsigned int number_cell_interior,
957  const typename ::Triangulation<dim>::cell_iterator &neighbor,
958  const unsigned int number_cell_exterior)
959  {
960  FaceToCellTopology<1> info;
961  info.cells_interior[0] = number_cell_interior;
962  info.cells_exterior[0] = number_cell_exterior;
963  info.interior_face_no = face_no;
964  if (cell->has_periodic_neighbor(face_no))
965  info.exterior_face_no = cell->periodic_neighbor_face_no(face_no);
966  else
967  info.exterior_face_no = cell->neighbor_face_no(face_no);
968 
969  info.subface_index = GeometryInfo<dim>::max_children_per_cell;
970  Assert(neighbor->level() <= cell->level(), ExcInternalError());
971  if (cell->level() > neighbor->level())
972  {
973  if (cell->has_periodic_neighbor(face_no))
974  info.subface_index =
975  cell->periodic_neighbor_of_coarser_periodic_neighbor(face_no)
976  .second;
977  else
978  info.subface_index =
979  cell->neighbor_of_coarser_neighbor(face_no).second;
980  }
981 
982  info.face_orientation = 0;
983  const unsigned int left_face_orientation =
984  !cell->face_orientation(face_no) + 2 * cell->face_flip(face_no) +
985  4 * cell->face_rotation(face_no);
986  const unsigned int right_face_orientation =
987  !neighbor->face_orientation(info.exterior_face_no) +
988  2 * neighbor->face_flip(info.exterior_face_no) +
989  4 * neighbor->face_rotation(info.exterior_face_no);
990  if (left_face_orientation != 0)
991  {
992  info.face_orientation = 8 + left_face_orientation;
993  Assert(right_face_orientation == 0,
994  ExcMessage(
995  "Face seems to be wrongly oriented from both sides"));
996  }
997  else
998  info.face_orientation = right_face_orientation;
999  return info;
1000  }
1001 
1002 
1003 
1010  inline bool
1011  compare_faces_for_vectorization(const FaceToCellTopology<1> &face1,
1012  const FaceToCellTopology<1> &face2)
1013  {
1014  if (face1.interior_face_no != face2.interior_face_no)
1015  return false;
1016  if (face1.exterior_face_no != face2.exterior_face_no)
1017  return false;
1018  if (face1.subface_index != face2.subface_index)
1019  return false;
1020  if (face1.face_orientation != face2.face_orientation)
1021  return false;
1022  return true;
1023  }
1024 
1025 
1026 
1033  template <int length>
1034  struct FaceComparator
1035  {
1036  bool
1037  operator()(const FaceToCellTopology<length> &face1,
1038  const FaceToCellTopology<length> &face2) const
1039  {
1040  for (unsigned int i = 0; i < length; ++i)
1041  if (face1.cells_interior[i] < face2.cells_interior[i])
1042  return true;
1043  else if (face1.cells_interior[i] > face2.cells_interior[i])
1044  return false;
1045  for (unsigned int i = 0; i < length; ++i)
1046  if (face1.cells_exterior[i] < face2.cells_exterior[i])
1047  return true;
1048  else if (face1.cells_exterior[i] > face2.cells_exterior[i])
1049  return false;
1050  if (face1.interior_face_no < face2.interior_face_no)
1051  return true;
1052  else if (face1.interior_face_no > face2.interior_face_no)
1053  return false;
1054  if (face1.exterior_face_no < face2.exterior_face_no)
1055  return true;
1056  else if (face1.exterior_face_no > face2.exterior_face_no)
1057  return false;
1058 
1059  // we do not need to check for subface_index and orientation because
1060  // those cannot be different if when all the other values are the
1061  // same.
1062  AssertDimension(face1.subface_index, face2.subface_index);
1063  AssertDimension(face1.face_orientation, face2.face_orientation);
1064 
1065  return false;
1066  }
1067  };
1068 
1069 
1070 
1071  template <int vectorization_width>
1072  void
1074  const std::vector<FaceToCellTopology<1>> &faces_in,
1075  const std::vector<bool> & hard_vectorization_boundary,
1076  std::vector<unsigned int> & face_partition_data,
1077  std::vector<FaceToCellTopology<vectorization_width>> &faces_out)
1078  {
1080  std::vector<std::vector<unsigned int>> faces_type;
1081 
1082  unsigned int face_start = face_partition_data[0],
1083  face_end = face_partition_data[0];
1084 
1085  face_partition_data[0] = faces_out.size();
1086  for (unsigned int partition = 0;
1087  partition < face_partition_data.size() - 1;
1088  ++partition)
1089  {
1090  std::vector<std::vector<unsigned int>> new_faces_type;
1091 
1092  // start with the end point for the last partition
1093  face_start = face_end;
1094  face_end = face_partition_data[partition + 1];
1095 
1096  // set the partitioner to the new vectorized lengths
1097  face_partition_data[partition + 1] = face_partition_data[partition];
1098 
1099  // loop over the faces in the current partition and reorder according
1100  // to the face type
1101  for (unsigned int face = face_start; face < face_end; ++face)
1102  {
1103  for (auto &face_type : faces_type)
1104  {
1105  // Compare current face with first face of type type
1106  if (compare_faces_for_vectorization(faces_in[face],
1107  faces_in[face_type[0]]))
1108  {
1109  face_type.push_back(face);
1110  goto face_found;
1111  }
1112  }
1113  faces_type.emplace_back(1, face);
1114  face_found:
1115  {}
1116  }
1117 
1118  // insert new faces in sorted list to get good data locality
1119  std::set<FaceToCellTopology<vectorization_width>,
1120  FaceComparator<vectorization_width>>
1121  new_faces;
1122  for (const auto &face_type : faces_type)
1123  {
1124  macro_face.interior_face_no =
1125  faces_in[face_type[0]].interior_face_no;
1126  macro_face.exterior_face_no =
1127  faces_in[face_type[0]].exterior_face_no;
1128  macro_face.subface_index = faces_in[face_type[0]].subface_index;
1129  macro_face.face_orientation =
1130  faces_in[face_type[0]].face_orientation;
1131  unsigned int no_faces = face_type.size();
1132  std::vector<unsigned char> touched(no_faces, 0);
1133 
1134  // do two passes through the data. The first is to identify
1135  // similar faces within the same index range as the cells which
1136  // will allow for vectorized read operations, the second picks up
1137  // all the rest
1138  unsigned int n_vectorized = 0;
1139  for (unsigned int f = 0; f < no_faces; ++f)
1140  if (faces_in[face_type[f]].cells_interior[0] %
1141  vectorization_width ==
1142  0)
1143  {
1144  bool is_contiguous = true;
1145  if (f + vectorization_width > no_faces)
1146  is_contiguous = false;
1147  else
1148  for (unsigned int v = 1; v < vectorization_width; ++v)
1149  if (faces_in[face_type[f + v]].cells_interior[0] !=
1150  faces_in[face_type[f]].cells_interior[0] + v)
1151  is_contiguous = false;
1152  if (is_contiguous)
1153  {
1154  AssertIndexRange(f,
1155  face_type.size() -
1156  vectorization_width + 1);
1157  for (unsigned int v = 0; v < vectorization_width; ++v)
1158  {
1159  macro_face.cells_interior[v] =
1160  faces_in[face_type[f + v]].cells_interior[0];
1161  macro_face.cells_exterior[v] =
1162  faces_in[face_type[f + v]].cells_exterior[0];
1163  touched[f + v] = 1;
1164  }
1165  new_faces.insert(macro_face);
1166  f += vectorization_width - 1;
1167  n_vectorized += vectorization_width;
1168  }
1169  }
1170 
1171  std::vector<unsigned int> untouched;
1172  untouched.reserve(no_faces - n_vectorized);
1173  for (unsigned int f = 0; f < no_faces; ++f)
1174  if (touched[f] == 0)
1175  untouched.push_back(f);
1176  unsigned int v = 0;
1177  for (const auto f : untouched)
1178  {
1179  macro_face.cells_interior[v] =
1180  faces_in[face_type[f]].cells_interior[0];
1181  macro_face.cells_exterior[v] =
1182  faces_in[face_type[f]].cells_exterior[0];
1183  ++v;
1184  if (v == vectorization_width)
1185  {
1186  new_faces.insert(macro_face);
1187  v = 0;
1188  }
1189  }
1190  if (v > 0 && v < vectorization_width)
1191  {
1192  // must add non-filled face
1193  if (hard_vectorization_boundary[partition + 1] ||
1194  partition == face_partition_data.size() - 2)
1195  {
1196  for (; v < vectorization_width; ++v)
1197  {
1198  // Dummy cell, not used
1199  macro_face.cells_interior[v] =
1201  macro_face.cells_exterior[v] =
1203  }
1204  new_faces.insert(macro_face);
1205  }
1206  else
1207  {
1208  // postpone to the next partition
1209  std::vector<unsigned int> untreated(v);
1210  for (unsigned int f = 0; f < v; ++f)
1211  untreated[f] = face_type[*(untouched.end() - 1 - f)];
1212  new_faces_type.push_back(untreated);
1213  }
1214  }
1215  }
1216 
1217  // insert sorted list to vector of faces
1218  for (auto it = new_faces.begin(); it != new_faces.end(); ++it)
1219  faces_out.push_back(*it);
1220  face_partition_data[partition + 1] += new_faces.size();
1221 
1222  // set the faces that were left over to faces_type for the next round
1223  faces_type = std::move(new_faces_type);
1224  }
1225 
1226 # ifdef DEBUG
1227  // final safety checks
1228  for (const auto &face_type : faces_type)
1229  AssertDimension(face_type.size(), 0U);
1230 
1231  AssertDimension(faces_out.size(), face_partition_data.back());
1232  unsigned int nfaces = 0;
1233  for (unsigned int i = face_partition_data[0];
1234  i < face_partition_data.back();
1235  ++i)
1236  for (unsigned int v = 0; v < vectorization_width; ++v)
1237  nfaces +=
1238  (faces_out[i].cells_interior[v] != numbers::invalid_unsigned_int);
1239  AssertDimension(nfaces, faces_in.size());
1240 
1241  std::vector<std::pair<unsigned int, unsigned int>> in_faces, out_faces;
1242  for (const auto &face_in : faces_in)
1243  in_faces.emplace_back(face_in.cells_interior[0],
1244  face_in.cells_exterior[0]);
1245  for (unsigned int i = face_partition_data[0];
1246  i < face_partition_data.back();
1247  ++i)
1248  for (unsigned int v = 0;
1249  v < vectorization_width &&
1250  faces_out[i].cells_interior[v] != numbers::invalid_unsigned_int;
1251  ++v)
1252  out_faces.emplace_back(faces_out[i].cells_interior[v],
1253  faces_out[i].cells_exterior[v]);
1254  std::sort(in_faces.begin(), in_faces.end());
1255  std::sort(out_faces.begin(), out_faces.end());
1256  AssertDimension(in_faces.size(), out_faces.size());
1257  for (unsigned int i = 0; i < in_faces.size(); ++i)
1258  {
1259  AssertDimension(in_faces[i].first, out_faces[i].first);
1260  AssertDimension(in_faces[i].second, out_faces[i].second);
1261  }
1262 # endif
1263  }
1264 
1265 #endif // ifndef DOXYGEN
1266 
1267  } // namespace MatrixFreeFunctions
1268 } // namespace internal
1269 
1270 
1272 
1273 #endif
internal::MatrixFreeFunctions::FaceSetup::FaceCategory
FaceCategory
Definition: face_setup_internal.h:123
internal::MatrixFreeFunctions::FaceSetup
Definition: face_setup_internal.h:74
internal::MatrixFreeFunctions::FaceToCellTopology::cells_interior
unsigned int cells_interior[vectorization_width]
Definition: face_info.h:63
tria_accessor.h
internal::MatrixFreeFunctions::FaceSetup::FaceCategory::ghosted
@ ghosted
internal::MatrixFreeFunctions::collect_faces_vectorization
void collect_faces_vectorization(const std::vector< FaceToCellTopology< 1 >> &faces_in, const std::vector< bool > &hard_vectorization_boundary, std::vector< unsigned int > &face_partition_data, std::vector< FaceToCellTopology< vectorization_width >> &faces_out)
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
tria.h
memory_consumption.h
utilities.h
internal::MatrixFreeFunctions::FaceToCellTopology< 1 >
LAPACKSupport::U
static const char U
Definition: lapack_support.h:167
GeometryInfo
Definition: geometry_info.h:1224
FaceToCellTopology< vectorization_width >
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
MPI_Comm
face_info.h
second
Point< 2 > second
Definition: grid_out.cc:4353
SparsityTools::partition
void partition(const SparsityPattern &sparsity_pattern, const unsigned int n_partitions, std::vector< unsigned int > &partition_indices, const Partitioner partitioner=Partitioner::metis)
Definition: sparsity_tools.cc:400
internal::MatrixFreeFunctions::TaskInfo
Definition: task_info.h:105
internal::MatrixFreeFunctions::FaceIdentifier::FaceIdentifier
FaceIdentifier()
Definition: face_setup_internal.h:51
internal::MatrixFreeFunctions::FaceSetup::inner_ghost_faces
std::vector< FaceToCellTopology< 1 > > inner_ghost_faces
Definition: face_setup_internal.h:136
internal::MatrixFreeFunctions::FaceSetup::FaceSetup
FaceSetup()
internal::MatrixFreeFunctions::FaceSetup::create_face
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)
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
CellId
Definition: cell_id.h:69
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
internal::MatrixFreeFunctions::FaceSetup::generate_faces
void generate_faces(const ::Triangulation< dim > &triangulation, const std::vector< std::pair< unsigned int, unsigned int >> &cell_levels, TaskInfo &task_info)
internal::MatrixFreeFunctions::FaceSetup::use_active_cells
bool use_active_cells
Definition: face_setup_internal.h:117
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
internal::MatrixFreeFunctions::FaceSetup::refinement_edge_faces
std::vector< FaceToCellTopology< 1 > > refinement_edge_faces
Definition: face_setup_internal.h:137
task_info.h
internal::MatrixFreeFunctions::FaceIdentifier::shared_faces
std::vector< std::pair< CellId, CellId > > shared_faces
Definition: face_setup_internal.h:56
internal::MatrixFreeFunctions::FaceSetup::face_is_owned
std::vector< FaceCategory > face_is_owned
Definition: face_setup_internal.h:132
internal::MatrixFreeFunctions::FaceSetup::initialize
void initialize(const ::Triangulation< dim > &triangulation, const MFAddData &additional_data, std::vector< std::pair< unsigned int, unsigned int >> &cell_levels)
parallel::Triangulation
TriangulationBase< dim, spacedim > Triangulation
Definition: tria_base.h:302
unsigned int
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
internal::MatrixFreeFunctions::FaceSetup::FaceCategory::locally_active_done_here
@ locally_active_done_here
internal::MatrixFreeFunctions::FaceSetup::FaceCategory::locally_active_done_elsewhere
@ locally_active_done_elsewhere
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
internal::MatrixFreeFunctions::FaceIdentifier::n_hanging_faces_larger_subdomain
unsigned int n_hanging_faces_larger_subdomain
Definition: face_setup_internal.h:58
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:191
internal::MatrixFreeFunctions::FaceSetup::boundary_faces
std::vector< FaceToCellTopology< 1 > > boundary_faces
Definition: face_setup_internal.h:135
config.h
numbers::invalid_subdomain_id
const types::subdomain_id invalid_subdomain_id
Definition: types.h:285
triangulation
const typename ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69
internal
Definition: aligned_vector.h:369
internal::MatrixFreeFunctions::FaceSetup::FaceCategory::locally_active_at_boundary
@ locally_active_at_boundary
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
internal::MatrixFreeFunctions::FaceIdentifier::n_hanging_faces_smaller_subdomain
unsigned int n_hanging_faces_smaller_subdomain
Definition: face_setup_internal.h:57
first
Point< 2 > first
Definition: grid_out.cc:4352
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
internal::MatrixFreeFunctions::FaceSetup::inner_faces
std::vector< FaceToCellTopology< 1 > > inner_faces
Definition: face_setup_internal.h:134
internal::MatrixFreeFunctions::FaceSetup::at_processor_boundary
std::vector< bool > at_processor_boundary
Definition: face_setup_internal.h:133
tria_base.h
internal::MatrixFreeFunctions::FaceIdentifier
Definition: face_setup_internal.h:49
internal::MatrixFreeFunctions::FaceSetup::FaceCategory::multigrid_refinement_edge
@ multigrid_refinement_edge