integration_info.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_integration_info_h
#define dealii_mesh_worker_integration_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/dof_info.h>
#include <deal.II/meshworker/local_results.h>
#include <deal.II/meshworker/vector_selector.h>
 
#include <memory>
 
DEAL_II_NAMESPACE_OPEN
 
namespace MeshWorker
{
  template <int dim, int spacedim = dim>
  class IntegrationInfo
  {
  private:
    std::vector<std::shared_ptr<FEValuesBase<dim, spacedim>>> fevalv;
 
  public:
    static constexpr unsigned int dimension       = dim;
    static constexpr unsigned int space_dimension = spacedim;
 
    IntegrationInfo();
 
    IntegrationInfo(const IntegrationInfo<dim, spacedim> &other);
 
    template <class FEVALUES>
    void
    initialize(const FiniteElement<dim, spacedim>             &el,
               const Mapping<dim, spacedim>                   &mapping,
               const Quadrature<FEVALUES::integral_dimension> &quadrature,
               const UpdateFlags                               flags,
               const BlockInfo *local_block_info = nullptr);
 
    void
    initialize_data(const std::shared_ptr<VectorDataBase<dim, spacedim>> &data);
 
    void
    clear();
 
    const FiniteElement<dim, spacedim> &
    finite_element() const;
 
    bool multigrid;
 
    const FEValuesBase<dim, spacedim> &
    fe_values() const;
 
 
    const FEValuesBase<dim, spacedim> &
    fe_values(const unsigned int i) const;
 
    std::vector<std::vector<std::vector<double>>> values;
 
    std::vector<std::vector<std::vector<Tensor<1, spacedim>>>> gradients;
 
    std::vector<std::vector<std::vector<Tensor<2, spacedim>>>> hessians;
 
    template <typename number>
    void
    reinit(const DoFInfo<dim, spacedim, number> &i);
 
    template <typename number>
    void
    fill_local_data(const DoFInfo<dim, spacedim, number> &info,
                    bool                                  split_fevalues);
 
    std::shared_ptr<VectorDataBase<dim, spacedim>> global_data;
 
    std::size_t
    memory_consumption() const;
 
  private:
    SmartPointer<const FiniteElement<dim, spacedim>,
                 IntegrationInfo<dim, spacedim>>
      fe_pointer;
 
    template <typename TYPE>
    void
    fill_local_data(std::vector<std::vector<std::vector<TYPE>>> &data,
                    VectorSelector                              &selector,
                    bool split_fevalues) const;
    unsigned int n_components;
  };
 
  template <int dim, int spacedim = dim>
  class IntegrationInfoBox
  {
  public:
    using CellInfo = IntegrationInfo<dim, spacedim>;
 
    IntegrationInfoBox();
 
    void
    initialize(const FiniteElement<dim, spacedim> &el,
               const Mapping<dim, spacedim>       &mapping,
               const BlockInfo                    *block_info = nullptr);
 
    template <typename VectorType>
    void
    initialize(const FiniteElement<dim, spacedim> &el,
               const Mapping<dim, spacedim>       &mapping,
               const AnyData                      &data,
               const VectorType                   &dummy,
               const BlockInfo                    *block_info = nullptr);
    template <typename VectorType>
    void
    initialize(const FiniteElement<dim, spacedim> &el,
               const Mapping<dim, spacedim>       &mapping,
               const AnyData                      &data,
               const MGLevelObject<VectorType>    &dummy,
               const BlockInfo                    *block_info = nullptr);
    void
    initialize_update_flags(bool neighbor_geometry = false);
 
    void
    add_update_flags_all(const UpdateFlags flags);
 
    void
    add_update_flags_cell(const UpdateFlags flags);
 
    void
    add_update_flags_boundary(const UpdateFlags flags);
 
    void
    add_update_flags_face(const UpdateFlags flags);
 
    void
    add_update_flags(const UpdateFlags flags,
                     const bool        cell     = true,
                     const bool        boundary = true,
                     const bool        face     = true,
                     const bool        neighbor = true);
 
    void
    initialize_gauss_quadrature(unsigned int n_cell_points,
                                unsigned int n_boundary_points,
                                unsigned int n_face_points,
                                const bool   force = true);
 
