Reference documentation for deal.II version GIT 29f9da0a34 2023-12-07 10:00:01+00:00
\(\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 - 2023 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 <typename CellIteratorType>
56  {
60  using type = CellIteratorType;
61  };
62 
70  template <typename 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 <typename 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
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 
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 <typename CellIteratorType,
277  class ScratchData,
278  class CopyData,
279  typename 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,
304 
305  const std_cxx20::type_identity_t<std::function<
306  void(const CellIteratorBaseType &, ScratchData &, CopyData &)>>
307  &cell_worker,
308  const std_cxx20::type_identity_t<std::function<void(const CopyData &)>>
309  &copier,
310 
311  const ScratchData &sample_scratch_data,
312  const CopyData &sample_copy_data,
313 
314  const AssembleFlags flags = assemble_own_cells,
315 
317  std::function<void(const CellIteratorBaseType &,
318  const unsigned int,
319  ScratchData &,
320  CopyData &)>> &boundary_worker =
321  std::function<void(const CellIteratorBaseType &,
322  const unsigned int,
323  ScratchData &,
324  CopyData &)>(),
325 
327  std::function<void(const CellIteratorBaseType &,
328  const unsigned int,
329  const unsigned int,
330  const CellIteratorBaseType &,
331  const unsigned int,
332  const unsigned int,
333  ScratchData &,
334  CopyData &)>> &face_worker =
335  std::function<void(const CellIteratorBaseType &,
336  const unsigned int,
337  const unsigned int,
338  const CellIteratorBaseType &,
339  const unsigned int,
340  const unsigned int,
341  ScratchData &,
342  CopyData &)>(),
343 
344  const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
345  const unsigned int chunk_size = 8
346 #endif
347  )
348  {
349  Assert(
350  (!cell_worker) == !(flags & work_on_cells),
351  ExcMessage(
352  "If you provide a cell worker function, you also need to request "
353  "that work should be done on cells by setting the 'work_on_cells' flag. "
354  "Conversely, if you don't provide a cell worker function, you "
355  "cannot set the 'work_on_cells' flag. One of these two "
356  "conditions is not satisfied."));
357 
362  ExcMessage(
363  "If you provide a face worker function, you also need to request "
364  "that work should be done on interior faces by setting either the "
365  "'assemble_own_interior_faces_once' flag or the "
366  "'assemble_own_interior_faces_both' flag. "
367  "Conversely, if you don't provide a face worker function, you "
368  "cannot set either of these two flags. One of these two "
369  "conditions is not satisfied."));
370 
373  ExcMessage(
374  "You can only 'specify assemble_ghost_faces_once' "
375  "OR 'assemble_ghost_faces_both', but not both of these flags."));
376 
377  Assert(
378  !(flags & cells_after_faces) ||
380  ExcMessage(
381  "The option 'cells_after_faces' only makes sense if you assemble on cells."));
382 
383  Assert(
384  (!face_worker) == !(flags & work_on_faces),
385  ExcMessage(
386  "If you provide a face worker function, you also need to request "
387  "that work should be done on faces by setting the 'work_on_faces' flag. "
388  "Conversely, if you don't provide a face worker function, you "
389  "cannot set the 'work_on_faces' flag. One of these two "
390  "conditions is not satisfied."));
391 
392  Assert(
393  (!boundary_worker) == !(flags & assemble_boundary_faces),
394  ExcMessage(
395  "If you provide a boundary face worker function, you also need to request "
396  "that work should be done on boundary faces by setting the 'assemble_boundary_faces' flag. "
397  "Conversely, if you don't provide a boundary face worker function, you "
398  "cannot set the 'assemble_boundary_faces' flag. One of these two "
399  "conditions is not satisfied."));
400 
401  auto cell_action = [&](const CellIteratorBaseType &cell,
402  ScratchData &scratch,
403  CopyData &copy) {
404  // First reset the CopyData class to the empty copy_data given by the
405  // user.
406  copy = sample_copy_data;
407 
408  // Store the dimension in which we are working for later use
409  const auto dim = cell->get_triangulation().dimension;
410 
411  const bool ignore_subdomain =
412  (cell->get_triangulation().locally_owned_subdomain() ==
414 
415  types::subdomain_id current_subdomain_id =
416  (cell->is_level_cell() ? cell->level_subdomain_id() :
417  cell->subdomain_id());
418 
419  const bool own_cell =
420  ignore_subdomain ||
421  (current_subdomain_id ==
422  cell->get_triangulation().locally_owned_subdomain());
423 
424  if ((!ignore_subdomain) &&
425  (current_subdomain_id == numbers::artificial_subdomain_id))
426  return;
427 
428  if (!(flags & (cells_after_faces)) &&
429  (((flags & (assemble_own_cells)) && own_cell) ||
430  ((flags & assemble_ghost_cells) && !own_cell)))
431  cell_worker(cell, scratch, copy);
432 
433  if (flags & (work_on_faces | work_on_boundary))
434  for (const unsigned int face_no : cell->face_indices())
435  {
436  if (cell->at_boundary(face_no) &&
437  !cell->has_periodic_neighbor(face_no))
438  {
439  // only integrate boundary faces of own cells
440  if ((flags & assemble_boundary_faces) && own_cell)
441  boundary_worker(cell, face_no, scratch, copy);
442  }
443  else
444  {
445  // interior face, potentially assemble
447  neighbor = cell->neighbor_or_periodic_neighbor(face_no);
448 
449  types::subdomain_id neighbor_subdomain_id =
451  if (neighbor->is_level_cell())
452  neighbor_subdomain_id = neighbor->level_subdomain_id();
453  // subdomain id is only valid for active cells
454  else if (neighbor->is_active())
455  neighbor_subdomain_id = neighbor->subdomain_id();
456 
457  const bool own_neighbor =
458  ignore_subdomain ||
459  (neighbor_subdomain_id ==
460  cell->get_triangulation().locally_owned_subdomain());
461 
462  // skip all faces between two ghost cells
463  if (!own_cell && !own_neighbor)
464  continue;
465 
466  // skip if the user doesn't want faces between own cells
467  if (own_cell && own_neighbor &&
470  continue;
471 
472  // skip face to ghost
473  if (own_cell != own_neighbor &&
474  !(flags &
476  continue;
477 
478  // Deal with refinement edges from the refined side. Assuming
479  // one-irregular meshes, this situation should only occur if
480  // both cells are active.
481  const bool periodic_neighbor =
482  cell->has_periodic_neighbor(face_no);
483 
484  if (dim > 1 && ((!periodic_neighbor &&
485  cell->neighbor_is_coarser(face_no) &&
486  neighbor->is_active()) ||
487  (periodic_neighbor &&
488  cell->periodic_neighbor_is_coarser(face_no) &&
489  neighbor->is_active())))
490  {
491  Assert(cell->is_active(), ExcInternalError());
492 
493  // skip if only one processor needs to assemble the face
494  // to a ghost cell and the fine cell is not ours.
495  if (!own_cell && (flags & assemble_ghost_faces_once))
496  continue;
497 
498  const std::pair<unsigned int, unsigned int>
499  neighbor_face_no =
500  periodic_neighbor ?
501  cell->periodic_neighbor_of_coarser_periodic_neighbor(
502  face_no) :
503  cell->neighbor_of_coarser_neighbor(face_no);
504 
505  face_worker(cell,
506  face_no,
508  neighbor,
509  neighbor_face_no.first,
510  neighbor_face_no.second,
511  scratch,
512  copy);
513 
515  {
516  // If own faces are to be assembled from both sides,
517  // call the faceworker again with swapped arguments.
518  // This is because we won't be looking at an adaptively
519  // refined edge coming from the other side.
520  face_worker(neighbor,
521  neighbor_face_no.first,
522  neighbor_face_no.second,
523  cell,
524  face_no,
526  scratch,
527  copy);
528  }
529  }
530  else if (dim == 1 && cell->level() > neighbor->level())
531  {
532  // In one dimension, there is no other check to do
533  const unsigned int neighbor_face_no =
534  periodic_neighbor ?
535  cell->periodic_neighbor_face_no(face_no) :
536  cell->neighbor_face_no(face_no);
537  Assert(periodic_neighbor ||
538  neighbor->face(neighbor_face_no) ==
539  cell->face(face_no),
540  ExcInternalError());
541 
542  face_worker(cell,
543  face_no,
545  neighbor,
546  neighbor_face_no,
548  scratch,
549  copy);
550 
552  {
553  // If own faces are to be assembled from both sides,
554  // call the faceworker again with swapped arguments.
555  face_worker(neighbor,
556  neighbor_face_no,
558  cell,
559  face_no,
561  scratch,
562  copy);
563  }
564  }
565  else
566  {
567  // If iterator is active and neighbor is refined, skip
568  // internal face.
569  if (::internal::is_active_iterator(cell) &&
570  neighbor->has_children())
571  continue;
572 
573  // Now neighbor is on the same refinement level.
574  // Double check.
575  Assert((!periodic_neighbor &&
576  !cell->neighbor_is_coarser(face_no)) ||
577  (periodic_neighbor &&
578  !cell->periodic_neighbor_is_coarser(face_no)),
579  ExcInternalError());
580 
581  // If we own both cells only do faces from one side (unless
582  // AssembleFlags says otherwise). Here, we rely on cell
583  // comparison that will look at cell->index().
584  if (own_cell && own_neighbor &&
586  (neighbor < cell))
587  continue;
588 
589  // We only look at faces to ghost on the same level once
590  // (only where own_cell=true and own_neighbor=false)
591  if (!own_cell)
592  continue;
593 
594  // now only one processor assembles faces_to_ghost. We let
595  // the processor with the smaller (level-)subdomain id
596  // assemble the face.
597  if (own_cell && !own_neighbor &&
598  (flags & assemble_ghost_faces_once) &&
599  (neighbor_subdomain_id < current_subdomain_id))
600  continue;
601 
602  const unsigned int neighbor_face_no =
603  periodic_neighbor ?
604  cell->periodic_neighbor_face_no(face_no) :
605  cell->neighbor_face_no(face_no);
606  Assert(periodic_neighbor ||
607  neighbor->face(neighbor_face_no) ==
608  cell->face(face_no),
609  ExcInternalError());
610 
611  face_worker(cell,
612  face_no,
614  neighbor,
615  neighbor_face_no,
617  scratch,
618  copy);
619  }
620  }
621  } // faces
622 
623  // Execute the cell_worker if faces are handled before cells
624  if ((flags & cells_after_faces) &&
625  (((flags & assemble_own_cells) && own_cell) ||
626  ((flags & assemble_ghost_cells) && !own_cell)))
627  cell_worker(cell, scratch, copy);
628  };
629 
630  // Submit to workstream
632  end,
633  cell_action,
634  copier,
635  sample_scratch_data,
636  sample_copy_data,
637  queue_length,
638  chunk_size);
639  }
640 
710  template <typename CellIteratorType,
711  class ScratchData,
712  class CopyData,
713  typename CellIteratorBaseType =
715  void
717  IteratorRange<CellIteratorType> iterator_range,
718  const std_cxx20::type_identity_t<std::function<
719  void(const CellIteratorBaseType &, ScratchData &, CopyData &)>>
720  &cell_worker,
721  const std_cxx20::type_identity_t<std::function<void(const CopyData &)>>
722  &copier,
723 
724  const ScratchData &sample_scratch_data,
725  const CopyData &sample_copy_data,
726 
727  const AssembleFlags flags = assemble_own_cells,
728 
730  std::function<void(const CellIteratorBaseType &,
731  const unsigned int,
732  ScratchData &,
733  CopyData &)>> &boundary_worker =
734  std::function<void(const CellIteratorBaseType &,
735  const unsigned int,
736  ScratchData &,
737  CopyData &)>(),
738 
740  std::function<void(const CellIteratorBaseType &,
741  const unsigned int,
742  const unsigned int,
743  const CellIteratorBaseType &,
744  const unsigned int,
745  const unsigned int,
746  ScratchData &,
747  CopyData &)>> &face_worker =
748  std::function<void(const CellIteratorBaseType &,
749  const unsigned int,
750  const unsigned int,
751  const CellIteratorBaseType &,
752  const unsigned int,
753  const unsigned int,
754  ScratchData &,
755  CopyData &)>(),
756 
757  const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
758  const unsigned int chunk_size = 8)
759  {
760  // Call the function above
761  mesh_loop<typename IteratorRange<CellIteratorType>::IteratorOverIterators,
762  ScratchData,
763  CopyData,
764  CellIteratorBaseType>(iterator_range.begin(),
765  iterator_range.end(),
766  cell_worker,
767  copier,
768  sample_scratch_data,
769  sample_copy_data,
770  flags,
771  boundary_worker,
772  face_worker,
773  queue_length,
774  chunk_size);
775  }
776 
836  template <typename CellIteratorType,
837  class ScratchData,
838  class CopyData,
839  class MainClass>
840  void
841  mesh_loop(const CellIteratorType &begin,
843  MainClass &main_class,
844  void (MainClass::*cell_worker)(const CellIteratorType &,
845  ScratchData &,
846  CopyData &),
847  void (MainClass::*copier)(const CopyData &),
848  const ScratchData &sample_scratch_data,
849  const CopyData &sample_copy_data,
850  const AssembleFlags flags = assemble_own_cells,
851  void (MainClass::*boundary_worker)(const CellIteratorType &,
852  const unsigned int,
853  ScratchData &,
854  CopyData &) = nullptr,
855  void (MainClass::*face_worker)(const CellIteratorType &,
856  const unsigned int,
857  const unsigned int,
858  const CellIteratorType &,
859  const unsigned int,
860  const unsigned int,
861  ScratchData &,
862  CopyData &) = nullptr,
863  const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
864  const unsigned int chunk_size = 8)
865  {
866  std::function<void(const CellIteratorType &, ScratchData &, CopyData &)>
867  f_cell_worker;
868 
869  std::function<void(
870  const CellIteratorType &, const unsigned int, ScratchData &, CopyData &)>
871  f_boundary_worker;
872 
873  std::function<void(const CellIteratorType &,
874  const unsigned int,
875  const unsigned int,
876  const CellIteratorType &,
877  const unsigned int,
878  const unsigned int,
879  ScratchData &,
880  CopyData &)>
881  f_face_worker;
882 
883  if (cell_worker != nullptr)
884  f_cell_worker = [&main_class,
885  cell_worker](const CellIteratorType &cell_iterator,
886  ScratchData &scratch_data,
887  CopyData &copy_data) {
888  (main_class.*cell_worker)(cell_iterator, scratch_data, copy_data);
889  };
890 
891  if (boundary_worker != nullptr)
892  f_boundary_worker =
893  [&main_class, boundary_worker](const CellIteratorType &cell_iterator,
894  const unsigned int face_no,
895  ScratchData &scratch_data,
896  CopyData &copy_data) {
897  (main_class.*
898  boundary_worker)(cell_iterator, face_no, scratch_data, copy_data);
899  };
900 
901  if (face_worker != nullptr)
902  f_face_worker = [&main_class,
903  face_worker](const CellIteratorType &cell_iterator_1,
904  const unsigned int face_index_1,
905  const unsigned int subface_index_1,
906  const CellIteratorType &cell_iterator_2,
907  const unsigned int face_index_2,
908  const unsigned int subface_index_2,
909  ScratchData &scratch_data,
910  CopyData &copy_data) {
911  (main_class.*face_worker)(cell_iterator_1,
912  face_index_1,
913  subface_index_1,
914  cell_iterator_2,
915  face_index_2,
916  subface_index_2,
917  scratch_data,
918  copy_data);
919  };
920 
921  mesh_loop(
922  begin,
923  end,
924  f_cell_worker,
925  [&main_class, copier](const CopyData &copy_data) {
926  (main_class.*copier)(copy_data);
927  },
928  sample_scratch_data,
929  sample_copy_data,
930  flags,
931  f_boundary_worker,
932  f_face_worker,
933  queue_length,
934  chunk_size);
935  }
936 
1016  template <typename CellIteratorType,
1017  class ScratchData,
1018  class CopyData,
1019  class MainClass,
1020  typename CellIteratorBaseType =
1022  void
1024  MainClass &main_class,
1025  void (MainClass::*cell_worker)(const CellIteratorBaseType &,
1026  ScratchData &,
1027  CopyData &),
1028  void (MainClass::*copier)(const CopyData &),
1029  const ScratchData &sample_scratch_data,
1030  const CopyData &sample_copy_data,
1031  const AssembleFlags flags = assemble_own_cells,
1032  void (MainClass::*boundary_worker)(const CellIteratorBaseType &,
1033  const unsigned int,
1034  ScratchData &,
1035  CopyData &) = nullptr,
1036  void (MainClass::*face_worker)(const CellIteratorBaseType &,
1037  const unsigned int,
1038  const unsigned int,
1039  const CellIteratorBaseType &,
1040  const unsigned int,
1041  const unsigned int,
1042  ScratchData &,
1043  CopyData &) = nullptr,
1044  const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
1045  const unsigned int chunk_size = 8)
1046  {
1047  // Call the function above
1048  mesh_loop<typename IteratorRange<CellIteratorType>::IteratorOverIterators,
1049  ScratchData,
1050  CopyData,
1051  MainClass,
1052  CellIteratorBaseType>(iterator_range.begin(),
1053  iterator_range.end(),
1054  main_class,
1055  cell_worker,
1056  copier,
1057  sample_scratch_data,
1058  sample_copy_data,
1059  flags,
1060  boundary_worker,
1061  face_worker,
1062  queue_length,
1063  chunk_size);
1064  }
1065 } // namespace MeshWorker
1066 
1068 
1069 #endif
IteratorOverIterators end() const
IteratorOverIterators begin()
static unsigned int n_threads()
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcMessage(std::string arg1)
void cell_action(IteratorType 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:195
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
constexpr int chunk_size
Definition: cuda_size.h:35
std::function< void(const CellIteratorBaseType &, const unsigned int, ScratchData &, CopyData &)> BoundaryWorkerFunctionType
Definition: mesh_loop.h:124
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
std::function< void(const CellIteratorBaseType &, ScratchData &, CopyData &)> CellWorkerFunctionType
Definition: mesh_loop.h:108
std::function< void(const CopyData &)> CopierFunctionType
Definition: mesh_loop.h:114
@ assemble_boundary_faces
@ assemble_own_interior_faces_once
@ assemble_ghost_faces_both
@ assemble_own_interior_faces_both
@ assemble_ghost_cells
@ assemble_ghost_faces_once
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
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:1270
void copy(const T *begin, const T *end, U *dest)
bool is_active_iterator(const DI &)
Definition: loop.h:49
const types::subdomain_id artificial_subdomain_id
Definition: types.h:363
const types::subdomain_id invalid_subdomain_id
Definition: types.h:342
static const unsigned int invalid_unsigned_int
Definition: types.h:221
typename type_identity< T >::type type_identity_t
Definition: type_traits.h:96
unsigned int subdomain_id
Definition: types.h:44
typename CellIteratorBaseType< CellIteratorType >::type type
Definition: mesh_loop.h:98
typename CellIteratorBaseType< CellIteratorType >::type type
Definition: mesh_loop.h:79