dof_info.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2009 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
 
#ifndef dealii_mesh_worker_dof_info_h
#define dealii_mesh_worker_dof_info_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/quadrature_lib.h>
 
#include <deal.II/dofs/block_info.h>
 
#include <deal.II/fe/fe_values.h>
 
#include <deal.II/meshworker/local_results.h>
#include <deal.II/meshworker/vector_selector.h>
 
#include <memory>
 
DEAL_II_NAMESPACE_OPEN
 
namespace MeshWorker
{
  // Forward declaration
#ifndef DOXYGEN
  template <int dim, class DOFINFO>
  class DoFInfoBox;
#endif

  template <int dim, int spacedim = dim, typename number = double>
  class DEAL_II_DEPRECATED_EARLY DoFInfo : public LocalResults<number>
  {
  public:
    typename Triangulation<dim, spacedim>::cell_iterator cell;

    typename Triangulation<dim, spacedim>::face_iterator face;

    unsigned int face_number;

    unsigned int sub_number;

    std::vector<types::global_dof_index> indices;

    std::vector<std::vector<types::global_dof_index>> indices_by_block;

    DoFInfo(const BlockInfo &block_info);

    DoFInfo(const DoFHandler<dim, spacedim> &dof_handler);

    template <class DHCellIterator>
    void
    reinit(const DHCellIterator &c);

    template <class DHCellIterator, class DHFaceIterator>
    void
    reinit(const DHCellIterator &c,
           const DHFaceIterator &f,
           const unsigned int    face_no);

    template <class DHCellIterator, class DHFaceIterator>
    void
    reinit(const DHCellIterator &c,
           const DHFaceIterator &f,
           const unsigned int    face_no,
           const unsigned int    subface_no);

    template <class DHFaceIterator>
    void
    set_face(const DHFaceIterator &f, const unsigned int face_no);

    template <class DHFaceIterator>
    void
    set_subface(const DHFaceIterator &f,
                const unsigned int    face_no,
                const unsigned int    subface_no);
 
    const BlockIndices &
    local_indices() const;
 

    ObserverPointer<const BlockInfo, DoFInfo<dim, spacedim>> block_info;

    bool level_cell;
 
  private:
    DoFInfo();

    void
    set_block_indices();

    template <class DHCellIterator>
    void
    get_indices(const DHCellIterator &c);

    std::vector<types::global_dof_index> indices_org;

    BlockIndices aux_local_indices;
 
    friend class DoFInfoBox<dim, DoFInfo<dim, spacedim, number>>;
  };
 

  template <int dim, class DOFINFO>
  class DEAL_II_DEPRECATED_EARLY DoFInfoBox
  {
  public:
    DoFInfoBox(const DOFINFO &seed);

    DoFInfoBox(const DoFInfoBox<dim, DOFINFO> &);

    DoFInfoBox &
    operator=(const DoFInfoBox<dim, DOFINFO> &);

    void
    reset();

    template <class ASSEMBLER>
    void
    assemble(ASSEMBLER &ass) const;
 

    DOFINFO cell;
    DOFINFO interior[GeometryInfo<dim>::faces_per_cell];
    DOFINFO exterior[GeometryInfo<dim>::faces_per_cell];

    bool interior_face_available[GeometryInfo<dim>::faces_per_cell];

    bool exterior_face_available[GeometryInfo<dim>::faces_per_cell];

    bool cell_valid;
  };
 
  //----------------------------------------------------------------------//
 
  template <int dim, int spacedim, typename number>
  DoFInfo<dim, spacedim, number>::DoFInfo()
    : face_number(numbers::invalid_unsigned_int)
    , sub_number(numbers::invalid_unsigned_int)
    , level_cell(false)
  {}
 
 
 
