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