loop.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2009 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
 
#ifndef dealii_mesh_worker_loop_h
#define dealii_mesh_worker_loop_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/template_constraints.h>
#include <deal.II/base/work_stream.h>
 
#include <deal.II/grid/filtered_iterator.h>
#include <deal.II/grid/tria.h>
 
#include <deal.II/meshworker/dof_info.h>
#include <deal.II/meshworker/integration_info.h>
#include <deal.II/meshworker/local_integrator.h>
 
#include <functional>
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declaration
#ifndef DOXYGEN
template <typename>
class TriaActiveIterator;
#endif
 
namespace internal
{
  template <class DI>
  inline bool
  is_active_iterator(const DI &)
  {
    return false;
  }
 
  template <typename AccessorType>
  inline bool
  is_active_iterator(const TriaActiveIterator<AccessorType> &)
  {
    return true;
  }
 
  template <typename AccessorType>
  inline bool
  is_active_iterator(
    const ::FilteredIterator<TriaActiveIterator<AccessorType>> &)
  {
    return true;
  }
 
  template <int dim, class DOFINFO, class A>
  void
  assemble(const MeshWorker::DoFInfoBox<dim, DOFINFO> &dinfo, A *assembler)
  {
    dinfo.assemble(*assembler);
  }
} // namespace internal
 
 
 
namespace MeshWorker
{
  class LoopControl
  {
  public:
    LoopControl()
      : own_cells(true)
      , ghost_cells(false)
      , faces_to_ghost(LoopControl::one)
      , own_faces(LoopControl::one)
      , cells_first(true)
    {}
 
    bool own_cells;
 
    bool ghost_cells;
 
    enum FaceOption
    {
      never,
      one,
      both
    };
 
    FaceOption faces_to_ghost;
 
    FaceOption own_faces;
 
    bool cells_first;
  };
 
 
 