  template <int dim, int spacedim, typename number>
  DoFInfo<dim, spacedim, number>::DoFInfo(
    const DoFHandler<dim, spacedim> &dof_handler)
    : face_number(numbers::invalid_unsigned_int)
    , sub_number(numbers::invalid_unsigned_int)
    , level_cell(false)
  {
    std::vector<types::global_dof_index> aux(1);
    aux[0] = dof_handler.get_fe().n_dofs_per_cell();
    aux_local_indices.reinit(aux);
  }
 
 
  template <int dim, int spacedim, typename number>
  template <class DHCellIterator>
  inline void
  DoFInfo<dim, spacedim, number>::get_indices(const DHCellIterator &c)
  {
    indices.resize(c->get_fe().n_dofs_per_cell());
    if (block_info == nullptr || block_info->local().size() == 0)
      c->get_active_or_mg_dof_indices(indices);
    else
      {
        indices_org.resize(c->get_fe().n_dofs_per_cell());
        c->get_active_or_mg_dof_indices(indices_org);
        set_block_indices();
      }
  }
 
 
  template <int dim, int spacedim, typename number>
  template <class DHCellIterator>
  inline void
  DoFInfo<dim, spacedim, number>::reinit(const DHCellIterator &c)
  {
    get_indices(c);
    level_cell = c->is_level_cell();
 
    cell        = typename Triangulation<dim, spacedim>::cell_iterator(*c);
    face_number = numbers::invalid_unsigned_int;
    sub_number  = numbers::invalid_unsigned_int;
    if (block_info)
      LocalResults<number>::reinit(block_info->local());
    else
      LocalResults<number>::reinit(aux_local_indices);
  }
 
 
  template <int dim, int spacedim, typename number>
  template <class DHFaceIterator>
  inline void
  DoFInfo<dim, spacedim, number>::set_face(const DHFaceIterator &f,
                                           const unsigned int    face_no)
  {
    face = static_cast<typename Triangulation<dim, spacedim>::face_iterator>(f);
    face_number = face_no;
    sub_number  = numbers::invalid_unsigned_int;
  }
 
 
  template <int dim, int spacedim, typename number>
  template <class DHCellIterator, class DHFaceIterator>
  inline void
  DoFInfo<dim, spacedim, number>::reinit(const DHCellIterator &c,
                                         const DHFaceIterator &f,
                                         const unsigned int    face_no)
  {
    if ((cell.state() != IteratorState::valid) ||
        cell != typename Triangulation<dim, spacedim>::cell_iterator(*c))
      get_indices(c);
    level_cell = c->is_level_cell();
 
    cell = typename Triangulation<dim, spacedim>::cell_iterator(*c);
    set_face(f, face_no);
 
    if (block_info)
      LocalResults<number>::reinit(block_info->local());
    else
      LocalResults<number>::reinit(aux_local_indices);
  }
 
 
  template <int dim, int spacedim, typename number>
  template <class DHFaceIterator>
  inline void
  DoFInfo<dim, spacedim, number>::set_subface(const DHFaceIterator &f,
                                              const unsigned int    face_no,
                                              const unsigned int    subface_no)
  {
    face = static_cast<typename Triangulation<dim, spacedim>::face_iterator>(f);
    face_number = face_no;
    sub_number  = subface_no;
  }
 
 
  template <int dim, int spacedim, typename number>
  template <class DHCellIterator, class DHFaceIterator>
  inline void
  DoFInfo<dim, spacedim, number>::reinit(const DHCellIterator &c,
                                         const DHFaceIterator &f,
                                         const unsigned int    face_no,
                                         const unsigned int    subface_no)
  {
    if (cell.state() != IteratorState::valid ||
        cell !=
          static_cast<typename Triangulation<dim, spacedim>::cell_iterator>(c))
      get_indices(c);
    level_cell = c->is_level_cell();
 
    cell = static_cast<typename Triangulation<dim, spacedim>::cell_iterator>(c);
    set_subface(f, face_no, subface_no);
 
    if (block_info)
      LocalResults<number>::reinit(block_info->local());
    else
      LocalResults<number>::reinit(aux_local_indices);
  }
 
 
  template <int dim, int spacedim, typename number>
  inline const BlockIndices &
  DoFInfo<dim, spacedim, number>::local_indices() const
  {
    if (block_info)
      return block_info->local();
    return aux_local_indices;
  }
 
