Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
mesh_loop.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 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_mesh_worker_mesh_loop_h
18 #define dealii_mesh_worker_mesh_loop_h
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/template_constraints.h>
23 #include <deal.II/base/work_stream.h>
24 
25 #include <deal.II/grid/filtered_iterator.h>
26 #include <deal.II/grid/tria.h>
27 
28 #include <deal.II/meshworker/assemble_flags.h>
29 #include <deal.II/meshworker/dof_info.h>
30 #include <deal.II/meshworker/integration_info.h>
31 #include <deal.II/meshworker/local_integrator.h>
32 #include <deal.II/meshworker/loop.h>
33 
34 #include <functional>
35 #include <type_traits>
36 
37 DEAL_II_NAMESPACE_OPEN
38 
39 template <typename>
40 class TriaActiveIterator;
41 
42 namespace MeshWorker
43 {
44  namespace internal
45  {
50  template <class CellIteratorType>
52  {
56  using type = CellIteratorType;
57  };
58 
66  template <class CellIteratorType>
67  struct CellIteratorBaseType<IteratorOverIterators<CellIteratorType>>
68  {
72  // Since we can have filtered iterators and the like as template
73  // arguments, we recursivelyremove the template layers to retrieve the
74  // underlying iterator type.
76  };
77 
86  template <class CellIteratorType>
87  struct CellIteratorBaseType<FilteredIterator<CellIteratorType>>
88  {
92  // Since we can have nested filtered iterators, we recursively
93  // remove the template layers to retrieve the underlying iterator type.
95  };
96  } // namespace internal
97 
173  template <class CellIteratorType,
174  class ScratchData,
175  class CopyData,
176  class CellIteratorBaseType =
178  void
180  const CellIteratorType & begin,
181  const typename identity<CellIteratorType>::type &end,
182 
183  const typename identity<std::function<
184  void(const CellIteratorBaseType &, ScratchData &, CopyData &)>>::type
185  &cell_worker,
186  const typename identity<std::function<void(const CopyData &)>>::type
187  &copier,
188 
189  const ScratchData &sample_scratch_data,
190  const CopyData & sample_copy_data,
191 
192  const AssembleFlags flags = assemble_own_cells,
193 
194  const typename identity<std::function<void(const CellIteratorBaseType &,
195  const unsigned int,
196  ScratchData &,
197  CopyData &)>>::type
198  &boundary_worker = std::function<void(const CellIteratorBaseType &,
199  const unsigned int,
200  ScratchData &,
201  CopyData &)>(),
202 
203  const typename identity<std::function<void(const CellIteratorBaseType &,
204  const unsigned int,
205  const unsigned int,
206  const CellIteratorBaseType &,
207  const unsigned int,
208  const unsigned int,
209  ScratchData &,
210  CopyData &)>>::type
211  &face_worker = std::function<void(const CellIteratorBaseType &,
212  const unsigned int,
213  const unsigned int,
214  const CellIteratorBaseType &,
215  const unsigned int,
216  const unsigned int,
217  ScratchData &,
218  CopyData &)>(),
219 
220  const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
221  const unsigned int chunk_size = 8)
222  {
223  Assert(
224  (!cell_worker) == !(flags & work_on_cells),
225  ExcMessage(
226  "If you specify a cell_worker, you need to set assemble_own_cells or assemble_ghost_cells."));
227 
228  Assert(
229  (flags &
232  ExcMessage(
233  "You can only specify assemble_own_interior_faces_once OR assemble_own_interior_faces_both."));
234 
235  Assert(
238  ExcMessage(
239  "You can only specify assemble_ghost_faces_once OR assemble_ghost_faces_both."));
240 
241  Assert(
242  !(flags & cells_after_faces) ||
244  ExcMessage(
245  "The option cells_after_faces only makes sense if you assemble on cells."));
246 
247  Assert((!face_worker) == !(flags & work_on_faces),
248  ExcMessage(
249  "If you specify a face_worker, assemble_face_* needs to be set."));
250 
251  Assert(
252  (!boundary_worker) == !(flags & assemble_boundary_faces),
253  ExcMessage(
254  "If you specify a boundary_worker, assemble_boundary_faces needs to be set."));
255 
256  auto cell_action = [&](const CellIteratorBaseType &cell,
257  ScratchData & scratch,
258  CopyData & copy) {
259  // First reset the CopyData class to the empty copy_data given by the
260  // user.
261  copy = sample_copy_data;
262 
263  const bool ignore_subdomain =
264  (cell->get_triangulation().locally_owned_subdomain() ==
266 
267  types::subdomain_id current_subdomain_id =
268  (cell->is_level_cell() ? cell->level_subdomain_id() :
269  cell->subdomain_id());
270 
271  const bool own_cell =
272  ignore_subdomain ||
273  (current_subdomain_id ==
274  cell->get_triangulation().locally_owned_subdomain());
275 
276  if ((!ignore_subdomain) &&
277  (current_subdomain_id == numbers::artificial_subdomain_id))
278  return;
279 
280  if (!(flags & (cells_after_faces)) &&
281  (((flags & (assemble_own_cells)) && own_cell) ||
282  ((flags & assemble_ghost_cells) && !own_cell)))
283  cell_worker(cell, scratch, copy);
284 
285  if (flags & (work_on_faces | work_on_boundary))
286  for (unsigned int face_no = 0;
287  face_no < GeometryInfo<CellIteratorBaseType::AccessorType::
288  Container::dimension>::faces_per_cell;
289  ++face_no)
290  {
291  if (cell->at_boundary(face_no) &&
292  !cell->has_periodic_neighbor(face_no))
293  {
294  // only integrate boundary faces of own cells
295  if ((flags & assemble_boundary_faces) && own_cell)
296  boundary_worker(cell, face_no, scratch, copy);
297  }
298  else
299  {
300  // interior face, potentially assemble
302  neighbor = cell->neighbor_or_periodic_neighbor(face_no);
303 
304  types::subdomain_id neighbor_subdomain_id =
306  if (neighbor->is_level_cell())
307  neighbor_subdomain_id = neighbor->level_subdomain_id();
308  // subdomain id is only valid for active cells
309  else if (neighbor->active())
310  neighbor_subdomain_id = neighbor->subdomain_id();
311 
312  const bool own_neighbor =
313  ignore_subdomain ||
314  (neighbor_subdomain_id ==
315  cell->get_triangulation().locally_owned_subdomain());
316 
317  // skip all faces between two ghost cells
318  if (!own_cell && !own_neighbor)
319  continue;
320 
321  // skip if the user doesn't want faces between own cells
322  if (own_cell && own_neighbor &&
325  continue;
326 
327  // skip face to ghost
328  if (own_cell != own_neighbor &&
329  !(flags &
331  continue;
332 
333  // Deal with refinement edges from the refined side. Assuming
334  // one-irregular meshes, this situation should only occur if
335  // both cells are active.
336  const bool periodic_neighbor =
337  cell->has_periodic_neighbor(face_no);
338 
339  if ((!periodic_neighbor &&
340  cell->neighbor_is_coarser(face_no)) ||
341  (periodic_neighbor &&
342  cell->periodic_neighbor_is_coarser(face_no)))
343  {
344  Assert(!cell->has_children(), ExcInternalError());
345  Assert(!neighbor->has_children(), ExcInternalError());
346 
347  // skip if only one processor needs to assemble the face
348  // to a ghost cell and the fine cell is not ours.
349  if (!own_cell && (flags & assemble_ghost_faces_once))
350  continue;
351 
352  const std::pair<unsigned int, unsigned int>
353  neighbor_face_no =
354  periodic_neighbor ?
355  cell->periodic_neighbor_of_coarser_periodic_neighbor(
356  face_no) :
357  cell->neighbor_of_coarser_neighbor(face_no);
358 
359  face_worker(cell,
360  face_no,
362  neighbor,
363  neighbor_face_no.first,
364  neighbor_face_no.second,
365  scratch,
366  copy);
367 
369  {
370  // If own faces are to be assembled from both sides,
371  // call the faceworker again with swapped arguments.
372  // This is because we won't be looking at an adaptively
373  // refined edge coming from the other side.
374  face_worker(neighbor,
375  neighbor_face_no.first,
376  neighbor_face_no.second,
377  cell,
378  face_no,
380  scratch,
381  copy);
382  }
383  }
384  else
385  {
386  // If iterator is active and neighbor is refined, skip
387  // internal face.
388  if (::internal::is_active_iterator(cell) &&
389  neighbor->has_children())
390  continue;
391 
392  // Now neighbor is on same level, double-check this:
393  Assert(cell->level() == neighbor->level(),
394  ExcInternalError());
395 
396  // If we own both cells only do faces from one side (unless
397  // AssembleFlags says otherwise). Here, we rely on cell
398  // comparison that will look at cell->index().
399  if (own_cell && own_neighbor &&
401  (neighbor < cell))
402  continue;
403 
404  // We only look at faces to ghost on the same level once
405  // (only where own_cell=true and own_neighbor=false)
406  if (!own_cell)
407  continue;
408 
409  // now only one processor assembles faces_to_ghost. We let
410  // the processor with the smaller (level-)subdomain id
411  // assemble the face.
412  if (own_cell && !own_neighbor &&
413  (flags & assemble_ghost_faces_once) &&
414  (neighbor_subdomain_id < current_subdomain_id))
415  continue;
416 
417  const unsigned int neighbor_face_no =
418  periodic_neighbor ?
419  cell->periodic_neighbor_face_no(face_no) :
420  cell->neighbor_face_no(face_no);
421  Assert(periodic_neighbor ||
422  neighbor->face(neighbor_face_no) ==
423  cell->face(face_no),
424  ExcInternalError());
425 
426  face_worker(cell,
427  face_no,
429  neighbor,
430  neighbor_face_no,
432  scratch,
433  copy);
434  }
435  }
436  } // faces
437 
438  // Execute the cell_worker if faces are handled before cells
439  if ((flags & cells_after_faces) &&
440  (((flags & assemble_own_cells) && own_cell) ||
441  ((flags & assemble_ghost_cells) && !own_cell)))
442  cell_worker(cell, scratch, copy);
443  };
444 
445  // Submit to workstream
446  WorkStream::run(begin,
447  end,
448  cell_action,
449  copier,
450  sample_scratch_data,
451  sample_copy_data,
452  queue_length,
453  chunk_size);
454  }
455 
525  template <class CellIteratorType,
526  class ScratchData,
527  class CopyData,
528  class CellIteratorBaseType =
530  void
532  IteratorRange<CellIteratorType> iterator_range,
533  const typename identity<std::function<
534  void(const CellIteratorBaseType &, ScratchData &, CopyData &)>>::type
535  &cell_worker,
536  const typename identity<std::function<void(const CopyData &)>>::type
537  &copier,
538 
539  const ScratchData &sample_scratch_data,
540  const CopyData & sample_copy_data,
541 
542  const AssembleFlags flags = assemble_own_cells,
543 
544  const typename identity<std::function<void(const CellIteratorBaseType &,
545  const unsigned int,
546  ScratchData &,
547  CopyData &)>>::type
548  &boundary_worker = std::function<void(const CellIteratorBaseType &,
549  const unsigned int,
550  ScratchData &,
551  CopyData &)>(),
552 
553  const typename identity<std::function<void(const CellIteratorBaseType &,
554  const unsigned int,
555  const unsigned int,
556  const CellIteratorBaseType &,
557  const unsigned int,
558  const unsigned int,
559  ScratchData &,
560  CopyData &)>>::type
561  &face_worker = std::function<void(const CellIteratorBaseType &,
562  const unsigned int,
563  const unsigned int,
564  const CellIteratorBaseType &,
565  const unsigned int,
566  const unsigned int,
567  ScratchData &,
568  CopyData &)>(),
569 
570  const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
571  const unsigned int chunk_size = 8)
572  {
573  // Call the function above
574  mesh_loop<typename IteratorRange<CellIteratorType>::IteratorOverIterators,
575  ScratchData,
576  CopyData,
577  CellIteratorBaseType>(iterator_range.begin(),
578  iterator_range.end(),
579  cell_worker,
580  copier,
581  sample_scratch_data,
582  sample_copy_data,
583  flags,
584  boundary_worker,
585  face_worker,
586  queue_length,
587  chunk_size);
588  }
589 
650  template <class CellIteratorType,
651  class ScratchData,
652  class CopyData,
653  class MainClass>
654  void
655  mesh_loop(const CellIteratorType & begin,
656  const typename identity<CellIteratorType>::type &end,
657  MainClass & main_class,
658  void (MainClass::*cell_worker)(const CellIteratorType &,
659  ScratchData &,
660  CopyData &),
661  void (MainClass::*copier)(const CopyData &),
662  const ScratchData & sample_scratch_data,
663  const CopyData & sample_copy_data,
664  const AssembleFlags flags = assemble_own_cells,
665  void (MainClass::*boundary_worker)(const CellIteratorType &,
666  const unsigned int,
667  ScratchData &,
668  CopyData &) = nullptr,
669  void (MainClass::*face_worker)(const CellIteratorType &,
670  const unsigned int,
671  const unsigned int,
672  const CellIteratorType &,
673  const unsigned int,
674  const unsigned int,
675  ScratchData &,
676  CopyData &) = nullptr,
677  const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
678  const unsigned int chunk_size = 8)
679  {
680  std::function<void(const CellIteratorType &, ScratchData &, CopyData &)>
681  f_cell_worker;
682 
683  std::function<void(
684  const CellIteratorType &, const unsigned int, ScratchData &, CopyData &)>
685  f_boundary_worker;
686 
687  std::function<void(const CellIteratorType &,
688  const unsigned int,
689  const unsigned int,
690  const CellIteratorType &,
691  const unsigned int,
692  const unsigned int,
693  ScratchData &,
694  CopyData &)>
695  f_face_worker;
696 
697  if (cell_worker != nullptr)
698  f_cell_worker = std::bind(cell_worker,
699  std::ref(main_class),
700  std::placeholders::_1,
701  std::placeholders::_2,
702  std::placeholders::_3);
703 
704  if (boundary_worker != nullptr)
705  f_boundary_worker = std::bind(boundary_worker,
706  std::ref(main_class),
707  std::placeholders::_1,
708  std::placeholders::_2,
709  std::placeholders::_3,
710  std::placeholders::_4);
711 
712  if (face_worker != nullptr)
713  f_face_worker = std::bind(face_worker,
714  std::ref(main_class),
715  std::placeholders::_1,
716  std::placeholders::_2,
717  std::placeholders::_3,
718  std::placeholders::_4,
719  std::placeholders::_5,
720  std::placeholders::_6,
721  std::placeholders::_7,
722  std::placeholders::_8);
723 
724  mesh_loop(begin,
725  end,
726  f_cell_worker,
727  std::bind(copier, main_class, std::placeholders::_1),
728  sample_scratch_data,
729  sample_copy_data,
730  flags,
731  f_boundary_worker,
732  f_face_worker,
733  queue_length,
734  chunk_size);
735  }
736 
816  template <class CellIteratorType,
817  class ScratchData,
818  class CopyData,
819  class MainClass,
820  class CellIteratorBaseType =
822  void
824  MainClass & main_class,
825  void (MainClass::*cell_worker)(const CellIteratorBaseType &,
826  ScratchData &,
827  CopyData &),
828  void (MainClass::*copier)(const CopyData &),
829  const ScratchData & sample_scratch_data,
830  const CopyData & sample_copy_data,
831  const AssembleFlags flags = assemble_own_cells,
832  void (MainClass::*boundary_worker)(const CellIteratorBaseType &,
833  const unsigned int,
834  ScratchData &,
835  CopyData &) = nullptr,
836  void (MainClass::*face_worker)(const CellIteratorBaseType &,
837  const unsigned int,
838  const unsigned int,
839  const CellIteratorBaseType &,
840  const unsigned int,
841  const unsigned int,
842  ScratchData &,
843  CopyData &) = nullptr,
844  const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
845  const unsigned int chunk_size = 8)
846  {
847  // Call the function above
848  mesh_loop<typename IteratorRange<CellIteratorType>::IteratorOverIterators,
849  ScratchData,
850  CopyData,
851  MainClass,
852  CellIteratorBaseType>(iterator_range.begin(),
853  iterator_range.end(),
854  main_class,
855  cell_worker,
856  copier,
857  sample_scratch_data,
858  sample_copy_data,
859  flags,
860  boundary_worker,
861  face_worker,
862  queue_length,
863  chunk_size);
864  }
865 } // namespace MeshWorker
866 
867 DEAL_II_NAMESPACE_CLOSE
868 
869 #endif
static const unsigned int invalid_unsigned_int
Definition: types.h:173
IteratorOverIterators begin()
const types::subdomain_id invalid_subdomain_id
Definition: types.h:258
typename CellIteratorBaseType< CellIteratorType >::type type
Definition: mesh_loop.h:94
void mesh_loop(const CellIteratorType &begin, const typename identity< CellIteratorType >::type &end, const typename identity< std::function< void(const CellIteratorBaseType &, ScratchData &, CopyData &)>>::type &cell_worker, const typename identity< std::function< void(const CopyData &)>>::type &copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const AssembleFlags flags=assemble_own_cells, const typename identity< std::function< void(const CellIteratorBaseType &, const unsigned int, ScratchData &, CopyData &)>>::type &boundary_worker=std::function< void(const CellIteratorBaseType &, const unsigned int, ScratchData &, CopyData &)>(), const typename identity< std::function< void(const CellIteratorBaseType &, const unsigned int, const unsigned int, const CellIteratorBaseType &, const unsigned int, const unsigned int, ScratchData &, CopyData &)>>::type &face_worker=std::function< void(const CellIteratorBaseType &, const unsigned int, const unsigned int, const CellIteratorBaseType &, const unsigned int, const unsigned int, ScratchData &, CopyData &)>(), const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: mesh_loop.h:179
typename CellIteratorBaseType< CellIteratorType >::type type
Definition: mesh_loop.h:75
static ::ExceptionBase & ExcMessage(std::string arg1)
unsigned int subdomain_id
Definition: types.h:43
#define Assert(cond, exc)
Definition: exceptions.h:1407
bool is_active_iterator(const DI &)
Definition: loop.h:46
const types::subdomain_id artificial_subdomain_id
Definition: types.h:275
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:1167
static unsigned int n_threads()
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:190
static ::ExceptionBase & ExcInternalError()
IteratorOverIterators end() const