Reference documentation for deal.II version Git 292c2606a1 2021-03-03 00:48:14 +0100
\(\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 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii_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 
274  template <class CellIteratorType,
275  class ScratchData,
276  class CopyData,
277  class CellIteratorBaseType =
279  void
281 #ifdef DOXYGEN
282  const CellIteratorType &begin,
283  const CellIteratorType &end,
284 
285  const CellWorkerFunctionType &cell_worker,
286  const CopierType & copier,
287 
288  const ScratchData &sample_scratch_data,
289  const CopyData & sample_copy_data,
290 
291  const AssembleFlags flags = assemble_own_cells,
292 
293  const BoundaryWorkerFunctionType &boundary_worker =
295 
296  const FaceWorkerFunctionType &face_worker = FaceWorkerFunctionType(),
297  const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
298  const unsigned int chunk_size = 8
299 #else
300  const CellIteratorType & begin,
301  const typename identity<CellIteratorType>::type &end,
302 
303  const typename identity<std::function<
304  void(const CellIteratorBaseType &, ScratchData &, CopyData &)>>::type
305  &cell_worker,
306  const typename identity<std::function<void(const CopyData &)>>::type
307  &copier,
308 
309  const ScratchData &sample_scratch_data,
310  const CopyData & sample_copy_data,
311 
312  const AssembleFlags flags = assemble_own_cells,
313 
314  const typename identity<std::function<void(const CellIteratorBaseType &,
315  const unsigned int,
316  ScratchData &,
317  CopyData &)>>::type
318  &boundary_worker = std::function<void(const CellIteratorBaseType &,
319  const unsigned int,
320  ScratchData &,
321  CopyData &)>(),
322 
323  const typename identity<std::function<void(const CellIteratorBaseType &,
324  const unsigned int,
325  const unsigned int,
326  const CellIteratorBaseType &,
327  const unsigned int,
328  const unsigned int,
329  ScratchData &,
330  CopyData &)>>::type
331  &face_worker = std::function<void(const CellIteratorBaseType &,
332  const unsigned int,
333  const unsigned int,
334  const CellIteratorBaseType &,
335  const unsigned int,
336  const unsigned int,
337  ScratchData &,
338  CopyData &)>(),
339 
340  const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
341  const unsigned int chunk_size = 8
342 #endif
343  )
344  {
345  Assert(
346  (!cell_worker) == !(flags & work_on_cells),
347  ExcMessage(
348  "If you provide a cell worker function, you also need to request "
349  "that work should be done on cells by setting the 'work_on_cells' flag. "
350  "Conversely, if you don't provide a cell worker function, you "
351  "cannot set the 'work_on_cells' flag. One of these two "
352  "conditions is not satisfied."));
353 
358  ExcMessage(
359  "If you provide a face worker function, you also need to request "
360  "that work should be done on interior faces by setting either the "
361  "'assemble_own_interior_faces_once' flag or the "
362  "'assemble_own_interior_faces_both' flag. "
363  "Conversely, if you don't provide a face worker function, you "
364  "cannot set either of these two flags. One of these two "
365  "conditions is not satisfied."));
366 
369  ExcMessage(
370  "You can only 'specify assemble_ghost_faces_once' "
371  "OR 'assemble_ghost_faces_both', but not both of these flags."));
372 
373  Assert(
374  !(flags & cells_after_faces) ||
376  ExcMessage(
377  "The option 'cells_after_faces' only makes sense if you assemble on cells."));
378 
379  Assert(
380  (!face_worker) == !(flags & work_on_faces),
381  ExcMessage(
382  "If you provide a face worker function, you also need to request "
383  "that work should be done on faces by setting the 'work_on_faces' flag. "
384  "Conversely, if you don't provide a face worker function, you "
385  "cannot set the 'work_on_faces' flag. One of these two "
386  "conditions is not satisfied."));
387 
388  Assert(
389  (!boundary_worker) == !(flags & assemble_boundary_faces),
390  ExcMessage(
391  "If you provide a boundary face worker function, you also need to request "
392  "that work should be done on boundary faces by setting the 'assemble_boundary_faces' flag. "
393  "Conversely, if you don't provide a boundary face worker function, you "
394  "cannot set the 'assemble_boundary_faces' flag. One of these two "
395  "conditions is not satisfied."));
396 
397  auto cell_action = [&](const CellIteratorBaseType &cell,
398  ScratchData & scratch,
399  CopyData & copy) {
400  // First reset the CopyData class to the empty copy_data given by the
401  // user.
402  copy = sample_copy_data;
403 
404  // Store the dimension in which we are working for later use
405  const auto dim = cell->get_triangulation().dimension;
406 
407  const bool ignore_subdomain =
408  (cell->get_triangulation().locally_owned_subdomain() ==
410 
411  types::subdomain_id current_subdomain_id =
412  (cell->is_level_cell() ? cell->level_subdomain_id() :
413  cell->subdomain_id());
414 
415  const bool own_cell =
416  ignore_subdomain ||
417  (current_subdomain_id ==
418  cell->get_triangulation().locally_owned_subdomain());
419 
420  if ((!ignore_subdomain) &&
421  (current_subdomain_id == numbers::artificial_subdomain_id))
422  return;
423 
424  if (!(flags & (cells_after_faces)) &&
425  (((flags & (assemble_own_cells)) && own_cell) ||
426  ((flags & assemble_ghost_cells) && !own_cell)))
427  cell_worker(cell, scratch, copy);
428 
429  if (flags & (work_on_faces | work_on_boundary))
430  for (const unsigned int face_no : cell->face_indices())
431  {
432  if (cell->at_boundary(face_no) &&
433  !cell->has_periodic_neighbor(face_no))
434  {
435  // only integrate boundary faces of own cells
436  if ((flags & assemble_boundary_faces) && own_cell)
437  boundary_worker(cell, face_no, scratch, copy);
438  }
439  else
440  {
441  // interior face, potentially assemble
443  neighbor = cell->neighbor_or_periodic_neighbor(face_no);
444 
445  types::subdomain_id neighbor_subdomain_id =
447  if (neighbor->is_level_cell())
448  neighbor_subdomain_id = neighbor->level_subdomain_id();
449  // subdomain id is only valid for active cells
450  else if (neighbor->is_active())
451  neighbor_subdomain_id = neighbor->subdomain_id();
452 
453  const bool own_neighbor =
454  ignore_subdomain ||
455  (neighbor_subdomain_id ==
456  cell->get_triangulation().locally_owned_subdomain());
457 
458  // skip all faces between two ghost cells
459  if (!own_cell && !own_neighbor)
460  continue;
461 
462  // skip if the user doesn't want faces between own cells
463  if (own_cell && own_neighbor &&
466  continue;
467 
468  // skip face to ghost
469  if (own_cell != own_neighbor &&
470  !(flags &
472  continue;
473 
474  // Deal with refinement edges from the refined side. Assuming
475  // one-irregular meshes, this situation should only occur if
476  // both cells are active.
477  const bool periodic_neighbor =
478  cell->has_periodic_neighbor(face_no);
479 
480  if (dim > 1 && ((!periodic_neighbor &&
481  cell->neighbor_is_coarser(face_no) &&
482  neighbor->is_active()) ||
483  (periodic_neighbor &&
484  cell->periodic_neighbor_is_coarser(face_no) &&
485  neighbor->is_active())))
486  {
487  Assert(cell->is_active(), ExcInternalError());
488 
489  // skip if only one processor needs to assemble the face
490  // to a ghost cell and the fine cell is not ours.
491  if (!own_cell && (flags & assemble_ghost_faces_once))
492  continue;
493 
494  const std::pair<unsigned int, unsigned int>
495  neighbor_face_no =
496  periodic_neighbor ?
497  cell->periodic_neighbor_of_coarser_periodic_neighbor(
498  face_no) :
499  cell->neighbor_of_coarser_neighbor(face_no);
500 
501  face_worker(cell,
502  face_no,
504  neighbor,
505  neighbor_face_no.first,
506  neighbor_face_no.second,
507  scratch,
508  copy);
509 
511  {
512  // If own faces are to be assembled from both sides,
513  // call the faceworker again with swapped arguments.
514  // This is because we won't be looking at an adaptively
515  // refined edge coming from the other side.
516  face_worker(neighbor,
517  neighbor_face_no.first,
518  neighbor_face_no.second,
519  cell,
520  face_no,
522  scratch,
523  copy);
524  }
525  }
526  else if (dim == 1 && cell->level() > neighbor->level())
527  {
528  // In one dimension, there is no other check to do
529  const unsigned int neighbor_face_no =
530  periodic_neighbor ?
531  cell->periodic_neighbor_face_no(face_no) :
532  cell->neighbor_face_no(face_no);
533  Assert(periodic_neighbor ||
534  neighbor->face(neighbor_face_no) ==
535  cell->face(face_no),
536  ExcInternalError());
537 
538  face_worker(cell,
539  face_no,
541  neighbor,
542  neighbor_face_no,
544  scratch,
545  copy);
546 
548  {
549  // If own faces are to be assembled from both sides,
550  // call the faceworker again with swapped arguments.
551  face_worker(neighbor,
552  neighbor_face_no,
554  cell,
555  face_no,
557  scratch,
558  copy);
559  }
560  }
561  else
562  {
563  // If iterator is active and neighbor is refined, skip
564  // internal face.
565  if (::internal::is_active_iterator(cell) &&
566  neighbor->has_children())
567  continue;
568 
569  // Now neighbor is on the same refinement level.
570  // Double check.
571  Assert(!cell->neighbor_is_coarser(face_no),
572  ExcInternalError());
573 
574  // If we own both cells only do faces from one side (unless
575  // AssembleFlags says otherwise). Here, we rely on cell
576  // comparison that will look at cell->index().
577  if (own_cell && own_neighbor &&
579  (neighbor < cell))
580  continue;
581 
582  // We only look at faces to ghost on the same level once
583  // (only where own_cell=true and own_neighbor=false)
584  if (!own_cell)
585  continue;
586 
587  // now only one processor assembles faces_to_ghost. We let
588  // the processor with the smaller (level-)subdomain id
589  // assemble the face.
590  if (own_cell && !own_neighbor &&
591  (flags & assemble_ghost_faces_once) &&
592  (neighbor_subdomain_id < current_subdomain_id))
593  continue;
594 
595  const unsigned int neighbor_face_no =
596  periodic_neighbor ?
597  cell->periodic_neighbor_face_no(face_no) :
598  cell->neighbor_face_no(face_no);
599  Assert(periodic_neighbor ||
600  neighbor->face(neighbor_face_no) ==
601  cell->face(face_no),
602  ExcInternalError());
603 
604  face_worker(cell,
605  face_no,
607  neighbor,
608  neighbor_face_no,
610  scratch,
611  copy);
612  }
613  }
614  } // faces
615 
616  // Execute the cell_worker if faces are handled before cells
617  if ((flags & cells_after_faces) &&
618  (((flags & assemble_own_cells) && own_cell) ||
619  ((flags & assemble_ghost_cells) && !own_cell)))
620  cell_worker(cell, scratch, copy);
621  };
622 
623  // Submit to workstream
624  WorkStream::run(begin,
625  end,
626  cell_action,
627  copier,
628  sample_scratch_data,
629  sample_copy_data,
630  queue_length,
631  chunk_size);
632  }
633 
703  template <class CellIteratorType,
704  class ScratchData,
705  class CopyData,
706  class CellIteratorBaseType =
708  void
710  IteratorRange<CellIteratorType> iterator_range,
711  const typename identity<std::function<
712  void(const CellIteratorBaseType &, ScratchData &, CopyData &)>>::type
713  &cell_worker,
714  const typename identity<std::function<void(const CopyData &)>>::type
715  &copier,
716 
717  const ScratchData &sample_scratch_data,
718  const CopyData & sample_copy_data,
719 
720  const AssembleFlags flags = assemble_own_cells,
721 
722  const typename identity<std::function<void(const CellIteratorBaseType &,
723  const unsigned int,
724  ScratchData &,
725  CopyData &)>>::type
726  &boundary_worker = std::function<void(const CellIteratorBaseType &,
727  const unsigned int,
728  ScratchData &,
729  CopyData &)>(),
730 
731  const typename identity<std::function<void(const CellIteratorBaseType &,
732  const unsigned int,
733  const unsigned int,
734  const CellIteratorBaseType &,
735  const unsigned int,
736  const unsigned int,
737  ScratchData &,
738  CopyData &)>>::type
739  &face_worker = std::function<void(const CellIteratorBaseType &,
740  const unsigned int,
741  const unsigned int,
742  const CellIteratorBaseType &,
743  const unsigned int,
744  const unsigned int,
745  ScratchData &,
746  CopyData &)>(),
747 
748  const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
749  const unsigned int chunk_size = 8)
750  {
751  // Call the function above
752  mesh_loop<typename IteratorRange<CellIteratorType>::IteratorOverIterators,
753  ScratchData,
754  CopyData,
755  CellIteratorBaseType>(iterator_range.begin(),
756  iterator_range.end(),
757  cell_worker,
758  copier,
759  sample_scratch_data,
760  sample_copy_data,
761  flags,
762  boundary_worker,
763  face_worker,
764  queue_length,
765  chunk_size);
766  }
767 
827  template <class CellIteratorType,
828  class ScratchData,
829  class CopyData,
830  class MainClass>
831  void
832  mesh_loop(const CellIteratorType & begin,
833  const typename identity<CellIteratorType>::type &end,
834  MainClass & main_class,
835  void (MainClass::*cell_worker)(const CellIteratorType &,
836  ScratchData &,
837  CopyData &),
838  void (MainClass::*copier)(const CopyData &),
839  const ScratchData & sample_scratch_data,
840  const CopyData & sample_copy_data,
841  const AssembleFlags flags = assemble_own_cells,
842  void (MainClass::*boundary_worker)(const CellIteratorType &,
843  const unsigned int,
844  ScratchData &,
845  CopyData &) = nullptr,
846  void (MainClass::*face_worker)(const CellIteratorType &,
847  const unsigned int,
848  const unsigned int,
849  const CellIteratorType &,
850  const unsigned int,
851  const unsigned int,
852  ScratchData &,
853  CopyData &) = nullptr,
854  const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
855  const unsigned int chunk_size = 8)
856  {
857  std::function<void(const CellIteratorType &, ScratchData &, CopyData &)>
858  f_cell_worker;
859 
860  std::function<void(
861  const CellIteratorType &, const unsigned int, ScratchData &, CopyData &)>
862  f_boundary_worker;
863 
864  std::function<void(const CellIteratorType &,
865  const unsigned int,
866  const unsigned int,
867  const CellIteratorType &,
868  const unsigned int,
869  const unsigned int,
870  ScratchData &,
871  CopyData &)>
872  f_face_worker;
873 
874  if (cell_worker != nullptr)
875  f_cell_worker = [&main_class,
876  cell_worker](const CellIteratorType &cell_iterator,
877  ScratchData & scratch_data,
878  CopyData & copy_data) {
879  (main_class.*cell_worker)(cell_iterator, scratch_data, copy_data);
880  };
881 
882  if (boundary_worker != nullptr)
883  f_boundary_worker =
884  [&main_class, boundary_worker](const CellIteratorType &cell_iterator,
885  const unsigned int face_no,
886  ScratchData & scratch_data,
887  CopyData & copy_data) {
888  (main_class.*
889  boundary_worker)(cell_iterator, face_no, scratch_data, copy_data);
890  };
891 
892  if (face_worker != nullptr)
893  f_face_worker = [&main_class,
894  face_worker](const CellIteratorType &cell_iterator_1,
895  const unsigned int face_index_1,
896  const unsigned int subface_index_1,
897  const CellIteratorType &cell_iterator_2,
898  const unsigned int face_index_2,
899  const unsigned int subface_index_2,
900  ScratchData & scratch_data,
901  CopyData & copy_data) {
902  (main_class.*face_worker)(cell_iterator_1,
903  face_index_1,
904  subface_index_1,
905  cell_iterator_2,
906  face_index_2,
907  subface_index_2,
908  scratch_data,
909  copy_data);
910  };
911 
912  mesh_loop(begin,
913  end,
914  f_cell_worker,
915  [&main_class, copier](const CopyData &copy_data) {
916  (main_class.*copier)(copy_data);
917  },
918  sample_scratch_data,
919  sample_copy_data,
920  flags,
921  f_boundary_worker,
922  f_face_worker,
923  queue_length,
924  chunk_size);
925  }
926 
1006  template <class CellIteratorType,
1007  class ScratchData,
1008  class CopyData,
1009  class MainClass,
1010  class CellIteratorBaseType =
1012  void
1014  MainClass & main_class,
1015  void (MainClass::*cell_worker)(const CellIteratorBaseType &,
1016  ScratchData &,
1017  CopyData &),
1018  void (MainClass::*copier)(const CopyData &),
1019  const ScratchData & sample_scratch_data,
1020  const CopyData & sample_copy_data,
1021  const AssembleFlags flags = assemble_own_cells,
1022  void (MainClass::*boundary_worker)(const CellIteratorBaseType &,
1023  const unsigned int,
1024  ScratchData &,
1025  CopyData &) = nullptr,
1026  void (MainClass::*face_worker)(const CellIteratorBaseType &,
1027  const unsigned int,
1028  const unsigned int,
1029  const CellIteratorBaseType &,
1030  const unsigned int,
1031  const unsigned int,
1032  ScratchData &,
1033  CopyData &) = nullptr,
1034  const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
1035  const unsigned int chunk_size = 8)
1036  {
1037  // Call the function above
1038  mesh_loop<typename IteratorRange<CellIteratorType>::IteratorOverIterators,
1039  ScratchData,
1040  CopyData,
1041  MainClass,
1042  CellIteratorBaseType>(iterator_range.begin(),
1043  iterator_range.end(),
1044  main_class,
1045  cell_worker,
1046  copier,
1047  sample_scratch_data,
1048  sample_copy_data,
1049  flags,
1050  boundary_worker,
1051  face_worker,
1052  queue_length,
1053  chunk_size);
1054  }
1055 } // namespace MeshWorker
1056 
1058 
1059 #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:280
typename CellIteratorBaseType< CellIteratorType >::type type
Definition: mesh_loop.h:79
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1466
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:394
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:393
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