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