Reference documentation for deal.II version Git a5ed68a04a 2019-09-22 06:50:58 -0600
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
loop.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2006 - 2018 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_loop_h
18 #define dealii_mesh_worker_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/dof_info.h>
29 #include <deal.II/meshworker/integration_info.h>
30 #include <deal.II/meshworker/local_integrator.h>
31 
32 #include <functional>
33 
34 DEAL_II_NAMESPACE_OPEN
35 
36 // Forward declaration
37 #ifndef DOXYGEN
38 template <typename>
39 class TriaActiveIterator;
40 #endif
41 
42 namespace internal
43 {
47  template <class DI>
48  inline bool
49  is_active_iterator(const DI &)
50  {
51  return false;
52  }
53 
54  template <class ACCESSOR>
55  inline bool
57  {
58  return true;
59  }
60 
61  template <class ACCESSOR>
62  inline bool
64  const ::FilteredIterator<TriaActiveIterator<ACCESSOR>> &)
65  {
66  return true;
67  }
68 
69  template <int dim, class DOFINFO, class A>
70  void
71  assemble(const MeshWorker::DoFInfoBox<dim, DOFINFO> &dinfo, A *assembler)
72  {
73  dinfo.assemble(*assembler);
74  }
75 } // namespace internal
76 
77 
78 
79 namespace MeshWorker
80 {
85  {
86  public:
91  : own_cells(true)
92  , ghost_cells(false)
93  , faces_to_ghost(LoopControl::one)
94  , own_faces(LoopControl::one)
95  , cells_first(true)
96  {}
97 
101  bool own_cells;
102 
108 
115  {
127  both
128  };
129 
143 
154 
160  };
161 
162 
163 
191  template <class INFOBOX, class DOFINFO, int dim, int spacedim, class ITERATOR>
192  void
194  ITERATOR cell,
195  DoFInfoBox<dim, DOFINFO> &dof_info,
196  INFOBOX & info,
197  const std::function<void(DOFINFO &, typename INFOBOX::CellInfo &)>
198  &cell_worker,
199  const std::function<void(DOFINFO &, typename INFOBOX::CellInfo &)>
200  & boundary_worker,
201  const std::function<void(DOFINFO &,
202  DOFINFO &,
203  typename INFOBOX::CellInfo &,
204  typename INFOBOX::CellInfo &)> &face_worker,
205  const LoopControl & loop_control)
206  {
207  const bool ignore_subdomain =
208  (cell->get_triangulation().locally_owned_subdomain() ==
210 
211  types::subdomain_id csid = (cell->is_level_cell()) ?
212  cell->level_subdomain_id() :
213  cell->subdomain_id();
214 
215  const bool own_cell =
216  ignore_subdomain ||
217  (csid == cell->get_triangulation().locally_owned_subdomain());
218 
219  dof_info.reset();
220 
221  if ((!ignore_subdomain) && (csid == numbers::artificial_subdomain_id))
222  return;
223 
224  dof_info.cell.reinit(cell);
225  dof_info.cell_valid = true;
226 
227  const bool integrate_cell = (cell_worker != nullptr);
228  const bool integrate_boundary = (boundary_worker != nullptr);
229  const bool integrate_interior_face = (face_worker != nullptr);
230 
231  if (integrate_cell)
232  info.cell.reinit(dof_info.cell);
233  // Execute this, if cells
234  // have to be dealt with
235  // before faces
236  if (integrate_cell && loop_control.cells_first &&
237  ((loop_control.own_cells && own_cell) ||
238  (loop_control.ghost_cells && !own_cell)))
239  cell_worker(dof_info.cell, info.cell);
240 
241  // Call the callback function in
242  // the info box to do
243  // computations between cell and
244  // face action.
245  info.post_cell(dof_info);
246 
247  if (integrate_interior_face || integrate_boundary)
248  for (unsigned int face_no = 0;
249  face_no <
250  GeometryInfo<
251  ITERATOR::AccessorType::Container::dimension>::faces_per_cell;
252  ++face_no)
253  {
254  typename ITERATOR::AccessorType::Container::face_iterator face =
255  cell->face(face_no);
256  if (cell->at_boundary(face_no) &&
257  !cell->has_periodic_neighbor(face_no))
258  {
259  // only integrate boundary faces of own cells
260  if (integrate_boundary && own_cell)
261  {
262  dof_info.interior_face_available[face_no] = true;
263  dof_info.interior[face_no].reinit(cell, face, face_no);
264  info.boundary.reinit(dof_info.interior[face_no]);
265  boundary_worker(dof_info.interior[face_no], info.boundary);
266  }
267  }
268  else if (integrate_interior_face)
269  {
270  // Interior face
272  cell->neighbor_or_periodic_neighbor(face_no);
273 
275  if (neighbor->is_level_cell())
276  neighbid = neighbor->level_subdomain_id();
277  // subdomain id is only valid for active cells
278  else if (neighbor->active())
279  neighbid = neighbor->subdomain_id();
280 
281  const bool own_neighbor =
282  ignore_subdomain ||
283  (neighbid ==
284  cell->get_triangulation().locally_owned_subdomain());
285 
286  // skip all faces between two ghost cells
287  if (!own_cell && !own_neighbor)
288  continue;
289 
290  // skip if the user doesn't want faces between own cells
291  if (own_cell && own_neighbor &&
292  loop_control.own_faces == LoopControl::never)
293  continue;
294 
295  // skip face to ghost
296  if (own_cell != own_neighbor &&
297  loop_control.faces_to_ghost == LoopControl::never)
298  continue;
299 
300  // Deal with refinement edges from the refined side. Assuming
301  // one-irregular meshes, this situation should only occur if both
302  // cells are active.
303  const bool periodic_neighbor =
304  cell->has_periodic_neighbor(face_no);
305 
306  if ((!periodic_neighbor && cell->neighbor_is_coarser(face_no)) ||
307  (periodic_neighbor &&
308  cell->periodic_neighbor_is_coarser(face_no)))
309  {
310  Assert(!cell->has_children(), ExcInternalError());
311  Assert(!neighbor->has_children(), ExcInternalError());
312 
313  // skip if only one processor needs to assemble the face
314  // to a ghost cell and the fine cell is not ours.
315  if (!own_cell &&
316  loop_control.faces_to_ghost == LoopControl::one)
317  continue;
318 
319  const std::pair<unsigned int, unsigned int> neighbor_face_no =
320  periodic_neighbor ?
321  cell->periodic_neighbor_of_coarser_periodic_neighbor(
322  face_no) :
323  cell->neighbor_of_coarser_neighbor(face_no);
324  const typename ITERATOR::AccessorType::Container::
325  face_iterator nface =
326  neighbor->face(neighbor_face_no.first);
327 
328  dof_info.interior_face_available[face_no] = true;
329  dof_info.exterior_face_available[face_no] = true;
330  dof_info.interior[face_no].reinit(cell, face, face_no);
331  info.face.reinit(dof_info.interior[face_no]);
332  dof_info.exterior[face_no].reinit(neighbor,
333  nface,
334  neighbor_face_no.first,
335  neighbor_face_no.second);
336  info.subface.reinit(dof_info.exterior[face_no]);
337 
338  face_worker(dof_info.interior[face_no],
339  dof_info.exterior[face_no],
340  info.face,
341  info.subface);
342  }
343  else
344  {
345  // If iterator is active and neighbor is refined, skip
346  // internal face.
347  if (::internal::is_active_iterator(cell) &&
348  neighbor->has_children())
349  {
350  Assert(
351  loop_control.own_faces != LoopControl::both,
352  ExcMessage(
353  "Assembling from both sides for own_faces is not "
354  "supported with hanging nodes!"));
355  continue;
356  }
357 
358  // Now neighbor is on same level, double-check this:
359  Assert(cell->level() == neighbor->level(),
360  ExcInternalError());
361 
362  // If we own both cells only do faces from one side (unless
363  // LoopControl says otherwise). Here, we rely on cell
364  // comparison that will look at cell->index().
365  if (own_cell && own_neighbor &&
366  loop_control.own_faces == LoopControl::one &&
367  (neighbor < cell))
368  continue;
369 
370  // independent of loop_control.faces_to_ghost,
371  // we only look at faces to ghost on the same level once
372  // (only where own_cell=true and own_neighbor=false)
373  if (!own_cell)
374  continue;
375 
376  // now only one processor assembles faces_to_ghost. We let the
377  // processor with the smaller (level-)subdomain id assemble
378  // the face.
379  if (own_cell && !own_neighbor &&
380  loop_control.faces_to_ghost == LoopControl::one &&
381  (neighbid < csid))
382  continue;
383 
384  const unsigned int neighbor_face_no =
385  periodic_neighbor ?
386  cell->periodic_neighbor_face_no(face_no) :
387  cell->neighbor_face_no(face_no);
388  Assert(periodic_neighbor ||
389  neighbor->face(neighbor_face_no) == face,
390  ExcInternalError());
391  // Regular interior face
392  dof_info.interior_face_available[face_no] = true;
393  dof_info.exterior_face_available[face_no] = true;
394  dof_info.interior[face_no].reinit(cell, face, face_no);
395  info.face.reinit(dof_info.interior[face_no]);
396  dof_info.exterior[face_no].reinit(neighbor,
397  neighbor->face(
398  neighbor_face_no),
399  neighbor_face_no);
400  info.neighbor.reinit(dof_info.exterior[face_no]);
401 
402  face_worker(dof_info.interior[face_no],
403  dof_info.exterior[face_no],
404  info.face,
405  info.neighbor);
406  }
407  }
408  } // faces
409  // Call the callback function in
410  // the info box to do
411  // computations between face and
412  // cell action.
413  info.post_faces(dof_info);
414 
415  // Execute this, if faces
416  // have to be handled first
417  if (integrate_cell && !loop_control.cells_first &&
418  ((loop_control.own_cells && own_cell) ||
419  (loop_control.ghost_cells && !own_cell)))
420  cell_worker(dof_info.cell, info.cell);
421  }
422 
423 
439  template <int dim,
440  int spacedim,
441  class DOFINFO,
442  class INFOBOX,
443  class ASSEMBLER,
444  class ITERATOR>
445  void
446  loop(ITERATOR begin,
447  typename identity<ITERATOR>::type end,
448  DOFINFO & dinfo,
449  INFOBOX & info,
450  const std::function<void(DOFINFO &, typename INFOBOX::CellInfo &)>
451  &cell_worker,
452  const std::function<void(DOFINFO &, typename INFOBOX::CellInfo &)>
453  & boundary_worker,
454  const std::function<void(DOFINFO &,
455  DOFINFO &,
456  typename INFOBOX::CellInfo &,
457  typename INFOBOX::CellInfo &)> &face_worker,
458  ASSEMBLER & assembler,
459  const LoopControl &lctrl = LoopControl())
460  {
461  DoFInfoBox<dim, DOFINFO> dof_info(dinfo);
462 
463  assembler.initialize_info(dof_info.cell, false);
464  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
465  {
466  assembler.initialize_info(dof_info.interior[i], true);
467  assembler.initialize_info(dof_info.exterior[i], true);
468  }
469 
470  // Loop over all cells
472  begin,
473  end,
474  std::bind(&cell_action<INFOBOX, DOFINFO, dim, spacedim, ITERATOR>,
475  std::placeholders::_1,
476  std::placeholders::_3,
477  std::placeholders::_2,
478  cell_worker,
479  boundary_worker,
480  face_worker,
481  lctrl),
482  std::bind(&::internal::assemble<dim, DOFINFO, ASSEMBLER>,
483  std::placeholders::_1,
484  &assembler),
485  info,
486  dof_info);
487  }
488 
489 
497  template <int dim, int spacedim, class ITERATOR, class ASSEMBLER>
498  void
499  integration_loop(ITERATOR begin,
500  typename identity<ITERATOR>::type end,
501  DoFInfo<dim, spacedim> & dof_info,
503  const LocalIntegrator<dim, spacedim> &integrator,
504  ASSEMBLER & assembler,
505  const LoopControl & lctrl = LoopControl())
506  {
507  std::function<void(DoFInfo<dim, spacedim> &,
509  cell_worker;
510  std::function<void(DoFInfo<dim, spacedim> &,
512  boundary_worker;
513  std::function<void(DoFInfo<dim, spacedim> &,
516  IntegrationInfo<dim, spacedim> &)>
517  face_worker;
518  if (integrator.use_cell)
519  cell_worker = std::bind(&LocalIntegrator<dim, spacedim>::cell,
520  &integrator,
521  std::placeholders::_1,
522  std::placeholders::_2);
523  if (integrator.use_boundary)
524  boundary_worker = std::bind(&LocalIntegrator<dim, spacedim>::boundary,
525  &integrator,
526  std::placeholders::_1,
527  std::placeholders::_2);
528  if (integrator.use_face)
529  face_worker = std::bind(&LocalIntegrator<dim, spacedim>::face,
530  &integrator,
531  std::placeholders::_1,
532  std::placeholders::_2,
533  std::placeholders::_3,
534  std::placeholders::_4);
535 
536  loop<dim, spacedim>(begin,
537  end,
538  dof_info,
539  box,
540  cell_worker,
541  boundary_worker,
542  face_worker,
543  assembler,
544  lctrl);
545  }
546 
547 } // namespace MeshWorker
548 
549 DEAL_II_NAMESPACE_CLOSE
550 
551 #endif
void loop(ITERATOR begin, typename identity< ITERATOR >::type end, DOFINFO &dinfo, 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, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
Definition: loop.h:446
const types::subdomain_id invalid_subdomain_id
Definition: types.h:279
FaceOption own_faces
Definition: loop.h:153
bool exterior_face_available[GeometryInfo< dim >::faces_per_cell]
Definition: dof_info.h:276
static ::ExceptionBase & ExcMessage(std::string arg1)
unsigned int subdomain_id
Definition: types.h:43
bool interior_face_available[GeometryInfo< dim >::faces_per_cell]
Definition: dof_info.h:270
#define Assert(cond, exc)
Definition: exceptions.h:1407
bool is_active_iterator(const DI &)
Definition: loop.h:49
void assemble(ASSEMBLER &ass) const
Definition: dof_info.h:472
const types::subdomain_id artificial_subdomain_id
Definition: types.h:296
DOFINFO interior[GeometryInfo< dim >::faces_per_cell]
Definition: dof_info.h:260
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
void integration_loop(ITERATOR begin, typename identity< ITERATOR >::type end, DoFInfo< dim, spacedim > &dof_info, IntegrationInfoBox< dim, spacedim > &box, const LocalIntegrator< dim, spacedim > &integrator, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
Definition: loop.h:499
DOFINFO exterior[GeometryInfo< dim >::faces_per_cell]
Definition: dof_info.h:264
FaceOption faces_to_ghost
Definition: loop.h:142
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:193
static ::ExceptionBase & ExcInternalError()