    std::size_t
    memory_consumption() const;
 
    UpdateFlags cell_flags;
    UpdateFlags boundary_flags;
 
    UpdateFlags face_flags;
 
    UpdateFlags neighbor_flags;
 
    Quadrature<dim> cell_quadrature;
 
    Quadrature<dim - 1> boundary_quadrature;
 
    Quadrature<dim - 1> face_quadrature;
    MeshWorker::VectorSelector cell_selector;
 
    MeshWorker::VectorSelector boundary_selector;
 
    MeshWorker::VectorSelector face_selector;
 
    std::shared_ptr<MeshWorker::VectorDataBase<dim, spacedim>> cell_data;
    std::shared_ptr<MeshWorker::VectorDataBase<dim, spacedim>> boundary_data;
    std::shared_ptr<MeshWorker::VectorDataBase<dim, spacedim>> face_data;
    template <class DOFINFO>
    void
    post_cell(const DoFInfoBox<dim, DOFINFO> &);
 
    template <class DOFINFO>
    void
    post_faces(const DoFInfoBox<dim, DOFINFO> &);
 
    CellInfo cell;
    CellInfo boundary;
    CellInfo face;
    CellInfo subface;
    CellInfo neighbor;
 
  };
 
 
  //----------------------------------------------------------------------//
 
  template <int dim, int sdim>
  inline IntegrationInfo<dim, sdim>::IntegrationInfo()
    : fevalv(0)
    , multigrid(false)
    , global_data(std::make_shared<VectorDataBase<dim, sdim>>())
    , n_components(numbers::invalid_unsigned_int)
  {}
 
#ifndef DOXYGEN
  template <int dim, int sdim>
  inline IntegrationInfo<dim, sdim>::IntegrationInfo(
    const IntegrationInfo<dim, sdim> &other)
    : multigrid(other.multigrid)
    , values(other.values)
    , gradients(other.gradients)
    , hessians(other.hessians)
    , global_data(other.global_data)
    , fe_pointer(other.fe_pointer)
    , n_components(other.n_components)
  {
    fevalv.resize(other.fevalv.size());
    for (unsigned int i = 0; i < other.fevalv.size(); ++i)
      {
        const FEValuesBase<dim, sdim> &p = *other.fevalv[i];
        const FEValues<dim, sdim>     *pc =
          dynamic_cast<const FEValues<dim, sdim> *>(&p);
        const FEFaceValues<dim, sdim> *pf =
          dynamic_cast<const FEFaceValues<dim, sdim> *>(&p);
        const FESubfaceValues<dim, sdim> *ps =
          dynamic_cast<const FESubfaceValues<dim, sdim> *>(&p);
 
        if (pc != nullptr)
          fevalv[i] =
            std::make_shared<FEValues<dim, sdim>>(pc->get_mapping(),
                                                  pc->get_fe(),
                                                  pc->get_quadrature(),
                                                  pc->get_update_flags());
        else if (pf != nullptr)
          fevalv[i] =
            std::make_shared<FEFaceValues<dim, sdim>>(pf->get_mapping(),
                                                      pf->get_fe(),
                                                      pf->get_quadrature(),
                                                      pf->get_update_flags());
        else if (ps != nullptr)
          fevalv[i] = std::make_shared<FESubfaceValues<dim, sdim>>(
            ps->get_mapping(),
            ps->get_fe(),
            ps->get_quadrature(),
            ps->get_update_flags());
        else
          DEAL_II_ASSERT_UNREACHABLE();
      }
  }
#endif
 
 
 
  template <int dim, int sdim>
  template <class FEVALUES>
  inline void
  IntegrationInfo<dim, sdim>::initialize(
    const FiniteElement<dim, sdim>                 &el,
    const Mapping<dim, sdim>                       &mapping,
    const Quadrature<FEVALUES::integral_dimension> &quadrature,
    const UpdateFlags                               flags,
    const BlockInfo                                *block_info)
  {
    fe_pointer = &el;
    if (block_info == nullptr || block_info->local().size() == 0)
      {
        fevalv.resize(1);
        fevalv[0] = std::make_shared<FEVALUES>(mapping, el, quadrature, flags);
      }
    else
      {
        fevalv.resize(el.n_base_elements());
        for (unsigned int i = 0; i < fevalv.size(); ++i)
          fevalv[i] = std::make_shared<FEVALUES>(mapping,
                                                 el.base_element(i),
                                                 quadrature,
                                                 flags);
      }
    n_components = el.n_components();
  }
 
 
  template <int dim, int spacedim>
  inline const FiniteElement<dim, spacedim> &
  IntegrationInfo<dim, spacedim>::finite_element() const
  {
    Assert(fe_pointer != nullptr, ExcNotInitialized());
    return *fe_pointer;
  }
 
