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