Reference documentation for deal.II version GIT 85919f3e70 2023-05-28 07:10:01+00:00
\(\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\}}\)
loop.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2006 - 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_loop_h
18 #define dealii_mesh_worker_loop_h
19 
20 #include <deal.II/base/config.h>
21 
24 
26 #include <deal.II/grid/tria.h>
27 
31 
32 #include <functional>
33 
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
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)
95  , cells_first(true)
96  {}
97 
101  bool own_cells;
102 
108 
115  {
127  both
128  };
129 
143 
154 
160  };
161 
162 
163 
189  template <class INFOBOX, class DOFINFO, int dim, int spacedim, class ITERATOR>
190  void
192  ITERATOR cell,
193  DoFInfoBox<dim, DOFINFO> &dof_info,
194  INFOBOX & info,
195  const std::function<void(DOFINFO &, typename INFOBOX::CellInfo &)>
196  &cell_worker,
197  const std::function<void(DOFINFO &, typename INFOBOX::CellInfo &)>
198  & boundary_worker,
199  const std::function<void(DOFINFO &,
200  DOFINFO &,
201  typename INFOBOX::CellInfo &,
202  typename INFOBOX::CellInfo &)> &face_worker,
203  const LoopControl & loop_control)
204  {
205  const bool ignore_subdomain =
206  (cell->get_triangulation().locally_owned_subdomain() ==
208 
209  types::subdomain_id csid = (cell->is_level_cell()) ?
210  cell->level_subdomain_id() :
211  cell->subdomain_id();
212 
213  const bool own_cell =
214  ignore_subdomain ||
215  (csid == cell->get_triangulation().locally_owned_subdomain());
216 
217  dof_info.reset();
218 
219  if ((!ignore_subdomain) && (csid == numbers::artificial_subdomain_id))
220  return;
221 
222  dof_info.cell.reinit(cell);
223  dof_info.cell_valid = true;
224 
225  const bool integrate_cell = (cell_worker != nullptr);
226  const bool integrate_boundary = (boundary_worker != nullptr);
227  const bool integrate_interior_face = (face_worker != nullptr);
228 
229  if (integrate_cell)
230  info.cell.reinit(dof_info.cell);
231  // Execute this, if cells
232  // have to be dealt with
233  // before faces
234  if (integrate_cell && loop_control.cells_first &&
235  ((loop_control.own_cells && own_cell) ||
236  (loop_control.ghost_cells && !own_cell)))
237  cell_worker(dof_info.cell, info.cell);
238 
239  // Call the callback function in
240  // the info box to do
241  // computations between cell and
242  // face action.
243  info.post_cell(dof_info);
244 
245  if (integrate_interior_face || integrate_boundary)
246  for (const unsigned int face_no : cell->face_indices())
247  {
248  typename ITERATOR::AccessorType::Container::face_iterator face =
249  cell->face(face_no);
250  if (cell->at_boundary(face_no) &&
251  !cell->has_periodic_neighbor(face_no))
252  {
253  // only integrate boundary faces of own cells
254  if (integrate_boundary && own_cell)
255  {
256  dof_info.interior_face_available[face_no] = true;
257  dof_info.interior[face_no].reinit(cell, face, face_no);
258  info.boundary.reinit(dof_info.interior[face_no]);
259  boundary_worker(dof_info.interior[face_no], info.boundary);
260  }
261  }
262  else if (integrate_interior_face)
263  {
264  // Interior face
266  cell->neighbor_or_periodic_neighbor(face_no);
267 
269  if (neighbor->is_level_cell())
270  neighbid = neighbor->level_subdomain_id();
271  // subdomain id is only valid for active cells
272  else if (neighbor->is_active())
273  neighbid = neighbor->subdomain_id();
274 
275  const bool own_neighbor =
276  ignore_subdomain ||
277  (neighbid ==
278  cell->get_triangulation().locally_owned_subdomain());
279 
280  // skip all faces between two ghost cells
281  if (!own_cell && !own_neighbor)
282  continue;
283 
284  // skip if the user doesn't want faces between own cells
285  if (own_cell && own_neighbor &&
286  loop_control.own_faces == LoopControl::never)
287  continue;
288 
289  // skip face to ghost
290  if (own_cell != own_neighbor &&
291  loop_control.faces_to_ghost == LoopControl::never)
292  continue;
293 
294  // Deal with refinement edges from the refined side. Assuming
295  // one-irregular meshes, this situation should only occur if both
296  // cells are active.
297  const bool periodic_neighbor =
298  cell->has_periodic_neighbor(face_no);
299 
300  if ((!periodic_neighbor && cell->neighbor_is_coarser(face_no)) ||
301  (periodic_neighbor &&
302  cell->periodic_neighbor_is_coarser(face_no)))
303  {
304  Assert(cell->is_active(), ExcInternalError());
305  Assert(neighbor->is_active(), ExcInternalError());
306 
307  // skip if only one processor needs to assemble the face
308  // to a ghost cell and the fine cell is not ours.
309  if (!own_cell &&
310  loop_control.faces_to_ghost == LoopControl::one)
311  continue;
312 
313  const std::pair<unsigned int, unsigned int> neighbor_face_no =
314  periodic_neighbor ?
315  cell->periodic_neighbor_of_coarser_periodic_neighbor(
316  face_no) :
317  cell->neighbor_of_coarser_neighbor(face_no);
318  const typename ITERATOR::AccessorType::Container::
319  face_iterator nface =
320  neighbor->face(neighbor_face_no.first);
321 
322  dof_info.interior_face_available[face_no] = true;
323  dof_info.exterior_face_available[face_no] = true;
324  dof_info.interior[face_no].reinit(cell, face, face_no);
325  info.face.reinit(dof_info.interior[face_no]);
326  dof_info.exterior[face_no].reinit(neighbor,
327  nface,
328  neighbor_face_no.first,
329  neighbor_face_no.second);
330  info.subface.reinit(dof_info.exterior[face_no]);
331 
332  face_worker(dof_info.interior[face_no],
333  dof_info.exterior[face_no],
334  info.face,
335  info.subface);
336  }
337  else
338  {
339  // If iterator is active and neighbor is refined, skip
340  // internal face.
341  if (::internal::is_active_iterator(cell) &&
342  neighbor->has_children())
343  {
344  Assert(
345  loop_control.own_faces != LoopControl::both,
346  ExcMessage(
347  "Assembling from both sides for own_faces is not "
348  "supported with hanging nodes!"));
349  continue;
350  }
351 
352  // Now neighbor is on same level, double-check this:
353  Assert(cell->level() == neighbor->level(),
354  ExcInternalError());
355 
356  // If we own both cells only do faces from one side (unless
357  // LoopControl says otherwise). Here, we rely on cell
358  // comparison that will look at cell->index().
359  if (own_cell && own_neighbor &&
360  loop_control.own_faces == LoopControl::one &&
361  (neighbor < cell))
362  continue;
363 
364  // independent of loop_control.faces_to_ghost,
365  // we only look at faces to ghost on the same level once
366  // (only where own_cell=true and own_neighbor=false)
367  if (!own_cell)
368  continue;
369 
370  // now only one processor assembles faces_to_ghost. We let the
371  // processor with the smaller (level-)subdomain id assemble
372  // the face.
373  if (own_cell && !own_neighbor &&
374  loop_control.faces_to_ghost == LoopControl::one &&
375  (neighbid < csid))
376  continue;
377 
378  const unsigned int neighbor_face_no =
379  periodic_neighbor ?
380  cell->periodic_neighbor_face_no(face_no) :
381  cell->neighbor_face_no(face_no);
382  Assert(periodic_neighbor ||
383  neighbor->face(neighbor_face_no) == face,
384  ExcInternalError());
385  // Regular interior face
386  dof_info.interior_face_available[face_no] = true;
387  dof_info.exterior_face_available[face_no] = true;
388  dof_info.interior[face_no].reinit(cell, face, face_no);
389  info.face.reinit(dof_info.interior[face_no]);
390  dof_info.exterior[face_no].reinit(neighbor,
391  neighbor->face(
392  neighbor_face_no),
393  neighbor_face_no);
394  info.neighbor.reinit(dof_info.exterior[face_no]);
395 
396  face_worker(dof_info.interior[face_no],
397  dof_info.exterior[face_no],
398  info.face,
399  info.neighbor);
400  }
401  }
402  } // faces
403  // Call the callback function in
404  // the info box to do
405  // computations between face and
406  // cell action.
407  info.post_faces(dof_info);
408 
409  // Execute this, if faces
410  // have to be handled first
411  if (integrate_cell && !loop_control.cells_first &&
412  ((loop_control.own_cells && own_cell) ||
413  (loop_control.ghost_cells && !own_cell)))
414  cell_worker(dof_info.cell, info.cell);
415  }
416 
417 
432  template <int dim,
433  int spacedim,
434  class DOFINFO,
435  class INFOBOX,
436  class ASSEMBLER,
437  class ITERATOR>
438  void
439  loop(ITERATOR begin,
441  DOFINFO & dinfo,
442  INFOBOX & info,
443  const std::function<void(DOFINFO &, typename INFOBOX::CellInfo &)>
444  &cell_worker,
445  const std::function<void(DOFINFO &, typename INFOBOX::CellInfo &)>
446  & boundary_worker,
447  const std::function<void(DOFINFO &,
448  DOFINFO &,
449  typename INFOBOX::CellInfo &,
450  typename INFOBOX::CellInfo &)> &face_worker,
451  ASSEMBLER & assembler,
452  const LoopControl &lctrl = LoopControl())
453  {
454  DoFInfoBox<dim, DOFINFO> dof_info(dinfo);
455 
456  assembler.initialize_info(dof_info.cell, false);
457  for (unsigned int i : GeometryInfo<dim>::face_indices())
458  {
459  assembler.initialize_info(dof_info.interior[i], true);
460  assembler.initialize_info(dof_info.exterior[i], true);
461  }
462 
463  // Loop over all cells
465  begin,
466  end,
467  [&cell_worker, &boundary_worker, &face_worker, &lctrl](
468  ITERATOR cell, INFOBOX &info, DoFInfoBox<dim, DOFINFO> &dof_info) {
469  cell_action<INFOBOX, DOFINFO, dim, spacedim, ITERATOR>(cell,
470  dof_info,
471  info,
472  cell_worker,
473  boundary_worker,
474  face_worker,
475  lctrl);
476  },
477  [&assembler](const MeshWorker::DoFInfoBox<dim, DOFINFO> &dinfo) {
478  ::internal::assemble<dim, DOFINFO, ASSEMBLER>(dinfo, &assembler);
479  },
480  info,
481  dof_info);
482  }
483 
484 
491  template <int dim, int spacedim, class ITERATOR, class ASSEMBLER>
492  void
495  DoFInfo<dim, spacedim> & dof_info,
497  const LocalIntegrator<dim, spacedim> &integrator,
498  ASSEMBLER & assembler,
499  const LoopControl & lctrl = LoopControl())
500  {
501  std::function<void(DoFInfo<dim, spacedim> &,
503  cell_worker;
504  std::function<void(DoFInfo<dim, spacedim> &,
506  boundary_worker;
507  std::function<void(DoFInfo<dim, spacedim> &,
511  face_worker;
512  if (integrator.use_cell)
513  cell_worker =
514  [&integrator](DoFInfo<dim, spacedim> & dof_info,
515  IntegrationInfo<dim, spacedim> &integration_info) {
516  integrator.cell(dof_info, integration_info);
517  };
518  if (integrator.use_boundary)
519  boundary_worker =
520  [&integrator](DoFInfo<dim, spacedim> & dof_info,
521  IntegrationInfo<dim, spacedim> &integration_info) {
522  integrator.boundary(dof_info, integration_info);
523  };
524  if (integrator.use_face)
525  face_worker =
526  [&integrator](DoFInfo<dim, spacedim> & dof_info_1,
527  DoFInfo<dim, spacedim> & dof_info_2,
528  IntegrationInfo<dim, spacedim> &integration_info_1,
529  IntegrationInfo<dim, spacedim> &integration_info_2) {
530  integrator.face(dof_info_1,
531  dof_info_2,
532  integration_info_1,
533  integration_info_2);
534  };
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 
550 
551 #endif
DOFINFO interior[GeometryInfo< dim >::faces_per_cell]
Definition: dof_info.h:264
void assemble(ASSEMBLER &ass) const
Definition: dof_info.h:493
bool exterior_face_available[GeometryInfo< dim >::faces_per_cell]
Definition: dof_info.h:280
bool interior_face_available[GeometryInfo< dim >::faces_per_cell]
Definition: dof_info.h:274
DOFINFO exterior[GeometryInfo< dim >::faces_per_cell]
Definition: dof_info.h:268
virtual void cell(DoFInfo< dim, spacedim, number > &dinfo, IntegrationInfo< dim, spacedim > &info) const
Definition: mesh_worker.cc:81
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
virtual void boundary(DoFInfo< dim, spacedim, number > &dinfo, IntegrationInfo< dim, spacedim > &info) const
Definition: mesh_worker.cc:91
FaceOption own_faces
Definition: loop.h:153
FaceOption faces_to_ghost
Definition: loop.h:142
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcMessage(std::string arg1)
void integration_loop(ITERATOR begin, std_cxx20::type_identity_t< ITERATOR > end, DoFInfo< dim, spacedim > &dof_info, IntegrationInfoBox< dim, spacedim > &box, const LocalIntegrator< dim, spacedim > &integrator, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
Definition: loop.h:493
void loop(ITERATOR begin, std_cxx20::type_identity_t< ITERATOR > 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:439
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 const char A
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
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:1270
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition: loop.h:71
bool is_active_iterator(const DI &)
Definition: loop.h:49
const types::subdomain_id artificial_subdomain_id
Definition: types.h:315
const types::subdomain_id invalid_subdomain_id
Definition: types.h:298
typename type_identity< T >::type type_identity_t
Definition: type_traits.h:96
unsigned int subdomain_id
Definition: types.h:44