  template <int dim, int spacedim>
  inline const FEValuesBase<dim, spacedim> &
  IntegrationInfo<dim, spacedim>::fe_values() const
  {
    AssertDimension(fevalv.size(), 1);
    return *fevalv[0];
  }
 
 
  template <int dim, int spacedim>
  inline const FEValuesBase<dim, spacedim> &
  IntegrationInfo<dim, spacedim>::fe_values(unsigned int i) const
  {
    AssertIndexRange(i, fevalv.size());
    return *fevalv[i];
  }
 
 
  template <int dim, int spacedim>
  template <typename number>
  inline void
  IntegrationInfo<dim, spacedim>::reinit(
    const DoFInfo<dim, spacedim, number> &info)
  {
    for (unsigned int i = 0; i < fevalv.size(); ++i)
      {
        FEValuesBase<dim, spacedim> &febase = *fevalv[i];
        if (info.sub_number != numbers::invalid_unsigned_int)
          {
            // This is a subface
            FESubfaceValues<dim, spacedim> &fe =
              dynamic_cast<FESubfaceValues<dim, spacedim> &>(febase);
            fe.reinit(info.cell, info.face_number, info.sub_number);
          }
        else if (info.face_number != numbers::invalid_unsigned_int)
          {
            // This is a face
            FEFaceValues<dim, spacedim> &fe =
              dynamic_cast<FEFaceValues<dim, spacedim> &>(febase);
            fe.reinit(info.cell, info.face_number);
          }
        else
          {
            // This is a cell
            FEValues<dim, spacedim> &fe =
              dynamic_cast<FEValues<dim, spacedim> &>(febase);
            fe.reinit(info.cell);
          }
      }
 
    const bool split_fevalues = info.block_info != nullptr;
    if (!global_data->empty())
      fill_local_data(info, split_fevalues);
  }
 
 
 
  //----------------------------------------------------------------------//
 
