deal.II version GIT relicensing-1989-gd7a2c90e4e 2024-10-14 01:50:00+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\}}\)
Loading...
Searching...
No Matches
loop.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2009 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
16#ifndef dealii_mesh_worker_loop_h
17#define dealii_mesh_worker_loop_h
18
19#include <deal.II/base/config.h>
20
23
25#include <deal.II/grid/tria.h>
26
30
31#include <functional>
32
34
35// Forward declaration
36#ifndef DOXYGEN
37template <typename>
39#endif
40
41namespace internal
42{
46 template <class DI>
47 inline bool
49 {
50 return false;
51 }
52
53 template <typename AccessorType>
54 inline bool
56 {
57 return true;
58 }
59
60 template <typename AccessorType>
61 inline bool
63 const ::FilteredIterator<TriaActiveIterator<AccessorType>> &)
64 {
65 return true;
66 }
67
68 template <int dim, class DOFINFO, class A>
69 void
71 {
72 dinfo.assemble(*assembler);
73 }
74} // namespace internal
75
76
77
200namespace MeshWorker
201{
206 {
207 public:
212 : own_cells(true)
213 , ghost_cells(false)
216 , cells_first(true)
217 {}
218
223
229
250
264
275
281 };
282
283
284
310 template <class INFOBOX,
311 class DOFINFO,
312 int dim,
313 int spacedim,
314 typename IteratorType>
315 void
317 IteratorType cell,
318 DoFInfoBox<dim, DOFINFO> &dof_info,
319 INFOBOX &info,
320 const std::function<void(DOFINFO &, typename INFOBOX::CellInfo &)>
321 &cell_worker,
322 const std::function<void(DOFINFO &, typename INFOBOX::CellInfo &)>
323 &boundary_worker,
324 const std::function<void(DOFINFO &,
325 DOFINFO &,
326 typename INFOBOX::CellInfo &,
327 typename INFOBOX::CellInfo &)> &face_worker,
328 const LoopControl &loop_control)
329 {
330 const bool ignore_subdomain =
331 (cell->get_triangulation().locally_owned_subdomain() ==
333
334 types::subdomain_id csid = (cell->is_level_cell()) ?
335 cell->level_subdomain_id() :
336 cell->subdomain_id();
337
338 const bool own_cell =
339 ignore_subdomain ||
340 (csid == cell->get_triangulation().locally_owned_subdomain());
341
342 dof_info.reset();
343
344 if ((!ignore_subdomain) && (csid == numbers::artificial_subdomain_id))
345 return;
346
347 dof_info.cell.reinit(cell);
348 dof_info.cell_valid = true;
349
350 const bool integrate_cell = (cell_worker != nullptr);
351 const bool integrate_boundary = (boundary_worker != nullptr);
352 const bool integrate_interior_face = (face_worker != nullptr);
353
354 if (integrate_cell)
355 info.cell.reinit(dof_info.cell);
356 // Execute this, if cells
357 // have to be dealt with
358 // before faces
359 if (integrate_cell && loop_control.cells_first &&
360 ((loop_control.own_cells && own_cell) ||
361 (loop_control.ghost_cells && !own_cell)))
362 cell_worker(dof_info.cell, info.cell);
363
364 // Call the callback function in
365 // the info box to do
366 // computations between cell and
367 // face action.
368 info.post_cell(dof_info);
369
370 if (integrate_interior_face || integrate_boundary)
371 for (const unsigned int face_no : cell->face_indices())
372 {
373 typename IteratorType::AccessorType::Container::face_iterator face =
374 cell->face(face_no);
375 if (cell->at_boundary(face_no) &&
376 !cell->has_periodic_neighbor(face_no))
377 {
378 // only integrate boundary faces of own cells
379 if (integrate_boundary && own_cell)
380 {
381 dof_info.interior_face_available[face_no] = true;
382 dof_info.interior[face_no].reinit(cell, face, face_no);
383 info.boundary.reinit(dof_info.interior[face_no]);
384 boundary_worker(dof_info.interior[face_no], info.boundary);
385 }
386 }
387 else if (integrate_interior_face)
388 {
389 // Interior face
391 cell->neighbor_or_periodic_neighbor(face_no);
392
394 if (neighbor->is_level_cell())
395 neighbid = neighbor->level_subdomain_id();
396 // subdomain id is only valid for active cells
397 else if (neighbor->is_active())
398 neighbid = neighbor->subdomain_id();
399
400 const bool own_neighbor =
401 ignore_subdomain ||
402 (neighbid ==
403 cell->get_triangulation().locally_owned_subdomain());
404
405 // skip all faces between two ghost cells
406 if (!own_cell && !own_neighbor)
407 continue;
408
409 // skip if the user doesn't want faces between own cells
410 if (own_cell && own_neighbor &&
411 loop_control.own_faces == LoopControl::never)
412 continue;
413
414 // skip face to ghost
415 if (own_cell != own_neighbor &&
416 loop_control.faces_to_ghost == LoopControl::never)
417 continue;
418
419 // Deal with refinement edges from the refined side. Assuming
420 // one-irregular meshes, this situation should only occur if both
421 // cells are active.
422 const bool periodic_neighbor =
423 cell->has_periodic_neighbor(face_no);
424
425 if ((!periodic_neighbor && cell->neighbor_is_coarser(face_no)) ||
426 (periodic_neighbor &&
427 cell->periodic_neighbor_is_coarser(face_no)))
428 {
429 Assert(cell->is_active(), ExcInternalError());
430 Assert(neighbor->is_active(), ExcInternalError());
431
432 // skip if only one processor needs to assemble the face
433 // to a ghost cell and the fine cell is not ours.
434 if (!own_cell &&
435 loop_control.faces_to_ghost == LoopControl::one)
436 continue;
437
438 const std::pair<unsigned int, unsigned int> neighbor_face_no =
439 periodic_neighbor ?
440 cell->periodic_neighbor_of_coarser_periodic_neighbor(
441 face_no) :
442 cell->neighbor_of_coarser_neighbor(face_no);
443 const typename IteratorType::AccessorType::Container::
444 face_iterator nface =
445 neighbor->face(neighbor_face_no.first);
446
447 dof_info.interior_face_available[face_no] = true;
448 dof_info.exterior_face_available[face_no] = true;
449 dof_info.interior[face_no].reinit(cell, face, face_no);
450 info.face.reinit(dof_info.interior[face_no]);
451 dof_info.exterior[face_no].reinit(neighbor,
452 nface,
453 neighbor_face_no.first,
454 neighbor_face_no.second);
455 info.subface.reinit(dof_info.exterior[face_no]);
456
457 face_worker(dof_info.interior[face_no],
458 dof_info.exterior[face_no],
459 info.face,
460 info.subface);
461 }
462 else
463 {
464 // If iterator is active and neighbor is refined, skip
465 // internal face.
467 neighbor->has_children())
468 {
469 Assert(
470 loop_control.own_faces != LoopControl::both,
472 "Assembling from both sides for own_faces is not "
473 "supported with hanging nodes!"));
474 continue;
475 }
476
477 // Now neighbor is on same level, double-check this:
478 Assert(cell->level() == neighbor->level(),
480
481 // If we own both cells only do faces from one side (unless
482 // LoopControl says otherwise). Here, we rely on cell
483 // comparison that will look at cell->index().
484 if (own_cell && own_neighbor &&
485 loop_control.own_faces == LoopControl::one &&
486 (neighbor < cell))
487 continue;
488
489 // independent of loop_control.faces_to_ghost,
490 // we only look at faces to ghost on the same level once
491 // (only where own_cell=true and own_neighbor=false)
492 if (!own_cell)
493 continue;
494
495 // now only one processor assembles faces_to_ghost. We let the
496 // processor with the smaller (level-)subdomain id assemble
497 // the face.
498 if (own_cell && !own_neighbor &&
499 loop_control.faces_to_ghost == LoopControl::one &&
500 (neighbid < csid))
501 continue;
502
503 const unsigned int neighbor_face_no =
504 periodic_neighbor ?
505 cell->periodic_neighbor_face_no(face_no) :
506 cell->neighbor_face_no(face_no);
507 Assert(periodic_neighbor ||
508 neighbor->face(neighbor_face_no) == face,
510 // Regular interior face
511 dof_info.interior_face_available[face_no] = true;
512 dof_info.exterior_face_available[face_no] = true;
513 dof_info.interior[face_no].reinit(cell, face, face_no);
514 info.face.reinit(dof_info.interior[face_no]);
515 dof_info.exterior[face_no].reinit(neighbor,
516 neighbor->face(
517 neighbor_face_no),
518 neighbor_face_no);
519 info.neighbor.reinit(dof_info.exterior[face_no]);
520
521 face_worker(dof_info.interior[face_no],
522 dof_info.exterior[face_no],
523 info.face,
524 info.neighbor);
525 }
526 }
527 } // faces
528 // Call the callback function in
529 // the info box to do
530 // computations between face and
531 // cell action.
532 info.post_faces(dof_info);
533
534 // Execute this, if faces
535 // have to be handled first
536 if (integrate_cell && !loop_control.cells_first &&
537 ((loop_control.own_cells && own_cell) ||
538 (loop_control.ghost_cells && !own_cell)))
539 cell_worker(dof_info.cell, info.cell);
540 }
541
542
557 template <int dim,
558 int spacedim,
559 class DOFINFO,
560 class INFOBOX,
561 typename AssemblerType,
562 typename IteratorType>
563 void
564 loop(IteratorType begin,
566 DOFINFO &dinfo,
567 INFOBOX &info,
568 const std::function<void(std_cxx20::type_identity_t<DOFINFO> &,
569 typename INFOBOX::CellInfo &)> &cell_worker,
570 const std::function<void(std_cxx20::type_identity_t<DOFINFO> &,
571 typename INFOBOX::CellInfo &)> &boundary_worker,
572 const std::function<void(std_cxx20::type_identity_t<DOFINFO> &,
574 typename INFOBOX::CellInfo &,
575 typename INFOBOX::CellInfo &)> &face_worker,
576 AssemblerType &assembler,
577 const LoopControl &lctrl = LoopControl())
578 {
579 DoFInfoBox<dim, DOFINFO> dof_info(dinfo);
580
581 assembler.initialize_info(dof_info.cell, false);
582 for (const unsigned int i : GeometryInfo<dim>::face_indices())
583 {
584 assembler.initialize_info(dof_info.interior[i], true);
585 assembler.initialize_info(dof_info.exterior[i], true);
586 }
587
588 // Loop over all cells
590 begin,
591 end,
592 [&cell_worker, &boundary_worker, &face_worker, &lctrl](
593 IteratorType cell, INFOBOX &info, DoFInfoBox<dim, DOFINFO> &dof_info) {
594 cell_action<INFOBOX, DOFINFO, dim, spacedim, IteratorType>(
595 cell,
596 dof_info,
597 info,
598 cell_worker,
599 boundary_worker,
600 face_worker,
601 lctrl);
602 },
603 [&assembler](const MeshWorker::DoFInfoBox<dim, DOFINFO> &dinfo) {
604 ::internal::assemble<dim, DOFINFO, AssemblerType>(dinfo,
605 &assembler);
606 },
607 info,
608 dof_info);
609 }
610
611
622 template <int dim,
623 int spacedim,
624 typename IteratorType,
625 typename AssemblerType>
627 integration_loop(IteratorType begin,
629 DoFInfo<dim, spacedim> &dof_info,
631 const LocalIntegrator<dim, spacedim> &integrator,
632 AssemblerType &assembler,
633 const LoopControl &lctrl = LoopControl())
634 {
635 std::function<void(DoFInfo<dim, spacedim> &,
637 cell_worker;
638 std::function<void(DoFInfo<dim, spacedim> &,
640 boundary_worker;
641 std::function<void(DoFInfo<dim, spacedim> &,
645 face_worker;
646 if (integrator.use_cell)
647 cell_worker =
648 [&integrator](DoFInfo<dim, spacedim> &dof_info,
649 IntegrationInfo<dim, spacedim> &integration_info) {
650 integrator.cell(dof_info, integration_info);
651 };
652 if (integrator.use_boundary)
653 boundary_worker =
654 [&integrator](DoFInfo<dim, spacedim> &dof_info,
655 IntegrationInfo<dim, spacedim> &integration_info) {
656 integrator.boundary(dof_info, integration_info);
657 };
658 if (integrator.use_face)
659 face_worker =
660 [&integrator](DoFInfo<dim, spacedim> &dof_info_1,
661 DoFInfo<dim, spacedim> &dof_info_2,
662 IntegrationInfo<dim, spacedim> &integration_info_1,
663 IntegrationInfo<dim, spacedim> &integration_info_2) {
664 integrator.face(dof_info_1,
665 dof_info_2,
666 integration_info_1,
667 integration_info_2);
668 };
669
670 loop<dim, spacedim>(begin,
671 end,
672 dof_info,
673 box,
674 cell_worker,
675 boundary_worker,
676 face_worker,
677 assembler,
678 lctrl);
679 }
680
681} // namespace MeshWorker
682
684
685#endif
DOFINFO interior[GeometryInfo< dim >::faces_per_cell]
Definition dof_info.h:263
void assemble(ASSEMBLER &ass) const
Definition dof_info.h:492
bool exterior_face_available[GeometryInfo< dim >::faces_per_cell]
Definition dof_info.h:279
bool interior_face_available[GeometryInfo< dim >::faces_per_cell]
Definition dof_info.h:273
DOFINFO exterior[GeometryInfo< dim >::faces_per_cell]
Definition dof_info.h:267
virtual void cell(DoFInfo< dim, spacedim, number > &dinfo, IntegrationInfo< dim, spacedim > &info) const
virtual void face(DoFInfo< dim, spacedim, number > &dinfo1, DoFInfo< dim, spacedim, number > &dinfo2, IntegrationInfo< dim, spacedim > &info1, IntegrationInfo< dim, spacedim > &info2) const
virtual void boundary(DoFInfo< dim, spacedim, number > &dinfo, IntegrationInfo< dim, spacedim > &info) const
FaceOption own_faces
Definition loop.h:274
FaceOption faces_to_ghost
Definition loop.h:263
#define DEAL_II_DEPRECATED
Definition config.h:205
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
void cell_action(IteratorType 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:316
void integration_loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DoFInfo< dim, spacedim > &dof_info, IntegrationInfoBox< dim, spacedim > &box, const LocalIntegrator< dim, spacedim > &integrator, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:627
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:564
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)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition loop.h:70
bool is_active_iterator(const DI &)
Definition loop.h:48
const types::subdomain_id artificial_subdomain_id
Definition types.h:362
const types::subdomain_id invalid_subdomain_id
Definition types.h:341
typename type_identity< T >::type type_identity_t
Definition type_traits.h:95
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()