  template <class INFOBOX,
            class DOFINFO,
            int dim,
            int spacedim,
            typename IteratorType>
  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)
  {
    const bool ignore_subdomain =
      (cell->get_triangulation().locally_owned_subdomain() ==
       numbers::invalid_subdomain_id);
 
    types::subdomain_id csid = (cell->is_level_cell()) ?
                                 cell->level_subdomain_id() :
                                 cell->subdomain_id();
 
    const bool own_cell =
      ignore_subdomain ||
      (csid == cell->get_triangulation().locally_owned_subdomain());
 
    dof_info.reset();
 
    if ((!ignore_subdomain) && (csid == numbers::artificial_subdomain_id))
      return;
 
    dof_info.cell.reinit(cell);
    dof_info.cell_valid = true;
 
    const bool integrate_cell          = (cell_worker != nullptr);
    const bool integrate_boundary      = (boundary_worker != nullptr);
    const bool integrate_interior_face = (face_worker != nullptr);
 
    if (integrate_cell)
      info.cell.reinit(dof_info.cell);
    // Execute this, if cells
    // have to be dealt with
    // before faces
    if (integrate_cell && loop_control.cells_first &&
        ((loop_control.own_cells && own_cell) ||
         (loop_control.ghost_cells && !own_cell)))
      cell_worker(dof_info.cell, info.cell);
 
    // Call the callback function in
    // the info box to do
    // computations between cell and
    // face action.
    info.post_cell(dof_info);
 
    if (integrate_interior_face || integrate_boundary)
      for (const unsigned int face_no : cell->face_indices())
        {
          typename IteratorType::AccessorType::Container::face_iterator face =
            cell->face(face_no);
          if (cell->at_boundary(face_no) &&
              !cell->has_periodic_neighbor(face_no))
            {
              // only integrate boundary faces of own cells
              if (integrate_boundary && own_cell)
                {
                  dof_info.interior_face_available[face_no] = true;
                  dof_info.interior[face_no].reinit(cell, face, face_no);
                  info.boundary.reinit(dof_info.interior[face_no]);
                  boundary_worker(dof_info.interior[face_no], info.boundary);
                }
            }
          else if (integrate_interior_face)
            {
              // Interior face
              TriaIterator<typename IteratorType::AccessorType> neighbor =
                cell->neighbor_or_periodic_neighbor(face_no);
 
              types::subdomain_id neighbid = numbers::artificial_subdomain_id;
              if (neighbor->is_level_cell())
                neighbid = neighbor->level_subdomain_id();
              // subdomain id is only valid for active cells
              else if (neighbor->is_active())
                neighbid = neighbor->subdomain_id();
 
              const bool own_neighbor =
                ignore_subdomain ||
                (neighbid ==
                 cell->get_triangulation().locally_owned_subdomain());
 
              // skip all faces between two ghost cells
              if (!own_cell && !own_neighbor)
                continue;
 
              // skip if the user doesn't want faces between own cells
              if (own_cell && own_neighbor &&
                  loop_control.own_faces == LoopControl::never)
                continue;
 
              // skip face to ghost
              if (own_cell != own_neighbor &&
                  loop_control.faces_to_ghost == LoopControl::never)
                continue;
 
              // Deal with refinement edges from the refined side. Assuming
              // one-irregular meshes, this situation should only occur if both
              // cells are active.
              const bool periodic_neighbor =
                cell->has_periodic_neighbor(face_no);
 
              if ((!periodic_neighbor && cell->neighbor_is_coarser(face_no)) ||
                  (periodic_neighbor &&
                   cell->periodic_neighbor_is_coarser(face_no)))
                {
                  Assert(cell->is_active(), ExcInternalError());
                  Assert(neighbor->is_active(), ExcInternalError());
 
                  // skip if only one processor needs to assemble the face
                  // to a ghost cell and the fine cell is not ours.
                  if (!own_cell &&
                      loop_control.faces_to_ghost == LoopControl::one)
                    continue;
 
                  const std::pair<unsigned int, unsigned int> neighbor_face_no =
                    periodic_neighbor ?
                      cell->periodic_neighbor_of_coarser_periodic_neighbor(
                        face_no) :
                      cell->neighbor_of_coarser_neighbor(face_no);
                  const typename IteratorType::AccessorType::Container::
                    face_iterator nface =
                      neighbor->face(neighbor_face_no.first);
 
                  dof_info.interior_face_available[face_no] = true;
                  dof_info.exterior_face_available[face_no] = true;
                  dof_info.interior[face_no].reinit(cell, face, face_no);
                  info.face.reinit(dof_info.interior[face_no]);
                  dof_info.exterior[face_no].reinit(neighbor,
                                                    nface,
                                                    neighbor_face_no.first,
                                                    neighbor_face_no.second);
                  info.subface.reinit(dof_info.exterior[face_no]);
 
                  face_worker(dof_info.interior[face_no],
                              dof_info.exterior[face_no],
                              info.face,
                              info.subface);
                }
              else
                {
                  // If iterator is active and neighbor is refined, skip
                  // internal face.
                  if (::internal::is_active_iterator(cell) &&
                      neighbor->has_children())
                    {
                      Assert(
                        loop_control.own_faces != LoopControl::both,
                        ExcMessage(
                          "Assembling from both sides for own_faces is not "
                          "supported with hanging nodes!"));
                      continue;
                    }
 
                  // Now neighbor is on same level, double-check this:
                  Assert(cell->level() == neighbor->level(),
                         ExcInternalError());
 
                  // If we own both cells only do faces from one side (unless
                  // LoopControl says otherwise). Here, we rely on cell
                  // comparison that will look at cell->index().
                  if (own_cell && own_neighbor &&
                      loop_control.own_faces == LoopControl::one &&
                      (neighbor < cell))
                    continue;
 
                  // independent of loop_control.faces_to_ghost,
                  // we only look at faces to ghost on the same level once
                  // (only where own_cell=true and own_neighbor=false)
                  if (!own_cell)
                    continue;
 
                  // now only one processor assembles faces_to_ghost. We let the
                  // processor with the smaller (level-)subdomain id assemble
                  // the face.
                  if (own_cell && !own_neighbor &&
                      loop_control.faces_to_ghost == LoopControl::one &&
                      (neighbid < csid))
                    continue;
 
                  const unsigned int neighbor_face_no =
                    periodic_neighbor ?
                      cell->periodic_neighbor_face_no(face_no) :
                      cell->neighbor_face_no(face_no);
                  Assert(periodic_neighbor ||
                           neighbor->face(neighbor_face_no) == face,
                         ExcInternalError());
                  // Regular interior face
                  dof_info.interior_face_available[face_no] = true;
                  dof_info.exterior_face_available[face_no] = true;
                  dof_info.interior[face_no].reinit(cell, face, face_no);
                  info.face.reinit(dof_info.interior[face_no]);
                  dof_info.exterior[face_no].reinit(neighbor,
                                                    neighbor->face(
                                                      neighbor_face_no),
                                                    neighbor_face_no);
                  info.neighbor.reinit(dof_info.exterior[face_no]);
 
                  face_worker(dof_info.interior[face_no],
                              dof_info.exterior[face_no],
                              info.face,
                              info.neighbor);
                }
            }
        } // faces
    // Call the callback function in
    // the info box to do
    // computations between face and
    // cell action.
    info.post_faces(dof_info);
 
    // Execute this, if faces
    // have to be handled first
    if (integrate_cell && !loop_control.cells_first &&
        ((loop_control.own_cells && own_cell) ||
         (loop_control.ghost_cells && !own_cell)))
      cell_worker(dof_info.cell, info.cell);
  }
 
 
  template <int dim,
            int spacedim,
            class DOFINFO,
            class INFOBOX,
            typename AssemblerType,
            typename IteratorType>
  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())
  {
    DoFInfoBox<dim, DOFINFO> dof_info(dinfo);
 
    assembler.initialize_info(dof_info.cell, false);
    for (const unsigned int i : GeometryInfo<dim>::face_indices())
      {
        assembler.initialize_info(dof_info.interior[i], true);
        assembler.initialize_info(dof_info.exterior[i], true);
      }
 
    // Loop over all cells
    WorkStream::run(
      begin,
      end,
      [&cell_worker, &boundary_worker, &face_worker, &lctrl](
        IteratorType cell, INFOBOX &info, DoFInfoBox<dim, DOFINFO> &dof_info) {
        cell_action<INFOBOX, DOFINFO, dim, spacedim, IteratorType>(
          cell,
          dof_info,
          info,
          cell_worker,
          boundary_worker,
          face_worker,
          lctrl);
      },
      [&assembler](const MeshWorker::DoFInfoBox<dim, DOFINFO> &dinfo) {
        ::internal::assemble<dim, DOFINFO, AssemblerType>(dinfo,
                                                                &assembler);
      },
      info,
      dof_info);
  }
 
 
  template <int dim,
            int spacedim,
            typename IteratorType,
            typename AssemblerType>
  DEAL_II_DEPRECATED_EARLY 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())
  {
    std::function<void(DoFInfo<dim, spacedim> &,
                       IntegrationInfo<dim, spacedim> &)>
      cell_worker;
    std::function<void(DoFInfo<dim, spacedim> &,
                       IntegrationInfo<dim, spacedim> &)>
      boundary_worker;
    std::function<void(DoFInfo<dim, spacedim> &,
                       DoFInfo<dim, spacedim> &,
                       IntegrationInfo<dim, spacedim> &,
                       IntegrationInfo<dim, spacedim> &)>
      face_worker;
    if (integrator.use_cell)
      cell_worker =
        [&integrator](DoFInfo<dim, spacedim>         &dof_info,
                      IntegrationInfo<dim, spacedim> &integration_info) {
          integrator.cell(dof_info, integration_info);
        };
    if (integrator.use_boundary)
      boundary_worker =
        [&integrator](DoFInfo<dim, spacedim>         &dof_info,
                      IntegrationInfo<dim, spacedim> &integration_info) {
          integrator.boundary(dof_info, integration_info);
        };
    if (integrator.use_face)
      face_worker =
        [&integrator](DoFInfo<dim, spacedim>         &dof_info_1,
                      DoFInfo<dim, spacedim>         &dof_info_2,
                      IntegrationInfo<dim, spacedim> &integration_info_1,
                      IntegrationInfo<dim, spacedim> &integration_info_2) {
          integrator.face(dof_info_1,
                          dof_info_2,
                          integration_info_1,
                          integration_info_2);
        };
 
    loop<dim, spacedim>(begin,
                        end,
                        dof_info,
                        box,
                        cell_worker,
                        boundary_worker,
                        face_worker,
                        assembler,
                        lctrl);
  }
 
} // namespace MeshWorker
 
DEAL_II_NAMESPACE_CLOSE
 
#endif