Reference documentation for deal.II version 9.3.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\}}\)
mesh_loop.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2021 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_mesh_worker_mesh_loop_h
18 #define dealii_mesh_worker_mesh_loop_h
19 
20 #include <deal.II/base/config.h>
21 
23 #include <deal.II/base/types.h>
25 
27 #include <deal.II/grid/tria.h>
28 
34 
35 #include <functional>
36 #include <type_traits>
37 
39 
40 // Forward declaration
41 #ifndef DOXYGEN
42 template <typename>
43 class TriaActiveIterator;
44 #endif
45 
46 namespace MeshWorker
47 {
48  namespace internal
49  {
54  template <class CellIteratorType>
56  {
60  using type = CellIteratorType;
61  };
62 
70  template <class CellIteratorType>
71  struct CellIteratorBaseType<IteratorOverIterators<CellIteratorType>>
72  {
76  // Since we can have filtered iterators and the like as template
77  // arguments, we recursivelyremove the template layers to retrieve the
78  // underlying iterator type.
80  };
81 
90  template <class CellIteratorType>
91  struct CellIteratorBaseType<FilteredIterator<CellIteratorType>>
92  {
96  // Since we can have nested filtered iterators, we recursively
97  // remove the template layers to retrieve the underlying iterator type.
99  };
100  } // namespace internal
101 
102 #ifdef DOXYGEN
103 
107  using CellWorkerFunctionType = std::function<
108  void(const CellIteratorBaseType &, ScratchData &, CopyData &)>;
109 
114  using CopierFunctionType = std::function<void(const CopyData &)>;
115 
121  std::function<void(const CellIteratorBaseType &,
122  const unsigned int,
123  ScratchData &,
124  CopyData &)>;
125 
130  using FaceWorkerFunctionType =
131  std::function<void(const CellIteratorBaseType &,
132  const unsigned int,
133  const unsigned int,
134  const CellIteratorBaseType &,
135  const unsigned int,
136  const unsigned int,
137  ScratchData &,
138  CopyData &)>;
139 #endif
140 
276  template <class CellIteratorType,
277  class ScratchData,
278  class CopyData,
279  class CellIteratorBaseType =
281  void
283 #ifdef DOXYGEN
284  const CellIteratorType &begin,
285  const CellIteratorType &end,
286 
287  const CellWorkerFunctionType &cell_worker,
288  const CopierType & copier,
289 
290  const ScratchData &sample_scratch_data,
291  const CopyData & sample_copy_data,
292 
293  const AssembleFlags flags = assemble_own_cells,
294 
295  const BoundaryWorkerFunctionType &boundary_worker =
297 
298  const FaceWorkerFunctionType &face_worker = FaceWorkerFunctionType(),
299  const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
300  const unsigned int chunk_size = 8
301 #else
302  const CellIteratorType & begin,
303  const typename identity<CellIteratorType>::type &end,
304 
305  const typename identity<std::function<
306  void(const CellIteratorBaseType &, ScratchData &, CopyData &)>>::type
307  &cell_worker,
308  const typename identity<std::function<void(const CopyData &)>>::type
309  &copier,
310 
311  const ScratchData &sample_scratch_data,
312  const CopyData & sample_copy_data,
313 
314  const AssembleFlags flags = assemble_own_cells,
315 
316  const typename identity<std::function<void(const CellIteratorBaseType &,
317  const unsigned int,
318  ScratchData &,
319  CopyData &)>>::type
320  &boundary_worker = std::function<void(const CellIteratorBaseType &,
321  const unsigned int,
322  ScratchData &,
323  CopyData &)>(),
324 
325  const typename identity<std::function<void(const CellIteratorBaseType &,
326  const unsigned int,
327  const unsigned int,
328  const CellIteratorBaseType &,
329  const unsigned int,
330  const unsigned int,
331  ScratchData &,
332  CopyData &)>>::type
333  &face_worker = std::function<void(const CellIteratorBaseType &,
334  const unsigned int,
335  const unsigned int,
336  const CellIteratorBaseType &,
337  const unsigned int,
338  const unsigned int,
339  ScratchData &,
340  CopyData &)>(),
341 
342  const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
343  const unsigned int chunk_size = 8
344 #endif
345  )
346  {
347  Assert(
348  (!cell_worker) == !(flags & work_on_cells),
349  ExcMessage(
350  "If you provide a cell worker function, you also need to request "
351  "that work should be done on cells by setting the 'work_on_cells' flag. "
352  "Conversely, if you don't provide a cell worker function, you "
353  "cannot set the 'work_on_cells' flag. One of these two "
354  "conditions is not satisfied."));
355 
360  ExcMessage(
361  "If you provide a face worker function, you also need to request "
362  "that work should be done on interior faces by setting either the "
363  "'assemble_own_interior_faces_once' flag or the "
364  "'assemble_own_interior_faces_both' flag. "
365  "Conversely, if you don't provide a face worker function, you "
366  "cannot set either of these two flags. One of these two "
367  "conditions is not satisfied."));
368 
371  ExcMessage(
372  "You can only 'specify assemble_ghost_faces_once' "
373  "OR 'assemble_ghost_faces_both', but not both of these flags."));
374 
375  Assert(
376  !(flags & cells_after_faces) ||
378  ExcMessage(
379  "The option 'cells_after_faces' only makes sense if you assemble on cells."));
380 
381  Assert(
382  (!face_worker) == !(flags & work_on_faces),
383  ExcMessage(
384  "If you provide a face worker function, you also need to request "
385  "that work should be done on faces by setting the 'work_on_faces' flag. "
386  "Conversely, if you don't provide a face worker function, you "
387  "cannot set the 'work_on_faces' flag. One of these two "
388  "conditions is not satisfied."));
389 
390  Assert(
391  (!boundary_worker) == !(flags & assemble_boundary_faces),
392  ExcMessage(
393  "If you provide a boundary face worker function, you also need to request "
394  "that work should be done on boundary faces by setting the 'assemble_boundary_faces' flag. "
395  "Conversely, if you don't provide a boundary face worker function, you "
396  "cannot set the 'assemble_boundary_faces' flag. One of these two "
397  "conditions is not satisfied."));
398 
399  auto cell_action = [&](const CellIteratorBaseType &cell,
400  ScratchData & scratch,
401  CopyData & copy) {
402  // First reset the CopyData class to the empty copy_data given by the
403  // user.
404  copy = sample_copy_data;
405 
406  // Store the dimension in which we are working for later use
407  const auto dim = cell->get_triangulation().dimension;
408 
409  const bool ignore_subdomain =
410  (cell->get_triangulation().locally_owned_subdomain() ==
412 
413  types::subdomain_id current_subdomain_id =
414  (cell->is_level_cell() ? cell->level_subdomain_id() :
415  cell->subdomain_id());
416 
417  const bool own_cell =
418  ignore_subdomain ||
419  (current_subdomain_id ==
420  cell->get_triangulation().locally_owned_subdomain());
421 
422  if ((!ignore_subdomain) &&
423  (current_subdomain_id == numbers::artificial_subdomain_id))
424  return;
425 
426  if (!(flags & (cells_after_faces)) &&
427  (((flags & (assemble_own_cells)) && own_cell) ||
428  ((flags & assemble_ghost_cells) && !own_cell)))
429  cell_worker(cell, scratch, copy);
430 
431  if (flags & (work_on_faces | work_on_boundary))
432  for (const unsigned int face_no : cell->face_indices())
433  {
434  if (cell->at_boundary(face_no) &&
435  !cell->has_periodic_neighbor(face_no))
436  {
437  // only integrate boundary faces of own cells
438  if ((flags & assemble_boundary_faces) && own_cell)
439  boundary_worker(cell, face_no, scratch, copy);
440  }
441  else
442  {
443  // interior face, potentially assemble
445  neighbor = cell->neighbor_or_periodic_neighbor(face_no);
446 
447  types::subdomain_id neighbor_subdomain_id =
449  if (neighbor->is_level_cell())
450  neighbor_subdomain_id = neighbor->level_subdomain_id();
451  // subdomain id is only valid for active cells
452  else if (neighbor->is_active())
453  neighbor_subdomain_id = neighbor->subdomain_id();
454 
455  const bool own_neighbor =
456  ignore_subdomain ||
457  (neighbor_subdomain_id ==
458  cell->get_triangulation().locally_owned_subdomain());
459 
460  // skip all faces between two ghost cells
461  if (!own_cell && !own_neighbor)
462  continue;
463 
464  // skip if the user doesn't want faces between own cells
465  if (own_cell && own_neighbor &&
468  continue;
469 
470  // skip face to ghost
471  if (own_cell != own_neighbor &&
472  !(flags &
474  continue;
475 
476  // Deal with refinement edges from the refined side. Assuming
477  // one-irregular meshes, this situation should only occur if
478  // both cells are active.
479  const bool periodic_neighbor =
480  cell->has_periodic_neighbor(face_no);
481 
482  if (dim > 1 && ((!periodic_neighbor &&
483  cell->neighbor_is_coarser(face_no) &&
484  neighbor->is_active()) ||
485  (periodic_neighbor &&
486  cell->periodic_neighbor_is_coarser(face_no) &&
487  neighbor->is_active())))
488  {
489  Assert(cell->is_active(), ExcInternalError());
490 
491  // skip if only one processor needs to assemble the face
492  // to a ghost cell and the fine cell is not ours.
493  if (!own_cell && (flags & assemble_ghost_faces_once))
494  continue;
495 
496  const std::pair<unsigned int, unsigned int>
497  neighbor_face_no =
498  periodic_neighbor ?
499  cell->periodic_neighbor_of_coarser_periodic_neighbor(
500  face_no) :
501  cell->neighbor_of_coarser_neighbor(face_no);
502 
503  face_worker(cell,
504  face_no,
506  neighbor,
507  neighbor_face_no.first,
508  neighbor_face_no.second,
509  scratch,
510  copy);
511 
513  {
514  // If own faces are to be assembled from both sides,
515  // call the faceworker again with swapped arguments.
516  // This is because we won't be looking at an adaptively
517  // refined edge coming from the other side.
518  face_worker(neighbor,
519  neighbor_face_no.first,
520  neighbor_face_no.second,
521  cell,
522  face_no,
524  scratch,
525  copy);
526  }
527  }
528  else if (dim == 1 && cell->level() > neighbor->level())
529  {
530  // In one dimension, there is no other check to do
531  const unsigned int neighbor_face_no =
532  periodic_neighbor ?
533  cell->periodic_neighbor_face_no(face_no) :
534  cell->neighbor_face_no(face_no);
535  Assert(periodic_neighbor ||
536  neighbor->face(neighbor_face_no) ==
537  cell->face(face_no),
538  ExcInternalError());
539 
540  face_worker(cell,
541  face_no,
543  neighbor,
544  neighbor_face_no,
546  scratch,
547  copy);
548 
550  {
551  // If own faces are to be assembled from both sides,
552  // call the faceworker again with swapped arguments.
553  face_worker(neighbor,
554  neighbor_face_no,
556  cell,
557  face_no,
559  scratch,
560  copy);
561  }
562  }
563  else
564  {
565  // If iterator is active and neighbor is refined, skip
566  // internal face.
567  if (::internal::is_active_iterator(cell) &&
568  neighbor->has_children())
569  continue;
570 
571  // Now neighbor is on the same refinement level.
572  // Double check.
573  Assert(!cell->neighbor_is_coarser(face_no),
574  ExcInternalError());
575 
576  // If we own both cells only do faces from one side (unless
577  // AssembleFlags says otherwise). Here, we rely on cell
578  // comparison that will look at cell->index().
579  if (own_cell && own_neighbor &&
581  (neighbor < cell))
582  continue;
583 
584  // We only look at faces to ghost on the same level once
585  // (only where own_cell=true and own_neighbor=false)
586  if (!own_cell)
587  continue;
588 
589  // now only one processor assembles faces_to_ghost. We let
590  // the processor with the smaller (level-)subdomain id
591  // assemble the face.
592  if (own_cell && !own_neighbor &&
593  (flags & assemble_ghost_faces_once) &&
594  (neighbor_subdomain_id < current_subdomain_id))
595  continue;
596 
597  const unsigned int neighbor_face_no =
598  periodic_neighbor ?
599  cell->periodic_neighbor_face_no(face_no) :
600  cell->neighbor_face_no(face_no);
601  Assert(periodic_neighbor ||
602  neighbor->face(neighbor_face_no) ==
603  cell->face(face_no),
604  ExcInternalError());
605 
606  face_worker(cell,
607  face_no,
609  neighbor,
610  neighbor_face_no,
612  scratch,
613  copy);
614  }
615  }
616  } // faces
617 
618  // Execute the cell_worker if faces are handled before cells
619  if ((flags & cells_after_faces) &&
620  (((flags & assemble_own_cells) && own_cell) ||
621  ((flags & assemble_ghost_cells) && !own_cell)))
622  cell_worker(cell, scratch, copy);
623  };
624 
625  // Submit to workstream
626  WorkStream::run(begin,
627  end,
628  cell_action,
629  copier,
630  sample_scratch_data,
631  sample_copy_data,
632  queue_length,
633  chunk_size);
634  }
635 
705  template <class CellIteratorType,
706  class ScratchData,
707  class CopyData,
708  class CellIteratorBaseType =
710  void
712  IteratorRange<CellIteratorType> iterator_range,
713  const typename identity<std::function<
714  void(const CellIteratorBaseType &, ScratchData &, CopyData &)>>::type
715  &cell_worker,
716  const typename identity<std::function<void(const CopyData &)>>::type
717  &copier,
718 
719  const ScratchData &sample_scratch_data,
720  const CopyData & sample_copy_data,
721 
722  const AssembleFlags flags = assemble_own_cells,
723 
724  const typename identity<std::function<void(const CellIteratorBaseType &,
725  const unsigned int,
726  ScratchData &,
727  CopyData &)>>::type
728  &boundary_worker = std::function<void(const CellIteratorBaseType &,
729  const unsigned int,
730  ScratchData &,
731  CopyData &)>(),
732 
733  const typename identity<std::function<void(const CellIteratorBaseType &,
734  const unsigned int,
735  const unsigned int,
736  const CellIteratorBaseType &,
737  const unsigned int,
738  const unsigned int,
739  ScratchData &,
740  CopyData &)>>::type
741  &face_worker = std::function<void(const CellIteratorBaseType &,
742  const unsigned int,
743  const unsigned int,
744  const CellIteratorBaseType &,
745  const unsigned int,
746  const unsigned int,
747  ScratchData &,
748  CopyData &)>(),
749 
750  const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
751  const unsigned int chunk_size = 8)
752  {
753  // Call the function above
754  mesh_loop<typename IteratorRange<CellIteratorType>::IteratorOverIterators,
755  ScratchData,
756  CopyData,
757  CellIteratorBaseType>(iterator_range.begin(),
758  iterator_range.end(),
759  cell_worker,
760  copier,
761  sample_scratch_data,
762  sample_copy_data,
763  flags,
764  boundary_worker,
765  face_worker,
766  queue_length,
767  chunk_size);
768  }
769 
829  template <class CellIteratorType,
830  class ScratchData,
831  class CopyData,
832  class MainClass>
833  void
834  mesh_loop(const CellIteratorType & begin,
835  const typename identity<CellIteratorType>::type &end,
836  MainClass & main_class,
837  void (MainClass::*cell_worker)(const CellIteratorType &,
838  ScratchData &,
839  CopyData &),
840  void (MainClass::*copier)(const CopyData &),
841  const ScratchData & sample_scratch_data,
842  const CopyData & sample_copy_data,
843  const AssembleFlags flags = assemble_own_cells,
844  void (MainClass::*boundary_worker)(const CellIteratorType &,
845  const unsigned int,
846  ScratchData &,
847  CopyData &) = nullptr,
848  void (MainClass::*face_worker)(const CellIteratorType &,
849  const unsigned int,
850  const unsigned int,
851  const CellIteratorType &,
852  const unsigned int,
853  const unsigned int,
854  ScratchData &,
855  CopyData &) = nullptr,
856  const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
857  const unsigned int chunk_size = 8)
858  {
859  std::function<void(const CellIteratorType &, ScratchData &, CopyData &)>
860  f_cell_worker;
861 
862  std::function<void(
863  const CellIteratorType &, const unsigned int, ScratchData &, CopyData &)>
864  f_boundary_worker;
865 
866  std::function<void(const CellIteratorType &,
867  const unsigned int,
868  const unsigned int,
869  const CellIteratorType &,
870  const unsigned int,
871  const unsigned int,
872  ScratchData &,
873  CopyData &)>
874  f_face_worker;
875 
876  if (cell_worker != nullptr)
877  f_cell_worker = [&main_class,
878  cell_worker](const CellIteratorType &cell_iterator,
879  ScratchData & scratch_data,
880  CopyData & copy_data) {
881  (main_class.*cell_worker)(cell_iterator, scratch_data, copy_data);
882  };
883 
884  if (boundary_worker != nullptr)
885  f_boundary_worker =
886  [&main_class, boundary_worker](const CellIteratorType &cell_iterator,
887  const unsigned int face_no,
888  ScratchData & scratch_data,
889  CopyData & copy_data) {
890  (main_class.*
891  boundary_worker)(cell_iterator, face_no, scratch_data, copy_data);
892  };
893 
894  if (face_worker != nullptr)
895  f_face_worker = [&main_class,
896  face_worker](const CellIteratorType &cell_iterator_1,
897  const unsigned int face_index_1,
898  const unsigned int subface_index_1,
899  const CellIteratorType &cell_iterator_2,
900  const unsigned int face_index_2,
901  const unsigned int subface_index_2,
902  ScratchData & scratch_data,
903  CopyData & copy_data) {
904  (main_class.*face_worker)(cell_iterator_1,
905  face_index_1,
906  subface_index_1,
907  cell_iterator_2,
908  face_index_2,
909  subface_index_2,
910  scratch_data,
911  copy_data);
912  };
913 
914  mesh_loop(begin,
915  end,
916  f_cell_worker,
917  [&main_class, copier](const CopyData &copy_data) {
918  (main_class.*copier)(copy_data);
919  },
920  sample_scratch_data,
921  sample_copy_data,
922  flags,
923  f_boundary_worker,
924  f_face_worker,
925  queue_length,
926  chunk_size);
927  }
928 
1008  template <class CellIteratorType,
1009  class ScratchData,
1010  class CopyData,
1011  class MainClass,
1012  class CellIteratorBaseType =
1014  void
1016  MainClass & main_class,
1017  void (MainClass::*cell_worker)(const CellIteratorBaseType &,
1018  ScratchData &,
1019  CopyData &),
1020  void (MainClass::*copier)(const CopyData &),
1021  const ScratchData & sample_scratch_data,
1022  const CopyData & sample_copy_data,
1023  const AssembleFlags flags = assemble_own_cells,
1024  void (MainClass::*boundary_worker)(const CellIteratorBaseType &,
1025  const unsigned int,
1026  ScratchData &,
1027  CopyData &) = nullptr,
1028  void (MainClass::*face_worker)(const CellIteratorBaseType &,
1029  const unsigned int,
1030  const unsigned int,
1031  const CellIteratorBaseType &,
1032  const unsigned int,
1033  const unsigned int,
1034  ScratchData &,
1035  CopyData &) = nullptr,
1036  const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
1037  const unsigned int chunk_size = 8)
1038  {
1039  // Call the function above
1040  mesh_loop<typename IteratorRange<CellIteratorType>::IteratorOverIterators,
1041  ScratchData,
1042  CopyData,
1043  MainClass,
1044  CellIteratorBaseType>(iterator_range.begin(),
1045  iterator_range.end(),
1046  main_class,
1047  cell_worker,
1048  copier,
1049  sample_scratch_data,
1050  sample_copy_data,
1051  flags,
1052  boundary_worker,
1053  face_worker,
1054  queue_length,
1055  chunk_size);
1056  }
1057 } // namespace MeshWorker
1058 
1060 
1061 #endif
static const unsigned int invalid_unsigned_int
Definition: types.h:196
IteratorOverIterators begin()
const types::subdomain_id invalid_subdomain_id
Definition: types.h:276
typename CellIteratorBaseType< CellIteratorType >::type type
Definition: mesh_loop.h:98
std::function< void(const CopyData &)> CopierFunctionType
Definition: mesh_loop.h:114
std::function< void(const CellIteratorBaseType &, const unsigned int, ScratchData &, CopyData &)> BoundaryWorkerFunctionType
Definition: mesh_loop.h:124
void mesh_loop(const CellIteratorType &begin, const CellIteratorType &end, const CellWorkerFunctionType &cell_worker, const CopierType &copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const AssembleFlags flags=assemble_own_cells, const BoundaryWorkerFunctionType &boundary_worker=BoundaryWorkerFunctionType(), const FaceWorkerFunctionType &face_worker=FaceWorkerFunctionType(), const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: mesh_loop.h:282
typename CellIteratorBaseType< CellIteratorType >::type type
Definition: mesh_loop.h:79
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1465
std::function< void(const CellIteratorBaseType &, const unsigned int, const unsigned int, const CellIteratorBaseType &, const unsigned int, const unsigned int, ScratchData &, CopyData &)> FaceWorkerFunctionType
Definition: mesh_loop.h:138
bool is_active_iterator(const DI &)
Definition: loop.h:49
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
VectorType::value_type * end(VectorType &V)
std::function< void(const CellIteratorBaseType &, ScratchData &, CopyData &)> CellWorkerFunctionType
Definition: mesh_loop.h:108
const types::subdomain_id artificial_subdomain_id
Definition: types.h:293
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
VectorType::value_type * begin(VectorType &V)
constexpr int chunk_size
Definition: cuda_size.h:35
void run(const std::vector< std::vector< Iterator >> &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1337
static unsigned int n_threads()
void copy(const T *begin, const T *end, U *dest)
void cell_action(ITERATOR cell, DoFInfoBox< dim, DOFINFO > &dof_info, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, const LoopControl &loop_control)
Definition: loop.h:191
static ::ExceptionBase & ExcInternalError()
IteratorOverIterators end() const