  template <int dim, int sdim>
  inline void
  IntegrationInfoBox<dim, sdim>::initialize_gauss_quadrature(unsigned int cp,
                                                             unsigned int bp,
                                                             unsigned int fp,
                                                             bool         force)
  {
    if (force || cell_quadrature.empty())
      cell_quadrature = QGauss<dim>(cp);
    if (force || boundary_quadrature.empty())
      boundary_quadrature = QGauss<dim - 1>(bp);
    if (force || face_quadrature.empty())
      face_quadrature = QGauss<dim - 1>(fp);
  }
 
 
  template <int dim, int sdim>
  inline void
  IntegrationInfoBox<dim, sdim>::add_update_flags_all(const UpdateFlags flags)
  {
    add_update_flags(flags, true, true, true, true);
  }
 
 
  template <int dim, int sdim>
  inline void
  IntegrationInfoBox<dim, sdim>::add_update_flags_cell(const UpdateFlags flags)
  {
    add_update_flags(flags, true, false, false, false);
  }
 
 
  template <int dim, int sdim>
  inline void
  IntegrationInfoBox<dim, sdim>::add_update_flags_boundary(
    const UpdateFlags flags)
  {
    add_update_flags(flags, false, true, false, false);
  }
 
 
  template <int dim, int sdim>
  inline void
  IntegrationInfoBox<dim, sdim>::add_update_flags_face(const UpdateFlags flags)
  {
    add_update_flags(flags, false, false, true, true);
  }
 
 
#ifndef DOXYGEN
  template <int dim, int sdim>
  inline void
  IntegrationInfoBox<dim, sdim>::initialize(const FiniteElement<dim, sdim> &el,
                                            const Mapping<dim, sdim> &mapping,
                                            const BlockInfo *block_info)
  {
    initialize_update_flags();
    initialize_gauss_quadrature((cell_flags & update_values) ?
                                  (el.tensor_degree() + 1) :
                                  el.tensor_degree(),
                                (boundary_flags & update_values) ?
                                  (el.tensor_degree() + 1) :
                                  el.tensor_degree(),
                                (face_flags & update_values) ?
                                  (el.tensor_degree() + 1) :
                                  el.tensor_degree(),
                                false);
 
    cell.template initialize<FEValues<dim, sdim>>(
      el, mapping, cell_quadrature, cell_flags, block_info);
    boundary.template initialize<FEFaceValues<dim, sdim>>(
      el, mapping, boundary_quadrature, boundary_flags, block_info);
    face.template initialize<FEFaceValues<dim, sdim>>(
      el, mapping, face_quadrature, face_flags, block_info);
    subface.template initialize<FESubfaceValues<dim, sdim>>(
      el, mapping, face_quadrature, face_flags, block_info);
    neighbor.template initialize<FEFaceValues<dim, sdim>>(
      el, mapping, face_quadrature, neighbor_flags, block_info);
  }
#endif
 
 
  template <int dim, int sdim>
  template <typename VectorType>
  void
  IntegrationInfoBox<dim, sdim>::initialize(const FiniteElement<dim, sdim> &el,
                                            const Mapping<dim, sdim> &mapping,
                                            const AnyData            &data,
                                            const VectorType &,
                                            const BlockInfo *block_info)
  {
    initialize(el, mapping, block_info);
    std::shared_ptr<VectorData<VectorType, dim, sdim>> p;
    VectorDataBase<dim, sdim>                         *pp;
 
    p = std::make_shared<VectorData<VectorType, dim, sdim>>(cell_selector);
    // Public member function of parent class was not found without
    // explicit cast
    pp = &*p;
    pp->initialize(data);
    cell_data = p;
    cell.initialize_data(p);
 
    p  = std::make_shared<VectorData<VectorType, dim, sdim>>(boundary_selector);
    pp = &*p;
    pp->initialize(data);
    boundary_data = p;
    boundary.initialize_data(p);
 
    p  = std::make_shared<VectorData<VectorType, dim, sdim>>(face_selector);
    pp = &*p;
    pp->initialize(data);
    face_data = p;
    face.initialize_data(p);
    subface.initialize_data(p);
    neighbor.initialize_data(p);
  }
 
  template <int dim, int sdim>
  template <typename VectorType>
  void
  IntegrationInfoBox<dim, sdim>::initialize(const FiniteElement<dim, sdim> &el,
                                            const Mapping<dim, sdim> &mapping,
                                            const AnyData            &data,
                                            const MGLevelObject<VectorType> &,
                                            const BlockInfo *block_info)
  {
    initialize(el, mapping, block_info);
    std::shared_ptr<MGVectorData<VectorType, dim, sdim>> p;
    VectorDataBase<dim, sdim>                           *pp;
 
    p = std::make_shared<MGVectorData<VectorType, dim, sdim>>(cell_selector);
    // Public member function of parent class was not found without
    // explicit cast
    pp = &*p;
    pp->initialize(data);
    cell_data = p;
    cell.initialize_data(p);
 
    p =
      std::make_shared<MGVectorData<VectorType, dim, sdim>>(boundary_selector);
    pp = &*p;
    pp->initialize(data);
    boundary_data = p;
    boundary.initialize_data(p);
 
    p  = std::make_shared<MGVectorData<VectorType, dim, sdim>>(face_selector);
    pp = &*p;
    pp->initialize(data);
    face_data = p;
    face.initialize_data(p);
    subface.initialize_data(p);
    neighbor.initialize_data(p);
  }
 
  template <int dim, int sdim>
  template <class DOFINFO>
  void
  IntegrationInfoBox<dim, sdim>::post_cell(const DoFInfoBox<dim, DOFINFO> &)
  {}
 
 
  template <int dim, int sdim>
  template <class DOFINFO>
  void
  IntegrationInfoBox<dim, sdim>::post_faces(const DoFInfoBox<dim, DOFINFO> &)
  {}
 
 
} // namespace MeshWorker
 
DEAL_II_NAMESPACE_CLOSE
 
#endif