Reference documentation for deal.II version Git d4ec8156ca 2019-11-18 14:03:00 -0700
\(\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  [&cell_worker, &boundary_worker, &face_worker, &lctrl](
475  ITERATOR cell, INFOBOX &info, DoFInfoBox<dim, DOFINFO> &dof_info) {
476  cell_action<INFOBOX, DOFINFO, dim, spacedim, ITERATOR>(cell,
477  dof_info,
478  info,
479  cell_worker,
480  boundary_worker,
481  face_worker,
482  lctrl);
483  },
484  [&assembler](const MeshWorker::DoFInfoBox<dim, DOFINFO> &dinfo) {
485  ::internal::assemble<dim, DOFINFO, ASSEMBLER>(dinfo, &assembler);
486  },
487  info,
488  dof_info);
489  }
490 
491 
499  template <int dim, int spacedim, class ITERATOR, class ASSEMBLER>
500  void
501  integration_loop(ITERATOR begin,
502  typename identity<ITERATOR>::type end,
503  DoFInfo<dim, spacedim> & dof_info,
505  const LocalIntegrator<dim, spacedim> &integrator,
506  ASSEMBLER & assembler,
507  const LoopControl & lctrl = LoopControl())
508  {
509  std::function<void(DoFInfo<dim, spacedim> &,
511  cell_worker;
512  std::function<void(DoFInfo<dim, spacedim> &,
514  boundary_worker;
515  std::function<void(DoFInfo<dim, spacedim> &,
518  IntegrationInfo<dim, spacedim> &)>
519  face_worker;
520  if (integrator.use_cell)
521  cell_worker =
522  [&integrator](DoFInfo<dim, spacedim> & dof_info,
523  IntegrationInfo<dim, spacedim> &integration_info) {
524  integrator.cell(dof_info, integration_info);
525  };
526  if (integrator.use_boundary)
527  boundary_worker =
528  [&integrator](DoFInfo<dim, spacedim> & dof_info,
529  IntegrationInfo<dim, spacedim> &integration_info) {
530  integrator.boundary(dof_info, integration_info);
531  };
532  if (integrator.use_face)
533  face_worker =
534  [&integrator](DoFInfo<dim, spacedim> & dof_info_1,
535  DoFInfo<dim, spacedim> & dof_info_2,
536  IntegrationInfo<dim, spacedim> &integration_info_1,
537  IntegrationInfo<dim, spacedim> &integration_info_2) {
538  integrator.face(dof_info_1,
539  dof_info_2,
540  integration_info_1,
541  integration_info_2);
542  };
543 
544  loop<dim, spacedim>(begin,
545  end,
546  dof_info,
547  box,
548  cell_worker,
549  boundary_worker,
550  face_worker,
551  assembler,
552  lctrl);
553  }
554 
555 } // namespace MeshWorker
556 
557 DEAL_II_NAMESPACE_CLOSE
558 
559 #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
virtual void cell(DoFInfo< dim, spacedim, number > &dinfo, IntegrationInfo< dim, spacedim > &info) const
Definition: mesh_worker.cc:81
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:1411
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:1177
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:501
DOFINFO exterior[GeometryInfo< dim >::faces_per_cell]
Definition: dof_info.h:264
virtual void face(DoFInfo< dim, spacedim, number > &dinfo1, DoFInfo< dim, spacedim, number > &dinfo2, IntegrationInfo< dim, spacedim > &info1, IntegrationInfo< dim, spacedim > &info2) const
Definition: mesh_worker.cc:101
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
virtual void boundary(DoFInfo< dim, spacedim, number > &dinfo, IntegrationInfo< dim, spacedim > &info) const
Definition: mesh_worker.cc:91
static ::ExceptionBase & ExcInternalError()