  //----------------------------------------------------------------------//
 
  template <int dim, int spacedim, typename number>
  inline DoFInfo<dim, spacedim, number>::DoFInfo(const BlockInfo &info)
    : face_number(numbers::invalid_unsigned_int)
    , sub_number(numbers::invalid_unsigned_int)
    , block_info(&info, typeid(*this).name())
    , level_cell(false)
  {
    indices_by_block.resize(info.local().size());
    for (unsigned int i = 0; i < indices_by_block.size(); ++i)
      indices_by_block[i].resize(info.local().block_size(i));
  }
 
 
 
  template <int dim, int spacedim, typename number>
  inline void
  DoFInfo<dim, spacedim, number>::set_block_indices()
  {
    for (unsigned int i = 0; i < indices.size(); ++i)
      {
        const std::pair<unsigned int, unsigned int> bi =
          block_info->local().global_to_local(this->block_info->renumber(i));
        indices_by_block[bi.first][bi.second] = indices_org[i];
      }
    // Remove this after
    // changing block codes
    for (unsigned int i = 0; i < indices.size(); ++i)
      indices[this->block_info->renumber(i)] = indices_org[i];
  }
 
 
 
  template <int dim, class DOFINFO>
  inline DoFInfoBox<dim, DOFINFO>::DoFInfoBox(const DOFINFO &seed)
    : cell(seed)
    , cell_valid(true)
  {
    for (const unsigned int i : GeometryInfo<dim>::face_indices())
      {
        exterior[i]                = seed;
        interior[i]                = seed;
        interior_face_available[i] = false;
        exterior_face_available[i] = false;
      }
  }
 
 
 
  template <int dim, class DOFINFO>
  inline DoFInfoBox<dim, DOFINFO>::DoFInfoBox(
    const DoFInfoBox<dim, DOFINFO> &other)
    : cell(other.cell)
    , cell_valid(other.cell_valid)
  {
    for (const unsigned int i : GeometryInfo<dim>::face_indices())
      {
        exterior[i]                = other.exterior[i];
        interior[i]                = other.interior[i];
        interior_face_available[i] = false;
        exterior_face_available[i] = false;
      }
  }
 
 
 
  template <int dim, class DOFINFO>
  inline DoFInfoBox<dim, DOFINFO> &
  DoFInfoBox<dim, DOFINFO>::operator=(const DoFInfoBox<dim, DOFINFO> &other)
  {
    cell       = other.cell;
    cell_valid = other.cell_valid;
    for (const unsigned int i : GeometryInfo<dim>::face_indices())
      {
        exterior[i]                = other.exterior[i];
        interior[i]                = other.interior[i];
        interior_face_available[i] = false;
        exterior_face_available[i] = false;
      }
    return *this;
  }
 
 
 
  template <int dim, class DOFINFO>
  inline void
  DoFInfoBox<dim, DOFINFO>::reset()
  {
    cell_valid = false;
    for (const unsigned int i : GeometryInfo<dim>::face_indices())
      {
        interior_face_available[i] = false;
        exterior_face_available[i] = false;
      }
  }
 
 
 
  template <int dim, class DOFINFO>
  template <class ASSEMBLER>
  inline void
  DoFInfoBox<dim, DOFINFO>::assemble(ASSEMBLER &assembler) const
  {
    if (!cell_valid)
      return;
 
    assembler.assemble(cell);
    for (const unsigned int i : GeometryInfo<dim>::face_indices())
      {
        // Only do something if data available
        if (interior_face_available[i])
          {
            // If both data
            // available, it is an
            // interior face
            if (exterior_face_available[i])
              assembler.assemble(interior[i], exterior[i]);
            else
              assembler.assemble(interior[i]);
          }
      }
  }
} // namespace MeshWorker
 
DEAL_II_NAMESPACE_CLOSE
 
#endif